¿Está disminuyendo la población de marcianos de ojos azules?

Supongamos que queremos probar la hipótesis de que la proporción de marcianos de ojos azules ha disminuido a lo largo del siglo XX. Desafortunadamente, la población marciana fluctúa mucho, por lo que cada década hay una gran diferencia en la población total [actualización: considere que la población marciana es constante en mil millones de marcianos. Los datos a continuación son muestras aleatorias en cada año]. El conjunto de datos (compuesto mientras escribo esto) podría verse así:

Year | Total martian population | Blue-eyed martians | Proportion
1910 | 400                      | 250                | 0.625
1920 | 2000                     | 1000               | 0.500
1930 | 70                       | 40                 | 0.571
1940 | 30                       | 14                 | 0.467
1950 | 10                       | 4                  | 0.400
1960 | 140                      | 52                 | 0.371
1970 | 50 000                   | 15 400             | 0.308
1980 | 70 000                   | 22 000             | 0.314
1990 | 1500                     | 80                 | 0.053
2000 | 5000                     | 800                | 0.160

Analizar años cuando la población marciana tiene menos de 100 años claramente no es tan significativo estadísticamente como cuando la población supera los 10 000, ya que en este último caso tenemos un conjunto de datos más amplio. Aún así, nos gustaría utilizar todos los datos disponibles para verificar nuestra hipótesis con un nivel de significación convencional del 95%.

¿Como procedemos? ¿Consideramos la importancia de cada año de acuerdo con el tamaño de la muestra en ese momento?

Más edición para adaptarse a las preocupaciones: la preocupación aquí es cómo ponderamos adecuadamente cada conjunto de datos teniendo en cuenta que son de tamaños tan diferentes. No hay sesgo de muestra ya que los datos se seleccionan aleatoriamente.

Esta respuesta describe tres formas de manejar adecuadamente los diferentes tamaños de muestra: un modelo lineal generalizado y dos regresiones de mínimos cuadrados ordinarios ponderados. En este caso, los tres funcionan bien. En general, cuando algunas proporciones están cerca de o , el GLM es mejor.$0$$1$

---

Debido a que los tamaños de muestra son tan pequeños en comparación con las poblaciones (menos del diez por ciento de ellos), en una aproximación excelente, la distribución de resultados de ojos azules y no azules en una muestra de tamaño es Binomial (porque las muestras son aleatorio). El otro parámetro binomial, , es la proporción verdadera (pero desconocida) de sujetos de ojos azules en la población. Por lo tanto, la posibilidad de observar personas de ojos azules es$n$$p$$k$

$$\binom{n}{k}p^k(1-p)^{n-k}.\tag{1}$$

Cada década sabemos y son los datos, pero no sabemos . Podemos estimarlo suponiendo que las probabilidades de registro correspondientes a varían por año linealmente (al menos para una buena aproximación). Esto significa que asumimos que hay números y tales que$n$$k$$p$$p$$\beta_0$$\beta_1$

$$\log(p) - \log(1-p) = \beta_0 + \beta_1 \times \text{Year}.$$

Equivalentemente

$$p = \frac{1}{1 + e^{-\beta_0-\beta_1\text{Year}}};\ 1-p = \frac{ e^{-\beta_0-\beta_1\text{Year}}}{1 + e^{-\beta_0-\beta_1\text{Year}}}.$$

Conectar esto a (1) da la posibilidad de observar de durante un año dado como$k$$n$$t$

$$\binom{n}{k} \frac{e^{-(\beta_0+\beta_1t)(n-k)}}{\left(1 + e^{-(\beta_0+\beta_1t)}\right)^n}.\tag{2}$$

Suponiendo que las muestras se obtienen de forma independiente en los años etc. y escribiendo los tamaños de muestra correspondientes y los recuentos de sujetos de ojos azules como y , la probabilidad de los datos es el producto de las probabilidades de los resultados individuales. Este producto es (por definición) la probabilidad de . Podemos estimar estos parámetros como los valores que maximizan la probabilidad; de manera equivalente, maximizan la probabilidad de registro$t_1, t_2,$ $n_i$$k_i$$(\beta_0, \beta_1)$$(\hat\beta_0, \hat\beta_1)$

$$\Lambda(\beta_0,\beta_1) = \sum_t \log\left(\binom{n}{k} \frac{e^{-(\beta_0+\beta_1t)(n-k)}}{\left(1 + e^{-(\beta_0+\beta_1t)}\right)^n}\right)\tag{3}$$

obtenido de .$(2)$

(Esto se simplifica considerablemente, usando reglas de logaritmos, que es una razón para expresar la relación de proporción de tiempo en términos de probabilidades de registro. Cuando todas las proporciones están entre y , aproximadamente, hay poca diferencia cualitativa entre usar probabilidades o su registro probabilidades: la curva ajustada será lineal o cercana a la lineal, respectivamente.)$0.2$$0.8$$p$

