Nothing
inst/doc/understanding-gs-designs.R
In GroupSeq: Group Sequential Design Probabilities - With Graphical User Interface
## ----knitr-setup, include = FALSE-----------------------------------------------------------------
library(knitr)
library("ggplot2")
knitr::opts_chunk$set(
comment = "#",
prompt = F,
tidy = FALSE,
cache = FALSE,
collapse = T,
echo = FALSE,
dpi = 300,
fig.width = 5, fig.height = 5
)
okabe_palette <- c('orange' = "#E69F00",
'sky blue' = "#56B4E9",
'bluish green' = "#009E73",
'yellow' = "#F0E442",
'blue' = "#0072B2",
'vermillion' = "#D55E00",
'reddish purple' = "#CC79A7")
old <- options(width = 100L, digits = 10)
## -------------------------------------------------------------------------------------------------
fill_colors = c(X1 = okabe_palette[["bluish green"]],
X2 = okabe_palette[["vermillion"]],
"X1+X2" = "#6B7E3A")
data = data.frame(x = c(1:10, 1:10, 1:20),
y = c(rep(2, 10), rep(1.5, 10), rep(1.5, 20)),
Sample = c(rep("X1", 10), rep("X2", 10), rep("X1+X2", 20))
)
## ---- out.width = "35%", message = FALSE, fig.height = 1, fig.width = 2.5, fig.align = "center", fig.cap = "Schematic view of a study sample - each circle represents some measurement."----
dat = data[1:10, ]
ggplot(dat, aes(x = x, y = y)) +
geom_point(size = 10, shape = 21, aes(fill = Sample)) +
scale_fill_manual(values = fill_colors[1]) +
lims(x = c(0.5, 10.5), y = c(2, 2)) +
theme_void() +
theme(legend.position = "bottom")
## ---- echo = FALSE, out.width = "65%", fig.cap = fig.cap, fig.height=3----------------------------
fig.cap = "Standard Normal Density with critial bounds marking 5% significance level."
blue = okabe_palette[["sky blue"]]
yellow = okabe_palette[["yellow"]]
crit = qnorm(1 - 0.025)
probLabel = paste0(pnorm(-crit)*100, "%")
x = seq(-3, 3, by = 0.01)
normalDensPlot =
ggplot(data.frame(x = x), aes(x = x)) +
stat_function(fun = dnorm, geom = "area", fill = blue, xlim = c(-4, -crit)) +
stat_function(fun = dnorm, geom = "area", fill = yellow, xlim = c(-crit, crit)) +
stat_function(fun = dnorm, geom = "area", fill = blue, xlim = c(crit, 4)) +
stat_function(fun = dnorm, color = okabe_palette[["bluish green"]], size = 1.5) +
geom_segment(x = 2.2, xend = 2.7, y = 0.02, yend = 0.07) +
geom_segment(x = -2.2, xend = -2.7, y = 0.02, yend = 0.07) +
annotate("text", label = probLabel, x = 2.9, y = 0.09, size = 3.5) +
annotate("text", label = probLabel, x = -2.9, y = 0.09, size = 3.5) +
lims(x = range(x)) +
labs(x = expression(x[1]),
y = expression(f(x[1]))) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.x = element_line())
normalDensPlot
## ---- echo = TRUE---------------------------------------------------------------------------------
qnorm(0.025)
## ---- out.width = "35%", message = FALSE, fig.height = 1.3, fig.width = 2.5, fig.align = "center"----
dat = data[1:20, ]
ggplot(dat, aes(x = x, y = y)) +
geom_point(size = 10, shape = 21, aes(fill = Sample)) +
scale_fill_manual(values = fill_colors[1:2]) +
lims(x = c(0.5, 10.5), y = c(1.2, 2.2)) +
theme_void() +
theme(legend.position = "bottom")
## ---- fig.width = 8, fig.height = 3, out.width = "85%"--------------------------------------------
normalDensPlotX2 = normalDensPlot +
labs(x = expression(x[2]), y = expression(f(x[2]))) +
stat_function(fun = dnorm, color = okabe_palette[["vermillion"]], size = 1.5)
gridExtra::grid.arrange(normalDensPlot, normalDensPlotX2, nrow = 1)
## ---- echo = TRUE---------------------------------------------------------------------------------
qnorm(0.025/2)
## ---- echo = FALSE, out.width = "65%", fig.cap = fig.cap, fig.height=3----------------------------
fig.cap = "Standard Normal Density with critial bounds marking 5% significance level."
blue = okabe_palette[["sky blue"]]
yellow = okabe_palette[["yellow"]]
crit = qnorm(1 - 0.025/2)
probLabel = paste0(pnorm(-crit)*100, "%")
x = seq(-3, 3, by = 0.01)
normalDensPlotAdj =
ggplot(data.frame(x = x), aes(x = x)) +
stat_function(fun = dnorm, geom = "area", fill = blue, xlim = c(-4, -crit)) +
stat_function(fun = dnorm, geom = "area", fill = yellow, xlim = c(-crit, crit)) +
stat_function(fun = dnorm, geom = "area", fill = blue, xlim = c(crit, 4)) +
stat_function(fun = dnorm, color = okabe_palette[["bluish green"]], size = 1.2) +
geom_segment(x = 2.2+.2, xend = 2.7+.2, y = 0.02, yend = 0.07) +
geom_segment(x = -2.2-.2, xend = -2.7-.2, y = 0.02, yend = 0.07) +
annotate("text", label = probLabel, x = 2.9, y = 0.09, size = 3.5) +
annotate("text", label = probLabel, x = -2.9, y = 0.09, size = 3.5) +
lims(x = range(x)) +
labs(x = expression(x[1]),
y = expression(f(x[1]))) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.x = element_line())
## ---- fig.width = 8, fig.height = 3, out.width = "85%"--------------------------------------------
normalDensPlotX2Adj = normalDensPlotAdj +
labs(x = expression(x[2]), y = expression(f(x[2]))) +
stat_function(fun = dnorm, color = okabe_palette[["vermillion"]], size = 1.5)
gridExtra::grid.arrange(normalDensPlotAdj, normalDensPlotX2Adj, nrow = 1)
## -------------------------------------------------------------------------------------------------
alpha = 0.05
alpha.Bonf = alpha / 2
alpha.Sidak = 1 - (1-alpha)^(1/2)
crit.Bonf = -rep(qnorm(alpha.Bonf / 2), 2)
crit.Sidak = -rep(qnorm(alpha.Sidak / 2), 2)
## ----prepare-2D-normal-density-plot, echo = FALSE, include = FALSE--------------------------------
mu = 0
Mean = c(mu, mu)
covariance = 0
Sigma0 = matrix(c(1, rep(covariance, 2), 1), 2)
width = 6.6
binwidth = 0.2
x1 = x2 = seq(from = mu - width/2, to = mu + width/2, by = binwidth)
grid = expand.grid(x1, x2)
z = mvtnorm::dmvnorm(grid, mean = Mean, sigma = Sigma0)
zmat = matrix(z, ncol = length(x1))
grid.col = expand.grid(x1[-1] - binwidth/2, x2[-1] - binwidth/2)
gc1 = grid.col[, 1]
gc2 = grid.col[, 2]
crit = crit.Sidak
get_color = Vectorize(function(x1, x2) {
if (x1 < -crit[1] || x1 > crit[1] || x2 < -crit[2] || x2 > crit[2])
blue else yellow
})
plot_normal_dens_2d = function(zmat, ...) {
persp(x1, x2, zmat,
expand = 0.5,
lphi = 0,
xlab = "x1",
ylab = "x2\n",
zlab = "f(x1, x2)",
shade = 0.2,
ticktype = "detailed",
...)
}
## ---- out.width = "55%"---------------------------------------------------------------------------
op = par(mar = c(0, 2, 0, 0))
plot_normal_dens_2d(zmat, phi = 45, col = get_color(gc1, gc2))
par(op)
## ---- include = FALSE-----------------------------------------------------------------------------
crit.Bonf = -rep(qnorm(alpha.Bonf / 2), 2)
crit.Sidak = -rep(qnorm(alpha.Sidak / 2), 2)
crit = crit.Bonf
crit = crit.Sidak
calc2Dprob = function(crit, sigma) {
calc_prob = function(lower, upper) {
as.numeric(mvtnorm::pmvnorm(lower = lower, upper = upper, sigma = sigma))
}
left_stripe = calc_prob(lower = c(-Inf, -Inf), upper = c(-crit[1], Inf))
lower_rectangle = calc_prob(lower = c(-crit[1], -Inf), upper = c(crit[1], -crit[2]))
2 * left_stripe + 2 * lower_rectangle # return two-sided prob
}
p.edge = as.numeric(mvtnorm::pmvnorm(lower = c(-Inf, -Inf),
upper = c(-crit[1], -crit[2]),
sigma = Sigma0))
p2D.Bonf = calc2Dprob(crit.Bonf, sigma = Sigma0)
p2D.Sidak = calc2Dprob(crit.Sidak, sigma = Sigma0)
## ---- out.width = "55%"---------------------------------------------------------------------------
orange = okabe_palette[["orange"]]
crit = crit.Sidak
get_color_edges = Vectorize(function(x1, x2) {
# Make edges grey
if (x1 < -crit[1] && x2 < -crit[2] ||
x1 < -crit[1] && x2 > crit[2] ||
x1 > crit[1] && x2 < -crit[2] ||
x1 > crit[1] && x2 > crit[2])
return("grey")
if (x1 < -crit[1] || x1 > crit[1] || x2 < -crit[2] || x2 > crit[2])
blue else yellow
})
op = par(mar = c(0, 2, 0, 0))
plot_normal_dens_2d(zmat, phi = 45, col = get_color_edges(gc1, gc2))
par(op)
## ---- out.width = "35%", message = FALSE, fig.height = 1.3, fig.width = 2.5, fig.align = "center"----
dat = data[1:20, ]
dat[11:20, "x"] = 11:20
dat = rbind(data.frame(x = 1:10, y = 1.5, Sample = "X1"), dat)
ggplot(dat, aes(x = x, y = y)) +
geom_point(size = 10, shape = 21, aes(fill = Sample)) +
scale_fill_manual(values = fill_colors[1:2]) +
lims(x = c(0.5, 20.5), y = c(1.2, 2.2)) +
theme_void() +
theme(legend.position = "bottom")
## ---- out.width = "35%", message = FALSE, fig.height = 1.3, fig.width = 2.5, fig.align = "center", fig.cap = "Schematic view of group sequential two-stage study."----
dat = data[c(1:10, 21:40), ]
ggplot(dat, aes(x = x, y = y)) +
geom_point(size = 10, shape = 21, aes(fill = Sample)) +
scale_fill_manual(values = fill_colors[c(1, 3)]) +
lims(x = c(0.5, 20.5), y = c(1.2, 2.2)) +
theme_void() +
theme(legend.position = "bottom")
## ---- fig.width = 8, fig.height = 3, out.width = "85%", fig.cap = "Normal densities for stage 1 and 2 with Bonferroni bounds."----
normalDensPlotX2Adj = normalDensPlotAdj +
labs(x = expression(x[2]), y = expression(f(x[2]))) +
stat_function(fun = dnorm, color = fill_colors[3], size = 1.5)
gridExtra::grid.arrange(normalDensPlotAdj, normalDensPlotX2Adj, nrow = 1)
## -------------------------------------------------------------------------------------------------
mu = 0
Mean = c(mu, mu)
covariance = sqrt(0.5)
Sigma = matrix(c(1, rep(covariance, 2), 1), 2)
z = mvtnorm::dmvnorm(grid, mean = Mean, sigma = Sigma)
zmat = matrix(z, ncol = length(x1))
## ---- fig.show="hold", out.width = "30%", fig.cap = "Bivariate normal density of two-stage design viewed from different angles."----
op = par(mar = rep(.1, 4))
plot_normal_dens_2d(zmat, phi = 40, col = get_color(gc1, gc2), box = FALSE)
plot_normal_dens_2d(zmat, phi = 65, col = get_color(gc1, gc2), box = FALSE)
plot_normal_dens_2d(zmat, phi = 90, col = get_color(gc1, gc2), box = FALSE)
par(op)
## -------------------------------------------------------------------------------------------------
p2D.Bonf = calc2Dprob(crit.Bonf, sigma = Sigma)
## ---- out.width = "40%", echo = FALSE-------------------------------------------------------------
include_graphics("figures/two-stage-Pocock-bounds.png")
## ---- out.width = "70%", echo = FALSE-------------------------------------------------------------
include_graphics("figures/two-stage-Pocock-bounds-result.png")
## -------------------------------------------------------------------------------------------------
crit.Pocock = rep(2.1783, 2)
p2D.Pocock = calc2Dprob(crit.Pocock, sigma = Sigma)
## ---- echo = FALSE, out.width = "50%", fig.cap = "Group sequential 2-stage Pocock design with one interim look at half of the total samples."----
n = 100
tt = c(.5, 1)
crit = crit.Pocock
x = seq(0.2, 1.0, length = n)
upperLine = data.frame(x = x, y = crit[1])
lowerLine = data.frame(x = x, y = -crit[1])
upperArea = data.frame(x = x, ymin = upperLine$y, ymax = 3, y = 3)
innerArea = data.frame(x = x, ymin = lowerLine$y, ymax = upperLine$y, y = 0)
lowerArea = data.frame(x = x, ymin = -3, ymax = lowerLine$y, y = -3)
p = ggplot(mapping = aes(x = x, y = y)) +
geom_ribbon(data = upperArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = blue) +
geom_ribbon(data = innerArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = yellow) +
geom_ribbon(data = lowerArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = blue) +
geom_line(data = upperLine) +
geom_line(data = lowerLine) +
geom_point(data = data.frame(x = tt, y = crit), size = 5, shape = 21, fill = fill_colors[c(1, 3)]) +
geom_point(data = data.frame(x = tt, y = -crit), size = 5, shape = 21, fill = fill_colors[c(1, 3)]) +
theme_minimal() +
labs(x = "Information Rate", y = "Critical Value") +
lims(y = c(-3, 3)) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.x = element_line())
p
## ---- echo = FALSE, out.width = "50%", fig.cap = "Group sequential 2-stage O'Brien-Fleming design with one interim look at 0.5."----
asOBF <- function(t, alpha = 0.05, side = 2) {
2 * (1 - stats::pnorm((stats::qnorm(1 - (alpha / side)/2)) / sqrt(t)))
}
n = 100
tt = c(.5, 1)
crit.OBF = -qnorm(asOBF(tt))
crit = crit.OBF
x = seq(0.2, 1.0, length = n)
upperLine = data.frame(x = x, y = -qnorm(asOBF(x)))
lowerLine = data.frame(x = x, y = qnorm(asOBF(x)))
upperArea = data.frame(x = x, ymin = upperLine$y, ymax = 6, y = 3)
innerArea = data.frame(x = x, ymin = lowerLine$y, ymax = upperLine$y, y = 0)
lowerArea = data.frame(x = x, ymin = -6, ymax = lowerLine$y, y = -3)
p = ggplot(mapping = aes(x = x, y = y)) +
geom_ribbon(data = upperArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = blue) +
geom_ribbon(data = innerArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = yellow) +
geom_ribbon(data = lowerArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = blue) +
geom_line(data = upperLine) +
geom_line(data = lowerLine) +
geom_point(data = data.frame(x = tt, y = crit), size = 5, shape = 21, fill = fill_colors[c(1, 3)]) +
geom_point(data = data.frame(x = tt, y = -crit), size = 5, shape = 21, fill = fill_colors[c(1, 3)]) +
theme_minimal() +
labs(x = "Information Rate", y = "Critical Value") +
lims(y = c(-6, 6)) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.x = element_line())
p
## ---- out.width = "70%", echo = FALSE-------------------------------------------------------------
include_graphics("figures/two-stage-OBF-bounds-result.png")
## ---- fig.show="hold", out.width = "45%", fig.cap = "Bivariate normal density with critical bounds of Pocock (left) and O'Brien-Fleming design (right) both with interim analysis after 50% of the sample."----
op = par(mar = rep(1.5, 4))
crit = crit.Pocock
plot_normal_dens_2d(zmat, phi = 80, col = get_color(gc1, gc2))
crit = crit.OBF
plot_normal_dens_2d(zmat, phi = 80, col = get_color(gc1, gc2))
par(op)
## ---- fig.show="hold", out.width = "30%", fig.cap = "Bivariate normal density of two-stage O'Brien-Fleming design with first stage after 50% (left), 70% (middle) and 90% (right) of all samples."----
crit = crit.OBF
mu = 0
Mean = c(mu, mu)
Sigma = matrix(c(1, rep(sqrt(0.5), 2), 1), 2)
z = mvtnorm::dmvnorm(grid, mean = Mean, sigma = Sigma)
zmat0.5 = matrix(z, ncol = length(x1))
Sigma = matrix(c(1, rep(sqrt(0.7), 2), 1), 2)
z = mvtnorm::dmvnorm(grid, mean = Mean, sigma = Sigma)
zmat0.7 = matrix(z, ncol = length(x1))
Sigma = matrix(c(1, rep(sqrt(0.9), 2), 1), 2)
z = mvtnorm::dmvnorm(grid, mean = Mean, sigma = Sigma)
zmat0.9 = matrix(z, ncol = length(x1))
op = par(mar = rep(.1, 4))
plot_normal_dens_2d(zmat0.5, phi = 90, col = get_color(gc1, gc2), box = FALSE)
plot_normal_dens_2d(zmat0.7, phi = 90, col = get_color(gc1, gc2), box = FALSE)
plot_normal_dens_2d(zmat0.9, phi = 90, col = get_color(gc1, gc2), box = FALSE)
par(op)
## ---- include = FALSE---------------------------------------------------------
options(old)
Try the GroupSeq package in your browser
Any scripts or data that you put into this service are public.
GroupSeq documentation built on Dec. 28, 2022, 1:23 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.
## ----knitr-setup, include = FALSE-----------------------------------------------------------------
library(knitr)
library("ggplot2")
knitr::opts_chunk$set(
comment = "#",
prompt = F,
tidy = FALSE,
cache = FALSE,
collapse = T,
echo = FALSE,
dpi = 300,
fig.width = 5, fig.height = 5
)
okabe_palette <- c('orange' = "#E69F00",
'sky blue' = "#56B4E9",
'bluish green' = "#009E73",
'yellow' = "#F0E442",
'blue' = "#0072B2",
'vermillion' = "#D55E00",
'reddish purple' = "#CC79A7")
old <- options(width = 100L, digits = 10)
## -------------------------------------------------------------------------------------------------
fill_colors = c(X1 = okabe_palette[["bluish green"]],
X2 = okabe_palette[["vermillion"]],
"X1+X2" = "#6B7E3A")
data = data.frame(x = c(1:10, 1:10, 1:20),
y = c(rep(2, 10), rep(1.5, 10), rep(1.5, 20)),
Sample = c(rep("X1", 10), rep("X2", 10), rep("X1+X2", 20))
)
## ---- out.width = "35%", message = FALSE, fig.height = 1, fig.width = 2.5, fig.align = "center", fig.cap = "Schematic view of a study sample - each circle represents some measurement."----
dat = data[1:10, ]
ggplot(dat, aes(x = x, y = y)) +
geom_point(size = 10, shape = 21, aes(fill = Sample)) +
scale_fill_manual(values = fill_colors[1]) +
lims(x = c(0.5, 10.5), y = c(2, 2)) +
theme_void() +
theme(legend.position = "bottom")
## ---- echo = FALSE, out.width = "65%", fig.cap = fig.cap, fig.height=3----------------------------
fig.cap = "Standard Normal Density with critial bounds marking 5% significance level."
blue = okabe_palette[["sky blue"]]
yellow = okabe_palette[["yellow"]]
crit = qnorm(1 - 0.025)
probLabel = paste0(pnorm(-crit)*100, "%")
x = seq(-3, 3, by = 0.01)
normalDensPlot =
ggplot(data.frame(x = x), aes(x = x)) +
stat_function(fun = dnorm, geom = "area", fill = blue, xlim = c(-4, -crit)) +
stat_function(fun = dnorm, geom = "area", fill = yellow, xlim = c(-crit, crit)) +
stat_function(fun = dnorm, geom = "area", fill = blue, xlim = c(crit, 4)) +
stat_function(fun = dnorm, color = okabe_palette[["bluish green"]], size = 1.5) +
geom_segment(x = 2.2, xend = 2.7, y = 0.02, yend = 0.07) +
geom_segment(x = -2.2, xend = -2.7, y = 0.02, yend = 0.07) +
annotate("text", label = probLabel, x = 2.9, y = 0.09, size = 3.5) +
annotate("text", label = probLabel, x = -2.9, y = 0.09, size = 3.5) +
lims(x = range(x)) +
labs(x = expression(x[1]),
y = expression(f(x[1]))) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.x = element_line())
normalDensPlot
## ---- echo = TRUE---------------------------------------------------------------------------------
qnorm(0.025)
## ---- out.width = "35%", message = FALSE, fig.height = 1.3, fig.width = 2.5, fig.align = "center"----
dat = data[1:20, ]
ggplot(dat, aes(x = x, y = y)) +
geom_point(size = 10, shape = 21, aes(fill = Sample)) +
scale_fill_manual(values = fill_colors[1:2]) +
lims(x = c(0.5, 10.5), y = c(1.2, 2.2)) +
theme_void() +
theme(legend.position = "bottom")
## ---- fig.width = 8, fig.height = 3, out.width = "85%"--------------------------------------------
normalDensPlotX2 = normalDensPlot +
labs(x = expression(x[2]), y = expression(f(x[2]))) +
stat_function(fun = dnorm, color = okabe_palette[["vermillion"]], size = 1.5)
gridExtra::grid.arrange(normalDensPlot, normalDensPlotX2, nrow = 1)
## ---- echo = TRUE---------------------------------------------------------------------------------
qnorm(0.025/2)
## ---- echo = FALSE, out.width = "65%", fig.cap = fig.cap, fig.height=3----------------------------
fig.cap = "Standard Normal Density with critial bounds marking 5% significance level."
blue = okabe_palette[["sky blue"]]
yellow = okabe_palette[["yellow"]]
crit = qnorm(1 - 0.025/2)
probLabel = paste0(pnorm(-crit)*100, "%")
x = seq(-3, 3, by = 0.01)
normalDensPlotAdj =
ggplot(data.frame(x = x), aes(x = x)) +
stat_function(fun = dnorm, geom = "area", fill = blue, xlim = c(-4, -crit)) +
stat_function(fun = dnorm, geom = "area", fill = yellow, xlim = c(-crit, crit)) +
stat_function(fun = dnorm, geom = "area", fill = blue, xlim = c(crit, 4)) +
stat_function(fun = dnorm, color = okabe_palette[["bluish green"]], size = 1.2) +
geom_segment(x = 2.2+.2, xend = 2.7+.2, y = 0.02, yend = 0.07) +
geom_segment(x = -2.2-.2, xend = -2.7-.2, y = 0.02, yend = 0.07) +
annotate("text", label = probLabel, x = 2.9, y = 0.09, size = 3.5) +
annotate("text", label = probLabel, x = -2.9, y = 0.09, size = 3.5) +
lims(x = range(x)) +
labs(x = expression(x[1]),
y = expression(f(x[1]))) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.x = element_line())
## ---- fig.width = 8, fig.height = 3, out.width = "85%"--------------------------------------------
normalDensPlotX2Adj = normalDensPlotAdj +
labs(x = expression(x[2]), y = expression(f(x[2]))) +
stat_function(fun = dnorm, color = okabe_palette[["vermillion"]], size = 1.5)
gridExtra::grid.arrange(normalDensPlotAdj, normalDensPlotX2Adj, nrow = 1)
## -------------------------------------------------------------------------------------------------
alpha = 0.05
alpha.Bonf = alpha / 2
alpha.Sidak = 1 - (1-alpha)^(1/2)
crit.Bonf = -rep(qnorm(alpha.Bonf / 2), 2)
crit.Sidak = -rep(qnorm(alpha.Sidak / 2), 2)
## ----prepare-2D-normal-density-plot, echo = FALSE, include = FALSE--------------------------------
mu = 0
Mean = c(mu, mu)
covariance = 0
Sigma0 = matrix(c(1, rep(covariance, 2), 1), 2)
width = 6.6
binwidth = 0.2
x1 = x2 = seq(from = mu - width/2, to = mu + width/2, by = binwidth)
grid = expand.grid(x1, x2)
z = mvtnorm::dmvnorm(grid, mean = Mean, sigma = Sigma0)
zmat = matrix(z, ncol = length(x1))
grid.col = expand.grid(x1[-1] - binwidth/2, x2[-1] - binwidth/2)
gc1 = grid.col[, 1]
gc2 = grid.col[, 2]
crit = crit.Sidak
get_color = Vectorize(function(x1, x2) {
if (x1 < -crit[1] || x1 > crit[1] || x2 < -crit[2] || x2 > crit[2])
blue else yellow
})
plot_normal_dens_2d = function(zmat, ...) {
persp(x1, x2, zmat,
expand = 0.5,
lphi = 0,
xlab = "x1",
ylab = "x2\n",
zlab = "f(x1, x2)",
shade = 0.2,
ticktype = "detailed",
...)
}
## ---- out.width = "55%"---------------------------------------------------------------------------
op = par(mar = c(0, 2, 0, 0))
plot_normal_dens_2d(zmat, phi = 45, col = get_color(gc1, gc2))
par(op)
## ---- include = FALSE-----------------------------------------------------------------------------
crit.Bonf = -rep(qnorm(alpha.Bonf / 2), 2)
crit.Sidak = -rep(qnorm(alpha.Sidak / 2), 2)
crit = crit.Bonf
crit = crit.Sidak
calc2Dprob = function(crit, sigma) {
calc_prob = function(lower, upper) {
as.numeric(mvtnorm::pmvnorm(lower = lower, upper = upper, sigma = sigma))
}
left_stripe = calc_prob(lower = c(-Inf, -Inf), upper = c(-crit[1], Inf))
lower_rectangle = calc_prob(lower = c(-crit[1], -Inf), upper = c(crit[1], -crit[2]))
2 * left_stripe + 2 * lower_rectangle # return two-sided prob
}
p.edge = as.numeric(mvtnorm::pmvnorm(lower = c(-Inf, -Inf),
upper = c(-crit[1], -crit[2]),
sigma = Sigma0))
p2D.Bonf = calc2Dprob(crit.Bonf, sigma = Sigma0)
p2D.Sidak = calc2Dprob(crit.Sidak, sigma = Sigma0)
## ---- out.width = "55%"---------------------------------------------------------------------------
orange = okabe_palette[["orange"]]
crit = crit.Sidak
get_color_edges = Vectorize(function(x1, x2) {
# Make edges grey
if (x1 < -crit[1] && x2 < -crit[2] ||
x1 < -crit[1] && x2 > crit[2] ||
x1 > crit[1] && x2 < -crit[2] ||
x1 > crit[1] && x2 > crit[2])
return("grey")
if (x1 < -crit[1] || x1 > crit[1] || x2 < -crit[2] || x2 > crit[2])
blue else yellow
})
op = par(mar = c(0, 2, 0, 0))
plot_normal_dens_2d(zmat, phi = 45, col = get_color_edges(gc1, gc2))
par(op)
## ---- out.width = "35%", message = FALSE, fig.height = 1.3, fig.width = 2.5, fig.align = "center"----
dat = data[1:20, ]
dat[11:20, "x"] = 11:20
dat = rbind(data.frame(x = 1:10, y = 1.5, Sample = "X1"), dat)
ggplot(dat, aes(x = x, y = y)) +
geom_point(size = 10, shape = 21, aes(fill = Sample)) +
scale_fill_manual(values = fill_colors[1:2]) +
lims(x = c(0.5, 20.5), y = c(1.2, 2.2)) +
theme_void() +
theme(legend.position = "bottom")
## ---- out.width = "35%", message = FALSE, fig.height = 1.3, fig.width = 2.5, fig.align = "center", fig.cap = "Schematic view of group sequential two-stage study."----
dat = data[c(1:10, 21:40), ]
ggplot(dat, aes(x = x, y = y)) +
geom_point(size = 10, shape = 21, aes(fill = Sample)) +
scale_fill_manual(values = fill_colors[c(1, 3)]) +
lims(x = c(0.5, 20.5), y = c(1.2, 2.2)) +
theme_void() +
theme(legend.position = "bottom")
## ---- fig.width = 8, fig.height = 3, out.width = "85%", fig.cap = "Normal densities for stage 1 and 2 with Bonferroni bounds."----
normalDensPlotX2Adj = normalDensPlotAdj +
labs(x = expression(x[2]), y = expression(f(x[2]))) +
stat_function(fun = dnorm, color = fill_colors[3], size = 1.5)
gridExtra::grid.arrange(normalDensPlotAdj, normalDensPlotX2Adj, nrow = 1)
## -------------------------------------------------------------------------------------------------
mu = 0
Mean = c(mu, mu)
covariance = sqrt(0.5)
Sigma = matrix(c(1, rep(covariance, 2), 1), 2)
z = mvtnorm::dmvnorm(grid, mean = Mean, sigma = Sigma)
zmat = matrix(z, ncol = length(x1))
## ---- fig.show="hold", out.width = "30%", fig.cap = "Bivariate normal density of two-stage design viewed from different angles."----
op = par(mar = rep(.1, 4))
plot_normal_dens_2d(zmat, phi = 40, col = get_color(gc1, gc2), box = FALSE)
plot_normal_dens_2d(zmat, phi = 65, col = get_color(gc1, gc2), box = FALSE)
plot_normal_dens_2d(zmat, phi = 90, col = get_color(gc1, gc2), box = FALSE)
par(op)
## -------------------------------------------------------------------------------------------------
p2D.Bonf = calc2Dprob(crit.Bonf, sigma = Sigma)
## ---- out.width = "40%", echo = FALSE-------------------------------------------------------------
include_graphics("figures/two-stage-Pocock-bounds.png")
## ---- out.width = "70%", echo = FALSE-------------------------------------------------------------
include_graphics("figures/two-stage-Pocock-bounds-result.png")
## -------------------------------------------------------------------------------------------------
crit.Pocock = rep(2.1783, 2)
p2D.Pocock = calc2Dprob(crit.Pocock, sigma = Sigma)
## ---- echo = FALSE, out.width = "50%", fig.cap = "Group sequential 2-stage Pocock design with one interim look at half of the total samples."----
n = 100
tt = c(.5, 1)
crit = crit.Pocock
x = seq(0.2, 1.0, length = n)
upperLine = data.frame(x = x, y = crit[1])
lowerLine = data.frame(x = x, y = -crit[1])
upperArea = data.frame(x = x, ymin = upperLine$y, ymax = 3, y = 3)
innerArea = data.frame(x = x, ymin = lowerLine$y, ymax = upperLine$y, y = 0)
lowerArea = data.frame(x = x, ymin = -3, ymax = lowerLine$y, y = -3)
p = ggplot(mapping = aes(x = x, y = y)) +
geom_ribbon(data = upperArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = blue) +
geom_ribbon(data = innerArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = yellow) +
geom_ribbon(data = lowerArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = blue) +
geom_line(data = upperLine) +
geom_line(data = lowerLine) +
geom_point(data = data.frame(x = tt, y = crit), size = 5, shape = 21, fill = fill_colors[c(1, 3)]) +
geom_point(data = data.frame(x = tt, y = -crit), size = 5, shape = 21, fill = fill_colors[c(1, 3)]) +
theme_minimal() +
labs(x = "Information Rate", y = "Critical Value") +
lims(y = c(-3, 3)) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.x = element_line())
p
## ---- echo = FALSE, out.width = "50%", fig.cap = "Group sequential 2-stage O'Brien-Fleming design with one interim look at 0.5."----
asOBF <- function(t, alpha = 0.05, side = 2) {
2 * (1 - stats::pnorm((stats::qnorm(1 - (alpha / side)/2)) / sqrt(t)))
}
n = 100
tt = c(.5, 1)
crit.OBF = -qnorm(asOBF(tt))
crit = crit.OBF
x = seq(0.2, 1.0, length = n)
upperLine = data.frame(x = x, y = -qnorm(asOBF(x)))
lowerLine = data.frame(x = x, y = qnorm(asOBF(x)))
upperArea = data.frame(x = x, ymin = upperLine$y, ymax = 6, y = 3)
innerArea = data.frame(x = x, ymin = lowerLine$y, ymax = upperLine$y, y = 0)
lowerArea = data.frame(x = x, ymin = -6, ymax = lowerLine$y, y = -3)
p = ggplot(mapping = aes(x = x, y = y)) +
geom_ribbon(data = upperArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = blue) +
geom_ribbon(data = innerArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = yellow) +
geom_ribbon(data = lowerArea, aes(x = x, ymin = ymin, ymax = ymax, y = y), fill = blue) +
geom_line(data = upperLine) +
geom_line(data = lowerLine) +
geom_point(data = data.frame(x = tt, y = crit), size = 5, shape = 21, fill = fill_colors[c(1, 3)]) +
geom_point(data = data.frame(x = tt, y = -crit), size = 5, shape = 21, fill = fill_colors[c(1, 3)]) +
theme_minimal() +
labs(x = "Information Rate", y = "Critical Value") +
lims(y = c(-6, 6)) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.x = element_line())
p
## ---- out.width = "70%", echo = FALSE-------------------------------------------------------------
include_graphics("figures/two-stage-OBF-bounds-result.png")
## ---- fig.show="hold", out.width = "45%", fig.cap = "Bivariate normal density with critical bounds of Pocock (left) and O'Brien-Fleming design (right) both with interim analysis after 50% of the sample."----
op = par(mar = rep(1.5, 4))
crit = crit.Pocock
plot_normal_dens_2d(zmat, phi = 80, col = get_color(gc1, gc2))
crit = crit.OBF
plot_normal_dens_2d(zmat, phi = 80, col = get_color(gc1, gc2))
par(op)
## ---- fig.show="hold", out.width = "30%", fig.cap = "Bivariate normal density of two-stage O'Brien-Fleming design with first stage after 50% (left), 70% (middle) and 90% (right) of all samples."----
crit = crit.OBF
mu = 0
Mean = c(mu, mu)
Sigma = matrix(c(1, rep(sqrt(0.5), 2), 1), 2)
z = mvtnorm::dmvnorm(grid, mean = Mean, sigma = Sigma)
zmat0.5 = matrix(z, ncol = length(x1))
Sigma = matrix(c(1, rep(sqrt(0.7), 2), 1), 2)
z = mvtnorm::dmvnorm(grid, mean = Mean, sigma = Sigma)
zmat0.7 = matrix(z, ncol = length(x1))
Sigma = matrix(c(1, rep(sqrt(0.9), 2), 1), 2)
z = mvtnorm::dmvnorm(grid, mean = Mean, sigma = Sigma)
zmat0.9 = matrix(z, ncol = length(x1))
op = par(mar = rep(.1, 4))
plot_normal_dens_2d(zmat0.5, phi = 90, col = get_color(gc1, gc2), box = FALSE)
plot_normal_dens_2d(zmat0.7, phi = 90, col = get_color(gc1, gc2), box = FALSE)
plot_normal_dens_2d(zmat0.9, phi = 90, col = get_color(gc1, gc2), box = FALSE)
par(op)
## ---- include = FALSE---------------------------------------------------------
options(old)
Try the GroupSeq package in your browser
Any scripts or data that you put into this service are public.
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.