Figures for "Power and precision" blog post

planning_precision.Rmd

---
title: "Power and precision"
author: "Richard D. Morey"
date: "11/06/2020"
output: 
  html_document:
    dev: png
    self_contained: no
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE,
                      fig.width = 6,
                      fig.height = 4,
                      dpi = 300)
```

```{r}
library(sysfonts)
library(showtext)
font_add_google("Roboto Condensed")
showtext_auto()

find_val_0 = Vectorize(function(power)
  uniroot(interval = c(0,1), function(p)
    pbinom(crit - 1, N, p, lower.tail = FALSE) - power)$root,
  "power")

find_val2 = function(power)
  c(find_val_0(power), find_val_1(power))

find_val_1 = Vectorize(function(power)
  uniroot(interval = c(0,1), function(p)
    pbinom(crit, N, p, lower.tail = TRUE) - power)$root,
  "power")

ci_alpha = Vectorize(function(alpha) 
  binom.test(crit, N, conf.level = 1 - 2*alpha)$conf.int,
"alpha")

find_val_t0 = Vectorize(function(power, N, crit)
  qlogis(uniroot(interval = c(0,1), function(q){
    delta = qlogis(q)
    suppressWarnings(
      pt(crit, N-1, ncp = delta * sqrt(N), lower.tail = FALSE) - power
    )}
  )$root),
  "power")

nb_lower_bound <- function(y, size, alpha = 0.025)
  qbeta(alpha, size, y + 1)


nb_upper_bound <- function(y, size, alpha = 0.025)
  qbeta(alpha, size, y, lower.tail = FALSE)



ci_nbinom <- Vectorize(function(y, size, alpha = .05){
  c(nb_lower_bound(y, size, alpha / 2), nb_upper_bound(y, size, alpha / 2) )
}, "y")
```

```{r}
N = 100
crit = 60
pow_vals = c(0, .01,.05,.1, .25)
col_vals = viridis::viridis(length(pow_vals), alpha = .3)
pow_vals_2 = c(pow_vals, rev(1 - pow_vals))
col_vals_2 = viridis::viridis(length(pow_vals_2), alpha = .5)

x_lim = c(.4,.75)
par(family = "Roboto Condensed", las = 1, xaxs = 'i',
    yaxs = 'i', mar =  c(5.1, 4.1, 4.1, 4.1))


curve(pbinom(crit - 1, N, x, lower.tail = FALSE), x_lim[1], x_lim[2], 1024,
      ylim = c(0, 1),
      axes=FALSE, 
      ylab = substitute(paste("Probability that ",X>=c), list(c = crit)),
      xlab = "True binomial probability parameter p",
      lwd = 2, xpd = TRUE, ty = 'n')
axis(1)
axis(2)

for(i in seq_along(pow_vals_2)[-1]){
  par = find_val_0(pow_vals_2[(i-1):i])
  rect(par[1], par()$usr[3], par[2], par()$usr[4], col = col_vals_2[i],
       border = NA)
  
}

for(i in seq_along(pow_vals_2)[c(-1,-length(pow_vals_2))]){
  pow = pow_vals_2[i]
  par = find_val_0(pow)
  text(par, .5, labels = pow, col = "#666666", srt = 90, adj=c(.5,.5))
}

curve(pbinom(crit - 1, N, x, lower.tail = FALSE), par()$usr[1], par()$usr[2], 1024,
      lwd = 2, xpd = TRUE, add = TRUE)


### Divide



curve(pbinom(crit - 1, N, x, lower.tail = FALSE), 0, 1, 1024,
      ylim = c(0, 1),
      axes=FALSE, 
      ylab = substitute(paste("Probability that ",X>=c), list(c = crit)),
      xlab = "True binomial probability parameter p",
      lwd = 2, xpd = TRUE, ty = 'l', col = rgb(0,0,0,.4))
axis(1)
axis(2)

#curve(pbinom(crit, N, x, lower.tail = TRUE), par()$usr[1], par()$usr[2], 1024,
#      lwd = 2, xpd = TRUE, add = TRUE, col = "#6666ffdd")
#mtext(substitute(paste("Probability that ",X<=c), list(c = crit)),
#      side = 4, las = 0, line = par()$mgp[1], col = "#6666ff")
#axis(4, las = 1)

