2

I'm wondering whether the pool() function from R's mice package correctly pools shared frailty models fitted on multiply imputed datasets using coxme(). See the following example:

# needed to tidy the coxme object
install.packages("remotes"); remotes::install_github("junkka/ehahelper")

# assuming the following packages are installed
library(survival)
library(coxme)
library(mice)
library(eha)
library(ehahelper)

set.seed(1234)

# prepare the dataframe by adding a cluster id
df <- stanford2
df$cid <- round(df$id / 10) + 1
head(df)

# impute, fit and pool
imp <- mice(df, print = FALSE)
fit <- with(imp, coxme(Surv(time, status) ~ age + t5 + (1 | cid)))
fit1 <- pool(fit)

# check the p values
summary(fit1, exponentiate = TRUE, conf.int = TRUE)
#  term estimate  std.error statistic        df   p.value        2.5 %       97.5 %
#1  age 1.029322 0.01064437 2.7150698 0.4995043 0.3856764 1.769147e-01 5.988783e+00
#2   t5 1.182635 0.18427536 0.9102967 0.4442556 0.6418694 1.210586e-26 1.155330e+26

The p values and confidence intervals seem overly large, and I suspect that pool() might not handle multiple coxme objects correctly (i.e. performs "Wald tests" with adjusted degrees of freedom). If I fill the pooled parameter estimates into a "template" coxme object and let the summary() function calculate the usual z and p values then the results look more sensible to me:

# create a template fit (e.g. by fitting the same frailty model on complete cases) 
fit0 <- coxme(Surv(time, status) ~ age + t5 + (1 | cid), data = df)

# copy parameter estimates from the pooled fit into the template fit object
vec <- fit1$pooled$estimate
names(vec) <- fit1$pooled$term
fit0$coefficients <- vec

# the same for the pooled estimates of the variances
vec <- fit1$pooled$t
names(vec) <- fit1$pooled$term
fit0$variance <- diag(vec)

# check the p values again, and compare with those above
summary(fit0)
#         coef exp(coef)   se(coef)    z      p
#age 0.0289002  1.029322 0.01064437 2.72 0.0066
#t5  0.1677452  1.182635 0.18427536 0.91 0.3600

Questions:

  1. Is it correct to perform Wald tests with pooled estimates and pooled standard errors?
  2. If not, what would be the correct way to pool coxme objects? Does anyone know of an R function or some code snippet to achieve this?

Thank you for your inputs!

schotti
  • 161
  • 1
  • 9
  • Aside: The advantage of letting ```summary()``` calculate z and p rather than doing it manually is that the model can be plotted using ```plot_model()``` from the ```sjPlot``` package. – schotti May 21 '21 at 10:03
  • 1
    Thanks for providing a simple, reproducible example. That made it a lot easier to identify the likely source of the difficulty. – EdM May 23 '21 at 17:29

1 Answers1

1

In the results from both fit1 and fit0, the coefficients (in the non-exponentiated scale), their standard errors, and the statistics calculated from their ratios are identical. So the pooling of point estimates and standard errors agrees between pool() and your method.

What's different is the degrees of freedom used to evaluate the significance of the statistics.

In your evaluation of fit0 you used a Z-test, or equivalently a t-test with infinite degrees of freedom, to evaluate whether the average coefficient values differ from 0. That is the limiting case if there were no error introduced by imputation and a very large number of degrees of freedom in the complete-data analysis. Yet there is some error being introduced by the imputation.

At the other extreme, if all the error is introduced by imputation, you only have 5 imputations (the default) among which to average, providing 4 degrees of freedom in a t-test to distinguish those mean coefficient estimates from 0. On that basis, the test statistics for age and t5 would have p-values of 0.053 and 0.041, respectively.

The proper calculation of the degrees of freedom, $\nu$, for multiple imputation is explained in Stef van Buuren's Flexible Imputation of Missing Data book online. That calculation takes into account the number of imputations and the variability introduced by the imputation.

When pool() estimated the degrees of freedom, it only got 0.50 and 0.44 for age and t5, respectively, as shown in summary(fit1). The reported p-values are consistent with the reported test statistics and t-distributions having those degrees of freedom.

