Multiple imputation using R

R has a long list of packages for multiple imputation. The main problem is integration: statistical procedures in other packages may or may not work with the imputation procedures. I have been using Amelia together with Zelig. Because they were written by the same group, they work well together. However, I have been having trouble with making multiple imputation to work with the plm package. After searching the internet, here comes the solution:

1. Impute the missing data using Amelia or Mice.
2. Estimate the model on each imputed data.
3. Use the mitools package to extract and combine results. 

For example, here is a simple example:

...

imp <- mice(d)

mydata <- imputationList(lapply(1:5, complete, x = imp))

fit <- lapply(mydata$imputations, function(x){

plm(cog3pl ~ oc + grade9 + boy + han + ruralbirth, data = x,

index = c("schids"), model = "pooling")})

betas <- MIextract(fit, fun = coef)

vars <- MIextract(fit, fun = vcov)

summary(MIcombine(betas, vars))

I bet this will work for most, if not all, estimation procedures in R.

Prediction, missing data, etc. in Stan

library(rstan)

N <- 1001

N_miss <- ceiling(N / 10)

N_obs <- N - N_miss

mu <- 3

sigma <- 2

y_obs <- rnorm(N_obs, mu, sigma)

missing_data_code <-

'

data {

  int N_obs;

  int N_miss;

  real y_obs[N_obs];

}

parameters {

  real mu;

  real sigma;

  real y_miss[N_miss];

}

model {

  // add prior on mu and sigma here if you want

  y_obs ~ normal(mu,sigma);

  y_miss ~ normal(mu,sigma);

}

generated quantities {

  real y_diff;

  y_diff <- y_miss[101] - y_miss[1];

}

'

results <- stan(model_code = missing_data_code,

                data = list(N_obs = N_obs, N_miss = N_miss, y_obs = y_obs))

y_diff <- apply(extract(results, c("y_miss[1]", "y_miss[101]")), 1:2, diff)

Handling missing data in R

Some estimation procedures, including MCMCglmm, does not handle missing values directly. Before estimating these models, missing values need to be excluded. From this post, the command to do this is:

newdata <- na.omit(mydata)