library(shiny)
ui <- fluidPage(
titlePanel("Simulating Statistical Power"),
tabsetPanel(
## INTRODUCTION
tabPanel("Introduction",
includeScript("../../../Matomo-tquant.js"),
sidebarLayout(
sidebarPanel(
HTML("<b>Purpose of the app</b><br>
This shiny app aims to help students to understand what is
statistical power and to estimate it through simulated data.
The app allows the user to perform simulations with commonly
used tests, such as one- and two-sample t test and one-way
ANOVA with and without repeated measures.")
),
mainPanel(
HTML("<br><b>Statistical power</b><br>
Power represents the ability of a statistical test to detect
an effect that is real.
It depends on four main elements:
<ul>
<li><b>the type I error</b> (i.e., the criterion which
determinates how unlikely a positive result must be
for the null hypothesis to be rejected)</li>
<li><b>the effect size</b> (i.e., the magnitude of the
effect of interest in the population)</li>
<li><b>the sample size used to detect the effect</b>
(i.e. which determinates the amount of sampling error inherent
in the test results)</li>
<li><b>the variability</b> (i.e., the standard deviation of
the sample)</li>
</ul>
<br>
<b>Why power matters</b>
<br>Estimating the statistical power (i.e., power analysis)
could be very useful <i>before</i> collecting the data.
<br>The power analysis allows you to decide:
<ul>
<li>how large a sample is needed to enable statistical
judgments that are accurate and reliable</li>
<li>how likely a statistical test will be to detect
effects of a given size in a particular situation.</li>
</ul>
<br>
<b>How to estimate power</b>
<br>Power can be determined analytically, by setting all
parameters that affect it.
Alternatively, it can be estimated by simulating data,
performing statistical tests on them, and calculating the
proportion of p-values that fall below the chosen significance
criterion.")
)
)
),
## WHAT IS POWER
tabPanel("What is power?",
sidebarLayout(
sidebarPanel(
HTML("<br> Change the value of Alpha and see how the power changes!
<br><br>"),
sliderInput("Alpha", "Alpha", .001,.999,0.05),
verbatimTextOutput("indices"),
HTML("<br> Let the red cross be your sample mean.
You want to know if its mean is different from a hypothetical
value.Thus, you perform a statistical test in order to make a
decision on rejecting or not rejecting the null hypothesis
(i.e., no difference).<br><br>")
),
mainPanel(
h3("What is power?"),
HTML("While we are always focusing on <b>Type I error</b>
(i.e., false positive), we frequently ignore that by
decreasing alpha the power becomes lower. The power is
determined by subtracting the <b>Type II error</b> (i.e.
false negative) from 1. It is the probability
to correctly reject H0 when it is false. In other
worlds, the power is the ability of a test to detect an
effect that is real.
<br>"),
plotOutput("powerplot"),
h3("Simulated power"),
HTML("The statistical power could be estimated by simulating data.
It requires the following steps:
<ol>
<li>Specify the effect size, the randomization procedure (e.g.,
normally distributed data), the test statistic, and the
significance level (i.e., alpha) that will be used.</li>
<li>Generate random samples from the distribution specified
in the previous step.</li>
<li>Calculate the test statistic from the simulated data
and determine if the null hypothesis is accepted or
rejected.</li>
<li>Repeat steps 2 and 3 several thousand times and count
the number of times the simulated data lead to a rejection of
the null hypothesis.</li>
</ol>
The power of the test is <b>the proportion of simulated
samples that lead to a rejection of H0</b>.
<br> Note that as the number of simulations increases, the
precision of the power estimate will increase as well.
Also note that if you set the hypothesized parameter equal to
the true parameter (i.e., an effect size of zero), the power
equals alpha."),
br(),
br(),
br()
)
)
),
## ONE SAMPLE T-TEST
tabPanel("One-sample t test",
sidebarLayout(
sidebarPanel(
h4("Sample size"),numericInput("n","", 20, min = 5, max = 200),
h4("Hypothesized mean"),sliderInput("mu1", "", 400, 420, 405),
h4("True mean"),sliderInput("mu2", "", 400, 420, 405),
h4("Standard deviation"),sliderInput("sd", "", 1, 20, 8),
h4("Alpha"),radioButtons("alpha", "", c(".01", ".05", ".10"),
selected = ".05", inline = TRUE),
h4("Number of replicates"),
radioButtons("nrep", "", c("500", "1000", "5000"),
selected = "1000", inline = TRUE),
actionButton("simulate", "Simulate power")
),
mainPanel(
h3("Power estimation in a one-sample t test"),
HTML("By setting the mean and the standard deviation of the
simulated samples, and the mean under the null hypotesis and
the significance criterion, you can estimate the power for a
one-sample t test.
<br><br> Here you can think of an experiment, where
participants compare two sounds with different frequencies. One
is given (hypothesized mean), the participant has to adjust
the other (true mean), so that they cannot hear the difference.
The app will simulate data for the participants, and calculate
the power of a one-sample t test, to compare the hypothesized
mean and the true mean."),
br(),
br(),
plotOutput("thistplot"),
br(),
br(),
br(),
plotOutput("tdensityplot"),
br(),
br()
)
), value = "One-sample t test"
),
## TWO SAMPLE T-TEST
tabPanel("Two-sample t test",
sidebarLayout(
sidebarPanel(
h4("Sample size"),
numericInput("tn1", "Group 1", 20, min = 5, max = 200),
numericInput("tn2", "Group 2", 20, min = 5, max = 200),
h4("Mean"),
sliderInput("tmu1", "Group 1", 400, 420, c(405)),
sliderInput("tmu2", "Group 2", 400, 420, c(405)),
h4("Standard deviation"),
sliderInput("tsd1", "Group 1", 1, 20, 8),
sliderInput("tsd2", "Group 2", 1, 20, 8),
h4("Alpha"),
(radioButtons("talpha", "", c(".01", ".05", ".10"),
selected = ".05", inline = TRUE)),
h4("Number of replicates"),
radioButtons("nrep","", c("500", "1000", "5000"),
selected = "1000", inline = TRUE),
actionButton("tsimulate", "Simulate power")
),
mainPanel(
h3("Power estimation in a two-sample t test"),
HTML("Here you want to estimate the power of your statistical
test in detecting true differences among the means of two
samples.
You need to set the size, the mean and the standard deviation
of both groups.
<br><br>
Now think of an experiment, where you compare two
different groups of participants. Maybe you have a group of
hearing impaired participants and one group of normal
participants. In your experiment all participants have to
adjust the frequency of a sound."),
br(),
br(),
plotOutput("t2histplot"),
br(),
br(),
br(),
plotOutput("t2densityplot"),
br(),
br()
)
)
),
## ONE WAY ANOVA
tabPanel("One-way ANOVA",
sidebarLayout(
sidebarPanel(
h4("Sample size"),
numericInput("an1", "Group 1", 20, min = 5, max = 200),
numericInput("an2", "Group 2", 20, min = 5, max = 200),
numericInput("an3", "Group 3", 20, min = 5, max = 200),
h4("Mean"),
sliderInput("amu1", "Group 1", 400, 420, c(405)),
sliderInput("amu2", "Group 2", 400, 420, c(405)),
sliderInput("amu3", "Group 3", 400, 420, c(405)),
h4("Standard deviation"),
sliderInput("s1", "Group 1", 1, 20, 8),
sliderInput("s2", "Group 2", 1, 20, 8),
sliderInput("s3", "Group 3", 1, 20, 8),
h4("Alpha"),
radioButtons("aalpha", "", c(".01", ".05", ".10"),
selected = ".05", inline = TRUE),
h4("Number of replicates"),
radioButtons("anrep","", c("500", "1000", "5000"),
selected = "1000", inline = TRUE),
actionButton("asimulate", "Simulate power")
),
mainPanel(
h3("Power estimation in a one-way analysis of variance"),
HTML("Finally, let's complicate the situation a bit! Now there
are three groups and you are performing a one-way ANOVA,
by testing the effect of a factor with three levels (e.g.,
treatment, placebo, nothing).
<br><br>
Think of an experiment, where you compare three
different groups of participants. Maybe you examine hearing
impaired participants, and one group gets a hearing aid, the
other group gets a cochlea implant, and the third groups does
not receive any support. All groups have to adjust the
frequency of the sound."),
br(),
br(),
plotOutput("ahistplot"),
br(),
br(),
br(),
plotOutput("adensityplot"),
br(),
br()
)
)
),
## ONE WAY REPEATED MEASURES ANOVA
tabPanel("One-way repeated measures ANOVA",
sidebarLayout(
sidebarPanel(
h4("Sample size"),
numericInput("rman", "", 20, min = 5, max = 200),
h4("Mean"),
sliderInput("rmamu1", "At time 1", 400, 420, c(405)),
sliderInput("rmamu2", "At time 2", 400, 420, c(405)),
sliderInput("rmamu3", "At time 3", 400, 420, c(405)),
h4("Standard deviation"),
sliderInput("rmaspers", "between subjects", 1, 20, 8),
sliderInput("rmaserror", "within subjects", 0, 10, 1),
h4("Alpha"),
radioButtons("rmaalpha", "", c(".01", ".05", ".10"),
selected = ".05", inline = TRUE),
h4("Number of replicates"),
radioButtons("rmanrep","", c("500", "1000", "5000"),
selected = "1000", inline = TRUE),
actionButton("rmasimulate", "Simulate power")
),
mainPanel(
h3("Power estimation in a one-way analysis of variance with repeated
measures"),
HTML("Let's think of a new situation. You observe just one group at
three different times (e.g., before (pre), during, and after
(post) a treatment).
As you want to test the effect of the time, you are
performing a one-way ANOVA with repeated measures.
<br><br>
Now think of an experiment, where you compare one
group of participants at three different times.
Maybe you have a training for hearing the frequency of a sound,
and you test your participants before, during and after the
training."),
br(),
br(),
plotOutput("rmahistplot", width = "100%"),
br(),
br(),
br(),
plotOutput("rmadensityplot"),
br(),
br()
)
)
)
) # end of TabSetPanel
) # end of UI
server <- function(input, output){
## INTRODUCTION
output$powerRangers_image <- renderImage({
outfile <- tempfile(fileext='powerRangers.png')
png(outfile, width=200, height=200)
dev.off()
list(src = outfile,
alt = "This is alternate text")
}, deleteFile = TRUE)
## VISUALIZATION
crit <- reactive({
qnorm(1-as.numeric(input$Alpha)) # setting the critical value
})
output$powerplot <- renderPlot({
par(mai=c(1,0,.5,0),mgp=c(0.5,0,0))
# distribution under H0
curve(dnorm(x,0,1),xlim=c(-3,6),
bty="n",xaxt="n",yaxt="n",lwd=2,
ylab='',xlab='')
# distribution under H1
curve(dnorm(x,3,1),lwd=2,add=TRUE)
# type II error area
polygon(c(0,seq(0,6,0.01),6), #
c(0,dnorm(mean=3,seq(0,6,0.01)),0),
col="#0080FF")
# power area
polygon(c(crit(),seq(crit(),10,0.01),10),
c(0,dnorm(mean=3,seq(crit(),10,0.01)),0),
col="#3ADF00")
# type I error area
polygon(c(crit(),seq(crit(),4,0.01),4),
c(0,dnorm(seq(crit(),4,0.01)),0),
col="#FF8000")
# rejecting line and sample point
abline(v=crit(),col="red",lwd=3)
text(x=crit()+1,y=0.1,labels="Rejecting H0",cex=1)
text(x=crit()-1.3,y=0.1,labels="Not rejecting H0",cex=1)
points(2,0,col="darkred",pch=4,cex=3,lwd=3)
# legend
legend(-3, 0.4,
legend=c(expression("Type I error"~(alpha)),
expression("Type II error"~(beta)),
expression("Power"~(1 - beta))),
fill=c("#FF8000","#0080FF","#3ADF00"),
cex=1.2, bty = "n")
})
output$indices <- renderText(paste("Alpha =",input$Alpha,
"\nBeta =", round(pnorm(crit(),mean=3),2),
"\nPower =",round(1-pnorm(crit(),mean=3),2),
"\nConclusion:\n",
if(crit()<2){"Rejecting H0"}
else{"Not rejecting H0"}))
## ONE SAMPLE T-TEST
# data
pval <- eventReactive(input$simulate, {
n1 <- input$n
hmu <- input$mu1; tmu <- input$mu2
s1 <- input$sd
nrep <- input$nrep
tp <- data.frame()
withProgress(message = "Simulating data", value = 0, {
for(i in 1:nrep){
y <- rnorm(n1, tmu, s1)
tp[i, 1] <- t.test(y, mu = hmu)$p.value
tp[i, 2] <- mean(y)
incProgress(1/as.numeric(nrep))
}
})
return(tp)
})
power <- reactive({
mean(pval()[,1] < as.numeric(input$alpha))
})
# distribution of p-values: plot
tpa <- eventReactive(input$simulate, {
par(mai = c(.6, .7, 1, .1), mgp = c(2, .7, 0))
hist(pval()[,1], col = "lightblue", border = "white", xlim = 0:1,
xlab = "p-value", main = "",
breaks = seq(0, 1, 0.05))
abline(v = input$alpha, col = "darkblue")
title("Distribution of p-values", adj = 0)
mtext(at = c(as.numeric(input$alpha) + 0.03,
par("usr")[4]- 0.1 * par("usr")[4]),
text = bquote(alpha == .(input$alpha)), cex = 1.2)
mtext(bquote(paste("Power =", " ",
frac("number of rejected hypothesis",
"total number of simulations"), " = ",
frac(.(sum(pval()[,1] < as.numeric(input$alpha))),
.(length(pval()[,1]))), " = ",
.(round(power(),3)))),
side = 3, line = 1, adj = 0.95)
box(bty = "l")
})
output$thistplot <- renderPlot({
tpa()
})
# Density of the mean values of simulated data: plot
tdp <- eventReactive(input$simulate, {
par(mai = c(.6, .7, 1, .1), mgp = c(2, .7, 0))
ymaxlab <- 2 / input$sd
plot(density(pval()[,2], adjust = 2), col = "blue", xlim = c(390, 430),
ylim = c(0, ymaxlab),
main = "",
ylab = "Density", lty = "dotted", bty = "n",
xlab = "Mean values of simulated data")
title("Distribution of the mean values of simulated data", adj = 0)
abline(v = input$mu1, col = "darkblue")
mtext(at = c(input$mu1 + 2, par("usr")[4] - 0.1 * par("usr")[4]),
text = bquote(mu[hypothesized] == .(input$mu1)),
cex = 1.2)
box(bty = "l")
})
output$tdensityplot <- renderPlot({
tdp()
})
## TWO SAMPLE T-TEST
# data
t2pval <- eventReactive(input$tsimulate, {
t2n1 <- input$tn1; t2n2 <- input$tn2
t2mu1 <- input$tmu1; t2mu2 <- input$tmu2
t2s1 <- input$tsd1; t2s2 <- input$tsd2
t2nrep <- input$nrep
t2p <- data.frame()
withProgress(message = "Simulating data", value = 0, {
for(i in 1:t2nrep){
mdata <- data.frame(y = c(rnorm(t2n1, t2mu1, t2s1),
rnorm(t2n2, t2mu2, t2s2)),
group = rep(c("one","two"),
times=c(t2n1, t2n2)))
t2p[i, 1] <- t.test(mdata[,1][mdata$group == "one"],
mdata[,1][mdata$group == "two"])$p.value
t2p[i, 2] <- mean(mdata[,1][mdata$group == "one"])
t2p[i, 3] <- mean(mdata[,1][mdata$group == "two"])
incProgress(1/as.numeric(t2nrep))
}
})
return(t2p)
})
t2power <- reactive({
mean(t2pval()[,1] < as.numeric(input$talpha))
})
# distribution of p-values: plot
t2pa <- eventReactive(input$tsimulate, {
par(mai = c(.6, .7, 1, .1), mgp = c(2, .7, 0))
hist(t2pval()[,1], col = "lightblue", border = "white", xlim = 0:1,
xlab = "p-value", main = "",
breaks = seq(0, 1, 0.05))
abline(v = input$talpha, col = "darkblue")
title("Distribution of p-values", adj = 0)
mtext(at = c(as.numeric(input$talpha) + 0.03,
par("usr")[4]- 0.1 * par("usr")[4]),
text = bquote(alpha == .(input$talpha)), cex = 1.2)
mtext(bquote(paste("Power =", " ",
frac("number of rejected hypothesis",
"total number of simulations"), " = ",
frac(.(sum(t2pval()[,1] < as.numeric(input$talpha))),
.(length(t2pval()[,1]))), " = ",
.(round(t2power(),3)))),
side = 3, line = 1, adj = 0.95)
box(bty = "l")
})
output$t2histplot <- renderPlot({
t2pa()
})
# Density of the mean values of simulated data: plot
t2dp <- eventReactive(input$tsimulate, {
par(mai = c(.6, .7, 1, .1), mgp = c(2, .7, 0))
ymaxlab <- 2 * max(c(1/input$tsd1, 1/input$tsd2))
gmu <- mean(c(input$tmu1, input$tmu2))
plot(density(t2pval()[,2], adjust = 2), col = "blue", xlim = c(390, 430),
ylim = c(0, ymaxlab),
main = "",
ylab = "Density", lty = "dotted", bty = "n",
xlab = "Mean values of simulated data")
title("Distribution of the mean values of simulated data", adj = 0)
lines(density(t2pval()[,3], adjust = 2), col = "red", lty = "dashed")
legend("topright", c("Group 1", "Group 2"),
col = c("blue", "red"),
lty = c("dotted", "dashed"), lwd = 2)
abline(v = gmu, col = "darkblue")
mtext(at = c(gmu + 2, par("usr")[4] - 0.1 * par("usr")[4]),
text = bquote(mu[total] == .(round(gmu,2))),
cex = 1.2)
box(bty = "l")
})
output$t2densityplot <- renderPlot({
t2dp()
})
## ONE WAY ANOVA
# data
apval <- eventReactive(input$asimulate, {
an1 <- input$an1; an2 <- input$an2; an3 <- input$an3
amu1 <- input$amu1; amu2 <- input$amu2; amu3 <- input$amu3
s1 <- input$s1; s2 <- input$s2; s3 <- input$s3
anrep <- input$anrep
ap <- data.frame()
withProgress(message = "Simulating data", value = 0, {
for(i in 1:anrep){
mdata <- data.frame(y = c(rnorm(an1, amu1, s1),
rnorm(an2, amu2, s2),
rnorm(an3, amu3, s3)),
group = rep(c("one","two","three"),
times=c(an1, an2, an3)))
ap[i, 1] <- summary(aov(y ~ group, data = mdata))[[1]]$"Pr(>F)"[1]
ap[i, 2] <- mean(mdata[,1][mdata$group == "one"])
ap[i, 3] <- mean(mdata[,1][mdata$group == "two"])
ap[i, 4] <- mean(mdata[,1][mdata$group == "three"])
incProgress(1/as.numeric(anrep))
}
})
return(ap)
})
apower <- reactive({
mean(apval()[,1] < as.numeric(input$aalpha))
})
# distribution of p-values: plot
hpa <- eventReactive(input$asimulate, {
par(mai = c(.6, .7, 1, .1), mgp = c(2, .7, 0))
hist(apval()[,1], col = "lightblue", border = "white", xlim = 0:1,
xlab = "p-value", main = "",
breaks = seq(0, 1, 0.05))
title("Distribution of p-values", adj = 0)
abline(v = input$aalpha, col = "darkblue")
mtext(at = c(as.numeric(input$aalpha) + 0.03,
par("usr")[4]- 0.1 * par("usr")[4]),
text = bquote(alpha == .(input$aalpha)), cex = 1.2)
mtext(bquote(paste("Power =", " ",
frac("number of rejected hypothesis",
"total number of simulations"), " = ",
frac(.(sum(apval()[,1] < as.numeric(input$aalpha))),
.(length(apval()[,1]))), " = ",
.(round(apower(),3)))),
side = 3, line = 1, adj = 0.95)
box(bty = "l")
})
output$ahistplot <- renderPlot({
hpa()
})
# Density of the mean values of simulated data: plot
adp <- eventReactive(input$asimulate, {
par(mai = c(.6, .7, 1, .1), mgp = c(2, .7, 0))
ymaxlab <- 2 * max(c(1/input$s1, 1/input$s2, 1/input$s3))
gmu <- mean(c(input$amu1, input$amu2, input$amu3))
plot(density(apval()[,2], adjust = 2), col = "blue", xlim = c(390, 430),
ylim = c(0, ymaxlab),
main = "",
ylab = "Density", lty = "dotted", bty = "n",
xlab = "Mean values of simulated data")
title("Distribution of the mean values of simulated data", adj = 0)
lines(density(apval()[,3], adjust = 2), col = "red", lty = "dashed")
lines(density(apval()[,4], adjust = 2), col = "green")
legend("topright", c("Group 1", "Group 2", "Group 3"),
col = c("blue", "red", "green"),
lty = c("dotted", "dashed", "solid"), lwd = 2)
abline(v = gmu, col = "darkblue")
mtext(at = c(gmu + 2, par("usr")[4] - 0.1 * par("usr")[4]),
text = bquote(mu[total] == .(round(gmu,2))),
cex = 1.2)
box(bty = "l")
})
output$adensityplot <- renderPlot({
adp()
})
## ONE WAY ANOVA WITH REPEATED MEASURES
#data
rmapval <- eventReactive(input$rmasimulate, {
rmanrep <- input$rmanrep
rman <- input$rman
rmamu1 <- input$rmamu1; rmamu2 <- input$rmamu2; rmamu3 <- input$rmamu3
rmaspers <- input$rmaspers; rmaserror <- input$rmaserror
rmamu <- c(rmamu1, rmamu2, rmamu3)
grandmu <- mean(rmamu)
rmap <- data.frame()
withProgress(message = "Simulating data", value = 0, {
for(i in 1:rmanrep){
dat <- replicate(rman, {
y <- grandmu + rnorm(3, rmamu-grandmu, rmaspers) +
rnorm(3, 0, rmaserror)
})
dat <- as.data.frame(t(dat))
dat$id <- factor(1:rman)
rmdata <- reshape(dat, varying = list(1:3), direction = "long")
colnames(rmdata) <- c("id", "time", "y")
rmdata$time <- as.factor(rmdata$time)
rmap[i, 1] <- summary(aov(y ~ time + Error(id/time),
data = rmdata))[[2]][[1]]$"Pr(>F)"[1]
rmap[i, 2] <- mean(dat[,1])
rmap[i, 3] <- mean(dat[,2])
rmap[i, 4] <- mean(dat[,3])
incProgress(1/as.numeric(rmanrep))
}
})
return(rmap)
})
rmapower <- reactive({
mean(rmapval()[,1] < as.numeric(input$rmaalpha))
})
# distribution of p-values: plot
hp <- eventReactive(input$rmasimulate, {
par(mai = c(.6, .7, 1, .1), mgp = c(2, .7, 0))
hist(rmapval()[,1], col = "lightblue", border = "white", xlim = 0:1,
xlab = "p-value", main = "",
breaks = seq(0, 1, 0.05))
title("Distribution of p-values", adj = 0)
abline(v = input$rmaalpha, col = "darkblue")
mtext(at = c(as.numeric(input$rmaalpha) + 0.03,
par("usr")[4]- 0.1 * par("usr")[4]),
text = bquote(alpha == .(input$rmaalpha)), cex = 1.2)
mtext(bquote(paste("Power =", " ",
frac("number of rejected hypothesis",
"total number of simulations"), " = ",
frac(.(sum(rmapval()[,1] < as.numeric(input$rmaalpha))),
.(length(rmapval()[,1]))), " = ",
.(round(rmapower(),3)))),
side = 3, line = 1, adj = 0.95)
box(bty = "l")
})
output$rmahistplot <- renderPlot({
hp()
})
# Density of the mean values of simulated data: plot
dp <- eventReactive(input$rmasimulate, {
par(mai = c(.6, .7, 1, .1), mgp = c(2, .7, 0))
ymaxlab <- (1 / (input$rmaspers + input$rmaserror)
+ 1 / (input$rmaspers * 2) + 1 / (2 * input$rmaserror + 1))
gmu <- mean(c(input$rmamu1, input$rmamu2, input$rmamu3))
plot(density(rmapval()[,2], adjust = 2), col = "blue",xlim = c(390, 430),
ylim = c(0, ymaxlab),
main = "",
ylab = "Density", lty = "dotted", bty = "n",
xlab = "Mean values of simulated data")
title("Distribution of the mean values of simulated data", adj = 0)
lines(density(rmapval()[,3], adjust = 2), col = "red", lty = "dashed")
lines(density(rmapval()[,4], adjust = 2), col = "green")
legend("topright", c("at time 1", "at time 2", "at time 3"),
col = c("blue", "red", "green"),
lty = c("dotted", "dashed", "solid"), lwd = 2)
abline(v = gmu, col = "darkblue")
mtext(at = c(gmu + 2, par("usr")[4] - 0.1 * par("usr")[4]),
text = bquote(mu[total] == .(round(gmu,2))),
cex = 1.2)
box(bty = "l")
})
output$rmadensityplot <- renderPlot({
dp()
})
}
shinyApp(ui, server)