inst/doc/Expertsurv-Vignette.R

In expertsurv: Incorporate Expert Opinion with Parametric Survival Models

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = 'center',
  out.width = '70%',
  error = TRUE,
  message = TRUE,
  warning = TRUE
)

## -----------------------------------------------------------------------------
#A param_expert object; which is a list of 
#length equal to the number of timepoints
param_expert_example1 <- list()

#If we have 1 timepoint and 2 experts
#dist is the names of the distributions
#wi is the weight assigned to each expert (usually 1)
#param1, param2, param3 are the parameters of the distribution
#e.g. for norm, param1 = mean, param2 = sd
#param3 is only used for the t-distribution and is the degress of freedom.
#We allow the following distributions:
#c("normal","t","gamma","lognormal","beta") 


param_expert_example1[[1]] <- data.frame(dist = c("norm","t"),
                                         wi = c(0.5,0.5), # Ensure Weights sum to 1
                                         param1 = c(0.1,0.12),
                                         param2 = c(0.005,0.005),
                                         param3 = c(NA,3))
param_expert_example1

#Naturally we will specify the timepoint for which these probabilities where elicited

timepoint_expert <- 14


#In case we wanted a second timepoint -- Just for illustration

# param_expert_example1[[2]] <- data.frame(dist = c("norm","norm"),
#                                          wi = c(1,1),
#                                          param1 = c(0.05,0.045),
#                                          param2 = c(0.005,0.005),
#                                          param3 = c(NA,NA))
# 
# timepoint_expert <- c(timepoint_expert,18)


## ----echo = FALSE, fig.cap = "Expert prior distributions"---------------------
img_temp <- "Vignette_Example_1_Expert_Opinion.png"
knitr::include_graphics(system.file(paste0("image/",img_temp), package = "expertsurv"))

## ----echo = FALSE, fig.cap = "Model Comparison"-------------------------------

img_temp <- "Vignette_Example_1_DIC.png"
knitr::include_graphics(system.file(paste0("image/",img_temp), package = "expertsurv"))


## ----echo = FALSE, fig.cap = "Survival function with Expert prior"------------
img_temp <- "Vignette_Example_1.png"
knitr::include_graphics(system.file(paste0("image/",img_temp), package = "expertsurv"))

## ----echo = FALSE, fig.cap = "Survival function with Expert prior (left) and Vague prior (right)"----

img_temp <- "Vignette_Example_2.png"
knitr::include_graphics(system.file(paste0("image/",img_temp), package = "expertsurv"))


## ----echo = FALSE, fig.cap = "Survival difference"----------------------------

img_temp <- "Vignette_Example_3.png"
knitr::include_graphics(system.file(paste0("image/",img_temp), package = "expertsurv"))


