2

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 ?

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 ?

utobi
  • 11,726
  • 1
    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. – Rune H Christensen May 23 '18 at 06:22
  • 1
    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. – Rune H Christensen May 23 '18 at 06:35
  • 1
    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? – Rune H Christensen May 23 '18 at 06:47
  • 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. – utobi May 23 '18 at 09:38
  • 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. – Rune H Christensen May 23 '18 at 11:16
  • So, about (1) I'm asking how to specify a corStruct and varClass such that the gls model fit is equal to that obtained by lme. In (2), I do have particular mixed-models in mind, such as random intercept, random slope and random interaction, with or without correlation between the ran.effs. But to keep it simple, I've asked here for the random interaction case. Do you have any idea how to trick gls in order to get a random-interaction mixed model from it, that is, to fit the specific Gaussian model that corresponds to this particular mixed-model? – utobi May 24 '18 at 14:57

0 Answers0