BookmarkSubscribeRSS Feed
Eagergrad
Calcite | Level 5

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:

CaseDateEventFieldIDRotationRotation_RepTillageTillage_RepIntIntRepPrecipLogNLogP
16/9/201511OR1NT1ORxNT12.21  
26/9/201512HV1NT2HVxNT12.211.437511-0.10297
36/9/201513OR2RT1ORxRT12.211.153039-0.31723
46/9/201514HV2RT2HVxRT12.21  
56/9/201515HV3NT3HVxNT22.211.418637-0.08575
66/9/201516OR3NT4ORxNT22.211.6716780.288355
76/9/201517NC1RT3NCxRT12.211.415636-0.46679
86/9/201518NC2NT5NCxNT12.210.465066-0.99305
96/9/201519HV4NT6HVxNT32.211.459139-0.41273
106/9/2015110OR4NT7ORxNT32.21  
116/9/2015111NC3NT8NCxNT22.211.166647-0.26135
126/9/2015112NC4RT4NCxRT22.211.362267-0.27624
136/9/2015113OR5RT5ORxRT22.21-0.51493-2.25524
146/9/2015114HV5RT6HVxRT22.21-0.25198-3.78837
156/9/2015115NC5RT7NCxRT32.210.398079-1.14898
166/9/2015116OR6RT8ORxRT32.21  
176/9/2015117HV6RT9HVxRT32.211.455399-0.70826
186/9/2015118NC6NT9NCxNT32.211.585474-0.42879
196/25/201521OR1NT1ORxNT118.82  
206/25/201522HV1NT2HVxNT118.822.1723070.61322
216/25/201523OR2RT1ORxRT118.821.7686240.465857
226/25/201524HV2RT2HVxRT118.822.1852210.240139
236/25/201525HV3NT3HVxNT218.822.2447720.682816
246/25/201526OR3NT4ORxNT218.821.8085260.512942
256/25/201527NC1RT3NCxRT118.821.9494050.189349
266/25/201528NC2NT5NCxNT118.821.139256-0.12419
276/25/201529HV4NT6HVxNT318.822.2784090.245731
286/25/2015210OR4NT7ORxNT318.821.1060740.106727
296/25/2015211NC3NT8NCxNT218.821.5918560.387852
306/25/2015212NC4RT4NCxRT218.822.2728110.548714
316/25/2015213OR5RT5ORxRT218.821.8572570.093303
326/25/2015214HV5RT6HVxRT218.821.989452-0.54507
336/25/2015215NC5RT7NCxRT318.821.22797-0.10663
346/25/2015216OR6RT8ORxRT318.822.480516-0.79576
356/25/2015217HV6RT9HVxRT318.822.1472310.321941
366/25/2015218NC6NT9NCxNT318.821.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.

hackathon24-white-horiz.png

The 2025 SAS Hackathon Kicks Off on June 11!

Watch the live Hackathon Kickoff to get all the essential information about the SAS Hackathon—including how to join, how to participate, and expert tips for success.

YouTube LinkedIn

What is ANOVA?

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.

Discussion stats
  • 0 replies
  • 1521 views
  • 0 likes
  • 1 in conversation