- Impute the missing data using Amelia or Mice.
- Estimate the model on each imputed data.
- 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.