# Applying Rubin's rule for combining multiply imputed datasets

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.

• Are you just wondering what the symbols in that formula mean? Do you know what $U$ stands for in the SPSS output? – gung - Reinstate Monica Jul 7 '15 at 21:58
• For the ANOVA output in particular (as part of the linear regression), F-statistic is reported along with sums of squares, df, and mean squared. I am specifically trying to figure out how to compute the within imputation variance component for the F-statistic. – user81715 Jul 7 '15 at 22:54

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

library(miceadds)


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)


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)


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)

• Thank you for your suggestion. I have in fact tried miceadds but one issue I ran into was pooling the F value for all 20 imputed datasets and not just the 5 that mice is programmed to read by default (I believe). I have tried entering all 20 f-values I obtained from the SPSS output and used mice to pool these estimates. Would this method be acceptable (i.e. do imputation in SPSS, pool these estimates in mice adds)? If there is a way in which I could have all 20 imputations read in MICE, I'd appreciate it. Thanks again! – user81715 Jul 7 '15 at 23:37
• I have updated my answer with an complete example from spss to miceadds. – Thomas Jul 7 '15 at 23:41
• Thank you again. I will try running this and hope to follow-up with the results! – user81715 Jul 7 '15 at 23:47
• By the way: mice is not restricted to any number of imputations. I often use mice with m > 20 imputations. Getting the estimates out of SPSS and pooling them in mice is difficult and tedious, because mice assumes by default that all steps are performed inside the package scope. If my answer has helped you, please tag it as correct or upvote! – Thomas Jul 7 '15 at 23:47
• > dat$IV1=factor(dat$IV1) Error in $<-.data.frame(*tmp*, "IV1", value = integer(0)) : replacement has 0 rows, data has 10> dat.mids=as.mids(dat) Error in mice(data[data[, .imp] == 0, -c(.imp, .id)], m = max(as.numeric(data[, : Data should contain at least two columns > fit=mi.anova(mi.res=dat.mids,formula="DV~IV1*IV2*INT",type=3)summary(fit) Error: unexpected symbol in "fit=mi.anova(mi.res=dat.mids,formula="DV~IV1*IV2*INT",type=3)summary" > – user81715 Jul 8 '15 at 0:32 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. • can you explain how you got to 4.23 from 3.8 (mean of variance estimates for each imputed dataset) and 0.35 (variance of the coefficient estimates for each of the five datasets)? I am missing a step. – llewmills Feb 6 '18 at 6:21 • @llewmills hmm the DF correction must have been$(1+1/m)$and not$(m-1)/(m+1)$. – AdamO Feb 6 '18 at 14:48 • thank you @AdamO. I assume you meant$\frac{m+1}{m}$, (rather than$\frac{1+1}{m}\$) because that delivered the exact result as the pool(fit). – llewmills Feb 6 '18 at 23:30
• I don't suppose you know the formula for pooling the p-value @AdamO? – llewmills Feb 7 '18 at 0:12
• @llewmills why would you calculate 5 p-values? You have a grand mean and a grand SE, now combine them and calculate one p-value for the multiply imputed analysis. – AdamO Feb 7 '18 at 3:40