can mice() handle crr()? Fine-Gray model

Question

My doubt is if it is possible to pool multiple imputation data set, from "mice()", on a fit model of Fine-Gray from "crr()", and if it is statistically correct...

example

library(survival)
library(mice)
library(cmprsk)

test1 <- as.data.frame(list(time=c(4,3,1,1,2,2,3,5,2,4,5,1, 4,3,1,1,2,2,3,5,2,4,5,1), 
                            status=c(1,1,1,0,2,2,0,0,1,1,2,0, 1,1,1,0,2,2,0,0,1,1,2,0),
                            x=c(0,2,1,1,NA,NA,0,1,1,2,0,1, 0,2,1,1,NA,NA,0,1,1,2,0,1),
                            sex=c(0,0,0,NA,1,1,1,1,NA,1,0,0, 0,0,0,NA,1,1,1,1,NA,1,0,0)))

dat <- mice(test1,m=10, seed=1982)

#Cox regression: cause 1

models.cox1 <- with(dat,coxph(Surv(time, status==1) ~ x +sex ))                 

summary(pool(models.cox1))

#Cox regression: cause 1 or 2

models.cox <- with(dat,coxph(Surv(time, status==1 | status==2) ~ x +sex ))                 
models.cox
summary(pool(models.cox))


#### crr()

#Fine-Gray model

models.FG<- with(dat,crr(ftime=time, fstatus=status,  cov1=test1[,c( "x","sex")], failcode=1, cencode=0, variance=TRUE))                 

summary(pool(models.FG))

#Error in pool(models.FG) : Object has no vcov() method.

models.FG

G

There is no vcov method for crr models which mice needs. Check by looking at one model. (m = models.FG$analyses[[1]]) ; vcov(m) . But we can access this with models.FG$analyses[[1]]$var. Check standard errors against returned values for m, against this sqrt(diag(models.FG$analyses[[1]]$var)). So maybe write own vcov method (also need coef method): vcov.crr <- function(object, ...) object$var ; coef.crr <- function(object, ...) object$coef . Then run again summary(pool(models.FG)) (i have no idea if it is statistically correct to pool values for this model type) – Danu

In situations where thare is no vcov method, one must ask the question whether the package authors want users to assume that further statistical analysis is valid beyond what they support. The email is easily obtainable with: maintainer("cmprsk") – Milquetoast

I got this answer: I think you can easily create one; eg coef.crr=function(object,...) object$coef vcov.crr=getS3method('vcov','coxph') I think then vcov(crrobject) should work – G

Now I should create a as.mira() list to pool() then... I am trying to do it but I'm doing something wrong... With as.mira(fitlist) I can not do a proper structure list list(). – G

@AndreuFerrero ; so the maintainer responded and just said to create the coef and vcov methods? When i did this in the comment above I was able to use pool without any further tweaking. Can you edit your question to show why you are having to use as.mira (ps fitlist is not in your question) – Danu

coef.crr=function(object,...) object$coef vcov.crr=getS3method('vcov','coxph') #Fine-Gray model models.FG<- with(dat,crr(ftime=time, fstatus=status, cov1=test1[,c( "x","sex")], failcode=1, cencode=0, variance=TRUE)) 8 cases omitted due to missing values #summary(pool(models.FG)) coef.crr.x<-coef.crr(models.FG) vcov.crr.x<-vcov.crr(models.FG) summary(pool(vcov.crr.x)) – G

when I use with() it was warning about: 8 cases omitted due to missing values in each m=... So to be sure about the results I was thinking to set one by one model in each m= and use as.mira() to then pool() them – G

When I summary() models.FG, it shows me the same results in ech m=, like with() didn't use different models from different imputations... I saw it when comparing results from STATA.12 – G

Yes, there is an issue with the with, likely because you are using cov1=test1 in the formula which has missing data. – Danu

vcov.crr <- function(object, ...) object$varor vcov.crr=getS3method('vcov','coxph') give me same SE, is that right? no differences in methods... if it is like this, when the results are iqual in Satat.12 in coef but not in vcov... So, caveat... – G

@AndreuFerrero ; Yes, it is easy to check - just look at the results of each function, however, I would expect differences re comment here. I dont use stata but there are different approaches towards imputation. – Danu

Please edit your question and move the additional information and clarifications you gave in the comments to the question. Currently, your question is only a piece of code without any further explanations - Thank you. – Treasurehouse

Answer 1

There are a couple of things that need to be done to get this to work.

Your initial data and imputation.

library(survival)
library(mice)
library(cmprsk)

test1 <- as.data.frame(list(time=c(4,3,1,1,2,2,3,5,2,4,5,1, 4,3,1,1,2,2,3,5,2,4,5,1), 
                            status=c(1,1,1,0,2,2,0,0,1,1,2,0, 1,1,1,0,2,2,0,0,1,1,2,0),
                            x=c(0,2,1,1,NA,NA,0,1,1,2,0,1, 0,2,1,1,NA,NA,0,1,1,2,0,1),
                            sex=c(0,0,0,NA,1,1,1,1,NA,1,0,0, 0,0,0,NA,1,1,1,1,NA,1,0,0)))