rect(par()$usr[1], par()$usr[3], 
     find_val_0(.01), par()$usr[4],
     border = NA, col = viridis::magma(10, alpha = .5)[3])

text(find_val_0(.01)/2, .9,
     "“small”", cex = 1.3, adj = c(.5,1))

rect(find_val_0(.99), par()$usr[3], 
     find_val_0(.01), par()$usr[4],
     border = NA, col = viridis::magma(10, alpha = .5)[7])

text(find_val_0(.5), .9,
     "“large”\n", cex = 1.3, adj = c(.5,1))


rect(par()$usr[2], par()$usr[3], 
     find_val_0(.99), par()$usr[4],
     border = NA, col = viridis::magma(10, alpha = .5)[8])

text(find_val_0(.99) + (1 - find_val_0(.99))/2, .9,
     "Reliably\ndetectably\n“large”", cex = 1.3, adj = c(.5,1))




#### Upper

curve(pbinom(crit, N, x, lower.tail = TRUE), x_lim[1], x_lim[2], 1024,
      ylim = c(0, 1),
      axes=FALSE, 
      ylab = "",
      xlab = "True binomial probability parameter p",
      lwd = 2, xpd = TRUE, ty = 'n')
axis(1)
mtext(substitute(paste("Probability that ",X<=c), list(c = crit)),
      side = 4, las = 0, line = par()$mgp[1], col = "#6666ff")
axis(4, las = 1)


for(i in seq_along(pow_vals_2)[-1]){
  par = find_val_1(pow_vals_2[(i-1):i])
  rect(par[1], par()$usr[3], par[2], par()$usr[4], col = col_vals_2[i],
       border = NA)
  
}

for(i in seq_along(pow_vals_2)[c(-1,-length(pow_vals_2))]){
  pow = pow_vals_2[i]
  par = find_val_1(pow)
  text(par, .5, labels = pow, col = "#6666ff", srt = 90, adj=c(.5,.5))
}

curve(pbinom(crit, N, x, lower.tail = TRUE), par()$usr[1], par()$usr[2], 1024,
      lwd = 2, xpd = TRUE, add = TRUE, col = "#6666ff", lty = 2)


#### Both

curve(pbinom(crit - 1, N, x, lower.tail = FALSE), x_lim[1], x_lim[2], 1024,
      ylim = c(0, 1),
      axes=FALSE, 
      ylab = substitute(paste("Probability that ",X>=c), list(c = crit)),
      xlab = "True binomial probability parameter p",
      lwd = 2, xpd = TRUE, ty = 'n')
axis(1)
axis(2)

for(i in seq_along(pow_vals)){
  pow = pow_vals[i]
  par = find_val2(pow)
  col = col_vals[i]
  rect(par[1], par()$usr[3], par[2], par()$usr[4], col = col,
       border = NA)
}

for(i in seq_along(pow_vals)[-1]){
  pow = pow_vals[i]
  par0 = find_val_0(pow)
  text(par0, .5, labels = pow, col = "#666666", srt = 90, adj=c(.5,.5))
  par1 = find_val_1(pow)
  text(par1, .5, labels = pow, col = "#6666ff", srt = 90, adj=c(.5,.5))
}


curve(pbinom(crit - 1, N, x, lower.tail = FALSE), par()$usr[1], par()$usr[2], 1024,
      lwd = 2, xpd = TRUE, add = TRUE)
curve(pbinom(crit, N, x, lower.tail = TRUE), par()$usr[1], par()$usr[2], 1024,
      lwd = 2, xpd = TRUE, add = TRUE, col = "#6666ff", lty = 2)

mtext(substitute(paste("Probability that ",X<=c), list(c = crit)),
      side = 4, las = 0, line = par()$mgp[1], col = "#6666ff")
axis(4, las = 1)


for(i in seq_along(pow_vals)[-1]){
  ci = ci_alpha(pow_vals[i])
  arrows(ci[1], pow_vals[i], ci[2], pow_vals[i], code = 3, angle = 90, 
         length = .03)
  }



```




```{r}

#### Normal

N = 100
crit = 2
x_lim = c(-.2,.6)

