repairing plotRisk

tagteam · Jun 22, 2024 · 9cabac2 · 9cabac2
1 parent 7a9f154
commit 9cabac2
Show file tree

Hide file tree

Showing 20 changed files with 40 additions and 2,996 deletions.
diff --git a/R/Score.R b/R/Score.R
@@ -575,10 +575,10 @@ Score.list <- function(object,
     if ((!missing(cutpoints)) && (!missing(split.method)) && tolower(split.method) != "none")
         if (verbose>0) warning("Calculation of sensitivities, specificities and predictive values only implemented for argument split.method='none'.\n
                 Argument 'cutpoints' is ignored.")
-    if ("Calibration" %in% plots) {
+    ## if ("Calibration" %in% plots) {
         ## add pseudo if needed
-        if (!("pseudo" %in% cens.method)) cens.method <- c(cens.method,"pseudo")
-    }
+        ## if (!("pseudo" %in% cens.method)) cens.method <- c(cens.method,"pseudo")
+    ## }
 
     # }}}
     # {{{ censoring model arguments
@@ -922,10 +922,10 @@ c.f., Chapter 7, Section 5 in Gerds & Kattan 2021. Medical risk prediction model
                 # need to communicate the censoring code of the event variable
                 # produced by Hist via model.frame in case of competing risks
                 if (response.type=="competing.risks"){
-                    censcode <- data[status==0,event[1]]
-                    margForm <- Hist(time,event,cens.code=censcode)~1
+                    censcode <- data[riskRegression_status==0,riskRegression_event[1]]
+                    margForm <- Hist(riskRegression_time,riskRegression_event,cens.code=censcode)~1
                 }else{
-                    censcode <- data[status==0,status[1]]
+                    censcode <- data[riskRegression_status==0,riskRegression_status[1]]
                     margForm <- update(formula,".~1")
                 }
                 margFit <- prodlim::prodlim(margForm,data=data)
@@ -935,6 +935,8 @@ c.f., Chapter 7, Section 5 in Gerds & Kattan 2021. Medical risk prediction model
                                    pseudovalue=c(prodlim::jackknife(margFit,cause=position.cause,times=times)))
                 if (response.type=="survival") jack[,pseudovalue:=1-pseudovalue]
             }
