Estimating Missing Data with aregImpute() {R}

Soil scientists routinely sample, characterize, and summarize patterns in soil properties in space, with depth, and through time. Invariably, some samples will be lost or sufficient funds required for complete characterization can run out. In these cases the scientist is left with a data table that contains holes (so to speak) in the rows/columns that are missing data. If the data are used within a regression, missing values in any of the predictor or the response variable result in row-wise deletion-- even if 9/10 variables are present. Furthermore, common multivariate methods (PCA, RDA, dissimilarity metrics, etc.) cannot effectively deal with missing data. The scientist is left with a couple options: 1) row-wise deletion of cases missing any variable, 2) re-sampling or re-characterizing the missing samples, or 3) estimating the missing values from other variables in the dataset. This last option is called missing data imputation. This is a broad topic with countless books and scientific papers written about it. Here is a fairly simple introduction to the topic of imputation. Fortunately for us non-experts, there is an excellent function (aregImpute()) in the Hmisc package for R.

Below is an example of filling missing data in a soil characterization database with the aregImpute function. For each missing value, 10 candidate multiple imputations are returned. Otherwise, the function is using default parameters-- there are a lot of options, so reading the manual page is highly recommended! From the example below, it looks like we are able to adequately predict missing observations in most variables-- R² values are all > 0.5 - 0.6. Note that we are using the aregImpute function to automatically find the "best model" for predicting missing values (for each variable).

Implementation

## impute missing data: with aregImpute
# updated version of methods used in transcan
# multiple impution, requesting 10 candidate values per NA
x.ar <- aregImpute(~ L + A + B + clay + silt + sand + ph + fe_d + fe_o + mn_d + mn_o + Fe + Ca + K + Al + Si + Ti + Zr + Rb + S + Zn, data=x, n.impute=10)


#
# R-squares for Predicting Non-Missing Values for Each Variable Using Last Imputations of Predictors
# not bad!
#
    L     A     B  clay  silt  sand    ph  fe_d  fe_o  mn_d  mn_o    Fe    Ca
0.949 0.933 0.934 1.000 1.000 1.000 0.567 0.950 0.597 0.906 0.902 0.913 0.844
    K    Al    Si    Ti    Zr    Rb     S    Zn
0.860 0.839 0.829 0.885 0.886 0.885 0.680 0.730

I am interested in replacing missing data with the mean of the multiple imputations for each case. The following code below demonstrates one possible approach. However, this is not the suggested approach for incorporating the imputed values into subsequent analysis! Regression models should be iteratively fit to data containing a single value of each multiple imputation, and model coefficients combined according to rules for mixture distributions. (Thanks for the tip Cyrus). There are functions within the Hmisc, rms, and Zelig packages for automating these procedures.

Implementation (slightly improper use of multiple imputation)

# get a list of those variables with imputed values
imp.vars <- names(x.ar$imputed)

# compute mean imputed value for each variable
# and extract the original indices of the missing data, by variable
imp.vars.mean <- lapply(x.ar$imputed, function(i) apply(i, 1, mean))
imp.vars.idx <- lapply(imp.vars.mean, function(i) as.integer(names(i)))

# copy origial data
x.no.na <- x

# loop over imputed variables
for(i in imp.vars)
        {
        print(i)
        # get the mean imputations for this variable
        imp.i <- imp.vars.mean[[i]]
        # get the original indices for NA
        idx.i <- imp.vars.idx[[i]]
        # replace original NA with imputed values
        x.no.na[idx.i, i] <- imp.i
        }

Links:

Customized Scatterplot Ideas

R: advanced statistical package

Exploration of Multivariate Data