I have five imputed datasets created with MICE in R, and am running run some post hoc analyses using the lsmeans package. Although MICE has great functions to easily pool and compare models (e.g. pool() and pool.compare()), they won't work here.

This leaves me in a bind regarding how to pool across lsmeans contrasts computed for each imputed dataset. I end up with a list of five lsmobj objects that output the following (each with slightly different values of course):

contrast          estimate       SE  df t.ratio p.value
 tpol.tpole.post -12.523126 4.021919 360  -3.114  0.0037
 tpol.tpole.fu    -9.557416 4.542045 360  -2.104  0.0621

P value adjustment: mvt method for 2 tests

contrast          estimate       SE  df t.ratio p.value

After some basic data manipulation to combine the respective rows from each imputation together into a list of data frames, I end up with:


         contrast  estimate       SE  df   t.ratio     p.value
1 tpol.tpole.post -12.52313 4.021919 360 -3.113719 0.003718103
2 tpol.tpole.post -12.66746 4.056950 360 -3.122411 0.003650298
3 tpol.tpole.post -12.74482 4.087365 360 -3.118101 0.003681073
4 tpol.tpole.post -13.32695 4.009558 360 -3.323796 0.001839164
5 tpol.tpole.post -12.38995 4.045162 360 -3.062906 0.004383206

       contrast   estimate       SE  df   t.ratio    p.value
1 tpol.tpole.fu  -9.557416 4.542045 360 -2.104210 0.06211271
2 tpol.tpole.fu  -9.526316 4.524851 360 -2.105333 0.06277749
3 tpol.tpole.fu  -8.240829 4.545330 360 -1.813032 0.11865308
4 tpol.tpole.fu -10.539075 4.508579 360 -2.337560 0.03500318
5 tpol.tpole.fu  -9.791866 4.490159 360 -2.180739 0.05187948

What would be an appropriate way to pool columns to get combined β estimates, SE, and significance tests?

Let's just tackle the first data frame above to keep it simple.

So far, my strategy is to use MICE's pool.scalar() function to combine estimates:

lsm.1 <- lsm.restruc[[1]] 
pooled.1 <- pool.scalar(lsm.1[,2], lsm.1[,3]^2, unique(lsm.1[,4]))

I have three questions about this approach:

1) pool.scalar needs a vector of variances as the second argument, so I square the SE estimates. Is this valid?

2) What should the sample size (third) parameter be based on here? Currently I'm entering df straight from the lsmeans analysis, but this is not sample size as requested by the function, so I suspect the pooled df will be incorrect. How do I determine the correct sample size? For example I had 183 participants in my study, randomised to three groups and assessed at three time points, contributing 549 individual data points altogether. Do I specify the total number of participants (183), the total number of individual data points (549), the number of participants compared in this particular contrast (123), the number of individual data points compared in this contrast (246), or something else... Clearly a bit confused here.

The other issue is how to pool the p-values. I assume it's best to calculate a new p-value from the pooled estimates, rather than trying to pool existing p-values. From pooled.1 I get an adjusted df, and can compute a new t-statistic from $\bar{Q} / \sqrt{T}$. Therefore to get a new p-value I run:

pooled.p <- 2 * pt(-(abs(pooled.1$qbar / sqrt(pooled.1$t))), pooled.1$df)

# Bonferroni correction to match the original correction applied by lsmeans 
pooled.p <- 2 * pooled.p

3) Are there any misconceptions or errors in the above logic or p-value calculations?

Any guidance or suggestions would be most welcome, particularly around the second question above.

  • $\begingroup$ I have never used the lsmeans package, so this is not the answer, but a couple of suggestions that have worked for me. Given that you have a list of results, have you tried pooling results with the mitools package? If not, you could try to use the MIextract and MIcombine functions on your list. Another package I have used to pool MI data is miceadds. Hope this helps! $\endgroup$
    – san
    Feb 3, 2016 at 5:17
  • $\begingroup$ mice also supports the "contr.treatment" function where you may need to manually execute each pairwise comparison (by changing "base="). I too am curious about the availability of a function or package that supports automated post-hoc testing for imputed data in R. $\endgroup$
    – user81715
    Mar 30, 2018 at 2:10

1 Answer 1


package emmeans supports pooling of model-estimated means:

cimputed_mice=mice(data=cmiss[,-1], method="rf", maxit = 10, seed=123) 

## perform a linear regression analysis on each imputed dataset
# and store the results

multi_lm = with(cimputed_mice, lm(cesd~MOSS+Income+decrease_inc+age+prescriptions+Chemo))

## pool the modeling results across imputations

## pooled summary    

##pooled R squared

## model-estimated means

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.