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