Created
June 12, 2020 16:30
-
-
Save richarddmorey/1ee3ebe3216aa1f0a68b931c2dbfa95b to your computer and use it in GitHub Desktop.
Figures for "Power and precision" blog post
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
--- | |
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