tests/testthat/helper.R

In bpcp: Beta Product Confidence Procedure for Right Censored Data

# leukemia data for testing

# create data leuk with first event missing
leuk_path <- system.file("data","leuk.rda",package="bpcp")
load(leuk_path)
leuk_na <- leuk
leuk_na[1,"status"] <- NA

# make character version of treatment column in leuk2
leuk2_path <- system.file("data","leuk2.rda",package="bpcp")
load(leuk2_path)
leuk2$treatment_char <- as.character(leuk2$treatment)

# load sclerosis data
sclerosis_path <- system.file("data","sclerosis.rda",package="bpcp")
load(sclerosis_path)

# dataset with covariates for testing
df_covars <- data.frame(time = c(1,5,10,15,20), status = c(0,1,1,0,1), group=c("T", "T", "T", "P", "P"), sex=c("M","F","M","F","M"))

df_nocens <- data.frame(time = c(1,1,3,5,7,9), status = c(1,1,1,1,1,1), group = c("P", "P", "P", "T", "T","T"))


# T/F whether estimate is within a given interval
contains_est <- function(lower,upper,est){
  if(est <= upper & est >= lower){
    return(TRUE)
  } else{
    return(FALSE)
  }
}

createTime<-function(x,n,testtime=1){
  # creates time (with no censoring) such that 
  # S(1)=x/n, where S is Kaplan-Meier
  if (x==n){
    time<- testtime + 1:n
  } else if (x==0){
    time<- testtime* (c(1:n)/(n+1))
  } else {
    time<- c( testtime*c((x+1):n)/(n+1), 
              testtime+ 1:x)
  }
  time
}

create_nocensoring_data <- function(x1,n1,x2,n2){
  Time<-c( createTime(x1,n1),
           createTime(x2,n2))
  Status<- rep(1,n1+n2)
  Group<- c(rep(1,n1),rep(2,n2))
  
  out <- data.frame(Time=Time,
                    Status=Status,
                    Group=Group)
  
  out
  
}

pvalue.CI.compatability<-function(x1,n1,x2,n2,Parmtype,Alternative,Nullparm, changeGroupOrder = F, func="bpcp",zero.one=F,method="standard"){
  # make sure that if the p-value for a test is p, that the appropriate 
  # (1-p) confidence limit equals the null parameter
  #
  # first create the data set (uncensored so that at Testtime=1 we get 
  # S1(1)=x1/n2 and S2(1)=x2/n2, where Sj(t) is the Kaplan-Meier at t) 
  Time<-c( createTime(x1,n1),
           createTime(x2,n2))
  # createTime requires Testtime to be 1....
  Testtime<- 1 
  Status<- rep(1,n1+n2)
  Group<- c(rep(1,n1),rep(2,n2))
  
  if(func=="bpcp"){
    pval <- bpcp2samp(time=Time,
                   status=Status,
                   group=Group, 
                   testtime=Testtime,
                   parmtype=Parmtype,
                   alternative=Alternative, 
                   nullparm=Nullparm,
                   conf.level=0.95, 
                   midp=FALSE,
                   changeGroupOrder = changeGroupOrder)$p.value
  } else{
    pval <- delta2samp(time=Time,
                    status=Status,
                    group=Group, 
                    testtime=Testtime,
                    parmtype=Parmtype,
                    alternative=Alternative, 
                    nullparm=Nullparm,
                    conf.level=0.95,
                    changeGroupOrder = changeGroupOrder,
                    zero.one.adjustment = zero.one,
                    method=method)$p.value
  }
  

  if (Alternative=="less"){
    ci.index<-2
  } else if (Alternative=="greater"){
    ci.index<-1
  } else stop("for this function Alternative must be either 'less' or 'greater' ")
  
  if(func=="bpcp"){
    confLimit<-bpcp2samp(time=Time,
                 status=Status,
                 group=Group, 
                 testtime=Testtime,
                 parmtype=Parmtype,
                 alternative=Alternative, 
                 nullparm=Nullparm,
                 conf.level=1-pval, 
                 midp=FALSE,
                 changeGroupOrder =changeGroupOrder)$conf.int[ci.index]
  } else{
    confLimit <- delta2samp(time=Time,
                            status=Status,
                            group=Group, 
                            testtime=Testtime,
                            parmtype=Parmtype,
                            alternative=Alternative, 
                            nullparm=Nullparm,
                            conf.level=1-pval, 
                            changeGroupOrder =changeGroupOrder,
                            zero.one.adjustment = zero.one,
                            method=method)$conf.int[ci.index]
  }
  
  unname(confLimit)
}

