# Is it possible to fit mixed-models via gls?

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),
data=Machines)
(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 ?


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

How would you fit it by gls ? Is it possible ?
• 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. May 23, 2018 at 6:22
• 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. May 23, 2018 at 6:35
• 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? May 23, 2018 at 6:47
• 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. May 23, 2018 at 11:16