Why does gls model without random effects yield a similar fit to mixed effects model?

Question

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))

Output

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) )

Output

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?

Answer 1

Random effects also model correlations. To explain this more formally, the model that both lme() and gls() are fitting is the following $$y_i = X_i\beta + \varepsilon_i, \quad \varepsilon_i \sim \mathcal N(0, V_i),$$ where $y_i$ denotes the outcome variable for the $i$-th group (i.e., follicles in your example, and the groups are indicated by the Mare variable), $X_i$ is the design matrix for the fixed effects $\beta$ (in your example $X_i$ has three columns, the intercept, sin(2 * pi * Time), and cos(2 * pi * Time)).

In your question, the focus is on $V_i$, which is the variance-covariance matrix for the error terms for the $i$-th group. You have the following options:

• When you use gls() you postulate that $V_i = \Sigma_i$, where the structure of $\Sigma_i$ is determined by the option you have selected in the correlation argument. I.e., you specified an ARMA model for it.
• When you use lme() in which you specify only the random argument but not the correlation argument, then you postulate that $V_i = Z_i D Z_i^\top + \sigma^2 \texttt{I}$, where $Z_i$ denotes the design matrix of the random effects, $D$ is the variance-covariance matrix of the random effects, and $\sigma^2$ is the variance of the within groups error terms with $\texttt I$ denoting the identity matrix. In your example, you specified that $Z_i$ has two columns (intercept and sin(2 * pi * Time)), and that the $2 \times 2$ covariance matrix of the random effects is diagonal with the covariance set to zero.
• When you use lme() in which you specify both the random and correlation argument, then you say that $V_i = Z_i D Z_i^\top + \Sigma_i$, where $\Sigma_i$ is the same matrix as in gls(). That is, you model the marginal covariance matrix $V_i$ using both random effects and a serial correlation structure. In you example, you combined the serial correlation structure mentioned in the first bullet with the random-effect structure mentioned in the second bullet.

The results of the likelihood ratio test between the two models suggest that the inclusion of random effects, i.e., the inclusion of the $Z_i D Z_i^\top$ term in the marginal covariance matrix does not improve the fit of the model over just using the ARMA serial correlation term $\Sigma_i$.