¿Por qué el uso del método de Newton para la optimización de la regresión logística se denomina mínimos cuadrados reponderados iterativos?

Question

¿Por qué el uso del método de Newton para la optimización de la regresión logística se denomina mínimos cuadrados reponderados iterativos?

No me parece claro porque la pérdida logística y la pérdida por mínimos cuadrados son cosas completamente distintas.

Answer 1

Resumen: Los MLG se ajustan mediante Marcador Fisher que, como señala Dimitriy V. Masterov, es Newton-Raphson con el hessiano esperado en su lugar (es decir, utilizamos una estimación de la información de Fisher en lugar de la información observada). Si utilizamos la función de enlace canónica, resulta que el hessiano observado es igual al hessiano esperado, por lo que la puntuación NR y Fisher son iguales en ese caso. En cualquier caso, veremos que la puntuación de Fisher en realidad está ajustando un modelo lineal de mínimos cuadrados ponderados, y las estimaciones de los coeficientes convergen* en un máximo de la probabilidad de regresión logística. Aparte de reducir el ajuste de una regresión logística a un problema ya resuelto, también obtenemos el beneficio de poder utilizar diagnósticos de regresión lineal en el ajuste WLS final para aprender sobre nuestra regresión logística.

Voy a centrarme en la regresión logística, pero para una perspectiva más general sobre la máxima verosimilitud en los MLG recomiendo la sección 15.3 de este capítulo que repasa esto y deriva IRLS en un entorno más general (creo que es de John Fox's Análisis de regresión aplicado y modelos lineales generalizados ).

$^*$ ver comentarios al final

---

La función de verosimilitud y puntuación

Ajustaremos nuestro MLG iterando algo de la forma $$ b^{(m+1)} = b^{(m)} - J^{-1}_{(m)}\nabla \ell(b^{(m)}) $$ donde $\ell$ es la probabilidad logarítmica y $J_{m}$ será el hessiano observado o esperado de la logverosimilitud.

Nuestra función de enlace es una función $g$ que asigna la media condicional $\mu_i = E(y_i | x_i)$ a nuestro predictor lineal, por lo que nuestro modelo para la media es $g(\mu_i) = x_i^T\beta$ . Sea $h$ sea la función de enlace inversa que asigna el predictor lineal a la media.

Para una regresión logística tenemos una probabilidad de Bernoulli con observaciones independientes, por lo que $$ \ell(b; y) = \sum_{i=1}^n y_i\log h(x_i^T b) + (1 - y_i) \log(1 - h(x_i^Tb)). $$ Tomar derivados, $$ \frac{\partial \ell}{\partial b_j} = \sum_{i=1}^n \frac{y_i}{h(x_i^T b)} h'(x_i^T b) x_{ij} - \frac{1 - y_i}{1 - h(x_i^T b)} h'(x_i^T b) x_{ij} $$ $$ = \sum_{i=1}^n x_{ij} h'(x_i^T b) \left(\frac{y_i}{h(x_i^T b)} - \frac{1 - y_i}{1 - h(x_i^T b)} \right) $$ $$ = \sum_i x_{ij} \frac{h'(x_i^T b)}{h(x_i^T b)(1 - h(x_i^T b))}(y_i - h(x_i^T b)). $$

---

Utilizar el enlace canónico

Supongamos ahora que utilizamos la función de enlace canónico $g_c = \text{logit}$ . Entonces $g^{-1}_c(x) := h_c(x) = \frac{1}{1+e^{-x}}$ así que $h_c' = h_c \cdot (1-h_c)$ lo que significa que se simplifica a $$ \frac{\partial \ell}{\partial b_j} = \sum_i x_{ij} (y_i - h_c(x_i^T b)) $$ así que $$ \nabla \ell (b; y) = X^T (y - \hat y). $$ Además, sigue utilizando $h_c$ , $$ \frac{\partial^2 \ell}{\partial b_k \partial b_j} = - \sum_i x_{ij} \frac{\partial}{\partial b_k} h_c(x_i^T b) = - \sum_i x_{ij}x_{ik} \left[h_c(x_i^T b) (1 - h_c(x_i^T b))\right]. $$

