Объединенные оценки для преобразованных множественных вмененных наборов данных

Question

У меня есть продольный набор данных, на котором я выполнил многократное продольное вменение (см. Здесь метод).

Я хочу рассчитать общий балл для каждого из m вмененных наборов данных, сложив вместе баллы за каждый период времени для каждого человека. Наконец, мне нужно получить оценку, t-значение и p-значение для анализа на основе результатов объединения для m наборов данных.

Вот данные игрушки с добавленными недостающими данными.

set.seed(345)
A0 <- rnorm(4,2,.5)
B0 <- rnorm(4,2+3,.5)
A1 <- rnorm(4,6,.5)
B1 <- rnorm(4,6+2,.5)
A2 <- rnorm(4,10,.5)
B2 <- rnorm(4,10+1,.5)
A3 <- rnorm(4,14,.5)
B3 <- rnorm(4,14+0,.5)
score <- c(A0,B0,A1,B1,A2,B2,A3,B3)
id <- rep(1:8,times = 4, length = 32)
timeFac <- rep(c(paste0("time",0:3)), each = 8, length = 32)
group <- rep(c("A","B"), times =2, each = 4, length = 32)
df <- data.frame(id = id, group = group, time = timeFac,  score = score)

# plots
(ggplot(df, aes(x = timeFac, y = score, group = group)) + 
    stat_summary(fun.y = "mean", geom = "line", aes(linetype = group)) +
    stat_summary(fun.y = "mean", geom = "point", aes(shape = group), size = 3) +
    coord_cartesian(ylim = c(0,18)))

# now place some NAs
df[sample(1:nrow(df), 10, replace = F),"score"] <- NA

df

Теперь для продольного многоуровневого вменения

library(mice)

ini <- mice(df, maxit=0)
pred <- ini$predictorMatrix
pred["score",] <- c(-2, 1, 2, 0)
imp <- mice(df, 
             method = c("", "", "", "2l.pan"),
             pred = pred, 
             maxit = 1, 
             seed = 71152)

Теперь мы можем извлечь 5 полных наборов данных из imp объект и рассчитать общую оценку для каждого номера идентификатора.

createDF <- function (i) {
d <- dcast(complete(imp,i), id + group ~ timeFac, value.var = "score")
d <- mutate(d, tot = time1 + time2 + time3) %>% select(-one_of("time1", "time2", "time3"))
return(d)
}

Создайте объект со всеми пятью полными наборами данных

dfList <- lapply(1:imp$m, function (i) createDF(i))

Который мы можем затем индексировать для проведения анализа

summary(lm(tot ~ group + time0, data = dfList[[1]]) 

Мой вопрос заключается в том, как мне вручную или с помощью пакета объединить эти анализы для получения одного коэффициента, t-значения и p-значения на основе пяти вмененных наборов данных? Я знаю, как объединить оценки перед выполнением операций над вмененными наборами данных (используя, например, pool(with(imp, lm(y~x)))) но не после.

В принципе, кто-нибудь может предложить лучший способ сделать то, что я пытаюсь сделать?

Я читал о правилах Рубина для объединения оценок, но пока не нашел простого проработанного примера того, как применять их вручную (если нет пакета, который может делать то, что мне нужно).

ps это может быть больше вопросом CV, но так как я надеюсь на пакет, который может делать то, что мне нужно, я выбрал SO