가치가있는 것을 위해, 나는 비교적 간단한 LMM에 대한 분산 추정의 안정성을 살펴보기 위해 약간의 시뮬레이션 연구를 수행했습니다 (를 sleepstudy통해 사용 가능한 데이터 세트 사용 lme4). 첫 번째 방법은 ngroups여러 주제에 대해 가능한 모든 주제 조합을 생성하고 가능한 각 조합에 대해 모델을 다시 맞추는 것입니다. 두 번째는 주제의 여러 무작위 하위 집합을 취합니다.

library(lme4)
library(ggplot2)
library(tidyr)

m0 <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy,
           control = lmerControl(optimizer = "nloptwrap"))
# set the number of factor levels
ngroups <- 3:18
# generate all possible combinations
combos <- lapply(X = ngroups,
                 FUN = function(x) combn(unique(sleepstudy$Subject), x))

# allocate output (sorry, this code is entirely un-optimized)
out <- list(matrix(NA, ncol(combos[[1]]), 1), matrix(NA, ncol(combos[[2]]), 1),
            matrix(NA, ncol(combos[[3]]), 1), matrix(NA, ncol(combos[[4]]), 1),
            matrix(NA, ncol(combos[[5]]), 1), matrix(NA, ncol(combos[[6]]), 1),
            matrix(NA, ncol(combos[[7]]), 1), matrix(NA, ncol(combos[[8]]), 1),
            matrix(NA, ncol(combos[[9]]), 1), matrix(NA, ncol(combos[[10]]), 1),
            matrix(NA, ncol(combos[[11]]), 1), matrix(NA, ncol(combos[[12]]), 1),
            matrix(NA, ncol(combos[[13]]), 1), matrix(NA, ncol(combos[[14]]), 1),
            matrix(NA, ncol(combos[[15]]), 1), matrix(NA, ncol(combos[[16]]), 1))
# took ~ 2.5 hrs on my laptop, commented out for safety
#system.time(for(ii in 1:length(combos)) {
#    for(jj in 1:ncol(combos[[ii]])) {
#    sls <- sleepstudy[sleepstudy$Subject %in% combos[[ii]][,jj],]
#    out[[ii]][jj] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
#        }
#    })

# pad with zeros, not all were equal
# from http://stackoverflow.com/questions/11148429/r-convert-asymmetric-list-to-matrix-number-of-elements-in-each-sub-list-diffe
max.len <- max(sapply(out, length))
corrected.list <- lapply(out, function(x) {c(x, rep(NA, max.len - length(x)))})
mat <- do.call(rbind, corrected.list)
mat <- data.frame(t(mat))
names(mat) <- paste0('s',3:18)
mat <- gather(mat, run, value)

ggplot(mat, aes(x = value, fill = run)) +
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash',
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) +
    theme_bw() +
    guides(fill = FALSE)

점선으로 표시된 검은 선은 분산의 원래 점 추정치이며, 패싯은 서로 다른 수의 주제를 나타냅니다 ( s3세 명의 대상 그룹, s4네 등).

그리고 다른 방법 :

ngroups <- 3:18
reps <- 500
out2<- matrix(NA, length(ngroups), reps)

for (ii in 1:length(ngroups)) {
    for(j in 1:reps) {
        sls <- sleepstudy[sleepstudy$Subject %in% sample(unique(sleepstudy$Subject), ngroups[i], replace = FALSE),]
        out2[i,j] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
    }
}
out2 <- data.frame(t(out2))
names(out2) <- paste0('s',3:18)
out2 <- gather(out2, run, value)

ggplot(out2, aes(x = value, fill = run)) +
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash',
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) +
    theme_bw() +
    guides(fill = FALSE)

(이 예제의 경우, 어쨌든) 차이는 나중에 14 명 이상이 아닌 경우 적어도 14 명의 피험자가있을 때까지 실제로 안정화되지 않는 것으로 보입니다.

답변