Specifying nested random effects in proc glimmix

Eagergrad — Tue, 21 Jan 2025 18:02:40 GMT

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.

topic Specifying nested random effects in proc glimmix in Statistical Procedures

Specifying nested random effects in proc glimmix