-1

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
Community
  • 1
  • 1
  • 1
    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) – user20650 Jan 22 '17 at 19:20
  • 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")` – IRTFM Jan 23 '17 at 05:28
  • 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 – Andreu Ferrero Jan 23 '17 at 15:40
  • 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()`. – Andreu Ferrero Jan 23 '17 at 16:10
  • @AndreuFerrero ; so the maintainer responded and just said to create the `coef` and `vcov` methods? When i did this in the [comment above](http://stackoverflow.com/questions/41794649/can-mice-handle-crr-fine-gray-model#comment70778357_41794649) 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) – user20650 Jan 23 '17 at 17:42
  • `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))` – Andreu Ferrero Jan 23 '17 at 19:35
  • 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 – Andreu Ferrero Jan 23 '17 at 19:39
  • 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 – Andreu Ferrero Jan 23 '17 at 19:47
  • Yes, there is an issue with the `with`, likely because you are using `cov1=test1` in the formula which has missing data. – user20650 Jan 23 '17 at 19:55
  • `vcov.crr <- function(object, ...) object$var`or `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... – Andreu Ferrero Jan 25 '17 at 17:53
  • @AndreuFerrero ; Yes, it is easy to check - just look at the results of each function, however, I would expect differences re comment [here](http://stackoverflow.com/questions/41794649/can-mice-handle-crr-fine-gray-model/41815070#comment70821529_41815070). I dont use stata but there are different approaches towards imputation. – user20650 Jan 25 '17 at 20:56
  • 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. – Uwe Jan 26 '17 at 09:22

1 Answers1

0

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

user20650
  • 24,654
  • 5
  • 56
  • 91
  • `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')` – Andreu Ferrero Jan 23 '17 at 22:02
  • 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) – user20650 Jan 23 '17 at 22:09
  • Thanks! I will try ! – Andreu Ferrero Jan 23 '17 at 22:13
  • 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) – user20650 Jan 23 '17 at 22:17
  • 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! – Andreu Ferrero Jan 28 '17 at 13:05
  • your very welcome - glad the results are similar; we must be doing something right ;) – user20650 Jan 28 '17 at 13:12
  • The first 5 decimals are just the same in coefficients and standard errors! – Andreu Ferrero Jan 28 '17 at 13:30
  • That's impressive, particularly as there is a stochastic element, so all good – user20650 Jan 28 '17 at 17:11
  • 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. – Andreu Ferrero Jan 28 '17 at 17:33
  • 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? – user20650 Jan 28 '17 at 17:43
  • The first 5 decimals again, without pooling. – Andreu Ferrero Jan 28 '17 at 18:43
  • `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?** – Andreu Ferrero Jul 15 '19 at 15:21
  • @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. – user20650 Jul 15 '19 at 22:07
  • 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` @user20650 – Andreu Ferrero Jul 16 '19 at 06:10
  • @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. – user20650 Jul 16 '19 at 21:40
  • 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.](https://www.ncbi.nlm.nih.gov/pubmed/23653257) – Andreu Ferrero Jul 17 '19 at 10:11