Verringert sich die Bevölkerung der blauäugigen Marsmenschen?

Nehmen wir an, wir wollen die Hypothese testen, dass der Anteil der blauäugigen Marsmenschen im Laufe des 20. Jahrhunderts abgenommen hat. Leider schwankt die Marsbevölkerung stark, so dass jedes Jahrzehnt ein großer Unterschied in der Gesamtbevölkerung besteht [Update: Betrachten Sie die Marsbevölkerung als konstant bei einer Milliarde Marsmenschen. Die folgenden Daten sind Zufallsstichproben für jedes Jahr. Der Datensatz (der sich beim Schreiben zusammensetzt) ​​könnte ungefähr so ​​aussehen:

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

Die Analyse von Jahren, in denen die Marsbevölkerung unter 100 Jahre alt ist, ist statistisch eindeutig nicht so aussagekräftig wie bei einer Bevölkerung über 10 000, da wir im letzteren Fall einen größeren Datensatz haben. Dennoch möchten wir alle verfügbaren Daten verwenden, um unsere Hypothese mit einem herkömmlichen Signifikanzniveau von 95% zu überprüfen.

Wie Fahren wir fort? Gewichten wir die Wichtigkeit eines jeden Jahres nach der Größe der Stichprobe zu der Zeit?

Weitere Änderungen, um den Bedenken gerecht zu werden: Hier geht es darum, wie wir jeden Datensatz angemessen gewichten, wobei zu berücksichtigen ist, dass er so unterschiedliche Größen hat. Es gibt keine Stichprobenverschiebung, da die Daten zufällig ausgewählt werden.

$0$$1$

---

$n$$p$$k$

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

Jedes Jahrzehnt kennen wir und diese sind durch die Daten gegeben - aber wir kennen . Wir können es schätzen, indem wir annehmen, dass die logarithmischen Quoten , die entsprechen, linear von Jahr zu Jahr variieren (zumindest in guter Näherung). Dies bedeutet, dass wir davon ausgehen, dass es die Zahlen und so dass$n$$k$$p$$p$$\beta_0$$\beta_1$

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

Gleichermaßen

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

Wenn Sie dies in (1) einstecken, haben Sie die Möglichkeit, aus während eines bestimmten Jahres als zu beobachten$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}$$

Unter der Annahme, dass die Proben unabhängig in den Jahren usw. erhalten werden und die entsprechenden und von blauäugigen Probanden als und , ist die Wahrscheinlichkeit der Daten das Produkt der Wahrscheinlichkeiten der einzelnen Ergebnisse. Dieses Produkt ist (per Definition) die Wahrscheinlichkeit von . Wir können diese Parameter als die Werte schätzen , die die Wahrscheinlichkeit maximieren; Entsprechend maximieren sie die Log-Wahrscheinlichkeit$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}$$

erhalten aus .$(2)$

(Dies vereinfacht sich erheblich, wenn Logarithmenregeln verwendet werden. Dies ist ein Grund, die Zeit-Verhältnis-Beziehung in Form von Log-Quoten auszudrücken. Wenn alle Proportionen ungefähr zwischen und , gibt es nur einen geringen qualitativen Unterschied zwischen der Verwendung von Wahrscheinlichkeiten oder ihrem Log Gewinnchancen: Die angepasste Kurve ist linear bzw. nahezu linear.)$0.2$$0.8$$p$

$(3)$ ist ein binomiales verallgemeinertes lineares Modell . Es muss durch numerische Minimierung von angepasst werden . Das Verfahren in (am Ende dieses Beitrags gezeigt) gibt die Lösung$\Lambda$glmR

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

Die Daten in dieser Figur sind mit Scheiben aufgetragen, deren Flächen proportional zu den Stichprobengrößen sind. Die GLM-Passform ist krummlinig. Zum Vergleich ist in grau die Zeile dargestellt, die wir erhalten würden, wenn wir die in der Frage angezeigten Daten in einen gewöhnlichen Löser für kleinste Quadrate ausgeben würden. Beide Passungen werden trotz der damals geringen Stichprobengröße durch die größeren Anteile in früheren Jahren beeinflusst. Die GLM-Anpassung kann jedoch die Anteile in den größten Proben, die 1970 und 1980 erhalten wurden, besser approximieren. Die gepunktete blaue Linie wird nachstehend beschrieben.$(\text{Year},\text{Proportion})$

Durch Hinzufügen eines quadratischen Terms können wir die Anpassungsgüte testen. Es verbessert die GLM-Anpassung erheblich (obwohl der Unterschied optisch nicht groß ist) und liefert den Beweis, dass dieses Modell die Variation der Ergebnisse nicht gut beschreibt. Ein Blick auf die Darstellung zeigt, dass das Ergebnis im Jahr 1990 viel niedriger war als vom Modell vorhergesagt.

---

Ein alternativer, aber vergleichbarer Ansatz besteht darin, für jedes Jahr einzeln zu schätzen , möglicherweise als (obwohl andere Schätzer möglich sind). Eine lineare Regression der logarithmischen Quoten dieser Schätzungen gegenüber dem Jahr, gewichtet mit den Stichprobengrößen oder der Regression der gewichteten kleinsten Quadrate, ergibt$p$$t_i$$k_i / n_i$$n_i$

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

Die Standardfehler dieser Schätzungen betragen bzw. , was darauf hinweist, dass sich die WLS-Schätzungen nicht wesentlich vom Binomial-GLM unterscheiden. (Die Standardfehler des GLM sind jedoch erheblich kleiner: Es "weiß", dass diese Stichprobengrößen ziemlich groß sind, während die lineare Regression überhaupt nichts über die Stichprobengrößen "weiß": Es gibt nur eine Folge von zehn separaten Beobachtungen.) Beachten Sie, dass dies der Fall ist Eine Alternative ist möglicherweise nicht verfügbar, wenn oder , es sei denn, ein anderer Schätzer der Wahrscheinlichkeiten wird verwendet (der keine Werte von oder ).$15.55$$0.00787$$k_i=n_i$$k_i=0$$0$$1$

Schließlich könnten wir einfach eine gewichtete Regression der kleinsten Quadrate der Rohwahrscheinlichkeitsschätzungen gegenüber dem Jahr durchführen, die umgekehrt durch eine Schätzung der Stichprobenvarianz gewichtet wird. Die Varianz einer Binomialverteilung Variable , Wieder ausgedrückt als Verhältnis ist . Das kann aus einer Stichprobe als geschätzt werden$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}.$$

Das Ergebnis erscheint in der Abbildung als gepunktete blaue Linie. In diesem Fall scheint es einen Kompromiss zwischen den GLM- und OLS-Anpassungen zu geben.

---

Der folgende RCode führte die Analysen durch und erzeugte die Figur.

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