Блог

Главная - Вопросы по программированию - Извлеките доверительный интервал для комбинации (суммы) оценок дисперсии из смешанной модели (lme4)

Извлеките доверительный интервал для комбинации (суммы) оценок дисперсии из смешанной модели (lme4)

#r #lme4 #mixed-models #random-effects

Вопрос:

Я использую смешанную модель, чтобы разделить дисперсию некоторых измерений на дисперсию между группами и внутри групп. Модель выглядит так:

 library(lme4)

k <- 100 
i <- 5

group_mean <- rnorm(k, 0, 1)
y <- rnorm(k*i, rep(group_mean, each = i), 0.5)

df <- data.frame(y, group = rep(1:k, each = i))
    
m <- lmer(y ~ (1|group), data = df)

Я могу извлечь между — и внутри группы SD и доверительные интервалы:

 VarCorr(m)
#>  Groups   Name        Std.Dev.
#>  group    (Intercept) 0.99958 
#>  Residual             0.50327
confint(m, oldNames = FALSE)
#> Computing profile confidence intervals ...
#>                           2.5 %     97.5 %
#> sd_(Intercept)|group  0.8640826 1.15724604
#> sigma                 0.4703110 0.54026452
#> (Intercept)          -0.3782316 0.02526364

Как я могу получить доверительные интервалы для суммы дисперсий ( sd_(Intercept)|group^2 sigma^2 )? Т. Е. оцененную дисперсию только y с одним наблюдением на каждое group .

Комментарии:

1. Вы могли бы его загрузить?

2. Это может быть решением проблемы. Я вижу, что у confint() есть возможность предоставить пользовательскую функцию загрузки.

3. Я бы изучил возможность использования bootMer .

4. Прочитав немного документации, я вижу , что confint() можно передать модель и пользовательскую функцию bootMer() , а затем рассчитать CI по полученным оценкам начальной загрузки ( confint(m, method = 'boot', FUN = get_total_sd, nsim = 10000) ). Работает довольно хорошо, хотя для 10000 повторных образцов этого простого набора данных требуется несколько минут. Хотя это не лишено смысла для моего использования.

5. Вам рекомендуется самостоятельно ответить на свой вопрос, основываясь на этом.

Ответ №1:

Как прокомментировал @Roland, это можно загрузить. lme4 обеспечивает хорошую основу для этого:

 library(lme4)
#> Loading required package: Matrix

k <- 100 
i <- 5

group_mean <- rnorm(k, 0, 1)
y <- rnorm(k*i, rep(group_mean, each = i), 0.5)

df <- data.frame(y, group = rep(1:k, each = i))

m <- lmer(y ~ (1|group), data = df)

get_total_sd <- function(mod) {
    varCorr_df <- as.data.frame(VarCorr(mod))
    sd_components <- varCorr_df[['sdcor']]
    names(sd_components) <- varCorr_df[['grp']]

    sd_total <- sqrt(sum(sd_components^2))

    c(sd_components, "Total" = sd_total)
}

confint(m, method = 'boot', FUN = get_total_sd, nsim = 10000)
#>              2.5 %    97.5 %
#> group    0.7655064 1.0329894
#> Residual 0.4921279 0.5657032
#> Total    0.9290245 1.1600250

confint(m, oldNames = FALSE)
#> Computing profile confidence intervals ...
#>                            2.5 %    97.5 %
#> sd_(Intercept)|group  0.77542074 1.0439971
#> sigma                 0.49426283 0.5677789
#> (Intercept)          -0.01914245 0.3472213

Создано 2021-06-03 пакетом reprex (v2.0.0)

Загруженный (процентиль) CIS для отдельных SDs аналогичны профилированным, но не равны. Я не уверен, что лучше. Это займет несколько минут, поэтому меня все еще интересуют другие подходы.

Метки: