In OakleyJ/SHELF: Tools to Support the Sheffield Elicitation Framework

knitr::opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE,
                      fig.pos = 'h',
                      fig.align = 'center')

knitr::opts_chunk$set(fig.height = 3,
                      fig.width = 4)

library(ggplot2)

makeSurvSummaryTable <- function(survDF, sf = 3, useWeights = FALSE){

  nTreatments <- length(levels(survDF$treatment))

  sTable <- array(NA, c(nTreatments, 5))
  rownames(sTable) <- levels(survDF$treatment)
  colnames(sTable) <- c("n", "events", "minimum", "median", "maximum")

  if(useWeights == FALSE){
  sv <- survival::survfit(survival::Surv(time, event) ~ treatment,
                          data = survDF)}else{
                            sv <- survival::survfit(survival::Surv(time, event) ~ treatment,
                                                    weights = weights,
                                                    data = survDF)

                          }
  table_output <- summary(sv)$table

  # Need to keep array format with column names if only one treatment group:
  if(nTreatments == 1){
    dnames <- names(table_output)
    dim(table_output) <- c(1, 9)
    colnames(table_output) <- dnames
  }

  sTable[, c("n", "events", "median")] <- table_output[, c("records", "events", "median")]
  sTable[, c("minimum")] <- tapply(survDF$time, survDF$treatment, min)
  sTable[, c("maximum")] <- tapply(survDF$time, survDF$treatment, max)

  sTable[, c("minimum", "maximum", "median") ] <- 
    signif( sTable[, c("minimum", "maximum", "median")  ], sf)

  sTable

}

Survival data

nTreatments <- length(levels(params$survivalDF$treatment))
sdf <- makeSurvSummaryTable(params$survivalDF, useWeights = params$useWeights)
knitr::kable(sdf, caption = "Summary statistics for the provided survival data")

 if(params$useWeights == TRUE){
        fit <- survival::survfit(survival::Surv(time, event) ~ treatment,
                                 weights = weights,
                                 data = params$survivalDF)}else{
        fit <- survival::survfit(survival::Surv(time, event) ~ treatment, data = params$survivalDF)
                                 }


      myplot<- survminer::ggsurvplot(fit, data = params$survivalDF, censor = TRUE,
                                      legend = "right",
                                      legend.title = "",
                            conf.int = TRUE,
                            legend.labs = levels(params$survivalDF$treatment),
                            xlim = c(0, params$targetTime),
                            xlab = paste0("Time (", params$timeUnit, ")"),
                            break.time.by = params$targetTime/8)
      myplot$plot +
        geom_vline(xintercept = params$targetTime, linetype="dotted")

knitr::kable(makeSurvivalTable(params$survivalDF,
                               params$breakTime,
                               params$truncationTime,
                               params$timeUnit,
                               dp = 2,
                               useWeights = params$useWeights),
             caption = "The survivor column is the proportion surviving after the end of the corresponding time interval. The hazard column is the proportion who do not survive to the end of the time interval, out of those who survived to the beginning of the time interval.")

r if(params$reportGroup1){paste0('# Results for treatment group: "', levels(params$survivalDF$treatment)[1], '"')}

r if(params$reportGroup1){"## Individual judgements"}

knitr::kable(params$individual[[1]])

if(params$inputMethod == "quartiles"){
  p1 <- plotQuartiles(vals = params$individual[[1]][2:4, ],
                              lower = params$individual[[1]][1, ],
                              upper = params$individual[[1]][5, ],
                              expertnames = colnames(params$individual[[1]]),
                              xlabel = "Surivivor proportion")
}

if(params$inputMethod == "tertiles"){
  p1 <- plotTertiles(vals = params$individual[[1]][2:4, ],
                              lower = params$individual[[1]][1, ],
                              upper = params$individual[[1]][5, ],
                              expertnames = colnames(params$individual[[1]]),
                              xlabel = "Surivivor proportion")
}
print(p1 + theme_bw()+ coord_flip())

r if(params$reportGroup1 & params$reportScenarioTest){"## Scenario Testing"}

sce <- survivalScenario(tLower = 0,
                        tUpper = min(params$truncationTime,
                                     max(params$survivalDF$time)),
                        expLower= params$expRange[1],
                        expUpper = params$expRange[2],
                        expGroup = levels(params$survivalDF$treatment)[1],
                        tTarget = params$targetTime,
                        survDf = params$survivalDF,
                        groups = levels(params$survivalDF$treatment),
                        xl = paste0("Time t (", params$timeUnit, ")"),
                        showPlot = FALSE,
                        useWeights = params$useWeights)
print(sce$KMplot)

r if(params$reportGroup1){"## RIO distribution"}

```{asis, echo = params$reportGroup1} The elicited RIO probabilities were as follows

```r
captionTextGroup1 <- paste0("The probability density function and cumulative distribution function fitted to RIO's probabilities for the Quantity of Interest: 
                            the proportion suriving for at least time ", params$targetTime, " ",
                            params$timeUnit, ' in the treatment group "',
                            levels(params$survivalDF$treatment)[1], '".')
expert <- 1
sf <- 3
mydf <- data.frame( params$allfits[[1]]$vals[expert, ], params$allfits[[1]]$probs[expert, ])
colnames(mydf) <- c( "$x$", "$P(X\\le x)$")
knitr::kable(mydf)

plotfit(params$allfits[[1]], xlab = "survival proportion", ylab = "density", d = params$dist1)
makeCDFPlot(lower = 0, upper = 1, v = params$allfits[[1]]$vals[1, ],
            p = params$allfits[[1]]$probs[1, ], dist = params$dist1,
            showFittedCDF = TRUE,  fit = params$allfits[[1]],
            xlab = "survival proportion", ylab = "cumulative probability")

