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
  • 1,013
  • 10
  • 24
  • 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