Sea $$ W = \text{diag}\left(h_c(x_1^T b)(1 - h_c(x_1^T b)), \dots, h_c(x_n^T b)(1 - h_c(x_n^T b))\right) = \text{diag}\left(\hat y_1(1 - \hat y_1), \dots, \hat y_n (1 - \hat y_n)\right). $$ Entonces tenemos $$ H = -X^TWX $$ y observe que esto no tiene ningún $y_i$ en ella, así que $E(H) = H$ (estamos viendo esto como una función de $b$ así que lo único aleatorio es $y$ sí mismo). Así, hemos demostrado que la puntuación de Fisher es equivalente a la de Newton-Raphson cuando utilizamos el enlace canónico en la regresión logística. También en virtud de $\hat y_i \in (0,1)$ $-X^TWX$ será siempre estrictamente negativa definida, aunque numéricamente si $\hat y_i$ se acerca demasiado a $0$ o $1$ entonces puede que tengamos pesos redondos para $0$ que puede hacer $H$ semidefinida negativa y, por tanto, computacionalmente singular.

Ahora cree el respuesta de trabajo $z = W^{-1}(y - \hat y)$ y observe que $$ \nabla \ell = X^T(y - \hat y) = X^T W z. $$

En conjunto, esto significa que podemos optimizar la logverosimilitud iterando $$ b^{(m+1)} = b^{(m)} + (X^T W_{(m)} X)^{-1}X^T W_{(m)} z_{(m)} $$ y $(X^T W_{(m)} X)^{-1}X^T W_{(m)} z_{(m)}$ es exactamente $\hat \beta$ para una regresión por mínimos cuadrados ponderados de $z_{(m)}$ en $X$ .

Comprobación 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 

y están de acuerdo.

---

Funciones de enlace no canónicas

Ahora bien, si no utilizamos el enlace canónico no obtendremos la simplificación de $\frac{h'}{h(1-h)} = 1$ en $\nabla \ell$ así que $H$ resulta mucho más complicado, por lo que vemos una diferencia notable al utilizar $E(H)$ en nuestra puntuación Fisher.

Esto es lo que haremos: ya hemos elaborado el plan general $\nabla \ell$ por lo que el hessiano será la principal dificultad. Necesitamos $$ \frac{\partial^2 \ell}{\partial b_k \partial b_j} = \sum_i x_{ij} \frac{\partial}{\partial b_k}h'(x_i^T b) \left(\frac{y_i}{h(x_i^T b)} - \frac{1 - y_i}{1 - h(x_i^T b)} \right) $$ $$ = \sum_i x_{ij}x_{ik} \left[h''(x_i^T b) \left(\frac{y_i}{h(x_i^T b)} - \frac{1 - y_i}{1 - h(x_i^T b)} \right) - h'(x_i^T b)^2\left(\frac{y_i}{h(x_i^T b)^2} + \frac{1-y_i}{(1-h(x_i^T b))^2} \right)\right] $$

