1

pbkrtest パッケージの関数 PBmodcomp を使用して、2 つの 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)

この問題を解決するにはどうすればよいですか?

ありがとう

4

1 に答える 1

1

ローランドが提案したこの問題の解決

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)
于 2014-08-04T22:58:33.690 に答える