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.
... View more
