Tutkimusten otoskoon ja voimakkuuden laskeminen

Jari Haukka
15-03-2013

library(Hmisc)

## Loading required package: survival

## Loading required package: splines

## Hmisc library by Frank E Harrell Jr
## 
## Type library(help='Hmisc'), ?Overview, or ?Hmisc.Overview') to see overall
## documentation.
## 
## NOTE:Hmisc no longer redefines [.factor to drop unused levels when
## subsetting.  To get the old behavior of Hmisc type dropUnusedLevels().

## Attaching package: 'Hmisc'

## The following object(s) are masked from 'package:survival':
## 
## untangle.specials

## The following object(s) are masked from 'package:base':
## 
## format.pval, round.POSIXt, trunc.POSIXt, units

library(xtable)

## Attaching package: 'xtable'

## The following object(s) are masked from 'package:Hmisc':
## 
## label, label<-

Vaarasuhteet ja otoskoko (n)

# Plot power vs. n for various odds ratios (base prob.=.1)
n <- seq(10, 1000, by = 10)
OR <- seq(0.2, 0.9, by = 0.1)

plot(0, 0, xlim = range(n), ylim = c(0, 1), xlab = "n", ylab = "Power", type = "n", 
    main = "Base probability 0.1")
abline(h = c(0.8, 0.9), lty = 5)
for (or in OR) {
    lines(n, bpower(0.1, odds.ratio = or, n = n))
    text(350, bpower(0.1, odds.ratio = or, n = 350) - 0.02, format(or))
}

plot(0, 0, xlim = range(n), ylim = c(0.01, 1), xlab = "n", ylab = "Power", type = "n", 
    main = "Base probability 0.05", log = "")
abline(h = c(0.8, 0.9), lty = 5)
for (or in OR) {
    lines(n, bpower(0.05, odds.ratio = or, n = n))
    text(350, bpower(0.05, odds.ratio = or, n = 350) - 0.02, format(or))
}

Miten jakaa ryhmät optimaalisesti kun halutaan verrata kahta ryhmää

Alla olevassa kuviossa on esitetty khdelle tutkimuskoolle (100 ja 200) miten ryhmät olisi jaettava optimaalisesti jotta saataisiin suurin mahdollinen tutkimuksen voimakkuus. Oletetaan että toisen ryhmän “onnistumisen tn” on p1=0.1 ja lasketaan ryhmien osuudet eri toisen ryhmän (p2) “onnistumisen” arvoille.

• Brittain E, Schlesselman JJ (1982): Optimal allocation for the comparison of proportions. Biometrics 38:1003–9.

tmp.p1 <- 0.1
tmp.p2 <- seq(0.2, 0.9, by = 0.05)
tmp.max.power.100 <- sapply(tmp.p2, function(x) ballocation(p1 = tmp.p1, p2 = x, 
    n = 100)[4])
tmp.max.power.200 <- sapply(tmp.p2, function(x) ballocation(p1 = tmp.p1, p2 = x, 
    n = 200)[4])



tmp.power.100 <- rep(NA, length(tmp.max.power.100))
for (i in seq(tmp.power.100)) tmp.power.100[i] <- bpower(p1 = tmp.p1, p2 = tmp.p2[i], 
    n1 = 100 * tmp.max.power.100[i], n2 = 100 * (1 - tmp.max.power.100[i]), 
    alpha = 0.05)
tmp.power.200 <- rep(NA, length(tmp.max.power.200))
for (i in seq(tmp.power.200)) tmp.power.200[i] <- bpower(p1 = tmp.p1, p2 = tmp.p2[i], 
    n1 = 200 * tmp.max.power.200[i], n2 = 200 * (1 - tmp.max.power.200[i]), 
    alpha = 0.05)

plot(tmp.p2, tmp.max.power.200, type = "l", xlab = "p2", ylab = "Fraction for group 1 (p1=0.1)", 
    ylim = c(0.4, 0.65))
lines(tmp.p2, tmp.max.power.100, col = "green", lty = 2)
text(0.2, tmp.max.power.200[1], "n=200", adj = 0)
text(0.2, tmp.max.power.100[1], "n=100", adj = 0)
abline(h = 0.5, lty = 5)
title(main = "Allocation of groups to get maximum power")

