Philippe Massicotte - 2 months ago 16
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 ...

Source (Stackoverflow)