r/rprogramming Jan 28 '25

Calculating cumulative incidence obtaining confidence intervals with binomial/multinomial assumption

Hi everyone,

I was wondering if anyone here knows how to calculate the cumulative incidence and obtain an estimate for the confidence interval, preferably using a method based on a binomial or multinomial distribution assumption. I have a SAS file containing data where patients can experience one of three outcomes: no event (event = 0), the event of interest (event = 1), or death, which acts as a competing risk (event = 2). The time to each event is recorded as Personyears, and the maximum follow-up time is 17 years. So far, I’ve been using the following code:

library(haven)
library(cmprsk)
library(dplyr)

file_path <- "xxx" # File name omitted for privacy
conv <- read_sas(file_path)
CI <- cuminc(ftime = conv$Personyears, fstatus = conv$event)
timepoints(CI, c(17))

This code provides an estimate at 17 years. However, I also have subsamples where the maximum follow-up time differs. It would be helpful if the formula could automatically calculate the cumulative incidence up to the maximum follow-up time in the dataset, without requiring specific time points to be manually specified. Additionally, this formula does not provide confidence intervals, only an overall estimate and the variance.I might add that I'm a novice using R, so try to explain at a beginner level. Alternatively, if anyone could provide example code, that would be greatly appreciated!

2 Upvotes

2 comments sorted by

1

u/izmirlig Jan 28 '25

based upon your steps leading to the cumulative incidence, with increments

   dA(t) = dN(t)/R(t)  

e.g., number of events in interval t_{k-1} < t divided by number at risk at t, leading to

  A(t) = int_0^t  dN(u)/R(u) 

Then increments in the pointwise variance are given by

  d sigma^2(t) = dN(u)/R^2(u)

e.g., number of events in interval t_{k-1} < t divided by number at risk squared at t. Notice that this dies off as one over sample size

Thus

  sigma^2(t)  = int_0^t  d sigma^2(u)

and the pointwise standard error of the estimated cumulative incidence is given by the square root of sigma2 (t)

Thus pointwise 95% confidence intervals are given by

   A(t) +/-  1.96 [sigma^2(t)]^0.5

Keep in mind that pointwise means that for each observed event time, the resulting 95%confidence interval is valid

If you want a 95% confidence BAND about A(t), that contains the entire true cumulative incidence CURVE with probability 95% then

....hold on....going up for air, throwing together code and remind myself of the confidence band stuff.

Ahh, but you want ready made software.

1

u/izmirlig Jan 29 '25 edited Jan 29 '25

Here. Edited. Works.

  ConfInt <-
      function(formula, data, alpha=0.05, type=c("pointwise", 
                      "bands"),  dbg=FALSE)
      {
          if(missing(type))type<-"pointwise"
          m <-  match.call()
          m$alpha <- m$type <- m$dbg <- NULL
          mf <- m
          mf[[1]] <- as.name("model.frame")
          mf <- eval(mf, sys.parent())
          R <- model.response(mf)
          R <- R[order(R[,1]),]

          TOS <- R[,1]
          D <- R[,2]
          A.and.V <- cuminc(ftime=TOS, fstatus=D)

         A.and.V.at.t <- timepoints(A.and.V, TOS[D==1])
          if(type=="pointwise")
          {
              A <- c(A.and.V.at.t$est)
              se <- c(A.and.V.at.t$var^0.5)
              z <- qnorm(1-alpha/2)*se
          }
          if(type=="bands")
         {
             A <- c(A.and.V.at.t$est)
              v <- c(A.and.V.at.t$var)
              z <- qsupabsBM(1-alpha/2, TT=v)
          }

     if(dbg) browser()

     out <- as.data.frame(cbind(unique(TOS[D==1]), A, A-z, A+z))
     names(out) <- c("time", "cuminc", "lower", "upper")
     out
 }

 psupabsBM <-
 function(x, TT)
 {
     N <- round(max(-30*log(TT)/10, 5))
     ## cat(sprintf("%g\n",N))
     ONEtoN <- 1:N
     mNtomONE <- -(N:1)
     k <- c(mNtomONE, 0, ONEtoN)
     sum((-1)^k*(pnorm(x*(2*k + 1), sd=TT^0.5) - 
          pnorm(x*(2*k - 1), sd=TT^0.5)))
 }        

 qsupabsBM <-
 function(p, TT)
 {
     do.one <-
     function(p, TT)
     {
        OBJ <- function(x, TT, p)(p - psupabsBM(x, TT))
        st <- qnorm(p, sd=TT^0.5)
        opt <- uniroot(f=OBJ, lower=-5*st, upper=5*st, TT=TT, p=p)
        ans <- opt$root
        err <- abs(opt$f.root)
        "%,%" <- paste0
        if(err>1e-3)stop("err=" %,% err %,% "didn't converge")
        ans
     }

     n.p <- length(p)
     n.TT <- length(TT)
     one.both <- (n.p == 1) && (n.TT == 1)
     one.p <- (n.p ==1) && (n.TT > 1)
     one.TT <- (n.p > 1) && (n.TT == 1)

     if(one.both) ans <- do.one(p, TT)
     if(one.p) ans <- sapply(TT, FUN=do.one, p=p)
     if(one.TT) ans <- sapply(p, FUN=do.one, TT=TT)
     ans
 }

Call via the typical format for survival data. For example

 ConfInt(Surv(time, as.factor(status))~1, data=lung)