bpcp2samp.nocensoring<-function(x1,n1,x2,n2,Parmtype,Alternative,Nullparm){
  Time<-c( createTime(x1,n1),
           createTime(x2,n2))
  Status<- rep(1,n1+n2)
  Group<- c(rep(1,n1),rep(2,n2))
  Testtime<- 1 
  
  
  out<-bpcp2samp(time=Time,
                 status=Status,
                 group=Group, 
                 testtime=Testtime,
                 parmtype=Parmtype,
                 alternative=Alternative, 
                 nullparm=Nullparm,
                 conf.level=0.95, 
                 midp=FALSE)
  out
}


pvalue.CI.compatability.fixtdelta<-function(x1,n1,x2,n2,Parmtype,Alternative,Nullparm, Adjustment){
  # make sure that if the p-value for a test is p, that the appropriate 
  # (1-p) confidence limit equals the null parameter
  #
  # first create the data set (uncensored so that at Testtime=1 we get 
  # S1(1)=x1/n2 and S2(1)=x2/n2, where Sj(t) is the Kaplan-Meier at t) 
  Time<-c( createTime(x1,n1),
           createTime(x2,n2))
  # createTime requires Testtime to be 1....
  Testtime<- 1 
  Status<- rep(1,n1+n2)
  Group<- c(rep(1,n1),rep(2,n2))
  
  # 
  pval<-fixtdelta(time=Time,
                  status=Status,
                  group=Group, 
                  testtime=Testtime,
                  parmtype=Parmtype,
                  alternative=Alternative, 
                  nullparm=Nullparm,
                  conf.level=0.95, 
                  adjustment=Adjustment)$p.value
  if (Alternative=="less"){
    ci.index<-2
  } else if (Alternative=="greater"){
    ci.index<-1
  } else stop("for this function Alternative must be either 'less' or 'greater' ")
  confLimit<-fixtdelta(time=Time,
                       status=Status,
                       group=Group, 
                       testtime=Testtime,
                       parmtype=Parmtype,
                       alternative=Alternative, 
                       nullparm=Nullparm,
                       conf.level=1-pval, 
                       adjustment=Adjustment)$conf.int[ci.index]
  confLimit
}


fixtdelta.nocensoring<-function(x1,n1,x2,n2,Parmtype,Alternative,Nullparm,CL=0.95,adjustment){
  Time<-c( createTime(x1,n1),
           createTime(x2,n2))
  Status<- rep(1,n1+n2)
  Group<- c(rep(1,n1),rep(2,n2))
  Testtime<- 1 
  
  
  out<-fixtdelta(time=Time,
                 status=Status,
                 group=Group, 
                 testtime=Testtime,
                 parmtype=Parmtype,
                 alternative=Alternative, 
                 nullparm=Nullparm,
                 conf.level=CL,
                 adjustment=adjustment)
  out
}

### One sample helper functions ###

# no events
no_events <- function(time,type=c("bpcp","km"),param="ci", arms=1){
  
  if(type=="bpcp"){
    
    res <- summary(bpcpfit(time = c(5,5,5,5,5), status = c(0,0,0,0,0)))
    
  } 
  
  if(param=="ci"){
    return(c(res[time,3],res[time,4]))
  } 
}

# no censoring

no_censoring <- function(time,type=c("bpcp","km"),param=c("ci","p"),arms=1,parmtype="difference",midp=F,conf.level=0.95, nullparm=NULL,  alternative = c("two.sided", "less", "greater"),changeGroupOrder=FALSE){
  
  if(type=="bpcp"){
    
    if(arms==1){
      
      res <- summary(bpcpfit(time = df_nocens$time, status = df_nocens$status, midp=midp,conf.level=conf.level))
      
      if(param=="ci"){
        
        out <- c(res[time+1,3],res[time+1,4])
      } 
      
    } else{
      
      testtime=unique(df_nocens$time)[time]
      
      res <- bpcp2samp(time = df_nocens$time, status = df_nocens$status, df_nocens$group, testtime = testtime,parmtype=parmtype,midp=midp,conf.level=conf.level, nullparm=nullparm, alternative=alternative,changeGroupOrder=changeGroupOrder)
      
      if(param=="ci"){
        out <- unname(c(res$conf.int[1],res$conf.int[2]))
        attr(out,"conf.level") <- 0.95
        
        
      } else{ out <- unname(res$p.value)}
      
      
    }
    }

    
  
  
  return(out)
}

# all censored

