Is it possible to fit multivariate Gaussian models implied by mixed-models through generalised least squares in R, by using, for instance, the gls function?

For instance, the random intercept model via lme is

mod.lme <- lme(score~Machine, random = ~1|Worker, data=Machines)

s2.lme <- as.numeric(VarCorr(mod.lme)[2,1]) #residual variance 

s2.ranef.lme <- as.numeric(VarCorr(mod.lme)[1,1]) # ran. eff. variance 

(tot.var.lme <- s2.lme + s2.ranef.lme)  # sum of variance components 

The corresponding Gaussian model can be fitted by gls as follows:

mod.gls <- gls(score~Machine, correlation =  corCompSymm(form = ~1| Worker), 
(tot.var.gls <- as.numeric(exp(attributes(mod.gls$apVar)$Pars[2]))^2)

mSt <- mod.gls$modelStruct 
cSt <- mSt$corStruct 
(rho <- coef(cSt, unconstrained = FALSE))

all.equal(tot.var.gls, tot.var.lme)  # total variances equal? 

s2.ranef.gls <- rho*tot.var.gls # get variance of the ran. eff from gls 

all.equal(as.numeric(s2.ranef.gls), s2.ranef.lme) # is equal to lme ? 

But what about this ?

mod2.lme <- lme(score~Machine, random = ~1|Worker/Machine, data = Machines)

How would you fit it by gls ? Is it possible ?

  • 1
    $\begingroup$ All LMMs correspond to a multivariate normal model (while the converse is not true) with a structured variance covariance matrix, so "all" you have to do is to work out the marginal variance covariance matrix for the nested random-effect model and fit that - whether gls is then able to parameterize that model is then the next question. I haven't worked out the math, so I don't know, but my guess is that you may have to write your own corStruct class if you need to use gls. $\endgroup$ May 23, 2018 at 6:22
  • 1
    $\begingroup$ Also beware of the difference in parameter spaces: the parameter space for the compound symmetry model is bigger than it is for the random intercept model. The random-effect variance is necessarily non-negative which leads to a non-negative corr but the corr in the compound symmetry model can also be negative (though not too much). So while two model fits can be equivalent (if $\rho \ge 0$) they need not be (if $\rho < 0$) and strictly speaking the underlying models are not the same. $\endgroup$ May 23, 2018 at 6:35
  • 1
    $\begingroup$ Last comment :-) Your last model indicates that you are interested in a model with nested random effects but in the Machines data it is actually not the case that Machine is nested in Worker - rather, these variables are crossed. In fact the last model "tricks" lme to fit a model with a random main effect of Worker and the random interaction for Worker:Machine. Was that intentional? $\endgroup$ May 23, 2018 at 6:47
  • $\begingroup$ I completely agree on your comments. Actually, I'm interested in fitting the particular multivariate normal that corresponds to a given LMM and yes the issue of the different parameter spaces must be taken into account. The example with Machines data is just for illustration. Thus, the point with this question is how to trick gls in order to deliver the same thing we get from lme. $\endgroup$
    – utobi
    May 23, 2018 at 9:38
  • $\begingroup$ Maybe I can ask for more clarification at this point: (1) Are you primarily asking how to derive the Gaussian model that corresponds to a mixed-effect model or are you asking how to write a corStruct class in order to fit it with gls? (2) If you have a particular and not any mixed-model in mind, please specify it. For a general (any) mixed model of the form $Y = X\beta + Z b + e$ with $e \sim MVN(0, R)$ and $b \sim MVN(0, G)$ we have $Cov(Y) = Z G Z' + R$ which does not simplify and therefore infeasible to fit with gls. Some structure is needed. $\endgroup$ May 23, 2018 at 11:16


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.