A população de marcianos de olhos azuis está diminuindo?

Suponha que desejamos testar a hipótese de que a proporção de marcianos de olhos azuis vem diminuindo ao longo do século XX. Infelizmente, a população marciana flutua muito; portanto, a cada década, há uma grande diferença na população total [atualização: considere a população marciana constante em um bilhão de marcianos. Os dados abaixo são amostras aleatórias a cada ano]. O conjunto de dados (composto no momento em que escrevo isso) pode ser algo como isto:

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

Analisar os anos em que a população marciana tem menos de 100 anos não é claramente tão significativo estatisticamente como quando a população está acima de 10.000, pois neste último caso, temos um conjunto de dados maior. Ainda assim, gostaríamos de usar todos os dados disponíveis para verificar nossa hipótese com um nível de significância convencional de 95%.

Como procedemos? Ponderamos a importância de cada ano de acordo com o tamanho da amostra da época?

Mais edições para atender às preocupações: a preocupação aqui é como ponderamos adequadamente cada conjunto de dados, tendo em mente que eles têm tamanhos tão diferentes. Não há viés de amostra, pois os dados são selecionados aleatoriamente.

Esta resposta descreve três maneiras de lidar com os diferentes tamanhos de amostra adequadamente: um Modelo Linear Generalizado e duas regressões de Mínimos Quadrados Ordinários ponderados. Nesse caso, os três funcionam bem. Em geral, quando algumas proporções estão próximas de ou , o GLM é melhor.$0$$1$

---

Como os tamanhos das amostras são muito pequenos em comparação com as populações (menos de dez por cento deles), para uma excelente aproximação, a distribuição dos resultados de olhos azuis e não de olhos azuis em uma amostra de tamanho é Binomial (porque as amostras são aleatória). O outro parâmetro binomial, , é a proporção verdadeira (mas desconhecida) de indivíduos de olhos azuis na população. Assim, a chance de observar pessoas de olhos azuis ép $n$$p$$k$

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

Cada década sabemos e --those são dadas pelos dados - mas não sabemos . Podemos estimar isso assumindo que as chances do log correspondente a variem linearmente por ano (pelo menos para uma boa aproximação). Isso significa que assumimos que existem números e tais 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}}}.$$

Ao conectar isso em (1), é possível observar de durante um determinado ano 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}$$

Assumindo que as amostras são independentemente obtido no ano etc e a escrever os correspondentes tamanhos da amostra e as contagens de sujeitos de olhos azuis como e , a probabilidade dos dados é o produto das probabilidades dos resultados individuais. Este produto é (por definição) a probabilidade de . Podemos estimar esses parâmetros como os valores que maximizam a probabilidade; equivalentemente, eles maximizam a probabilidade do log$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}$$

obtido de .$(2)$

(Isso simplifica consideravelmente o uso de regras de logaritmos, que é um dos motivos para expressar a relação de proporção do tempo em termos de probabilidades de log. Quando todas as proporções estão entre e , aproximadamente, há pouca diferença qualitativa entre o uso de probabilidades ou seu log. odds: a curva ajustada será linear ou próxima de linear, respectivamente.)$0.2$$0.8$$p$

$(3)$ é um modelo linear generalizado binomial . Ele deve ser ajustado minimizando numericamente . O procedimento em (mostrado no final deste post) fornece a solução$\Lambda$glmR

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

Os dados nesta figura são plotados com discos cujas áreas são proporcionais aos tamanhos da amostra. O ajuste do GLM é curvilíneo. Mostrada para comparação, em cinza, é a linha que obteríamos apenas despejando os dados mostrados na pergunta em um solucionador de Mínimos Quadrados Ordinários. Ambos os ajustes são influenciados pelas maiores proporções nos anos anteriores, apesar dos pequenos tamanhos de amostra. No entanto, o ajuste GLM faz um trabalho melhor na aproximação das proporções nas maiores amostras obtidas em 1970 e 1980. A linha azul pontilhada é descrita abaixo.$(\text{Year},\text{Proportion})$

Ao adicionar um termo quadrático, podemos testar a qualidade do ajuste. Melhora significativamente o ajuste do GLM (embora visualmente a diferença não seja grande), fornecendo evidências de que este modelo não descreve bem a variação nos resultados. Observar o gráfico indica que o resultado em 1990 foi muito menor do que o modelo prevê.

---

Uma abordagem alternativa, mas comparável, é estimar individualmente para cada ano , talvez como (embora outros estimadores sejam possíveis). Uma regressão linear das chances logarítmicas dessas estimativas em relação ao ano, ponderada pelo tamanho da amostra ou regressão dos mínimos quadrados ponderados, gera$p$$t_i$$k_i / n_i$$n_i$

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

Os erros padrão dessas estimativas são e , respectivamente, indicando que as estimativas do WLS não são significativamente diferentes do GLM binomial. (Porém, os erros padrão do GLM são consideravelmente menores: "sabe" que esses tamanhos de amostra são bastante grandes, enquanto a regressão linear "não sabe" nada sobre o tamanho das amostras: possui apenas uma sequência de dez observações separadas.) Observe que isso alternativa pode não estar disponível se ou , a menos que um estimador diferente de probabilidades seja usado (o que não produz valores de ou ).$15.55$$0.00787$$k_i=n_i$$k_i=0$$0$$1$

Finalmente, podemos simplesmente executar uma regressão ponderada dos mínimos quadrados das estimativas brutas de probabilidade relação ao ano, inversamente ponderadas por uma estimativa da variação da amostra. A variação de uma variável binomial , re-expressa como uma proporção é . Isso pode ser estimado a partir de uma amostra 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}.$$

Seu resultado aparece na figura como uma linha azul pontilhada. Nesse caso, parece haver comprometimento entre os ajustes GLM e OLS.

---

O Rcódigo a seguir executou as análises e produziu a 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