curve(pt(crit, N - 1, x * sqrt(N), lower.tail = FALSE), x_lim[1], x_lim[2], 1024,
      ylim = c(0, 1),
      axes=FALSE, 
      ylab = substitute(paste("Probability that ",t>=c), list(c = crit)),
      xlab = expression(paste("True (standardized) effect size ", delta)),
      lwd = 2, xpd = TRUE, ty = 'n')
axis(1)
axis(2)

for(i in seq_along(pow_vals_2)[-1]){
  par = find_val_t0(pow_vals_2[(i-1):i], N = N, crit = crit)
  par[par == Inf] = par()$usr[2]
  par[par == -Inf] = par()$usr[1]
  rect(par[1], par()$usr[3], par[2], par()$usr[4], col = col_vals_2[i],
       border = NA)
  
}

for(i in seq_along(pow_vals_2)[c(-1,-length(pow_vals_2))]){
  pow = pow_vals_2[i]
  par = find_val_t0(pow, N = N, crit = crit)
  text(par, .5, labels = pow, col = "#666666", srt = 90, adj=c(.5,.5))
}


curve(pt(crit, N - 1, x * sqrt(N), lower.tail = FALSE), par()$usr[1], par()$usr[2], 1024,
      lwd = 2, xpd = TRUE, add = TRUE)

ci_level = 90
ci = effsize::cohen.d(d = scale(rnorm(N)) + crit/sqrt(N), f = NA, 
                      conf.level = ci_level/100, noncentral = TRUE)$conf.int
arrows(ci[1], par()$usr[4], ci[2], par()$usr[4], xpd = TRUE,
       code = 3, angle = 90, length = .03)
text(crit/sqrt(N), par()$usr[4], paste0(ci_level,"% ","CI, if t=",crit), adj=c(.5,-0.2), xpd = TRUE)


### Normal 2

par(mar =  c(5.1, 4.1, 5.1, 4.1))

alpha = .025
N = c(100, 200, 400, 800)
crit = -qt(alpha, N - 1)
x_lim = c(-.2,.6)
col_vals = viridis::viridis(length(N))
ci_level = 95
vadj = .1 
vspc = .1


curve(pt(crit, N[1] - 1, x * sqrt(N[1]), lower.tail = FALSE), x_lim[1], x_lim[2], 1024,
      ylim = c(0, 1),
      axes=FALSE, 
      ylab = substitute(paste("Probability that ",t>=phantom(),"criterion"), list(c = crit)),
      xlab = expression(paste("True (standardized) effect size ", delta)),
      lwd = 2, xpd = TRUE, ty = 'n')
axis(1)
axis(2)

abline(h = c(alpha, 1-alpha), lty = 2, col = "gray")
abline(v = 0,lty = 2, col = "gray")

for(i in seq_along(N)){
  curve(pt(crit[i], N[i] - 1, x * sqrt(N[i]), lower.tail = FALSE), par()$usr[1], par()$usr[2], 1024,
        lwd = 2, xpd = TRUE, add = TRUE, col = col_vals[i])
  ci = effsize::cohen.d(d = scale(rnorm(N[i])) + crit[i]/sqrt(N[i]), f = NA, 
                        conf.level = ci_level/100, noncentral = TRUE)$conf.int
  arrows(ci[1], par()$usr[4] + vadj + vspc*(i-1), 
         ci[2], par()$usr[4] + vadj + vspc*(i-1), xpd = TRUE,
         code = 3, angle = 90, length = .03, col = col_vals[i])
  text(0, par()$usr[4] + vadj + vspc*(i-1), 
       paste0("N=",N[i]), adj = 1.1, xpd = TRUE)
}

for(i in seq_along(N)){
  points(find_val_t0(1-alpha, N=N[i], crit=crit[i]), 1-alpha,
         pch = 19, col = col_vals[i],xpd = TRUE)
}

points(0, alpha)

text(
  x = -.15,
  y = par()$usr[4] + vadj,
  label = paste0(ci_level,"% ","CI"),
  xpd = TRUE,
  srt = 90,
  adj = c(-.1,.5)
  )

legend("bottomright", legend = paste0("N=",N), col = col_vals, lwd = 2,
       lty = 1, bty = 'n')

