Skip to content

Instantly share code, notes, and snippets.

Created August 11, 2011 10:05
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save even4void/1139310 to your computer and use it in GitHub Desktop.
A simple animated demonstration of Power Analysis
mha <- 1 # mean under the alternative
es <- 2 # observed effect (deviation from mean under H0)
x <- seq(-6, 6, length=1000)
dh0 <- dnorm(x, 0, 1) <- function(es, mha, verbose=FALSE) {
dh1 <- dnorm(x, mha, 1)
plot.window(xlim=range(x), ylim=c(0,.6))
## show the two bell-shaped curves
draw.curve <- function(x, y, col) {
polygon(c(rev(x), x), c(rep(0, 1000), y), col=paste(col, "50", sep=""),
lines(x, y, lwd=1, col=col)
draw.curve(x, dh0, col="#2A4480")
draw.curve(x, dh1, col="#FFE973")
mtext(expression(mu[H0]), 1, -1.25, at=0)
mtext(expression(mu[H1]), 1, -1.25, at=mha)
## show the observed effect
abline(v=es, lty=2)
## highlight Type I and 1 - Type II error
polygon(c(rev(x[x>=es]), x[x>=es]), c(rep(0, sum(x>=es)), dh0[x>=es]),
col="#BF303075", border=NA)
lines(x[x>=es], dh0[x>=es], col="#BF3030", lwd=1)
polygon(c(rev(x[x>=es & dh0<dh1]), x[x>=es & dh0<dh1]),
c(rev(dh0[x>=es & dh0<dh1]), dh1[x>=es & dh0<dh1]),
col="#1D737375", border=NA)
lines(x[x>=es & dh0<dh1], dh1[x>=es & dh0<dh1], col="#1D7373", lwd=1)
alpha <- pnorm(es, 0, 1, lower.tail=FALSE)
beta <- pnorm(es, mha, 1)
labs.h0 <- "Type I Error (Alpha) "
labs.h1a <- "Type II Error (Beta) "
labs.h1b <- "1 - Type II Error (Power) "
if (verbose) {
labs.h0 <- paste(labs.h0, round(alpha*100, 1), "%", sep="")
labs.h1a <- paste(labs.h1a, round(beta*100, 1), "%", sep="")
labs.h1b <- paste(labs.h1b, round((1-beta)*100, 1), "%", sep="")
legend("topleft", c(labs.h0, labs.h1a, labs.h1b), pch=15,
col=c("#BF3030", "#FFE973", "#1D7373"), bty="n")
} <- function() {
update.display <- function(...) { <- slider(no=1)
s.mha <- slider(no=2)
s.verbose <- slider("show.numbers"), s.mha, s.verbose)
slider(update.display, sl.names=c("Effect", "Mean Under HA"),
sl.mins=c(0.5, 0.5), sl.maxs=c(4, 3),
sl.deltas=c(.05,.1), sl.defaults=c(es, mha),
function(...) {slider("show.numbers", obj.value=TRUE);
function(...) {slider("show.numbers", obj.value=FALSE);
but.names=c("On", "Off"))
slider("show.numbers", obj.value=FALSE)
}, mha, verbose=FALSE)
Copy link

A little interactive demo of power analysis. This is largely inspired from the Java applet found here: Here is a sample screenshot.

Copy link

maahiir commented Oct 30, 2022

I tried running the same code but I am getting error when I try to implement I get the following error

Error in tkconfigure(sc, variable = slider1) :
could not find function "tkconfigure"

Could you help me with this?

Copy link

Error in tkconfigure(sc, variable = slider1) :
could not find function "tkconfigure"

You likely need to install the tcltk2 package in R, e.g. install.packages("tcltk2", dependencies = TRUE).

Copy link

maahiir commented Oct 31, 2022

Thanks :)
This worked ::)


Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment