R/elicitSurvivalExtrapolation.R
In OakleyJ/SHELF: Tools to Support the Sheffield Elicitation Framework
Defines functions
checkPlot
makeSurvSummaryTable
elicitSurvivalExtrapolation
Documented in
elicitSurvivalExtrapolation
#' Elicitation for survival extrapolation
#'
#' Opens up a web browser in which you can implement the SHELF protocol for survival extrapolation.
#' Start with uploading a .csv file of individual patient survival data (time, event to indicate censoring, and
#' treatment group). Then elicit individual judgements, perform scenario testing as required, and elicit a RIO distribution.
#' Judgements for two treatment groups can be elicited in the same session.
#'
#'
#' @author Jeremy Oakley <j.oakley@@sheffield.ac.uk>
#' @examples
#'
#' \dontrun{
#'
#' # make a suitable csv file using a built in data set from the survival package
#' sdf <- survival::veteran[, c("time", "status", "trt")]
#' colnames(sdf) <- c("time", "event", "treatment")
#' sdf$treatment <- factor(sdf$treatment, labels = c("standard", "test"))
#'
#' # write the data frame sdf to a .csv file in the current working directory
#' write.csv(sdf, file = "testFile.csv", row.names = FALSE)
#'
#' # Run the app and upload testFile.csv in the first tab, and change unit of time to "days"
#'
#' elicitSurvivalExtrapolation()
#'
#' }
#' @import shiny
#' @import survival
#' @importFrom survminer ggsurvplot
#' @export
elicitSurvivalExtrapolation<- function(){
runApp(list(
ui = shinyUI(fluidPage(
# Application title
titlePanel("SHELF: survival extrapolation"),
# sidebarLayout(
mainPanel(tags$style(type="text/css",
".shiny-output-error { visibility: hidden; }",
".shiny-output-error:before { visibility: hidden; }"
),
tabsetPanel(
tabPanel("Upload data",
sidebarLayout(
sidebarPanel(
wellPanel(
fileInput("file", label = "Survival data upload"),
uiOutput("targetTime"),
uiOutput("truncationTime"),
uiOutput("timeUnit"),
uiOutput("caseWeights"),
uiOutput("fontSize")
)
)
,
mainPanel(
helpText('Upload a .csv file with three columns:
"time", "event" and "treatment" (in that order).
Values in the "event" column should be 0 for
a censored observation, and 1 otherwise. The
"treatment" column should be included even
if there is only one treatment group. If the
data include case weights for each observation,
include a fourth column "weights"'),
plotOutput("KMbasic"),
tableOutput("survivalSummary"),
)
)
),
tabPanel("KM table",
sidebarLayout(
sidebarPanel(
wellPanel(
uiOutput("breakTime"),
)),
mainPanel(
helpText('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'),
tableOutput("KMdata")
)
)
),
tabPanel("Individual judgements",
sidebarLayout(
sidebarPanel(
wellPanel(
uiOutput("treatmentGroup"),
numericInput("nExperts", label = h5("Number of experts"),
value = 2, min = 1),
radioButtons("elicMethod", "Elicitation method",
c("Quartiles" = "quartiles",
"Tertiles" = "tertiles")),
textInput("indAxis", label = h5("Axis limits"),
value = "0, 1"),
hr(style = "border-top: 1px solid #000000;"),
h5("Scenario testing"),
helpText("Use next tab before displaying here"),
checkboxInput("showScenario",
label = h5("Show scenario interval"))
)),
mainPanel(
h4("Quantity of interest:"),
textOutput("defn"),
hr(),
helpText('Enter plausible limits and quantiles as values
between 0 and 1.'),
uiOutput("EnterQuartilesGp1"),
uiOutput("EnterQuartilesGp2"),
plotOutput("individualQuartiles")
)
)
),
tabPanel("Scenario testing",
sidebarLayout(
sidebarPanel(
wellPanel(
uiOutput("scenarioGroup"),
uiOutput("scenarioTruncationTime"),
uiOutput("expRange")
)),
mainPanel(
helpText('In this tab you can consider a scenario in which the hazard
remains constant after a specified time point. The shaded
bar indicates a point-wise approximate 95% credible interval for
the survivor function based on this assumption. Intervals can be
added to the "Individual judgements" tab.'),
plotOutput("KMplot"),
htmlOutput("reportInterval")
)
)
),
tabPanel("RIO elicitation",
sidebarLayout(
sidebarPanel(
wellPanel(
uiOutput("RIOtreatmentGroup"),
uiOutput("RIOvaluesGp1"),
uiOutput("RIOprobsGp1"),
uiOutput("RIOvaluesGp2"),
uiOutput("RIOprobsGp2"),
uiOutput("RIOdist1"),
uiOutput("RIOdist2"),
hr(style = "border-top: 1px solid #000000;"),
textInput("RIOfq", label = h5("Feedback quantiles"),
value = "0.1, 0.9"),
textOutput("RIOfq"),
h5("Feedback probability"),
fluidRow(
column(width = 7,
"Probability QoI <= "
),
column(width = 5,
numericInput("RIOfeedbackProbability", label = NULL,
value = 0.5, min = 0, max = 1)
)
),
textOutput("RIOfp"),
hr(style = "border-top: 1px solid #000000;"),
uiOutput("xaxis1"),
uiOutput("xaxis2")
)),
mainPanel(
h4("Quantity of Interest (QoI):"),
textOutput("RIOdefn"),
hr(),
plotOutput("pdfPlot"),
plotOutput("cdfPlot")
)
)
),
tabPanel("Extrapolation plot",
sidebarLayout(
sidebarPanel(
wellPanel(
numericInput("alpha", "Credible Interval Size",
0.95, min = 0, max = 1),
)),
mainPanel(
plotOutput("extrapolationPlot")
)
)
)
,
tabPanel("Report",
wellPanel(
h5("Report content"),
checkboxInput("reportData", "KM plot and summary data",
value = TRUE, width = NULL),
checkboxInput("reportTable", "Survivor and hazard table",
value = TRUE, width = NULL),
uiOutput("reportGroup1"),
uiOutput("reportGroup2"),
checkboxInput("reportIndividual", "Individual judgements",
value = TRUE, width = NULL),
checkboxInput("reportScenarioTest", "Scenario Test results",
value = TRUE, width = NULL),
checkboxInput("reportExtrapolation", "Plot extrapolation",
value = TRUE, width = NULL),
checkboxInput("reportDistributions", "All fitted distributions",
value = TRUE, width = NULL),
selectInput("outFormat", label = "Report format",
choices = list('html' = "html_document",
'pdf' = "pdf_document",
'Word' = "word_document"),
width = "30%"),
downloadButton("report", "Download report")
)
)
)
)
)
),
server = function(input, output) {
# File uploading ----
survivalDF <- reactive({
req(input$file)
inFile <- input$file
sdf <- utils::read.csv(inFile$datapath)
colnames(sdf)[1:3] <- c("time", "event", "treatment")
sdf$treatment <- as.factor(sdf$treatment)
if(length(levels(sdf$treatment)) > 2){
showNotification(
"This app will only with work with one or two treatment groups.",
duration = 60,
type = "error"
)
sdf <- NULL
}
sdf
})
output$KMbasic <- renderPlot({
req(survivalDF(),input$targetTime, input$truncationTime)
# Truncate data frame for plotting
index <- survivalDF()$time > input$truncationTime
truncatedDf <- survivalDF()
truncatedDf$event[index] <- 0
truncatedDf$time[index] <- input$truncationTime
if(caseWeight_value() == TRUE){
fit <- survival::survfit(survival::Surv(time, event) ~ treatment,
weights = weights,
data = truncatedDf)}else{
fit <- survival::survfit(survival::Surv(time, event) ~ treatment, data = truncatedDf)
}
myplot<- survminer::ggsurvplot(fit, data = truncatedDf, censor = TRUE,
legend = "right",
legend.title = "",
conf.int = TRUE,
legend.labs = levels(truncatedDf$treatment),
xlim = c(0, input$targetTime),
xlab = paste0("Time t (", input$timeUnit, ")"),
ylab = "S(t)",
break.time.by = input$targetTime/8)
myplot$plot +
geom_vline(xintercept = input$targetTime, linetype="dotted") +
theme_bw(base_size = input$fs)
})
output$survivalSummary <-renderTable({
req(survivalDF())
makeSurvSummaryTable(survivalDF(),
useWeights = caseWeight_value())},
rownames = TRUE)
output$targetTime <- renderUI({
req(survivalDF())
maxTime <- signif(max(survivalDF()$time), 1)
numericInput("targetTime", "Target extrapolation time", value = 2* maxTime)
})
output$truncationTime <- renderUI({
req(survivalDF())
maxTime <- signif(max(survivalDF()$time), 1)
numericInput("truncationTime", "Data truncation time", value = 2* maxTime)
})
output$timeUnit <- renderUI({
req(survivalDF())
selectInput("timeUnit", "Unit of time",
c("days", "weeks", "months", "years"),
selected = "years")
})
output$caseWeights <- renderUI({
req(survivalDF())
if("weights" %in% colnames(survivalDF())){
checkboxInput("useWeights", "Use case weights", value = TRUE)
}
})
caseWeight_value <- reactive({
if (is.null(input$useWeights)) {
return(FALSE)
} else {
return(input$useWeights)
}
})
output$fontSize <- renderUI({
req(survivalDF())
numericInput("fs", label = "Font size (all plots)", value = 16)
})
# KM table tab ----
output$breakTime <- renderUI({
req(survivalDF(), input$timeUnit)
maxTime <- min(max(survivalDF()$time), input$truncationTime)
numericInput("breakTime", paste0("Time interval duration (",
input$timeUnit, ")"),
value = signif(maxTime/4, 1))
})
output$KMdata <- renderTable({
req(survivalDF(), input$breakTime,
input$truncationTime)
makeSurvivalTable(survivalDF(), input$breakTime,
input$truncationTime, input$timeUnit, dp = 2,
useWeights = caseWeight_value())
})
# Scenario testing tab ----
output$scenarioGroup <- renderUI({
selectInput("scenarioGroup", "Treatment group",
choices = levels(survivalDF()$treatment))
})
output$expRange <- renderUI({
req(survivalDF())
maxTime <- min(max(survivalDF()$time), input$truncationTime)
textInput("expRange", label = h5("Constant hazard fitting period"),
value = paste0(0, ", ", signif(maxTime, 1) ))
})
expRange <- reactive({
tryCatch(eval(parse(text = paste("c(", input$expRange, ")"))),
error = function(e){NULL})
})
scenario <- reactiveValues()
observe({
req(survivalDF(), expRange(), input$scenarioGroup,
input$targetTime,
input$truncationTime)
sce <- survivalScenario(tLower = 0,
tUpper = min(input$truncationTime,
max(survivalDF()$time)),
expLower= expRange()[1],
expUpper = expRange()[2],
expGroup = input$scenarioGroup,
tTarget = input$targetTime,
survDf = survivalDF(),
groups = levels(survivalDF()$treatment),
xl = paste0("Time t (", input$timeUnit, ")"),
showPlot = FALSE,
fontsize = input$fs,
useWeights = caseWeight_value())
scenario$plot <- sce$KMplot
scenario$interval <- sce$interval
})
output$KMplot <- renderPlot({
# req(survivalDF(), expRange(), input$scenarioGroup,
# input$targetTime,
# input$truncationTime)
# scenario <- survivalScenario(tLower = 0,
# tUpper = min(input$truncationTime,
# max(survivalDF()$time)),
# expLower= expRange()[1],
# expUpper = expRange()[2],
# expGroup = input$scenarioGroup,
# tTarget = input$targetTime,
# survDf = survivalDF(),
# groups = levels(survivalDF()$treatment),
# xl = paste0("Time (", input$timeUnit, ")"),
# showPlot = FALSE,
# fontsize = input$fs)
# scenario$KMplot
req(scenario$plot)
scenario$plot
})
scenarioIntervalGp1 <- reactive({
req(survivalDF(), expRange(),input$scenarioGroup,
input$targetTime,
input$truncationTime)
sce <- survivalScenario(tLower = 0,
tUpper = min(input$truncationTime,
max(survivalDF()$time)),
expLower= expRange()[1],
expUpper = expRange()[2],
expGroup = levels(survivalDF()$treatment)[1],
tTarget = input$targetTime,
survDf = survivalDF(),
groups = levels(survivalDF()$treatment),
xl = paste0("Time (", input$timeUnit, ")"),
showPlot = FALSE,
fontsize = input$fs,
useWeights = caseWeight_value())
sce$interval
})
scenarioIntervalGp2 <- reactive({
req(survivalDF(), expRange(),input$scenarioGroup,
input$targetTime,
input$truncationTime)
interval <- NULL
if(length(levels(survivalDF()$treatment)) > 1){
interval <- survivalScenario(tLower = 0,
tUpper = min(input$truncationTime,
max(survivalDF()$time)),
expLower= expRange()[1],
expUpper = expRange()[2],
expGroup = levels(survivalDF()$treatment)[2],
tTarget = input$targetTime,
survDf = survivalDF(),
groups = levels(survivalDF()$treatment),
xl = paste0("Time (", input$timeUnit, ")"),
showPlot = FALSE,
fontsize = input$fs,
useWeights = caseWeight_value())$interval
}
interval
})
output$reportInterval <- renderUI({
# req(input$scenarioGroup, survivalDF())
# txt <- NULL
# if(input$scenarioGroup == levels(survivalDF()$treatment)[1]){
# req(scenarioIntervalGp1())
# txt <- paste0("(", signif(scenarioIntervalGp1()[1], 2)
# , ", ", signif(scenarioIntervalGp1()[2], 2), ")" )
# }
#
# if(input$scenarioGroup != levels(survivalDF()$treatment)[1]){
# req(scenarioIntervalGp2())
# txt <- paste0("(", signif(scenarioIntervalGp2()[1], 2), ", ",
# signif(scenarioIntervalGp2()[2], 2), ")" )
# }
req(scenario$interval)
txt <- paste0("(", signif(scenario$interval[1], 2),
", ", signif(scenario$interval[2], 2), ")" )
HTML(paste(hr(),
h5(paste0("Approximate 95% credible interval at the target time: ",
txt))))
})
# Individual judgements tab ----
initialVals1 <- reactive({
initialdf <- matrix(c(0, 0.25, 0.5, 0.75, 1), nrow = 5, ncol = input$nExperts
)
if(input$elicMethod == "quartiles"){
rownames(initialdf) <- c("L", "Q1", "M", "Q3", "U")
}
if(input$elicMethod == "tertiles"){
rownames(initialdf) <- c("L", "T1", "M", "T2", "U")
}
colnames(initialdf) <- LETTERS[1:input$nExperts]
initialdf
})
initialVals2 <- reactive({
initialdf <- matrix(0.1*(1:5), nrow = 5, ncol = input$nExperts
)
if(input$elicMethod == "quartiles"){
rownames(initialdf) <- c("L", "Q1", "M", "Q3", "U")
}
if(input$elicMethod == "tertiles"){
rownames(initialdf) <- c("L", "T1", "M", "T2", "U")
}
colnames(initialdf) <- LETTERS[1:input$nExperts]
initialdf
})
output$treatmentGroup <- renderUI({
selectInput("treatmentGroup", "Treatment group",
choices = levels(survivalDF()$treatment))
})
output$tgpNumber <- reactive({
which(input$treatmentGroup == levels(survivalDF()$treatment))
})
outputOptions(output, 'tgpNumber', suspendWhenHidden = FALSE)
indAxis <- reactive({
tryCatch(eval(parse(text = paste("c(", input$indAxis, ")"))),
error = function(e){NULL})
})
output$defn <- renderText({
paste0("Proportion surviving for at least ", input$targetTime,
" ", input$timeUnit,
' in the treatment group "',
input$treatmentGroup,'".')
})
output$EnterQuartilesGp1 <- renderUI({
conditionalPanel(
condition = "output.tgpNumber == 1",
shinyMatrix::matrixInput(
inputId = "myvals1", value = initialVals1(),
class = "numeric",
cols = list(names = TRUE,
editableNames = TRUE),
rows = list(names = TRUE))
)
})
output$EnterQuartilesGp2 <- renderUI({
conditionalPanel(
condition = "output.tgpNumber == 2",
shinyMatrix::matrixInput(inputId = "myvals2", value = initialVals2(),
class = "numeric",
cols = list(names = TRUE,
editableNames = TRUE),
rows = list(names = TRUE))
)
})
output$individualQuartiles <- renderPlot({
req(input$treatmentGroup, survivalDF()$treatment,
indAxis())
p1 <- NULL
if(input$treatmentGroup == levels(survivalDF()$treatment)[1]){
req(input$myvals1)
if(checkPlot(input$myvals1)== TRUE){
if(input$elicMethod == "quartiles"){
p1 <- plotQuartiles(vals = input$myvals1[2:4, ],
lower = input$myvals1[1, ],
upper = input$myvals1[5, ],
expertnames = colnames(input$myvals1),
xlabel = paste0("S(T=",input$targetTime,")"))
}
if(input$elicMethod == "tertiles"){
p1 <- plotTertiles(vals = input$myvals1[2:4, ],
lower = input$myvals1[1, ],
upper = input$myvals1[5, ],
expertnames = colnames(input$myvals1),
xlabel = paste0("S(T=",input$targetTime,")"))
}
}}
if(input$treatmentGroup != levels(survivalDF()$treatment)[1]){
req(input$myvals2)
if(checkPlot(input$myvals2)== TRUE){
if(input$elicMethod == "quartiles"){
p1 <- plotQuartiles(vals = input$myvals2[2:4, ],
lower = input$myvals2[1, ],
upper = input$myvals2[5, ],
expertnames = colnames(input$myvals2),
xlabel = paste0("S(T=",input$targetTime,")"))
}
if(input$elicMethod == "tertiles"){
p1 <- plotTertiles(vals = input$myvals2[2:4, ],
lower = input$myvals2[1, ],
upper = input$myvals2[5, ],
expertnames = colnames(input$myvals2),
xlabel = paste0("S(T=",input$targetTime,")"))
}
}
}
req(p1)
if(input$showScenario == TRUE){
if(input$treatmentGroup == levels(survivalDF()$treatment)[1]){
interval1 <- tryCatch(eval(scenarioIntervalGp1()), error = function(e){NULL})
if(!is.null(interval1)){
p1 <- p1 + geom_vline(xintercept = scenarioIntervalGp1(),
linetype = "dashed")
}
}
if(input$treatmentGroup != levels(survivalDF()$treatment)[1]){
interval2 <- tryCatch(eval(scenarioIntervalGp1()), error = function(e){NULL})
if(!is.null(interval2)){
p1 <- p1 + geom_vline(xintercept = scenarioIntervalGp2(),
linetype = "dashed")
}
}
}
print(p1 + theme_bw(base_size = input$fs)+
xlim(indAxis()[1], indAxis()[2]) + coord_flip())
})
# RIO tab ----
output$RIOdefn <- renderText({
paste0("Proportion surviving for at least ", input$targetTime,
" ", input$timeUnit,
' in the treatment group "',
input$RIOtreatmentGroup,'".')
})
output$RIOtreatmentGroup <- renderUI({
selectInput("RIOtreatmentGroup", "Treatment group",
choices = levels(survivalDF()$treatment))
})
output$RIOtgpNumber <- reactive({
which(input$RIOtreatmentGroup == levels(survivalDF()$treatment))
})
outputOptions(output, 'RIOtgpNumber', suspendWhenHidden = FALSE)
output$RIOvaluesGp1 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 1",
textInput("values1", label = h5("QoI values"),
value = "0.25, 0.5, 0.75")
)
})
output$xaxis1 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 1",
textInput("xaxis1", label = h5(paste0("x-axis limits (",
levels(survivalDF()$treatment)[1], ")")),
value = "0, 1")
)
})
output$xaxis2 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 2",
textInput("xaxis2", label = h5(paste0("x-axis limits (",
levels(survivalDF()$treatment)[2], ")")),
value = "0, 1")
)
})
output$RIOprobsGp1 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 1",
textInput("probs1", label = h5("QoI cumulative probabilities"),
value = "0.25, 0.5, 0.75")
)
})
output$RIOvaluesGp2 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 2",
textInput("values2", label = h5("QoI values"),
value = "0.25, 0.5, 0.75")
)
})
output$RIOprobsGp2 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 2",
textInput("probs2", label = h5("QoI cumulative probabilities"),
value = "0.25, 0.5, 0.75")
)
})
output$RIOdist1 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 1",
selectInput("RIOdist1", label = h5("QoI distribution"),
choices = list(Histogram = "hist",
Normal = "normal",
'Student-t' = "t",
'Skew normal' = "skewnormal",
Gamma = "gamma",
'Log normal' = "lognormal",
'Log Student-t' = "logt",
Beta = "beta",
'Mirror gamma' = "mirrorgamma",
'Mirror log normal' = "mirrorlognormal",
'Mirror log Student-t' = "mirrorlogt",
'Best fitting' = "best"),
#choiceValues = 1:8,
selected = 1)
)
})
output$RIOdist2 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 2",
selectInput("RIOdist2", label = h5("Distribution"),
choices = list(Histogram = "hist",
Normal = "normal",
'Student-t' = "t",
'Skew normal' = "skewnormal",
Gamma = "gamma",
'Log normal' = "lognormal",
'Log Student-t' = "logt",
Beta = "beta",
'Mirror gamma' = "mirrorgamma",
'Mirror log normal' = "mirrorlognormal",
'Mirror log Student-t' = "mirrorlogt",
'Best fitting' = "best"),
#choiceValues = 1:8,
selected = 1)
)
})
# Hack to avoid CRAN check NOTE
X1 <- X2 <- xpos <- ypos <- hjustvar <- vjustvar <- annotateText <- NULL
p1 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$probs1, ")"))),
error = function(e){NULL})
})
p2 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$probs2, ")"))),
error = function(e){NULL})
})
v1 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$values1, ")"))),
error = function(e){NULL})
})
v2 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$values2, ")"))),
error = function(e){NULL})
})
fq <- reactive({
tryCatch(eval(parse(text = paste("c(", input$RIOfq, ")"))),
error = function(e){NULL})
})
xaxis1 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$xaxis1, ")"))),
error = function(e){NULL})
})
xaxis2 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$xaxis2, ")"))),
error = function(e){NULL})
})
myfit1 <- reactive({
req(v1(), p1())
check <- checkJudgementsValid(probs = p1(), vals = v1(),
tdf = 10,
lower = 0,
upper= 1)
if(check$valid == TRUE){
fitdist(vals = v1(), probs = p1(), lower = 0,
upper = 1,
tdf = 10)
}
})
myfit2 <- reactive({
req( v2(), p2())
check <- checkJudgementsValid(probs = p2(), vals = v2(),
tdf = 10,
lower = 0,
upper= 1)
if(check$valid == TRUE){
fitdist(vals = v2(), probs = p2(), lower = 0,
upper = 1,
tdf = 10)
}
})
# RIO pdf plots ----
pdfPlot1 <- reactive({
req(myfit1())
plotfit(myfit1(), d = input$RIOdist1,
ql = fq()[1], qu = fq()[2],
xl = xaxis1()[1], xu = xaxis1()[2],
fs = input$fs,
returnPlot = TRUE,
showPlot = FALSE,
ylab = expression(f[S(T)](x)))
})
pdfPlot2 <- reactive({
req(myfit2())
plotfit(myfit2(), d = input$RIOdist2,
ql = fq()[1], qu = fq()[2],
xl = xaxis2()[1], xu = xaxis2()[2],
fs = input$fs,
returnPlot = TRUE,
showPlot = FALSE,
ylab = expression(f[S(T)](x)))
})
output$pdfPlot <- renderPlot({
req(survivalDF(), input$RIOtreatmentGroup)
if(input$RIOtreatmentGroup == levels(survivalDF()$treatment)[1]){
req(pdfPlot1())
print(pdfPlot1())
}
if(input$RIOtreatmentGroup != levels(survivalDF()$treatment)[1]){
req(pdfPlot2())
print(pdfPlot2())
}
})
# RIO cdf plots ----
cdfPlot1 <- reactive({
req(myfit1(), v1(), p1(), input$RIOdist1)
makeCDFPlot(lower = 0,
v = v1(),
p = p1(),
upper = 1,
fit = myfit1(),
dist = input$RIOdist1,
showFittedCDF = TRUE,
fontsize = input$fs,
xaxisLower = xaxis1()[1], xaxisUpper = xaxis1()[2],
ylab = expression(P(S(T) <= x)))
})
cdfPlot2 <- reactive({
req(myfit2(), v2(), p2())
makeCDFPlot(lower = 0,
v = v2(),
p = p2(),
upper = 1,
fit = myfit2(),
dist = input$RIOdist2,
showFittedCDF = TRUE,
fontsize = input$fs,
xaxisLower = xaxis2()[1], xaxisUpper = xaxis2()[2],
ylab = expression(P(S(T) <= x)))
})
output$cdfPlot <- renderPlot({
req(survivalDF(), input$RIOtreatmentGroup)
if(input$RIOtreatmentGroup == levels(survivalDF()$treatment)[1]){
req(cdfPlot1())
print(cdfPlot1())
}
if(input$RIOtreatmentGroup != levels(survivalDF()$treatment)[1]){
req(cdfPlot2())
print(cdfPlot2())
}
})
quantileValues <- reactive({
req(fq(), input$RIOtreatmentGroup, survivalDF())
if(input$RIOtreatmentGroup == levels(survivalDF()$treatment)[1]){
req(myfit1())
myfit <- myfit1()
dist <- input$RIOdist1
}
if(input$RIOtreatmentGroup != levels(survivalDF()$treatment)[1]){
req(myfit2())
myfit <- myfit2()
dist <- input$RIOdist2
}
if(min(fq())<=0 | max(fq())>=1){
return(NULL)
}else{
FB <- feedback(myfit,
quantiles = fq(),
ex = 1)
if(dist == "best"){
values <- FB$fitted.quantiles[,
as.character(myfit$best.fitting[1,
1])]
}else{
values <- FB$fitted.quantiles[, dist]
}
return(values)
}
})
probabilityValues <- reactive({
req(survivalDF(), input$RIOtreatmentGroup,
input$RIOfeedbackProbability)
if(input$RIOtreatmentGroup == levels(survivalDF()$treatment)[1]){
req(myfit1())
myfit <- myfit1()
dist <- input$RIOdist1
}
if(input$RIOtreatmentGroup != levels(survivalDF()$treatment)[1]){
req(myfit2())
myfit <- myfit2()
dist <- input$RIOdist2
}
FB <- feedback(myfit,
values = input$RIOfeedbackProbability,
ex = 1)
if(dist == "best"){
probs <- FB$fitted.probabilities[,
as.character(myfit$best.fitting[1,
1])]
}else{
probs <- FB$fitted.probabilities[, dist]
}
return(probs)
})
# ...and display on the PDF tab...
output$RIOfq <- renderText({
req(quantileValues())
paste0("Fitted quantiles: ", quantileValues()[1], ", ", quantileValues()[2])
})
output$RIOfp <- renderText({
req(probabilityValues())
#paste0("P(QoI <= ", input$RIOfeedbackProbability,") = ", probabilityValues())
paste0(" = ", probabilityValues())
})
# Extrapolation plot tab ----
output$extrapolationPlot <- renderPlot({
req(survivalDF(),
input$targetTime, myfit1(), input$truncationTime,
input$alpha > 0, input$alpha < 1)
if(input$RIOdist1 == "best"){
fqRIO1 <- myfit1()$best.fitting[1, 1]
}else{
fqRIO1 <- input$RIOdist1
}
if(length(levels(survivalDF()$treatment)) > 1){
group2RIO <- myfit2()
if(input$RIOdist2 == "best"){
fqRIO2 <- myfit2()$best.fitting[1, 1]
}else{
fqRIO2 <- input$RIOdist2
}
}else{
group2RIO <- NULL
fqRIO2 <- NULL
}
survivalExtrapolatePlot(survivalDF(),
myfit1 = myfit1(),
myfit2 = group2RIO,
fqDist1 = fqRIO1,
fqDist2 = fqRIO2,
tTruncate = input$truncationTime,
tTarget = input$targetTime,
alpha = input$alpha,
useWeights = caseWeight_value(),
groups = levels(survivalDF()$treatment),
xl = paste0("Time t (", input$timeUnit, ")"),
fontsize = input$fs,
breakTime = input$targetTime/8,
showPlot = TRUE,
returnPlot = FALSE)
})
observeEvent(input$exit, {
stopApp()
})
# Report writing ----
output$reportGroup1 <- renderUI({
req(survivalDF())
checkboxInput("reportGroup1", paste0('Results for treatment group "',
levels(survivalDF()$treatment)[1], '"'), value = TRUE)
})
output$reportGroup2 <- renderUI({
req(survivalDF())
if(length(levels(survivalDF()$treatment)) >1){
checkboxInput("reportGroup2", paste0('Results for treatment group "',
levels(survivalDF()$treatment)[2], '"'), value = TRUE)
}
})
output$report <- downloadHandler(
filename = function(){switch(input$outFormat,
html_document = "extrapolation-report.html",
pdf_document = "extrapolation-report.pdf",
word_document = "extrapolation-report.docx")},
content = function(file) {
# Copy the report file to a temporary directory before processing it, in
# case we don't have write permissions to the current working dir (which
# can happen when deployed).
tempReport <- file.path(tempdir(), "elicitationShinySummarySurvivalExtrapolation.Rmd")
file.copy(system.file("shinyAppFiles", "elicitationShinySummarySurvivalExtrapolation.Rmd",
package="SHELF"),
tempReport, overwrite = TRUE)
# Set up parameters to pass to Rmd document
dist1 <- input$RIOdist1
dist2 <- NULL
reportGroup1 <- input$reportGroup1
reportGroup2 <- FALSE
if(length(levels(survivalDF()$treatment)) == 1){
allfits <- list(myfit1())
individual <- list(input$myvals1)
}
if(length(levels(survivalDF()$treatment)) == 2){
allfits <- list(myfit1(), myfit2())
individual <- list(input$myvals1, input$myvals2)
reportGroup2 <- input$reportGroup2
dist2 <- input$RIOdist2
}
names(allfits) <- names(individual) <- levels(survivalDF()$treatment)
params <- list(allfits = allfits, individual = individual,
survivalDF = survivalDF(),
targetTime = input$targetTime,
timeUnit = input$timeUnit,
truncationTime = input$truncationTime,
breakTime = input$breakTime,
reportGroup1 = reportGroup1,
reportGroup2 = reportGroup2,
reportDistributions = input$reportDistributions,
dist1 = dist1,
dist2 = dist2,
inputMethod = input$elicMethod,
reportScenarioTest = input$reportScenarioTest,
expRange = expRange(),
reportExtrapolation = input$reportExtrapolation,
useWeights = caseWeight_value(),
alpha = input$alpha)
# Knit the document, passing in the `params` list, and eval it in a
# child of the global environment (this isolates the code in the document
# from the code in this app).
rmarkdown::render(tempReport, output_file = file,
params = params,
output_format = input$outFormat,
envir = new.env(parent = globalenv())
)
}
)
}
), launch.browser = TRUE)
}
# quantileValues1 <- reactive({
# ssq <- myfit1()$ssq[1, is.na(myfit1()$ssq[1,])==F]
# best.index <- which(ssq == min(ssq))[1]
#
# ex <- 1
# pl <- limits1()[1]
# pu <- limits1()[2]
# if(as.numeric(input$radio1)==8){index<-best.index}else{index<-as.numeric(input$radio1) - 1}
# if(as.numeric(input$radio1)==1){
# if(pl == -Inf & myfit1()$limits[ex,1] > -Inf){pl <- myfit1()$limits[ex,1]}
# if(pu == Inf & myfit1()$limits[ex,2] < Inf){pu <- myfit1()$limits[ex,2] }
# if(pl == -Inf & myfit1()$limits[ex,1] == -Inf){
# pl <- qnorm(0.001, myfit1()$Normal[ex,1], myfit1()$Normal[ex,2])}
# if(pu == Inf & myfit1()$limits[ex,2] == Inf){
# pu <- qnorm(0.999, myfit1()$Normal[ex,1], myfit1()$Normal[ex,2])}
# p <- c(0, myfit1()$probs[ex,], 1)
# x <- c(pl, myfit1()$vals[ex,], pu)
# values <- qhist(c(0.05, 0.95), x, p)
# }
#
# if(as.numeric(input$radio1)>1){
# temp<-feedback(myfit1(), quantiles=c(0.05, 0.95), ex=1)
# values=temp$fitted.quantiles[,index]
# }
# data.frame(quantiles=c(0.05, 0.95), values=values)
#
# })
#
#
# quantileValues2 <- reactive({
# ssq <- myfit2()$ssq[1, is.na(myfit2()$ssq[1,])==F]
# best.index <- which(ssq == min(ssq))[1]
#
# ex <- 1
# pl <- limits1()[1]
# pu <- limits1()[2]
# if(as.numeric(input$radio2)==8){index<-best.index}else{index<-as.numeric(input$radio2) - 1}
# if(as.numeric(input$radio2)==1){
# if(pl == -Inf & myfit2()$limits[ex,1] > -Inf){pl <- myfit2()$limits[ex,1]}
# if(pu == Inf & myfit2()$limits[ex,2] < Inf){pu <- myfit2()$limits[ex,2] }
# if(pl == -Inf & myfit2()$limits[ex,1] == -Inf){
# pl <- qnorm(0.001, myfit2()$Normal[ex,1], myfit2()$Normal[ex,2])}
# if(pu == Inf & myfit2()$limits[ex,2] == Inf){
# pu <- qnorm(0.999, myfit2()$Normal[ex,1], myfit2()$Normal[ex,2])}
# p <- c(0, myfit2()$probs[ex,], 1)
# x <- c(pl, myfit2()$vals[ex,], pu)
# values <- qhist(c(0.05, 0.95), x, p)
# }
#
# if(as.numeric(input$radio1)>1){
# temp<-feedback(myfit1(), quantiles=c(0.05, 0.95), ex=1)
# values=temp$fitted.quantiles[,index]
# }
# data.frame(quantiles=c(0.05, 0.95), values=values)
#
# })
#
#
# output$valuesPDF1 <- renderTable({
# quantileValues1()
# })
#
# output$valuesPDF2 <- renderTable({
# quantileValues2()
# })
# Helper function KMextrapolate ----
# Helper function survivalTable ----
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
}
# Helper function check quartiles can be plotted ----
checkPlot <- function(myvals){
if(anyNA(myvals)){
return(FALSE)
}
if(any(diff(myvals)<=0)){
return(FALSE)
}
return(TRUE)
}
OakleyJ/SHELF documentation built on June 9, 2025, 11:05 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions checkPlot makeSurvSummaryTable elicitSurvivalExtrapolation
Documented in elicitSurvivalExtrapolation
#' Elicitation for survival extrapolation
#'
#' Opens up a web browser in which you can implement the SHELF protocol for survival extrapolation.
#' Start with uploading a .csv file of individual patient survival data (time, event to indicate censoring, and
#' treatment group). Then elicit individual judgements, perform scenario testing as required, and elicit a RIO distribution.
#' Judgements for two treatment groups can be elicited in the same session.
#'
#'
#' @author Jeremy Oakley <j.oakley@@sheffield.ac.uk>
#' @examples
#'
#' \dontrun{
#'
#' # make a suitable csv file using a built in data set from the survival package
#' sdf <- survival::veteran[, c("time", "status", "trt")]
#' colnames(sdf) <- c("time", "event", "treatment")
#' sdf$treatment <- factor(sdf$treatment, labels = c("standard", "test"))
#'
#' # write the data frame sdf to a .csv file in the current working directory
#' write.csv(sdf, file = "testFile.csv", row.names = FALSE)
#'
#' # Run the app and upload testFile.csv in the first tab, and change unit of time to "days"
#'
#' elicitSurvivalExtrapolation()
#'
#' }
#' @import shiny
#' @import survival
#' @importFrom survminer ggsurvplot
#' @export
elicitSurvivalExtrapolation<- function(){
runApp(list(
ui = shinyUI(fluidPage(
# Application title
titlePanel("SHELF: survival extrapolation"),
# sidebarLayout(
mainPanel(tags$style(type="text/css",
".shiny-output-error { visibility: hidden; }",
".shiny-output-error:before { visibility: hidden; }"
),
tabsetPanel(
tabPanel("Upload data",
sidebarLayout(
sidebarPanel(
wellPanel(
fileInput("file", label = "Survival data upload"),
uiOutput("targetTime"),
uiOutput("truncationTime"),
uiOutput("timeUnit"),
uiOutput("caseWeights"),
uiOutput("fontSize")
)
)
,
mainPanel(
helpText('Upload a .csv file with three columns:
"time", "event" and "treatment" (in that order).
Values in the "event" column should be 0 for
a censored observation, and 1 otherwise. The
"treatment" column should be included even
if there is only one treatment group. If the
data include case weights for each observation,
include a fourth column "weights"'),
plotOutput("KMbasic"),
tableOutput("survivalSummary"),
)
)
),
tabPanel("KM table",
sidebarLayout(
sidebarPanel(
wellPanel(
uiOutput("breakTime"),
)),
mainPanel(
helpText('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'),
tableOutput("KMdata")
)
)
),
tabPanel("Individual judgements",
sidebarLayout(
sidebarPanel(
wellPanel(
uiOutput("treatmentGroup"),
numericInput("nExperts", label = h5("Number of experts"),
value = 2, min = 1),
radioButtons("elicMethod", "Elicitation method",
c("Quartiles" = "quartiles",
"Tertiles" = "tertiles")),
textInput("indAxis", label = h5("Axis limits"),
value = "0, 1"),
hr(style = "border-top: 1px solid #000000;"),
h5("Scenario testing"),
helpText("Use next tab before displaying here"),
checkboxInput("showScenario",
label = h5("Show scenario interval"))
)),
mainPanel(
h4("Quantity of interest:"),
textOutput("defn"),
hr(),
helpText('Enter plausible limits and quantiles as values
between 0 and 1.'),
uiOutput("EnterQuartilesGp1"),
uiOutput("EnterQuartilesGp2"),
plotOutput("individualQuartiles")
)
)
),
tabPanel("Scenario testing",
sidebarLayout(
sidebarPanel(
wellPanel(
uiOutput("scenarioGroup"),
uiOutput("scenarioTruncationTime"),
uiOutput("expRange")
)),
mainPanel(
helpText('In this tab you can consider a scenario in which the hazard
remains constant after a specified time point. The shaded
bar indicates a point-wise approximate 95% credible interval for
the survivor function based on this assumption. Intervals can be
added to the "Individual judgements" tab.'),
plotOutput("KMplot"),
htmlOutput("reportInterval")
)
)
),
tabPanel("RIO elicitation",
sidebarLayout(
sidebarPanel(
wellPanel(
uiOutput("RIOtreatmentGroup"),
uiOutput("RIOvaluesGp1"),
uiOutput("RIOprobsGp1"),
uiOutput("RIOvaluesGp2"),
uiOutput("RIOprobsGp2"),
uiOutput("RIOdist1"),
uiOutput("RIOdist2"),
hr(style = "border-top: 1px solid #000000;"),
textInput("RIOfq", label = h5("Feedback quantiles"),
value = "0.1, 0.9"),
textOutput("RIOfq"),
h5("Feedback probability"),
fluidRow(
column(width = 7,
"Probability QoI <= "
),
column(width = 5,
numericInput("RIOfeedbackProbability", label = NULL,
value = 0.5, min = 0, max = 1)
)
),
textOutput("RIOfp"),
hr(style = "border-top: 1px solid #000000;"),
uiOutput("xaxis1"),
uiOutput("xaxis2")
)),
mainPanel(
h4("Quantity of Interest (QoI):"),
textOutput("RIOdefn"),
hr(),
plotOutput("pdfPlot"),
plotOutput("cdfPlot")
)
)
),
tabPanel("Extrapolation plot",
sidebarLayout(
sidebarPanel(
wellPanel(
numericInput("alpha", "Credible Interval Size",
0.95, min = 0, max = 1),
)),
mainPanel(
plotOutput("extrapolationPlot")
)
)
)
,
tabPanel("Report",
wellPanel(
h5("Report content"),
checkboxInput("reportData", "KM plot and summary data",
value = TRUE, width = NULL),
checkboxInput("reportTable", "Survivor and hazard table",
value = TRUE, width = NULL),
uiOutput("reportGroup1"),
uiOutput("reportGroup2"),
checkboxInput("reportIndividual", "Individual judgements",
value = TRUE, width = NULL),
checkboxInput("reportScenarioTest", "Scenario Test results",
value = TRUE, width = NULL),
checkboxInput("reportExtrapolation", "Plot extrapolation",
value = TRUE, width = NULL),
checkboxInput("reportDistributions", "All fitted distributions",
value = TRUE, width = NULL),
selectInput("outFormat", label = "Report format",
choices = list('html' = "html_document",
'pdf' = "pdf_document",
'Word' = "word_document"),
width = "30%"),
downloadButton("report", "Download report")
)
)
)
)
)
),
server = function(input, output) {
# File uploading ----
survivalDF <- reactive({
req(input$file)
inFile <- input$file
sdf <- utils::read.csv(inFile$datapath)
colnames(sdf)[1:3] <- c("time", "event", "treatment")
sdf$treatment <- as.factor(sdf$treatment)
if(length(levels(sdf$treatment)) > 2){
showNotification(
"This app will only with work with one or two treatment groups.",
duration = 60,
type = "error"
)
sdf <- NULL
}
sdf
})
output$KMbasic <- renderPlot({
req(survivalDF(),input$targetTime, input$truncationTime)
# Truncate data frame for plotting
index <- survivalDF()$time > input$truncationTime
truncatedDf <- survivalDF()
truncatedDf$event[index] <- 0
truncatedDf$time[index] <- input$truncationTime
if(caseWeight_value() == TRUE){
fit <- survival::survfit(survival::Surv(time, event) ~ treatment,
weights = weights,
data = truncatedDf)}else{
fit <- survival::survfit(survival::Surv(time, event) ~ treatment, data = truncatedDf)
}
myplot<- survminer::ggsurvplot(fit, data = truncatedDf, censor = TRUE,
legend = "right",
legend.title = "",
conf.int = TRUE,
legend.labs = levels(truncatedDf$treatment),
xlim = c(0, input$targetTime),
xlab = paste0("Time t (", input$timeUnit, ")"),
ylab = "S(t)",
break.time.by = input$targetTime/8)
myplot$plot +
geom_vline(xintercept = input$targetTime, linetype="dotted") +
theme_bw(base_size = input$fs)
})
output$survivalSummary <-renderTable({
req(survivalDF())
makeSurvSummaryTable(survivalDF(),
useWeights = caseWeight_value())},
rownames = TRUE)
output$targetTime <- renderUI({
req(survivalDF())
maxTime <- signif(max(survivalDF()$time), 1)
numericInput("targetTime", "Target extrapolation time", value = 2* maxTime)
})
output$truncationTime <- renderUI({
req(survivalDF())
maxTime <- signif(max(survivalDF()$time), 1)
numericInput("truncationTime", "Data truncation time", value = 2* maxTime)
})
output$timeUnit <- renderUI({
req(survivalDF())
selectInput("timeUnit", "Unit of time",
c("days", "weeks", "months", "years"),
selected = "years")
})
output$caseWeights <- renderUI({
req(survivalDF())
if("weights" %in% colnames(survivalDF())){
checkboxInput("useWeights", "Use case weights", value = TRUE)
}
})
caseWeight_value <- reactive({
if (is.null(input$useWeights)) {
return(FALSE)
} else {
return(input$useWeights)
}
})
output$fontSize <- renderUI({
req(survivalDF())
numericInput("fs", label = "Font size (all plots)", value = 16)
})
# KM table tab ----
output$breakTime <- renderUI({
req(survivalDF(), input$timeUnit)
maxTime <- min(max(survivalDF()$time), input$truncationTime)
numericInput("breakTime", paste0("Time interval duration (",
input$timeUnit, ")"),
value = signif(maxTime/4, 1))
})
output$KMdata <- renderTable({
req(survivalDF(), input$breakTime,
input$truncationTime)
makeSurvivalTable(survivalDF(), input$breakTime,
input$truncationTime, input$timeUnit, dp = 2,
useWeights = caseWeight_value())
})
# Scenario testing tab ----
output$scenarioGroup <- renderUI({
selectInput("scenarioGroup", "Treatment group",
choices = levels(survivalDF()$treatment))
})
output$expRange <- renderUI({
req(survivalDF())
maxTime <- min(max(survivalDF()$time), input$truncationTime)
textInput("expRange", label = h5("Constant hazard fitting period"),
value = paste0(0, ", ", signif(maxTime, 1) ))
})
expRange <- reactive({
tryCatch(eval(parse(text = paste("c(", input$expRange, ")"))),
error = function(e){NULL})
})
scenario <- reactiveValues()
observe({
req(survivalDF(), expRange(), input$scenarioGroup,
input$targetTime,
input$truncationTime)
sce <- survivalScenario(tLower = 0,
tUpper = min(input$truncationTime,
max(survivalDF()$time)),
expLower= expRange()[1],
expUpper = expRange()[2],
expGroup = input$scenarioGroup,
tTarget = input$targetTime,
survDf = survivalDF(),
groups = levels(survivalDF()$treatment),
xl = paste0("Time t (", input$timeUnit, ")"),
showPlot = FALSE,
fontsize = input$fs,
useWeights = caseWeight_value())
scenario$plot <- sce$KMplot
scenario$interval <- sce$interval
})
output$KMplot <- renderPlot({
# req(survivalDF(), expRange(), input$scenarioGroup,
# input$targetTime,
# input$truncationTime)
# scenario <- survivalScenario(tLower = 0,
# tUpper = min(input$truncationTime,
# max(survivalDF()$time)),
# expLower= expRange()[1],
# expUpper = expRange()[2],
# expGroup = input$scenarioGroup,
# tTarget = input$targetTime,
# survDf = survivalDF(),
# groups = levels(survivalDF()$treatment),
# xl = paste0("Time (", input$timeUnit, ")"),
# showPlot = FALSE,
# fontsize = input$fs)
# scenario$KMplot
req(scenario$plot)
scenario$plot
})
scenarioIntervalGp1 <- reactive({
req(survivalDF(), expRange(),input$scenarioGroup,
input$targetTime,
input$truncationTime)
sce <- survivalScenario(tLower = 0,
tUpper = min(input$truncationTime,
max(survivalDF()$time)),
expLower= expRange()[1],
expUpper = expRange()[2],
expGroup = levels(survivalDF()$treatment)[1],
tTarget = input$targetTime,
survDf = survivalDF(),
groups = levels(survivalDF()$treatment),
xl = paste0("Time (", input$timeUnit, ")"),
showPlot = FALSE,
fontsize = input$fs,
useWeights = caseWeight_value())
sce$interval
})
scenarioIntervalGp2 <- reactive({
req(survivalDF(), expRange(),input$scenarioGroup,
input$targetTime,
input$truncationTime)
interval <- NULL
if(length(levels(survivalDF()$treatment)) > 1){
interval <- survivalScenario(tLower = 0,
tUpper = min(input$truncationTime,
max(survivalDF()$time)),
expLower= expRange()[1],
expUpper = expRange()[2],
expGroup = levels(survivalDF()$treatment)[2],
tTarget = input$targetTime,
survDf = survivalDF(),
groups = levels(survivalDF()$treatment),
xl = paste0("Time (", input$timeUnit, ")"),
showPlot = FALSE,
fontsize = input$fs,
useWeights = caseWeight_value())$interval
}
interval
})
output$reportInterval <- renderUI({
# req(input$scenarioGroup, survivalDF())
# txt <- NULL
# if(input$scenarioGroup == levels(survivalDF()$treatment)[1]){
# req(scenarioIntervalGp1())
# txt <- paste0("(", signif(scenarioIntervalGp1()[1], 2)
# , ", ", signif(scenarioIntervalGp1()[2], 2), ")" )
# }
#
# if(input$scenarioGroup != levels(survivalDF()$treatment)[1]){
# req(scenarioIntervalGp2())
# txt <- paste0("(", signif(scenarioIntervalGp2()[1], 2), ", ",
# signif(scenarioIntervalGp2()[2], 2), ")" )
# }
req(scenario$interval)
txt <- paste0("(", signif(scenario$interval[1], 2),
", ", signif(scenario$interval[2], 2), ")" )
HTML(paste(hr(),
h5(paste0("Approximate 95% credible interval at the target time: ",
txt))))
})
# Individual judgements tab ----
initialVals1 <- reactive({
initialdf <- matrix(c(0, 0.25, 0.5, 0.75, 1), nrow = 5, ncol = input$nExperts
)
if(input$elicMethod == "quartiles"){
rownames(initialdf) <- c("L", "Q1", "M", "Q3", "U")
}
if(input$elicMethod == "tertiles"){
rownames(initialdf) <- c("L", "T1", "M", "T2", "U")
}
colnames(initialdf) <- LETTERS[1:input$nExperts]
initialdf
})
initialVals2 <- reactive({
initialdf <- matrix(0.1*(1:5), nrow = 5, ncol = input$nExperts
)
if(input$elicMethod == "quartiles"){
rownames(initialdf) <- c("L", "Q1", "M", "Q3", "U")
}
if(input$elicMethod == "tertiles"){
rownames(initialdf) <- c("L", "T1", "M", "T2", "U")
}
colnames(initialdf) <- LETTERS[1:input$nExperts]
initialdf
})
output$treatmentGroup <- renderUI({
selectInput("treatmentGroup", "Treatment group",
choices = levels(survivalDF()$treatment))
})
output$tgpNumber <- reactive({
which(input$treatmentGroup == levels(survivalDF()$treatment))
})
outputOptions(output, 'tgpNumber', suspendWhenHidden = FALSE)
indAxis <- reactive({
tryCatch(eval(parse(text = paste("c(", input$indAxis, ")"))),
error = function(e){NULL})
})
output$defn <- renderText({
paste0("Proportion surviving for at least ", input$targetTime,
" ", input$timeUnit,
' in the treatment group "',
input$treatmentGroup,'".')
})
output$EnterQuartilesGp1 <- renderUI({
conditionalPanel(
condition = "output.tgpNumber == 1",
shinyMatrix::matrixInput(
inputId = "myvals1", value = initialVals1(),
class = "numeric",
cols = list(names = TRUE,
editableNames = TRUE),
rows = list(names = TRUE))
)
})
output$EnterQuartilesGp2 <- renderUI({
conditionalPanel(
condition = "output.tgpNumber == 2",
shinyMatrix::matrixInput(inputId = "myvals2", value = initialVals2(),
class = "numeric",
cols = list(names = TRUE,
editableNames = TRUE),
rows = list(names = TRUE))
)
})
output$individualQuartiles <- renderPlot({
req(input$treatmentGroup, survivalDF()$treatment,
indAxis())
p1 <- NULL
if(input$treatmentGroup == levels(survivalDF()$treatment)[1]){
req(input$myvals1)
if(checkPlot(input$myvals1)== TRUE){
if(input$elicMethod == "quartiles"){
p1 <- plotQuartiles(vals = input$myvals1[2:4, ],
lower = input$myvals1[1, ],
upper = input$myvals1[5, ],
expertnames = colnames(input$myvals1),
xlabel = paste0("S(T=",input$targetTime,")"))
}
if(input$elicMethod == "tertiles"){
p1 <- plotTertiles(vals = input$myvals1[2:4, ],
lower = input$myvals1[1, ],
upper = input$myvals1[5, ],
expertnames = colnames(input$myvals1),
xlabel = paste0("S(T=",input$targetTime,")"))
}
}}
if(input$treatmentGroup != levels(survivalDF()$treatment)[1]){
req(input$myvals2)
if(checkPlot(input$myvals2)== TRUE){
if(input$elicMethod == "quartiles"){
p1 <- plotQuartiles(vals = input$myvals2[2:4, ],
lower = input$myvals2[1, ],
upper = input$myvals2[5, ],
expertnames = colnames(input$myvals2),
xlabel = paste0("S(T=",input$targetTime,")"))
}
if(input$elicMethod == "tertiles"){
p1 <- plotTertiles(vals = input$myvals2[2:4, ],
lower = input$myvals2[1, ],
upper = input$myvals2[5, ],
expertnames = colnames(input$myvals2),
xlabel = paste0("S(T=",input$targetTime,")"))
}
}
}
req(p1)
if(input$showScenario == TRUE){
if(input$treatmentGroup == levels(survivalDF()$treatment)[1]){
interval1 <- tryCatch(eval(scenarioIntervalGp1()), error = function(e){NULL})
if(!is.null(interval1)){
p1 <- p1 + geom_vline(xintercept = scenarioIntervalGp1(),
linetype = "dashed")
}
}
if(input$treatmentGroup != levels(survivalDF()$treatment)[1]){
interval2 <- tryCatch(eval(scenarioIntervalGp1()), error = function(e){NULL})
if(!is.null(interval2)){
p1 <- p1 + geom_vline(xintercept = scenarioIntervalGp2(),
linetype = "dashed")
}
}
}
print(p1 + theme_bw(base_size = input$fs)+
xlim(indAxis()[1], indAxis()[2]) + coord_flip())
})
# RIO tab ----
output$RIOdefn <- renderText({
paste0("Proportion surviving for at least ", input$targetTime,
" ", input$timeUnit,
' in the treatment group "',
input$RIOtreatmentGroup,'".')
})
output$RIOtreatmentGroup <- renderUI({
selectInput("RIOtreatmentGroup", "Treatment group",
choices = levels(survivalDF()$treatment))
})
output$RIOtgpNumber <- reactive({
which(input$RIOtreatmentGroup == levels(survivalDF()$treatment))
})
outputOptions(output, 'RIOtgpNumber', suspendWhenHidden = FALSE)
output$RIOvaluesGp1 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 1",
textInput("values1", label = h5("QoI values"),
value = "0.25, 0.5, 0.75")
)
})
output$xaxis1 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 1",
textInput("xaxis1", label = h5(paste0("x-axis limits (",
levels(survivalDF()$treatment)[1], ")")),
value = "0, 1")
)
})
output$xaxis2 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 2",
textInput("xaxis2", label = h5(paste0("x-axis limits (",
levels(survivalDF()$treatment)[2], ")")),
value = "0, 1")
)
})
output$RIOprobsGp1 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 1",
textInput("probs1", label = h5("QoI cumulative probabilities"),
value = "0.25, 0.5, 0.75")
)
})
output$RIOvaluesGp2 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 2",
textInput("values2", label = h5("QoI values"),
value = "0.25, 0.5, 0.75")
)
})
output$RIOprobsGp2 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 2",
textInput("probs2", label = h5("QoI cumulative probabilities"),
value = "0.25, 0.5, 0.75")
)
})
output$RIOdist1 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 1",
selectInput("RIOdist1", label = h5("QoI distribution"),
choices = list(Histogram = "hist",
Normal = "normal",
'Student-t' = "t",
'Skew normal' = "skewnormal",
Gamma = "gamma",
'Log normal' = "lognormal",
'Log Student-t' = "logt",
Beta = "beta",
'Mirror gamma' = "mirrorgamma",
'Mirror log normal' = "mirrorlognormal",
'Mirror log Student-t' = "mirrorlogt",
'Best fitting' = "best"),
#choiceValues = 1:8,
selected = 1)
)
})
output$RIOdist2 <- renderUI({
conditionalPanel(
condition = "output.RIOtgpNumber == 2",
selectInput("RIOdist2", label = h5("Distribution"),
choices = list(Histogram = "hist",
Normal = "normal",
'Student-t' = "t",
'Skew normal' = "skewnormal",
Gamma = "gamma",
'Log normal' = "lognormal",
'Log Student-t' = "logt",
Beta = "beta",
'Mirror gamma' = "mirrorgamma",
'Mirror log normal' = "mirrorlognormal",
'Mirror log Student-t' = "mirrorlogt",
'Best fitting' = "best"),
#choiceValues = 1:8,
selected = 1)
)
})
# Hack to avoid CRAN check NOTE
X1 <- X2 <- xpos <- ypos <- hjustvar <- vjustvar <- annotateText <- NULL
p1 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$probs1, ")"))),
error = function(e){NULL})
})
p2 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$probs2, ")"))),
error = function(e){NULL})
})
v1 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$values1, ")"))),
error = function(e){NULL})
})
v2 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$values2, ")"))),
error = function(e){NULL})
})
fq <- reactive({
tryCatch(eval(parse(text = paste("c(", input$RIOfq, ")"))),
error = function(e){NULL})
})
xaxis1 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$xaxis1, ")"))),
error = function(e){NULL})
})
xaxis2 <- reactive({
tryCatch(eval(parse(text = paste("c(", input$xaxis2, ")"))),
error = function(e){NULL})
})
myfit1 <- reactive({
req(v1(), p1())
check <- checkJudgementsValid(probs = p1(), vals = v1(),
tdf = 10,
lower = 0,
upper= 1)
if(check$valid == TRUE){
fitdist(vals = v1(), probs = p1(), lower = 0,
upper = 1,
tdf = 10)
}
})
myfit2 <- reactive({
req( v2(), p2())
check <- checkJudgementsValid(probs = p2(), vals = v2(),
tdf = 10,
lower = 0,
upper= 1)
if(check$valid == TRUE){
fitdist(vals = v2(), probs = p2(), lower = 0,
upper = 1,
tdf = 10)
}
})
# RIO pdf plots ----
pdfPlot1 <- reactive({
req(myfit1())
plotfit(myfit1(), d = input$RIOdist1,
ql = fq()[1], qu = fq()[2],
xl = xaxis1()[1], xu = xaxis1()[2],
fs = input$fs,
returnPlot = TRUE,
showPlot = FALSE,
ylab = expression(f[S(T)](x)))
})
pdfPlot2 <- reactive({
req(myfit2())
plotfit(myfit2(), d = input$RIOdist2,
ql = fq()[1], qu = fq()[2],
xl = xaxis2()[1], xu = xaxis2()[2],
fs = input$fs,
returnPlot = TRUE,
showPlot = FALSE,
ylab = expression(f[S(T)](x)))
})
output$pdfPlot <- renderPlot({
req(survivalDF(), input$RIOtreatmentGroup)
if(input$RIOtreatmentGroup == levels(survivalDF()$treatment)[1]){
req(pdfPlot1())
print(pdfPlot1())
}
if(input$RIOtreatmentGroup != levels(survivalDF()$treatment)[1]){
req(pdfPlot2())
print(pdfPlot2())
}
})
# RIO cdf plots ----
cdfPlot1 <- reactive({
req(myfit1(), v1(), p1(), input$RIOdist1)
makeCDFPlot(lower = 0,
v = v1(),
p = p1(),
upper = 1,
fit = myfit1(),
dist = input$RIOdist1,
showFittedCDF = TRUE,
fontsize = input$fs,
xaxisLower = xaxis1()[1], xaxisUpper = xaxis1()[2],
ylab = expression(P(S(T) <= x)))
})
cdfPlot2 <- reactive({
req(myfit2(), v2(), p2())
makeCDFPlot(lower = 0,
v = v2(),
p = p2(),
upper = 1,
fit = myfit2(),
dist = input$RIOdist2,
showFittedCDF = TRUE,
fontsize = input$fs,
xaxisLower = xaxis2()[1], xaxisUpper = xaxis2()[2],
ylab = expression(P(S(T) <= x)))
})
output$cdfPlot <- renderPlot({
req(survivalDF(), input$RIOtreatmentGroup)
if(input$RIOtreatmentGroup == levels(survivalDF()$treatment)[1]){
req(cdfPlot1())
print(cdfPlot1())
}
if(input$RIOtreatmentGroup != levels(survivalDF()$treatment)[1]){
req(cdfPlot2())
print(cdfPlot2())
}
})
quantileValues <- reactive({
req(fq(), input$RIOtreatmentGroup, survivalDF())
if(input$RIOtreatmentGroup == levels(survivalDF()$treatment)[1]){
req(myfit1())
myfit <- myfit1()
dist <- input$RIOdist1
}
if(input$RIOtreatmentGroup != levels(survivalDF()$treatment)[1]){
req(myfit2())
myfit <- myfit2()
dist <- input$RIOdist2
}
if(min(fq())<=0 | max(fq())>=1){
return(NULL)
}else{
FB <- feedback(myfit,
quantiles = fq(),
ex = 1)
if(dist == "best"){
values <- FB$fitted.quantiles[,
as.character(myfit$best.fitting[1,
1])]
}else{
values <- FB$fitted.quantiles[, dist]
}
return(values)
}
})
probabilityValues <- reactive({
req(survivalDF(), input$RIOtreatmentGroup,
input$RIOfeedbackProbability)
if(input$RIOtreatmentGroup == levels(survivalDF()$treatment)[1]){
req(myfit1())
myfit <- myfit1()
dist <- input$RIOdist1
}
if(input$RIOtreatmentGroup != levels(survivalDF()$treatment)[1]){
req(myfit2())
myfit <- myfit2()
dist <- input$RIOdist2
}
FB <- feedback(myfit,
values = input$RIOfeedbackProbability,
ex = 1)
if(dist == "best"){
probs <- FB$fitted.probabilities[,
as.character(myfit$best.fitting[1,
1])]
}else{
probs <- FB$fitted.probabilities[, dist]
}
return(probs)
})
# ...and display on the PDF tab...
output$RIOfq <- renderText({
req(quantileValues())
paste0("Fitted quantiles: ", quantileValues()[1], ", ", quantileValues()[2])
})
output$RIOfp <- renderText({
req(probabilityValues())
#paste0("P(QoI <= ", input$RIOfeedbackProbability,") = ", probabilityValues())
paste0(" = ", probabilityValues())
})
# Extrapolation plot tab ----
output$extrapolationPlot <- renderPlot({
req(survivalDF(),
input$targetTime, myfit1(), input$truncationTime,
input$alpha > 0, input$alpha < 1)
if(input$RIOdist1 == "best"){
fqRIO1 <- myfit1()$best.fitting[1, 1]
}else{
fqRIO1 <- input$RIOdist1
}
if(length(levels(survivalDF()$treatment)) > 1){
group2RIO <- myfit2()
if(input$RIOdist2 == "best"){
fqRIO2 <- myfit2()$best.fitting[1, 1]
}else{
fqRIO2 <- input$RIOdist2
}
}else{
group2RIO <- NULL
fqRIO2 <- NULL
}
survivalExtrapolatePlot(survivalDF(),
myfit1 = myfit1(),
myfit2 = group2RIO,
fqDist1 = fqRIO1,
fqDist2 = fqRIO2,
tTruncate = input$truncationTime,
tTarget = input$targetTime,
alpha = input$alpha,
useWeights = caseWeight_value(),
groups = levels(survivalDF()$treatment),
xl = paste0("Time t (", input$timeUnit, ")"),
fontsize = input$fs,
breakTime = input$targetTime/8,
showPlot = TRUE,
returnPlot = FALSE)
})
observeEvent(input$exit, {
stopApp()
})
# Report writing ----
output$reportGroup1 <- renderUI({
req(survivalDF())
checkboxInput("reportGroup1", paste0('Results for treatment group "',
levels(survivalDF()$treatment)[1], '"'), value = TRUE)
})
output$reportGroup2 <- renderUI({
req(survivalDF())
if(length(levels(survivalDF()$treatment)) >1){
checkboxInput("reportGroup2", paste0('Results for treatment group "',
levels(survivalDF()$treatment)[2], '"'), value = TRUE)
}
})
output$report <- downloadHandler(
filename = function(){switch(input$outFormat,
html_document = "extrapolation-report.html",
pdf_document = "extrapolation-report.pdf",
word_document = "extrapolation-report.docx")},
content = function(file) {
# Copy the report file to a temporary directory before processing it, in
# case we don't have write permissions to the current working dir (which
# can happen when deployed).
tempReport <- file.path(tempdir(), "elicitationShinySummarySurvivalExtrapolation.Rmd")
file.copy(system.file("shinyAppFiles", "elicitationShinySummarySurvivalExtrapolation.Rmd",
package="SHELF"),
tempReport, overwrite = TRUE)
# Set up parameters to pass to Rmd document
dist1 <- input$RIOdist1
dist2 <- NULL
reportGroup1 <- input$reportGroup1
reportGroup2 <- FALSE
if(length(levels(survivalDF()$treatment)) == 1){
allfits <- list(myfit1())
individual <- list(input$myvals1)
}
if(length(levels(survivalDF()$treatment)) == 2){
allfits <- list(myfit1(), myfit2())
individual <- list(input$myvals1, input$myvals2)
reportGroup2 <- input$reportGroup2
dist2 <- input$RIOdist2
}
names(allfits) <- names(individual) <- levels(survivalDF()$treatment)
params <- list(allfits = allfits, individual = individual,
survivalDF = survivalDF(),
targetTime = input$targetTime,
timeUnit = input$timeUnit,
truncationTime = input$truncationTime,
breakTime = input$breakTime,
reportGroup1 = reportGroup1,
reportGroup2 = reportGroup2,
reportDistributions = input$reportDistributions,
dist1 = dist1,
dist2 = dist2,
inputMethod = input$elicMethod,
reportScenarioTest = input$reportScenarioTest,
expRange = expRange(),
reportExtrapolation = input$reportExtrapolation,
useWeights = caseWeight_value(),
alpha = input$alpha)
# Knit the document, passing in the `params` list, and eval it in a
# child of the global environment (this isolates the code in the document
# from the code in this app).
rmarkdown::render(tempReport, output_file = file,
params = params,
output_format = input$outFormat,
envir = new.env(parent = globalenv())
)
}
)
}
), launch.browser = TRUE)
}
# quantileValues1 <- reactive({
# ssq <- myfit1()$ssq[1, is.na(myfit1()$ssq[1,])==F]
# best.index <- which(ssq == min(ssq))[1]
#
# ex <- 1
# pl <- limits1()[1]
# pu <- limits1()[2]
# if(as.numeric(input$radio1)==8){index<-best.index}else{index<-as.numeric(input$radio1) - 1}
# if(as.numeric(input$radio1)==1){
# if(pl == -Inf & myfit1()$limits[ex,1] > -Inf){pl <- myfit1()$limits[ex,1]}
# if(pu == Inf & myfit1()$limits[ex,2] < Inf){pu <- myfit1()$limits[ex,2] }
# if(pl == -Inf & myfit1()$limits[ex,1] == -Inf){
# pl <- qnorm(0.001, myfit1()$Normal[ex,1], myfit1()$Normal[ex,2])}
# if(pu == Inf & myfit1()$limits[ex,2] == Inf){
# pu <- qnorm(0.999, myfit1()$Normal[ex,1], myfit1()$Normal[ex,2])}
# p <- c(0, myfit1()$probs[ex,], 1)
# x <- c(pl, myfit1()$vals[ex,], pu)
# values <- qhist(c(0.05, 0.95), x, p)
# }
#
# if(as.numeric(input$radio1)>1){
# temp<-feedback(myfit1(), quantiles=c(0.05, 0.95), ex=1)
# values=temp$fitted.quantiles[,index]
# }
# data.frame(quantiles=c(0.05, 0.95), values=values)
#
# })
#
#
# quantileValues2 <- reactive({
# ssq <- myfit2()$ssq[1, is.na(myfit2()$ssq[1,])==F]
# best.index <- which(ssq == min(ssq))[1]
#
# ex <- 1
# pl <- limits1()[1]
# pu <- limits1()[2]
# if(as.numeric(input$radio2)==8){index<-best.index}else{index<-as.numeric(input$radio2) - 1}
# if(as.numeric(input$radio2)==1){
# if(pl == -Inf & myfit2()$limits[ex,1] > -Inf){pl <- myfit2()$limits[ex,1]}
# if(pu == Inf & myfit2()$limits[ex,2] < Inf){pu <- myfit2()$limits[ex,2] }
# if(pl == -Inf & myfit2()$limits[ex,1] == -Inf){
# pl <- qnorm(0.001, myfit2()$Normal[ex,1], myfit2()$Normal[ex,2])}
# if(pu == Inf & myfit2()$limits[ex,2] == Inf){
# pu <- qnorm(0.999, myfit2()$Normal[ex,1], myfit2()$Normal[ex,2])}
# p <- c(0, myfit2()$probs[ex,], 1)
# x <- c(pl, myfit2()$vals[ex,], pu)
# values <- qhist(c(0.05, 0.95), x, p)
# }
#
# if(as.numeric(input$radio1)>1){
# temp<-feedback(myfit1(), quantiles=c(0.05, 0.95), ex=1)
# values=temp$fitted.quantiles[,index]
# }
# data.frame(quantiles=c(0.05, 0.95), values=values)
#
# })
#
#
# output$valuesPDF1 <- renderTable({
# quantileValues1()
# })
#
# output$valuesPDF2 <- renderTable({
# quantileValues2()
# })
# Helper function KMextrapolate ----
# Helper function survivalTable ----
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
}
# Helper function check quartiles can be plotted ----
checkPlot <- function(myvals){
if(anyNA(myvals)){
return(FALSE)
}
if(any(diff(myvals)<=0)){
return(FALSE)
}
return(TRUE)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.