tableCaption1 <-  paste0("Percentiles from the distribution fitted to RIO's probabilities for the Quantity of Interest: the proportion suriving for at least time ", params$targetTime, " ",
                            params$timeUnit, ' in the treatment group "',
                            levels(params$survivalDF$treatment)[1], '".')
fb1 <- feedback(params$allfits[[1]],
                quantiles = c(0.01, 0.1, 0.5, 0.9, 0.99))$fitted.quantiles[, params$dist1]
fb1 <- matrix(fb1, nrow = 1)
colnames(fb1) <- paste0(c(1, 10, 50, 90, 99), "%")
knitr::kable(fb1, caption = tableCaption1)

r if(params$reportGroup2){paste0('# Results for treatment group: "', levels(params$survivalDF$treatment)[2], '"')}

r if(params$reportGroup2){"## Individual judgements"}

knitr::kable(params$individual[[2]])

if(params$inputMethod == "quartiles"){
  p2 <- plotQuartiles(vals = params$individual[[2]][2:4, ],
                              lower = params$individual[[2]][1, ],
                              upper = params$individual[[2]][5, ],
                              expertnames = colnames(params$individual[[2]]),
                              xlabel = "Surivivor proportion")
}

if(params$inputMethod == "tertiles"){
  p2 <- plotTertiles(vals = params$individual[[2]][2:4, ],
                              lower = params$individual[[2]][1, ],
                              upper = params$individual[[2]][5, ],
                              expertnames = colnames(params$individual[[2]]),
                              xlabel = "Surivivor proportion")
}
print(p2 + theme_bw()+ coord_flip())

r if(params$reportGroup2 & params$reportScenarioTest){"## Scenario Testing"}

sce <- survivalScenario(tLower = 0,
                        tUpper = min(params$truncationTime,
                                     max(params$survivalDF$time)),
                        expLower= params$expRange[1],
                        expUpper = params$expRange[2],
                        expGroup = levels(params$survivalDF$treatment)[2],
                        tTarget = params$targetTime,
                        survDf = params$survivalDF,
                        groups = levels(params$survivalDF$treatment),
                        xl = paste0("Time t (", params$timeUnit, ")"),
                        showPlot = FALSE,
                        useWeights = params$useWeights)
print(sce$KMplot)

r if(params$reportGroup2){"## RIO distribution"}

```{asis, echo = params$reportGroup2} The elicited RIO probabilities were as follows

```r
captionTextGroup2 <- paste0("The probability density function and cumulative distribution function fitted to RIO's probabilities for the Quantity of Interest: 
                            the proportion suriving for at least time ", params$targetTime, " ",
                            params$timeUnit, ' in the treatment group "',
                            levels(params$survivalDF$treatment)[2], '".')
expert <- 1
sf <- 3
mydf <- data.frame( params$allfits[[2]]$vals[expert, ], params$allfits[[1]]$probs[expert, ])
colnames(mydf) <- c( "$x$", "$P(X\\le x)$")
knitr::kable(mydf)

plotfit(params$allfits[[2]], xlab = "survival proportion", ylab = "density", d = params$dist2)
makeCDFPlot(lower = 0, upper = 1, v = params$allfits[[2]]$vals[1, ],
            p = params$allfits[[2]]$probs[1, ], dist = params$dist2,
            showFittedCDF = TRUE, fit = params$allfits[[2]],
            xlab = "survival proportion", ylab = "cumulative probability")

tableCaption2 <-  paste0("Percentiles from the distribution fitted to RIO's probabilities for the Quantity of Interest: the proportion suriving for at least time ", params$targetTime, " ",
                            params$timeUnit, ' in the treatment group "',
                            levels(params$survivalDF$treatment)[2], '".')
fb2 <- feedback(params$allfits[[2]],
                quantiles = c(0.01, 0.1, 0.5, 0.9, 0.99))$fitted.quantiles[, params$dist2]
fb2 <- matrix(fb2, nrow = 1)
colnames(fb2) <- paste0(c(1, 10, 50, 90, 99), "%")
knitr::kable(fb2, caption = tableCaption2)

r if(params$reportExtrapolation){"# Extrapolation results"}

if(length(params$allfits) > 1){
  group2RIO <- params$allfits[[2]]
  fqRIO2 <- params$dist2}else{
    group2RIO <- NULL
    fqRIO2 <- NULL
}

survivalExtrapolatePlot(params$survivalDF,
                        myfit1 = params$allfits[[1]],
                                          myfit2 = group2RIO,
                                          fqDist1 = params$dist1,
                                          fqDist2 = fqRIO2,
                                          tTruncate = params$truncationTime,
                                          tTarget = params$targetTime,
                                          alpha = params$alpha,
                                          useWeights = params$useWeights,
                                          groups = levels(params$survivalDF$treatment),
                                          xl = paste0("Time t (", params$timeUnit, ")"),
                                          breakTime = params$targetTime/8,
                                          showPlot = TRUE,
                                          returnPlot = FALSE)

r if(params$reportDistributions){"# Appendix: all fitted distributions"} r if(params$reportDistributions){paste0('## Treatment group: "', levels(params$survivalDF$treatment)[1], '"')}

fit <- params$allfits[[1]]
bin.left <- NA
bin.right <- NA
chips <- NA
roulette <- FALSE
filename <- system.file("shinyAppFiles", "distributionsChild.Rmd", package="SHELF")

r if(params$reportDistributions & params$reportGroup2 ){paste0('## Treatment group: "', levels(params$survivalDF$treatment)[2], '"')}

fit <- params$allfits[[2]]
bin.left <- NA
bin.right <- NA
chips <- NA
roulette <- FALSE
filename <- system.file("shinyAppFiles", "distributionsChild.Rmd", package="SHELF")