diff --git a/DESCRIPTION b/DESCRIPTION index 6a1bb03..f285d57 100644 --- a/DESCRIPTION +++ b/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 = "tag@biostat.ku.dk"), person("Johan Sebastian", "Ohlendorff", role = "aut"), diff --git a/NAMESPACE b/NAMESPACE index 5a9e448..828350c 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS b/NEWS index d587596..2c0f663 100644 --- a/NEWS +++ b/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. diff --git a/R/RcppExports.R b/R/RcppExports.R index 9ab7985..fc567b4 100644 --- a/R/RcppExports.R +++ b/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) { diff --git a/R/calcSeCox.R b/R/calcSeCox.R index a489f36..af615ae 100644 --- a/R/calcSeCox.R +++ b/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, diff --git a/R/iidCox.R b/R/iidCox.R index cd907b7..1156176 100644 --- a/R/iidCox.R +++ b/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 diff --git a/R/predictCox.R b/R/predictCox.R index 82bad89..0aba6ee 100644 --- a/R/predictCox.R +++ b/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)) diff --git a/R/riskRegression-package.R b/R/riskRegression-package.R index 0c607b2..bfedba2 100644 --- a/R/riskRegression-package.R +++ b/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" diff --git a/R/wglm.R b/R/wglm.R index 83d2b5d..36d2261 100644 --- a/R/wglm.R +++ b/R/wglm.R @@ -3,9 +3,9 @@ ## Author: Brice Ozenne ## Created: sep 1 2020 (14:58) ## Version: -## Last-Updated: Oct 17 2024 (09:51) +## Last-Updated: Oct 20 2024 (13:59) ## By: Brice Ozenne -## Update #: 733 +## Update #: 774 ##---------------------------------------------------------------------- ## ### Commentary: @@ -29,7 +29,9 @@ #' @param ties [character] method used to handle ties when using a Cox model (\code{"breslow"} or \code{"efron"}). #' Ignored if fitter equals to \code{"prodlim"}. #' @param product.limit [logical] if \code{TRUE} the survival is computed using the product limit estimator. -#' @param store [vector of length 2] Whether prediction should only be computed for unique covariate sets and mapped back to the original dataset (\code{data="minimal"}) and whether the influence function should be stored in a memory efficient way (\code{iid="minimal"}). Otherwise use \code{data="full"} and/or \code{iid="full"}. +#' @param iid [logical] should the influence function of the logistic regression parameters be computed, accounting for the uncertainty of the weights. This can be computationally and memory intensive. +#' @param se [logical] should the variance-covariance matrix of the logistic regression parameters be stored, accounting for the uncertainty of the weights. This can be computationally and memory intensive. +#' @param store [vector] when evaluating the iid, should prediction be only computed for unique covariate sets and mapped back to the original dataset (\code{data="minimal"}). Otherwise use \code{data="full"}. #' #' @details First, a Cox model is fitted (argument formula.censor) #' and the censoring probabilities are computed relative to each timepoint (argument times) to obtain the censoring weights. @@ -82,7 +84,8 @@ #' @export wglm <- function(formula.event, times, data, formula.censor = ~1, cause = NA, - fitter = NULL, ties = NULL, product.limit = NULL, store = NULL){ + fitter = NULL, ties = NULL, product.limit = NULL, + iid = FALSE, se = TRUE, store = NULL){ tol <- 1e-12 mycall <- match.call() @@ -123,10 +126,12 @@ wglm <- function(formula.event, times, data, formula.censor = ~1, cause = NA, stop("Mismatch between argument \'formula.event\' and argument \'data\'. \n", "Could not find the variable(s) \'",paste(setdiff(c(varSurv$time,varSurv$status), names(data)), collapse = "\', \'"),"\' \n.") } + if(is.logical(data[[varSurv$status]])){ + data[[varSurv$status]] <- as.numeric(data[[varSurv$status]]) + } if(any(is.na(data))){ warning("Argument \'data\' contains missing values. \n") } - init.Hist <- Hist(time = data[[varSurv$time]], event = data[[varSurv$status]]) allStates <- c(attr(init.Hist, "cens.code"),attr(init.Hist, "states")) code.cens <- attr(init.Hist, "cens.code") @@ -262,7 +267,7 @@ wglm <- function(formula.event, times, data, formula.censor = ~1, cause = NA, iIndex <- which(data[[obs.newname[iTime]]]>0) if(n.censor[iTime]>0){ iPred <- predictRisk(object.censor, diag = TRUE, newdata = data[iIndex,,drop=FALSE], times = pmin(data.time[iIndex], times[iTime]) - tol, - type = "survival", product.limit = product.limit)[,1] + type = "survival", product.limit = product.limit, store = store)[,1] }else{ iPred <- rep(1, length(iIndex)) } @@ -295,6 +300,15 @@ wglm <- function(formula.event, times, data, formula.censor = ~1, cause = NA, out$product.limit <- product.limit out$fitter <- fitter class(out) <- append("wglm",class(out)) + if(iid || se){ + out$iid <- iid(out, simplify = FALSE) + } + if(se){ + out$vcov <- stats::vcov(out) + if(iid == FALSE){ + out$iid <- NULL + } + } return(out) } @@ -319,6 +333,13 @@ formula.wglm <- function(x, ...){ #' @param ... Not used. #' @export coef.wglm <- function(object, times = NULL, simplify = TRUE, ...){ + + ## ** normalize user input + dots <- list(...) + if(length(dots)>0){ + stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n") + } + if(is.null(times)){ times <- object$times }else{ @@ -330,12 +351,90 @@ coef.wglm <- function(object, times = NULL, simplify = TRUE, ...){ } n.times <- length(times) + ## ** extract coefficients M.coef <- do.call(rbind,setNames(lapply(object$fit, coef),object$time)) if(simplify && length(times)==1){ out <- stats::setNames(M.coef[object$times %in% times,], colnames(M.coef)) }else{ out <- M.coef[object$times %in% times,,drop=FALSE] } + + ## ** export + return(out) +} + +## * confint.wglm +#' @title Confidence intervals for Estimate from IPCW Logistic Regressions +#' @description Display the confidence intervals w.r.t. the estimated regression parameters from IPCW logistic regressions. +#' +#' @param object a wglm object. +#' @param parm not used. For compatibility with the generic method. +#' @param level [numeric, 0-1] Level of confidence. +#' @param times [numeric vector] time points at which the estimates should be output. +#' @param type [character] should the robust variance-covariance matrix be computing accounting for the uncertainty of the IPCW (\code{"robust"}) +#' or ignoring the uncertainty of the IPCW (\code{"robust-wknown"}), or model-based ignoring the uncertainty of the IPCW (\code{"model-wknown"})? +#' @param ... Not used. +#' @export +confint.wglm <- function(object, parm = NULL, level = 0.95, times = NULL, type = "robust", ...){ + + ## ** extract coefficients + dots <- list(...) + if(length(dots)>0){ + stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n") + } + + outTable <- model.tables(object, times = times, type = type, level = level) + + ## ** export + out <- outTable[,c("estimate","lower","upper")] + rownames(out) <- paste0(outTable$name,"(t=",outTable$time,")") + return(out) +} + +## * model.tables.wglm +#' @title Statistical Inference for Estimate from IPCW Logistic Regressions +#' @description Export estimated regression parameters from IPCW logistic regressions with their uncertainty (standard errors, confidence intervals and p-values). +#' +#' @param x a wglm object. +#' @param times [numeric vector] time points at which the estimates should be output. +#' @param type [character] should the robust variance-covariance matrix be computing accounting for the uncertainty of the IPCW (\code{"robust"}) +#' or ignoring the uncertainty of the IPCW (\code{"robust-wknown"}), or model-based ignoring the uncertainty of the IPCW (\code{"model-wknown"})? +#' @param level [numeric, 0-1] Level of confidence. +#' @param ... Not used. +#' @export +model.tables.wglm <- function(x, times = NULL, type = "robust", level = 0.95, ...){ + + ## ** normalize user input + dots <- list(...) + if(length(dots)>0){ + stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n") + } + + if(is.null(times)){ + times <- x$times + }else{ + if(any(times %in% x$times == FALSE)){ + stop("Incorrect specification of argument \'times\' \n", + "Should be one of \"",paste0(times,collapse="\" \""),"\" \n") + + } + } + n.times <- length(times) + + ## ** extract coefficients and se + Mcoef <- coef(x, times = times, simplify = FALSE) + Mse <- lapply(stats::vcov(x, times = times, type = type, simplify = FALSE), function(iM){sqrt(diag(iM))}) + ls.out <- lapply(1:length(Mse), function(iTime){data.frame(name = names(Mse[[iTime]]), time = times[iTime], estimate = Mcoef[iTime,], se = Mse[[iTime]])}) + out <- do.call(rbind,ls.out) + rownames(out) <- NULL + + ## ** evaluate ci and p.value + out$statistic <- out$estimate/out$se + out$lower <- out$estimate + qnorm((1-level)/2)*out$se + out$upper <- out$estimate + qnorm(1-(1-level)/2)*out$se + out$p.value <- 2*(1-pnorm(abs(out$statistic))) + + ## ** export return(out) } @@ -345,16 +444,53 @@ coef.wglm <- function(object, times = NULL, simplify = TRUE, ...){ #' #' @param object a wglm object. #' @param times [numeric vector] time points at which the variance-covariance matrix should be output. +#' @param type [character] should the robust variance-covariance matrix be computing accounting for the uncertainty of the IPCW (\code{"robust"}) +#' or ignoring the uncertainty of the IPCW (\code{"robust-wknown"}), or model-based ignoring the uncertainty of the IPCW (\code{"model-wknown"})? #' @param simplify [logical] should the ouput be converted to a matrix when only one timepoint is requested. Otherwise will always return a list. #' @param ... Not used. #' @export -vcov.wglm <- function(object, times = NULL, simplify = TRUE, ...){ - ls.iid <- iid(object, times = times, simplify = simplify) - if(is.list(ls.iid)){ - out <- lapply(ls.iid, crossprod) +vcov.wglm <- function(object, times = NULL, type = "robust", simplify = TRUE, ...){ + + ## ** normalize user input + dots <- list(...) + if(length(dots)>0){ + stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n") + } + + type <- match.arg(type, c("robust","robust-wknown","model-wknown")) + if(is.null(times)){ + times <- object$times }else{ - out <- crossprod(ls.iid) + if(any(times %in% object$times == FALSE)){ + stop("Incorrect specification of argument \'times\' \n", + "Should be one of \"",paste0(times,collapse="\" \""),"\" \n") + + } + } + + ## ** evaluate variance-covariance matrix + if(type == "robust"){ + if(!is.null(object$vcov)){ + out <- object$vcov[match(times, object$times)] + }else{ + if(!is.null(object$iid)){ + ls.iid <- object$iid[match(times, object$times)] + }else{ + ls.iid <- iid(object, times = times, simplify = simplify) + } + if(is.list(ls.iid)){ + out <- lapply(ls.iid, crossprod) + }else{ + out <- crossprod(ls.iid) + } + } + }else if(type == "model-wknown"){ + out <- lapply(object$fit[match(times, object$times)],stats::vcov) + }else if(type == "robust-wknown"){ + out <- lapply(object$fit[match(times, object$times)], function(iO){crossprod(lava::iid(iO))}) } + + ## ** export return(out) } @@ -373,11 +509,7 @@ summary.wglm <- function(object, print = TRUE, se = "robust", times = NULL, ...) } } n.times <- length(times) - if(se == "robust"){ - object.iid <- lava::iid(object, simplify = FALSE, times = times) - }else if(se == "robust-wknown"){ - object.iid <- lapply(object$fit[match(times, object$times)],lava::iid) - } + object.vcov <- stats::vcov(object, simplify = FALSE, times = times, type = se) out <- setNames(vector(mode = "list", length = n.times), times) if(print){ @@ -387,7 +519,7 @@ summary.wglm <- function(object, print = TRUE, se = "robust", times = NULL, ...) iTime2 <- which(object$times == times[iTime]) suppressWarnings(out[[iTime]] <- summary(object$fit[[iTime2]])$coef) if(se %in% c("robust","robust-wknown")){ - out[[iTime]][,"Std. Error"] <- sqrt(diag(crossprod(object.iid[[iTime]]))) + out[[iTime]][,"Std. Error"] <- sqrt(diag(object.vcov[[iTime]])) out[[iTime]][,3] <- out[[iTime]][,"Estimate"]/out[[iTime]][,"Std. Error"] ## name change from z to t stat when using quasibinomial instead binomial out[[iTime]][,4] <- 2*(1-pnorm(abs(out[[iTime]][,3]))) } @@ -412,6 +544,11 @@ summary.wglm <- function(object, print = TRUE, se = "robust", times = NULL, ...) print.wglm <- function(x, times = NULL, short = FALSE, ...){ ## ** prepare + dots <- list(...) + if(length(dots)>0){ + stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n") + } + if(is.null(times)){ times <- x$times }else{ @@ -476,6 +613,13 @@ print.wglm <- function(x, times = NULL, short = FALSE, ...){ #' @param ... Not used. #' @export score.wglm <- function(x, indiv = FALSE, times = NULL, simplify = TRUE, ...){ + + ## ** noramlize user input + dots <- list(...) + if(length(dots)>0){ + stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n") + } + if(inherits(x,"glm") && !inherits(x,"wglm")){ x <- list(fit = list("1" = x), times = "1") @@ -491,6 +635,7 @@ score.wglm <- function(x, indiv = FALSE, times = NULL, simplify = TRUE, ...){ } n.times <- length(times) + ## ** compute score out <- setNames(vector(mode = "list", length = n.times), times) for(iTime in 1:n.times){ iTime2 <- which(x$times == times[iTime]) @@ -505,6 +650,7 @@ score.wglm <- function(x, indiv = FALSE, times = NULL, simplify = TRUE, ...){ if(indiv==FALSE){out[[iTime]] <- colSums(out[[iTime]])} } + ## ** export if(n.times == 1 && simplify){ return(out[[1]]) }else{ @@ -522,6 +668,13 @@ score.wglm <- function(x, indiv = FALSE, times = NULL, simplify = TRUE, ...){ #' @param ... Not used. #' @export information.wglm <- function(x, times = NULL, simplify = TRUE, ...){ + + ## ** normalize user input + dots <- list(...) + if(length(dots)>0){ + stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n") + } + if(inherits(x,"glm") && !inherits(x,"wglm")){ x <- list(fit = list("1" = x), times = "1") @@ -537,6 +690,7 @@ information.wglm <- function(x, times = NULL, simplify = TRUE, ...){ } n.times <- length(times) + ## ** compute information out <- setNames(vector(mode = "list", length = n.times), times) for(iTime in 1:n.times){ iTime2 <- which(x$times == times[iTime]) @@ -549,6 +703,7 @@ information.wglm <- function(x, times = NULL, simplify = TRUE, ...){ rownames(out[[iTime]]) <- colnames(out[[iTime]]) } + ## ** export if(n.times == 1 && simplify){ return(out[[1]]) }else{ @@ -564,9 +719,17 @@ information.wglm <- function(x, times = NULL, simplify = TRUE, ...){ #' @param x a wglm object. #' @param times [numeric vector] time points at which the iid should be output. #' @param simplify [logical] should the ouput be converted to a matrix when only one timepoint is requested. Otherwise will always return a list. +#' @param store [vector] when evaluating the iid, should prediction be only computed for unique covariate sets and mapped back to the original dataset (\code{data="minimal"}). Otherwise use \code{data="full"}. #' @param ... Not used. #' @export -iid.wglm <- function(x, times = NULL, simplify = TRUE, ...){ +iid.wglm <- function(x, times = NULL, simplify = TRUE, store = NULL, ...){ + + ## ** normalize user input + dots <- list(...) + if(length(dots)>0){ + stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n") + } + if(inherits(x,"glm") && !inherits(x,"wglm")){ x <- list(fit = list("1" = x), times = "1", @@ -582,16 +745,19 @@ iid.wglm <- function(x, times = NULL, simplify = TRUE, ...){ } } + + ## ** prepare (glm score and information) n.times <- length(times) ls.score <- lava::score(x, times = times, simplify = FALSE, indiv = TRUE) ls.info <- lava::information(x, times = times, simplify = FALSE) + ## ** evaluate iid out <- setNames(vector(mode = "list", length = n.times), times) - for(iTime in 1:n.times){ - iTime2 <- which(x$times == times[iTime]) - iObject <- x$fit[[iTime2]] + for(iT in 1:n.times){ ## iT <- 1 + iT2 <- which(x$times == times[iT]) + iObject <- x$fit[[iT2]] - ## ** compute the uncertainty related to the weights + ## *** compute the uncertainty related to the weights ## S score, I information matrix, H hessian, X design matrix, Y outcome, \pi fitted probabilities ## \beta regresssion parameters (logit) \eta nuisance parameters (censoring) ## The estimating equation of \beta is @@ -611,7 +777,7 @@ iid.wglm <- function(x, times = NULL, simplify = TRUE, ...){ ## = (S_i + n*AIF_{W,X*(Y-pi)}(O_i) / I ## NOTE: IF_{W_O}(O_i) = dW(O,\eta)/d\eta \IF_\eta(O_i) - if(x$n.censor[iTime]>0){ + if(x$n.censor[iT2]>0){ n.obs <- sum(coxN(iObject)) W2 <- iObject$prior.weights2 X <- stats::model.matrix(iObject) @@ -621,25 +787,25 @@ iid.wglm <- function(x, times = NULL, simplify = TRUE, ...){ attr(factor,"factor") <- lapply(apply(colMultiply_cpp(X, -W2*(Y - pi)), 2, list), function(iVec){cbind(iVec[[1]])}) iPred <- predictRisk(x$model.censor, diag = TRUE, newdata = x$data, times = iObject$time.prior.weights, - type = "survival", product.limit = x$product.limit, average.iid = factor) + type = "survival", product.limit = x$product.limit, average.iid = factor, store = store) } - ## ** assemble uncertainty + ## *** assemble uncertainty ## (S+dS/dW)/I - if(x$n.censor[iTime]>0){ - out[[iTime]] <- (ls.score[[iTime]] + do.call(cbind,attr(iPred, "average.iid"))*NROW(x$data)) %*% solve(ls.info[[iTime]]) + if(x$n.censor[iT2]>0){ + out[[iT]] <- (ls.score[[iT]] + do.call(cbind,attr(iPred, "average.iid"))*NROW(x$data)) %*% solve(ls.info[[iT]]) }else{ - out[[iTime]] <- ls.score[[iTime]] %*% solve(ls.info[[iTime]]) + out[[iT]] <- ls.score[[iT]] %*% solve(ls.info[[iT]]) } ## ## WRONG VERSION ## factor1 <- colMultiply_cpp(X, (Y - pi)) %*% iVcov - ## factor2 <- - ls.score[[iTime]] %*% crossprod(iVcov) %*% t(colMultiply_cpp(X, scale = -pi*(1-pi))) %*% X + ## factor2 <- - ls.score[[iT]] %*% crossprod(iVcov) %*% t(colMultiply_cpp(X, scale = -pi*(1-pi))) %*% X ## factor <- TRUE ## attr(factor,"factor") <- lapply(apply(factor1 + factor2, 2, list), function(iVec){cbind(iVec[[1]])}) ## iPred <- predictRisk(x$model.censor, diag = TRUE, newdata = x$data, times = iObject$time.prior.weights, ## type = "survival", product.limit = FALSE, average.iid = factor) - ## out[[iTime]] <- out[[iTime]] - do.call(cbind,attr(iPred, "average.iid"))*NROW(x$data) + ## out[[iT]] <- out[[iT]] - do.call(cbind,attr(iPred, "average.iid"))*NROW(x$data) } @@ -664,15 +830,25 @@ iid.wglm <- function(x, times = NULL, simplify = TRUE, ...){ #' weights.wglm <- function(object, times = NULL, simplify = TRUE, ...){ - out <- object$data[object$name.IPCW] + ## ** normalize user input + dots <- list(...) + if(length(dots)>0){ + stop("Unknown argument(s) \'",paste(names(dots),collapse="\' \'"),"\'. \n") + } if(!is.null(times)){ if(any(times %in% object$time == FALSE)){ stop("Unknown timepoint ",paste(setdiff(times,object$time), collapse = ", ")," in argument \'times\'. \n", "Valid timepoints: ",paste(object$time, collapse = ", "),". \n") } + } + + ## ** extract weights + out <- object$data[object$name.IPCW] + if(!is.null(times)){ out <- out[match(times,object$time)] } + ## ** export if(simplify && NCOL(out)==1){ return(out[[1]]) }else{ diff --git a/man/confint.wglm.Rd b/man/confint.wglm.Rd new file mode 100644 index 0000000..d403e54 --- /dev/null +++ b/man/confint.wglm.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/wglm.R +\name{confint.wglm} +\alias{confint.wglm} +\title{Confidence intervals for Estimate from IPCW Logistic Regressions} +\usage{ +\method{confint}{wglm}(object, parm = NULL, level = 0.95, times = NULL, type = "robust", ...) +} +\arguments{ +\item{object}{a wglm object.} + +\item{parm}{not used. For compatibility with the generic method.} + +\item{level}{[numeric, 0-1] Level of confidence.} + +\item{times}{[numeric vector] time points at which the estimates should be output.} + +\item{type}{[character] should the robust variance-covariance matrix be computing accounting for the uncertainty of the IPCW (\code{"robust"}) +or ignoring the uncertainty of the IPCW (\code{"robust-wknown"}), or model-based ignoring the uncertainty of the IPCW (\code{"model-wknown"})?} + +\item{...}{Not used.} +} +\description{ +Display the confidence intervals w.r.t. the estimated regression parameters from IPCW logistic regressions. +} diff --git a/man/iid.wglm.Rd b/man/iid.wglm.Rd index 72cc5e5..b321c2c 100644 --- a/man/iid.wglm.Rd +++ b/man/iid.wglm.Rd @@ -4,7 +4,7 @@ \alias{iid.wglm} \title{IID for IPCW Logistic Regressions} \usage{ -\method{iid}{wglm}(x, times = NULL, simplify = TRUE, ...) +\method{iid}{wglm}(x, times = NULL, simplify = TRUE, store = NULL, ...) } \arguments{ \item{x}{a wglm object.} @@ -13,6 +13,8 @@ \item{simplify}{[logical] should the ouput be converted to a matrix when only one timepoint is requested. Otherwise will always return a list.} +\item{store}{[vector] when evaluating the iid, should prediction be only computed for unique covariate sets and mapped back to the original dataset (\code{data="minimal"}). Otherwise use \code{data="full"}.} + \item{...}{Not used.} } \description{ diff --git a/man/model.tables.wglm.Rd b/man/model.tables.wglm.Rd new file mode 100644 index 0000000..387d050 --- /dev/null +++ b/man/model.tables.wglm.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/wglm.R +\name{model.tables.wglm} +\alias{model.tables.wglm} +\title{Statistical Inference for Estimate from IPCW Logistic Regressions} +\usage{ +\method{model.tables}{wglm}(x, times = NULL, type = "robust", level = 0.95, ...) +} +\arguments{ +\item{x}{a wglm object.} + +\item{times}{[numeric vector] time points at which the estimates should be output.} + +\item{type}{[character] should the robust variance-covariance matrix be computing accounting for the uncertainty of the IPCW (\code{"robust"}) +or ignoring the uncertainty of the IPCW (\code{"robust-wknown"}), or model-based ignoring the uncertainty of the IPCW (\code{"model-wknown"})?} + +\item{level}{[numeric, 0-1] Level of confidence.} + +\item{...}{Not used.} +} +\description{ +Export estimated regression parameters from IPCW logistic regressions with their uncertainty (standard errors, confidence intervals and p-values). +} diff --git a/man/vcov.wglm.Rd b/man/vcov.wglm.Rd index 0488cf9..dfd8417 100644 --- a/man/vcov.wglm.Rd +++ b/man/vcov.wglm.Rd @@ -4,13 +4,16 @@ \alias{vcov.wglm} \title{Variance-covariance for IPCW Logistic Regressions} \usage{ -\method{vcov}{wglm}(object, times = NULL, simplify = TRUE, ...) +\method{vcov}{wglm}(object, times = NULL, type = "robust", simplify = TRUE, ...) } \arguments{ \item{object}{a wglm object.} \item{times}{[numeric vector] time points at which the variance-covariance matrix should be output.} +\item{type}{[character] should the robust variance-covariance matrix be computing accounting for the uncertainty of the IPCW (\code{"robust"}) +or ignoring the uncertainty of the IPCW (\code{"robust-wknown"}), or model-based ignoring the uncertainty of the IPCW (\code{"model-wknown"})?} + \item{simplify}{[logical] should the ouput be converted to a matrix when only one timepoint is requested. Otherwise will always return a list.} \item{...}{Not used.} diff --git a/man/wglm.Rd b/man/wglm.Rd index 807f118..7a2e4c0 100644 --- a/man/wglm.Rd +++ b/man/wglm.Rd @@ -13,6 +13,8 @@ wglm( fitter = NULL, ties = NULL, product.limit = NULL, + iid = FALSE, + se = TRUE, store = NULL ) } @@ -34,7 +36,11 @@ Ignored if fitter equals to \code{"prodlim"}.} \item{product.limit}{[logical] if \code{TRUE} the survival is computed using the product limit estimator.} -\item{store}{[vector of length 2] Whether prediction should only be computed for unique covariate sets and mapped back to the original dataset (\code{data="minimal"}) and whether the influence function should be stored in a memory efficient way (\code{iid="minimal"}). Otherwise use \code{data="full"} and/or \code{iid="full"}.} +\item{iid}{[logical] should the influence function of the logistic regression parameters be computed, accounting for the uncertainty of the weights. This can be computationally and memory intensive.} + +\item{se}{[logical] should the variance-covariance matrix of the logistic regression parameters be stored, accounting for the uncertainty of the weights. This can be computationally and memory intensive.} + +\item{store}{[vector] when evaluating the iid, should prediction be only computed for unique covariate sets and mapped back to the original dataset (\code{data="minimal"}). Otherwise use \code{data="full"}.} } \value{ an object of class \code{"wglm"}. diff --git a/slowtests/test-predictCSC.R b/slowtests/test-predictCSC.R index 020c0ef..cad188c 100644 --- a/slowtests/test-predictCSC.R +++ b/slowtests/test-predictCSC.R @@ -3,9 +3,9 @@ ## author: Brice Ozenne ## created: maj 18 2017 (09:23) ## Version: -## last-updated: Oct 17 2024 (10:21) +## last-updated: Oct 20 2024 (21:33) ## By: Brice Ozenne -## Update #: 350 +## Update #: 351 #---------------------------------------------------------------------- ## ### Commentary: @@ -463,7 +463,7 @@ test_that("predict.CSC (no covariates): compare to mstate",{ }) ## ** With covariates test_that("predict.CSC (covariates): compare to mstate",{ - for(iX in 0:1){ + for(iX in 0:1){ ## iX <- 0 newdata <- data.frame(X1 = iX, X2 = 0, X16 = 0) newdata.L <- data.frame(X1.1 = c(iX, 0), X1.2 = c(0, iX), X2.1 = c(0, 0), X2.2 = c(0, 0), diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 45ae154..0b595c0 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -95,7 +95,7 @@ BEGIN_RCPP END_RCPP } // calcSeCif2_cpp -List calcSeCif2_cpp(const std::vector& ls_IFbeta, const std::vector& ls_X, const std::vector& ls_cumhazard, const arma::mat& ls_hazard, const arma::mat& survival, const arma::mat& cif, const std::vector< std::vector >& ls_IFcumhazard, const std::vector& ls_IFhazard, const NumericMatrix& eXb, int nJumpTime, const NumericVector& JumpMax, const NumericVector& tau, const arma::vec& tauIndex, int nTau, int nObs, int theCause, int nCause, bool hazardType, arma::vec nVar, int nNewObs, NumericMatrix strata, bool exportSE, bool exportIF, bool exportIFsum, bool diag); +List calcSeCif2_cpp(const std::vector& ls_IFbeta, const std::vector& ls_X, const std::vector& ls_cumhazard, const arma::mat& ls_hazard, const arma::mat& survival, const arma::mat& cif, const std::vector< std::vector >& ls_IFcumhazard, const std::vector& ls_IFhazard, const arma::mat& eXb, int nJumpTime, const NumericVector& JumpMax, const NumericVector& tau, const arma::vec& tauIndex, int nTau, int nObs, int theCause, int nCause, bool hazardType, arma::vec nVar, int nNewObs, arma::mat strata, bool exportSE, bool exportIF, bool exportIFsum, bool diag); RcppExport SEXP _riskRegression_calcSeCif2_cpp(SEXP ls_IFbetaSEXP, SEXP ls_XSEXP, SEXP ls_cumhazardSEXP, SEXP ls_hazardSEXP, SEXP survivalSEXP, SEXP cifSEXP, SEXP ls_IFcumhazardSEXP, SEXP ls_IFhazardSEXP, SEXP eXbSEXP, SEXP nJumpTimeSEXP, SEXP JumpMaxSEXP, SEXP tauSEXP, SEXP tauIndexSEXP, SEXP nTauSEXP, SEXP nObsSEXP, SEXP theCauseSEXP, SEXP nCauseSEXP, SEXP hazardTypeSEXP, SEXP nVarSEXP, SEXP nNewObsSEXP, SEXP strataSEXP, SEXP exportSESEXP, SEXP exportIFSEXP, SEXP exportIFsumSEXP, SEXP diagSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -108,7 +108,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::mat& >::type cif(cifSEXP); Rcpp::traits::input_parameter< const std::vector< std::vector >& >::type ls_IFcumhazard(ls_IFcumhazardSEXP); Rcpp::traits::input_parameter< const std::vector& >::type ls_IFhazard(ls_IFhazardSEXP); - Rcpp::traits::input_parameter< const NumericMatrix& >::type eXb(eXbSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type eXb(eXbSEXP); Rcpp::traits::input_parameter< int >::type nJumpTime(nJumpTimeSEXP); Rcpp::traits::input_parameter< const NumericVector& >::type JumpMax(JumpMaxSEXP); Rcpp::traits::input_parameter< const NumericVector& >::type tau(tauSEXP); @@ -120,7 +120,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< bool >::type hazardType(hazardTypeSEXP); Rcpp::traits::input_parameter< arma::vec >::type nVar(nVarSEXP); Rcpp::traits::input_parameter< int >::type nNewObs(nNewObsSEXP); - Rcpp::traits::input_parameter< NumericMatrix >::type strata(strataSEXP); + Rcpp::traits::input_parameter< arma::mat >::type strata(strataSEXP); Rcpp::traits::input_parameter< bool >::type exportSE(exportSESEXP); Rcpp::traits::input_parameter< bool >::type exportIF(exportIFSEXP); Rcpp::traits::input_parameter< bool >::type exportIFsum(exportIFsumSEXP); @@ -173,8 +173,8 @@ BEGIN_RCPP END_RCPP } // calcAIFsurv_cpp -std::vector< std::vector > calcAIFsurv_cpp(const std::vector& ls_IFcumhazard, const arma::mat& IFbeta, const std::vector& cumhazard0, const arma::mat& survival, const arma::colvec& eXb, const arma::mat& X, const NumericVector& prevStrata, const std::vector& ls_indexStrata, const std::vector& factor, int nTimes, int nObs, int nStrata, int nVar, int diag, bool exportCumHazard, bool exportSurvival); -RcppExport SEXP _riskRegression_calcAIFsurv_cpp(SEXP ls_IFcumhazardSEXP, SEXP IFbetaSEXP, SEXP cumhazard0SEXP, SEXP survivalSEXP, SEXP eXbSEXP, SEXP XSEXP, SEXP prevStrataSEXP, SEXP ls_indexStrataSEXP, SEXP factorSEXP, SEXP nTimesSEXP, SEXP nObsSEXP, SEXP nStrataSEXP, SEXP nVarSEXP, SEXP diagSEXP, SEXP exportCumHazardSEXP, SEXP exportSurvivalSEXP) { +std::vector< std::vector > calcAIFsurv_cpp(const std::vector& ls_IFcumhazard, const arma::mat& IFbeta, const std::vector& cumhazard0, const arma::mat& survival, const arma::colvec& eXb, const arma::mat& X, const NumericVector& prevStrata, const std::vector& ls_indexStrata, const std::vector& ls_indexStrataTime, const std::vector& factor, int nTimes, int nObs, int nStrata, int nVar, int diag, bool exportCumHazard, bool exportSurvival); +RcppExport SEXP _riskRegression_calcAIFsurv_cpp(SEXP ls_IFcumhazardSEXP, SEXP IFbetaSEXP, SEXP cumhazard0SEXP, SEXP survivalSEXP, SEXP eXbSEXP, SEXP XSEXP, SEXP prevStrataSEXP, SEXP ls_indexStrataSEXP, SEXP ls_indexStrataTimeSEXP, SEXP factorSEXP, SEXP nTimesSEXP, SEXP nObsSEXP, SEXP nStrataSEXP, SEXP nVarSEXP, SEXP diagSEXP, SEXP exportCumHazardSEXP, SEXP exportSurvivalSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -186,6 +186,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); Rcpp::traits::input_parameter< const NumericVector& >::type prevStrata(prevStrataSEXP); Rcpp::traits::input_parameter< const std::vector& >::type ls_indexStrata(ls_indexStrataSEXP); + Rcpp::traits::input_parameter< const std::vector& >::type ls_indexStrataTime(ls_indexStrataTimeSEXP); Rcpp::traits::input_parameter< const std::vector& >::type factor(factorSEXP); Rcpp::traits::input_parameter< int >::type nTimes(nTimesSEXP); Rcpp::traits::input_parameter< int >::type nObs(nObsSEXP); @@ -194,7 +195,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type diag(diagSEXP); Rcpp::traits::input_parameter< bool >::type exportCumHazard(exportCumHazardSEXP); Rcpp::traits::input_parameter< bool >::type exportSurvival(exportSurvivalSEXP); - rcpp_result_gen = Rcpp::wrap(calcAIFsurv_cpp(ls_IFcumhazard, IFbeta, cumhazard0, survival, eXb, X, prevStrata, ls_indexStrata, factor, nTimes, nObs, nStrata, nVar, diag, exportCumHazard, exportSurvival)); + rcpp_result_gen = Rcpp::wrap(calcAIFsurv_cpp(ls_IFcumhazard, IFbeta, cumhazard0, survival, eXb, X, prevStrata, ls_indexStrata, ls_indexStrataTime, factor, nTimes, nObs, nStrata, nVar, diag, exportCumHazard, exportSurvival)); return rcpp_result_gen; END_RCPP } @@ -552,7 +553,7 @@ static const R_CallMethodDef CallEntries[] = { {"_riskRegression_calcSeMinimalCSC_cpp", (DL_FUNC) &_riskRegression_calcSeMinimalCSC_cpp, 36}, {"_riskRegression_calcSeCif2_cpp", (DL_FUNC) &_riskRegression_calcSeCif2_cpp, 25}, {"_riskRegression_calcSeMinimalCox_cpp", (DL_FUNC) &_riskRegression_calcSeMinimalCox_cpp, 33}, - {"_riskRegression_calcAIFsurv_cpp", (DL_FUNC) &_riskRegression_calcAIFsurv_cpp, 16}, + {"_riskRegression_calcAIFsurv_cpp", (DL_FUNC) &_riskRegression_calcAIFsurv_cpp, 17}, {"_riskRegression_calculateDelongCovarianceFast", (DL_FUNC) &_riskRegression_calculateDelongCovarianceFast, 2}, {"_riskRegression_colCumSum", (DL_FUNC) &_riskRegression_colCumSum, 1}, {"_riskRegression_quantileProcess_cpp", (DL_FUNC) &_riskRegression_quantileProcess_cpp, 7}, diff --git a/src/calcSeCSC.cpp b/src/calcSeCSC.cpp index 5ccd7a1..e764b9a 100644 --- a/src/calcSeCSC.cpp +++ b/src/calcSeCSC.cpp @@ -348,14 +348,14 @@ List calcSeMinimalCSC_cpp(const arma::vec& seqTau, // horizon time for the predi // * calcSeCif2_cpp: compute IF for the absolute risk (method 2) // [[Rcpp::export]] List calcSeCif2_cpp(const std::vector& ls_IFbeta, const std::vector& ls_X, - const std::vector& ls_cumhazard, const arma::mat& ls_hazard, const arma::mat& survival, const arma::mat& cif, - const std::vector< std::vector >& ls_IFcumhazard, const std::vector& ls_IFhazard, - const NumericMatrix& eXb, - int nJumpTime, const NumericVector& JumpMax, + const std::vector& ls_cumhazard, const arma::mat& ls_hazard, const arma::mat& survival, const arma::mat& cif, + const std::vector< std::vector >& ls_IFcumhazard, const std::vector& ls_IFhazard, + const arma::mat& eXb, + int nJumpTime, const NumericVector& JumpMax, const NumericVector& tau, const arma::vec& tauIndex, int nTau, int nObs, int theCause, int nCause, bool hazardType, arma::vec nVar, - int nNewObs, NumericMatrix strata, + int nNewObs, arma::mat strata, bool exportSE, bool exportIF, bool exportIFsum, bool diag){ arma::mat X_IFbeta; @@ -460,51 +460,51 @@ List calcSeCif2_cpp(const std::vector& ls_IFbeta, const std::vector0){ - X_IFbeta = ls_IFbeta[iCause] * (ls_X[iCause].row(iNewObs)).t(); - ieXb = eXb(iNewObs,iCause); + X_IFbeta = ls_IFbeta[iCause] * (ls_X[iCause].row(iNewObs)).t(); + ieXb = eXb(iNewObs,iCause); } // Rcout << "3 "; if(hazardType || (iCause != theCause)){ - if(nVar[iCause]>0){ - if(diag){ - IFcumhazard += ieXb * (ls_IFcumhazard[iCause][iStrataCause].cols(iUvec_linspace) + X_IFbeta * ls_tcumhazard[iCause].submat(iUvec_strata,iUvec_linspace)); - }else{ - IFcumhazard += ieXb * (ls_IFcumhazard[iCause][iStrataCause] + X_IFbeta * ls_tcumhazard[iCause].row(iStrataCause)); - } - }else{ - if(diag){ - IFcumhazard += ls_IFcumhazard[iCause][iStrataCause].cols(iUvec_linspace); - }else{ - IFcumhazard += ls_IFcumhazard[iCause][iStrataCause]; - } - } + if(nVar[iCause]>0){ + if(diag){ + IFcumhazard += ieXb * (ls_IFcumhazard[iCause][iStrataCause].cols(iUvec_linspace) + X_IFbeta * ls_tcumhazard[iCause].submat(iUvec_strata,iUvec_linspace)); + }else{ + IFcumhazard += ieXb * (ls_IFcumhazard[iCause][iStrataCause] + X_IFbeta * ls_tcumhazard[iCause].row(iStrataCause)); + } + }else{ + if(diag){ + IFcumhazard += ls_IFcumhazard[iCause][iStrataCause].cols(iUvec_linspace); + }else{ + IFcumhazard += ls_IFcumhazard[iCause][iStrataCause]; + } + } } // Rcout << "4 "; if(iCause == theCause){ - if(nVar[iCause] > 0){ - if(diag){ - hazard = ieXb * ls_hazard.submat(iUvec_linspace,iUvec_strata); - IFhazard = ieXb * (ls_IFhazard[iStrataCause].cols(iUvec_linspace) + X_IFbeta * ls_thazard.submat(iUvec_strata,iUvec_linspace)); - }else{ - hazard = ieXb * ls_hazard.col(iStrataCause); - IFhazard = ieXb * (ls_IFhazard[iStrataCause] + X_IFbeta * ls_thazard.row(iStrataCause)); - } - }else{ - if(diag){ - hazard = ls_hazard.submat(iUvec_linspace,iUvec_strata); - IFhazard = ls_IFhazard[iStrataCause].cols(iUvec_linspace); - }else{ - hazard = ls_hazard.col(iStrataCause); - IFhazard = ls_IFhazard[iStrataCause]; - } - } + if(nVar[iCause] > 0){ + if(diag){ + hazard = ieXb * ls_hazard.submat(iUvec_linspace,iUvec_strata); + IFhazard = ieXb * (ls_IFhazard[iStrataCause].cols(iUvec_linspace) + X_IFbeta * ls_thazard.submat(iUvec_strata,iUvec_linspace)); + }else{ + hazard = ieXb * ls_hazard.col(iStrataCause); + IFhazard = ieXb * (ls_IFhazard[iStrataCause] + X_IFbeta * ls_thazard.row(iStrataCause)); + } + }else{ + if(diag){ + hazard = ls_hazard.submat(iUvec_linspace,iUvec_strata); + IFhazard = ls_IFhazard[iStrataCause].cols(iUvec_linspace); + }else{ + hazard = ls_hazard.col(iStrataCause); + IFhazard = ls_IFhazard[iStrataCause]; + } + } } } @@ -516,13 +516,13 @@ List calcSeCif2_cpp(const std::vector& ls_IFbeta, const std::vector0){ - if(iJump==0){ - IF_tempo = IFhazard.col(iJump); - }else{ - // survival is evaluated just before the jump - IF_tempo = (IFhazard.col(iJump) - IFcumhazard.col(iJump-1) * hazard[iJump]) * survival(iNewObs,iJump); - } - cumIF_tempo = cumIF_tempo + IF_tempo; + if(iJump==0){ + IF_tempo = IFhazard.col(iJump); + }else{ + // survival is evaluated just before the jump + IF_tempo = (IFhazard.col(iJump) - IFcumhazard.col(iJump-1) * hazard[iJump]) * survival(iNewObs,iJump); + } + cumIF_tempo = cumIF_tempo + IF_tempo; } @@ -531,32 +531,32 @@ List calcSeCif2_cpp(const std::vector& ls_IFbeta, const std::vector& ls_IFbeta, const std::vector& ls_IFbeta, const std::vector > calcAIFsurv_cpp(const std::vector& ls_IFcumhazard, const arma::mat& IFbeta, @@ -367,6 +370,7 @@ std::vector< std::vector > calcAIFsurv_cpp(const std::vector& ls_indexStrata, + const std::vector& ls_indexStrataTime, const std::vector& factor, int nTimes, int nObs, @@ -435,13 +439,13 @@ std::vector< std::vector > calcAIFsurv_cpp(const std::vector