+        }else{
+            jack <- NULL
         }
     } else{
         Weights <- NULL

diff --git 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  5 2024 (07:27) 
+## Last-Updated: Jun 17 2024 (10:15) 
 ##           By: Thomas Alexander Gerds
-##     Update #: 23
+##     Update #: 24
 #----------------------------------------------------------------------
 ##
 ### Commentary:
@@ -84,12 +84,10 @@ computePerformance <- function(DT,
     # }}}
     # {{{ collect data for calibration plots
     if ("Calibration" %in% plots){
-        if (response.type[[1]]=="binary" || cens.type[[1]]=="uncensored")
-            out[["Calibration"]]$plotframe <- DT[model!=0]
-        else{
-            out[["Calibration"]]$plotframe <- merge(jack,DT[model!=0],by=c("riskRegression_ID","times"))
-        }
+        out[["Calibration"]]$plotframe <- DT[model!=0]
         out[["Calibration"]]$plotframe[,model:=factor(model,levels=models$levels,labels=models$labels)]
+        if (length(jack)>0)
+            out[["Calibration"]]$plotframe <- merge(jack,DT[model!=0],by=c("riskRegression_ID","times"))
     }
     # }}}
     ## make sure that Brier score comes first, so that we can remove the null.model afterwards

diff --git 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 (19:12) 
+## Last-Updated: Jun 22 2024 (08:35) 
 ##           By: Thomas Alexander Gerds
-##     Update #: 243
+##     Update #: 249
 #----------------------------------------------------------------------
 ## 
 ### Commentary: 
@@ -89,7 +89,7 @@ plotRisk <- function(x,
                      preclipse=0,
                      preclipse.shade=FALSE,
                      ...){
-    model = riskRegression_event = riskRegression_time = risk = status = event = cause = NULL
+    model = riskRegression_event = riskRegression_time = risk = riskRegression_status = 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]
@@ -112,9 +112,9 @@ plotRisk <- function(x,
     if (missing(models)){
         models <- pframe[,levels(model)[1:2]]
     }else{
-        if (any(!(models %in% unique(pframe$models)))){
+        if (any(!(models %in% unique(pframe$model)))){
             stop(paste0("Fitted object does not contain models named: ",paste0(models,collapse = ", "),
-                        "\nAvailable models are named: ",paste0(unique(pframe$models),collapse = ", ")))
+                        "\nAvailable models are named: ",paste0(unique(pframe$model),collapse = ", ")))
         }
         if (is.na(models[2])) stop("Need two models for a scatterplot of predicted risks")
         fitted.models <- x$models
@@ -145,7 +145,7 @@ plotRisk <- function(x,
         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 <- factor(riskRegression_status,levels = c(-1,0,1),labels = c("No event","Censored","Event"))
                 r[riskRegression_time>times] <- "No event"
                 r
             }]
@@ -154,17 +154,16 @@ plotRisk <- function(x,
             # sort such that event-free, censored, current cause, cr1, cr2, ...
             Rfactor <- pframe[model==modelnames[1],{
                 # need to use x$states as object$states is always sorted no matter the cause of interest
-                if (any(status == 0)){
-                    r = factor(as.character(event),levels = 0:(nCR+2),labels = c("No event",x$states,"Censored"))
+                if (any(riskRegression_status == 0)){
+                    r = factor(as.character(riskRegression_event),levels = 0:(nCR+2),labels = c("No event",x$states,"Censored"))
                 } else{
-                    r = factor(as.character(event),levels = 0:(nCR+1),labels = c("No event",x$states))
+                    r = factor(as.character(riskRegression_event),levels = 0:(nCR+1),labels = c("No event",x$states))
                 }
-                ## r[status == 0] = "Censored"
                 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){
-                        if (any(status == 0))
+                        if (any(riskRegression_status == 0))
                             r = factor(r,
                                        levels = c("No event","Censored",states),
                                        labels = c("No event","Censored","Event","Competing risk"))
@@ -173,7 +172,7 @@ plotRisk <- function(x,
                                        levels = c("No event",states),
                                        labels = c("No event","Event","Competing risk"))
                     }else{
-                        if (any(status == 0))
+                        if (any(riskRegression_status == 0))
                             r = factor(r,
                                        levels = c("No event","Censored",states),
                                        labels = c("No event","Censored","Event",paste0("Competing risk ",1:nCR)))
@@ -246,8 +245,14 @@ plotRisk <- function(x,
                 paste(paste(this.legend),collapse = "\n"))
     }
     if (x$response.type == "competing.risks") {
-        message("Counts of events and censored are evaluated at the prediction time horizon (times=",times,"):\n\n",
-                paste(this.legend,collapse = "\n"),"\n\nwhere the event type values (",paste(states,collapse = ", "),") in the data correspond to labels:\n Event, ",ifelse(nCR == 1,"Competing risk",paste0("Competing risk ",1:nCR,collapse = ", "))
+        message("Counts of events and censored are evaluated at the prediction time horizon (times=",
+                times,
+                "):\n\n",
+                paste(this.legend,collapse = "\n"),
+                "\n\nwhere the event type values (",
+                paste(states,collapse = ", "),
+                ") in the data correspond to labels:\n Event, ",
+                ifelse(nCR == 1,"Competing risk",paste0("Competing risk ",1:nCR,collapse = ", "))
                 )
     }
     legend.DefaultArgs <- list(legend=this.legend,

diff --git a/slowtests/slowtest-score.R b/slowtests/slowtest-score.R
@@ -3,9 +3,9 @@
 ## Author: Thomas Alexander Gerds
 ## Created: Dec  6 2020 (09:25) 
 ## Version: 
-## Last-Updated: Jun  5 2024 (18:13) 
+## Last-Updated: Jun 17 2024 (10:58) 
 ##           By: Thomas Alexander Gerds
-##     Update #: 13
+##     Update #: 15
 #----------------------------------------------------------------------
 ## 
 ### Commentary: 
@@ -83,7 +83,6 @@ test_that("competing risks outcome: robustness against order of data set",{
     expect_equal(ignore_attr=TRUE,s3$AUC,s3b$AUC)
 })
 # }}}
-
 # {{{ "survival outcome: Brier Score pec vs Score"
 if (requireNamespace("pec",quietly=TRUE)){
     test_that("survival outcome: Brier Score pec vs Score",{
@@ -100,7 +99,6 @@ if (requireNamespace("pec",quietly=TRUE)){
     })
 }
 # }}}
-
 # {{{ "survival outcome: matrix input"
 test_that("survival outcome: matrix input",{
     set.seed(112)
@@ -145,7 +143,7 @@ test_that("Number of models and time points", {
         ## expect_equal(ignore_attr=TRUE,r1$AUC$score[model=="a"],r2$AUC$score[model=="a" & times==1000])
     }
 })
-
+# }}}
 # {{{ "Bootstrap cross validation
 test_that("Number of models and time points", {
     if (!requireNamespace("pec",quietly=TRUE)){
@@ -180,7 +178,6 @@ test_that("Number of models and time points", {
     }
 })
 # }}}
-
 # {{{ "LOOB: Number of models and time points"
 test_that("LOOB: Number of models and time points", {   
     library(testthat)
@@ -228,9 +225,7 @@ test_that("LOOB: Number of models and time points", {
     expect_equal(ignore_attr=TRUE,A1$Brier$score[model=="fit1"&times==365.25*4],
                  A$Brier$score[model=="fit1"&times==365.25*4])
 })
-
-
-
+# }}}
 # {{{ "vcov AUC"
 test_that("vcov AUC",{
     set.seed(112)
@@ -248,7 +243,7 @@ test_that("vcov AUC",{
     expect_equal(ignore_attr=TRUE,dim(test$AUC$vcov),c(4,4))
 })
 # }}}
-
+# {{{ LOOB binary
 ## GIVES WARNING 
 ## already exporting variable(s): data, split.method, Weights, N, trainseeds
 test_that("loob binary",{
@@ -275,7 +270,8 @@ test_that("loob binary",{
     expect_equal(ignore_attr=TRUE,loob.se0$AUC$contrasts$delta,loob.se1$AUC$contrasts$delta)
     expect_equal(ignore_attr=TRUE,loob.se0$Brier$contrasts$delta,loob.se1$Brier$contrasts$delta)
 })
-
+# }}}
+# {{{ bootcv binary
 test_that("bootcv binary (multi.state.test)",{
     learndat=sampleData(200,outcome="binary")
     lr1a = glm(Y~X6,data=learndat,family=binomial)
@@ -302,7 +298,7 @@ test_that("bootcv binary (multi.state.test)",{
     ## expect_equal(ignore_attr=TRUE,bootcv[[3]][[m]]$contrasts[,.(p)],bootcv[[4]][[m]]$contrasts[,.(p)])
     ## }
 })
-
+# }}}
 # {{{ "Brier score"
 if (requireNamespace("pec",quietly=TRUE)){
     test_that("Brier score",{

diff --git a/tests/testthat/test-AUC.R b/tests/testthat/test-AUC.R
diff --git a/tests/testthat/test-BinomialRegression.R b/tests/testthat/test-BinomialRegression.R
diff --git a/tests/testthat/test-BrierScore.R b/tests/testthat/test-BrierScore.R