Skip to content

Commit

Permalink
clean up protected names
Browse files Browse the repository at this point in the history
  • Loading branch information
@tagteam
tagteam committed Jun 4, 2024
1 parent 34184b3 commit 56bd408
Show file tree
Hide file tree
Showing 18 changed files with 263 additions and 191 deletions.
29 changes: 13 additions & 16 deletions R/AUC.binary.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,9 @@
## Author: Thomas Alexander Gerds
## Created: Jan 11 2022 (17:04)
## Version:
## Last-Updated: Jun 4 2024 (07:21)
## Last-Updated: Jun 4 2024 (15:24)
## By: Thomas Alexander Gerds
## Update #: 6
## Update #: 16
#----------------------------------------------------------------------
##
### Commentary:
Expand All @@ -14,7 +14,6 @@
#----------------------------------------------------------------------
##
### Code:

AUC.binary <- function(DT,
breaks=NULL,
se.fit,
Expand All @@ -31,7 +30,7 @@ AUC.binary <- function(DT,
ROC,
cutpoints,
...){
PPV=NPV=Prisks=Prisks2=model=risk=ReSpOnSe=FPR=TPR=riskRegression_ID=NULL
PPV=NPV=Prisks=Prisks2=model=risk=riskRegression_event=FPR=TPR=riskRegression_ID=NULL
## do not want to depend on Daim as they turn marker to ensure auc > 0.5
delongtest <- function(risk,
score,
Expand All @@ -54,7 +53,6 @@ AUC.binary <- function(DT,
#riskcontrols <- as.matrix(risk[Controls,])
riskcontrols <- as.matrix(risk[!Cases,])
riskcases <- as.matrix(risk[Cases,])

# new method, uses a fast implementation of delongs covariance matrix
# Fast Implementation of DeLong’s Algorithm for Comparing the Areas Under Correlated Receiver Operating Characteristic Curves
# article can be found here:
Expand Down Expand Up @@ -124,7 +122,6 @@ AUC.binary <- function(DT,
}
out
}

auRoc.numeric <- function(X,D,breaks,ROC,cutpoints=NULL){
if (is.null(breaks)) breaks <- rev(sort(unique(X))) ## need to reverse when high X is concordant with {response=1}
TPR <- c(prodlim::sindex(jump.times=X[D==1],eval.times=breaks,comp="greater",strict=FALSE)/sum(D==1))
Expand Down Expand Up @@ -162,10 +159,10 @@ AUC.binary <- function(DT,
ROC <- TRUE
}
if (is.factor(DT[["risk"]])){
score <- aucDT[,auRoc.factor(risk,ReSpOnSe,ROC=ROC),by=list(model)]
score <- aucDT[,auRoc.factor(risk,riskRegression_event,ROC=ROC),by=list(model)]
}
else{
score <- aucDT[,auRoc.numeric(risk,ReSpOnSe,breaks=breaks,ROC=ROC,cutpoints=cutpoints),by=list(model)]
score <- aucDT[,auRoc.numeric(risk,riskRegression_event,breaks=breaks,ROC=ROC,cutpoints=cutpoints),by=list(model)]
}
if (ROC==FALSE){
setnames(score,"V1","AUC")
Expand All @@ -184,20 +181,20 @@ AUC.binary <- function(DT,
for (i in 1:length(cutpoints)){
temp.TPR <- subset(temp.TPR.ic,cutpoints==cutpoints[i])
aucDT.temp <- merge(aucDT,temp.TPR)
some.fun <- function(ReSpOnSe,risk,TPR,FPR,PPV,NPV,Prisks,Prisks2,cut,N){
meanY <- mean(ReSpOnSe)
some.fun <- function(riskRegression_event,risk,TPR,FPR,PPV,NPV,Prisks,Prisks2,cut,N){
meanY <- mean(riskRegression_event)
out <- list(TPR = TPR[1],
SE.TPR = sd(ReSpOnSe/meanY * ((risk > cut)-TPR))/sqrt(N),
SE.TPR = sd(riskRegression_event/meanY * ((risk > cut)-TPR))/sqrt(N),
FPR = FPR[1],
SE.FPR = sd((1-ReSpOnSe)/(1-meanY) * ((risk > cut)-FPR))/sqrt(N),
SE.FPR = sd((1-riskRegression_event)/(1-meanY) * ((risk > cut)-FPR))/sqrt(N),
PPV = PPV[1],
SE.PPV = sd((risk > cut)/Prisks[1] * (ReSpOnSe - PPV))/sqrt(N),
SE.PPV = sd((risk > cut)/Prisks[1] * (riskRegression_event - PPV))/sqrt(N),
NPV = NPV[1],
SE.NPV = sd((risk <= cut)/Prisks2[1] * ((1-ReSpOnSe) - NPV))/sqrt(N),
SE.NPV = sd((risk <= cut)/Prisks2[1] * ((1-riskRegression_event) - NPV))/sqrt(N),
cutpoint = cut)
out
}
res.cut[[i]] <- aucDT.temp[,some.fun(ReSpOnSe,risk,TPR,FPR,PPV,NPV,Prisks,Prisks2,cutpoints[i],N), by = list(model)]
res.cut[[i]] <- aucDT.temp[,some.fun(riskRegression_event,risk,TPR,FPR,PPV,NPV,Prisks,Prisks2,cutpoints[i],N), by = list(model)]
}
output <- list(score=AUC,ROC=ROC, res.cut=do.call("rbind",res.cut))
}
Expand All @@ -210,7 +207,7 @@ AUC.binary <- function(DT,
delong.res <- delongtest(risk=xRisk,
score=output$score,
dolist=dolist,
response=aucDT[model==model[1],ReSpOnSe],
response=aucDT[model==model[1],riskRegression_event],
cause="1",
alpha=alpha,
multi.split.test=multi.split.test,
Expand Down
99 changes: 80 additions & 19 deletions R/AUC.competing.risks.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,9 @@
## Author: Thomas Alexander Gerds
## Created: Jan 11 2022 (17:06)
## Version:
## Last-Updated: Jun 4 2024 (11:56)
## Last-Updated: Jun 4 2024 (14:35)
## By: Thomas Alexander Gerds
## Update #: 29
## Update #: 36
#----------------------------------------------------------------------
##
### Commentary:
Expand Down Expand Up @@ -62,11 +62,27 @@ AUC.competing.risks <- function(DT,
if (!is.null(cutpoints)){
breaks <- sort(cutpoints,decreasing = TRUE)
aucDT[,nth.times:=as.numeric(factor(times))]
cutpoint.helper.fun <- function(FPR,TPR,risk,ipcwCases,ipcwControls1,ipcwControls2, N, riskRegression_time,times,riskRegression_event,cens.model,nth.times,conservative, IC.G, cutpoints,se.fit){
cutpoint.helper.fun <- function(FPR,
TPR,
risk,
ipcwCases,
ipcwControls1,
ipcwControls2,
N,
riskRegression_time,
times,
riskRegression_event,
cens.model,
nth.times,
conservative,
IC.G,
cutpoints,
se.fit){
den_TPR<-sum(ipcwCases) ## estimate the cumulative incidence via IPCW
den_FPR<-sum(ipcwControls1+ipcwControls2)
indeces <- sindex(risk,cutpoints,comp = "greater",TRUE)
res <- list()
# FIXME
ordered <- order(riskRegression_time) ## can probably move this outside to improve computation time, for now keep it
for (i in 1:length(cutpoints)){
den_PPV <- sum(ipcwCases[risk > cutpoints[i]]+ipcwControls1[risk > cutpoints[i]] + ipcwControls2[risk > cutpoints[i]])
Expand All @@ -77,8 +93,24 @@ AUC.competing.risks <- function(DT,
if (se.fit){
IC0.TPR <- ipcwCases*N*((risk > cutpoints[i])-TPRi)/den_TPR
IC0.FPR <- (ipcwControls1+ipcwControls2)*N*((risk > cutpoints[i])-FPRi)/(1-den_TPR)
SE.TPR <- sd(getInfluenceCurve.Brier(times,riskRegression_time[ordered],IC0.TPR[ordered],IC0.TPR[ordered],IC.G,cens.model,nth.times,conservative,riskRegression_event[ordered]))/sqrt(N)
SE.FPR <- sd(getInfluenceCurve.Brier(times,riskRegression_time[ordered],IC0.FPR[ordered],IC0.FPR[ordered],IC.G,cens.model,nth.times,conservative,riskRegression_event[ordered]))/sqrt(N)
SE.TPR <- sd(getInfluenceCurve.Brier(t = times,
time = riskRegression_time[ordered],
IC0 = IC0.TPR[ordered],
residuals = IC0.TPR[ordered],
IC.G = IC.G,
cens.model = cens.model,
nth.times = nth.times,
conservative = conservative,
event = riskRegression_event[ordered]))/sqrt(N)
SE.FPR <- sd(getInfluenceCurve.Brier(t = times,
time = riskRegression_time[ordered],
IC0 = IC0.FPR[ordered],
residuals = IC0.FPR[ordered],
IC.G = IC.G,
cens.model = cens.model,
nth.times = nth.times,
conservative = conservative,
event = riskRegression_event[ordered]))/sqrt(N)
}
else {
SE.TPR <- SE.FPR <- NA
Expand All @@ -92,16 +124,19 @@ AUC.competing.risks <- function(DT,
PPV <- (TPRi*den_TPR)/den_PPV
if (se.fit){
IC0.PPV <- (risk > cutpoints[i])/den_PPV*(((ipcwCases+ipcwControls2)*N)*(1*(riskRegression_event==1)-1*(riskRegression_event!=0)*PPV)-ipcwControls1*N*PPV) #OBS, check other causes, paul's implementation
SE.PPV <- sd(getInfluenceCurve.Brier(times,riskRegression_time[ordered],IC0.PPV[ordered],IC0.PPV[ordered],IC.G,cens.model,nth.times,conservative,riskRegression_event[ordered]))/sqrt(N)
SE.PPV <- sd(getInfluenceCurve.Brier(t = times,
time = riskRegression_time[ordered],
IC0 = IC0.PPV[ordered],
residuals = IC0.PPV[ordered],
IC.G = IC.G,
cens.model = cens.model,
nth.times = nth.times,
conservative = conservative,
event = riskRegression_event[ordered]))/sqrt(N)
}
else {
SE.PPV <- NA
}
## Alternative implementation
# IC0.PPV <- (risk > cutpoints[i])/den_PPV*(ipcwCases*N - PPV)
# weights.PPV <- ((risk > cutpoints[i])/(den_PPV))*((1-PPV)*ipcwCases*N - PPV*(N*ipcwControls1+N*ipcwControls2)) #include censoring from denominator
# weights.PPV <- ((risk > cutpoints[i])/(den_PPV))*(ipcwCases*N) #exclude censoring from denominator
# SE.PPV <- sd(getInfluenceCurve.Brier(times,riskRegression_time[ordered],IC0.PPV[ordered],weights.PPV[ordered],IC.G,cens.model,nth.times,conservative,riskRegression_event[ordered]))/sqrt(N)
}
else {
PPV <- NA
Expand All @@ -110,25 +145,51 @@ AUC.competing.risks <- function(DT,
NPV <- ((1-FPRi)*den_FPR)/den_NPV
if (se.fit){
IC0.NPV <- (risk <= cutpoints[i])/den_NPV*(((ipcwCases+ipcwControls2)*N)*(1*(riskRegression_event!=1 & riskRegression_event!=0)-1*(riskRegression_event!=0)*NPV)+ipcwControls1*N*(1-NPV)) #OBS, check other causes, paul's implementation
SE.NPV <- sd(getInfluenceCurve.Brier(times,riskRegression_time[ordered],IC0.NPV[ordered],IC0.NPV[ordered],IC.G,cens.model,nth.times,conservative,riskRegression_event[ordered]))/sqrt(N)
SE.NPV <- sd(getInfluenceCurve.Brier(t = times,
time = riskRegression_time[ordered],
IC0 = IC0.NPV[ordered],
residuals = IC0.NPV[ordered],
IC.G = IC.G,
cens.model = cens.model,
nth.times = nth.times,
conservative = conservative,
event = riskRegression_event[ordered]))/sqrt(N)
}
else {
SE.NPV <- NA
}

## Alternative implementation
# IC0.NPV <- (risk <= cutpoints[i])/den_NPV*((ipcwControls1+ipcwControls2)*N - NPV)
# weights.NPV <- ((risk <= cutpoints[i])/(den_NPV))*((ipcwControls1+ipcwControls2)*N)
# SE.NPV <- sd(getInfluenceCurve.Brier(times,riskRegression_time[ordered],IC0.NPV[ordered],weights.NPV[ordered],IC.G,cens.model,nth.times,conservative,riskRegression_event[ordered]))/sqrt(N)
}
else {
NPV <- NA
}
res[[i]] <- data.table(risk = cutpoints[i], TPR=TPRi, SE.TPR=SE.TPR,FPR=FPRi, SE.FPR=SE.FPR,PPV=PPV,SE.PPV=SE.PPV,NPV=NPV,SE.NPV=SE.NPV)
res[[i]] <- data.table(risk = cutpoints[i],
TPR=TPRi,
SE.TPR=SE.TPR,
FPR=FPRi,
SE.FPR=SE.FPR,
PPV=PPV,
SE.PPV=SE.PPV,
NPV=NPV,
SE.NPV=SE.NPV)
}
do.call("rbind",res)
}
output <- list(res.cut=aucDT[, cutpoint.helper.fun(FPR,TPR,risk,ipcwCases,ipcwControls1,ipcwControls2, N, riskRegression_time,times[1],riskRegression_status*riskRegression_event,cens.model,nth.times[1],conservative, MC, cutpoints,se.fit),by=list(model,times)])
output <- list(res.cut=aucDT[, cutpoint.helper.fun(FPR = FPR,
TPR = TPR,
risk = risk,
ipcwCases = ipcwCases,
ipcwControls1 = ipcwControls1,
ipcwControls2 = ipcwControls2,
N = N,
riskRegression_time = riskRegression_time,
times = times[1],
riskRegression_event = riskRegression_status*riskRegression_event,
cens.model = cens.model,
nth.times = nth.times[1],
conservative = conservative,
IC.G = MC,
cutpoints = cutpoints,
se.fit = se.fit),by=list(model,times)])
}
else if (ROC) {
if (is.null(breaks)){
Expand Down
Loading

0 comments on commit 56bd408

Please sign in to comment.