Calculating predicted means (or predicted probabilities) and SE after multiple imputation in R

Question

I want to calculate predicted values and standard errors, but I can't simply use predict(), as I’m using 15 multiply imputed datasets (Amelia package generated). I run regression models on each dataset. Afterwards, results are combined into a single set of model coefficients and standard errors using the Amelia function mi.meld(), which uses Rubin’s rule.

Sample data and code:

dd<-list()
for (i in 1:15){
dd[[i]] <- data.frame(
  Age=runif(50,20,90),
  Cat=factor(sample(0:4, 50, replace=T)),
  Outcome = sample(0:1, 50, replace=T)
)}

b.out<-NULL
se.out<-NULL
for(i in 1:15) {
  ols.out<-glm(Outcome~Age+factor(Cat), data=dd[[i]],family="binomial")
  b.out <- rbind(b.out, ols.out$coef)
  se.out <- rbind(se.out, coef(summary(ols.out))[,2])}
mod0 <- mi.meld(q = b.out, se = se.out)

> mod0
$q.mi
     (Intercept)         Age factor(Cat)1 factor(Cat)2 factor(Cat)3 factor(Cat)4
[1,]   0.0466825 -0.00577106    0.5291908  -0.09760264    0.4058684    0.3125109

$se.mi
     (Intercept)        Age factor(Cat)1 factor(Cat)2 factor(Cat)3 
factor(Cat)4
[1,]    1.863276 0.02596468     1.604759     1.398322     1.414589     
1.332743

Now comes the problematic part. I want to calculate predicted values (in this case, predicted probabilities) and standard errors for the following set of predictor values:

data.predict <- data.frame(Cat=as.factor(c(0:4)), Age=53.6)
print(data.predict)
  Cat  Age
1   0 53.6
2   1 53.6
3   2 53.6
4   3 53.6
5   4 53.6

If I had fitted a model on 1 dataset, I would simply do this:

prediction<- predict(mod1, data.predict, type="response",se.fit=T)

However, I don't have a model object, I just have coefficients stored.. Now I have looked into two ways of solving this, the first being manually calculating the predictions in this way: predict() with arbitrary coefficients in r However I don't know how to get to the standard errors.. The other idea I had was to create a fake model object, like this function creates: https://gist.github.com/MrFlick/ae299d8f3760f02de6bf and to use it in predict(), but as the standard errors of the model aren't used, there is also no way of calculating the standard errors of the prediction..

Does anyone have a suggestion on how to solve this? I have tried to explain my problem clearly and with sample code, but if my question somehow isn't clear, please let me know so I can provide additional information. Thank you for your help!

Answer 1

I'm sure you don't need this answer after nearly years, but I've just bee working on a similar problem and I thought I would put the answer here for posterity.

Andrew Heiss put this solution on gisthub here - https://gist.github.com/andrewheiss/a3134085e92c6607db39c5b14e1b879e

I've modified it a tiny bit (partly because I think the default behaviour of 'nest' might have been altered in tidyverse since he wrote this ?)

The code (the hard work!) is almost entirely from Andrew Heiss. The annotation here is a mix of mine and his.

This is using the Amelia Africa dataset, for my real life problem I had a different dataset (obviously) and did the first few steps a bit differently, hwich was all fine.

library(tidyverse)
library(Amelia)
library(broom)

# Use the africa dataset from Amelia
data(africa)
set.seed(1234)
imp_amelia <- amelia(x = africa, m = 5, cs = "country", ts = "year", logs = "gdp_pc", p2s = 0) # do the imputations -- for me, it was fine to do this bit in 'mice'

# Gather all the imputed datasets into one data frame and run a model on each
models_imputed_df <- bind_rows(unclass(imp_amelia$imputations), .id = "m") %>%
  group_by(m) %>%
  nest() %>% 
  mutate(model = data %>% map(~ lm(gdp_pc ~ trade + civlib, data = .)))

# again - for my real life problem the models looked very different to this, and used rms - and this was also totally fine.

models_imputed_df
#> # A tibble: 5 x 3
#>   m     data               model   
#>   <chr> <list>             <list>  
#> 1 imp1  <tibble [120 × 7]> <S3: lm>
#> 2 imp2  <tibble [120 × 7]> <S3: lm>
#> 3 imp3  <tibble [120 × 7]> <S3: lm>
#> 4 imp4  <tibble [120 × 7]> <S3: lm>
#> 5 imp5  <tibble [120 × 7]> <S3: lm>


# We want to see how GDP per capita varies with changes in civil liberties, so
# we create a new data frame with values for each of the covariates in the
# model. We include the full range of civil liberties (from 0 to 1) and the mean
# of trade.

# ie. this is a 'skelton' data frame of all your variables that you want to make predictions over.

new_data <- data_frame(civlib = seq(0, 1, 0.1), 
                       trade = mean(africa$trade, na.rm = TRUE))
new_data
#> # A tibble: 11 x 2
#>    civlib trade
#>     <dbl> <dbl>
#>  1  0.     62.6
#>  2  0.100  62.6
#>  3  0.200  62.6
#>  4  0.300  62.6
#>  5  0.400  62.6
#>  6  0.500  62.6
#>  7  0.600  62.6
#>  8  0.700  62.6
#>  9  0.800  62.6
#> 10  0.900  62.6
#> 11  1.00   62.6

# write a function to meld predictions

meld_predictions <- function(x) {
  # x is a data frame with m rows and two columns:
  #
  # m  .fitted  .se.fit
  # 1  1.05     0.34
  # 2  1.09     0.28
  # x  ...      ...

  # Meld the fitted values using Rubin's rules
  x_melded <- mi.meld(matrix(x$.fitted), matrix(x$.se.fit))

  data_frame(.fitted = as.numeric(x_melded$q.mi),
             .se.fit = as.numeric(x_melded$se.mi))
}

# We augment/predict using new_data in each of the imputed models, then we group
# by each of the values of civil liberties (so each value, like 0.1 and 0.2 has
# 5 values, 1 from each of the imputed models), and then we meld those 5
# predicted values into a single value with meld_predictions()

predict_melded <- data_frame(models = models_imputed_df$model) %>%
  mutate(m = 1:n(),
         fitted = models %>% map(~ augment(., newdata = new_data))) %>% 
  unnest(fitted) %>% 
  dplyr::select(-models) %>% #### I needed to add this row to make the code work, once you've used the models to get the fit you don't need them in the data object anymore.  I took this line out because it was slowing everything down, then realised the code only works with this line... not sure why?
  group_by(civlib) %>%  
  nest(data=c(m, .fitted, .se.fit)) %>%  # needed to change this here from gisthub to make the nested 'data' have all the imputations in it, not just estimates from one of the imputations.
  mutate(fitted_melded = data %>% map(~ meld_predictions(.))) %>% 
  unnest(fitted_melded) %>% 
  mutate(ymin = .fitted + (qnorm(0.025) * .se.fit),
         ymax = .fitted + (qnorm(0.975) * .se.fit))


## NB. this is still on the link scale -- you'd need to write an extra few lines to exponentiate everything and get your predictions and se's on the response scale
# Plot!
ggplot(predict_melded, aes(x = civlib, y = .fitted)) +
  geom_line(color = "blue") +
  geom_ribbon(aes(ymin = ymin, ymax = ymax), alpha = 0.2, fill = "blue")

Hope that's at least a little bit helpful to anyone else stumbling across this.