Skip to content

Commit

Permalink
... restaure tests (!) fixing slow tests ...
Browse files Browse the repository at this point in the history
  • Loading branch information
@bozenne
bozenne committed Oct 14, 2024
1 parent bb3088b commit fe6aa16
Show file tree
Hide file tree
Showing 37 changed files with 7,395 additions and 117 deletions.
4 changes: 2 additions & 2 deletions R/ate-iid.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,9 @@
## Author: Brice Ozenne
## Created: apr 5 2018 (17:01)
## Version:
## Last-Updated: sep 10 2024 (09:54)
## Last-Updated: Oct 14 2024 (10:01)
## By: Brice Ozenne
## Update #: 1341
## Update #: 1352
##----------------------------------------------------------------------
##
### Commentary:
Expand Down
108 changes: 53 additions & 55 deletions R/ate-pointEstimate.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,9 @@
## Author: Brice Ozenne
## Created: jun 27 2019 (10:43)
## Version:
## Last-Updated: sep 9 2024 (17:15)
## Last-Updated: Oct 14 2024 (09:49)
## By: Brice Ozenne
## Update #: 993
## Update #: 1015
##----------------------------------------------------------------------
##
### Commentary:
Expand Down Expand Up @@ -112,6 +112,7 @@ ATE_TI <- function(object.event,
index.lastjumpC <- max(index.obsSINDEXjumpC)
time.jumpC <- time.jumpC[1:index.lastjumpC]
}

