Fixing Variance Values in Lme4

Fixing variance values in lme4

This is possible, if a little hacky. Here's a reproducible example:

Fit the original model:

library(lme4)
set.seed(101)
ss <- sleepstudy[sample(nrow(sleepstudy),size=round(0.9*nrow(sleepstudy))),]
m1 <- lmer(Reaction~Days+(1|Subject)+(0+Days|Subject),ss)
fixef(m1)
## (Intercept)        Days 
##   251.55172    10.37874 

Recover the deviance (in this case REML-criterion) function:

dd <- as.function(m1)

I'm going to set the standard deviations to zero so that I have something to compare with, i.e. the coefficients of the regular linear model. (The parameter vector for dd is a vector containing the column-wise, lower-triangular, concatenated Cholesky factors for the scaled random effects terms in the model. Luckily, if all you have are scalar/intercept-only random effects (e.g. (1|x)), then these correspond to the random-effects standard deviations, scaled by the model standard deviation).

(ff <- dd(c(0,0)))  ## new REML: 1704.708
environment(dd)$pp$beta(1)  ## new parameters
##    [1] 251.11920  10.56979

Matches:

coef(lm(Reaction~Days,ss))
## (Intercept)        Days 
##   251.11920    10.56979 

If you want to construct a new merMod object you can do it as follows ...

opt <- list(par=c(0,0),fval=ff,conv=0)
lmod <- lFormula(Reaction~Days+(1|Subject)+(0+Days|Subject),ss)
m1X <- mkMerMod(environment(dd), opt, lmod$reTrms, fr = lmod$fr,
         mc = quote(hacked_lmer()))

Now suppose we want to set the variances to particular non-zero value (say (700,30)). This will be a little bit tricky because of the scaling by the residual standard deviation ...

newvar <- c(700,30)
ff2 <- dd(sqrt(newvar)/sigma(m1))
opt2 <- list(par=c(0,0),fval=ff,conv=0)
m2X <- mkMerMod(environment(dd), opt, lmod$reTrms, fr = lmod$fr,
         mc = quote(hacked_lmer()))
VarCorr(m2X)
unlist(VarCorr(m2X))
##   Subject Subject.1 
## 710.89304  30.46684

So this doesn't get us quite where we wanted (because the residual variance changes ...)

buildMM <- function(theta) {
   dd <- as.function(m1)
   ff <- dd(theta)
   opt <- list(par=c(0,0),fval=ff,conv=0)
   mm <- mkMerMod(environment(dd), opt, lmod$reTrms, fr = lmod$fr,
         mc = quote(hacked_lmer()))
   return(mm)
}

objfun <- function(x,target=c(700,30)) {
   mm <- buildMM(sqrt(x))
   return(sum((unlist(VarCorr(mm))-target)^2))
}
s0 <- c(700,30)/sigma(m1)^2
opt <- optim(fn=objfun,par=s0)
mm_final <- buildMM(sqrt(opt$par))
summary(mm_final)
##  Random effects:
##  Groups    Name        Variance Std.Dev.
##  Subject   (Intercept) 700      26.458  
##  Subject.1 Days         30       5.477  
##  Residual              700      26.458  
## Number of obs: 162, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  251.580      7.330   34.32
## Days          10.378      1.479    7.02

By the way, it's not generally recommended to use random effects when the grouping variables have a very small number (e.g. <5 or 6) levels: see here ...

texreg with lmer and lme objects -- variances differ

It seems pretty clear that it's reporting variance for model 1 and sigma for model 2. Variance is sigma squared, so the answer to your question is "both"

Equivalence of a mixed model fitted by lme and lmer

I can tell you what's going on and what should in principle fix it, but at the moment the fix doesn't work ...

The residuals being plotted take all of the random effects into account, which in the case of the lmer fit includes the individual-level random effects (the (0+dummy(Var4,"1")|Var5) term), which leads to weird residuals for the Var4==1 group. To illustrate this:

plot(ModLMER, col = Dataone$Var4+1)

i.e., you can see that the weird residuals are exactly the ones in red == those for which Var4==1.

In theory we should be able to get the same residuals via:

res <- Dataone$Var1 - predict(ModLMER, re.form = ~(1|Var3))

i.e., ignore the group-specific observation-level random effect term. However, it looks like there is a bug at the moment ("contrasts can be applied only to factors with 2 or more levels").

An extremely hacky solution is to construct the random-effect predictions without the observation-level term yourself:

## fixed-effect predictions
p0 <- predict(ModLMER, re.form = NA)
## construct RE prediction, Var3 term only:
Z <- getME(ModLMER, "Z")
b <- drop(getME(ModLMER, "b"))
## zero out observation-level components
b[101:200] <- 0
## add RE predictions to fixed predictions
p1 <- drop(p0 + Z %*% b)
## plot fitted vs residual
plot(p1, Dataone$Var1 - p1)

For what it's worth, this also works:

library(glmmTMB)
ModGLMMTMB <- glmmTMB(Var1~I(Var2)+I(Var2^2)+(1|Var3),
                      dispformula = ~factor(Var4),
                      REML = TRUE,
                      data = Dataone)

---

Related Topics