Multiple Imputation with MICE in R

This is my code from the Feb 27 meeting on how to deal with missing data. I wouldn't say that all the methods I used here are foolproof, but feel free to take a look if you're wanting to give MICE a try! Happy to answer any questions.

Code:

install.packages("foreign")
install.packages("mice")
install.packages("VIM")
library("foreign")
library("mice")
library("VIM")
data<-read.spss(file = "/Users/kirsten/Box Sync/Methods/Missing Data Presentation/PALS missing data - vars of interest.sav",
                  to.data.frame=TRUE)
str(data) #Checking variable types, formats, values
md.pattern(data) #mice command for calling missing data pattern - this is a little overwhelming with this many vars/cases
#this next command calls a plot comparing missing data on two variables
marginplot(data[, c("alc015_age21", "alc015_age22")], col = mdc(1:2), cex = 1.2, cex.lab = 1.2, cex.numbers = 1.3, pch = 19)

#Testing imputation with default prediction methods (pmm, logreg, polyreg)
imp <- mice(data, seed = 23109)
#Looking at results
print(imp)
#Checking for implausible values - not really sure if there are possible implausible values here...
imp$imp$marstatDbinary.21
imp$imp$alc015_age22 #This has labels, not numerical values...wonder how to change that setting
#calling the first imputation - this pulls the original values with missing data filled in 
#by the first set of imputed values
complete(imp)
#calling the second imputation
complete(imp, 2)
#Inspecting distributions of original data and imputed data
require(lattice)
stripplot(imp, pch = 20, cex = 1.2)
#looking at scatterplot of original and imputed data
xyplot(imp, alc015_age21 ~ alc015_age22 | .imp, pch = 20, cex = 1.4)
#checking convergence
plot(imp, c("highed", "ptsuTOTX7_age21", "ptsuTOTX7_age22"))
#looking at kernal density plots 
densityplot(imp, scales = list(x = list(relation = "free")), layout = c(5, 1))

###SETTING UP ACTUAL IMPUTATION###

#Dropping a few unnecessary variables
str(data)
data1 <- data[c(-18,-19)]
str(data1)
#Recoding marital status to factor
data1$marstatDbinary.21 <- factor(data1$marstatDbinary.21)
data1$marstatDbinary.22 <- factor(data1$marstatDbinary.22)
#Recoding ethnicity to white/nonwhite
summary(data1$ethnic)
data1$ethnic <- as.numeric(data$ethnic)
data1$ethnic <- ifelse(data1$ethnic > 1, c("2"), c("1")) 
#Recoding ordinal variables to continuous
data1$ethnic <- as.numeric(data$ethnic)
data1$alc015_age21 <-as.numeric(data$alc015_age21)
data1$alc015_age22 <- as.numeric(data$alc015_age22)
data1$highed <- as.numeric(data$highed)
str(data1)
#Removing idnum from predictor set
imp <- mice(data1, print = FALSE)
imp$predictorMatrix
pred1 <- imp$predictorMatrix
pred1[,"idnum"] <- 0
pred1
#Test imputation with default settings, but incorporating specified predictor set
imp2 <- mice(data1, pred = pred1, seed = 23109)
#Looking at results
print(imp2)
#Checks out
#Running full imputation with 40 imputations with default settings
imp3 <- mice(data1, m = 40, pred = pred1, seed = 23109)


#looking at distributions of original and imputed data
print(imp3)
require(lattice)
stripplot(imp3, pch = 20, cex = 1.2)
#checking convergence
plot(imp3, c("highed", "ptsuTOTX7_age21", "ptsuTOTX7_age22"))
#Looking at kernel density plots of continuous variables
densityplot(imp3, scales = list(x = list(relation = "free")), layout = c(5, 1))

#Running analysis
fit <- with(imp3, lm(alc015_age22 ~ psuTOTX6_age21 + alc015_age21+ currented.22 
                     + marstatDbinary.22 + gender + ethnic + highed))

#pooling and printing results
print(pool(fit))
round(summary(pool(fit)), 2)

#Running the analysis with complete case data in original dataset
fit2 <- lm(alc015_age22 ~ psuTOTX6_age21 + alc015_age21+ currented.22 + marstatDbinary.22 + gender + ethnic + highed, data=data1)
summary(fit2)