I suspect my issue is related to my random effects statement within the glimmix proc, but I don't know how to specify the structure/statement properly. Here's some background on the project where the variable names are in parentheses following their description.
We are investigating the amount of nitrate leaching (logN) under different crop rotations and tillage systems. Precipitation is associated with how much nitrate leaches from the soil, so I've included this as a quantitative predictor (precip). There are 3 levels of crop rotations (rotations) and two levls of tillage (tillage). The combination of these two effects yields 6 treatments (trt), and there are 3 replicates of each of these treatments (intrep), yielding a total of 18 plots (fieldID). When there is sufficient rain to get samples from at least 6 of the 18 plots, it is considered a rain event (event). Due to field conditions, we don't get samples from all 18 plots every time we sample, so there are missing data. I've read that proc glimmix can handle missing data better than proc mixed, so I've opted to use proc glimmix. I'm considering "event" a random effect because we are sampling the same plots throughout the study. Additionally, I'd like to considering "intrep" a random effect nested within "event" since there is variability across plots , and each time an "event" occurs, we get samples from the plots. I'm interested in determining if and how much of a difference there is in nitrate leaching between the levels of rotation, tillage, and their interaction.
Here is a sample of my data:
Case | Date | Event | FieldID | Rotation | Rotation_Rep | Tillage | Tillage_Rep | Int | IntRep | Precip | LogN | LogP |
1 | 6/9/2015 | 1 | 1 | OR | 1 | NT | 1 | ORxNT | 1 | 2.21 | ||
2 | 6/9/2015 | 1 | 2 | HV | 1 | NT | 2 | HVxNT | 1 | 2.21 | 1.437511 | -0.10297 |
3 | 6/9/2015 | 1 | 3 | OR | 2 | RT | 1 | ORxRT | 1 | 2.21 | 1.153039 | -0.31723 |
4 | 6/9/2015 | 1 | 4 | HV | 2 | RT | 2 | HVxRT | 1 | 2.21 | ||
5 | 6/9/2015 | 1 | 5 | HV | 3 | NT | 3 | HVxNT | 2 | 2.21 | 1.418637 | -0.08575 |
6 | 6/9/2015 | 1 | 6 | OR | 3 | NT | 4 | ORxNT | 2 | 2.21 | 1.671678 | 0.288355 |
7 | 6/9/2015 | 1 | 7 | NC | 1 | RT | 3 | NCxRT | 1 | 2.21 | 1.415636 | -0.46679 |
8 | 6/9/2015 | 1 | 8 | NC | 2 | NT | 5 | NCxNT | 1 | 2.21 | 0.465066 | -0.99305 |
9 | 6/9/2015 | 1 | 9 | HV | 4 | NT | 6 | HVxNT | 3 | 2.21 | 1.459139 | -0.41273 |
10 | 6/9/2015 | 1 | 10 | OR | 4 | NT | 7 | ORxNT | 3 | 2.21 | ||
11 | 6/9/2015 | 1 | 11 | NC | 3 | NT | 8 | NCxNT | 2 | 2.21 | 1.166647 | -0.26135 |
12 | 6/9/2015 | 1 | 12 | NC | 4 | RT | 4 | NCxRT | 2 | 2.21 | 1.362267 | -0.27624 |
13 | 6/9/2015 | 1 | 13 | OR | 5 | RT | 5 | ORxRT | 2 | 2.21 | -0.51493 | -2.25524 |
14 | 6/9/2015 | 1 | 14 | HV | 5 | RT | 6 | HVxRT | 2 | 2.21 | -0.25198 | -3.78837 |
15 | 6/9/2015 | 1 | 15 | NC | 5 | RT | 7 | NCxRT | 3 | 2.21 | 0.398079 | -1.14898 |
16 | 6/9/2015 | 1 | 16 | OR | 6 | RT | 8 | ORxRT | 3 | 2.21 | ||
17 | 6/9/2015 | 1 | 17 | HV | 6 | RT | 9 | HVxRT | 3 | 2.21 | 1.455399 | -0.70826 |
18 | 6/9/2015 | 1 | 18 | NC | 6 | NT | 9 | NCxNT | 3 | 2.21 | 1.585474 | -0.42879 |
19 | 6/25/2015 | 2 | 1 | OR | 1 | NT | 1 | ORxNT | 1 | 18.82 | ||
20 | 6/25/2015 | 2 | 2 | HV | 1 | NT | 2 | HVxNT | 1 | 18.82 | 2.172307 | 0.61322 |
21 | 6/25/2015 | 2 | 3 | OR | 2 | RT | 1 | ORxRT | 1 | 18.82 | 1.768624 | 0.465857 |
22 | 6/25/2015 | 2 | 4 | HV | 2 | RT | 2 | HVxRT | 1 | 18.82 | 2.185221 | 0.240139 |
23 | 6/25/2015 | 2 | 5 | HV | 3 | NT | 3 | HVxNT | 2 | 18.82 | 2.244772 | 0.682816 |
24 | 6/25/2015 | 2 | 6 | OR | 3 | NT | 4 | ORxNT | 2 | 18.82 | 1.808526 | 0.512942 |
25 | 6/25/2015 | 2 | 7 | NC | 1 | RT | 3 | NCxRT | 1 | 18.82 | 1.949405 | 0.189349 |
26 | 6/25/2015 | 2 | 8 | NC | 2 | NT | 5 | NCxNT | 1 | 18.82 | 1.139256 | -0.12419 |
27 | 6/25/2015 | 2 | 9 | HV | 4 | NT | 6 | HVxNT | 3 | 18.82 | 2.278409 | 0.245731 |
28 | 6/25/2015 | 2 | 10 | OR | 4 | NT | 7 | ORxNT | 3 | 18.82 | 1.106074 | 0.106727 |
29 | 6/25/2015 | 2 | 11 | NC | 3 | NT | 8 | NCxNT | 2 | 18.82 | 1.591856 | 0.387852 |
30 | 6/25/2015 | 2 | 12 | NC | 4 | RT | 4 | NCxRT | 2 | 18.82 | 2.272811 | 0.548714 |
31 | 6/25/2015 | 2 | 13 | OR | 5 | RT | 5 | ORxRT | 2 | 18.82 | 1.857257 | 0.093303 |
32 | 6/25/2015 | 2 | 14 | HV | 5 | RT | 6 | HVxRT | 2 | 18.82 | 1.989452 | -0.54507 |
33 | 6/25/2015 | 2 | 15 | NC | 5 | RT | 7 | NCxRT | 3 | 18.82 | 1.22797 | -0.10663 |
34 | 6/25/2015 | 2 | 16 | OR | 6 | RT | 8 | ORxRT | 3 | 18.82 | 2.480516 | -0.79576 |
35 | 6/25/2015 | 2 | 17 | HV | 6 | RT | 9 | HVxRT | 3 | 18.82 | 2.147231 | 0.321941 |
36 | 6/25/2015 | 2 | 18 | NC | 6 | NT | 9 | NCxNT | 3 | 18.82 | 1.634183 | -0.1134 |
I've provided (what I think is) the best iteration of the model below:
proc glimmix data=x.working plots=studentpanel;
class rotation tillage event intrep fieldID;
model logn = precip rotation tillage rotation*tillage /dist=normal link=identity zeta=1E-8;
random event;
lsmeans rotation|tillage / adjust=tukey lines;
lsmeans rotation|tillage / plots=meanplot (sliceby=tillage join);
run;
While this model doesn't produce any errors, warnings, or notes that I was getting with other models (e.g., estimated G matrix NPD, MIVQUE0 estimate of profiled variance is linearly related to other covariance parameters, Obtaining minimum variance quadratic unbiased estimates as starting values for the covariance parameters failed., etc.), I don't think this model captures the study design and complexity.
I have minimal experience using SAS and while I've read through a variety of SAS community posts and SAS publications to provide clarification, I find myself more confused as the alternative models I develop yield more errors than I initially had.
Any guidance or suggestions are appreciated.
ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.
Find more tutorials on the SAS Users YouTube channel.