Applying Rubin's rule for combining multiply imputed datasets

Question

I am hoping to pool the results of a pretty basic set of analysis performed on a multiply imputed data (e.g. multiple regression, ANOVA). Multiple imputation and the analyses have been completed in SPSS but SPSS does not provide pooled results for a few statistics including the F-value, covariance matrix, R-Squared etc.

I have made a few attempts to address this issue by venturing into R or trying macros that are available and have not successfully resolved the problem (e.g. with running into issues with pooling the stats for more than 5 imputations in Mice, for example).

At this point, I would like to try computing these by hand, applying Rubin's rule, using the output that SPSS generates. However, I am not sure how I can derive the within-imputation variance ($\bar U = \frac 1 m\sum_{j=1}^mU_j$) based on the output SPSS generates.

I would really appreciate a detailed instruction on this.

Answer 1

Rubin's rules can only be applied to parameters following a normal distribution. For parameters with a F or Chi Square distribution a different set of formulas is needed:

• Allison, P. D. (2002). Missing data. Newbury Park, CA: Sage.

For performing an ANOVA on multiple imputed datasets you could use the R package miceadds (pdf; miceadds::mi.anova).

Update 1

Here is a complete example:

1. Export your data from SPSS to R. In Spss save your dataset as .csv

2. Read in your dataset:

   library(miceadds)   
   dat <– read.csv(file='your-dataset.csv')
   

   Lets assume, that $reading$ is your dependent variable and that you have two factors

   • gender, with male = 0 and female = 1
   • treatment, with control = 0 and 'received treatment' = 1

   Now lets convert them to factors:

   dat$gender    <- factor(dat$gender)
   dat$treatment <- factor(dat$treatment)
   

3. Convert your dataset to a mids object, wehere we assume, that the first variable holds the imputation number (Imputation_ in SPSS):

   dat.mids <- as.mids(dat)
   

4. Now you can perform an ANOVA:

   fit <- mi.anova(mi.res=dat.mids, formula="reading~gender*treatment", type=3)
   summary(fit)
   

Update 2 This is a reply to your second comment:

What you describe here is a data import/export related problem between SPSS and R. You could try to import the .sav file directly into R and there are a bunch of dedicated packages for that: foreign, rio, gdata, Hmisc, etc. I prefer the csv-way, but that's a matter of taste and/or depends on the nature of your problem. Maybe you should also check some tutorials on youtube or other sources on the internet.

library(foreign)
dat <- read.spss(file='path-to-sav', use.value.labels=F, to.data.frame=T)

Update 3 This is a reply to your first comment:

Yes, you can do your analysis in SPSS and pool the F values in miceadds (please note this example is taken from the miceadds::micombine.F help page):

library(miceadds)
Fvalues <- c(6.76 , 4.54 , 4.23 , 5.45 , 4.78, 6.76 , 4.54 , 4.23 , 5.45 , 4.78, 
             6.76 , 4.54 , 4.23 , 5.45 , 4.78, 6.76 , 4.54 , 4.23 , 5.45 , 4.78 )
micombine(Fvalues, df1=4)

Answer 2

You correctly wrote down the pooled estimator:

$$ \bar{U} = \frac{1}{m} \sum_{i=1}^m U_i$$

Where $U_i$ represents the analytic results from the $i$-th imputed dataset. Normally, analytic results have some normal approximating distribution from which we draw inference or create confidence bounds. This is mainly done using the mean value ($U_i$) and its standard error. T-tests, linear regressions, logistic regressions, and basically most analyses can be adequately summarized in terms of that value $U_i$ and its standard error $\text{se}(U_i)$.

Rubin's Rules uses the law of total variance to write down the variance as the sum of a between and within imputation variance:

$$\text{var}(\bar{U}) = E[\text{var}(\bar{U}|U_i)] + \mbox{var}\left(E[\bar{U}|U_i]\right)$$

The first term is the within-variance such that $E[\text{var}(\bar{U}|U_i) = \frac{1}{m}\sum_{i=1}^m V_i$ where $V_i$ is the variance of the analysis result from the $i$-th complete or imputed dataset. The latter term is the between-imputation variance: $ \mbox{var}\left(E[\bar{U}|U_i]\right) = \frac{M+1}{M-1} \sum_{i=1}^m\left(U_i - \bar{U}\right)^2$. I've never quite grasped the DF correction here, but this is basically the accepted approach.

Anyway, since the recommended number of imputations is small (Rubin suggests as few as 5), it is typically possible to compute this number by hand fitting each analysis. A by-hand example is listed below:

require(mice)
set.seed(123)
nhimp <- mice(nhanes)
sapply(1:5, function(i) {
  fit <- lm(chl ~ bmi, data=complete(nhimp, i))
  print(c('coef'=coef(fit)[2], 'var'=vcov(fit)[2, 2]))
})

Gives the following output:

coef.bmi      var 
2.123417 4.542842 
3.295818 3.801829 
2.866338 3.034773 
1.994418 4.124130 
3.153911 3.531536

So the within variance is the average of the imputation specific point estimate variances: 3.8 (average of second column). The between variance is 0.35 variance of the first column). Using the DF correction we get variance 4.23. This agrees with the pool command given in the mice package.

> fit <- with(data=nhimp,exp=lm(chl~bmi))
> summary(pool(fit))
                  est        se        t       df   Pr(>|t|)     lo 95      hi 95 nmis       fmi     lambda
(Intercept) 119.03466 54.716451 2.175482 19.12944 0.04233303  4.564233 233.505080   NA 0.1580941 0.07444487
bmi           2.68678  2.057294 1.305978 18.21792 0.20781073 -1.631731   7.005291    9 0.1853028 0.10051760

which shows the SE = 2.057 for the model coefficient, (Variance = SE**2 = 4.23).

I fail to see how increasing the number of imputed datasets creates any particular issue. If you cannot supply an example of the error, I don't know how to be more helpful. But by-hand combination is certain to accommodate a variety of modeling strategies.

This paper discusses other ways that the law of total variance can derive other estimates of the variance of the pooled estimate. In particular, the authors point out (correctly) that the necessary assumption for Rubin's Rules is not normality of the point estimates but something called congeniality. WRT normality, most point estimates that come from regression models have rapid convergence under the central limit theorem, and the bootstrap can show you this.