## ----include = FALSE, eval=FALSE----------------------------------------------
#  library("expertsurv"); library("dplyr"); library("ggplot2")
#  data(bc)
#  set.seed(236236)
#  ## Age at diagnosis is correlated with survival time. A longer survival time
#  ## gives a younger mean age
#  bc$age <- rnorm(dim(bc)[1], mean = 65 - scale(bc$recyrs, scale=F), sd = 5)
#  ## Create age at diagnosis in days - used later for matching to expected rates
#  bc$agedays <- floor(bc$age * 365.25)
#  ## Create some random diagnosis dates between 01/01/1984 and 31/12/1989
#  bc$diag <- as.Date(floor(runif(dim(bc)[1], as.Date("01/01/1984", "%d/%m/%Y"),
#                                 as.Date("31/12/1989", "%d/%m/%Y"))),
#                     origin="1970-01-01")
#  ## Create sex (assume all are female)
#  bc$sex <- factor("female")
#  ## 2-level prognostic variable
#  bc$group2 <- ifelse(bc$group=="Good", "Good", "Medium/Poor")
#  head(bc)
#  
#  
#  ## reshape US lifetable to be a tidy data.frame, and convert rates to per person-year as flexsurv regression is in years
#  survexp.us.df <- as.data.frame.table(survexp.us, responseName = "exprate") %>%
#    mutate(exprate = 365.25 * exprate)
#  survexp.us.df$age <- as.numeric(as.character(survexp.us.df$age))
#  survexp.us.df$year <- as.numeric(as.character(survexp.us.df$year))
#  
#  ## Obtain attained age and attained calendar year in (whole) years
#  bc <- bc %>% mutate(attained.age.yr = floor(age + recyrs),
#                      attained.year = lubridate::year(diag + rectime))
#  
#  ## merge in (left join) expected rates at event time
#  bc <- bc %>% left_join(survexp.us.df, by = c("attained.age.yr"="age",
#                                               "attained.year"="year",
#                                               "sex"="sex"))
#  
#  param_expert_example_gpm <- list()
#  param_expert_example_gpm[[1]] <- data.frame(dist = c("norm"),
#                                           wi = c(1), # Ensure Weights sum to 1
#                                           param1 = c(0.40),
#                                           param2 = c(0.05),
#                                           param3 = c(NA))
#  
#  timepoint_expert <- 20
#  
#  example_gpm  <- expertsurv::fit.models.expert(formula=Surv(recyrs,censrec)~as.factor(group2),data=bc,
#                                 distr=c("gomp"),
#                                 method="mle",
#                                 opinion_type = "survival",
#                                 times_expert = timepoint_expert,
#                                 param_expert = param_expert_example_gpm,
#                                 id_St = min(which(bc$group2 =="Good")))
#  
#  expert_opinion_flex <- example_gpm$misc$flex_expert_opinion[[1]]
#  expert_opinion_flex$bhazard_par <- c(0.5)
#  
#  model.gomp.sep.rs <- expertsurv:::flexsurvreg(Surv(recyrs, censrec)~as.factor(group2),
#                                       data=bc, dist="gompertz",
#                                       anc = list(shape = ~ as.factor(group2)),
#                                       bhazard=exprate,expert_opinion = expert_opinion_flex)
#  
#  
#  
#  ## All-cause survival
#  ss.gomp.sep.rs.surv <- flexsurv::standsurv(model.gomp.sep.rs,
#                                      type = "survival",
#                                      at = list(list(group2 = "Good"),
#                                                list(group2 = "Medium/Poor")),
#                                      t = seq(0,30,length=50),
#                                      rmap=list(sex = sex,
#                                                year = diag,
#                                                age = agedays
#                                      ),
#                                      ratetable = survexp.us,
#                                      scale.ratetable = 365.25,
#                                      newdata = bc)
#  #> Marginal all-cause survival will be calculated
#  #> Calculating marginal expected survival and hazard
#  
#  # All-cause hazard
#  ss.gomp.sep.rs.haz <- flexsurv::standsurv(model.gomp.sep.rs,
#                                     type = "hazard",
#                                     at = list(list(group2 = "Good"),
#                                               list(group2 = "Medium/Poor")),
#                                     t = seq(0,30,length=50),
#                                     rmap=list(sex = sex,
#                                               year = diag,
#                                               age = agedays
#                                     ),
#                                     ratetable = survexp.us,
#                                     scale.ratetable = 365.25,
#                                     newdata = bc)
#  
#  plot(example_gpm, add.km = T, t = 0:50,plot_opinion  = TRUE)
#  ggsave("inst/image/Vignette_Example_4.png", units = "in", width = 8, height = 6)
#  
#  
#  plot(ss.gomp.sep.rs.surv, expected = T)+theme_bw()+ylab("Survival")+labs(color = "")+
#    scale_color_manual( values = c("Expected" = "Black", "group2=Good" = "red","group2=Medium/Poor" = "blue"),
#      labels = c("Expected" = "Expected", "group2=Good" ="Group - Good" , "group2=Medium/Poor"="Group - Medium/Poor" ))
#  ggsave("inst/image/Vignette_Example_5.png", units = "in", width = 8, height = 6)
#  
#  
#  plot(ss.gomp.sep.rs.haz, expected = T)+theme_bw()+ylab("Hazard")+labs(color = "")+
#    scale_color_manual( values = c("Expected" = "Black", "group2=Good" = "red","group2=Medium/Poor" = "blue"),
#      labels = c("Expected" = "Expected", "group2=Good" ="Group - Good" , "group2=Medium/Poor"="Group - Medium/Poor" ))
#  ggsave("inst/image/Vignette_Example_6.png", units = "in", width = 8, height = 6)
#  
#  

## ----echo = FALSE, fig.cap = "Survival Estimates without General Population Mortality"----

img_temp <- "Vignette_Example_4.png"
knitr::include_graphics(system.file(paste0("image/",img_temp), package = "expertsurv"))


## ----echo = FALSE, fig.cap = "All-cause Survival Estimates including General Population Mortality"----
img_temp <- "Vignette_Example_5.png"
knitr::include_graphics(system.file(paste0("image/",img_temp), package = "expertsurv"))

## ----echo = FALSE, fig.cap = "All-cause Hazard Functions including General Population Mortality"----
img_temp <- "Vignette_Example_6.png"
knitr::include_graphics(system.file(paste0("image/",img_temp), package = "expertsurv"))

## ----eval=FALSE, include=FALSE------------------------------------------------
#  
#  #Uniform prior on lambda implies non uniform St*
#  
#  t_expert <- 10
#  n <- 1000000
#  a <- 0
#  b <- .3
#  lambda<-runif(n,a,b)
#  St <-exp(-lambda*t_expert)
#  St_vec <- seq(0,1, by = 0.001)
#  
#  St_1 <- exp(-Inf*t_expert)
#  St_2<- exp(-b*t_expert)
#  k <- -log(St_2)/t_expert
#  
#  plot(density(St))
#  #lines(St_vec, -log(St_vec)/t_expert)
#  lines(St_vec, (1/(St_vec*t_expert))/k, col = "red")
#  
#  
#  
#  #Uniform St* implies non uniform lambda
#  
#  t_expert <- 10
#  n <- 1000000
#  a <- 0.01
#  b <- 0.999
#  St<-runif(n,a,b)
#  lambda <- -log(St)/t_expert
#  lambda_vec <- seq(0,1, by = 0.001)
#  
#  lambda_1 <- -log(b)/t_expert
#  lambda_2<- -log(a)/t_expert
#  
#  k <- -log(St_2)/t_expert
#  
#  plot(density(lambda))
#  #lines(St_vec, -log(St_vec)/t_expert)
#  lines(lambda_vec, t_expert*exp(-lambda_vec*t_expert), col = "red")
#  
#