I am trying to answer a question from Pinhiero and Bates Mixed Effects Models in S and S-Plus, explaining how random effects fail to confer any benefit over a gls model that has mixed effects.
This model contains random effects
library(nlme) fm1Ovar.lme <- lme( follicles ~ 1 + sin(2*pi*Time) + cos(2*pi*Time), data = Ovary, random = pdDiag(~sin(2*pi*Time)), corr = corARMA(p = 1, q = 1))
fm1Ovar.lme Linear mixed-effects model fit by REML Data: Ovary Log-restricted-likelihood: -771.9471 Fixed: follicles ~ 1 + sin(2 * pi * Time) + cos(2 * pi * Time) (Intercept) sin(2 * pi * Time) cos(2 * pi * Time) 12.1248728 -2.9198264 -0.8487095 Random effects: Formula: ~sin(2 * pi * Time) | Mare Structure: Diagonal (Intercept) sin(2 * pi * Time) Residual StdDev: 2.614435 1.004898 3.733423 Correlation Structure: ARMA(1,1) Formula: ~1 | Mare Parameter estimate(s): Phi1 Theta1 0.7868896 -0.2793591 Number of Observations: 308 Number of Groups: 11
Whereas this model does not
fm1Ovar.gls <- gls( follicles ~ 1 + sin(2*pi*Time) + cos(2*pi*Time), data = Ovary, corr = corARMA(form = ~ 1 | Mare, p = 1, q = 1) )
fm1Ovar.gls Generalized least squares fit by REML Model: follicles ~ 1 + sin(2 * pi * Time) + cos(2 * pi * Time) Data: Ovary Log-restricted-likelihood: -773.3402 Coefficients: (Intercept) sin(2 * pi * Time) cos(2 * pi * Time) 12.0587080 -2.8832412 -0.8035591 Correlation Structure: ARMA(1,1) Formula: ~1 | Mare Parameter estimate(s): Phi1 Theta1 0.8908103 -0.3496050 Degrees of freedom: 308 total; 305 residual Residual standard error: 4.597197
Comparing log likelihood ratios
anova(fm1Ovar.lme, fm1Ovar.gls) Model df AIC BIC logLik Test L.Ratio p-value fm1Ovar.lme 1 8 1559.894 1589.657 -771.9471 fm1Ovar.gls 2 6 1558.680 1581.002 -773.3402 1 vs 2 2.786262 0.2483
Indicates that the
lme model with the random effects offers no significant benefit over the
gls model without. So I guess we would favour the more parsimonious
gls model with fewer parameters.
I was wondering how one would explain this? Presumably the absence of any benefit from including random effects in the model says something about those random effects, but what? Is it that the between-Mare (i.e. between-subject) variation is very low?
I suppose a broader issue is that I don't understand what random effects actually are in this context. I get that with the
gls model we are modeling correlation between observations within groups (Mares in this case), similar to what occurs in a mixed effects model. So what extra information does the random effect modeling give us that modeling correlation structure does not?