sensitivity analysis

paper/rsabsence.Rnw

@@ -52,6 +52,7 @@ Reproducibility <- TRUE
 
 ## packages
 library(ggplot2) # plotting
+library(gridExtra) # combining ggplots
 library(dplyr) # data manipulation
 library(reporttools) # reporting of p-values
 
@@ -700,6 +701,109 @@ which is available in our git repository.% The effect estimates and standard
 
 \bibliography{bibliography}
 
+\appendix
+\begin{appendixbox}
+% \label{first:sensitivity}
+% \section{Sensitivity analyses}
+
+
+\begin{center}
+<< "sensitivity", fig.height = 8 >>=
+## compute number of successful replications as a function of the equivalence margin
+marginseq <- seq(0.01, 4.5, 0.01)
+alphaseq <- c(0.005, 0.05, 0.1)
+sensitivityGrid <- expand.grid(m = marginseq, a = alphaseq)
+equivalenceDF <- lapply(X = seq(1, nrow(sensitivityGrid)), FUN = function(i) {
+    m <- sensitivityGrid$m[i]
+    a <- sensitivityGrid$a[i]
+    rpcbNull$ptosto <- with(rpcbNull, pmax(pnorm(q = smdo, mean = m, sd = so,
+                                                 lower.tail = TRUE),
+                                           pnorm(q = smdo, mean = -m, sd = so,
+                                                 lower.tail = FALSE)))
+    rpcbNull$ptostr <- with(rpcbNull, pmax(pnorm(q = smdr, mean = m, sd = sr,
+                                                 lower.tail = TRUE),
+                                           pnorm(q = smdr, mean = -m, sd = sr,
+                                                 lower.tail = FALSE)))
+    successes <- sum(rpcbNull$ptosto <= a & rpcbNull$ptostr <= a)
+    data.frame(margin = m, alpha = a,
+               successes = successes, proportion = successes/nrow(rpcbNull))
+}) %>%
+    bind_rows()
+
+## plot number of successes as a function of margin
+nmax <- nrow(rpcbNull)
+bks <- seq(0, nmax, round(nmax/5))
+labs <- paste0(bks, " (", bks/nmax*100, "%)")
+plotA <- ggplot(data = equivalenceDF,
+                aes(x = margin, y = successes,
+                    color = factor(alpha, ordered = TRUE))) +
+    facet_wrap(~ "Equivalence test") +
+    geom_vline(xintercept = margin, lty = 2, alpha = 0.4) +
+    geom_step(alpha = 0.9, linewidth = 0.8) +
+    scale_y_continuous(breaks = bks, labels = labs) +
+    ## scale_y_continuous(labels = scales::percent) +
+    guides(color = guide_legend(reverse = TRUE)) +
+    labs(x = bquote("Equivalence margin" ~ Delta),
+         y = "Successful replications",
+         color = bquote(italic("p")["TOST"] ~ "threshold" ~ alpha)) +
+    theme_bw() +
+    theme(panel.grid.minor = element_blank(),
+          panel.grid.major = element_blank(),
+          strip.background = element_rect(fill = alpha("tan", 0.4)),
+          strip.text = element_text(size = 12),
+          legend.position = c(0.85, 0.25),
+          plot.background = element_rect(fill = "transparent", color = NA),
+          ## axis.text.y = element_text(hjust = 0),
+          legend.box.background = element_rect(fill = "transparent", colour = NA))
+
+## compute number of successful replications as a function of the prior scale
+priorsdseq <- seq(0, 40, 0.1)
+bfThreshseq <- c(3, 6, 10)
+sensitivityGrid2 <- expand.grid(s = priorsdseq, thresh = bfThreshseq)
+bfDF <- lapply(X = seq(1, nrow(sensitivityGrid2)), FUN = function(i) {
+    priorsd <- sensitivityGrid2$s[i]
+    thresh <- sensitivityGrid2$thresh[i]
+    rpcbNull$BForig <- with(rpcbNull, BF01(estimate = smdo, se = so, unitvar = priorsd^2))
+    rpcbNull$BFrep <- with(rpcbNull, BF01(estimate = smdr, se = sr, unitvar = priorsd^2))
+    successes <- sum(rpcbNull$BForig >= thresh & rpcbNull$BFrep >= thresh)
+    data.frame(priorsd = priorsd, thresh = thresh,
+               successes = successes, proportion = successes/nrow(rpcbNull))
+}) %>%
+    bind_rows()
+
+## plot number of successes as a function of prior sd
+plotB <- ggplot(data = bfDF,
+                aes(x = priorsd, y = successes, color = factor(thresh, ordered = TRUE))) +
+    facet_wrap(~ "Bayes factor") +
+    geom_vline(xintercept = 4, lty = 2, alpha = 0.4) +
+    geom_step(alpha = 0.9, linewidth = 0.8) +
+    scale_y_continuous(breaks = bks, labels = labs, limits = c(0, nmax)) +
+    ## scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
+    guides(color = guide_legend(reverse = TRUE)) +
+    labs(x = "Prior distribution scale",
+         y = "Successful replications ",
+         color = bquote("BF"["01"] ~ "threshold" ~ gamma)) +
+    theme_bw() +
+    theme(panel.grid.minor = element_blank(),
+          panel.grid.major = element_blank(),
+          strip.background = element_rect(fill = alpha("tan", 0.4)),
+          strip.text = element_text(size = 12),
+          legend.position = c(0.85, 0.25),
+          plot.background = element_rect(fill = "transparent", color = NA),
+          ## axis.text.y = element_text(hjust = 0),
+          legend.box.background = element_rect(fill = "transparent", colour = NA))
+grid.arrange(plotA, plotB, ncol = 1)
+@
+
+\captionof{figure}{Number of successful replications of original null results in
+  the RPCB as a function of the margin $\Delta$ of equivalence test
+  ($p_{\text{TOST}} \leq \alpha$ in both studies) or the scale of the prior
+  distribution for the effect under the alternative $H_{1}$ of the Bayes
+  factor ($\BF_{01} \geq \gamma$ in both studies).}
+\end{center}
+
+\end{appendixbox}
+
 
 << "sessionInfo1", eval = Reproducibility, results = "asis" >>=
 ## print R sessionInfo to see system information and package versions