# NOT RUN {
# get the pbc data from the survival package
require(survival)
data(pbc, package="survival")
# generate categorical versions of some of the baseline covariates
pbc$ageg[!is.na(pbc$age)] <-
ifelse(pbc$age[!is.na(pbc$age)] <= median(pbc$age, na.rm=TRUE), "Low", "High")
pbc$albuming[!is.na(pbc$albumin)]<-
ifelse(pbc$albumin[!is.na(pbc$albumin)] <= median(pbc$albumin, na.rm=TRUE), "Low", "High")
pbc$phosg[!is.na(pbc$alk.phos)] <-
ifelse(pbc$alk.phos[!is.na(pbc$alk.phos)]<= median(pbc$alk.phos,na.rm=TRUE), "Low", "High")
pbc$astg[!is.na(pbc$ast)] <-
ifelse(pbc$ast[!is.na(pbc$ast)] <= median(pbc$ast, na.rm=TRUE), "Low", "High")
pbc$bilig[!is.na(pbc$bili)] <-
ifelse(pbc$bili[!is.na(pbc$bili)] <= median(pbc$bili, na.rm=TRUE), "Low", "High")
pbc$cholg[!is.na(pbc$chol)] <-
ifelse(pbc$chol[!is.na(pbc$chol)] <= median(pbc$chol, na.rm=TRUE), "Low", "High")
pbc$copperg[!is.na(pbc$copper)] <-
ifelse(pbc$copper[!is.na(pbc$copper)] <= median(pbc$copper, na.rm=TRUE), "Low", "High")
pbc$ageg[is.na(pbc$age)] <- "No Data"
pbc$albuming[is.na(pbc$albumin)] <- "No Data"
pbc$phosg[is.na(pbc$alk.phos)] <- "No Data"
pbc$astg[is.na(pbc$ast)] <- "No Data"
pbc$bilig[is.na(pbc$bili)] <- "No Data"
pbc$cholg[is.na(pbc$chol)] <- "No Data"
pbc$copperg[is.na(pbc$copper)] <- "No Data"
#eliminate treatment NAs
pbcdat <- pbc[!is.na(pbc$trt), ]
# PFS and OS endpoints
set.seed(2006)
pbcdat$'event.pfs' <- sample(c(0,1),dim(pbcdat)[1],replace=TRUE)
pbcdat$'timepfs' <- sample(1:5000,dim(pbcdat)[1],replace=TRUE)
pbcdat$'event.os' <- pbcdat$event
pbcdat$'timeos' <- pbcdat$time
#variable importance for OS for the created categorical variables
#(higher is more important, also works for numeric variables)
varnames <- c('ageg', 'sex', 'bilig', 'cholg', 'astg', 'albuming', 'phosg')
# define function the eval_function()
# Attention: The eval_function ALWAYS needs to return a dataframe with one row.
# Include exception handling, like if(N1>0 && N2>0) hr <- exp(coxph(...) )
# to avoid program abort due to errors
hazardratio <- function(D) {
HRpfs <- tryCatch(exp(coxph(Surv(D$timepfs, D$event.pfs) ~ D$trt )$coefficients[[1]]),
warning=function(w) {NA})
HRpfs <- 1/HRpfs
HR.pfs <- round(HRpfs, 2)
HR.pfs[HR.pfs > 10] <- 10
HR.pfs[HR.pfs < 0.00001] <- 0.00001
HRos <- tryCatch(exp(coxph(Surv(D$timeos, D$event.os) ~ D$trt )$coefficients[[1]]),
warning=function(w) {NA})
HRos <- 1/HRos
HR.os <- round(HRos, 2)
HR.os[HR.os > 10] <- 10
HR.os[HR.os < 0.00001] <- 0.00001
data.frame( HR.pfs, HR.os#, N.of.subjects,N1 ,N2
)
}
# run subscreen
# }
# NOT RUN {
results <- subscreencalc(data=pbcdat,
eval_function=hazardratio,
endpoints = c("timepfs" , "event.pfs", "timeos", "event.os"),
treat="trt",
subjectid = "id",
factors=c("ageg", "sex", "bilig", "cholg", "copperg"),
use_complement = FALSE,
factorial = FALSE)
# visualize the results of the subgroup screening with a Shiny app
subscreenshow(results)
# }
Run the code above in your browser using DataLab