Philippe Massicotte - 1 year ago 81
R Question

# Fixing variance values in lme4

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

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

``````(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()))
``````

If you want to set the variances to (700,30), then follow the same procedure but start with `dd(sqrt(700),sqrt(30))`. Calling the deviance function sets the variances internally and re-estimates the fixed-effect parameters internally ...

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

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download