Podsumowanie: GLM są odpowiednie dzięki punktacji Fishera, która, jak zauważa Dimitriy V. Masterov, jest Newtonem-Raphsonem z oczekiwanym Hesjanem (tj. Zamiast szacowanej informacji używamy szacunku informacji Fishera). Jeśli używamy kanonicznej funkcji łącza, okazuje się, że obserwowany Hesjan jest równy spodziewanemu Hesji, więc NR i Fisher są w tym przypadku takie same. Tak czy inaczej, zobaczymy, że punktacja Fishera faktycznie pasuje do ważonego modelu liniowego metodą najmniejszych kwadratów, a oszacowania współczynników z tego zbiegają się * na maksymalnym prawdopodobieństwie regresji logistycznej. Oprócz zmniejszenia dopasowania regresji logistycznej do już rozwiązanego problemu, możemy również skorzystać z możliwości zastosowania diagnostyki regresji liniowej przy ostatecznym dopasowaniu WLS, aby dowiedzieć się o naszej regresji logistycznej.
Skupię się na regresji logistycznej, ale dla bardziej ogólnego spojrzenia na maksymalne prawdopodobieństwo w GLM polecam sekcję 15.3 tego rozdziału, która to omawia i wyprowadza IRLS w bardziej ogólnym otoczeniu (myślę, że pochodzi z Applied Johna Foxa Analiza regresji i uogólnione modele liniowe ).
∗ patrz komentarze na końcu
Funkcja prawdopodobieństwa i wyniku
Będziemy dopasowywać nasz GLM przez iterację czegoś w postaci
b(m+1)=b(m)−J−1(m)∇ℓ(b(m))
gdzie
ℓ jest prawdopodobieństwem logarytmicznym i
Jm będzie albo zaobserwowany lub oczekiwany Hesjan o prawdopodobieństwie dziennika.
Nasza funkcja połączenia jest funkcją która odwzorowuje średnią warunkową μ i = E ( y i | x i ) na nasz predyktor liniowy, więc naszym modelem średniej jest g ( μ i ) = x T i β . Niech h będzie funkcją odwrotnego łącza odwzorowującą predyktor liniowy na średnią.gμi=E(yi|xi)g(μi)=xTiβh
Dla regresji logistycznej istnieje prawdopodobieństwo Bernoulliego z niezależnymi obserwacjami, więc
Biorąc pochodne,
ℓ(b;y)=∑i=1nyilogh(xTib)+(1−yi)log(1−h(xTib)).
= n ∑ i=1xijh′(x T i b)( y i∂ℓ∂bj=∑i=1nyih(xTib)h′(xTib)xij−1−yi1−h(xTib)h′(xTib)xij
=∑ixijh′(x T i b)=∑i=1nxijh′(xTib)(yih(xTib)−1−yi1−h(xTib))
=∑ixijh′(xTib)h(xTib)(1−h(xTib))(yi−h(xTib)).
Za pomocą łącza kanonicznego
Załóżmy teraz, że używamy kanonicznej funkcji łącza . Następnie g - 1 c ( x ) : = h c ( x ) =gc=logit soh ′ c =hc⋅(1-hc),co oznacza, że upraszcza to
∂g−1c(x):=hc(x)=11+e−xh′c=hc⋅(1−hc)
, tak
∇£ -l(b,y)=X, T(y - y ).
Ponadto nadal używahc,
∂2
∂ℓ∂bj=∑ixij(yi−hc(xTib))
∇ℓ(b;y)=XT(y−y^).
hc∂2ℓ∂bk∂bj=−∑ixij∂∂bkhc(xTib)=−∑ixijxik[hc(xTib)(1−hc(xTib))].
Niech
W=diag(hc(xT1b)(1−hc(xT1b)),…,hc(xTnb)(1−hc(xTnb)))=diag(y^1(1−y^1),…,y^n(1−y^n)).
Then we have
H=−XTWX
and note how this doesn't have any
yi in it anymore, so
E(H)=H (we're viewing this as a function of
b so the only random thing is
y itself). Thus we've shown that Fisher scoring is equivalent to Newton-Raphson when we use the canonical link in logistic regression. Also by virtue of
y^i∈(0,1) −XTWXy^i010H
z=W−1(y−y^) and note that
∇ℓ=XT(y−y^)=XTWz.
All together this means that we can optimize the log likelihood by iterating
b(m+1)=b(m)+(XTW(m)X)−1XTW(m)z(m)
and
(XTW(m)X)−1XTW(m)z(m) is exactly
β^ for a weighted least squares regression of
z(m) on
X.
Checking this in R
:
set.seed(123)
p <- 5
n <- 500
x <- matrix(rnorm(n * p), n, p)
betas <- runif(p, -2, 2)
hc <- function(x) 1 /(1 + exp(-x)) # inverse canonical link
p.true <- hc(x %*% betas)
y <- rbinom(n, 1, p.true)
# fitting with our procedure
my_IRLS_canonical <- function(x, y, b.init, hc, tol=1e-8) {
change <- Inf
b.old <- b.init
while(change > tol) {
eta <- x %*% b.old # linear predictor
y.hat <- hc(eta)
h.prime_eta <- y.hat * (1 - y.hat)
z <- (y - y.hat) / h.prime_eta
b.new <- b.old + lm(z ~ x - 1, weights = h.prime_eta)$coef # WLS regression
change <- sqrt(sum((b.new - b.old)^2))
b.old <- b.new
}
b.new
}
my_IRLS_canonical(x, y, rep(1,p), hc)
# x1 x2 x3 x4 x5
# -1.1149687 2.1897992 1.0271298 0.8702975 -1.2074851
glm(y ~ x - 1, family=binomial())$coef
# x1 x2 x3 x4 x5
# -1.1149687 2.1897992 1.0271298 0.8702975 -1.2074851
and they agree.
Non-canonical link functions
Now if we're not using the canonical link we don't get the simplification of h′h(1−h)=1 in ∇ℓ so H becomes much more complicated, and we therefore see a noticeable difference by using E(H) in our Fisher scoring.
Here's how this will go: we already worked out the general ∇ℓ so the Hessian will be the main difficulty. We need
∂2ℓ∂bk∂bj=∑ixij∂∂bkh′(xTib)(yih(xTib)−1−yi1−h(xTib))
=∑ixijxik[h′′(xTib)(yih(xTib)−1−yi1−h(xTib))−h′(xTib)2(yih(xTib)2+1−yi(1−h(xTib))2)]
Via the linearity of expectation all we need to do to get E(H) is replace each occurrence of yi with its mean under our model which is μi=h(xTiβ). Each term in the summand will therefore contain a factor of the form
h′′(xTib)(h(xTiβ)h(xTib)−1−h(xTiβ)1−h(xTib))−h′(xTib)2(h(xTiβ)h(xTib)2+1−h(xTiβ)(1−h(xTib))2).
But to actually do our optimization we'll need to estimate each
β, and at step
m b(m) is the best guess we have. This means that this will reduce to
h′′(xTib)(h(xTib)h(xTib)−1−h(xTib)1−h(xTib))−h′(xTib)2(h(xTib)h(xTib)2+1−h(xTib)(1−h(xTib))2)
=−h′(xTib)2(1h(xTib)+11−h(xTib))
=−h′(xTib)2h(xTib)(1−h(xTib)).
This means we will use
J with
Jjk=−∑ixijxikh′(xTib)2h(xTib)(1−h(xTib)).
Now let
W∗=diag(h′(xT1b)2h(xT1b)(1−h(xT1b)),…,h′(xTnb)2h(xTnb)(1−h(xTnb)))
and note how under the canonical link
h′c=hc⋅(1−hc) reduces
W∗ to
W from the previous section. This lets us write
J=−XTW∗X
except this is now
E^(H) rather than necessarily being
H itself, so this can differ from Newton-Raphson. For all
i W∗ii>0 so aside from numerical issues
J will be negative definite.
We have
∂ℓ∂bj=∑ixijh′(xTib)h(xTib)(1−h(xTib))(yi−h(xTib))
so letting our new working response be
z∗=D−1(y−y^) with
D=diag(h′(xT1b),…,h′(xTnb)), we have
∇ℓ=XTW∗z∗.
All together we are iterating
b(m+1)=b(m)+(XTW∗(m)X)−1XTW∗(m)z∗(m)
so this is still a sequence of WLS regressions except now it's not necessarily Newton-Raphson.
I've written it out this way to emphasize the connection to Newton-Raphson, but frequently people will factor the updates so that each new point b(m+1) is itself the WLS solution, rather than a WLS solution added to the current point b(m). If we wanted to do this, we can do the following:
b(m+1)=b(m)+(XTW∗(m)X)−1XTW∗(m)z∗(m)
=(XTW∗(m)X)−1(XTW∗(m)Xb(m)+XTW∗(m)z∗(m))
=(XTW∗(m)X)−1XTW∗(m)(Xb(m)+z∗(m))
so if we're going this way you'll see the working response take the form
η(m)+D−1(m)(y−y^(m)), but it's the same thing.
Let's confirm that this works by using it to perform a probit regression on the same simulated data as before (and this is not the canonical link, so we need this more general form of IRLS).
my_IRLS_general <- function(x, y, b.init, h, h.prime, tol=1e-8) {
change <- Inf
b.old <- b.init
while(change > tol) {
eta <- x %*% b.old # linear predictor
y.hat <- h(eta)
h.prime_eta <- h.prime(eta)
w_star <- h.prime_eta^2 / (y.hat * (1 - y.hat))
z_star <- (y - y.hat) / h.prime_eta
b.new <- b.old + lm(z_star ~ x - 1, weights = w_star)$coef # WLS
change <- sqrt(sum((b.new - b.old)^2))
b.old <- b.new
}
b.new
}
# probit inverse link and derivative
h_probit <- function(x) pnorm(x, 0, 1)
h.prime_probit <- function(x) dnorm(x, 0, 1)
my_IRLS_general(x, y, rep(0,p), h_probit, h.prime_probit)
# x1 x2 x3 x4 x5
# -0.6456508 1.2520266 0.5820856 0.4982678 -0.6768585
glm(y~x-1, family=binomial(link="probit"))$coef
# x1 x2 x3 x4 x5
# -0.6456490 1.2520241 0.5820835 0.4982663 -0.6768581
and again the two agree.
Comments on convergence
Finally, a few quick comments on convergence (I'll keep this brief as this is getting really long and I'm no expert at optimization). Even though theoretically each J(m) is negative definite, bad initial conditions can still prevent this algorithm from converging. In the probit example above, changing the initial conditions to b.init=rep(1,p)
results in this, and that doesn't even look like a suspicious initial condition. If you step through the IRLS procedure with that initialization and these simulated data, by the second time through the loop there are some y^i that round to exactly 1 and so the weights become undefined. If we're using the canonical link in the algorithm I gave we won't ever be dividing by y^i(1−y^i) to get undefined weights, but if we've got a situation where some y^i are approaching 0 or 1, such as in the case of perfect separation, then we'll still get non-convergence as the gradient dies without us reaching anything.