$(3)$ es un modelo lineal generalizado binomial . Debe ajustarse minimizando numéricamente . El procedimiento en (que se muestra al final de esta publicación) da la solución$\Lambda$glmR

$$(\hat\beta_0, \hat\beta_1)_\text{GLM} = (31.498711, -0.0163568).$$

Los datos en esta figura se trazan con discos cuyas áreas son proporcionales a los tamaños de muestra. El ajuste GLM es curvilíneo. En comparación, en gris, se muestra la línea que obtendríamos simplemente volcando los datos que se muestran en la pregunta en un solucionador de mínimos cuadrados ordinarios. Ambos ajustes están influenciados por las mayores proporciones en años anteriores, a pesar de los pequeños tamaños de muestra en ese momento. Sin embargo, el ajuste GLM hace un mejor trabajo al aproximar las proporciones en las muestras más grandes obtenidas en 1970 y 1980. La línea azul punteada se describe a continuación.$(\text{Year},\text{Proportion})$

Al agregar un término cuadrático, podemos probar la bondad del ajuste. Mejora significativamente el ajuste GLM (aunque visualmente la diferencia no es grande), proporcionando evidencia de que este modelo no describe bien la variación en los resultados. Mirar la gráfica indica que el resultado en 1990 fue mucho más bajo de lo que predice el modelo.

---

Un enfoque alternativo, pero comparable, es estimar individualmente para cada año , tal vez como (aunque son posibles otros estimadores). Una regresión lineal de las probabilidades de registro de estas estimaciones con respecto al año, ponderada por los tamaños de muestra , o regresión de mínimos cuadrados ponderados, arroja$p$$t_i$$k_i / n_i$$n_i$

$$(\hat\beta_0, \hat\beta_1)_\text{WLS} = (36.12744, -0.018706).$$

Los errores estándar de estas estimaciones son y , respectivamente, lo que indica que las estimaciones de WLS no son significativamente diferentes del GLM binomial. (Sin embargo, los errores estándar de GLM son considerablemente menores: "sabe" que estos tamaños de muestra son bastante grandes, mientras que la regresión lineal "no sabe" nada sobre los tamaños de muestra: solo tiene una secuencia de diez observaciones separadas). Tenga en cuenta que esto La alternativa podría no estar disponible si o , a menos que se use un estimador diferente de las probabilidades (que no produce valores de o ).$15.55$$0.00787$$k_i=n_i$$k_i=0$$0$$1$

Finalmente, podríamos simplemente realizar una regresión ponderada de mínimos cuadrados de las estimaciones de probabilidad bruta contra el año, inversamente ponderada por una estimación de la varianza de la muestra. La varianza de una variable binomial , reexpresada como una proporción es . Eso puede estimarse a partir de una muestra como$k/n$$(n,p)$$X$$X/n$$p(1-p)/n$

$$p(1-p)n \approx \frac{k}{n}\frac{n-k}{n}/n = \frac{k(n-k)}{n^3}.$$

Su resultado aparece en la figura como una línea azul punteada. En este caso, parece haber un compromiso entre los ajustes GLM y OLS.

---

El siguiente Rcódigo realizó los análisis y produjo la figura.

year <- seq(1910, 2000, by=10)
total <- c(40, 200, 7, 3, 1, 14, 5000, 7000, 150, 500) * 10
blue <- c(250, 1000, 40, 14, 4, 52, 15400, 22000, 80, 800)
X <- data.frame(Year=year, Success=blue, Failure=total-blue,
                Proportion=blue/total, Total=total)
#
# GLM
#
fit <- glm(cbind(Success, Failure) ~ Year, X, family="binomial")
summary(fit)
#
# WLS of the log odds (an alternative)
#
fit.WLS <- lm(log(Success/Failure) ~ Year, X, weights=Total)
summary(fit.WLS)
#
# Plot the results.
#
X.more <- data.frame(Year=1901:2010)
X.more$Prediction <- predict(fit, X.more, type="response")
plot(X$Year, X$Proportion, ylim=0:1,
     type="p", pch=21, bg="Red", cex=sqrt(X$Total/2000),
     xlab="Year", ylab="Proportion",
     main="GLM and OLS Fits", sub="GLM: solid line; OLS: dotted line")
lines(X.more, lwd=2)
abline(lm(Proportion ~ Year, X), 
       lty=3, lwd=3, col="Gray") #The OLS fit
abline(lm(Proportion ~ Year, X, weights=Total^3/(Success*Failure)), 
       lty=3, lwd=3, col="Blue") #The weighted OLS fit to the proportions