(2024.20.10) better handling duplicated times predictCox with diag

DESCRIPTION

@@ -2,7 +2,7 @@ Package: riskRegression
 Type: Package
 Title: Risk Regression Models and Prediction Scores for Survival Analysis with
     Competing Risks
-Version: 2024.16.10
+Version: 2024.20.10
 Authors@R: c(person("Thomas Alexander", "Gerds", role = c("aut", "cre"),
                     email = "[email protected]"),
              person("Johan Sebastian", "Ohlendorff", role = "aut"),

NAMESPACE

@@ -24,6 +24,7 @@ S3method(confint,ate)
 S3method(confint,influenceTest)
 S3method(confint,predictCSC)
 S3method(confint,predictCox)
+S3method(confint,wglm)
 S3method(coxBaseEstimator,coxph)
 S3method(coxBaseEstimator,phreg)
 S3method(coxBaseEstimator,prodlim)
@@ -89,6 +90,7 @@ S3method(is.iidCox,default)
 S3method(is.iidCox,phreg)
 S3method(is.iidCox,prodlim)
 S3method(model.tables,ate)
+S3method(model.tables,wglm)
 S3method(nobs,CauseSpecificCox)
 S3method(nobs,multinom)
 S3method(nobs,wglm)
@@ -370,6 +372,7 @@ importFrom(stats,uniroot)
 importFrom(stats,update)
 importFrom(stats,update.formula)
 importFrom(stats,var)
+importFrom(stats,vcov)
 importFrom(stats,weights)
 importFrom(stats,wilcox.test)
 importFrom(survival,Surv)

NEWS

@@ -1,3 +1,9 @@
+Version: 2024.18.10 (Brice - 2024/18/10)
+- [feature] add confint and model.tables method to wglm for extracting coefficients with the confidence intervals and p-values.
+- [internal] wglm automatically compute the variance covariance matrix instead of doing it on the fly when calling summary. Can be disabled with argument se.
+- [internal] speed-up for iidCox with non-unique times: it internals works on unique times and expand the results at the end to match the user-request.
+- [internal] speed-up for predictCox with non-unique times when diag=TRUE: call iidCox with unique times to reduce memory usage.
+
 Version: 2024.09.10 (Brice - 2024/09/10)
 - add argument product.limit in predictCox and remove remove predictCoxPL. Side effect: while the influence function for the hazard is unchanged it the influence function relative to the survival scaled by the product limit survival (instead of exponential approximation).
 - fix bug in how ties are handled in ate when using the double robust estimator.

R/RcppExports.R

@@ -39,8 +39,8 @@ calcSeMinimalCox_cpp <- function(seqTau, newSurvival, hazard0, cumhazard0, newX,
     .Call(`_riskRegression_calcSeMinimalCox_cpp`, seqTau, newSurvival, hazard0, cumhazard0, newX, neweXb, IFbeta, Ehazard0, cumEhazard0, hazard_iS0, cumhazard_iS0, delta_iS0, sample_eXb, sample_time, indexJumpSample_time, jump_time, indexJumpTau, lastSampleTime, newdata_index, factor, nTau, nNewObs, nSample, nStrata, p, diag, exportSE, exportIF, exportIFmean, exportHazard, exportCumhazard, exportSurvival, debug)
 }
 
-calcAIFsurv_cpp <- function(ls_IFcumhazard, IFbeta, cumhazard0, survival, eXb, X, prevStrata, ls_indexStrata, factor, nTimes, nObs, nStrata, nVar, diag, exportCumHazard, exportSurvival) {
-    .Call(`_riskRegression_calcAIFsurv_cpp`, ls_IFcumhazard, IFbeta, cumhazard0, survival, eXb, X, prevStrata, ls_indexStrata, factor, nTimes, nObs, nStrata, nVar, diag, exportCumHazard, exportSurvival)
+calcAIFsurv_cpp <- function(ls_IFcumhazard, IFbeta, cumhazard0, survival, eXb, X, prevStrata, ls_indexStrata, ls_indexStrataTime, factor, nTimes, nObs, nStrata, nVar, diag, exportCumHazard, exportSurvival) {
+    .Call(`_riskRegression_calcAIFsurv_cpp`, ls_IFcumhazard, IFbeta, cumhazard0, survival, eXb, X, prevStrata, ls_indexStrata, ls_indexStrataTime, factor, nTimes, nObs, nStrata, nVar, diag, exportCumHazard, exportSurvival)
 }
 
 calculateDelongCovarianceFast <- function(Xs, Ys) {

R/calcSeCox.R

@@ -3,9 +3,9 @@
 ## author: Brice Ozenne
 ## created: maj 27 2017 (11:46) 
 ## Version: 
-## last-updated: sep 10 2024 (10:31) 
+## last-updated: Oct 20 2024 (13:50) 
 ##           By: Brice Ozenne
-##     Update #: 903
+##     Update #: 950
 #----------------------------------------------------------------------
 ## 
 ### Commentary: 
@@ -74,9 +74,21 @@ calcSeCox <- function(object, times, nTimes, type, diag,
     ## ** Computation of the influence function
     if(is.iidCox(object)){
         store.iid <- object$iid$store.iid
-        iid.object <- selectJump(object$iid, times = times, type = type)
+        if(diag && any(duplicated(times)) && store.iid[[1]] != "minimal"){
+            iid.object <- selectJump(object$iid, times = sort(unique(times)), type = type)
+            compress.time <- TRUE
+        }else{
+            iid.object <- selectJump(object$iid, times = times, type = type)
+            compress.time <- FALSE
+        }
     }else{
-        iid.object <- iidCox(object, tau.hazard = times, store.iid = store.iid, return.object = FALSE)
+        if(diag && any(duplicated(times)) && store.iid[[1]] != "minimal"){
+            iid.object <- iidCox(object, tau.hazard = sort(unique(times)), store.iid = store.iid, return.object = FALSE)
+            compress.time <- TRUE
+        }else{
+            iid.object <- iidCox(object, tau.hazard = times, store.iid = store.iid, return.object = FALSE)
+            compress.time <- FALSE
+        }
     }
 
     ## ** Prepare arguments
@@ -211,25 +223,28 @@ calcSeCox <- function(object, times, nTimes, type, diag,
 
             for(iStrata in 1:nStrata){ ## iStrata <- 1
                 indexStrata <- which(new.strata==iStrata)
+                if(compress.time){
+                    indexStrata.time <- match(times[indexStrata], iid.object$time[[iStrata]])                        
+                }else{
+                    indexStrata.time <- indexStrata
+                }
                 if(length(indexStrata)==0){next}
                 iPrevalence <- length(indexStrata)/new.n
 
                 ## compute iid
                 if("hazard" %in% type){
                     if(nVar.lp==0){
-                        iIFhazard <- iid.object$IFhazard[[iStrata]][,indexStrata,drop=FALSE]
+                        iIFhazard <- iid.object$IFhazard[[iStrata]][,indexStrata.time,drop=FALSE]
                     }else{
-                        iIFhazard <- rowMultiply_cpp(iid.object$IFhazard[[iStrata]][,indexStrata,drop=FALSE] + rowMultiply_cpp(X_IFbeta_mat[,indexStrata,drop=FALSE],
-                                                                                                                               scale = Lambda0$hazard[[iStrata]][indexStrata]),
+                        iIFhazard <- rowMultiply_cpp(iid.object$IFhazard[[iStrata]][,indexStrata.time,drop=FALSE] + rowMultiply_cpp(X_IFbeta_mat[,indexStrata,drop=FALSE], scale = Lambda0$hazard[[iStrata]][indexStrata]),
                                                      scale = new.eXb[indexStrata])
                     }
                 }
                 if("cumhazard" %in% type || "survival" %in% type){
                     if(nVar.lp==0){
-                        iIFcumhazard <- iid.object$IFcumhazard[[iStrata]][,indexStrata,drop=FALSE]
+                        iIFcumhazard <- iid.object$IFcumhazard[[iStrata]][,indexStrata.time,drop=FALSE]
                     }else{
-                        iIFcumhazard <- rowMultiply_cpp(iid.object$IFcumhazard[[iStrata]][,indexStrata,drop=FALSE] + rowMultiply_cpp(X_IFbeta_mat[,indexStrata,drop=FALSE],
-                                                                                                                                     scale = Lambda0$cumhazard[[iStrata]][indexStrata]),
+                        iIFcumhazard <- rowMultiply_cpp(iid.object$IFcumhazard[[iStrata]][,indexStrata.time,drop=FALSE] + rowMultiply_cpp(X_IFbeta_mat[,indexStrata,drop=FALSE], scale = Lambda0$cumhazard[[iStrata]][indexStrata]),
                                                         scale = new.eXb[indexStrata])
                     }
                     if("survival" %in% type && ("iid" %in% export || "average.iid" %in% export)){
@@ -337,6 +352,11 @@ calcSeCox <- function(object, times, nTimes, type, diag,
         if(is.null(new.survival)){
             new.survival <- matrix()
         }
+        if(compress.time){
+            new.indexStrataTime <- lapply(1:nStrata, function(iS){match(times[new.indexStrata[[iS]]+1], iid.object$time[[iS]])-1})
+        }else{
+            new.indexStrataTime <- new.indexStrata
+        }
 
         ## C++
         if("hazard" %in% type){
@@ -348,6 +368,7 @@ calcSeCox <- function(object, times, nTimes, type, diag,
                                               X = new.LPdata,
                                               prevStrata = new.prevStrata,
                                               ls_indexStrata = new.indexStrata,
+                                              ls_indexStrataTime = new.indexStrataTime,
                                               factor = factor,
                                               nTimes = nTimes,
                                               nObs = object.n,
@@ -366,6 +387,7 @@ calcSeCox <- function(object, times, nTimes, type, diag,
                                                  X = new.LPdata,
                                                  prevStrata = new.prevStrata,
                                                  ls_indexStrata = new.indexStrata,
+                                                 ls_indexStrataTime = new.indexStrataTime,
                                                  factor = factor,
                                                  nTimes = nTimes,
                                                  nObs = object.n,

R/iidCox.R

@@ -150,23 +150,30 @@ iidCox.coxph <- function(object, newdata = NULL,
                  "each element being the vector of times for each strata \n")
         }
 
-        need.order <- FALSE
+        need.order <- vector(mode = "logical", length = nStrata)
         tau.oorder <- vector(mode = "list", length = nStrata)
         etimes.max <- vector(mode = "numeric", length = nStrata)
-        for(iStrata in 1:nStrata){
+
+        for(iStrata in 1:nStrata){ ## iStrata <- 1
             etimes.max[iStrata] <- attr(tau.hazard[[iStrata]],"etimes.max")
-            need.order <- need.order + is.unsorted(tau.hazard[[is.unsorted]])
-            tau.oorder[[iStrata]] <- order(order(tau.hazard[[iStrata]]))
-            tau.hazard[[iStrata]] <- sort(tau.hazard[[iStrata]])
+            need.order[iStrata] <- is.unsorted(tau.hazard[[is.unsorted]])  || any(duplicated(tau.hazard[[is.unsorted]]))
+            if(need.order){
+                Utau.hazard[[iStrata]] <- sort(unique(tau.hazard[[iStrata]]))
+            }else{
+                Utau.hazard[[iStrata]] <- tau.hazard[[iStrata]]
+            }
+            tau.oorder[[iStrata]] <- match(tau.hazard[[iStrata]],Utau.hazard[[iStrata]])
         }
-        need.order <- (need.order>0)
-
     }else if(!is.null(tau.hazard)){
-        etimes.max <- attr(tau.hazard,"etimes.max")
+        etimes.max <- attr(tau.hazard,"etimes.max")        
 
-        need.order <- is.unsorted(tau.hazard)
-        tau.oorder <- lapply(1:nStrata, function(iS){order(order(tau.hazard))})
-        tau.hazard <- sort(tau.hazard)
+        need.order <- is.unsorted(tau.hazard) || any(duplicated(tau.hazard))
+        if(need.order){
+            Utau.hazard <- sort(unique(tau.hazard))
+        }else{
+            Utau.hazard <- tau.hazard
+        }
+        tau.oorder <- match(tau.hazard,Utau.hazard)
     }else{
         etimes.max <- NULL        
         need.order <- FALSE
@@ -367,15 +374,17 @@ iidCox.coxph <- function(object, newdata = NULL,
             }
 
             ## tau.hazard
-            if(is.null(tau.hazard)){
+            if(is.list(tau.hazard)){
+                iTau.oorder <- tau.oorder[[iStrata]]
+                tau.hazard_strata <- Utau.hazard[[nStrata]]
+            }else if(!is.null(tau.hazard)){
+                iTau.oorder <- tau.oorder
+                tau.hazard_strata <- Utau.hazard
+            }else{
                 tau.hazard_strata <- unique(object.time_strata[[iStrata]][object.status_strata[[iStrata]] == 1])
                 if(!is.null(tau.max)){
                     tau.hazard_strata <- tau.hazard_strata[tau.hazard_strata<=tau.max]
                 }
-            }else if(is.list(tau.hazard)){
-                tau.hazard_strata <- tau.hazard[[nStrata]]
-            }else{
-                tau.hazard_strata <- tau.hazard
             }
 
             ## E
@@ -418,14 +427,14 @@ iidCox.coxph <- function(object, newdata = NULL,
                 tau.hazard_strata <- 0
             }
             if(need.order){
-                out$time[[iStrata]] <- tau.hazard_strata[tau.oorder[[iStrata]]]
+                out$time[[iStrata]] <- tau.hazard_strata[iTau.oorder]
             }else{
                 out$time[[iStrata]] <- tau.hazard_strata
             }
             if(store.iid=="minimal"){
                 if(need.order && nVar.lp>0){
-                    out$calcIFhazard$Elambda0[[iStrata]] <- IFlambda_res$Elambda0[,tau.oorder[[iStrata]],drop=FALSE]
-                    out$calcIFhazard$cumElambda0[[iStrata]] <- IFlambda_res$cumElambda0[,tau.oorder[[iStrata]],drop=FALSE]
+                    out$calcIFhazard$Elambda0[[iStrata]] <- IFlambda_res$Elambda0[,iTau.oorder,drop=FALSE]
+                    out$calcIFhazard$cumElambda0[[iStrata]] <- IFlambda_res$cumElambda0[,iTau.oorder,drop=FALSE]
                 }else{
                     out$calcIFhazard$Elambda0[[iStrata]] <- IFlambda_res$Elambda0
                     out$calcIFhazard$cumElambda0[[iStrata]] <- IFlambda_res$cumElambda0
@@ -441,8 +450,8 @@ iidCox.coxph <- function(object, newdata = NULL,
                     colnames(IFlambda_res$cumhazard) <- tau.hazard_strata
                 } 
                 if(need.order){
-                    out$IFhazard[[iStrata]] <- IFlambda_res$hazard[,tau.oorder[[iStrata]],drop=FALSE]
-                    out$IFcumhazard[[iStrata]] <- IFlambda_res$cumhazard[,tau.oorder[[iStrata]],drop=FALSE]
+                    out$IFhazard[[iStrata]] <- IFlambda_res$hazard[,iTau.oorder,drop=FALSE]
+                    out$IFcumhazard[[iStrata]] <- IFlambda_res$cumhazard[,iTau.oorder,drop=FALSE]
                 }else{
                     out$IFhazard[[iStrata]] <- IFlambda_res$hazard
                     out$IFcumhazard[[iStrata]] <- IFlambda_res$cumhazard

R/predictCox.R

@@ -234,6 +234,10 @@ predictCox <- function(object,
         nTimes <- 0
         times <- numeric(0)
     }else{
+        if(!is.numeric(times) && !is.integer(times)){
+            stop("Argument \'times\' should be a numeric vector. \n",
+                 "Provide class: ",class(times)[1],". \n")
+        }
         nTimes <- length(times)
     }
     needOrder <- (nTimes[1]>0 && is.unsorted(times))

R/riskRegression-package.R

@@ -99,6 +99,6 @@ NULL
 #' @importFrom grDevices col2rgb gray
 #' @importFrom graphics bxp abline axis box legend lines mtext par plot points segments text title polygon par boxplot
 #' @importFrom utils capture.output find head select.list setTxtProgressBar tail txtProgressBar
-#' @importFrom stats confint cov as.formula coef delete.response drop.terms family formula get_all_vars lm glm median model.frame model.matrix model.response model.tables na.fail na.omit nobs optim pnorm predict qnorm quantile rbinom reformulate rexp runif sd setNames smooth terms terms.formula time uniroot update update.formula var weights wilcox.test
+#' @importFrom stats confint cov as.formula coef delete.response drop.terms family formula get_all_vars lm glm median model.frame model.matrix model.response model.tables na.fail na.omit nobs optim pnorm predict qnorm quantile rbinom reformulate rexp runif sd setNames smooth terms terms.formula time uniroot update update.formula var vcov weights wilcox.test
 "_PACKAGE"