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:


#  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
          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

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)
  • $\begingroup$ 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. $\endgroup$ – gregmacfarlane May 13 '14 at 20:30
  • 1
    $\begingroup$ @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. $\endgroup$ – Roland May 14 '14 at 8:16
  • $\begingroup$ Is there anything else besides the bounty to increase the quality of my question and/or increase the interest? $\endgroup$ – Roland May 16 '14 at 13:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.