dat <- mice(test1,m=10, print=FALSE)

There is no vcov method for crr models which mice requires, however, we can access the covariance matrix using the model$var returned value.

So write own vcov method to extract, and also need a coef method.

vcov.crr <- function(object, ...) object$var # or getS3method('vcov','coxph')
coef.crr <- function(object, ...) object$coef

There is also an error in how the model is passed to with.mids: your code has cov1=test1[,c( "x","sex")], but really you want cov1 to use the imputed data. I am not sure how to correctly write this as an expression due to the cov1 requiring a matrix with relevant variables, but you can easily hard code a function.

# This function comes from mice:::with.mids
Andreus_with <- 
function (data, ...) {
    call <- match.call()
    if (!is.mids(data)) 
        stop("The data must have class mids")
    analyses <- as.list(1:data$m)
    for (i in 1:data$m) {
        data.i <- complete(data, i)
        analyses[[i]] <- crr(ftime=data.i[,'time'], fstatus=data.i[,'status'],  
                         cov1=data.i[,c( "x","sex")], 
                         failcode=1, cencode=0, variance=TRUE)
    }
    object <- list(call = call, call1 = data$call, nmis = data$nmis, 
        analyses = analyses)
    oldClass(object) <- c("mira", "matrix")
    return(object)
}

---

EDIT:

The mice internals have changed since this answer; it now uses the broom package to extract elements from the fitted crr model. So tidy and glance methods for crr models are required:

tidy.crr <- function(x, ...) {
  co = coef(x)
  data.frame(term = names(co), 
             estimate = unname(co), 
             std.error=sqrt(diag(x$var)), 
             stringsAsFactors = FALSE)
}

glance.crr <- function(x, ...){ }

---

The above code then allows the data to be pooled.

models.FG <- Andreus_with(dat)                 
summary(pool(models.FG))

Note that this gives warnings over df.residual not being defined, and so large samples are assumed. I'm not familiar with crr so a more sensible value can perhaps be extracted -- this would then be added to the tidy method. (mice version ‘3.6.0’)

Danu

vcov.crr <- getS3method('vcov','coxph') coef.crr <- function(object, ...) object$coef Just like this for de other vcov method??? Both get same result... It is really similar to STATA.12 results, the coef are the same, still the vcov... Could you please, code for the: vcov.crr=getS3method('vcov','coxph') – G

Andreu, the only difference i nthe vcov functions (between mine and the maintainers) is that theirs add names to the matrix, so it wont effect the results. ( to see the code type getS3method('vcov','coxph') in to the terminal, you will see it uses object$var, as does mine) – Danu

Thanks! I will try ! – G

ps I would expect some differences in results as stata and R will be using different imputed values : you could try running it for longer with more chains (btw remember to set the seed in R to make your results reproducible) – Danu

Using same imputed data set in STATA.12 and the code you gave me the resuls now are more similar when I set none robust option for standard error!!! Thanks for you answer! – G

your very welcome - glad the results are similar; we must be doing something right ;) – Danu

The first 5 decimals are just the same in coefficients and standard errors! – G

That's impressive, particularly as there is a stochastic element, so all good – Danu

I imputed the missing data in R with mice() and then I pool the results with your code in R. Finally I import the imputed data set from R to STATA and pool again the fine-gray model in STATA to see what happen with the coefficients and standard errors. Results were really similar if I set standard errors without robust option. I am thinking to do all the process only in STATA to see differences... I guess that if I find more different the pooled results it will be because the imputation method in STATA too. – G

ahh right okay, you used the same data.. well anyway equality up to 5 decimal places is pretty good. Are the results from the crr models exactly the same without pooling? – Danu

The first 5 decimals again, without pooling. – G

summary(pool(models.FG)) Error: No glance method for objects of class crr for package ‘mice’ version 3.5.0 ¿Should I go back to an old version? – G

@AndreuFerrero ; mice has had a major rewrite -- it now uses methods from broom to summarise the models. You'll need to step through the relevant mice function, similar to above, to see what needs to be tweaked - probably be as simple as writing a quick glance function for ccr model. Sorry, I dont have time just now. – Danu

I appreciate so much the effort you have done. I have used the function you wrote many times. I read something about broom. Could you give me a [link] about how write glance function for ccr model? I mean, an example of something similar. I will figure out how to do it, I am improving my ´code´. @user20650 @Danu – G

@AndreuFerrero ; please see update -- hopefully enough to get you started. If you have a ref or method on how to deal with the residuals please feel free to edit. – Danu

Wow! I just found this paper by now. Not easy. Trying about residuals. "...obtained the variance of the MI estimator from the usual MI variance estimator, combining between-imputation and within-imputation variances." Moreno-Betancur M, Latouche A. Regression modeling of the cumulative incidence function with missing causes of failure using pseudo-values. Stat Med. 2013 Aug 15;32(18):3206-23. – G