Computing Statistical Power
[This article was first published on Coffee and Econometrics in the Morning, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Today’s task is to compute a statistical power function in R. For a good definition of power, go here. Thinking about the power of a testing procedure is an often-overlooked aspect of hypothesis testing.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
I also explain what power (in a statistical sense) is in the video, as well as use my shadenorm() function to show what power looks like on the classical hypothesis testing graph. Along the way, I demonstrate how writing a function can make your life easier. Here’s the video:
Here is the script file I walk through in the video (copy into notepad for a better view):
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
## ----------------------- ## | |
## Computing and Graphing ## | |
## Statistical Power ## | |
## ----------------------- ## | |
## Parameters of the Problem ## | |
mu0 = 6 ## Null-hypothesized value | |
alpha = 0.05 ## Significance Level | |
s.sq = 16 ## Variance from pilot study | |
n1 = 64 ## Sample Size | |
alt = 6.5 | |
se = sqrt(s.sq/n1) | |
## -------------------------- ## | |
## What is Statistical Power? ## | |
## -------------------------- ## | |
cuts = c(1-alpha/2, alpha/2) | |
crits = qnorm(cuts, mu0, sqrt(s.sq/n1)) ## Critical Values Based on Null | |
shadenorm(mu=mu0, sig = se, outside=crits) ## My own function | |
shadenorm(mu = 6.5, sig = se, lines=TRUE, outside=crits, col="blue") | |
## ------------------------------ ## | |
## Write a Function to Compute it ## | |
## ------------------------------ ## | |
power = function(theta, mu, var, n, alpha=0.05){ | |
crit.l = qnorm(alpha/2, mu, sqrt(var/ n)) ## Critical Value Based on Null | |
crit.h = qnorm(1-alpha/2, mu, sqrt(var/ n)) ## Critical Value Based on Null | |
pr.high = pnorm(crit.h, theta, sd = sqrt(var/ n),lower.tail=FALSE) ## Prob Reject High | |
pr.low = pnorm(crit.l, theta, sd = sqrt(var/ n)) ## Prob Reject Low | |
pow = pr.low+pr.high | |
pow | |
} | |
power(thet=6.5, mu=6, var= 16, n = 64) | |
thet = seq(3,9, by=0.01) | |
pow = power(thet, mu0, s.sq, n1) | |
plot(thet, pow, type="l", ylim=c(0,1), main="My Power Plot") | |
n2 = 256 ## New Larger Sample Size | |
pow2 = power(thet, mu0, s.sq,n2) | |
lines(thet, pow2, lty="dashed", col="blue") | |
abline(h=0.05, lty="dotted") |
My script file uses a new version of my shadenorm() command that is a little easier to use (only change: now, there are no justbelow or justabove arguments. If you just specify below, that’s justbelow=TRUE…).
This is available here:
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
## -------------------------- ## | |
## An R Function that Shades ## | |
## under a normal density. ## | |
## ## | |
## This is a convenience ## | |
## function for polygon() ## | |
## -------------------------- ## | |
shadenorm = function(below=NULL, above=NULL, pcts = c(0.025,0.975), mu=0, sig=1, numpts = 500, color = "gray", dens = 40, | |
lines=FALSE,between=NULL,outside=NULL){ | |
if(is.null(between)){ | |
bothnull = is.null(below) & is.null(above) | |
if(bothnull==TRUE){ | |
below = ifelse(is.null(below), qnorm(pcts[1],mu,sig), below) | |
above = ifelse(is.null(above), qnorm(pcts[2],mu,sig), above) | |
} | |
} | |
if(is.null(outside)==FALSE){ | |
if(is.numeric(outside)==FALSE){if(outside==TRUE){outside=qnorm(pcts,mu,sig)}} | |
below = min(outside) | |
above = max(outside) | |
} | |
lowlim = mu - 4*sig | |
uplim = mu + 4*sig | |
x.grid = seq(lowlim,uplim, length= numpts) | |
dens.all = dnorm(x.grid,mean=mu, sd = sig) | |
if(lines==FALSE){ | |
plot(x.grid, dens.all, type="l", xlab="X", ylab="Density") | |
} | |
if(lines==TRUE){ | |
lines(x.grid,dens.all) | |
} | |
if(is.null(below)==FALSE){ | |
x.below = x.grid[x.grid<below] | |
dens.below = dens.all[x.grid<below] | |
polygon(c(x.below,rev(x.below)),c(rep(0,length(x.below)),rev(dens.below)),col=color,density=dens) | |
} | |
if(is.null(above)==FALSE){ | |
x.above = x.grid[x.grid>above] | |
dens.above = dens.all[x.grid>above] | |
polygon(c(x.above,rev(x.above)),c(rep(0,length(x.above)),rev(dens.above)),col=color,density=dens) | |
} | |
if(is.null(between)==FALSE){ | |
if(is.numeric(between)==FALSE){if(between==TRUE){between=qnorm(pcts,mu,sig)}} | |
from = min(between) | |
to = max(between) | |
x.between = x.grid[x.grid>from&x.grid<to] | |
dens.between = dens.all[x.grid>from&x.grid<to] | |
polygon(c(x.between,rev(x.between)),c(rep(0,length(x.between)),rev(dens.between)),col=color,density=dens) | |
} | |
} |
To leave a comment for the author, please follow the link and comment on their blog: Coffee and Econometrics in the Morning.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.