Могу ли я объединить вмененные оценки модели случайных эффектов, используя пакет mi?
Похоже, что mi В последние несколько лет пакет несколько переписан.
"Старый" способ работы хорошо описан в следующем руководстве: http://thomasleeper.com/Rcourse/Tutorials/mi.html
"Новый" способ сделать что-то (придерживаясь демонстрационной модели Leeper) выглядит примерно так:
#load mi
library(mi)
#set seed
set.seed(10)
#simulate some data (with some observations missing)
x1 <- runif(100, 0, 5)
x2 <- rnorm(100)
y <- 2*x1 + 20*x2 + rnorm(100)
mydf <- cbind.data.frame(x1, x2, y)
mydf$x1[sample(1:nrow(mydf), 20, FALSE)] <- NA
mydf$x2[sample(1:nrow(mydf), 10, FALSE)] <- NA
# Convert to a missing_data.frame
mydf_mdf <- missing_data.frame(mydf)
# impute
mydf_imp <- mi(mydf_mdf)
Хотя имена функций изменились, это на самом деле очень похоже на "старый" способ ведения дел.
Самым большим изменением (с моей точки зрения) является замена следующих "старых" функций
lm.mi(formula, mi.object, ...)
glm.mi(formula, mi.object, family = gaussian, ...)
bayesglm.mi(formula, mi.object, family = gaussian, ...)
polr.mi(formula, mi.object, ...)
bayespolr.mi(formula, mi.object, ...)
lmer.mi(formula, mi.object, rescale=FALSE, ...)
glmer.mi(formula, mi.object, family = gaussian, rescale=FALSE, ...),
Ранее пользователь мог вычислить модель для каждого вмененного набора данных, используя одну из этих функций, а затем объединить результаты, используя mi.pooled() (или же coef.mi() если мы будем следовать примеру Липера).
В текущей версии mi (У меня установлен v1.0), эти последние шаги, кажется, были объединены в одну функцию, pool(), pool() Похоже, что функция считывает семейство и функцию связи, которые были назначены переменной в процессе импутации выше, а затем оценивает модель с bayesglm используя указанную формулу, как показано ниже.
# run models on imputed data and pool the results
summary(pool(y ~ x1 + x2, mydf_imp))
##
## Call:
## pool(formula = y ~ x1 + x2, data = mydf_imp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.98754 -0.40923 0.03393 0.46734 2.13848
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.34711 0.25979 -1.336 0.215
## x1 2.07806 0.08738 23.783 1.46e-13 ***
## x2 19.90544 0.11068 179.844 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.7896688)
##
## Null deviance: 38594.916 on 99 degrees of freedom
## Residual deviance: 76.598 on 97 degrees of freedom
## AIC: 264.74
##
## Number of Fisher Scoring iterations: 7
Похоже, мы приближаемся к восстановлению наших смоделированных бета-значений (2 и 20). Другими словами, он ведет себя как ожидалось.
Давайте возьмем немного больший набор данных с наивно смоделированным случайным эффектом только для получения группирующей переменной.
mydf2 <- data.frame(x1 = rep(runif(100, 0, 5), 20)
,x2 = rep(rnorm(100, 0, 2.5), 20)
,group_var = rep(1:20, each = 100)
,noise = rep(rnorm(100), 20))
mydf2$y <- 2*mydf2$x1 + 20*mydf2$x2 + mydf2$noise
mydf2$x1[sample(1:nrow(mydf2), 200, FALSE)] <- NA
mydf2$x2[sample(1:nrow(mydf2), 100, FALSE)] <- NA
# Convert to a missing_data.frame
mydf2_mdf <- missing_data.frame(mydf2)
show(mydf2_mdf)
## Object of class missing_data.frame with 2000 observations on 5 variables
##
## There are 4 missing data patterns
##
## Append '@patterns' to this missing_data.frame to access the corresponding pattern for every observation or perhaps use table()
##
## type missing method model
## x1 continuous 200 ppd linear
## x2 continuous 100 ppd linear
## group_var continuous 0 <NA> <NA>
## noise continuous 0 <NA> <NA>
## y continuous 0 <NA> <NA>
##
## family link transformation
## x1 gaussian identity standardize
## x2 gaussian identity standardize
## group_var <NA> <NA> standardize
## noise <NA> <NA> standardize
## y <NA> <NA> standardize
поскольку missing_data.frame() кажется, интерпретирует group_var как непрерывный, я использую change() функция от mi переназначить "un" для "неупорядоченного категориального", а затем действуйте, как указано выше.
mydf2_mdf <- change(mydf2_mdf, y = "group_var", what = "type", to = "un" )
# impute
mydf2_imp <- mi(mydf2_mdf)
Теперь, если версия 1.0 mi удалил функциональность предыдущих версий (т.е. функциональность, доступную с lmer.mi а также glmer.mi), Я бы предположил, что добавление случайного эффекта в формулу должно указывать pool() к соответствующему lme4 функция. Тем не менее, первоначальное сообщение об ошибке предполагает, что это не так.
# run models on imputed data and pool the results
summary(pool(y ~ x1 + x2 + (1|group_var), mydf2_imp))
## Warning in Ops.factor(1, group_var): '|' not meaningful for factors
## Warning in Ops.factor(1, group_var): '|' not meaningful for factors
## Error in if (prior.scale[j] < min.prior.scale) {: missing value where TRUE/FALSE needed
После моего предупреждающего сообщения и извлечения целых чисел из моего фактора получается оценка, но результаты показывают, что pool() все еще оценивает модель с фиксированным эффектом bayesglm и удерживая мою попытку постоянного случайного эффекта.
summary(pool(y ~ x1 + x2 + (1|as.numeric(as.character(group_var))), mydf2_imp))
##
## Call:
## pool(formula = y ~ x1 + x2 + (1 | as.numeric(as.character(group_var))),
## data = mydf2_imp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.93633 -0.69923 0.01073 0.56752 2.12167
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.383e-01 2.596e+02 0.001
## x1 1.995e+00 1.463e-02 136.288
## x2 2.000e+01 8.004e-03 2499.077
## 1 | as.numeric(as.character(group_var))TRUE -3.105e-08 2.596e+02 0.000
## Pr(>|t|)
## (Intercept) 1
## x1 <2e-16 ***
## x2 <2e-16 ***
## 1 | as.numeric(as.character(group_var))TRUE 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.8586836)
##
## Null deviance: 5384205.2 on 1999 degrees of freedom
## Residual deviance: 1713.9 on 1996 degrees of freedom
## AIC: 5377
##
## Number of Fisher Scoring iterations: 4
Мои вопросы:
- Можно ли легко сгруппировать оценки случайных эффектов, используя
mi?, а также - Если да, то как?
2 ответа
Вы можете указать FUN аргумент pool() Функция для изменения оценки. В вашем случае это было бы summary(pool(y ~ x1 + x2 + (1|as.numeric(as.character(group_var))), data = mydf2_imp, FUN = lmer)), Это может или не может на самом деле работать, но это законный синтаксис. Если это не удается, то вы можете использовать complete функция для создания завершенных data.frames, вызов lmer на каждом и усредните результаты самостоятельно, что было бы что-то вроде
dfs <- complete(mydf2_imp)
estimates <- lapply(dfs, FUN = lme4, formula = y ~ x1 + x2 +
(1|as.numeric(as.character(group_var))))
rowMeans(sapply(estimates, FUN = fixef))
Просто, чтобы предоставить альтернативу, есть пакет, который фокусируется на MI для моделей со смешанными эффектами, а также объединяет полученные результаты (mitml найди это здесь).
Использовать пакет довольно просто. Полагается на пакеты pan а также jomo для вменения, но он также может обрабатывать ввод из других пакетов MI (?as.mitml.list).
Пул оценки из модели смешанных эффектов в основном автоматизированы и включены в testEstimates функция.
require(mitml)
require(lme4)
data(studentratings)
# impute example data using 'pan'
fml <- ReadDis + SES ~ ReadAchiev + (1|ID)
imp <- panImpute(studentratings, formula=fml, n.burn=1000, n.iter=100, m=5)
implist <- mitmlComplete(imp, print=1:5)
# fit model using lme4
fit.lmer <- with(implist, lmer(SES ~ (1|ID)))
# pool results using 'Rubin's rules'
testEstimates(fit.lmer, var.comp=TRUE)
Выход:
# Call:
# testEstimates(model = fit.lmer, var.comp = TRUE)
# Final parameter estimates and inferences obtained from 5 imputed data sets.
# Estimate Std.Error t.value df p.value RIV FMI
# (Intercept) 46.988 1.119 41.997 801.800 0.000 0.076 0.073
# Estimate
# Intercept~~Intercept|ID 38.272
# Residual~~Residual 298.446
# ICC|ID 0.114
# Unadjusted hypothesis test as appropriate in larger samples.