R: получить коэффициенты и CI из результатов начальной загрузки модели со смешанными эффектами.

Рабочие данные выглядят так:

set.seed(1234)
df <- data.frame(y = rnorm(1:30), 
                 fac1 = as.factor(sample(c("A","B","C","D","E"),30, replace = T)),
                 fac2 = as.factor(sample(c("NY","NC","CA"),30,replace = T)),
                 x = rnorm(1:30))

Модель lme устанавливается как:

library(lme4)
mixed <- lmer(y ~ x + (1|fac1) + (1|fac2), data = df)

Я использовал bootMer для запуска параметрической начальной загрузки, и я могу успешно получить коэффициенты (перехват) и SE для фиксированных и случайных эффектов:

mixed_boot_sum <- function(data){s <- sigma(data)
c(beta = getME(data, "fixef"), theta = getME(data, "theta"), sigma = s)}

mixed_boot <- bootMer(mixed, FUN = mixed_boot_sum, nsim = 100, type = "parametric", use.u = FALSE)

Мой первый вопрос: как получить коэффициенты (наклон) каждого отдельного уровня двух случайных эффектов из результатов начальной загрузки mixed_boot?

У меня нет проблем с извлечением коэффициентов (наклона) из модели mixed с помощью функции augment из пакета broom, см. ниже:

library(broom)
mixed.coef <- augment(mixed, df)

Однако похоже, что broom не может работать с объектом класса boot. Я не могу использовать вышеуказанные функции непосредственно на mixed_boot.

Я также пытался изменить mixed_boot_sum, добавив mmList (я думал, что это то, что я ищу), но R жалуется:

Error in bootMer(mixed, FUN = mixed_boot_sum, nsim = 100, type = "parametric",  : 
  bootMer currently only handles functions that return numeric vectors

Кроме того, можно ли получить CI как для фиксированных, так и для случайных эффектов, указав также FUN?

Теперь я очень запутался в правильных спецификациях для FUN, чтобы удовлетворить мои потребности. Любая помощь в отношении моего вопроса будет принята с благодарностью!


r lme4 statistics-bootstrap coefficients
person Chuan    schedule 06.09.2016    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Мой первый вопрос: как получить коэффициенты (наклон) каждого отдельного уровня двух случайных эффектов из результатов начальной загрузки Mixed_boot?

Я не уверен, что вы подразумеваете под «коэффициентами (наклоном) каждого отдельного уровня». broom::augment(mixed, df) дает прогнозы (остатки и т. д.) для каждого наблюдения. Если вам нужны прогнозируемые коэффициенты на каждом уровне, я бы попробовал

mixed_boot_coefs <- function(fit){
   unlist(coef(fit))
}

что для оригинальной модели дает

mixed_boot_coefs(mixed)
## fac1.(Intercept)1 fac1.(Intercept)2 fac1.(Intercept)3 fac1.(Intercept)4 
##        -0.4973925        -0.1210432        -0.3260958         0.2645979 
## fac1.(Intercept)5           fac1.x1           fac1.x2           fac1.x3 
##        -0.6288728         0.2187408         0.2187408         0.2187408 
##           fac1.x4           fac1.x5 fac2.(Intercept)1 fac2.(Intercept)2 
##         0.2187408         0.2187408        -0.2617613        -0.2617613 
##  ...

Если вы хотите, чтобы результирующий объект был более четко назван, вы можете использовать:

flatten <- function(cc) setNames(unlist(cc),
                                outer(rownames(cc),colnames(cc),
                                      function(x,y) paste0(y,x)))

mixed_boot_coefs <- function(fit){
   unlist(lapply(coef(fit),flatten))
}

При прогоне через bootMer/confint/boot::boot.ci эти функции дадут доверительные интервалы для каждого из этих значений (обратите внимание, что все наклоны facW.xZ идентичны для разных групп, поскольку модель предполагает случайные изменения только в точке пересечения). Другими словами, любую информацию, которую вы знаете, как извлечь из подобранной модели (условные режимы/BLUP [ranef], предсказанные точки пересечения и наклоны для каждого уровня группирующей переменной [coef], оценки параметров [fixef, getME], случайные эффекты дисперсии [VarCorr], прогнозы при определенных условиях [predict]...) можно использовать в аргументе FUN bootMer, при условии, что его структуру можно преобразовать в простой числовой вектор.

person Ben Bolker    schedule 06.09.2016
comment
Большое спасибо! @Бен Болкер. Это хорошо работает. Еще один вопрос, можно ли восстановить исходное имя каждого уровня группирующей переменной при использовании coef ? Я вижу, в ваших результатах R выдает числовые значения уровней вместо исходных имен. - person Chuan; 07.09.2016