lme4中的固定差异值

Question

我正在使用函数使用lme4 R包创建线性混合模型。在这个模型中，我有四个随机效果和一个固定效果（截距）。我的问题是关于随机效应的估计方差。是否有可能以类似的方式为SAS和PARMS参数指定协方差参数的初始值。

在以下例子中，所估计的方差是：

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 

Answer 1

这是可能的，如果有点hacky。这里有一个重复的例子：

配合原有的模式：

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-标准）功能：

dd <- as.function(m1) 

我要设定标准偏差归零，以便我可以比较一下，即正则线性模型的系数。 （dd的参数向量是一个向量，包含模型中随机效应项的列方向，下三角，级联Cholesky因子。幸运的是，如果您拥有的只有标量/截距随机效应（例如， (1|x)），那么它们对应于由模型标准偏差缩放的随机效应标准偏差）。

(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对象，你可以如下做到这一点...

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））。这将是一个有点棘手，因为残留的标准偏差缩放的...

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 

所以这并不让我们相当，我们希望（因为剩余方差的变化......）

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）水平：看here ...