tagteam / riskRegression

R package for risk regression and prediction with censored data
44 stars 17 forks source link

Error when using predictRisk.aalen #39

Closed RobinDenz1 closed 2 years ago

RobinDenz1 commented 2 years ago

Dear package maintainers,

first of all I would like to thank you again for creating this very useful R-package. In particular, I am a big fan of the predictRisk function and have incorporated it into two of my own R-packages.

When playing around with different models I noticed that the predictRisk.aalen method of the predictRisk function currently does not work. It only returns the error FIXME, no matter what kind of inputs are passed to it. From the source code (and the error message) it is apparent that this is intentional (There is an unconditional stop() function right at the start of the function body).

Since I would like to use this feature I was wondering what exactly the problem with the current implementation is. Is there something fundamentally wrong with it or is it just a coding issue? If it's the latter I could try to fix it.

With kind regards, Robin Denz

tagteam commented 2 years ago

Hej Robin,

at some point these functions (predictRisk.aalen and predictRisk.cox.aalen) stopped working and hence the FIXME's. I do not use either of them and will remove them from the next version of riskRegression unless of course someone (e.g., you) spends the time to repair these functions. The current broken version tries to pull c-code from the almost obsolete package pec. But, all is opensource and you can start with what was done. However, mayhaps, the timereg package which fits these objects has developed own predict functions which could be adapted? This would probably be better, because then the predictRisk functions would not fail next time timereg is updated. best Thomas

RobinDenz1 commented 2 years ago

Thank you for the fast reply. There is indeed a predict method in the timereg package now. I have written three short wrapper functions around them which could be used to replace the faulty predictRisk methods. Here is the code I am proposing:

## * predictRisk.aalen
##' @export
##' @rdname predictRisk
##' @method predictRisk aalen
predictRisk.aalen <- function(object, newdata, times, ...) {

  if (is.null(object$B.iid)) {
    stop("Need resample.iid=1 when fitting the aalen model to make",
         " predictions. Please re-fit the aalen model with resample.idd=1.")
  }

  out <- timereg::predict.aalen(object=object,
                                newdata=newdata,
                                times=times,
                                se=FALSE,
                                ...)$S0
  return(1 - out)
}

## * predictRisk.cox.aalen
##' @export
##' @rdname predictRisk
##' @method predictRisk cox.aalen
predictRisk.cox.aalen <- function(object, newdata, times, ...) {

  out <- timereg::predict.cox.aalen(object=object,
                                    newdata=newdata,
                                    times=times,
                                    se=FALSE,
                                    ...)$S0
  return(1 - out)
}

## * predictRisk.comp.risk
##' @export
##' @rdname predictRisk
##' @method predictRisk comp.risk
predictRisk.comp.risk <- function(object, newdata, times, ...) {

  out <- timereg::predict.comprisk(object=object,
                                   newdata=newdata,
                                   times=times,
                                   se=FALSE,
                                   ...)$P1
  return(out)
}

And here are a few examples to show that they work as intended:

# aalen model
data(sTRACE)
out <- aalen(Surv(time, status==9) ~ sex + diabetes + chf + vf,
             data=sTRACE, max.time=7, n.sim=0, resample.iid=1)

predictRisk.aalen(object=out, newdata=sTRACE[1:5,], times=c(1, 2, 3))

# cox.aalen model
out <- cox.aalen(Surv(time,status==9) ~ prop(age) + prop(sex) +
                 prop(diabetes) + chf + vf,
                 data=sTRACE, max.time=7, n.sim=0, resample.iid=1)

predictRisk.cox.aalen(object=out, newdata=sTRACE[1:5,], times=c(1, 2, 3))

# comp.risk model
data(bmt) 

add <- comp.risk(Event(time, cause) ~ platelet + age + tcell, data=bmt, cause=1)
ndata <- data.frame(platelet=c(1, 0, 0), age=c(0, 1, 0), tcell=c(0, 0, 1))

predictRisk.comp.risk(object=add, newdata=ndata, times=c(1, 2, 3))

If you see any problems with these functions, please let me know and I will do my best to improve them accordingly.

tagteam commented 2 years ago

Excellent. Thank you very much! on github now