if(attr(estimator,"integral")){
## *** jump time index of the event time stopped at tau
index.obsSINDEXjumpC.int <- do.call(cbind,lapply(times, function(tau){
Expand Down Expand Up @@ -214,58 +215,58 @@ ATE_TI <- function(object.event,

## ** Compute augmentation term
if(attr(estimator,"integral")){
augTerm <- matrix(0, nrow = n.obs, ncol = n.times)

## integral is sum over individuals from 0 to min(T_i,\tau) of dM_i^C. It is therefore non-0 only if:
## - some of the requested times are after the first censoring time
## - some of the censoring jump times are before the event times (always the case?)
if(index.lastjumpC>0 && any(beforeEvent.jumpC)){

## absolute risk at event times
predTempo <- predictRisk(object.event, newdata = mydata, times = c(times, time.jumpC), cause = cause, product.limit = product.limit,
iid = (method.iid==2)*return.iid.nuisance)
F1.tau <- predTempo[,1:n.times,drop=FALSE]
F1.jump <- predTempo[,n.times + (1:index.lastjumpC),drop=FALSE]
if((method.iid==2)*return.iid.nuisance){
out$store$iid.nuisance.outcome <- attr(predTempo,"iid")
}
## - some of the censoring jump times are before the event times
## this would correspond to index.lastjumpC>0 && any(beforeEvent.jumpC)
## the first should be automatically enforced in ate_initArgs when creating attr(estimator,"integral")

## survival at all jump of the censoring process
S.jump <- predictRisk(object.event, type = "survival", newdata = mydata, times = time.jumpC-tol, product.limit = product.limit,
iid = (method.iid==2)*return.iid.nuisance)
if((method.iid==2)*return.iid.nuisance){
out$store$iid.nuisance.survival <- attr(S.jump,"iid")
attr(S.jump,"iid") <- NULL
}
augTerm <- matrix(0, nrow = n.obs, ncol = n.times)

## martingale for the censoring process
## at all times of jump of the censoring process
G.jump <- 1-predictRisk(object.censor, newdata = mydata, times = if(index.lastjumpC>1){c(0,time.jumpC[1:(index.lastjumpC-1)])}else{0},
product.limit = product.limit, iid = (method.iid==2)*return.iid.nuisance)
## absolute risk at event times and jump times of the censoring process
predTempo <- predictRisk(object.event, newdata = mydata, times = c(times, time.jumpC), cause = cause, product.limit = product.limit,
iid = (method.iid==2)*return.iid.nuisance)
F1.tau <- predTempo[,1:n.times,drop=FALSE]
F1.jump <- predTempo[,n.times + (1:index.lastjumpC),drop=FALSE]
if((method.iid==2)*return.iid.nuisance){
out$store$iid.nuisance.outcome <- attr(predTempo,"iid")
}

if(return.iid.nuisance && (method.iid==2)){
out$store$iid.nuisance.censoring <- -attr(G.jump,"iid")
attr(G.jump,"iid") <- NULL
}
dN.jump <- do.call(cbind,lapply(time.jumpC, function(iJump){(mydata[[eventVar.time]] == iJump)*(mydata[[eventVar.status]] == level.censoring)}))
dLambda.jump <- predictCox(object.censor, newdata = mydata, times = time.jumpC, type = "hazard", iid = (method.iid==2)*return.iid.nuisance)
if((method.iid==2)*return.iid.nuisance){
out$store$iid.nuisance.martingale <- dLambda.jump$hazard.iid
}
dM.jump <- dN.jump - dLambda.jump$hazard
## integral = \int_0^min(T_i,\tau) (F1(\tau|A_i,W_i)-F1(t|A_i,W_i)) / S(t|A_i,W_i) * dM_i^C(t)/Gc(t|A_i,W_i)
## = F1(\tau|A_i,W_i) \int_0^min(T_i,\tauf) dM_i^C(t) / (S(t|A_i,W_i) * Gc(t|A_i,W_i)) - \int_0^min(T_i,\tauf) F1(t|A_i,W_i) dM_i^C(t) / (S(t|A_i,W_i) * Gc(t|A_i,W_i))
integrand <- matrix(0, nrow = n.obs, ncol = index.lastjumpC)
integrand2 <- matrix(0, nrow = n.obs, ncol = index.lastjumpC)
index.beforeEvent.jumpC <- which(beforeEvent.jumpC) ## identify which increment to put in the integral, i.e. when the individual was still at risk of being censored
integrand[index.beforeEvent.jumpC] <- dM.jump[index.beforeEvent.jumpC] / (G.jump[index.beforeEvent.jumpC] * S.jump[index.beforeEvent.jumpC])
integrand2[index.beforeEvent.jumpC] <- F1.jump[index.beforeEvent.jumpC] * integrand[index.beforeEvent.jumpC]
integral <- rowCumSum(integrand)
integral2 <- rowCumSum(integrand2)
augTerm[,beforeTau.nJumpC!=0] <- F1.tau[,beforeTau.nJumpC!=0,drop=FALSE] * integral[,beforeTau.nJumpC.n0,drop=FALSE] - integral2[,beforeTau.nJumpC.n0,drop=FALSE]
## survival at all jump of the censoring process
S.jump <- predictRisk(object.event, type = "survival", newdata = mydata, times = time.jumpC-tol, product.limit = product.limit,
iid = (method.iid==2)*return.iid.nuisance)
if((method.iid==2)*return.iid.nuisance){
out$store$iid.nuisance.survival <- attr(S.jump,"iid")
attr(S.jump,"iid") <- NULL
}

## martingale for the censoring process
## at all times of jump of the censoring process
G.jump <- 1-predictRisk(object.censor, newdata = mydata, times = if(index.lastjumpC>1){c(0,time.jumpC[1:(index.lastjumpC-1)])}else{0},
product.limit = product.limit, iid = (method.iid==2)*return.iid.nuisance)

if(return.iid.nuisance && (method.iid==2)){
out$store$iid.nuisance.censoring <- -attr(G.jump,"iid")
attr(G.jump,"iid") <- NULL
}
dN.jump <- do.call(cbind,lapply(time.jumpC, function(iJump){(mydata[[eventVar.time]] == iJump)*(mydata[[eventVar.status]] == level.censoring)}))
dLambda.jump <- predictCox(object.censor, newdata = mydata, times = time.jumpC, type = "hazard", iid = (method.iid==2)*return.iid.nuisance)
if((method.iid==2)*return.iid.nuisance){
out$store$iid.nuisance.martingale <- dLambda.jump$hazard.iid
}
dM.jump <- dN.jump - dLambda.jump$hazard
## integral = \int_0^min(T_i,\tau) (F1(\tau|A_i,W_i)-F1(t|A_i,W_i)) / S(t|A_i,W_i) * dM_i^C(t)/Gc(t|A_i,W_i)
## = F1(\tau|A_i,W_i) \int_0^min(T_i,\tauf) dM_i^C(t) / (S(t|A_i,W_i) * Gc(t|A_i,W_i)) - \int_0^min(T_i,\tauf) F1(t|A_i,W_i) dM_i^C(t) / (S(t|A_i,W_i) * Gc(t|A_i,W_i))
integrand <- matrix(0, nrow = n.obs, ncol = index.lastjumpC)
integrand2 <- matrix(0, nrow = n.obs, ncol = index.lastjumpC)
index.beforeEvent.jumpC <- which(beforeEvent.jumpC) ## identify which increment to put in the integral, i.e. when the individual was still at risk of being censored
integrand[index.beforeEvent.jumpC] <- dM.jump[index.beforeEvent.jumpC] / (G.jump[index.beforeEvent.jumpC] * S.jump[index.beforeEvent.jumpC])
integrand2[index.beforeEvent.jumpC] <- F1.jump[index.beforeEvent.jumpC] * integrand[index.beforeEvent.jumpC]
integral <- rowCumSum(integrand)
integral2 <- rowCumSum(integrand2)
augTerm[,beforeTau.nJumpC!=0] <- F1.tau[,beforeTau.nJumpC!=0,drop=FALSE] * integral[,beforeTau.nJumpC.n0,drop=FALSE] - integral2[,beforeTau.nJumpC.n0,drop=FALSE]
}

## ** Compute individual contribution to the ATE + influence function for the G-formula
for(iC in 1:n.contrasts){ ## iC <- 1
if(attr(estimator,"export.GFORMULA")){
Expand Down Expand Up @@ -297,12 +298,10 @@ ATE_TI <- function(object.event,
out$store$iid.IPTW[[iC]][data.index,] <- out$store$iid.IPTW[[iC]][data.index,] + rowCenter_cpp(iIID.ate, center = iATE)/n.obs
}
if(attr(estimator,"export.AIPTW")){
if(attr(estimator,"IPCW")){
if(inherits(object.event,"wglm")){
iIID.ate <- F1.ctf.tau[[iC]] + colMultiply_cpp(iW.IPCW * Y.tau - F1.ctf.tau[[iC]], scale = iW.IPTW[,iC])
}else{
iIID.ate <- F1.ctf.tau[[iC]] + colMultiply_cpp(iW.IPCW * Y.tau - F1.ctf.tau[[iC]] + augTerm, scale = iW.IPTW[,iC])
}
if(attr(estimator,"integral")){ ## full augmentation
iIID.ate <- F1.ctf.tau[[iC]] + colMultiply_cpp(iW.IPCW * Y.tau - F1.ctf.tau[[iC]] + augTerm, scale = iW.IPTW[,iC])
}else if(inherits(object.event,"wglm") && attr(estimator,"IPCW")){ ## not full augmentation
iIID.ate <- F1.ctf.tau[[iC]] + colMultiply_cpp(iW.IPCW * Y.tau - F1.ctf.tau[[iC]], scale = iW.IPTW[,iC])
}else{
iIID.ate <- F1.ctf.tau[[iC]] + colMultiply_cpp(Y.tau - F1.ctf.tau[[iC]], scale = iW.IPTW[,iC])
}
Expand All @@ -315,8 +314,7 @@ ATE_TI <- function(object.event,
out$store$iid.AIPTW[[iC]][data.index,] <- out$store$iid.AIPTW[[iC]][data.index,] + rowCenter_cpp(iIID.ate, center = iATE)/n.obs
}
}



## ** save quantities useful for the calculation of iid.nuisance
if(return.iid.nuisance){
out$store$n.obs <- n.obs
Expand Down
13 changes: 6 additions & 7 deletions R/ate.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,9 @@
## author: Thomas Alexander Gerds
## created: Oct 23 2016 (08:53)
## Version:
## last-updated: sep 11 2024 (18:32)
## last-updated: Oct 14 2024 (09:36)
## By: Brice Ozenne
## Update #: 2360
## Update #: 2368
#----------------------------------------------------------------------
##
### Commentary:
Expand Down Expand Up @@ -836,7 +836,7 @@ ate_initArgs <- function(object.event,
if(inherits(model.censor,"coxph") || inherits(model.censor,"cph") || inherits(model.censor,"phreg")){
censoringMF <- coxModelFrame(model.censor)
test.censor <- censoringMF$status == 1
n.censor <- sapply(times, function(t){sum(test.censor * (censoringMF$stop <= t))})
n.censor <- sapply(times, function(t){sum(test.censor * (censoringMF$stop < t))})

info.censor <- SurvResponseVar(coxFormula(model.censor))
censorVar.status <- info.censor$status
Expand All @@ -845,16 +845,15 @@ ate_initArgs <- function(object.event,
}else{ ## G-formula or IPTW (no censoring)
censorVar.status <- NA
censorVar.time <- NA



if(inherits(model.event,"CauseSpecificCox")){
test.censor <- model.event$response[,"status"] == 0
n.censor <- sapply(times, function(t){sum(test.censor * (model.event$response[,"time"] <= t))})
n.censor <- sapply(times, function(t){sum(test.censor * (model.event$response[,"time"] < t))})
level.censoring <- attr(model.event$response,"cens.code")
}else if(inherits(model.event,"coxph") || inherits(model.event,"cph") || inherits(model.event,"phreg")){
censoringMF <- coxModelFrame(model.event)
test.censor <- censoringMF$status == 0
n.censor <- sapply(times, function(t){sum(test.censor * (censoringMF$stop <= t))})
n.censor <- sapply(times, function(t){sum(test.censor * (censoringMF$stop < t))})
level.censoring <- 0
}else if(inherits(model.event,"wglm")){
n.censor <- object.event$n.censor
Expand Down
26 changes: 24 additions & 2 deletions R/calcSeCSC.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,9 @@
## author: Brice Ozenne
## created: maj 27 2017 (21:23)
## Version:
## last-updated: mar 10 2023 (14:27)
## last-updated: Oct 14 2024 (12:59)
## By: Brice Ozenne
## Update #: 1099
## Update #: 1110
#----------------------------------------------------------------------
##
### Commentary:
Expand Down Expand Up @@ -183,6 +183,28 @@ calcSeCSC <- function(object, cif, hazard, cumhazard, survival, object.time, obj
if("iid" %in% export){
out$iid <- aperm(out$iid, c(2,3,1))
}
if("average.iid" %in% export && !is.null(attr(export,"factor"))){
## calcSeCif2_cpp average ignoring the weights defined in attr(export,"factor")
factor <- attr(export, "factor")
if(diag){
## Not necessary due to error message
## Attribute "factor" of argument 'average.iid' not available when 'iid' is TRUE
}else{
## May occur when compressing data as it bypass the error message
n.newobs <- dim(out$iid)[3]
out$average.iid <- lapply(1:length(factor), function(iF){ ## iObs <- 1
Reduce("+",lapply(1:n.newobs, function(iObs){
if(nTimes==1){
return(cbind(out$iid[,,iObs]*factor[[1]][iObs,]))
}else{
return(rowMultiply_cpp(out$iid[,,iObs], factor[[1]][iObs,]))
}
}))/n.newobs})
if(length(factor)==1){
out$average.iid <- out$average.iid[[1]]
}
}
}

}else if("average.iid" %in% export){

Expand Down
Loading

0 comments on commit fe6aa16

Please sign in to comment.