tmp.ulos <- cbind(tmp.max.power.100, tmp.power.100, tmp.max.power.200, tmp.power.200)

dimnames(tmp.ulos) <- list(format(tmp.p2), c("Fraction 1/n=100", "Power/n=100", 
    "Fraction 1/n=200", "Power/n=200"))

print(xtable(tmp.ulos, digits = c(0, 2, 4, 2, 4), caption = "Allocation fraction and power for n=100 and n=200"), 
    type = "html")

Tapahtuma tutkimuksen otoskoko

# In this example, 4 plots are drawn on one page, one plot for each
# combination of noncompliance percentage.  Within a plot, the 5-year
# mortality % in the control group is on the x-axis, and separate curves
# are drawn for several % reductions in mortality with the intervention.
# The accrual period is 1.5y, with all patients followed at least 5y and
# some 6.5y.


# par(mfrow=c(2,2),oma=c(3,0,3,0)) par(mfrow=c(2,2),oma=c(1,0,1,0))
par(mfrow = c(1, 1))


morts <- seq(10, 25, length = 50)
red <- c(10, 15, 20, 25)

# Power vs Non-Adherence for Main Comparison alpha=.05, 2-tailed, Total
# N=14000

for (noncomp in c(0, 10, 15, -1)) {
    if (noncomp >= 0) 
        nc.i <- nc.c <- noncomp else {
        nc.i <- 25
        nc.c <- 15
    }
    z <- paste("Drop-in ", nc.c, "%, Non-adherence ", nc.i, "%", sep = "")
    plot(0, 0, xlim = range(morts), ylim = c(0, 1), xlab = "5-year Mortality in Control Patients (%)", 
        ylab = "Power", type = "n")
    title(z)
    cat(z, "\n")
    lty <- 0
    for (r in red) {
        lty <- lty + 1
        power <- morts
        i <- 0
        for (m in morts) {
            i <- i + 1
            power[i] <- cpower(5, 14000, m/100, r, 1.5, 5, nc.c, nc.i, pr = FALSE)
        }
        lines(morts, power, lty = lty)
    }
    if (noncomp == 0) 
        legend(18, 0.55, rev(paste(red, "% reduction", sep = "")), lty = 4:1, 
            bty = "n")
}

## Drop-in 0%, Non-adherence 0%

## Drop-in 10%, Non-adherence 10%

## Drop-in 15%, Non-adherence 15%

## Drop-in 15%, Non-adherence 25%

Jatkuva vaste (normaalijakauma)

tmp.N <- seq(10, 200, by = 10)

tmp.1 <- sapply(tmp.N, function(x) power.t.test(n = x, delta = 1, sd = 1)$power)
tmp.2 <- sapply(tmp.N, function(x) power.t.test(n = x, delta = 1, sd = 2)$power)
tmp.4 <- sapply(tmp.N, function(x) power.t.test(n = x, delta = 1, sd = 4)$power)
tmp.1.2 <- sapply(tmp.N, function(x) power.t.test(n = x, delta = 1, sd = 1/2)$power)



plot(range(tmp.N), range(c(tmp.1.2, tmp.4)), type = "n", xlab = "N", ylab = "Power", 
    log = "")
abline(h = 0.8, lty = 5, col = "green", v = tmp.N)
lines(tmp.N, tmp.1)
lines(tmp.N, tmp.2)
lines(tmp.N, tmp.4)
lines(tmp.N, tmp.1.2)

text(rep(10, 5), c(tmp.1[1], tmp.2[1], tmp.4[1], tmp.1.2[1]), c("1", "2", "4", 
    ".5"), adj = 0)

library(xtable)
tmp.ulos <- cbind(tmp.1, tmp.2, tmp.4, tmp.1.2)
dimnames(tmp.ulos) <- list(tmp.N, c("1", "2", "4", "0.5"))
print(xtable(tmp.ulos, digits = c(0, 4, 4, 4, 4), caption = "Power. N and SD as proportion of delta (100%, 200%, 400%, 50%))"), 
    type = "html")