```


```{r fig.height = 10}
ci_alpha = .05
sizes = c(1, 25)

par(mfrow = c(2,1))
par(family = "Roboto Condensed", las = 1, xaxs = 'i',
    yaxs = 'i', mar =  c(1, 4.1, 5.1, 4.1))

size = sizes[1]

plot(0,0,ylim = c(0,50), xlim = c(0,1), 
     xlab = "",
     ylab = "Observed failures", ty='n', axes=FALSE)
axis(2)
axis(3)
mtext("True neg. binomial prob. parameter p",3,line = par()$mgp[1])

x_vals = seq(floor(par()$usr[3]), ceiling(par()$usr[4]))
ci = ci_nbinom(y = x_vals, size = size, alpha = ci_alpha)
segments(ci[1,],x_vals,ci[2,],x_vals, lwd = 2, xpd = TRUE)

text(par()$usr[2], par()$usr[4]/2, paste0("size=",size), 
     cex = 1.3, adj = .9, xpd = TRUE)

size = sizes[2]

par(mar =  c(5.1, 4.1, 1, 4.1))

plot(0,0,ylim = c(0,50), xlim = c(0,1), 
     xlab = "True neg. binomial prob. parameter p",
     ylab = "Observed failures", ty='n', axes=FALSE)
axis(1)
axis(2)

x_vals = seq(floor(par()$usr[3]), ceiling(par()$usr[4]))
ci = ci_nbinom(y = x_vals, size = size, alpha = ci_alpha)
segments(ci[1,],x_vals,ci[2,],x_vals, lwd = 2, xpd = TRUE)

text(par()$usr[2], par()$usr[4]/2, paste0("size=",size), 
     cex = 1.3, adj = .9, xpd = TRUE)


```


```{r}
vadj = .1
vspc = .1

sizes = c(1,10,25,50)

p_null = 1/3

par(family = "Roboto Condensed", las = 1, xaxs = 'i',
    yaxs = 'i', mfrow = c(1,1), mar = c(5.1, 4.1, 5.1, 4.1))

crits = qnbinom(1-ci_alpha/2, sizes, p_null)
col_vals = viridis::viridis(length(crits))
  
curve(pnbinom(crits[1]-1, sizes[1], x, lower.tail = FALSE), 
      0, .5, 1024, ylim = c(0,1),
      xlab = "True neg. binomial prob. parameter p",
      ylab = substitute(paste("Probability that ",X>=phantom(),"criterion")),
      axes=FALSE, xpd = TRUE, lwd = 2, typ='n')
abline(v = p_null,col = "gray", lty = 2)
abline(h = c(ci_alpha/2, 1-ci_alpha/2),
       col = "gray", lty = 2)
axis(1)
axis(2)

for(i in seq_along(sizes)){
  curve(pnbinom(crits[i]-1, sizes[i], x, lower.tail = FALSE),
        par()$usr[1], par()$usr[2], 1024, 
        add = TRUE, xpd = TRUE, lwd = 2, col = col_vals[i])
  ci = ci_nbinom(crits[i], size = sizes[i],alpha = ci_alpha)
  arrows(ci[1], par()$usr[4] + vadj + vspc*(i-1), 
         ci[2], par()$usr[4] + vadj + vspc*(i-1), xpd = TRUE,
         code = 3, angle = 90, length = .03, col = col_vals[i])
  text(p_null, par()$usr[4] + vadj + vspc*(i-1), 
       paste0("size=",sizes[i]), adj = -.2, xpd = TRUE)

}

points(p_null, ci_alpha/2)

for(i in seq_along(N)){
  points(ci_nbinom(crits[i], size = sizes[i], alpha = ci_alpha)[1],
         1-ci_alpha/2,
         pch = 19, col = col_vals[i], xpd = TRUE)
}

legend("topright", legend = paste0("size=",sizes), col = col_vals,
       lwd = 2, bty = 'n')

text(
  x = par()$usr[2],
  y = par()$usr[4] + vadj,
  label = paste0((1-ci_alpha)*100,"% ","CI"),
  xpd = TRUE,
  srt = 90,
  adj = c(-.1,.5)
  )


```