inst/doc/ReplicationSuccess.R

In ReplicationSuccess: Design and Analysis of Replication Studies

## ----"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)