There thus might be an error in the calculation of degrees of freedom because you started with a coxme model. By default, the pool() function tries to get the complete-data degrees of freedom (typically, the number of cases minus the number of fitted coefficients) from the models themselves. The degrees of freedom easily extracted from the coxme objects (e.g. fit$analyses[[4]]$df) are very small and probably aren't what's appropriate for calculating the degrees of freedom for multiple imputation. The pool() function allows for a dfcom argument where you can specify the degrees of freedom in the complete-data analysis. Provide a reasonable value for that argument in your call to pool()* and you should be OK. My guess is that the results will be much closer to your estimates based on the Z-test.

*I'm not sure how many degrees of freedom to associate with the fixed-effect estimates in a coxme model. With 113 events in this case, 2 fixed coefficients, and a single random effect based on 19 clusters, something on the order of 100 would probably be OK.

EdM
  • 57,766
  • 7
  • 66
  • 187
  • Thanks a lot also for this answer, EdM. It's the dfs, right! What I don't quite understand is, ```coxme``` on complete cases does Wald tests equivalent to t tests with infinite dfs. (Accordingly, I can't find any dfs in e.g. the ```coxme``` object fit0 (except for the $\chi^2$ tests)). Doesn't this mean that the dfs after multiply imputing should also be infinite, and the Wald tests should still be appropriate? With van Buuren: $\nu_\mbox{com} =$ dfs in the complete data $= \infty$, therefore $\nu_\mbox{obs} = \infty$ according to formula (2.31), and $\nu = \nu_\mbox{old} = \infty$ (2.32)? – schotti May 23 '21 at 17:40
  • @schotti yes, the fixed-effect tests in `coxme` are based on a Wald test: a simple test relying on asymptotic normality that doesn't always hold with small data samples. There isn't, however, a near-infinite amount of information going into the coefficient estimates, as there would be if thousands of cases led to an effectively infinite df in a t-test on coefficients. At some level, imputation needs to take the ratio of cases to predictors into account. That's a general issue with models fit by MLE, but I'm not much of an expert; consult van Buuren and others on pooling such models. – EdM May 23 '21 at 18:30
  • According to Barnard and Rubin, for large $\nu_\mbox{com}$ we have $\nu_\mbox{obs} \approx \nu_\mbox{com} \cdot (1-\lambda)$, and $\nu$ is half the harmonic mean of $\nu_\mbox{obs}$ and $\nu_\mbox{old}$. This basically reduces dfs and increases p values. But this magic cannot happen if $\nu_\mbox{com}$ is assumed to be $\infty$ (as do the fixed-effect tests in ```coxme```) because then $\nu_\mbox{obs}$ and $\nu$ are $\infty$ as well. In other words, ```coxme```’s somewhat crude approach to testing makes it impossible to account for the additional variation introduced by imputation. – schotti May 23 '21 at 23:52
  • I assume ```coxme``` does so based on the rationale (or "hope"?) that mixed effect cox regression analyses usually include enough observations so that additional cases wouldn't make a big difference anymore. And as long as I'm sure that adjusting the dfs for imputation doesn't make them drop too far, say below 100 or so, the same approximation should work for me, since qnorm(0.975) and qt(0.975, df = 100) are still fairly close. To check this, I can set the dfcom argument in the ```pool()``` function, and if the dfs drop lower, I can still replace the z by a t test. Great, thanks again! – schotti May 23 '21 at 23:52
  • 1
    @schotti Cox models typically report Wald z-tests for coefficients; that's not specific to mixed modeling with `coxme`. The degrees of freedom for that asymptotic test don't represent the actual number of cases behind the results. Mixed models have an additional problem, as it's not alway clear how to calculate degrees of freedom even when the coefficient tests are t-tests, as in a standard [linear mixed model](https://stats.stackexchange.com/q/264437/28500). So `coxme` faces both challenges. A reasonable choice for the dfcom argument seems like the best work-around to this non-expert. – EdM May 24 '21 at 13:06