Curvas características de operación

Hasta hace muy pocos años la forma más sencilla de calcular el tamaño de la muestra necesario para cumplir con ciertos requisitos en el contraste de hipótesis era utilizar las curvas características de operación. Estas curvas dibujadas en los apéndices de ciertos libros se utilizan del mismo modo que las tablas de las distribuciones. Fijados unos parámetros se busca el valor deseado de la respuesta, en este caso el tamaño de la muestra.

Los gràficos del apéndice V del libro de Douglas C. Montgomery (2012) siempre me han fascinado. Pues bien, el objetivo de este artículo es reproducir esos gráficos con R. En realidad es un ejercicio puramente gráfico donde se utilizan profusamente los trucos gráficos y las funciones de bajo nivel de R. Ahora hay funciones en los programas estadísticos, también en R, para el cálculo de la potencia de un test que determinan el tamaño de la muestra con precisión y no es necesario recurrir a las curvas características de operación. Justamente para dibujar estos gráficos utilizaré una de estas funciones, en concreto la función pwr.anova.test() del paquete pwr.

En concreto vamos a centrarnos en el caso del ANOVA de un factor y, más exactamente, con \(a=4\) niveles. El problema radica en que los parámetros de la función de potencia y los de las curvas características de operación son distintos. Los detalles se pueden estudiar en el documento de Francesc Carmona (2015).

Esencialmente, para calcular el tamaño de la muestra se necesita fijar el nivel de significación \(\alpha\), el tamaño del efecto y la potencia deseada. En el caso de las curvas de operación los parámetros son: el nivel de significación, los grados de libertad \(\nu_1\) (que depende del número de niveles \(a\)) y \(\nu_2\) (que depende del número de réplicas \(n\)) y del parámetro \(\Phi\) de no centralidad de la \(F\). Las relaciones entre unos parámetros y otros son: \[\nu_1=a-1\] \[\nu_2=a(n-1)\] \[\Phi^2 = n f^2\] Con esta información ya estamos capacitados para escribir el código R que calcula las curvas características de operación de este caso. Observemos que en el mismo gráfico se representan las curvas para los dos niveles de significación habituales y que el eje de las Y está en escala logarítmica. Además hay un cierto trabajo de detalle para situar las flechas explicativas y las etiquetas en sus oportunos lugares.

El código

library(pwr)  # pwr.anova.test
library(plotrix)  # boxed.labels

nu1 <- 3  # a-1
a <- nu1 + 1
nu2 <- c(6, 7, 8, 9, 10, 12, 15, 20, 30, 60, 500)  # a*(n-1)
n <- nu2/a + 1
p <- c(seq(0.01, 0.08, 0.01), seq(0.1, 0.8, 0.1), 1)

oldpar <- par(mar = c(4, 4, 5, 2) + 0.1, cex = 0.8)  # bottom, left, top, right
plot(c(1, 6), c(1, 0), ylim = c(0, 1), type = "n", axes = F, xlab = "", 
    ylab = "Probability of accepting the hypothesis", cex.lab = 0.8)
alpha <- 0.05
phi <- seq(1, 4, 0.01)
for (i in n) {
    beta <- pwr.anova.test(k = a, n = i, f = phi/sqrt(i), sig.level = alpha, 
        power = NULL)$power
    lines(phi, log10((1 - beta) * 100)/2, lwd = 1.5, col = "blue")
}
alpha <- 0.01
phi <- seq(1, 5, 0.01)
for (i in n) {
    beta <- pwr.anova.test(k = a, n = i, f = phi/sqrt(i), sig.level = alpha, 
        power = NULL)$power
    lines(phi + 1, log10((1 - beta) * 100)/2, lwd = 1.5, col = "blue")
}
polygon(c(1, 6, 6, 1), c(-0.3, -0.3, 0, 0), col = "white", border = F)
# axis
axis(2, at = log10(p * 100)/2, label = format(p, nsmall = 2), las = 1)
axis(1, at = 1:6, label = c(1:3, expression(paste(Phi, " (for ", 
    alpha == 0.05, ")")), "", ""))
axis(3, at = 1:6, label = c(expression(paste(Phi, " (for ", alpha == 
    0.01, ")")), 1:5))
# lines
for (i in p) segments(1, log10(i * 100)/2, 6, log10(i * 100)/2)
for (i in 1:6) segments(i, 0, i, 1)
# title
title("Operating Characteristic Curves for the Fixed Effects Model Analysis of Variance", 
    line = 3)
# labels
text(1.5, 0.97, expression(nu[1] == 3))
text(1.5, 0.55, expression(alpha == 0.05))
text(4.5, 0.68, expression(alpha == 0.01))
# arrows
arrows(3, 0.33, 5.1, 0.33, length = 0.1, code = 3)  # 6
arrows(2.77, 0.365, 4.7, 0.365, length = 0.1, code = 3)  # 7
arrows(2.61, 0.4, 4.42, 0.4, length = 0.1, code = 3)  # 8
arrows(2.48, 0.435, 4.2, 0.435, length = 0.1, code = 3)  # 9
arrows(2.38, 0.47, 4.03, 0.47, length = 0.1, code = 3)  # 10
arrows(2.22, 0.505, 3.8, 0.505, length = 0.1, code = 3)  # 12
arrows(2.09, 0.54, 3.6, 0.54, length = 0.1, code = 3)  # 15
arrows(1.97, 0.575, 3.4, 0.575, length = 0.1, code = 3)  # 20
arrows(1.85, 0.61, 3.23, 0.61, length = 0.1, code = 3)  # 30
arrows(1.72, 0.645, 3.05, 0.645, length = 0.1, code = 3)  # 60
arrows(1.62, 0.68, 2.9, 0.68, length = 0.1, code = 3)  # infty
# labels
boxed.labels(3.3, 0.33, expression(nu[2] == 6), bg = "white", border = F, 
    cex = 0.8)
s <- 3.2
r <- 0.365
for (i in nu2[2:10]) {
    boxed.labels(s, r, i, bg = "white", border = F, cex = 0.8)
    s <- s - 0.07
    r <- r + 0.035
}
boxed.labels(s, 0.68, expression(infinity), bg = "white", border = F, 
    cex = 1)
par(oldpar)

Por último reproduzco la imagen original: