fix CRAN checks

bozenne · Jul 18, 2023 · 0447c42 · 0447c42
1 parent 81e9bd5
commit 0447c42
Show file tree

Hide file tree

Showing 23 changed files with 197 additions and 108 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -2,7 +2,7 @@ Package: BuyseTest
 Type: Package
 Title: Generalized Pairwise Comparisons
 Version: 3.0.0
-Date: 2023-07-17
+Date: 2023-07-18
 Authors@R: c(
     person("Brice", "Ozenne", role = c("aut", "cre"), email = "brice.mh.ozenne@gmail.com", comment = c(ORCID = "0000-0001-9694-2956")),
     person("Julien", "Peron", role = "ctb"),

diff --git a/NEWS b/NEWS
@@ -12,6 +12,7 @@
   function, now the seed is used to generate sample or simulation
   specific seeds. Results are now identical whether one or more CPUs
   are used, and reproducible when multiple CPUs are used.
+- add some support for paired data via strata argument in BuyseTest
 
 ** Internal change
 - more uniform effect of the strata argument

diff --git a/R/BuyseTTEM.R b/R/BuyseTTEM.R
@@ -3,9 +3,9 @@
 ## Author: Brice Ozenne
 ## Created: nov 18 2020 (12:15) 
 ## Version: 
-## Last-Updated: jun 27 2023 (13:54) 
+## Last-Updated: jul 18 2023 (10:18) 
 ##           By: Brice Ozenne
-##     Update #: 639
+##     Update #: 665
 ##----------------------------------------------------------------------
 ## 
 ### Commentary: 
@@ -127,16 +127,22 @@ BuyseTTEM.prodlim <- function(object, treatment, iid, iid.surv = "exp", ...){
     level.strata <- object$peron$level.strata
     n.strata <- object$peron$n.strata
     strata.var <- object$peron$strata.var
-
+    object.stratified <- NCOL(object$X)>1
+
     ## ** compute start/stop indexes per strata
     iIndexX.C <- which(X[[treatment]]==0)
     iIndexX.T <- which(X[[treatment]]==1)
 
     ## position in the results of the pair (treatment,strata)
     iMindexX <- do.call(rbind,lapply(1:n.strata, function(iStrata){
-        iIndexS <- which(X[["..strata.."]]==iStrata)
-        return(c(C = intersect(iIndexX.C,iIndexS),
-                 T = intersect(iIndexX.T,iIndexS)))
+        if(object.stratified){
+            iIndexS <- which(X[["..strata.."]]==iStrata)
+            return(c(C = intersect(iIndexX.C,iIndexS),
+                     T = intersect(iIndexX.T,iIndexS)))
+        }else{
+            return(c(C = iIndexX.C,
+                     T = iIndexX.T))
+        }        
     }))
 
     index.start <- matrix(NA, nrow = n.strata, ncol = 2, dimnames = list(NULL,level.treatment))
@@ -210,7 +216,6 @@ BuyseTTEM.prodlim <- function(object, treatment, iid, iid.surv = "exp", ...){
     ## cif <- object$peron$cif[[1]][[1]][[1]]
     ## (cif[cif[,"index.cif.after"]+1,"cif"] - cif[cif[,"index.cif.before"]+1,"cif"]) - cif[,"dcif"]
 
-
     ## ** table iid
     if(iid){
         n.obs <- NROW(object$model.response)
@@ -244,8 +249,13 @@ BuyseTTEM.prodlim <- function(object, treatment, iid, iid.surv = "exp", ...){
                 object$peron$iid.hazard[[iStrata]][[iTreat]] <- vector(mode = "list", length=n.CR)
                 object$peron$iid.cif[[iStrata]][[iTreat]] <- vector(mode = "list", length=n.CR)
 
-                iIndStrata <- intersect(which(model.matrix[[treatment]]==level.treatment[iTreat]),
-                                        which(model.matrix[["..strata.."]]==iStrata))
+                if(object.stratified){    
+                    iIndStrata <- intersect(which(model.matrix[[treatment]]==level.treatment[iTreat]),
+                                            which(model.matrix[["..strata.."]]==iStrata))
+                }else{
+                    iIndStrata <- which(model.matrix[[treatment]]==level.treatment[iTreat])
+                }
+
                 iIndex.jump <- object$peron$jumpSurvHaz[[iStrata]][[iTreat]]$index.jump
                 iN.jump <- length(iIndex.jump)
 
@@ -332,6 +342,7 @@ BuyseTTEM.survreg <- function(object, treatment, n.grid = 1e3, iid, ...){
     treatment <- object$peron$treatment
     level.treatment <- object$peron$level.treatment
     level.strata <- object$peron$level.strata
+    object.stratified <- NCOL(object$X)>1
 
     if(is.null(object$xlevels[[treatment]])){
         mf[[treatment]] <- factor(mf[[treatment]], levels = sort(unique(mf[[treatment]])), labels = level.treatment)
@@ -347,7 +358,6 @@ BuyseTTEM.survreg <- function(object, treatment, n.grid = 1e3, iid, ...){
     object$peron$n.CR <- 1 ## survival case (one type of event)
 
     ## ** prepare for iid
-
     if(iid){
 
         ## *** extract information
@@ -412,11 +422,16 @@ BuyseTTEM.survreg <- function(object, treatment, n.grid = 1e3, iid, ...){
     for(iStrata in 1:n.strata){ ## iStrata <- 1
         for(iTreat in level.treatment){ ## iTreat <- 1
 
-            iIndex.obs <- intersect(
-                intersect(which(mf[[treatment]]==iTreat),
-                          which(object$peron$X[,"..strata.."]==iStrata)),
-                which(mf[,1][,2]==1)
-            )
+            if(object.stratified){
+                iIndex.obs <- intersect(
+                    intersect(which(mf[[treatment]]==iTreat),
+                              which(object$peron$X[,"..strata.."]==iStrata)),
+                    which(mf[,1][,2]==1)
+                )
+            }else{
+                iIndex.obs <- intersect(which(mf[[treatment]]==iTreat),
+                                        which(mf[,1][,2]==1))
+            }
 
             ## jump time  in this strata
             iNewdata <- mf[iIndex.obs[1],,drop=FALSE]

diff --git a/R/BuyseTest-Peron.R b/R/BuyseTest-Peron.R
@@ -3,9 +3,9 @@
 ## Author: Brice Ozenne
 ## Created: okt 12 2020 (11:10) 
 ## Version: 
-## Last-Updated: jun 27 2023 (12:58) 
+## Last-Updated: jul 18 2023 (12:03) 
 ##           By: Brice Ozenne
-##     Update #: 527
+##     Update #: 548
 ##----------------------------------------------------------------------
 ## 
 ### Commentary: 
@@ -19,7 +19,7 @@
 #' @noRd
 calcPeron <- function(data,
                       model.tte, fitter, args,
-                      method.score,
+                      method.score, paired,
                       treatment,
                       level.treatment,
                       endpoint,
@@ -46,17 +46,22 @@ calcPeron <- function(data,
     test.CR <- setNames(vector(mode = "logical", length = D.UTTE), endpoint.UTTE)
 
     ## prepare formula
-    if(is.null(model.tte)){        
+    if(is.null(model.tte)){
         model.tte <- vector(length = D.UTTE, mode = "list")
         names(model.tte) <- endpoint.UTTE
         tofit <- TRUE
         if(length(args)==0){args <- NULL}
-
-        if(any(fitter=="prodlim")){
-            txt.modelUTTE <- paste0("prodlim::Hist(",endpoint.UTTE,",",status.UTTE,") ~ ",treatment," + ..strata..")
-        }else if(any(fitter=="survreg")){
-            txt.modelUTTE <- paste0("survival::Surv(",endpoint.UTTE,",",status.UTTE,") ~ ",treatment," + ..strata..")
+
+        txt.fitter <- sapply(fitter, switch,
+                             "prodlim" = "prodlim::Hist",
+                             "survreg" = "survival::Surv",
+                             NA)
+        if(paired){
+            txt.modelUTTE <- paste0(txt.fitter,"(",endpoint.UTTE,",",status.UTTE,") ~ ",treatment)        
+        }else{
+            txt.modelUTTE <- paste0(txt.fitter,"(",endpoint.UTTE,",",status.UTTE,") ~ ",treatment," + ..strata..")        
         }
+
     }else{
         tofit <- FALSE
     }
@@ -79,7 +84,6 @@ calcPeron <- function(data,
             }
 
         }
-
         model.tte[[iUTTE]] <- BuyseTTEM(model.tte[[iUTTE]], treatment = treatment, level.treatment = level.treatment, level.strata = level.strata, iid = iidNuisance)
     }
 

diff --git a/R/BuyseTest-check.R b/R/BuyseTest-check.R
@@ -3,9 +3,9 @@
 ## Author: Brice Ozenne
 ## Created: apr 27 2018 (23:32) 
 ## Version: 
-## Last-Updated: jul  4 2023 (17:03) 
+## Last-Updated: jul 18 2023 (09:34) 
 ##           By: Brice Ozenne
-##     Update #: 342
+##     Update #: 344
 ##----------------------------------------------------------------------
 ## 
 ### Commentary: 
@@ -151,8 +151,8 @@ testArgs <- function(name.call,
     }else{
         strata.tempo <- data[[strata]]
         if(is.factor(strata.tempo)){strata.tempo <- droplevels(strata.tempo)} ## otherwise the next tapply statement generates NA when there are empty levels which leads to an error
-
-        if(any(sapply(Ustatus.TTE, function(iS){tapply(data[[iS]], strata.tempo, function(iVec){sum(iVec!=0)})})==0)){
+        ## if non-paired data (i.e. more than 2 obs per strata)
+        if(any(table(strata.tempo)>2) && any(sapply(Ustatus.TTE, function(iS){tapply(data[[iS]], strata.tempo, function(iVec){sum(iVec!=0)})})==0)){
             warning("BuyseTest: time to event variables with only censored events in at least one strata \n")
         }
     }

diff --git a/R/BuyseTest-initialization.R b/R/BuyseTest-initialization.R
@@ -65,7 +65,6 @@ initializeArgs <- function(status,
     if(is.null(cpus)){ cpus <- option$cpus }
     if(is.null(keep.pairScore)){ keep.pairScore <- option$keep.pairScore }
     if(is.null(scoring.rule)){ scoring.rule <- option$scoring.rule }
-    if(is.null(pool.strata)){ pool.strata <- option$pool.strata }
     if(is.null(hierarchical)){ hierarchical <- option$hierarchical }
     if(is.null(correction.uninf)){ correction.uninf <- option$correction.uninf }
     if(is.null(method.inference)){ method.inference <- option$method.inference }
@@ -220,7 +219,22 @@ initializeArgs <- function(status,
     ## ** pool.strata
     if(is.null(strata)){
         pool.strata <- 0
-        attr(pool.strata,"original") <- "none"
+        attr(pool.strata,"type") <- "none"
+        attr(pool.strata,"original") <- NA
+    }else if(is.null(pool.strata)){
+        pool.strata <- switch(tolower(option$pool.strata),
+                              "buyse" = 0,
+                              "cmh" = 1,
+                              "equal" = 2,
+                              "var-favorable" = 3.1,
+                              "var-unfavorable" = 3.2,
+                              "var-netbenefit" = 3.3,
+                              "var-winratio" = 3.4,
+                              NA
+                              )
+        attr(pool.strata,"type") <- option$pool.strata   
+        attr(pool.strata,"original") <- NA   
+
     }else if(is.character(pool.strata)){
         pool.strata_save <- tolower(pool.strata)
         pool.strata <- switch(pool.strata_save,
@@ -233,6 +247,7 @@ initializeArgs <- function(status,
                                "var-winratio" = 3.4,
                                NA
                               )
+        attr(pool.strata,"type") <- pool.strata_save
         attr(pool.strata,"original") <- pool.strata_save
     }else{
         pool.strata <- NA
@@ -365,7 +380,7 @@ initializeArgs <- function(status,
 }
 
 ## * initializeData
-initializeData <- function(data, type, endpoint, Uendpoint, D, scoring.rule, status, Ustatus, method.inference, censoring, strata, treatment, hierarchical, copy,
+initializeData <- function(data, type, endpoint, Uendpoint, D, scoring.rule, status, Ustatus, method.inference, censoring, strata, pool.strata, treatment, hierarchical, copy,
                            keep.pairScore, endpoint.TTE, status.TTE, iidNuisance, weightEndpoint, weightObs){
 
     if (!data.table::is.data.table(data)) {
@@ -443,7 +458,7 @@ initializeData <- function(data, type, endpoint, Uendpoint, D, scoring.rule, sta
         D.UTTE <- 0
         index.UTTE <- rep(-100, D)
     }
-    
+
     ## ** scoring method for each endpoint
     ## check if status
     n.CR <- sapply(status, function(iC){max(data[[iC]])})
@@ -471,6 +486,7 @@ initializeData <- function(data, type, endpoint, Uendpoint, D, scoring.rule, sta
     })
     attr(method.score,"test.censoring") <- test.censoring
     attr(method.score,"test.CR") <- test.CR
+    paired <- all(n.obsStrata==2)
 
     ## ** previously analyzed distinct TTE endpoints
     if((scoring.rule==1) && hierarchical){ ## only relevant when using Peron scoring rule with hierarchical GPC
@@ -528,6 +544,17 @@ initializeData <- function(data, type, endpoint, Uendpoint, D, scoring.rule, sta
         weightObs <- data$..weight..
     }
 
+    ## ** pool.strata
+    ## set default pool.strata to Buyse for paired data
+    ## otherwise pooling will do something strange
+    if(paired && pool.strata !=0){
+        if(is.na(attr(pool.strata,"original"))){
+            pool.strata[] <- 0
+        }else{
+            warning("Weights from the \"buyse\" pooling scheme (argument \'pool.strata\') are recommended for paired data. \n")
+        }
+    }
+
     ## ** export
     keep.cols <- union(c(treatment, "..strata.."),
                        na.omit(attr(method.inference,"resampling-strata")))
@@ -540,8 +567,8 @@ initializeData <- function(data, type, endpoint, Uendpoint, D, scoring.rule, sta
                 weightObs = weightObs,
                 index.strata = tapply(data[["..rowIndex.."]], data[["..strata.."]], list),
                 level.treatment = level.treatment,
-                level.strata = level.strata,
-                method.score = method.score,
+                level.strata = level.strata, pool.strata = pool.strata,
+                method.score = method.score, paired = paired,
                 n.strata = n.strata,
                 n.obs = n.obs,
                 n.obsStrata = n.obsStrata,
@@ -554,7 +581,7 @@ initializeData <- function(data, type, endpoint, Uendpoint, D, scoring.rule, sta
                 endpoint.UTTE = endpoint.UTTE,
                 status.UTTE = status.UTTE,
                 D.UTTE = D.UTTE,
-                index.UTTE = index.UTTE,
+                index.UTTE = index.UTTE,                
                 keep.pairScore = keep.pairScore
                 ))
 }

diff --git a/R/BuyseTest-print.R b/R/BuyseTest-print.R
@@ -17,6 +17,7 @@ printGeneral <- function(status,
                          level.treatment,
                          scoring.rule,
                          M.status,
+                         method.score, paired,
                          neutral.as.uninf,
                          correction.uninf,
                          operator,
@@ -29,7 +30,7 @@ printGeneral <- function(status,
                          weightEndpoint,
                          Wscheme,
                          ...){
-
+    
     if(!is.null(strata)){
         n.strata <- length(level.strata)
     }else{
@@ -94,8 +95,12 @@ printGeneral <- function(status,
     cat("   - 2 groups  ",if(D>1){" "},": Control = ",level.treatment[1]," and Treatment = ",level.treatment[2],"\n", sep = "")
     cat("   - ",D," endpoint",if(D>1){"s"},": \n", sep = "")
     print(df.endpoint, row.names = FALSE, quote = FALSE, right = FALSE)
-
-    if(n.strata>1){
+    if(paired){
+        txt.variable <- switch(as.character(length(strata)),
+                               "1" = "variable",
+                               "variables")        
+        cat("   - ", n.strata, " pairs (",txt.variable,": ",paste(strata, collapse = " "),") \n", sep = "")
+    }else if(n.strata>1){
         txt.variable <- switch(as.character(length(strata)),
                                "1" = "variable",
                                "variables")