Как получить среднее значение по группе после множественного вменения?
Предположим, мы используем следующий пример из пакета mice:
imp <- mice(nhanes, m=5, print=FALSE, seed=55152)
Я хочу получить средства nhanes$chl, сгруппированные по nhanes$hyp. Без использования множественного вменения можно было бы сделать:
library(mice)
library(dplyr)
nhanes %>%
group_by(hyp) %>%
summarise(chl=mean(chl, na.rm=TRUE))
Как я могу получить такие объединенные оценки средних значений chl на группу hyp после многократного вменения с помощью mice?
Ответы 1
Без вменения используются только строки с полными вариантами для chl и hyp, что называется удалением по списку. Используя этот метод, мы получаем эти mean по группам:
> tapply(X=nhanes$chl, INDEX=nhanes$hyp, FUN=mean, na.rm=TRUE)
1 2
181.8182 229.0000
Множественное вменение,
> m. <- 5
> imp <- mice::mice(nhanes, m=m., print=FALSE, seed=55152)
дает нам m вмененных наборов данных, которые мы обычно объединяем в длинном формате. (Обратите внимание, что рекомендуется m > 20.1)
> long <- mice::complete(imp, action='long')
> with(long, table(.imp)) ## nrows by .imp
.imp
1 2 3 4 5
25 25 25 25 25
> nrow(nhanes) ## nrows original
[1] 25
Каждый вмененный набор данных будет иметь свою собственную статистику (т. е. в данном случае группу mean). Например, .imp == 1 будет иметь:
> tapply(long[long$.imp == 1, ]$chl, long[long$.imp == 1, ]$hyp, mean, na.rm=TRUE)
1 2
185.8421 212.5000
Согласно правилам Рубина, для получения оценки нам нужно среднее значение этой вмененной статистики. Соответственно, мы рассчитываем статистику для каждого вменения,
> (r0 <- do.call('rbind', by(long, ~ .imp, \(x) tapply(x$chl, x$hyp, mean))))
1 2
1 185.8421 212.5000
2 186.9000 231.0000
3 187.9000 237.8000
4 189.0000 232.5000
5 188.1053 215.8333
и возьмем средства этих результатов.
> colMeans(r0)
1 2
187.5495 225.9267
Мы могли бы сообщать оценки со стандартными ошибками, которые мы вычисляем соответствующим образом:
> ## compute mean and variance for each imputation
> R1 <- by(long, ~ .imp, \(x)
+ sapply(c(mean=mean, var=\(x) var(x)/length(x)), \(f)
+ tapply(x$chl, x$hyp, f))) |>
+ simplify2array()
>
> Q <- rowMeans(R1[, 1, ]) ## calculate mean estimates
> U <- rowMeans(R1[, 2, ]) ## calculate within variances
>
> B <- rowSums(((R1[, 1, ] - Q)^2))/(m. - 1) ## calculate between variances
>
> T <- U + (1 + 1/m.)*B ## calculate total variances
>
> cbind(Estimate=Q, 'Std. Error'=sqrt(T))
Estimate Std. Error
1 187.5495 9.175414
2 225.9267 21.404731
Данные:
> nhanes <- mice::nhanes
Обратите внимание, что я исправил расчет стандартных ошибок, в котором неправильно использовалась дисперсия вместо дисперсии выборочного среднего.