all_censored <- function(time,type=c("bpcp","km"),param="ci",arms=1){
  
  if(type=="bpcp"){
    
    res <- summary(bpcpfit(time = c(1,5,10,15,20), status = c(0,0,0,0,0)))
    
  } 
  
  if(param=="ci"){
    return(c(res[time,3],res[time,4]))
  } 
}

# last censored

last_censored  <- function(time,param=c("ci","est"),arms=1){
    
    res <- summary(bpcpfit(time = c(1,5,10,15,20), status = c(1,1,1,1,0)))
    
    if(param=="ci"){
      return(c(res[time,3],res[time,4]))
    }
    
    if(param=="est"){
      return(c(res[time,2]))
    }
  
   
}


# data.frame of Borkowf's Table 7
borkowf_tb7 <- data.frame(
  k = c(0, 1, 1, 2, 2, 3, 3, 3, 4, 5, 5, 5, 5, 6, 7, 7, 7, 7, 7),
  Time = c("0", "6", "6+", "7", "9+", "10", "10+", "11+", "13", "16", "17+", "19+", "20+", "22", "23", "25+", "32+", "34+", "35"),
  Interval = c(
    "[0,6)", "[6,6]", "(6,7)", "[7,9]", "(9,10)","[10,10]", "(10,11]", "(11,13)", 
    "[13,16)", "[16,17]", "(17,19]", "(19,20]", "(20,22)", "[22,23)", "[23,25]", 
    "(25,32]", "(32,34]", "(34,35]", "(35,∞)"
  ),
  
  # Method (a) Lower, Upper
  a_Lower = c(1.000, 0.707, 0.707, 0.636, 0.636, 0.564, 0.564, 0.564, 0.481, 0.404, 0.404, 0.404, 0.404, 0.286, 0.184, 0.184, 0.184, 0.184, 0.184),
  a_Upper = c(1.000, 1.000, 1.000, 0.977, 0.977, 0.942, 0.942, 0.942, 0.900, 0.851, 0.851, 0.851, 0.851, 0.789, 0.712, 0.712, 0.712, 0.712, 0.712),
  
  # Method (c) Lower, Upper
  c_Lower = c(1.000, 0.707, 0.704, 0.634, 0.629, 0.559, 0.554, 0.548, 0.470, 0.398, 0.391, 0.383, 0.374, 0.277, 0.186, 0.176, 0.153, 0.138, 0.122),
  c_Upper = c(1.000, 1.000, 1.000, 0.980, 0.984, 0.947, 0.952, 0.958, 0.910, 0.857, 0.864, 0.872, 0.881, 0.799, 0.710, 0.720, 0.744, 0.758, 0.775),
  
  # Method (d) Lower, Upper
  d_Lower = c(0.935, 0.700, 0.697, 0.629, 0.624, 0.556, 0.551, 0.545, 0.469, 0.397, 0.390, 0.382, 0.373, 0.277, 0.186, 0.176, 0.153, 0.138, 0.122),
  d_Upper = c(1.000, 1.000, 1.000, 0.985, 0.989, 0.950, 0.955, 0.961, 0.912, 0.858, 0.865, 0.873, 0.882, 0.799, 0.710, 0.720, 0.744, 0.758, 0.775),
  
  # Method (e) Lower, Upper
  e_Lower = c(0.911, 0.683, 0.680, 0.614, 0.610, 0.544, 0.538, 0.533, 0.460, 0.391, 0.384, 0.376, 0.367, 0.275, 0.189, 0.179, 0.155, 0.141, 0.124),
  e_Upper = c(1.000, 0.997, 1.000, 0.970, 0.975, 0.938, 0.943, 0.949, 0.903, 0.852, 0.859, 0.867, 0.875, 0.797, 0.713, 0.722, 0.746, 0.761, 0.777)
)

# Cancer data used in simulation
data(cancer,package="survival")

# example with 100% failure rate in observed
d_failure <- subset(colon,etype==1 & extent==1 & rx !="Lev")

# change time to years
d_failure$time<- d_failure$time/365.25

# refactor rx
d_failure$rx <- factor(d_failure$rx, levels=c("Obs","Lev+5FU"))

# data from Borkowf's paper
d<-leuk2[leuk2$treatment=="6-MP",]
d<-d[order(d$time),]


# make dataset with 100% survival in 1 arm, ending on 0% in the other
d_failure_2 <- d_failure
d_failure_2[which(d_failure_2$rx=="Lev+5FU"),]$status <- 1

# make dataset with 0% survival in 1 arm, variable survival in the other
d_failure_3 <- d_failure
d_failure_3[which(d_failure_3$rx=="Obs"),]$status <- 1