3
$\begingroup$

I want the mathematical formula the statistical programs are using to compute alpha and beta. I typically use G-Power to estimate the necessary sample size for a given power (1-beta) for an a-priori analysis. The input it requires for this calculation is: Type of analysis (t-test, difference between two independent samples) the number of tails (2) alpha (0.05) power (0.8) effect size (0.2) = estimated necessary sample size (788) I have looked everywhere I know to find out what the function/formula is that is being used to compute the sample size value from these variables.

I know that:

  • alpha=type I error=FP/(FP+TN)
  • beta=type II error= FN/(TP+FN)
  • where FP=false positive
  • TP=true negative
  • FN=false negative
  • TN=true negative

But I don't know how to compute false positives or true negatives a-priori.

​Unless Stevens (2009) is wrong, power (the inverse of beta) is a function of alpha therefore β ∝ α (e.g. if, {α = .1, β = .37, power = .63}, {α = .05, β = .52, power = .48}, {α = .01, β = .78, power = .22}, see p.4). I am simply interested in finding the function that describes the relationship between alpha and beta/power. My question: Do you know where I can find the function to describe this relationship?

The reason? To actually understand what I am doing every time I run a statistical analysis.

when I perused the source code for an R package pwr I found:

What looks like the n is obtained numerically using a root finding algorithm computed with the highlighted part of the function below (uniroot looks like a general function to find where a supplied function equals 0).

> pwr.t.test
function (n = NULL, d = NULL, sig.level = 0.05, power = NULL, 
    type = c("two.sample", "one.sample", "paired"), alternative = c("two.sided", 
        "less", "greater")) 
{
    if (sum(sapply(list(n, d, power, sig.level), is.null)) != 
        1) 
        stop("exactly one of n, d, power, and sig.level must be NULL")
    if (!is.null(sig.level) && !is.numeric(sig.level) || any(0 > 
        sig.level | sig.level > 1)) 
        stop(sQuote("sig.level"), " must be numeric in [0, 1]")
    if (!is.null(power) && !is.numeric(power) || any(0 > power | 
        power > 1)) 
        stop(sQuote("power"), " must be numeric in [0, 1]")
    type <- match.arg(type)
    alternative <- match.arg(alternative)
    tsample <- switch(type, one.sample = 1, two.sample = 2, paired = 1)
    ttside <- switch(alternative, less = 1, two.sided = 2, greater = 3)
    tside <- switch(alternative, less = 1, two.sided = 2, greater = 1)
    if (tside == 2 && !is.null(d)) 
        d <- abs(d)
    if (ttside == 1) {
        p.body <- quote({
            nu <- (n - 1) * tsample
            pt(qt(sig.level/tside, nu, lower = TRUE), nu, ncp = sqrt(n/tsample) * 
                d, lower = TRUE)
        })
    }
    if (ttside == 2) {
        p.body <- quote({
            nu <- (n - 1) * tsample
            qu <- qt(sig.level/tside, nu, lower = FALSE)
            pt(qu, nu, ncp = sqrt(n/tsample) * d, lower = FALSE) + 
                pt(-qu, nu, ncp = sqrt(n/tsample) * d, lower = TRUE)
        })
    }
    if (ttside == 3) {
        p.body <- quote({
            nu <- (n - 1) * tsample
            pt(qt(sig.level/tside, nu, lower = FALSE), nu, ncp = sqrt(n/tsample) * 
                d, lower = FALSE)
        })
    }
    if (is.null(power)) 
        power <- eval(p.body)
    else if (is.null(n)) 
        n <- uniroot(function(n) eval(p.body) - power, c(2 + 
            1e-10, 1e+07))$root
        else if (is.null(d)) {
            if (ttside == 2) {
                d <- uniroot(function(d) eval(p.body) - power, c(1e-07, 
                    10))$root
            }
            if (ttside == 1) {
                d <- uniroot(function(d) eval(p.body) - power, c(-10, 
                    5))$root
            }
            if (ttside == 3) {
                d <- uniroot(function(d) eval(p.body) - power, c(-5, 
                    10))$root
            }
        }
        else if (is.null(sig.level)) 
            sig.level <- uniroot(function(sig.level) eval(p.body) - 
                power, c(1e-10, 1 - 1e-10))$root
    else stop("internal error")
    NOTE <- switch(type, paired = "n is number of *pairs*", two.sample = "n is number in *each* group", 
        NULL)
    METHOD <- paste(switch(type, one.sample = "One-sample", two.sample = "Two-sample", 
        paired = "Paired"), "t test power calculation")
    structure(list(n = n, d = d, sig.level = sig.level, power = power, 
        alternative = alternative, note = NOTE, method = METHOD), 
        class = "power.htest")
}
<environment: namespace:pwr>

I'm ​just ​not savvy enough with the notation to transform this into something prettier e.g. f(x)=...

$\endgroup$
  • $\begingroup$ For a general relationship between alpha and beta, did you take a look at stats.stackexchange.com/questions/59202/… ? $\endgroup$ – Roland Jan 28 '14 at 20:04
  • $\begingroup$ @roland Yes I've seen that post, i specifically want the mathematical function though. How would you suggest i change my question to make that more clear? $\endgroup$ – timothy.s.lau Jan 28 '14 at 20:06
1
$\begingroup$

Check out Introduction to sample size determination and power analysis for clinical trials by John M. Lachin. Equations (6) and (7) look like they will be useful to you. $$N = \frac{4\sigma^2(Z_\alpha+Z_\beta)^2}{\mu_1^2}$$ $$Z_\beta=\frac{\left|\mu_1\right|\sqrt{N}-2Z_\alpha\sigma}{2\sigma} $$

This is for equal sample sizes, but a form allowing for unequal sample sizes is found in the linked paper.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.