La population de Martiens aux yeux bleus diminue-t-elle?

Supposons que nous voulons tester l'hypothèse selon laquelle la proportion de martiens aux yeux bleus a diminué tout au long du 20e siècle. Malheureusement, la population martienne fluctue fortement, donc chaque décennie il y a une grande différence dans la population totale [mise à jour: considérez la population martienne comme constante à un milliard de martiens. Les données ci-dessous sont des échantillons aléatoires chaque année]. L'ensemble de données (composé au moment où j'écris ceci) pourrait ressembler à ceci:

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

Analyser les années où la population martienne a moins de 100 ans n'est clairement pas aussi significatif statistiquement que lorsque la population est supérieure à 10 000, car dans ce dernier cas, nous avons un ensemble de données plus important. Néanmoins, nous aimerions utiliser toutes les données disponibles pour vérifier notre hypothèse avec un niveau de signification conventionnel de 95%.

Comment nous procédons? Pesons-nous l'importance de chaque année en fonction de la taille de l'échantillon à l'époque?

Plus de modifications pour répondre aux préoccupations: la préoccupation ici est de savoir comment pondérer adéquatement chaque ensemble de données en gardant à l'esprit qu'ils sont de tailles différentes. Il n'y a pas de biais d'échantillon car les données sont sélectionnées au hasard.

Cette réponse décrit trois façons de gérer les différentes tailles d'échantillon de manière appropriée: un modèle linéaire généralisé et deux régressions des moindres carrés ordinaires pondérés. Dans ce cas, les trois fonctionnent bien. En général, lorsque certaines proportions sont proches$0$ ou $1$, le GLM est meilleur.

---

Parce que les tailles d'échantillon sont si petites par rapport aux populations (moins de dix pour cent d'entre elles), pour une excellente approximation, la distribution des yeux bleus et non bleus donne un échantillon de taille $n$est binomial (car les échantillons sont aléatoires). L'autre paramètre binomial,$p$, est la proportion réelle (mais inconnue) de sujets aux yeux bleus dans la population. Ainsi, la chance d'observer$k$ les gens aux yeux bleus est

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

Chaque décennie, nous savons que et ceux-ci sont donnés par les données - mais nous ne savons pas . On peut l'estimer en supposant que la cote de log correspondant à varie linéairement d'une année à l'autre (au moins à une bonne approximation). Cela signifie que nous supposons qu'il existe des nombres et tels que$n$$k$$p$$p$$\beta_0$$\beta_1$

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

De manière équivalente,

$$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}}}.$$

Le brancher sur (1) donne la possibilité d'observer sur pendant une année donnée comme$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}$$

En supposant que les échantillons sont obtenus indépendamment aux années etc. et en écrivant les tailles et les nombres d'échantillons correspondants de sujets aux yeux bleus comme et , la probabilité des données est le produit des probabilités des résultats individuels. Ce produit est (par définition) la probabilité de . Nous pouvons estimer ces paramètres comme les valeurs qui maximisent la probabilité; de façon équivalente, ils maximisent la probabilité logarithmique$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}$$

obtenu à partir de .$(2)$

(Cela simplifie considérablement, en utilisant des règles de logarithmes, ce qui est une raison pour exprimer la relation temps-proportion en termes de cotes de log. Lorsque toutes les proportions sont comprises entre et , environ, il y a peu de différence qualitative entre l'utilisation des probabilités ou leur log cotes: la courbe ajustée sera respectivement linéaire ou proche de linéaire.)$0.2$$0.8$$p$

$(3)$ est un modèle linéaire généralisé binomial . Il doit être ajusté en minimisant numériquement . La procédure en (indiquée à la fin de ce post) donne la solution$\Lambda$glmR

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

Les données de cette figure sont tracées avec des disques dont les zones sont proportionnelles à la taille des échantillons. L'ajustement GLM est curviligne. Pour la comparaison, en gris, la ligne que nous obtiendrions simplement en vidant les données affichées dans la question dans un solveur Ordinaire des moindres carrés. Les deux ajustements sont influencés par les proportions plus importantes des années précédentes, malgré la petite taille des échantillons à l'époque. Cependant, l'ajustement GLM fait un meilleur travail d'approximation des proportions dans les plus grands échantillons obtenus en 1970 et 1980. La ligne bleue en pointillés est décrite ci-dessous.$(\text{Year},\text{Proportion})$

En ajoutant un terme quadratique, nous pouvons tester la qualité de l'ajustement. Il améliore considérablement l'ajustement GLM (bien que visuellement la différence ne soit pas grande), ce qui prouve que ce modèle ne décrit pas bien la variation des résultats. L'examen du graphique indique que le résultat en 1990 était bien inférieur à celui prévu par le modèle.

---

Une approche alternative, mais comparable, consiste à estimer individuellement pour chaque année , peut-être comme (bien que d'autres estimateurs soient possibles). Une régression linéaire des cotes logarithmiques de ces estimations par rapport à l'année, pondérée par les tailles d'échantillon , ou régression des moindres carrés pondérés, donne$p$$t_i$$k_i / n_i$$n_i$

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

Les erreurs-types de ces estimations sont respectivement de et , ce qui indique que les estimations WLS ne diffèrent pas significativement du GLM binomial. (Les erreurs standard du GLM sont cependant beaucoup plus petites: il "sait" que ces tailles d'échantillon sont assez grandes tandis que la régression linéaire "ne sait" rien du tout sur les tailles d'échantillon: il n'a qu'une séquence de dix observations distinctes.) Notez que cette l'alternative pourrait ne pas être disponible si ou , sauf si un estimateur différent des probabilités est utilisé (qui ne produit pas de valeurs de ou ).$15.55$$0.00787$$k_i=n_i$$k_i=0$$0$$1$

Enfin, nous pourrions simplement effectuer une régression des moindres carrés pondérés des estimations de probabilité brutes rapport à l'année, inversement pondérée par une estimation de la variance de l'échantillon. La variance d'une variable binomiale , ré-exprimée en proportion est . Cela peut être estimé à partir d'un échantillon comme$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}.$$

Son résultat apparaît sur la figure sous la forme d'une ligne bleue en pointillés. Dans ce cas, il semble qu'il y ait un compromis entre les ajustements GLM et OLS.

---

Le Rcode suivant a effectué les analyses et produit la figure.

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