Nothing
## ----"knitr-options", echo = FALSE--------------------------------------------
library(knitr)
opts_chunk$set(size = "small",
fig.height = 4,
fig.align = "center",
cache = FALSE,
message = FALSE,
warning = FALSE)
## ----"data-loading"-----------------------------------------------------------
library(ReplicationSuccess)
data("RProjects")
str(RProjects, width = 80, strict.width = "cut")
## computing zo, zr, c
RProjects$zo <- with(RProjects, fiso/se_fiso)
RProjects$zr <- with(RProjects, fisr/se_fisr)
RProjects$c <- with(RProjects, se_fiso^2/se_fisr^2)
## ----"plot-projects", fig.height = 6------------------------------------------
## plots of effect estimates
par(mfrow = c(2, 2), las = 1, pty = "s",
mar = c(4, 2.5, 2.5, 1))
for (p in unique(RProjects$project)) {
data_project <- subset(RProjects, project == p)
significant <- ifelse(data_project$pr1 < 0.025, "#DF536BCC", "#00000080")
plot(rr ~ ro, data = data_project, ylim = c(-0.5, 1), col = significant,
xlim = c(-0.5, 1), main = p, xlab = expression(italic(r)[o]),
cex = 0.7, pch = 19, ylab = expression(italic(r)[r]))
legend("topleft", legend = "replication significant", cex = 0.8, pch = 20,
col = "#DF536BCC", bty = "n")
abline(h = 0, lty = 2)
abline(a = 0, b = 1, col = "grey")
}
## ----forestPlots, fig.height = 5----------------------------------------------
data("protzko2020")
## forestplots of effect estimates
parOld <- par(mar = c(5, 8, 4, 2), mfrow = c(2, 2))
experiments <- unique(protzko2020$experiment)[1:4]
for (ex in experiments) {
## compute CIs
dat <- subset(protzko2020, experiment == ex)
za <- qnorm(p = 0.975)
plotDF <- data.frame(lower = dat$smd - za*dat$se,
est = dat$smd,
upper = dat$smd + za*dat$se)
colpalette <- c("#000000", "#1B9E77", "#D95F02")
cols <- colpalette[dat$type]
yseq <- seq(1, nrow(dat))
## forestplot
plot(x = plotDF$est, y = yseq, xlim = c(-0.15, 0.8),
ylim = c(0.8*min(yseq), 1.05*max(yseq)), type = "n",
yaxt = "n", xlab = "Effect estimate (SMD)", ylab = "")
abline(v = 0, col = "#0000004D")
arrows(x0 = plotDF$lower, x1 = plotDF$upper, y0 = yseq, angle = 90,
code = 3, length = 0.05, col = cols)
points(y = yseq, x = plotDF$est, pch = 20, lwd = 2, col = cols)
axis(side = 2, at = yseq, las = 1, labels = dat$type, cex.axis = 0.85)
title(main = ex)
}
## -----------------------------------------------------------------------------
for (p in unique(RProjects$project)) {
data_project <- subset(RProjects, project == p)
significant_O <- data_project$po1 < 0.025
significant_R <- data_project$pr1 < 0.025
success <- significant_O & significant_R
cat(paste0(p, ": \n"))
cat(paste0(round(mean(significant_O)*100, 1), "% original studies significant (",
sum(significant_O), "/", length(significant_O), ")\n"))
cat(paste0(round(mean(significant_R)*100, 1), "% replications significant (",
sum(significant_R), "/", length(significant_R), ")\n"))
cat(paste0(round(mean(success)*100, 1),
"% both studies significant in the same direction (",
sum(success), "/", length(success), ")\n \n"))
}
## -----------------------------------------------------------------------------
sampleSizeSignificance(zo = 2.5, power = 0.8, level = 0.05, designPrior = "conditional")
sampleSizeSignificance(zo = 2.5, power = 0.8, level = 0.05, designPrior = "predictive")
sampleSizeSignificance(zo = 2.5, power = 0.8, level = 0.05, designPrior = "conditional",
shrinkage = 0.25)
## ----"plot-powerSignificance", echo = FALSE, fig.height = 4-------------------
po <- seq(0.0001, 0.05, 0.0001)/2
## plot power
plot(po, powerSignificance(zo = p2z(po, alternative = "one.sided"),
designPrior = "conditional")*100,
type = "l", ylim = c(0, 100), lwd = 1.5, ylab = "Power (%)",
xlab = expression(italic(p)[o]), las = 1, yaxt = "n")
axis(side = 2, at = seq(0, 100, 25), las = 1)
axis(side = 3, at = seq(0.0, 0.025, by = 0.005),
labels = c("", round(p2z(p = seq(0.005, 0.025, by = 0.005),
alternative = "one.sided"), 2)))
mtext(text = expression(paste(italic(z)[o])), side = 3, line = 2)
abline(h = seq(0, 100, 25), lty = 1, col = adjustcolor("lightgrey", alpha.f = 0.3))
# abline(h = 50, lty = 1, lwd = 1.5, col = adjustcolor("lightgrey", alpha.f = 0.4))
# abline(h = 50, col = "#333333B3", lty = 3)
lines(po, powerSignificance(zo = p2z(po, alternative = "one.sided"),
designPrior = "predictive")*100, lwd = 2, lty = 2)
legend("topright", legend = c("conditional", "predictive"),
title = "Design prior", lty = c(1, 2), lwd = 1.5, bty = "n")
## ----"plot-sampleSizeSignificance", echo = FALSE, fig.height = 4--------------
## plot sample size
plot(po, sampleSizeSignificance(zo = p2z(po, alternative = "one.sided"),
designPrior = "conditional", power = 0.8),
type = "l", ylim = c(0.5, 8), log = "y", lwd = 1.5,
ylab = expression(paste("Relative sample size ", n[r]/n[o])),
xlab = expression(italic(p)[o]), las = 1, yaxt = "n")
axis(side = 3, at = seq(0.0, 0.025, by = 0.005),
labels = c("", round(p2z(p = seq(0.005, 0.025, by = 0.005),
alternative = "one.sided"), 2)))
axis(side = 2, las = 1, at = c(0.5, 1, 2, 4, 8, 16, 32),
labels = c("1/2", "1", "2", "4", "8", "16", "32"))
mtext(text = expression(paste(italic(z)[o])), side = 3, line = 2)
abline(h = c(0.5, 1, 2, 4, 8), lty = 1, col = adjustcolor("lightgrey", alpha.f = 0.3))
# abline(h = 1, col = "#333333B3", lty = 3)
abline(h = 1, lty = 1, lwd = 1.5, col = adjustcolor("lightgrey", alpha.f = 0.4))
lines(po, sampleSizeSignificance(zo = p2z(po, alternative = "one.sided"),
designPrior = "predictive", power = 0.8),
lwd = 2, lty = 2)
legend("topleft", legend = c("conditional", "predictive"),
title = "Design prior", lty = c(1, 2), lwd = 1.5, bty = "n")
## ----"plot-predictionInterval", fig.height = 6--------------------------------
## compute prediction intervals for replication projects
par(mfrow = c(2, 2), las = 1, mai = rep(0.65, 4))
for (p in unique(RProjects$project)) {
data_project <- subset(RProjects, project == p)
pQ <- Qtest(thetao = data_project$fiso,
thetar = data_project$fisr,
seo = data_project$se_fiso,
ser = data_project$se_fisr)
PI <- predictionInterval(thetao = data_project$fiso,
seo = data_project$se_fiso,
ser = data_project$se_fisr)
## transforming back to correlation scale
PI <- tanh(PI)
incompatible <- pQ < 0.05 ## incompatible at 5% level
color <- ifelse(incompatible == FALSE, "#00000099", "#DF536BCC")
study <- seq(1, nrow(data_project))
plot(data_project$rr, study, col = color, pch = 20, cex = 0.5,
xlim = c(-0.5, 1), xlab = expression(italic(r)[r]), ylab = "Study",
main = paste0(p, ": ", round(mean(incompatible)*100, 0), "% incompatible"))
arrows(PI$lower, study, PI$upper, study, length = 0.02,
angle = 90, code = 3, col = color)
abline(v = 0, lty = 3)
}
## ----example-morewedge, echo = FALSE------------------------------------------
## data from study
morewedge <- subset(RProjects, study == "Morewedge et al. (2010), Science")
## Functions
ssv <- function(zo, so, a) {
tau2 <- so^2/(zo^2/qnorm(p = a/2, lower.tail = FALSE)^2 - 1)
return(tau2)
}
## Parameters and computations
optionsOld <- options(scipen = 5)
theta_o <- morewedge$fiso
theta_r <- morewedge$fisr
so <- morewedge$se_fiso
sr <- morewedge$se_fisr
c <- so^2/sr^2
zo <- theta_o/so
zr <- theta_r/sr
alpha <- 0.05
za <- qnorm(p = alpha/2, lower.tail = FALSE)
ps <- signif(pSceptical(zo = zo, zr = zr, c = c, alternative = "two.sided", type = "nominal"), 2)
ps1 <- signif(pSceptical(zo = zo, zr = zr, c = c, alternative = "one.sided", type = "nominal"), 2)
ps1r <- signif(pSceptical(zo = zo, zr = zr, c = c, alternative = "one.sided", type = "golden"), 2)
tau <- sqrt(ssv(zo, so, alpha)) # sufficiently sceptical prior sd at alpha = 0.05
s2_p <- 1/(1/so^2 + 1/tau^2) # posterior variance
mu_p <- s2_p*theta_o/so^2 # posterior mean
## ----"plot-pSceptical", echo = FALSE, fig.height = 4--------------------------
## Plot
df <- data.frame(estimate = factor(c("Original Study",
# "Posterior",
# "Sceptical Prior",
"Replication Study"),
levels = c("Original Study",
# "Posterior",
# "Sceptical Prior",
"Replication Study")),
x = c(1, 2),
col = c(1, 4),
theta = c(theta_o,
theta_r),
lower = c(theta_o - za*so,
theta_r - za*sr),
upper = c(theta_o + za*so,
theta_r + za*sr),
p = c(signif(2*pnorm(q = zo, lower.tail = FALSE), 1),
signif(2*pnorm(q = zr, lower.tail = FALSE), 2)))
ylims <- c(-0.1, 1)
xlims <- c(0.5, 2.5)
cex <- 1
plot(x = df$x, y = df$theta, pch = " ",
xlim = xlims, ylim = ylims,
xlab = "", ylab = "Effect Size", xaxt = "n", las = 1)
axis(side = 1, at = df$x, labels = df$estimate,
cex.axis = 0.9, mgp = c(3, 2, 0))
abline(h = seq(-2, 2, 0.25), lty = 1, col = adjustcolor("lightgrey", alpha.f = 0.4))
abline(h = 0, lty = 2)
# ## one-sided CIs
# arrows(x0 = df$x[-3], x1 = df$x[-3], y0 = df$theta[-3] + 0.05, y1 = df$upper[-3],
# col = df$col[-3], code = 2, angle = 90, length = 0, lwd = 1.5)
# points(x = df$x[-3], y = df$upper[-3], pch = 17, cex = 1.3*cex, col = df$col[-3])
## two-sided CI
arrows(x0 = df$x, x1 = df$x, y0 = df$theta + 0.05, y1 = df$upper,
col = df$col, code = 2, angle = 90, length = 0.1, lwd = 1.5)
arrows(x0 = df$x, x1 = df$x, y0 = df$theta - 0.05, y1 = df$lower,
col = df$col, code = 2, angle = 90, length = 0.1, lwd = 1.5)
## arrows and text
points(x = df$x, y = df$theta, pch = 20, cex = cex*1.5, col = df$col)
text(x = df$x[1] - 0.2, y = df$theta[1],
labels = bquote(paste(hat(theta)[o], " = ", .(round(df$theta[1], 2)))))
text(x = df$x[2] + 0.2, y = df$theta[2],
labels = bquote(paste(hat(theta)[r], " = ", .(round(df$theta[2], 2)))))
text(x = df$x[1], y = df$upper[1] + 0.1,
labels = bquote(paste(p[o], " = ", .(round(morewedge$po1, 4)))))
text(x = df$x[2], df$upper[2] + 0.1,
labels = bquote(paste(p[r], " = ", .(round(morewedge$pr1, 4)))))
text(x = 1.5, y = 0.3,
labels = bquote(paste(p[S], " = ", .(ps1))))
# arrows(x0 = 2.1, x1 = 2.9, y0 = 0.9, y1 = 0.9, length = 0.05, lwd = 1.2, col = "#000000B2")
# text(x = 2.5, y = 0.98, labels = "reverse-Bayes", cex = 0.75, col = "#000000B2")
# arrows(x0 = 3.15, x1 = 4.5, y0 = 0.9, y1 = 0.9, length = 0.05, lwd = 1.2, col = "#000000B2")
# text(x = 3.8, y = 0.98, labels = "prior-data conflict assessment", cex = 0.75, col = "#000000B2")
# points(x = 2, y = 0, pch = 1, cex = 2)
# arrows(x0 = 1.85, x1 = 1.98, y0 = -0.2, y1 = -0.03, length = 0.05, col = "#000000B2")
# text(x = 1.85, y = -0.25, labels = "fixed at zero", col = "#000000B2", cex = 0.75)
box()
## ----morewedge, echo = TRUE---------------------------------------------------
morewedge <- subset(RProjects, study == "Morewedge et al. (2010), Science")
## ----echo = FALSE-------------------------------------------------------------
ps_gold <- signif(pSceptical(zo = zo, zr = zr, c = c, alternative = "one.sided", type = "golden"), 2)
ps_contr <- signif(pSceptical(zo = zo, zr = zr, c = c, alternative = "one.sided", type = "controlled"), 2)
## ----echo = TRUE--------------------------------------------------------------
print(pS_nominal <- pSceptical(zo = morewedge$zo, zr = morewedge$zr,
c = morewedge$c, alternative = "one.sided",
type = "nominal"))
print(pS_golden <- pSceptical(zo = morewedge$zo, zr = morewedge$zr,
c = morewedge$c, alternative = "one.sided",
type = "golden"))
print(pS_controlled <- pSceptical(zo = morewedge$zo, zr = morewedge$zr,
c = morewedge$c, alternative = "one.sided",
type = "controlled"))
## ----"plot-pSceptical-projects", fig.height = 4-------------------------------
## computing one-sided golden and controlled sceptical p-value for replication projects
RProjects$psG <- with(RProjects,
pSceptical(zo = zo, zr = zr, c = c,
alternative = "one.sided", type = "golden"))
RProjects$psC <- with(RProjects,
pSceptical(zo = zo, zr = zr, c = c,
alternative = "one.sided", type = "controlled"))
## ----echo = FALSE-------------------------------------------------------------
for (p in unique(RProjects$project)) {
data_project <- subset(RProjects, project == p)
cat(paste0(p, ": \n"))
success_sceptG <- (data_project$psG < 0.025)
cat(paste0(round(mean(success_sceptG)*100, 2),
"% smaller than 0.025 (one-sided golden sceptical p-value) \n"))
success_sceptC <- (data_project$psC < 0.025)
cat(paste0(round(mean(success_sceptC)*100, 2),
"% smaller than 0.025 (one-sided controlled sceptical p-value) \n"))
success_tradit <- (data_project$po1 < 0.025) & (data_project$pr1 < 0.025)
cat(paste0(round(mean(success_tradit)*100, 2),
"% smaller than 0.025 (both one-sided traditional p-values) \n"))
if(sum(success_sceptG != success_tradit) > 0){
discrep <- data_project[(success_sceptG != success_tradit),
c("ro", "rr", "c", "po1", "pr1", "psG")]
## print effect estimates, 1sided p-values, and c of discrepant studies
cat("Discrepant studies (golden vs significance): \n")
print(signif(discrep, 2), row.names = FALSE)
}
if(sum(success_sceptC != success_tradit) > 0){
discrep <- data_project[(success_sceptC != success_tradit),
c("ro", "rr", "c", "po1", "pr1", "psC")]
## print effect estimates, 1sided p-values, and c of discrepant studies
cat("Discrepant studies (controlled vs significance): \n")
print(signif(discrep, 2), row.names = FALSE)
}
if(sum(success_sceptC != success_sceptG) > 0){
discrep <- data_project[(success_sceptC != success_sceptG),
c("ro", "rr", "c", "po1", "pr1", "psG", "psC")]
## print effect estimates, 1sided p-values, and c of discrepant studies
cat("Discrepant studies (golden vs controlled): \n")
print(signif(discrep, 2), row.names = FALSE)
}
cat("\n \n")
}
## ----golden-vs-controlled, fig.height = 6-------------------------------------
par(mfrow = c(2, 2), las = 1, pty = "s",
mar = c(4, 2.5, 2.5, 1))
myaxis <- c(10e-5, 0.001, 0.01, 0.1, 1)
for (p in unique(RProjects$project)) {
data_project <- subset(RProjects, project == p)
plot(psG ~ psC, data = data_project, ylim = c(10e-5, 1),
xlim = c(10e-5, 1), main = p, xlab = expression(paste("controlled ", p[S])),
ylab = expression(paste("golden ", p[S])),
pch = 19,
col = "#00000099",
cex = 1.3,
axes = F,
log = "xy")
abline(h = 0, lty = 2)
abline(a = 0, b = 1, col = "grey")
abline(h = 0.025, lty = 2, col = "grey")
abline(v = 0.025, lty = 2, col = "grey")
axis(1, at = myaxis, labels = myaxis)
axis(2,at = myaxis, labels = myaxis)
box()
}
## ----"threshold-p-sceptical2", echo = FALSE-----------------------------------
## computing nominal, controlled, and golden replication success levels
## for one-sided uncalibrated sceptical p-value
thresh_gol <- levelSceptical(level = 0.025, alternative = "one.sided",
type = "golden")
thresh_contr <- levelSceptical(level = 0.025, alternative = "one.sided",
type = "controlled", c = c(1, 2))
thresh_nom <- levelSceptical(level = 0.025, alternative = "one.sided",
type = "nominal")
## ----"threshold-p-sceptical"--------------------------------------------------
## computing nominal, golden and controlled replication success levels
## for one-sided uncalibrated sceptical p-value
print(rs_level_nom <- levelSceptical(level = 0.025, alternative = "one.sided",
type = "nominal"))
print(rs_level_gol <- levelSceptical(level = 0.025, alternative = "one.sided",
type = "golden"))
print(rs_level_contr <- levelSceptical(level = 0.025, alternative = "one.sided",
type = "controlled", c = c(1, 2)))
## -----------------------------------------------------------------------------
sampleSizeReplicationSuccess(zo = 2.5, power = 0.8, level = 0.025,
alternative = "one.sided",
designPrior = "conditional",
type = c("golden", "controlled"))
sampleSizeReplicationSuccess(zo = 2.5, power = 0.8, level = 0.025,
alternative = "one.sided",
designPrior = "predictive",
type = c("golden", "controlled"))
## ----"plot-powerReplicationSuccess", echo = FALSE, fig.height = 4-------------
## plot power
po <- seq(0.001, 0.05, length.out = 100)/2
plot(po, powerReplicationSuccess(zo = p2z(po, alternative = "one.sided"),
designPrior = "conditional",
level = 0.025,
alternative = "one.sided",
type = "golden")*100,
type = "l", ylim = c(0, 100), lwd = 1.5, ylab = "Power (%)", las = 1,
xlab = expression(italic(p)[o]), yaxt = "n")
abline(h = seq(0, 100, 25), lty = 1, col = adjustcolor("lightgrey", alpha.f = 0.3))
# abline(h = 50, lty = 1, lwd = 1.5, col = adjustcolor("lightgrey", alpha.f = 0.4))
axis(side = 2, at = seq(0, 100, 25), las = 1)
axis(side = 3, at = seq(0.0, 0.025, by = 0.005),
labels = c("", round(p2z(p = seq(0.005, 0.025, by = 0.005),
alternative = "one.sided"), 2)))
mtext(text = expression(paste( italic(z)[o])), side = 3, line = 2)
lines(po, powerReplicationSuccess(zo = p2z(po, alternative = "one.sided"),
designPrior = "predictive",
level = 0.025,
alternative = "one.sided",
type = "golden")*100,
lwd = 2, lty = 2)
lines(po, powerReplicationSuccess(zo = p2z(po, alternative = "one.sided"),
designPrior = "conditional",
level = 0.025,
alternative = "one.sided",
type = "controlled")*100,
col = "red", lwd = 2)
lines(po, powerReplicationSuccess(zo = p2z(po, alternative = "one.sided"),
designPrior = "predictive",
level = 0.025,
alternative = "one.sided",
type = "controlled")*100,
col = "red", lty = 2, lwd = 2)
legend("topright", legend = c("conditional", "predictive"),
title = "Design prior", lty = c(1, 2), lwd = 1.5, bty = "n")
# abline(h = 50, lty = 3)
legend("bottomleft",
legend = c(expression(paste("golden ", p[S])),
expression(paste("controlled ", p[S]))),
col = c("black", "red"),
lty = 1, lwd = 2,
bty = "n")
## ----"plot-sampleSizeReplicationSuccess", echo = FALSE, fig.height = 4--------
# po <- seq(0.0001, 0.05, 0.0001)
po <- seq(0.001, 0.05, length.out = 100)/2
## plot sample size
plot(po,
sampleSizeReplicationSuccess(zo = p2z(po, alternative = "one.sided"),
power = 0.8,
level = 0.025,
designPrior = "conditional",
alternative = "one.sided",
type = "golden"),
type = "l", log = "y", lwd = 1.5, las = 1, ylim = c(0.5, 20),
ylab = expression(paste("Relative sample size ", n[r]/n[o])),
xlab = expression(italic(p)[o]), yaxt = "n")
axis(side = 2, las = 1, at = c(0.5, 1, 2, 4, 8, 16, 32),
labels = c("1/2", "1", "2", "4", "8", "16", "32"))
# abline(h = 1, lty = 3)
abline(h = c(0.5, 1, 2, 4, 8, 16, 32), lty = 1, col = adjustcolor("lightgrey", alpha.f = 0.3))
abline(h = 1, lty = 1, lwd = 1.5, col = adjustcolor("lightgrey", alpha.f = 0.4))
axis(side = 3, at = seq(0.0, 0.025, by = 0.005),
labels = c("", round(p2z(p = seq(0.005, 0.025, by = 0.005),
alternative = "one.sided"), 2)))
mtext(text = expression(paste(italic(z)[o])), side = 3, line = 2)
suppressWarnings({
lines(po, sampleSizeReplicationSuccess(zo = p2z(po, alternative = "one.sided"),
power = 0.8, level = 0.025,
designPrior = "predictive",
alternative = "one.sided",
type = "golden"),
lwd = 2, lty = 2)
lines(po,
sampleSizeReplicationSuccess(zo = p2z(po, alternative = "one.sided"),
power = 0.8,
level = 0.025,
designPrior = "conditional",
alternative = "one.sided",
type = "controlled"),
col = "red", lwd = 2)
lines(po,
sampleSizeReplicationSuccess(zo = p2z(po, alternative = "one.sided"),
power = 0.8,
level = 0.025,
designPrior = "predictive",
alternative = "one.sided",
type = "controlled"),
col = "red", lwd = 2, lty = 2)
})
legend("topleft", legend = c("conditional", "predictive"),
title = "Design prior", lty = c(1, 2), lwd = 1.5, bty = "n")
legend("bottomright",
legend = c(expression(paste("golden ", p[S])),
expression(paste("controlled ", p[S]))),
col = c("black", "red"),
lty = 1, lwd = 2,
bty = "n")
## ----echo = FALSE-------------------------------------------------------------
d <- morewedge$fisr/morewedge$fiso
dminrs <- effectSizeReplicationSuccess(zo = morewedge$zo, c = morewedge$c)
dminsign <- effectSizeSignificance(zo = morewedge$zo, c = morewedge$c)
## ----echo = FALSE, fig.height = 3.5, eval = FALSE-----------------------------
# data("SSRP")
# intpow_cond <- with(SSRP, powerSignificanceInterim(zo = fiso/se_fiso, zi = fisi/se_fisi, c = nr/no,
# f = ni/nr, designPrior = "conditional"))
# intpow_pred <- with(SSRP, powerSignificanceInterim(zo = fiso/se_fiso, zi = fisi/se_fisi, c = nr/no,
# f = ni/nr, designPrior = "predictive"))
# plot(intpow_cond*100, intpow_pred*100,
# xlab = "Conditional power (in %)",
# ylab = "Predictive power (in %)",
# pch = 20,
# cex = 1.2,
# xlim = c(80, 100),
# ylim = c(0, 100),
# yaxt = "n", las = 1,
# col = "#00000099")
# axis(side = 2, at = seq(0, 100, 25), las = 1)
# abline(a = 0, b = 1, col = "grey")
## ----"shrinkage", echo = FALSE, fig.height = 3.5------------------------------
zo <- seq(0, 4, 0.01)
s <- pmax(1 - 1/zo^2, 0)
shrinkage <- (1 - s)*100
plot(zo, shrinkage, type = "l", ylim = c(0, 100), las = 1,
xlab = expression(paste(italic(z)[o])), ylab = "Shrinkage (%)",
yaxt = "n")
axis(side = 2, at = seq(0, 100, 25), las = 1)
axis(side = 3, at = seq(0, 4, by = 1),
labels = c(signif(z2p(seq(0, 3, by = 1), alternative = "one.sided"), 2),
signif(z2p(4, alternative = "one.sided"), 1)))
abline(h = seq(0, 100, 25), lty = 1, col = adjustcolor("lightgrey", alpha.f = 0.3))
mtext(text = expression(italic(p)[o]), side = 3, line = 2)
## ----include=FALSE------------------------------------------------------------
options(optionsOld)
par(parOld)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.