Skip to content

Instantly share code, notes, and snippets.

@richarddmorey

richarddmorey/planning_precision.Rmd

Created June 12, 2020 16:30
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save richarddmorey/1ee3ebe3216aa1f0a68b931c2dbfa95b to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save richarddmorey/1ee3ebe3216aa1f0a68b931c2dbfa95b to your computer and use it in GitHub Desktop.
Download ZIP
Figures for "Power and precision" blog post
Raw
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)
)
```
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment