Twój model zakłada, że sukces gniazda może być postrzegany jako hazard: Bóg przewraca załadowaną monetę z bokami oznaczonymi „sukces” i „porażka”. Wynik rzutu dla jednego gniazda jest niezależny od wyniku rzutu dla każdego innego gniazda.
Ptaki mają jednak coś dla siebie: moneta może w niektórych temperaturach znacznie sprzyjać sukcesowi w porównaniu do innych. Kiedy więc masz szansę obserwować gniazda w danej temperaturze, liczba sukcesów jest równa liczbie udanych rzutów tej samej monety - tej dla tej temperatury. Odpowiedni rozkład dwumianowy opisuje szanse na sukces. Oznacza to, że ustala prawdopodobieństwo zerowego sukcesu, jednego, dwóch, ... itd. Poprzez liczbę gniazd.
Jednym rozsądnym oszacowaniem związku między temperaturą a tym, jak Bóg ładuje monety, jest proporcja sukcesów zaobserwowanych w tej temperaturze. Jest to oszacowanie maksymalnego prawdopodobieństwa (MLE).
71033/7.3/73
5,10,15,200,3,2,32,7,5,3
Górny rząd rysunku pokazuje MLE w każdej z czterech zaobserwowanych temperatur. Czerwona krzywa na panelu „Dopasuj” pokazuje, w jaki sposób ładowana jest moneta, w zależności od temperatury. Z założenia ślad ten przechodzi przez każdy z punktów danych. (To, co robi w temperaturach pośrednich, jest nieznane; z grubsza połączyłem wartości, aby podkreślić ten punkt).
Ten „nasycony” model nie jest zbyt przydatny, właśnie dlatego, że nie daje nam podstaw do oszacowania, jak Bóg załaduje monety w temperaturach pośrednich. Aby to zrobić, musimy założyć, że istnieje jakaś krzywa „trendu”, która wiąże ładunki monet z temperaturą.
Dolny rząd figury wpisuje się w taki trend. Tendencja jest ograniczona w zakresie tego, co może zrobić: po wykreśleniu w odpowiednich współrzędnych („log odds”), jak pokazano w panelach „Logit Response” po lewej stronie, może ona podążać tylko po linii prostej. Każda taka prosta określa obciążenie monety we wszystkich temperaturach, jak pokazano za pomocą odpowiedniej zakrzywionej linii w panelach „Dopasuj”. To obciążenie z kolei określa rozkłady dwumianowe we wszystkich temperaturach. Dolny rząd przedstawia te rozkłady temperatur, w których zaobserwowano gniazda. (Czarne przerywane linie oznaczają oczekiwane wartości rozkładów, pomagając w ich dość dokładnej identyfikacji. Nie widać tych linii w górnym rzędzie rysunku, ponieważ pokrywają się one z czerwonymi segmentami.)
Teraz trzeba dokonać kompromisu: linia może przechodzić blisko niektórych punktów danych, tylko po to, aby oddalić się od innych. Powoduje to, że odpowiedni rozkład dwumianowy przypisuje niższe prawdopodobieństwo większości obserwowanych wartości niż wcześniej. Widać to wyraźnie przy 10 stopniach i 15 stopniach: prawdopodobieństwo zaobserwowanych wartości nie jest najwyższym możliwym prawdopodobieństwem, ani nie jest zbliżone do wartości przypisanych w górnym rzędzie.
Regresja logistyczna przesuwa się i porusza możliwymi liniami wokół (w układzie współrzędnych stosowanym przez panele „Logit Response”), przekształca ich wysokości w prawdopodobieństwa dwumianowe (panele „Fit”), ocenia szanse przypisane do obserwacji (cztery prawe panele ) i wybiera linię, która daje najlepszą kombinację tych szans.
Co jest „najlepsze”? Po prostu to, że łączne prawdopodobieństwo wszystkich danych jest tak duże, jak to możliwe. W ten sposób żadne pojedyncze prawdopodobieństwo (czerwone segmenty) nie może być naprawdę małe, ale zwykle większość prawdopodobieństw nie będzie tak wysoka, jak w modelu nasyconym.
Oto jedna iteracja wyszukiwania regresji logistycznej, w której linia została obrócona w dół:
1015stopnie, ale okropne zadanie dopasowania innych danych. (Przy 5 i 20 stopniach dwumianowe prawdopodobieństwa przypisane do danych są tak małe, że nawet nie widać czerwonych segmentów.) Ogólnie rzecz biorąc, jest to znacznie gorsze dopasowanie niż te pokazane na pierwszym rysunku.
Mam nadzieję, że ta dyskusja pomogła ci w opracowaniu mentalnego obrazu prawdopodobieństw dwumianowych zmieniających się wraz ze zmianą linii, przy jednoczesnym zachowaniu tych samych danych. Dopasowanie linii przez regresję logistyczną próbuje sprawić, że te czerwone słupki ogólnie są tak wysokie, jak to możliwe. Tak więc związek między regresją logistyczną a rodziną rozkładów dwumianowych jest głęboki i bliski.
Dodatek: R
kod do tworzenia cyfr
#
# Create example data.
#
X <- data.frame(temperature=c(5,10,15,20),
nests=c(2,7,5,3),
successes=c(0,3,2,3))
#
# A function to plot a Binomial(n,p) distribution and highlight the value `k0`.
#
plot.binom <- function(n, p, k0, highlight="#f02020", ...) {
plot(0:n, dbinom(0:n, n, p), type="h", yaxt="n",
xlab="Trials", ylab="Probability", ...)
abline(v = p*n, lty=3, lwd=2)
if(!missing(k0)) lines(rep(k0,2), c(0, dbinom(k0,n,p)), lwd=2, col=highlight)
}
#
# A function to convert from probability to log odds.
#
logit <- function(p) log(p) - log(1-p)
#
# Fit a saturated model, then the intended model.
#
# Ordinarily the formula for the saturated model would be in the form
# `... ~ factor(temperature)`, but the following method makes it possible to
# plot the predicted values in a visually effective way.
#
fit.0 <- glm(cbind(successes, nests-successes) ~ factor(round(temperature/5)),
data=X, family=binomial)
summary(fit.0)
fit <- glm(cbind(successes, nests-successes) ~ temperature,
data=X, family=binomial)
summary(fit)
#
# Plot both fits, one per row.
#
lfits <- list(fit.0, fit)
par.old <- par(mfrow=c(length(lfits), nrow(X)+2))
for (fit in lfits) {
#
# Construct arrays of plotting points.
#
X$p.hat <- predict(fit, type="response")
Y <- data.frame(temperature = seq(min(X$temperature), max(X$temperature),
length.out=101))
Y$p.hat <- predict(fit, type="response", newdata=Y) # Probability
Y$lambda.hat <- predict(fit, type="link", newdata=Y) # Log odds
#
# Plot the fit in terms of log odds.
#
with(Y, plot(temperature, lambda.hat, type="n",
yaxt="n", bty="n", main="Logit Response",
ylab=expression(hat(lambda))))
if (isTRUE(diff(range(Y$lambda.hat)) < 6)) {
# Draw gridlines and y-axis labels
p <- c( .10, .25, .5, .75, .9)
q <- logit(p)
suppressWarnings(rug(q, side=2))
abline(h=q, col="#d0d0d0")
mtext(signif(p, 2), at=q, side=2, cex=0.6)
}
with(Y, lines(temperature, lambda.hat, lwd=2, col="#f02020"))
#
# Plot the data and the fit in terms of probability.
#
with(X, plot(temperature, successes/nests, ylim=0:1,
cex=sqrt(nests), pch=21, bg="Gray",
main="Fit"))
with(Y, lines(temperature, p.hat, col="#f02020", lwd=2))
#
# Plot the Binomial distributions associated with each row of the data.
#
apply(X, 1, function(x) plot.binom(x[2], x[4], x[3], bty="n", lwd=2, col="Gray",
main=paste(x[1], "Degrees")))
}
par(mfrow=par.old)