# app.R
#
# Input: texts.txt
#
# by: Paul Wellingerhof, 24th August 2018
library(shiny)
library(shinydashboard)
# options(width=10000)
##################
### load files ###
##################
# this file already includes a bit of HTML-formatting to improve this R-code's
# readability
txt <- strsplit(paste(readLines("texts.txt"), collapse=" "),
"<!-- [[:alnum:]]* -->")[[1]][-1]
##########
### ui ###
##########
### ui-input
# --- Demonstrations --- #
# Publication Bias Test
inSetSize <- sliderInput("nrepl", "How big is the set of studies?",
1, 10, 8, 1)
inPropSig <- sliderInput("propsig",
"How many of the 8 studies are significant?", 0, 8, 7, 1)
inBetaMode <- sliderInput("pmax", "Mode of the Beta distribution",
0, 1, .5, .05)
inShapeSum <- sliderInput("shapesum", paste("Sum of the two shape parameters",
"of the Beta distribution (dispersion)"), 10, 200, 10, 1)
# Identify Publication Bias
shoPow <- checkboxInput("shopow", "Show power values", FALSE)
setSelOne <- selectInput("setsel1", label="", choices=c("", "1. Set",
"2. Set", "3. Set"))
setSelTwo <- selectInput("setsel2", label="", choices=c("", "1. Set",
"2. Set", "3. Set"))
setSelThree <- selectInput("setsel3", label="", choices=c("", "1. Set",
"2. Set", "3. Set"))
resetBut <- actionButton("resetbut", "Reset")
submitBut <- actionButton("submitbut", "Submit!",
style="color:#ffffff;background-color:#337ab7;border-color:#2e6da4")
# --- Data Input --- #
inPropSigDI <- sliderInput("propsigDI",
"How many of the studies were significant?", 0, 10, 5, 1)
pbtBut <- actionButton("pbtbut", "Publication Bias Test!",
style="color:#fff; background-color:#ea8100;border-color:#ce7100")
### ui-page
ui <- dashboardPage(
# HEADER
dashboardHeader(title="Excess of Success"),
# SIDEBAR
dashboardSidebar(
sidebarMenu(id="tabs",
# --- Welcome --- #
menuItem("Welcome!", tabName="welcome", icon=icon("info-circle")),
# --- Introduction --- #
menuItem("Introduction", tabName="intro", icon=icon("align-left")),
# --- Demonstrations --- #
menuItem("Demonstrations", tabName="demons", icon=icon("child"),
# Publication Bias Test
menuSubItem("Publication Bias Test", tabName="powsums"),
# Identify Publication Bias
menuSubItem("Identify Publication Bias", tabName="sets")
),
# --- Data Input --- #
menuItem("Data Input", tabName="userin", icon=icon("file"))
)
),
# BODY
dashboardBody(
includeScript("../../../Matomo-tquant.js"),
# make tick green and cross red (Identify Publication Bias)
tags$style(".fa-check {color:#15ff00}"),
tags$style(".fa-times {color:#ff0000}"),
tabItems(
# --- Welcome --- #
tabItem(tabName="welcome",
fluidRow(
box(
title="Welcome!", width=4, htmlOutput("welcometxt")
)
)
),
# --- Introduction --- #
tabItem(tabName="intro",
fluidRow(
column(width=4,
box(
width=NULL, htmlOutput("intro1")
),
box(
width=NULL, htmlOutput("intro2")
)
),
column(width=4,
box(
width=NULL, htmlOutput("intro3")
)
),
column(width=4,
box(
width=NULL, htmlOutput("intro4")
),
box(
title="References", width=NULL, htmlOutput("refs")
)
)
)
),
# --- Demonstrations --- #
# Publication Bias Test #
tabItem(tabName="powsums",
fluidRow(
column(width=4,
box(width=NULL,
htmlOutput("PBT1")
),
box(title="Study set attributes", width=NULL,
inSetSize, inPropSig
),
box(title="Power sample generation", width=NULL,
inBetaMode, inShapeSum, plotOutput("betadist")
)
),
column(width=8,
box(title="Power sample", width=NULL,
column(width=6, htmlOutput("info")),
column(width=6, plotOutput("sampledist"))
),
box(title=HTML("Computing the <i>p</i>-value"), width=NULL,
column(width=9, htmlOutput("combmat")),
column(width=3, htmlOutput("pval"))
)
)
)
),
# Identify Publication Bias
tabItem(tabName="sets",
fluidRow(
column(width=4,
box(width=NULL,
htmlOutput("setsintro")
)
),
column(width=8,
box(width=NULL,
shoPow,
column(width=4, tableOutput("settab1"), setSelOne,
uiOutput("rslt1"), htmlOutput("sol1"), resetBut),
column(width=4, tableOutput("settab2"), setSelTwo,
uiOutput("rslt2"), htmlOutput("sol2")),
column(width=4, tableOutput("settab3"), setSelThree,
uiOutput("rslt3"), htmlOutput("sol3"), submitBut)
)
)
)
),
# --- Data Input --- #
tabItem(tabName="userin",
fluidRow(
column(width=4,
box(width=NULL,
htmlOutput("DIintro"), fileInput("powerdat", ""),
inPropSigDI, pbtBut
)
),
column(width=8,
box(width=NULL,
column(width=6, htmlOutput("combmatDI")),
column(width=6, htmlOutput("pvalDI"),
plotOutput("sampledistDI"), htmlOutput("powsumtext"))
)
)
)
)
)
)
)
##############
### server ###
##############
server <- function(input, output, session){
# --------- Welcome --------- #
# text: WELCOME
output$welcometxt <- renderUI({
HTML(txt[1])
})
# --------- Introduction --------- #
# text: INTRO1
output$intro1 <- renderUI({
HTML(txt[2])
})
# text: INTRO2
output$intro2 <- renderUI({
HTML(txt[3])
})
# text: INTRO3
output$intro3 <- renderUI({
HTML(txt[4])
})
#table: INTRO4
output$intro4 <- renderUI({
HTML(txt[5])
})
# text: REFERENCES
output$refs <- renderUI({
HTML(txt[6])
})
# --------- Demonstrations --------- #
########## Publication Bias Test ###########
# output: SMALL INTRO TEXT FOR PBT
output$PBT1 <- renderUI({
HTML(txt[7])
})
# observeEvent: DEPENDENT PROPORTION SIGNIFICANT SLIDER
observeEvent(input$nrepl, {
if(input$propsig > input$nrepl) newPropsig <- 1
else newPropsig <- input$propsig
updateSliderInput(session, inputId="propsig", max=input$nrepl,
value=newPropsig, label=paste("How many of the", input$nrepl,
"studies are significant?"), step=1)
})
# reactive: POWERSAMPLE & COMBINATION-MATRIX
powlist <- reactive({
# sampling
vec <- round(rbeta(input$nrepl,
# alpha-parameter
shape1=input$pmax*(input$shapesum - 2) + 1,
# beta-parameter
shape2=input$shapesum - (input$pmax*(input$shapesum - 2) + 1)),
3
)
# distribution
ps <- sapply(0:length(vec),
function(i) dbinom(i, length(vec), mean(vec)))
# combination-matrix
pow <- vec
beta <- 1 - pow
allCombs <- as.matrix(expand.grid(rep(list(0:1), input$nrepl)))
combsPow <- allCombs[sapply(1:nrow(allCombs),
function(i) sum(allCombs[i, ])) >= input$propsig, , drop=FALSE]
combsBeta <- 1 - combsPow
if(input$nrepl == 1 && input$propsig == 0){
mPB <- t(t(sapply(1:nrow(combsPow), function(i) combsPow[i, ]*pow)))
mB <- t(t(sapply(1:nrow(combsBeta), function(i) combsBeta[i, ]*beta)))
} else {
mPB <- t(sapply(1:nrow(combsPow), function(i) combsPow[i, ]*pow))
mB <- t(sapply(1:nrow(combsBeta), function(i) combsBeta[i, ]*beta))
}
mPB[mPB == 0] <- mB[mB != 0]
mPB <- cbind(mPB, apply(mPB, 1, prod))
attr(mPB, "dimnames") <- list(Combination=1:nrow(mPB),
Study=c(1:input$nrepl, "Product"))
# p-value
p <- sum(mPB[, ncol(mPB)])
# list output
list(vec=vec, ps=ps, max=which.max(ps) - 1, mPB=mPB, p=p)
})
# output: BETA DISTRIBUTION
output$betadist <- renderPlot({
par(mai=c(.8, .8, .5, .2), mgp=c(2, .7, 0))
curve(dbeta(x,
input$pmax*(input$shapesum - 2) + 1,
input$shapesum - (input$pmax*(input$shapesum - 2) + 1)),
xlab="Power value", ylab="Density",
main="Beta Probability-Density-Function")
})
# output: SAMPLE DISTRIBUTION
output$sampledist <- renderPlot({
par(mai=c(.8, .8, .5, .2), mgp=c(2, .7, 0))
plot(0:input$nrepl, powlist()$ps, type="h", col="blue", ylim=c(0, 1),
axes=FALSE, xlab="Number of successful replications (n)",
ylab="P(X = n)", xlim=c(-.25, length(powlist()$vec) + .25),
main="Sample Probability Distribution")
if (input$nrepl <= 10) axis(1, at=0:input$nrepl)
else axis(1)
axis(2)
abline(v=sum(powlist()$vec), col="red", lty=2)
legend("topleft", legend="Sum of Power values", col="red", lty=2, bty="n")
graphics::box()
})
# output: TEXT (powersample & its statistics)
output$info <- renderUI({
# labels
labels <- c(
"Power vector",
"Sum of Power values ~ n, where P(X = n) maximal",
"Sample Mean, Sample Standard Deviation"
)
labels <- gsub("$", "</b></br>", gsub("^", "<b>", labels))
# values
vals <- c(
gsub("(.{55,}?)\\s", "\\1<br/>",
paste(format(powlist()$vec), collapse=" ")),
paste(sum(powlist()$vec), "~", powlist()$max),
paste0("<i>M</i> = ", round(mean(powlist()$vec), 3), ", <i>SD</i> = ",
round(sd(powlist()$vec), 3))
)
vals <- gsub("$", "</pre>",
gsub("^", "<pre>", vals))
HTML(paste(labels, vals))
})
# output: COMBINATION-MATRIX
output$combmat <- renderUI({
HTML(paste("<pre style='height:20em'>",
paste(gsub("$", "</br>", capture.output(powlist()$mPB)), collapse=""),
"</pre>", collapse="")
)
})
# output: P-VALUE
output$pval <- renderUI({
withMathJax(
sprintf("The result of summing up the product column is $$p = %.03f.$$",
powlist()$p
)
)
})
########## Identify Publication Bias ###########
# output: INTRO TEXT FOR SETS
output$setsintro <- renderUI({
HTML(txt[8])
})
# reactive: TABLE (set 1)
tabs <- reactive({
input$resetbut
isolate({
# set 1
set1 <- data.frame(study=1:20, n=10, g.est=0, alpha=.05, t=0, p=1,
sig=0)
for(s in 1:nrow(set1)){
samp1 <- rnorm(set1$n[s])
samp2 <- rnorm(set1$n[s])
while(set1$sig[s] == 0){
pval <- t.test(samp1, samp2, var.equal=TRUE)$p.value
if (pval < set1$alpha[s]){
set1$g.est[s] <- -diff(t.test(samp1, samp2,
var.equal=TRUE)$estimate)
set1$t[s] <- t.test(samp1, samp2, var.equal=TRUE)$statistic
set1$p[s] <- pval
set1$sig[s] <- 1
} else-if (set1$n[s] == 30) {
set1$g.est[s] <- -diff(t.test(samp1, samp2,
var.equal=TRUE)$estimate)
set1$t[s] <- t.test(samp1, samp2,var.equal=TRUE)$statistic
set1$p[s] <- pval
break
} else {
set1$n[s] <- set1$n[s] + 1
samp1[set1$n[s]] <- rnorm(1)
samp2[set1$n[s]] <- rnorm(1)
}
}
}
idx1 <- sort(as.numeric(row.names(set1[order(set1$p), ][1:5, ])))
set1 <- set1Pow <- round(set1[idx1, c(2, 3, 5:7)], 3)
set1Pow$power <- sapply(1:nrow(set1),
function(i) round(power.t.test(set1$n[i], set1$g.est[i])$power, 3))
# set 2
set2 <- data.frame(study=1:100, n=sample(10:30, 100, TRUE), g.est=.1,
alpha=.05, t=0, p=1, sig=0)
for(s in 1:nrow(set2)){
tres <- t.test(rnorm(set2$n[s]), rnorm(set2$n[s], .1), var.equal=TRUE)
set2$g.est[s] <- -diff(tres$estimate)
set2$t[s] <- tres$statistic
set2$p[s] <- tres$p.value
if(set2$p[s] < set2$alpha[s]) set2$sig <- 1
}
idx2 <- sort(as.numeric(row.names(set2[order(set2$p), ][1:5, ])))
set2 <- set2Pow <- round(set2[idx2, c(2, 3, 5:7)], 3)
set2Pow$power <- sapply(1:nrow(set2),
function(i) round(power.t.test(set2$n[i], set2$g.est[i])$power, 3))
# set 3
set3 <- data.frame(study=1:5, n=sample(10:30, 5, TRUE), g.est=.8,
alpha=.05, t=0, p=1, sig=0)
for(s in 1:nrow(set3)){
tres <- t.test(rnorm(set3$n[s]), rnorm(set3$n[s], .8), var.equal=TRUE)
set3$g.est[s] <- -diff(tres$estimate)
set3$t[s] <- tres$statistic
set3$p[s] <- tres$p.value
if(set3$p[s] < set3$alpha[s]) set3$sig <- 1
}
set3 <- set3Pow <- round(set3[, c(2, 3, 5:7)], 3)
set3Pow$power <- sapply(1:nrow(set3),
function(i) round(power.t.test(set3$n[i], set3$g.est[i])$power, 3))
set1Pow$n <- as.integer(set1Pow$n)
set2Pow$n <- as.integer(set2Pow$n)
set3Pow$n <- as.integer(set3Pow$n)
list(tab1=set1Pow, tab2=set2Pow, tab3=set3Pow)
})
})
setSamp <- reactive({
input$resetbut
isolate({
sample(1:3)
})
})
# output: SET1
output$settab1 <- renderTable(digits=3, {
if(input$shopow) tabs()[[setSamp()[1]]][, -5]
else tabs()[[setSamp()[1]]][, -c(5, 6)]
})
# output: SET2
output$settab2 <- renderTable(digits=3, {
if(input$shopow) tabs()[[setSamp()[2]]][, -5]
else tabs()[[setSamp()[2]]][, -c(5, 6)]
})
# output: SET3
output$settab3 <- renderTable(digits=3, {
if(input$shopow) tabs()[[setSamp()[3]]][, -5]
else tabs()[[setSamp()[3]]][, -c(5, 6)]
})
# observeEvent: RESET SELECTIVE INPUT
observeEvent(input$resetbut, {
updateSelectInput(session, "setsel1", selected="")
updateSelectInput(session, "setsel2", selected="")
updateSelectInput(session, "setsel3", selected="")
output$rslt1 <- renderUI({})
output$rslt2 <- renderUI({})
output$rslt3 <- renderUI({})
output$sol1 <- renderUI({})
output$sol2 <- renderUI({})
output$sol3 <- renderUI({})
})
# reactive: RESULTS
observeEvent(input$submitbut, {
# output: ICON1
output$rslt1 <- renderUI({
isolate(
if(substr(input$setsel1, 1, 1) == as.character(setSamp()[1])){
icon("check")
} else {
icon("times")
}
)
})
# output: ICON2
output$rslt2 <- renderUI({
isolate(
if(substr(input$setsel2, 1, 1) == as.character(setSamp()[2])){
icon("check")
} else {
icon("times")
}
)
})
# output: ICON1
output$rslt3 <- renderUI({
isolate(
if(substr(input$setsel3, 1, 1) == as.character(setSamp()[3])){
icon("check")
} else {
icon("times")
}
)
})
# output: SOLULTION1
output$sol1 <- renderUI({
HTML(paste0("Solution: This is the ", setSamp()[1], ". set."))
})
# output: SOLUTION2
output$sol2 <- renderUI({
HTML(paste0("Solution: This is the ", setSamp()[2], ". set."))
})
# output: SOLUTION3
output$sol3 <- renderUI({
HTML(paste0("Solution: This is the ", setSamp()[3], ". set."))
})
})
# --------- Data Input --------- #
# output: DI-INTRO
output$DIintro <- renderUI({
HTML(txt[9])
})
# reactive: POWER DATAFRAME
powdat <- reactive({
file <- input$powerdat
if(is.null(file)) return()
dat <- read.table(file$datapath, header=TRUE)
# if(ncol(dat) > 1 | class(dat[, 1]) != "numeric") return()
dat
})
# observeEvent: PROPORTION SIGNIFICANT DI
observeEvent(input$powerdat, {
updateSliderInput(session, "propsigDI", max=nrow(powdat()))
})
# reactive: COMBINATION MATRIX, P-VALUE
pbtDI <- reactive({
pow <- powdat()[, 1]
beta <- 1 - pow
allCombs <- as.matrix(expand.grid(rep(list(0:1), length(pow))))
combsPow <- allCombs[sapply(1:nrow(allCombs),
function(i) sum(allCombs[i, ])) >= input$propsigDI, , drop=FALSE]
combsBeta <- 1 - combsPow
if(length(pow) == 1 && input$propsigDI == 0){
mPB <- t(t(sapply(1:nrow(combsPow), function(i) combsPow[i, ]*pow)))
mB <- t(t(sapply(1:nrow(combsBeta), function(i) combsBeta[i, ]*beta)))
} else {
mPB <- t(sapply(1:nrow(combsPow), function(i) combsPow[i, ]*pow))
mB <- t(sapply(1:nrow(combsBeta), function(i) combsBeta[i, ]*beta))
}
mPB[mPB == 0] <- mB[mB != 0]
mPB <- cbind(mPB, apply(mPB, 1, prod))
attr(mPB, "dimnames") <- list(Combination=1:nrow(mPB),
Study=c(1:length(pow), "Product"))
# distribution
ps <- dbinom(0:length(pow), length(pow), mean(pow))
# p-value
p <- sum(mPB[, ncol(mPB)])
list(mPB=mPB, p=p, ps=ps, pow=pow, max=which.max(ps) - 1)
})
# observeEvent: ALL DI PBT OUTPUTS
observeEvent(input$pbtbut, {
# output: COMBINATION MATRIX DI
output$combmatDI <- renderUI({
HTML(paste("<pre style='height:40em'>",
paste(gsub("$", "</br>", capture.output(pbtDI()$mPB)), collapse=""),
"</pre>", collapse="")
)
})
# output: P-VALUE DI
output$pvalDI <- renderUI({
withMathJax(
sprintf(
"The result of summing up the product column is $$p = %.03f.$$",
pbtDI()$p
)
)
})
# output: SAMPLE DISTRIBUTION
output$sampledistDI <- renderPlot({
par(mai=c(.8, .8, .5, .2), mgp=c(2, .7, 0))
plot(0:length(pbtDI()$pow), pbtDI()$ps, type="h", col="blue",
ylim=c(0, 1), axes=FALSE,
xlab="Number of successful replications (n)",
ylab="P(X = n)", xlim=c(-.25, length(pbtDI()$pow) + .25),
main="Sample Probability Distribution")
if (length(pbtDI()$pow) <= 10) axis(1, at=0:length(pbtDI()$pow))
else axis(1)
axis(2)
abline(v=sum(pbtDI()$pow), col="red", lty=2)
legend("topleft", legend="Sum of Power values", col="red", lty=2,
bty="n")
graphics::box()
})
# output: POWSUM TEXT
output$powsumtext <- renderUI({
HTML(paste(tags$b("Sum of Power values ~ n, where P(X = n) maximal"),
"<pre>", sum(pbtDI()$pow), " ~ ", pbtDI()$max, "</pre>", collapse=""))
})
})
}
################
### shinyApp ###
################
shinyApp(ui, server)