I’m new to using multiple imputations and I would like an opinion on using it with survival analysis in R.

I would like to perform a multiple imputation on data with missing values (mice package) and perform a survival time analysis (coxph). I want to perform a likelihood ratio test for a model with a spline(rms::rcs) in the covariates and a model without.

In this link (https://www.rdocumentation.org/packages/mice/versions/3.13.0/topics/pool.compare), it seems that 'glm / lm' models can be compared by using pool.compare (or D1/D3 function) corporated in the "mice" package.

#load package and data  


#Create dummy data /  "time to death" and "death" category to apply coxph 
time_death<- as.integer(runif(nrow(nhanes2), min = 100, max = 1000))
death  <- as.integer((runif(nrow(nhanes2), min = 0, max = 2)))
nhanes2 <- nhanes2 %>% mutate(death=death , time_death=time_death)

# mice
imp <- mice(nhanes2, m = 10, seed = 1, print = FALSE)

# produce the models to compare, a full model and
# an intercept only restricted model  
fit.imputed.full <- with(imp, lm(bmi ~ age + chl))
fit.imputed.res <- with(imp, lm(bmi ~ 1))  

# compare models using pool.compare()
pooled.comparison <- pool.compare(fit.imputed.full, fit.imputed.res) #Same to D1 or D3

#pooled.comparison <- D1(fit.imputed.full, fit.imputed.res)
pooled.comparison$pvalue  #<-- success

For survival time analysis, we have confirmed that "coxph function" can do the similar analysis.

# produce the cox proportional hazard models (after mice) to compare, a full model and control
fit.imputed.full_coxph <- with(imp, coxph(Surv(time_death,death)~ bmi + age + hyp+ chl))
fit.imputed.full_ctrl <- with(imp, coxph(Surv(time_death,death)~ bmi + chl))

pooled.comparison <- pool.compare(fit.imputed.full_coxph, fit.imputed.full_ctrl) #Same to D1 or D3
pooled.comparison$pvalue #<-success

But for models with restricted cubic spline, is it not possible to compare between the two models (model with rcs and without rcs)?

#Can we apply restricted cubic spline (rcs) to pool.compare (D1)??
# produce the cox proportional hazard models (after mice) to compare, a full model and control
fit.imputed.full_coxph_rcs <- with(imp, rms::cph(Surv(time_death,death)~ bmi + age + hyp+ chl))
fit.imputed.full_ctrl_rcs <- with(imp, rms::cph(Surv(time_death,death)~ rcs(bmi,k=3) + chl))

ddist <- datadist(nhanes2)

pooled.comparison <- pool.compare(fit.imputed.full_coxph_rcs, fit.imputed.full_ctrl_rcs) #Same to D1 or D3


Error $ operator is invalid for atomic vectors

Please see attached code. Thank you very much for your cooperation. Could you provide us with the best solution or suggestions?

  • $\begingroup$ Note that your fit.imputed.full_coxph_rcs and fit.imputed.full_ctrl_rcs aren't nested, as the former includes predictors not in the second while the second includes rcs() terms not in the first. I suspect that pool.compare() requires nested models, like those in your first 2 examples, where the predictors in one model are a subset of those in the other. $\endgroup$
    – EdM
    Sep 23, 2021 at 14:59
  • $\begingroup$ @EdM thank you for your answering and sincere reply. I understand the nested data in the model plays an important role to 'pool.compare'. However, could we perform to make nested data in rcs()? Do you know the examples? $\endgroup$ Sep 23, 2021 at 22:01

1 Answer 1


Some of this question seems to be specific to software implementation, which is off-topic on this site, but there are some statistical issues that should be addressed.

First, the attempt to compare fit.imputed.full_coxph_rcs and fit.imputed.full_ctrl_rcs via a likelihood-ratio test, as shown in the question, isn't valid. Such tests can only done on nested models in which the predictors of one model form a subset of the other model. That's not the case for those two models as shown, as the former includes predictors not in the second while the second includes rcs() terms not in the first. I can't say for sure that explains the error message, but check with truly nested models first in which the only difference in predictors has to do with rcs() terms.

Second, it's not clear that the pool.compare() function is the best way to do such a likelihood-ratio test. According to the documentation you cite, the function has been deprecated in favor of a D3() function that implements the methods for such tests described by Van Buuren. Combining likelihood-ratio tests among multiply imputed data sets requires some care.

Third, you might consider doing Wald tests on the rcs() terms instead. In the rms package many tests are implemented as Wald tests on the coefficient estimates and their covariances. In particular, that is how anova() works on cph objects; if such a model includes rcs() terms, the test for its nonlinear terms is a Wald test. Some notes in your code suggest that you have specified a D1()-type function for the model comparison. That's a Wald test, not a likelihood-ratio test.

A few hints for getting beyond this problem.

First, once the rms package is loaded the rcs() function should be accessible to the standard coxph() function. See if that works.

Second, the rms package inherits its own imputation function, aregImpute(), from the Hmisc package. That's designed to work directly with rms functions like cph(). You could thus consider working completely within the rms package for this.

Third, the functions used by the mice package can be quite picky about the way that model objects are constructed. Make sure to read the documentation very carefully; make sure that all required packages are installed and working properly and that the model objects you feed to functions in mice have the values and corresponding names that the functions are looking for. This could be something as simple as your code setting the datadist option after you did the cph() models; perhaps without the datadist set first the cph objects lack some critical value needed to work with mice.

Finally, think about why you really need to do this comparison. If you have a large enough data set and the cubic spline fits OK, how important is it to see if that does any better than a simple linear fit?

  • $\begingroup$ Thank you for your sincere reply! Your reply helps me lots! $\endgroup$ Oct 13, 2021 at 14:16

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.