Is the following procedure to measure the quality of an imputation ok?

I'd like to compare different kinds of imputation techniques, i.e. methods which allow to fill missing data fields in a data frame.

For now, I'm only using the R package mice, which uses multiple imputation with chained equations. Below, I'm describing the steps I am following, and I'd appreciate comments if this is a reasonable approach.

1. The data I am currently using is an excerpt from the boys dataset (which is included in mice), where I make sure that the region variable reg is always there.

df <- boys[,c("age","hgt","wgt","bmi","hc","reg")]
df <- df[which(!(df[,factor.variable.name] %in% NA)),]

2. To be able to make the imputation results more comparable, I log-transform and normalize/centralize df column-wise.

3. An extra column contains a numeric value numeric.class as an indicator of the reg value is added.

4. I'm imputing the values of the 5 numeric variables "age", "hgt","wgt","bmi","hc" with a two-level linear model which uses numeric.class as class variable, which I do by

 predMatr[1:5,1:5] <- (matrix(data=1, 5, ncol=5) -diag(5))
predMatr[,"numeric.class"] <- -2
predMatr[,"reg"] <- 0
method <- c(rep("2l.norm",5),"", "")


and calling imp <- mice(df, meth=method, pred=predMatr, maxit=30, m=1).

5. To measure the performance of this, I actually replace the df in the mice call from the last step by a validation dataframe which has been preprocessed in the following way.

6. Setting the validation percentage to, say, 10%, and going column by column, I randomly select and replace 10% of all non-NA entries in a given column by NA. I picked this stratified approach since, usually, there are not complete rows missing, but rather only single (or possibly multiple) entries.

7. For the entries where I introduced extra NAs, I compute the RMSE column-wise, i.e. the mean of the square of the difference between the imputed values and the acual values.

8. I repeat steps 5-7, replacing in each column 10% of the entries by NA, where I make sure that these entries were not replaced earlier. This way, I get 10 folds of cross-validation data frames and 10 RMSEs for each numerical variable.

I have the following questions:

1. Are there any methodical flaws in this approach?

2. In particular, did I implement the two-level linear model in a correct way? (Should I have introduced dummy variables?)

3. Is the RMSE per numerical feature, averaged over the 10 folds, a reasonable way to measure the quality of the imputation?

4. As I describe it, I only compare the imputed values of the first imputation (m=1) with the actual values. What's the most reasonable way to include multiple imputations in the RMSE?

Edit: Below you find my R script:

library("mice")
library("lattice")
library("mi")

#  determine the number of folds.
k <- 10

# reduce boys dataset to 5 features and 1 classifier
df <- boys[,c("age", "hgt"  ,"wgt", "bmi", "hc" ,"reg")]

number.features <- 5
factor.variable.name <- "reg"

numeric.feature.names <- c("age","hgt","wgt","bmi","hc")

# make sure region classifier is always there
df <- df[which(!(df[,factor.variable.name] %in% NA)),]

# # introduce numeric class variable (needed for mice)
l=unique(as.character(df[,factor.variable.name])) # determine unique levels
df[,"numeric.class"] <- as.numeric(factor(df[,factor.variable.name], levels=l)) # write numeric levels as new feature

# z-transform numeric features

df[,numeric.feature.names]<- log(df[,numeric.feature.names])
df[,numeric.feature.names]<- data.frame(lapply(df[,numeric.feature.names], FUN=scale))

# take the feature-wise non-NA indices to create CV sets.

idx.list <- list(which(!(df$age %in% NA)),which(!(df$hgt %in% NA)), which(!(df$wgt %in% NA)), bmi.NNA.idx <- which(!(df$bmi %in% NA)),which(!(df$hc %in% NA)) ) idx.meta.df <- data.frame(length = integer(number.features), remainder = integer(number.features), largerfolds=integer(number.features)) # store the number of the NA-features for (i in (1:number.features)){ idx.meta.df$length[i] <- length(idx.list[[i]])}

# divide number of NA features by number of folds; residual for making first few folds larger
idx.meta.df$foldsize <- idx.meta.df$length %/% k
idx.meta.df$largerfolds <- idx.meta.df$length %% k

# shuffle each index list

idx.list<- lapply(idx.list,sample)

#initialize empty list of k lists of indices
fold.idx <- rep(list((list())),k)

# initialize CV indices
for (i in (1:k)){
for (j in 1:number.features){
# needs special treatment for foldsize > length

# size of block of indices for fold number i of feature j is given by idx.meta.df$foldsize[j] # except for the first idx.meta.df$largerfolds[j] elements. They are 1 element larger.
if (i <= idx.meta.df$largerfolds[j]){ size <- idx.meta.df$foldsize[j] + 1
}
else{
size <- idx.meta.df$foldsize[j] } # assing block of 'size' indices fold.idx[[i]][[j]] <- idx.list[[j]][((i-1)*size +1):((i)*size)] } } folds.df <- rep(list((data.frame()),k)) # create CV folds for (i in (1:k)){ with.na.df <- df for (j in (1:number.features)){ with.na.df[fold.idx[[i]][[j]],numeric.feature.names[j]]<- NA } folds.df[[i]]<-with.na.df } imputed.folds.df <- rep(list((data.frame()),k)) ## prepare multilevel imputation for mice dry.imp <- mice(df, print=FALSE) predMatr <- dry.imp$pred
## use all numeric features to predict every other features (even those which don't have NAs in the original data)
predMatr[1:number.features,1:number.features] <- (matrix(data=1, nrow=number.features, ncol=number.features) -diag(number.features))
predMatr[,"numeric.class"] <- -2
predMatr[,"reg"] <- 0
method <- c(rep("2l.norm",number.features),"", "")

# # make imputation
for (i in (1:k)){
# impute with mice, using 2 factor linear model which as numeric.class as class indicator
cat("two-level MICE for fold # ", i, "\n")
imp <- mice(folds.df[[i]], meth=method, pred=predMatr, maxit=30, m=1)
imputed.folds.df.mice[[i]] <- complete(imp)

}

# calculate RMSE and other statistics per fold
CV.results <- data.frame(RMSE1mice=numeric(0), RMSE2mice=numeric(0), RMSE3mice=numeric(0), RMSE4mice=numeric(0), RMSE5mice=numeric(0))

for (i in (1:k)){
for (j in (1:number.features)){

CV.results[i,j]<- mean((imputed.folds.df.mice[[i]][fold.idx[[i]][[j]],numeric.feature.names[j]] - df[fold.idx[[i]][[j]],numeric.feature.names[j]])^2)
}
}

• Would it be possible for you to put more of your code in the question? As it is, it doesn't really qualify as a MWE. – gregmacfarlane May 13 '14 at 20:30
• @gmacfarlane: I can do that. Most of the questions I have are not concerned with the implementation, but the methodical approach. I'll paste a minimal version of my R script if it helps. – Roland May 14 '14 at 8:16
• Is there anything else besides the bounty to increase the quality of my question and/or increase the interest? – Roland May 16 '14 at 13:14