在LME4中修复方差值 [英] Fixing variance values in lme4

问题描述

我正在使用lme4 R包使用lmer()函数创建线性混合模型.在此模型中，我有四个随机效应和一个固定效应(拦截).我的问题是关于随机效应的估计方差.是否可以像在SAS中使用PARMS参数一样，为协方差参数指定初始值.

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)

我想修复这些(例如)

c(0.09902947, 0.02460464, 0.05848691, 0.06093686)

所以没有估算值.

    > 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

推荐答案

这可能，即使有点hacky.这是一个可重现的示例:

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

适合原始模型:

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 

恢复偏差(在这种情况下为REML标准)功能:

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

dd <- as.function(m1)

我将标准偏差设置为零，以便可以比较一些东西，即常规线性模型的系数. (dd的参数向量是一个向量，其中包含模型中缩放的随机效应项的按列，下三角，级联的Cholesky因子.幸运的是，如果您拥有的都是标量/仅截取随机效应(例如(1|x))，则它们对应于随机效应标准差，由模型标准差缩放).

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

比赛:

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

如果要构造一个新的merMod对象，可以按以下步骤进行操作...

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

现在假设我们要将方差​​设置为特定的非零值(例如(700,30)).这会有点棘手，因为要按残留标准偏差进行缩放...

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

顺便说一句，通常不建议在分组变量的数量非常少(例如< 5或6)级别时使用随机效果:请参见

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

这篇关于在LME4中修复方差值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！