在 lmer 对象上运行参数引导程序 (PBmodcomp) 时出错答案

【问题标题】：Error runnning parametric bootstrap (PBmodcomp) on lmer objects在 lmer 对象上运行参数引导程序 (PBmodcomp) 时出错
【发布时间】：2014-08-04 04:50:13
【问题描述】：

我正在尝试使用 pbkrtest 包中的函数 PBmodcomp 对两个 lmer 模型进行模型比较。但是我收到以下错误。

Error in `[<-.data.frame`(`*tmp*`, , rcol, value = c(0.318337014579985,  : 
  replacement has 4080 rows, data has 4458

我的数据在这里可用：https://www.dropbox.com/s/oweyw767qtpbqot/Data.txt

head(dat)
        Subject time        age  cognition gender
60002.1   60002    1  0.4898039 -0.6915897      2
60002.2   60002    2  4.4898039 -0.8999999      2
60002.3   60002    3  8.4898039 -1.1619855      2
60008.1   60008    1  2.4898039 -0.2106083      2
60008.2   60008    2  6.4898039  0.3355440      2
60008.3   60008    3 10.4898039 -0.7309111      2

我运行的代码是

library(lme4)
library(pbkrtest)
m2 <- lmer(cognition ~ age + (age | Subject), data = dat, REML = FALSE, na.action = na.omit, control = lmerControl(optimizer = "Nelder_Mead"))
m3 <- lmer(cognition ~ age + gender + (age | Subject), data = dat, REML = FALSE, na.action = na.omit, control = lmerControl(optimizer = "Nelder_Mead"))
pb <- PBmodcomp(m3, m2)

我该如何解决这个问题？

谢谢

【问题讨论】：

你有没有试过traceback()或debug()这个问题？
在拟合模型之前您是否尝试过移除 NA？
我在 lmer 模型中指定了 na.action = na.omit，所以我认为应该删除 NA。
是的，他们是。但是这可能导致m2 和m3 的数据不同。因此，我建议在创建模型之前尝试删除它们。
感谢 Roland，这解决了这个特定问题。

标签： r lme4 lmer statistics-bootstrap

【解决方案1】：

按照 Roland 的建议解决此问题

dat2 <-dat[complete.cases(dat[,3:4]),] #remove cases with missing values

m2 <- lmer(cognition ~ age + (age | Subject), data = dat2, REML = FALSE, na.action = na.omit, control = lmerControl(optimizer = "Nelder_Mead"))
m3 <- lmer(cognition ~ age + gender + (age | Subject), data = dat2, REML = FALSE, na.action = na.omit, control = lmerControl(optimizer = "Nelder_Mead"))

pb <- PBmodcomp(m3, m2, nsim = 10)

pb
Parametric bootstrap test; time: 3.93 sec; samples: 10 extremes: 0;
large : cognition ~ age + gender + (age | Subject)
small : cognition ~ age + (age | Subject)
         stat df   p.value    
LRT    15.629  1 7.706e-05 ***
PBtest 15.629      0.09091 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

编辑

由于我的实际数据正在检查多个响应变量，因此我在循环中运行了这些函数，因此我也粘贴了下面的循环代码。

#first loop with the reduced model
fitlmer.strings.01 <- function(X){
    s <- colnames(X)
    lapply(s, function(.s){
    y <- dat[, .s]
    index <- complete.cases(cbind(y, dat[, "age"]))
    dat <- data.frame(dat[index, ])
    y <- y[index]
    out <- with(dat, lmer(y ~ age + (age | Subject), REML = FALSE, control = lmerControl(optimizer = "Nelder_Mead")))    
    out  
    })
    } 
output.01 <- fitlmer.strings.01(dat[4:5])

#second loop with the full model 
fitlmer.strings.02 <- function(X){
    s <- colnames(X)
    lapply(s, function(.s){
    y <- dat[, .s]
    index <- complete.cases(cbind(y, dat[, "age"]))
    dat <- data.frame(dat[index, ])
    y <- y[index]
    out <- with(dat, lmer(y ~ age + gender + (age | Subject), REML = FALSE, control = lmerControl(optimizer = "Nelder_Mead")))    
    out  
    })
    } 

output.02 <- fitlmer.strings.02(dat[4:5])
pb <- PBmodcomp(output.02[[1]], output.01[[1]], nsim = 10)

【讨论】：