How to pool results from post hoc lsmeans analysis across multiple imputations with MICE

Question

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):

[[1]]
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

[[2]]
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:

lsm.restruc

[[1]]
         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

[[2]]
       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.