r - Fixing variance values in lme4

Question

I am using the lme4 R package to create a linear mixed model using the lmer() function. In this model I have four random effects and one fixed effect (intercept). My question is about the estimated variances of the random effects. Is it possible to specify initial values for the covariance parameters in a similar way as it can be done in SAS with the PARMS argument.

In the following example, the estimated variances are:

c(0.00000, 0.03716, 0.00000, 0.02306)

I would like to fix these to (for example)

c(0.09902947, 0.02460464, 0.05848691, 0.06093686)

so there are not estimated.

    > summary(mod1)
    Linear mixed model fit by maximum likelihood  ['lmerMod']
    Formula: log_cumcover_mod ~ (1 | kildestationsnavn) + (1 | year) + (1 |  
        kildestationsnavn:year) + (1 | proevetager)
       Data: res

         AIC      BIC   logLik deviance df.resid 
       109.9    122.9    -48.9     97.9       59 

    Scaled residuals: 
        Min      1Q  Median      3Q     Max 
    -2.1056 -0.6831  0.2094  0.8204  1.7574 

    Random effects:
     Groups                 Name        Variance Std.Dev.
     kildestationsnavn:year (Intercept) 0.00000  0.0000  
     kildestationsnavn      (Intercept) 0.03716  0.1928  
     proevetager            (Intercept) 0.00000  0.0000  
     year                   (Intercept) 0.02306  0.1518  
     Residual                           0.23975  0.4896  
    Number of obs: 65, groups:  
    kildestationsnavn:year, 6; kildestationsnavn, 3; proevetager, 2; year, 2

    Fixed effects:
                Estimate Std. Error t value
    (Intercept)   4.9379     0.1672   29.54

Answer 1

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 ...