I am trying to fit a mixed model in r and get a variance grouped by a factor.
The initial model in SAS is:
proc glimmix data=Marsh.KISR_quad_mean (where=(eval=1)) plots=all ; *initial eval=0;
by eval;
class Site Plot Season Trt;
model Perc_HYAM= Site Season Site*Season
Trt Season*Trt Site*Trt Site*Trt*Season / ddfm=kr ;
random Plot*Site ;
random Plot*Site*Season ;
lsmeans Site Season Site*season Trt Season*Trt site*trt site*trt*season / cl ;
run;
The R version that works is:
Eval_Site_Fixed_NI_AOV <-function(dvar) paste(dvar,
"~ SITE + (1|BLOCK:SITE) +",
"SEASON + SITE:SEASON + (1|BLOCK:SITE:SEASON) +",
"TRT + SITE:TRT + SEASON:TRT + SITE:SEASON:TRT", sep=" ")
lmer(Eval_Site_Fixed_NI_AOV("perc_HYAM"),
data=KISR_Quad_means_E1)
The model with grouped variance in SAS is:
proc glimmix data=Marsh.KISR_quad_mean (where=(eval=1)) plots=all; *initial eval=0;
by eval;
class Site Plot Season Trt;
model Perc_HYAM= Site Season Site*Season
Trt Season*Trt Site*Trt Site*Trt*Season / ddfm=kr ;
random Plot*Site ;
**random Plot*Site*Season /group=Trt ; **
lsmeans Site Season Site*season Trt Season*Trt site*trt site*trt*season / ilink cl ;
run;
The random PlotSiteSeason is grouped by factor=Trt. My attempt, which is clearly wrong, was to use:
Eval_Site_TGroup_AOV <-function(dvar) paste(dvar,
"~ SITE + (1|BLOCK:SITE) +",
"SEASON + SITE:SEASON + (TRT|BLOCK:SITE:SEASON) +",
"TRT + SITE:TRT + SEASON:TRT + SITE:SEASON:TRT", sep=" ")
I have searched a few other packages that may allow specification of random effects but have not found a good example of how to equate what I do with SAS to r.
The random effects statement in r does not produce an separate estimate of variance for 2 treatment levels but a large matrix that exceeds the number of observations.