A través de la linealidad de la expectativa todo lo que necesitamos hacer para conseguir $E(H)$ es sustituir cada aparición de $y_i$ con su media según nuestro modelo, que es $\mu_i=h(x_i^T\beta)$ . Por tanto, cada término del sumando contendrá un factor de la forma $$ h''(x_i^T b) \left(\frac{h(x_i^T \beta)}{h(x_i^T b)} - \frac{1 - h(x_i^T \beta)}{1 - h(x_i^T b)} \right) - h'(x_i^T b)^2\left(\frac{h(x_i^T \beta)}{h(x_i^T b)^2} + \frac{1-h(x_i^T \beta)}{(1-h(x_i^T b))^2} \right). $$ Pero para hacer realmente nuestra optimización necesitaremos estimar cada $\beta$ y en el paso $m$ $b^{(m)}$ es la mejor suposición que tenemos. Esto significa que se reducirá a $$ h''(x_i^T b) \left(\frac{h(x_i^T b)}{h(x_i^T b)} - \frac{1 - h(x_i^T b)}{1 - h(x_i^T b)} \right) - h'(x_i^T b)^2\left(\frac{h(x_i^T b)}{h(x_i^T b)^2} + \frac{1-h(x_i^T b)}{(1-h(x_i^T b))^2} \right) $$ $$ = - h'(x_i^T b)^2\left(\frac{1}{h(x_i^T b)} + \frac{1}{1-h(x_i^T b)} \right) $$ $$ = -\frac{h'(x_i^T b)^2}{h(x_i^T b)(1-h(x_i^T b))}. $$ Esto significa que utilizaremos $J$ con $$ J_{jk} = -\sum_i x_{ij}x_{ik} \frac{h'(x_i^T b)^2}{h(x_i^T b)(1-h(x_i^T b))}. $$

Ahora dejemos que $$ W^* = \text{diag}\left(\frac{h'(x_1^T b)^2}{h(x_1^T b)(1-h(x_1^T b))} ,\dots, \frac{h'(x_n^T b)^2}{h(x_n^T b)(1-h(x_n^T b))}\right) $$ y observe cómo bajo el enlace canónico $h_c' = h_c \cdot (1-h_c)$ reduce $W^*$ a $W$ de la sección anterior. Esto nos permite escribir $$ J = -X^TW^*X $$ excepto que esto es ahora $\hat E(H)$ en lugar de ser necesariamente $H$ por lo que puede diferir de Newton-Raphson. Para todos los $i$ $W_{ii}^* > 0$ así que aparte de las cuestiones numéricas $J$ será definida negativa.

Tenemos $$ \frac{\partial \ell}{\partial b_j} = \sum_i x_{ij} \frac{h'(x_i^T b)}{h(x_i^T b)(1 - h(x_i^T b))}(y_i - h(x_i^T b)) $$ así que dejemos que nuestra nueva respuesta de trabajo sea $z^* = D^{-1}(y-\hat y)$ con $D=\text{diag}\left(h'(x_1^T b), \dots, h'(x_n^T b)\right)$ tenemos $\nabla \ell = X^TW^*z^*$ .

Todos juntos estamos iterando $$ b^{(m+1)} = b^{(m)} + (X^T W_{(m)}^* X)^{-1}X^T W_{(m)}^* z_{(m)}^* $$ así que esto sigue siendo una secuencia de regresiones WLS excepto que ahora no es necesariamente Newton-Raphson.

---

Lo he escrito así para enfatizar la conexión con Newton-Raphson, pero con frecuencia la gente factorizará las actualizaciones de modo que cada nuevo punto $b^{(m+1)}$ es la propia solución WLS, en lugar de una solución WLS añadida al punto actual $b^{(m)}$ . Si quisiéramos hacer esto, podemos hacer lo siguiente: $$ b^{(m+1)} = b^{(m)} + (X^T W_{(m)}^* X)^{-1}X^T W_{(m)}^* z_{(m)}^* $$ $$ = (X^T W_{(m)}^* X)^{-1}\left(X^T W_{(m)}^* Xb^{(m)}+ X^TW^*_{(m)}z_{(m)}^* \right) $$ $$ = (X^T W_{(m)}^* X)^{-1}X^TW_{(m)}^*\left(Xb^{(m)}+ z_{(m)}^* \right) $$ así que si vamos por este camino verás que la respuesta de trabajo toma la forma $\eta^{(m)} + D^{-1}_{(m)}(y - \hat y^{(m)})$ pero es lo mismo.

---

Confirmemos que esto funciona utilizándolo para realizar una regresión probit sobre los mismos datos simulados que antes (y esto no es el enlace canónico, por lo que necesitamos esta forma más general de 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 

y de nuevo los dos se ponen de acuerdo.

---

Comentarios sobre la convergencia

Por último, algunos comentarios rápidos sobre la convergencia (seré breve porque esto se está alargando mucho y no soy un experto en optimización). Aunque teóricamente cada $J_{(m)}$ es negativo definido, las malas condiciones iniciales pueden impedir que este algoritmo converja. En el ejemplo probit anterior, cambiar las condiciones iniciales a b.init=rep(1,p) resulta en esto, y ni siquiera parece una condición inicial sospechosa. Si usted paso a través del procedimiento IRLS con que la inicialización y estos datos simulados, por segunda vez a través del bucle hay algunos $\hat y_i$ que redondean exactamente $1$ y así los pesos se vuelven indefinidos. Si estamos utilizando el enlace canónico en el algoritmo que he dado nunca vamos a estar dividiendo por $\hat y_i (1 - \hat y_i)$ para obtener pesos indefinidos, pero si tenemos una situación en la que algunos $\hat y_i$ se acercan $0$ o $1$ como en el caso de separación perfecta, entonces seguiremos obteniendo no-convergencia ya que el gradiente muere sin que alcancemos nada.