diff --git a/R/AUC.binary.R b/R/AUC.binary.R index b88ec1b..b62a9d9 100644 --- a/R/AUC.binary.R +++ b/R/AUC.binary.R @@ -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: @@ -14,7 +14,6 @@ #---------------------------------------------------------------------- ## ### Code: - AUC.binary <- function(DT, breaks=NULL, se.fit, @@ -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, @@ -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: @@ -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)) @@ -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") @@ -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)) } @@ -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, diff --git a/R/AUC.competing.risks.R b/R/AUC.competing.risks.R index 9f2b8f9..0673ebf 100644 --- a/R/AUC.competing.risks.R +++ b/R/AUC.competing.risks.R @@ -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: @@ -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]]) @@ -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 @@ -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 @@ -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)){ diff --git a/R/AUC.survival.R b/R/AUC.survival.R index 69ac1f8..2c205b3 100644 --- a/R/AUC.survival.R +++ b/R/AUC.survival.R @@ -3,9 +3,9 @@ ## Author: Thomas Alexander Gerds ## Created: Jan 11 2022 (17:06) ## Version: -## Last-Updated: Jun 4 2024 (11:51) +## Last-Updated: Jun 4 2024 (14:34) ## By: Thomas Alexander Gerds -## Update #: 16 +## Update #: 29 #---------------------------------------------------------------------- ## ### Commentary: @@ -44,7 +44,7 @@ AUC.survival <- function(DT, ## order data data.table::setorder(aucDT,model,times,-risk) ## identify cases and controls - aucDT[,Cases:=(riskRegression_time <= times & status==cause)] + aucDT[,Cases:=(riskRegression_time <= times & riskRegression_status==cause)] aucDT[,Controls:=(riskRegression_time > times)] ## prepare Weights aucDT[Cases==0,ipcwCases:=0] @@ -61,11 +61,26 @@ AUC.survival <- 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,ipcwControls, N, riskRegression_time,times,event,cens.model,nth.times,conservative, IC.G, cutpoints,se.fit){ + cutpoint.helper.fun <- function(FPR, + TPR, + risk, + ipcwCases, + ipcwControls, + 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(ipcwControls) indeces <- sindex(risk,cutpoints,comp = "greater",TRUE) res <- list() + ## browser() ordered <- order(riskRegression_time) ## can probably move this outside to improve computation time, for now keep it SE.TPR <- SE.FPR <- SE.PPV <- SE.NPV <- NA for (i in 1:length(cutpoints)){ @@ -76,9 +91,25 @@ AUC.survival <- function(DT, FPRi <- FPR[indeces[i]] if (se.fit){ IC0.TPR <- ipcwCases*N*((risk > cutpoints[i])-TPRi)/den_TPR - SE.TPR <- sd(getInfluenceCurve.Brier(times,riskRegression_time[ordered],IC0.TPR[ordered],IC0.TPR[ordered],IC.G,cens.model,nth.times,conservative,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) IC0.FPR <- (ipcwControls)*N*((risk > cutpoints[i])-FPRi)/(1-den_TPR) - SE.FPR <- sd(getInfluenceCurve.Brier(times,riskRegression_time[ordered],IC0.FPR[ordered],IC0.FPR[ordered],IC.G,cens.model,nth.times,conservative,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 { @@ -87,15 +118,18 @@ AUC.survival <- function(DT, if (den_PPV > 1e-10){ PPV <- (TPRi*den_TPR)/den_PPV if (se.fit){ - IC0.PPV <- (risk > cutpoints[i])/den_PPV*(((ipcwCases)*N)*(1*(event==1)-1*(event!=0)*PPV)-ipcwControls*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,event[ordered]))/sqrt(N) + #OBS, check other causes, paul's implementation + IC0.PPV <- (risk > cutpoints[i])/den_PPV*(((ipcwCases)*N)*(1*(riskRegression_event==1)-1*(riskRegression_event!=0)*PPV)-ipcwControls*N*PPV) + 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) } - - ## 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,event[ordered]))/sqrt(N) } else { PPV <- NA @@ -103,14 +137,17 @@ AUC.survival <- function(DT, if (den_NPV > 1e-10){ NPV <- ((1-FPRi)*den_FPR)/den_NPV if (se.fit){ - IC0.NPV <- (risk <= cutpoints[i])/den_NPV*(((ipcwCases)*N)*(1*(event!=1 & event!=0)-1*(event!=0)*NPV)+ipcwControls*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,event[ordered]))/sqrt(N) + IC0.NPV <- (risk <= cutpoints[i])/den_NPV*(((ipcwCases)*N)*(1*(riskRegression_event!=1 & riskRegression_event!=0)-1*(riskRegression_event!=0)*NPV)+ipcwControls*N*(1-NPV)) #OBS, check other causes, paul's implementation + 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) } - - ## 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,event[ordered]))/sqrt(N) } else { NPV <- NA @@ -119,7 +156,22 @@ AUC.survival <- function(DT, } do.call("rbind",res) } - output <- list(res.cut=aucDT[, cutpoint.helper.fun(FPR,TPR,risk,ipcwCases,ipcwControls, N, riskRegression_time,times[1],status,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, + ipcwControls = ipcwControls, + N = N, + riskRegression_time = riskRegression_time, + times = times[1], + riskRegression_event = riskRegression_status, + 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==TRUE) { if (is.null(breaks)){ @@ -146,7 +198,7 @@ AUC.survival <- function(DT, data.table::setorder(aucDT,model,times,riskRegression_ID) aucDT[,IF.AUC:=getInfluenceCurve.AUC(t = times[1], time = riskRegression_time, - status = riskRegression_status, + event = riskRegression_status, WTi = WTi, Wt = Wt, risk = risk, diff --git a/R/Brier.binary.R b/R/Brier.binary.R index 16a1c54..e36312a 100644 --- a/R/Brier.binary.R +++ b/R/Brier.binary.R @@ -3,9 +3,9 @@ ## Author: Thomas Alexander Gerds ## Created: Jan 11 2022 (17:03) ## Version: -## Last-Updated: Jun 4 2024 (07:21) +## Last-Updated: Jun 4 2024 (15:08) ## By: Thomas Alexander Gerds -## Update #: 6 +## Update #: 7 #---------------------------------------------------------------------- ## ### Commentary: @@ -29,8 +29,8 @@ Brier.binary <- function(DT, dolist, keep.residuals=FALSE, ...){ - residuals=Brier=risk=model=ReSpOnSe=lower=upper=se=riskRegression_ID=IF.Brier=NULL - DT[,residuals:=(ReSpOnSe-risk)^2] + residuals=Brier=risk=model=riskRegression_event=lower=upper=se=riskRegression_ID=IF.Brier=NULL + DT[,residuals:=(riskRegression_event-risk)^2] if (se.fit==1L){ data.table::setkey(DT,model,riskRegression_ID) score <- DT[,data.table::data.table(Brier=sum(residuals)/N,se=sd(residuals)/sqrt(N)),by=list(model)] @@ -86,7 +86,7 @@ Brier.binary <- function(DT, list(iid.decomp = DT[,data.table::data.table(riskRegression_ID,model,IF.Brier)])) } if (keep.residuals[[1]] == TRUE) { - output <- c(output,list(residuals=DT[,data.table::data.table(riskRegression_ID,ReSpOnSe,model,risk,residuals)])) + output <- c(output,list(residuals=DT[,data.table::data.table(riskRegression_ID,riskRegression_event,model,risk,residuals)])) } output } diff --git a/R/Brier.competing.risks.R b/R/Brier.competing.risks.R index b110469..1f029a8 100644 --- a/R/Brier.competing.risks.R +++ b/R/Brier.competing.risks.R @@ -3,9 +3,9 @@ ## Author: Thomas Alexander Gerds ## Created: Jan 11 2022 (17:04) ## Version: -## Last-Updated: Jun 4 2024 (11:55) +## Last-Updated: Jun 4 2024 (19:08) ## By: Thomas Alexander Gerds -## Update #: 7 +## Update #: 8 #---------------------------------------------------------------------- ## ### Commentary: @@ -56,7 +56,7 @@ Brier.competing.risks <- function(DT, cens.model=cens.model, conservative = conservative, nth.times=nth.times[1], - riskRegression_event = riskRegression_status*riskRegression_event),by=list(model,times)] + event = riskRegression_status*riskRegression_event),by=list(model,times)] score <- DT[,data.table(Brier=sum(residuals)/N, se=sd(IF.Brier)/sqrt(N)),by=list(model,times)] if (se.fit==TRUE){ diff --git a/R/Brier.survival.R b/R/Brier.survival.R index cfb88ee..16fb398 100644 --- a/R/Brier.survival.R +++ b/R/Brier.survival.R @@ -3,9 +3,9 @@ ## Author: Thomas Alexander Gerds ## Created: Jan 11 2022 (17:04) ## Version: -## Last-Updated: Jun 4 2024 (11:50) +## Last-Updated: Jun 4 2024 (14:11) ## By: Thomas Alexander Gerds -## Update #: 7 +## Update #: 8 #---------------------------------------------------------------------- ## ### Commentary: @@ -31,7 +31,7 @@ Brier.survival <- function(DT, keep.residuals=FALSE, IC.data, ...){ - IC0=IPCW=nth.times=riskRegression_ID=times=raw.Residuals=risk=Brier=residuals=WTi=Wt=status=setorder=model=IF.Brier=data.table=sd=lower=qnorm=se=upper=NULL + IC0=IPCW=nth.times=riskRegression_ID = riskRegression_time = riskRegression_status=times=raw.Residuals=risk=Brier=residuals=WTi=Wt=status=setorder=model=IF.Brier=data.table=sd=lower=qnorm=se=upper=NULL ## compute 0/1 outcome: DT[riskRegression_time<=times & status==1,residuals:=(1-risk)^2/WTi] DT[riskRegression_time<=times & status==0,residuals:=0] diff --git a/R/Score.R b/R/Score.R index f272796..6095e4a 100644 --- a/R/Score.R +++ b/R/Score.R @@ -212,7 +212,6 @@ ##' \code{N=300} ##' \code{mean(sapply(1:B, function(b){match(1,sample(1:N,size=N,replace=TRUE),nomatch=0)}))} ##' -##' ##' ## Bootstrap without replacement (training size set to be 70 percent of data) ##' B=10, M=.7 ##' @@ -523,7 +522,7 @@ Score.list <- function(object, roc.grid=switch(roc.method,"vertical"=seq(0,1,.01),"horizontal"=seq(1,0,-.01)), cutpoints = NULL, ...){ - se.conservative=IPCW=IF.AUC.conservative=IF.AUC0=IF.AUC=IC0=Brier=AUC=casecontrol=se=nth.times=time=status=riskRegression_ID=WTi=risk=IF.Brier=lower=upper=crossval=b=time=status=model=reference=p=model=pseudovalue=ReSpOnSe=residuals=event=j=NULL + se.conservative=IPCW=IF.AUC.conservative=IF.AUC0=IF.AUC=IC0=Brier=AUC=casecontrol=se=nth.times=riskRegression_time=riskRegression_time_status=riskRegression_ID=WTi=risk=IF.Brier=lower=upper=crossval=b=time=status=model=reference=p=model=pseudovalue=riskRegression_event=residuals=event=j=NULL # }}} theCall <- match.call() @@ -648,20 +647,16 @@ c.f., Chapter 7, Section 5 in Gerds & Kattan 2021. Medical risk prediction model if (length(grep("^riskRegression_",names(data)))>0){ stop("Internal variable name clash. Cannot have variables with names starting with riskRegression_ in the data.") } - ## put ReSpOnSe for binary and (time, event, status) in the first column(s) - if (response.type=="binary"){ - data[,event:=factor(ReSpOnSe, - levels=1:(length(states)+1), - labels=c(states,"event-free"))] - } + ## put change names of outcome if (response.type %in% c("survival","competing.risks")){ colnames(response) <- paste0("riskRegression_",colnames(response)) - # FIXME: is order okay? + # FIXME: check if order still correct data <- cbind(response,data) N <- as.numeric(NROW(data)) neworder <- data[,order(riskRegression_time,-riskRegression_status)] data.table::setorder(data,riskRegression_time,-riskRegression_status) }else{ + # add variable riskRegression_event data <- cbind(response,data) N <- as.numeric(NROW(data)) neworder <- 1:N diff --git a/R/computePerformance.R b/R/computePerformance.R index 11a3520..c2bf437 100644 --- a/R/computePerformance.R +++ b/R/computePerformance.R @@ -3,9 +3,9 @@ ## Author: Thomas Alexander Gerds ## Created: Feb 27 2022 (09:12) ## Version: -## Last-Updated: Jun 4 2024 (11:41) +## Last-Updated: Jun 4 2024 (12:06) ## By: Thomas Alexander Gerds -## Update #: 21 +## Update #: 22 #---------------------------------------------------------------------- ## ### Commentary: @@ -62,7 +62,7 @@ computePerformance <- function(DT, keep.residuals=keep.residuals, keep.vcov=keep.vcov, keep.iid=keep.iid, - dolist=dolist,Q=probs,ROC=FALSE,MC=MC,IC.data=IC.data,breaks=breaks,cutpoints=cutpoints) ## will break survival + dolist=dolist,Q=probs,ROC=FALSE,MC=MC,IC.data=IC.data,breaks=breaks,cutpoints=cutpoints) if (response.type=="competing.risks") { input <- c(input,list(cause=cause,states=states)) } diff --git a/R/crossvalPerf.loob.AUC.R b/R/crossvalPerf.loob.AUC.R index 3f2f60b..1ce2b08 100644 --- a/R/crossvalPerf.loob.AUC.R +++ b/R/crossvalPerf.loob.AUC.R @@ -3,9 +3,9 @@ ## Author: Thomas Alexander Gerds ## Created: Jun 4 2024 (09:20) ## Version: -## Last-Updated: Jun 4 2024 (11:24) +## Last-Updated: Jun 4 2024 (15:08) ## By: Thomas Alexander Gerds -## Update #: 8 +## Update #: 11 #---------------------------------------------------------------------- ## ### Commentary: @@ -37,7 +37,7 @@ crossvalPerf.loob.AUC <- function(times, # crossvalPerf.loob.Brier has more arguments ...) { # initializing output - bfold <- fold <- AUC <- ReSpOnSe <- model <- b <- risk <- casecontrol <- IF.AUC <- IF.AUC0 <- se <- IF.AUC.conservative <- se.conservative <- lower <- upper <- NF <- reference <- event <- status0 <- NULL + bfold <- fold <- AUC <- riskRegression_event <- model <- b <- risk <- casecontrol <- IF.AUC <- IF.AUC0 <- se <- IF.AUC.conservative <- se.conservative <- lower <- upper <- NF <- reference <- riskRegression_event <- riskRegression_time <- riskRegression_status0 <- riskRegression_ID <- NULL if (response.type=="binary") { auc.loob <- data.table(expand.grid(times=0,model=mlevs)) #add times to auc.loob; now we can write less code for the same thing! times <- 0 @@ -72,7 +72,7 @@ crossvalPerf.loob.AUC <- function(times, if (response.type=="binary"){ t <- 0 ## the following indices have to be logical!!! - cases.index <- Response[,ReSpOnSe==1] + cases.index <- Response[,riskRegression_event==1] controls.index <- !cases.index }else{ t <- times[s] @@ -161,7 +161,6 @@ crossvalPerf.loob.AUC <- function(times, muCase = muCase, muControls = muControls, nu = nu1tauPm, - # FIXME HERE firsthit = first_hits[s], cases = cases.index, controls =controls.index, diff --git a/R/crossvalPerf.loob.Brier.R b/R/crossvalPerf.loob.Brier.R index 5617400..a0c5f0b 100644 --- a/R/crossvalPerf.loob.Brier.R +++ b/R/crossvalPerf.loob.Brier.R @@ -3,9 +3,9 @@ ## Author: Thomas Alexander Gerds ## Created: Jun 4 2024 (09:16) ## Version: -## Last-Updated: Jun 4 2024 (09:21) +## Last-Updated: Jun 4 2024 (15:08) ## By: Thomas Alexander Gerds -## Update #: 5 +## Update #: 9 #---------------------------------------------------------------------- ## ### Commentary: @@ -37,25 +37,25 @@ crossvalPerf.loob.Brier <- function(times, conservative, cens.model, cause){ - status <- residuals <- risk <- WTi <- event <- ReSpOnSe <- Brier <- IC0 <- nth.times <- IF.Brier <- lower <- se <- upper <- model <- NF <- IPCW <- response <- reference <- status0 <- IBS <- NULL + riskRegression_status = riskRegression_time <- residuals <- risk <- WTi <- riskRegression_event <- riskRegression_event <- Brier <- IC0 <- nth.times <- IF.Brier <- lower <- se <- upper <- model <- NF <- IPCW <- response <- reference <- riskRegression_status0 <- IBS <- NULL ## sum across bootstrap samples where subject i is out of bag if (cens.type=="rightCensored"){ if (response.type=="survival"){ ## event of interest before times - DT.B[time<=times & status==1,residuals:=(1-risk)^2/WTi] + DT.B[riskRegression_time<=times & riskRegression_status==1,residuals:=(1-risk)^2/WTi] } else{ ## competing risks ## event of interest before times - DT.B[time<=times & status==1 & event==cause,residuals:=(1-risk)^2/WTi] + DT.B[riskRegression_time<=times & riskRegression_status==1 & riskRegression_event==cause,residuals:=(1-risk)^2/WTi] ## competing event before times - DT.B[time<=times & status==1 &event!=cause,residuals:=(0-risk)^2/WTi] + DT.B[riskRegression_time<=times & riskRegression_status==1 &riskRegression_event!=cause,residuals:=(0-risk)^2/WTi] } ## right censored before times - DT.B[time<=times & status==0,residuals:=0] + DT.B[riskRegression_time<=times & riskRegression_status==0,residuals:=0] ## no event at times - DT.B[time>times,residuals:=(risk)^2/Wt] + DT.B[riskRegression_time>times,residuals:=(risk)^2/Wt] }else{ - DT.B[,residuals:=(ReSpOnSe-risk)^2] + DT.B[,residuals:=(riskRegression_event-risk)^2] } ## for each individual sum the residuals of the bootstraps where this individual is out-of-bag ## divide by number of times out-off-bag later @@ -94,7 +94,7 @@ crossvalPerf.loob.Brier <- function(times, if (cens.type[1]=="rightCensored" && !conservative){ # merge with ## FIXME HERE - rr_vars <- grep("^riskRegression_",names(testdata)) + rr_vars <- grep("^riskRegression_",names(data)) DT.B <- data[,rr_vars,with=FALSE][DT.B,,on="riskRegression_ID"] DT.B[,nth.times:=as.numeric(factor(times))] WW <- data.table(riskRegression_ID=1:N,WTi=Weights$IPCW.subject.times,key="riskRegression_ID") @@ -123,20 +123,20 @@ crossvalPerf.loob.Brier <- function(times, #amount of observations in the data, does not fulfill this. #the calculations in getInfluenceCurve.Brier cannot accomodate this (for now). if (response.type == "survival"){ - DT.B[,status0:=status] + DT.B[,riskRegression_status0:=riskRegression_status] } else { - DT.B[,status0:=status*event] + DT.B[,riskRegression_status0:=riskRegression_status*riskRegression_event] } DT.B[,IF.Brier:=getInfluenceCurve.Brier(t=times[1], - time=time, - IC0, + time=riskRegression_time, + IC0 = IC0, residuals=residuals, IC.G=Weights$IC, cens.model=cens.model, conservative = conservative, nth.times=nth.times[1], - event = status0),by=list(model,times)] + event = riskRegression_status0),by=list(model,times)] score.loob <- DT.B[,data.table(Brier=sum(residuals)/N, se=sd(IF.Brier)/sqrt(N)),by=byvars] } @@ -230,8 +230,8 @@ crossvalPerf.loob.Brier <- function(times, DT.B[,model:=factor(model,levels=mlevs,mlabels)] if (all(c("Wt","WTi")%in%names(DT.B))){ DT.B[,IPCW:=1/WTi] - DT.B[time>=times,IPCW:=1/Wt] - DT.B[time=times,IPCW:=1/Wt] + DT.B[riskRegression_time 2){ wdata[,event2 := as.numeric(event)] - wdata[status == 0,event2 := 0] + wdata[riskRegression_status == 0,event2 := 0] wdata[,event2 := NULL] } # else { # fitCSC <- NULL # } - # wdata[,status:=1-status] + # wdata[,riskRegression_status:=1-riskRegression_status] Y <- data[["riskRegression_time"]] - status <- data[["status"]] + riskRegression_status <- data[["riskRegression_status"]] ## fit Cox model for censoring times args <- list(x=TRUE,y=TRUE,eps=0.000001,linear.predictors=TRUE) args$surv <- TRUE @@ -65,7 +65,7 @@ getCensoringWeights <- function(formula, incrementsOfTime <- diff(sort(Y)) minimumIncrement <- min(incrementsOfTime[incrementsOfTime > 0]) #need > 0 in case of ties IPCW.subject.times <- as.numeric(rms::survest(fit,newdata=wdata,times=Y-minimumIncrement/2,what='parallel',se.fit=FALSE)) #Y-minimumIncrement/2 - # subject.times <- Y[(((Y<=max(times))*status)==1)] + # subject.times <- Y[(((Y<=max(times))*riskRegression_status)==1)] out <- list(IPCW.times=IPCW.times,IPCW.subject.times=IPCW.subject.times,method=cens.model) if (influence.curve==TRUE){ @@ -88,50 +88,50 @@ getCensoringWeights <- function(formula, ic.weights <- matrix(0,N,N) ## (i,j)'th entry is f_j(tilde{T}_i-;X_i)/G(tilde{T}_i|X_i) when delta_i != 0 and time_i <= tau; otherwise it is f_j(tau-;X_i)/G(tau | X_i) # k=0 ## counts subject-times with event before t for (i in 1:N){ - if (Y[i]<=times[t.ind] && status[i] != 0){ ## min(T,C)<=t, note that (residuals==0) => (status==0) + if (Y[i]<=times[t.ind] && riskRegression_status[i] != 0){ ## min(T,C)<=t, note that (residuals==0) => (status==0) # k=k+1 ic.weights[i,] <- IC.subject[,1,i] }else if (Y[i] > times[t.ind]){## min(T,C)>t - ic.weights[i,] <- IC.times[,t.ind,i] + ic.weights[i,] <- IC.times[,t.ind,i] } } IC.weights[[t.ind]] <- ic.weights } - IC <- list(censoring.save.memory = censoring.save.memory, IC.weights = IC.weights) + IC <- list(censoring.save.memory = censoring.save.memory, IC.weights = IC.weights) } else { - IC <- list(censoring.save.memory = censoring.save.memory,fit=fit,wdata=wdata, TiMinus = TiMinus) + IC <- list(censoring.save.memory = censoring.save.memory,fit=fit,wdata=wdata, TiMinus = TiMinus) } out <- c(out,list(IC=IC)) } }, "discrete"={ - stop("Use stratification with Cox instead. ") + stop("Use stratification with Cox instead. ") }, { warning("Using other models (than Cox) for getting the censoring weights is under construction. ") vv <- all.vars(formula(delete.response(terms(formula)))) new.formula<-as.formula(paste0("Surv(riskRegression_time,riskRegression_status)",paste0("~",paste0(paste(vv,collapse = "+"))))) wdata <- copy(data) - wdata[,status:=1-status] + wdata[,riskRegression_status:=1-riskRegression_status] input <- list(formula=new.formula,data=wdata) message("Fitting censoring model to data ...", appendLF = FALSE) fit <- do.call(cens.model,input) message("done!") if (influence.curve){ - if (is.null(data[["event"]])){ - new.formula<-as.formula(paste0("Surv(riskRegression_time,riskRegression_status==1)",paste0("~",paste0(paste(vv,collapse = "+"))))) - } - else { - new.formula<-as.formula(paste0("Surv(riskRegression_time,riskRegression_event==1)",paste0("~",paste0(paste(vv,collapse = "+"))))) - } - input <- list(formula=new.formula,data=wdata) - message("Fitting time model to data ...", appendLF = FALSE) - fit.time <- do.call(cens.model,input) - message("done!") + if (is.null(data[["event"]])){ + new.formula<-as.formula(paste0("Surv(riskRegression_time,riskRegression_status==1)",paste0("~",paste0(paste(vv,collapse = "+"))))) + } + else { + new.formula<-as.formula(paste0("Surv(riskRegression_time,riskRegression_event==1)",paste0("~",paste0(paste(vv,collapse = "+"))))) + } + input <- list(formula=new.formula,data=wdata) + message("Fitting time model to data ...", appendLF = FALSE) + fit.time <- do.call(cens.model,input) + message("done!") } else { - fit.time <- NULL + fit.time <- NULL } # IC.data <- list(Stimes = diag(1-predictRisk(fit.time,wdata,wdata$time,1)), Gtimes=diag(1-predictRisk(fit,wdata,wdata$time,1))) diff --git a/R/getInfluenceCurve.R b/R/getInfluenceCurve.R deleted file mode 100644 index 8ded0f0..0000000 --- a/R/getInfluenceCurve.R +++ /dev/null @@ -1,34 +0,0 @@ -getInfluenceCurve.AUC <- function(t, - time, - event, - WTi, - Wt, - risk, - MC, - auc, - nth.times, - conservative, - cens.model){ - ## assert that time is sorted in ascending order with stop - if (is.unsorted(time)){ - stop("Internal error. Time is not sorted in ascending order. ") - } - - conservativeIFcalculation <- getIC0AUC(time,event,t,risk,WTi,Wt,auc) - if (conservative[[1]] || cens.model[[1]] == "none"){ - conservativeIFcalculation[["ic0"]] - } - else { - conservativeIFcalculation[["ic0"]]+getInfluenceFunction.AUC.censoring.term(time = time, - event = event, - t = t, - IFcalculationList = conservativeIFcalculation, - MC = MC, - cens.model = cens.model, - Wt = Wt, - auc = auc, - nth.times = nth.times) - } -} - - diff --git a/R/getResponse.R b/R/getResponse.R index 5a3a14d..cccb32a 100644 --- a/R/getResponse.R +++ b/R/getResponse.R @@ -28,7 +28,7 @@ getResponse <- function(formula,cause,data,vars){ } } ## coercing to 0/1 variable - response <- data.table(ReSpOnSe=as.numeric(response==cause)) + response <- data.table(riskRegression_event=as.numeric(response==cause)) data.table::setattr(response,"states",c("0","1")) data.table::setattr(response,"event","1") data.table::setattr(response,"model","binary") @@ -39,6 +39,7 @@ getResponse <- function(formula,cause,data,vars){ } } else{ + warning("Methods for continuous outcomes are either not\n (not yet) implemented or not well tested.") data.table::setattr(response,"model","continuous") attr(response,"event") <- NULL } diff --git a/R/plotCalibration.R b/R/plotCalibration.R index 13bf65c..5e43b49 100644 --- a/R/plotCalibration.R +++ b/R/plotCalibration.R @@ -3,9 +3,9 @@ ## author: Thomas Alexander Gerds ## created: Feb 23 2017 (11:15) ## Version: -## last-updated: Jun 4 2024 (08:07) +## last-updated: Jun 4 2024 (15:08) ## By: Thomas Alexander Gerds -## Update #: 397 +## Update #: 398 #---------------------------------------------------------------------- ## ### Commentary: @@ -202,7 +202,7 @@ plotCalibration <- function(x, pframe <- x$Calibration$plotframe if (is.null(pframe)) stop("Object has no information for calibration plot.\nYou should call the function \"riskRegression::Score\" with plots=\"calibration\".") - Rvar <- grep("^(ReSpOnSe|pseudovalue)$",names(pframe),value=TRUE) + Rvar <- grep("^(riskRegression_event|pseudovalue)$",names(pframe),value=TRUE) if (!missing(models)){ fitted.models <- pframe[,unique(model)] if (all(models%in%fitted.models)){ diff --git a/R/plotRisk.R b/R/plotRisk.R index 604eb01..ece40b1 100644 --- a/R/plotRisk.R +++ b/R/plotRisk.R @@ -3,9 +3,9 @@ ## Author: Thomas Alexander Gerds ## Created: Mar 13 2017 (16:53) ## Version: -## Last-Updated: Jun 4 2024 (07:22) +## Last-Updated: Jun 4 2024 (19:12) ## By: Thomas Alexander Gerds -## Update #: 241 +## Update #: 243 #---------------------------------------------------------------------- ## ### Commentary: @@ -89,7 +89,7 @@ plotRisk <- function(x, preclipse=0, preclipse.shade=FALSE, ...){ - model = ReSpOnSe = risk = status = event = cause = NULL + model = riskRegression_event = riskRegression_time = risk = status = event = cause = NULL if (is.null(x$risks$score)) stop("No predicted risks in object. You should set summary='risks' when calling Score.") if (!is.null(x$null.model)){ pframe <- x$risks$score[model!=x$null.model] @@ -139,14 +139,14 @@ plotRisk <- function(x, # order according to current cause of interest states <- c(cause,states[cause != states]) if (x$response.type=="binary"){ - Rfactor <- factor(pframe[model==modelnames[1],ReSpOnSe],levels = rev(states),labels = c("No event","Event")) + Rfactor <- factor(pframe[model==modelnames[1],riskRegression_event],levels = rev(states),labels = c("No event","Event")) } else{ if (x$response.type=="survival"){ Rfactor <- pframe[model==modelnames[1], { r <- factor(status,levels = c(-1,0,1),labels = c("No event","Censored","Event")) - r[time>times] <- "No event" + r[riskRegression_time>times] <- "No event" r }] }else{ ## competing risks @@ -160,7 +160,7 @@ plotRisk <- function(x, r = factor(as.character(event),levels = 0:(nCR+1),labels = c("No event",x$states)) } ## r[status == 0] = "Censored" - r[time>times] = "No event" + r[riskRegression_time>times] = "No event" # check if all states are anonymous numbers if (suppressWarnings(sum(is.na(as.numeric(states))) == 0)){ if (nCR == 1){ diff --git a/R/riskQuantile.R b/R/riskQuantile.R index 3a676a2..b1b0bb9 100644 --- a/R/riskQuantile.R +++ b/R/riskQuantile.R @@ -3,9 +3,9 @@ ## author: Thomas Alexander Gerds ## created: Jan 9 2016 (19:31) ## Version: -## last-updated: Jun 4 2024 (07:21) +## last-updated: Jun 4 2024 (15:49) ## By: Thomas Alexander Gerds -## Update #: 307 +## Update #: 309 #---------------------------------------------------------------------- ## ### Commentary: @@ -24,13 +24,13 @@ getQuantile <- function(x,Fx,Q){ } riskQuantile.binary <- function(DT,N,NT,NF,dolist,Q,...){ - reference=model=ReSpOnSe=risk=cause=X=riskRegression_ID=NULL + reference=model=riskRegression_event=risk=cause=X=riskRegression_ID=NULL models <- unique(DT[,model]) if (missing(Q)) Q <- c(0.05,0.25,0.5,0.75,0.95) else Q <- sort(Q) - score.event <- DT[ReSpOnSe==1,data.table(t(quantile(risk,probs=Q))),by=list(model)] + score.event <- DT[riskRegression_event==1,data.table(t(quantile(risk,probs=Q))),by=list(model)] score.event[,cause:="event"] - score.eventfree <- DT[ReSpOnSe==0,data.table(t(quantile(risk,probs=Q))),by=list(model)] + score.eventfree <- DT[riskRegression_event==0,data.table(t(quantile(risk,probs=Q))),by=list(model)] score.eventfree[,cause:="event-free"] score.overall <- DT[,data.table(t(quantile(risk,probs=Q))),by=list(model)] score.overall[,cause:="overall"] @@ -48,9 +48,9 @@ riskQuantile.binary <- function(DT,N,NT,NF,dolist,Q,...){ N <- NROW(DTdiff) Xrange <- DTdiff[,range(X)] Xmed <- DTdiff[,median(X)] - changedist.event <- DTdiff[ReSpOnSe==1,data.table(t(quantile(X,probs=Q))),by=list(model)] + changedist.event <- DTdiff[riskRegression_event==1,data.table(t(quantile(X,probs=Q))),by=list(model)] changedist.event[,cause:="event"] - changedist.eventfree <- DTdiff[ReSpOnSe==0,data.table(t(quantile(X,probs=Q))),by=list(model)] + changedist.eventfree <- DTdiff[riskRegression_event==0,data.table(t(quantile(X,probs=Q))),by=list(model)] changedist.eventfree[,cause:="event-free"] changedist.overall <- DTdiff[,data.table(t(quantile(X,probs=Q))),by=list(model)] changedist.overall[,cause:="overall"] @@ -105,7 +105,7 @@ riskQuantile.survival <- function(DT,N,NT,NF,dolist,Q,...){ ## ## For 'event-free analyses' P(X<=x|T>t) is estimated by P(T>t|X<=x) P(X<=x)/P(T>t) ####### - surv <- DT[model==models[[1]],data.table::data.table("surv"=(1/N*sum((time>times)/Wt))),by=times] + surv <- DT[model==models[[1]],data.table::data.table("surv"=(1/N*sum((riskRegression_time>times)/Wt))),by=times] surv[,cuminc:=1-surv] getQ.event <- function(Q,tp,X,time,status,WTi,surv){ uX <- sort(unique(X[time<=tp & status==1])) @@ -125,11 +125,11 @@ riskQuantile.survival <- function(DT,N,NT,NF,dolist,Q,...){ } ## a <- DT[model==1] ## system.time(a[,getQ.eventFree(Q=Q,tp=times[1],X=risk,time=time,Wt=Wt,surv=surv)]) - score.eventfree <- DT[,getQ.eventFree(Q=Q,tp=times,X=risk,time=time,Wt=Wt,surv=surv),by=list(model,times)] + score.eventfree <- DT[,getQ.eventFree(Q=Q,tp=times,X=risk,time=riskRegression_time,Wt=Wt,surv=surv),by=list(model,times)] ## setkey(DT,model,times) ## save(surv,file="~/tmp/surv.rda") ## save(DT,file="~/tmp/DT.rda") - score.event <- DT[,getQ.event(Q=Q,tp=times,X=risk,time=time,status=status,WTi=WTi,surv=surv),by=list(model,times)] + score.event <- DT[,getQ.event(Q=Q,tp=times,X=risk,time=riskRegression_time,status=riskRegression_status,WTi=WTi,surv=surv),by=list(model,times)] score.overall <- DT[,data.table(t(quantile(risk,probs=Q))),by=list(model,times)] score.overall[,cause:="overall"] colnames(score.overall) <- colnames(score.event) @@ -147,8 +147,8 @@ riskQuantile.survival <- function(DT,N,NT,NF,dolist,Q,...){ N <- NROW(DTdiff) Xrange <- DTdiff[,range(X)] Xmed <- DTdiff[,median(X)] - changedist.eventfree <- DTdiff[,getQ.eventFree(Q=Q,tp=times,X=X,time=time,Wt=Wt,surv=surv),by=list(model,times)] - changedist.event <- DTdiff[,getQ.event(Q=Q,tp=times,X=X,time=time,status=status,WTi=WTi,surv=surv),by=list(model,times)] + changedist.eventfree <- DTdiff[,getQ.eventFree(Q=Q,tp=times,X=X,time=riskRegression_time,Wt=Wt,surv=surv),by=list(model,times)] + changedist.event <- DTdiff[,getQ.event(Q=Q,tp=times,X=X,time=riskRegression_time,status=riskRegression_status,WTi=WTi,surv=surv),by=list(model,times)] changedist.overall <- DTdiff[,data.table(t(quantile(X,probs=Q))),by=list(model,times)] changedist.overall[,cause:="overall"] colnames(changedist.overall) <- colnames(changedist.event) @@ -211,8 +211,8 @@ riskQuantile.competing.risks <- function(DT,N,NT,NF,dolist,cause,states,Q,...){ ## ## For 'event-free analyses' P(X<=x|T>t) is estimated by P(T>t|X<=x) P(X<=x)/P(T>t) ####### - surv <- DT[model==models[[1]],data.table::data.table("surv"=(1/N*sum((time>times)/Wt))),by=times] - cuminc <- lapply(states.code,function(cc){DT[model==models[[1]],data.table::data.table("cuminc"=1/N*sum((event==cc & time<=times)/WTi)),by=times]}) + surv <- DT[model==models[[1]],data.table::data.table("surv"=(1/N*sum((riskRegression_time>times)/Wt))),by=times] + cuminc <- lapply(states.code,function(cc){DT[model==models[[1]],data.table::data.table("cuminc"=1/N*sum((riskRegression_event==cc & riskRegression_time<=times)/WTi)),by=times]}) names(cuminc) <- states getQ.states <- function(Q,tp,X,time,event,WTi,cuminc,states.code){ uX <- sort(unique(X)) @@ -234,8 +234,8 @@ riskQuantile.competing.risks <- function(DT,N,NT,NF,dolist,cause,states,Q,...){ qR[,cause:="event-free"] qR } - score.eventfree <- DT[,getQ.eventFree(Q=Q,tp=times,X=risk,time=time,Wt=Wt,surv=surv),by=list(model,times)] - score.states <- DT[,getQ.states(Q=Q,tp=times,X=risk,time=time,event=event,WTi=WTi,cuminc=cuminc,states.code=states.code),by=list(model,times)] + score.eventfree <- DT[,getQ.eventFree(Q=Q,tp=times,X=risk,time=riskRegression_time,Wt=Wt,surv=surv),by=list(model,times)] + score.states <- DT[,getQ.states(Q=Q,tp=times,X=risk,time=riskRegression_time,event=riskRegression_event,WTi=WTi,cuminc=cuminc,states.code=states.code),by=list(model,times)] score.overall <- DT[,data.table(t(quantile(risk,probs=Q))),by=list(model,times)] score.overall[,cause:="overall"] colnames(score.overall) <- colnames(score.states) @@ -253,8 +253,8 @@ riskQuantile.competing.risks <- function(DT,N,NT,NF,dolist,cause,states,Q,...){ N <- NROW(DTdiff) Xrange <- DTdiff[,range(X)] Xmed <- DTdiff[,median(X)] - changedist.eventfree <- DTdiff[,getQ.eventFree(Q=Q,tp=times,X=X,time=time,Wt=Wt,surv=surv),by=list(model,times)] - changedist.states <- DTdiff[,getQ.states(Q=Q,tp=times,X=X,time=time,event=event,WTi=WTi,cuminc=cuminc,states.code=states.code),by=list(model,times)] + changedist.eventfree <- DTdiff[,getQ.eventFree(Q=Q,tp=times,X=X,time=riskRegression_time,Wt=Wt,surv=surv),by=list(model,times)] + changedist.states <- DTdiff[,getQ.states(Q=Q,tp=times,X=X,time=riskRegression_time,event=riskRegression_event,WTi=WTi,cuminc=cuminc,states.code=states.code),by=list(model,times)] changedist.overall <- DTdiff[,data.table(t(quantile(X,probs=Q))),by=list(model,times)] changedist.overall[,cause:="overall"] colnames(changedist.overall) <- colnames(changedist.states) diff --git a/man/Score.Rd b/man/Score.Rd index 629ce96..f59bb21 100644 --- a/man/Score.Rd +++ b/man/Score.Rd @@ -286,7 +286,6 @@ replacement from the data, which does not depend on the sample size: \code{N=300} \code{mean(sapply(1:B, function(b){match(1,sample(1:N,size=N,replace=TRUE),nomatch=0)}))} - ## Bootstrap without replacement (training size set to be 70 percent of data) B=10, M=.7 @@ -317,7 +316,9 @@ testdat <- sampleData(40,outcome="binary") ## score logistic regression models lr1 = glm(Y~X1+X2+X7+X9,data=learndat,family=binomial) lr2 = glm(Y~X3+X5,data=learndat,family=binomial) -Score(list("LR(X1+X2+X7+X9)"=lr1,"LR(X3+X5)"=lr2),formula=Y~1,data=testdat) +x=Score(list("LR(X1+X2+X7+X9)"=lr1,"LR(X3+X5)"=lr2),formula=Y~1,data=testdat) +print(x) +summary(x) ## ROC curve and calibration plot xb=Score(list("LR(X1+X2+X7+X9)"=lr1,"LR(X3+X5+X6)"=lr2),formula=Y~1, diff --git a/man/plotCalibration.Rd b/man/plotCalibration.Rd index b2b5034..b72c651 100644 --- a/man/plotCalibration.Rd +++ b/man/plotCalibration.Rd @@ -10,7 +10,7 @@ plotCalibration( times, method = "nne", cens.method, - round = TRUE, + round = 2, bandwidth = NULL, q = 10, bars = FALSE,