Tanto en el modelo de regresión lineal como en el de regresión logística, asumimos una cierta distribución condicional para \(Y|\XX\) que depende de un conjunto de parámetros \(\bftheta\).
En regresión lineal, suponemos \(Y|\XX=\xx\sim\mathcal{N}(\bftheta^\top\xx,\sigma^2)\).
En regresión logística, que recordemos está pensado para problemas de clasificación, suponemos \(Y|\XX=\xx\sim\text{Binom}(1,g(\bftheta^\top\xx))\), donde \(g(z)=1/(1+e^{-z})\).
En esta sección mostraremos que ambos métodos son casos especiales de una familia más amplia de modelos, llamados Modelos Lineales Generalizados (GLM) . También mostraremos cómo otros modelos en GLM pueden derivarse y aplicarse a problemas de regresión y clasificación.
Para una variable \(\mathbf{Z}\in\RR^q\), denotamos \(\media\left[\mathbf{Z}\right]\) su esperanza, dada por \[
\media\left[\bfZ\right]:=\int \mathbf{z}\,d\prob_{\bfZ}(\zz)\in\RR^q,
\]
mientras que \(\cov\left[\bfZ\right]\) es su covarianza definida por \[
\cov\left[\bfZ\right]:=\media\left[(\bfZ-\media[\bfZ])(\bfZ-\media[\bfZ])^\top\right]\in\RR^{q\times q}.
\]
La familia exponencial
Definición 1. (Familia exponencial) Decimos que una clase de distribuciones pertenece a la familia exponencial si se puede escribir en la forma
con \(h:\RR\to\RR\), \(\bfeta\in\RR^q\), \(T:\RR\to\RR^q\) y \(A:\RR^q\to\RR\).
En la expresión (1) se distinguen los siguientes elementos:
\(\bfeta\in\RR^q\) es el parámetro natural.
\(T:\RR\to\RR^q\) es el estadístico suficiente.
\(h:\RR\to\RR\) es la función base que depende únicamente de \(y\).
\(A:\RR^q\to\RR\) es la función de partición que normaliza la expresión, garantizando que la distribución integre a 1.
Observar que (1) puede expresarse también como \[
p(y|\bfeta)=\frac{h(y)\exp\left\{\bfeta^\top T(y)\right\}}{\exp\left\{A(\bfeta)\right\}}=\frac{h(y)\exp\left\{\bfeta^\top T(y)\right\}}{\displaystyle\int h(y)\exp\left\{\bfeta^\top T(y)\right\}\,dy}.
\tag{2}
\]
Para un mismo \(h\) y \(T\), la familia exponencial consiste en todas las distribuciones que se obtienen al variar el parámetro natural \(\bfeta\).
Ejemplo 1 Distribución Bernoulli
Sea \(Y\in\{0,1\}\) tal que \(Y\sim\text{Binom}(1,\phi)\), con \(\phi:=\prob[Y=1]\). Entonces
y el estadístico suficiente es \[
T(y)=(y,y^2)\in\RR^2.
\]
El resto de los elementos son: \[
b(y)=\frac{1}{\sqrt{2\pi}}, \qquad A(\bfeta)=\frac{\mu^2}{2\sigma^2}+\log\sigma=-\frac{\eta_1^2}{4\eta_2}+\frac{1}{2}\log(-2\eta_2).\]
📝
Exprese la distribución normal univariada con varianza conocida como familia exponencial.
Existen muchas otras distribuciones que son miembros de la familia exponencial: la Multinomial (que veremos más adelante), la Poisson (para modelar datos de conteo), la Gamma y la Exponencial (para modelar variables continuas no negativas, como intervalos de tiempo); la Beta y la Dirichlet (para distribuciones sobre probabilidades); y muchas más.
Importante
La función de partición \(A:\RR^q\to\RR\) queda definida a partir de la siguiente expresión:
Además, también se puede obtener la Hessiana: \[
\nabla^2 A(\bfeta)=\cov\left[T(Y)\right]\in\RR^{q\times q}.
\]
En consecuencia, \(\nabla^2 A(\bfeta)\succeq 0\), lo cual implica que \(A\) es una función convexa.
📝
¿Cuándo \(A\) es estrictamente convexa?
Construcción de GLMs
Supongamos que deseas construir un modelo para estimar el número \(y\) de clientes que llegan a tu tienda (o el número de visitas a páginas web en tu sitio) en una hora determinada, basándote en ciertas características como promociones en la tienda, publicidad reciente, clima, día de la semana, etc. Sabemos que la distribución Poisson usualmente proporciona un buen modelo para el número de visitantes. Sabiendo esto, ¿cómo podemos proponer un modelo para nuestro problema? Afortunadamente, la Poisson es una distribución de la familia exponencial, por lo que podemos aplicar un Modelo Lineal Generalizado (GLM). En esta sección, describiremos un método para construir modelos GLM para problemas como este.
Más generalmente, consideremos un problema de regresión o clasificación o regresión donde queremos predecir el valor de alguna variable aleatoria \(Y\in\RR\) a partir de un conjunto de predictoras \(\XX\in\RR^p\). Un GLM está definido a partir del siguiente supuesto:
Supuesto del modelo
\(Y|\XX=\xx\) se distribuye según una distribución de la familia exponencial con parámetro natural de la forma
El objetivo es predecir \[
h_{\bfTheta}(\xx):=\media\left[T(Y)\mid\xx;\bfTheta\right]
\]
📝
¿Cómo sería el GLM para el caso del ejemplo inicial, en el que \(Y\) representa el número de clientes que llegan a una tienda, suponiendo que \(Y|X\) sigue una distribución Poisson?
Los GLM son una clase de algoritmos de aprendizaje que tienen muchas propiedades deseables, como la facilidad de aprendizaje. Además, los modelos resultantes son frecuentemente muy efectivos. Generalmente \(T(y)=y\) y, en tal caso, GLM permite predecir \[
h_{\bftheta}(\xx)=\media[Y\mid\xx;\bfTheta].\]
Esto permite ver a los modelos de regresión lineal y regresión logística como casos particulares de GLM.
Por otra parte, la conexión entre el parámetro natural \(\bfeta\) (predictor lineal) y la media \(\media[T(Y)\mid\xx,\bfTheta]\) se establece mediante una función de enlace canónica. Por ejemplo, en regresión logística \[
\eta = \bftheta^\top \xx = \text{logit}(\underbrace{\media[Y\mid\xx,\bftheta]}_{\sigma(\bftheta^\top\xx)}),
\] donde \(\text{logit}:(0,1)\to\RR\) es la inversa de la función logística \(\sigma\); esto es \[
\text{logit}\,p=\log\frac{p}{1-p}.
\]
Ajuste de parámetros
Vamos a abordar el problema de estimación de los parámetros en un GLM de la manera habitual: definiendo el problema de optimización a partir del principio de máxima verosimilitud y luego, de ser posible, aplicando la condición de optimalidad.
Función de verosimilitud
Por simplicidad, vamos a considerar \(\bftheta\in\RR^p\) tal que \(\eta:=\bftheta^\top\xx\in\RR\).
Sea \(\{\xx_i,y_i\}_{i=1}^n\) un conjunto de datos i.i.d. Para todo \(i=1,\ldots,n\), se tiene que
Entonces \[
\ell(\tilde{\bfeta})=\sum_{i=1}^n\left[\eta_i\, T(y_i)-A(\eta_i)\right]+\text{constante},
\tag{3}
\]
Optimalidad
De (3), el problema de optimización asociado a la máxima verosimilitud es \[
\begin{array}{ll}
\text{minimizar } & \sum_{i=1}^n\left[A(\eta_i)-\eta_i\, T(y_i)\right]
\end{array}
\]
¡Cuidado! La variable de optimización es \(\bftheta\), que determina cada parámetro natural a partir de la expresión \(\eta_i=\bftheta^\top \xx_i\). Dado que tenemos convexidad, podemos plantear \(\nabla J(\bftheta)=0\) (como antes, llamamos \(J\) a la función objetivo). Tenemos \[
\frac{\partial J}{\partial \theta_k}=\sum_{i=1}^n \left[A'(\eta_i)-T(y_i)\right]x_{ik}
\]
En forma vectorial: \[
X^\top(A'(\tilde{\bfeta})-T(\yy))=\bfzero,
\tag{4}
\] donde \[
A'(\tilde{\bfeta})=(A'(\eta_1),\ldots,A'(\eta_n))^\top.
\]
Ahora bien, sabemos que \[
A'(\eta_i)=\media\left[T(Y)\mid\XX=\xx_i\right]\qquad\forall i=1,\ldots,n.
\]
Por lo tanto, la expresión (4) puede reescribirse como \[
\sum_{i=1}^n(T(y_i)-\media\left[T(Y)\mid \XX=\xx_i\right])\,\xx_i=\bfzero.
\]
Observar que lo del paréntesis es el residuo entre el valor observado de \(T(y_i)\) y el valor esperado según el modelo. La condición de optimalidad establece que la suma de estos residuos ponderados por los predictores debe ser cero. Por otra parte, no existe una solución explícita para el punto óptimo \(\hat{\bftheta}\) en general; por supuesto se recurre a métodos iterativos para aproximarlo.
IRLS (Iterative Reweigthed Least Squares)
La aplicación del método de Newton al problema de estimación de parámetros en un GLM puede reformularse como un problema de mínimos cuadrados ponderados (WLS), en el que cada iteración resuelve un sistema lineal con pesos actualizados según la varianza del modelo.
Importante
El problema de mínimos cuadrados ponderados (WLS) es \[
\begin{array}{ll}
\text{minimizar } & \frac{1}{2}\|W^{1/2}(\yy-X\bftheta)\|_2^2,
\end{array}
\] donde \(W:=\text{diag}\,(w_1,\ldots,w_n)\) contiene los pesos de cada observación. Si \(X\) es de rango completo, la solución cerrada es \[
\hat{\bfbeta}^{\text{WLS}}:=(X^\top W X)^{-1}X^\top W \yy.
\]
Recordemos que el método de Newton actualiza los parámetros mediante \[
\bftheta_{t+1}:=\bftheta_t-\left[\nabla^2 J(\bftheta_t)\right]^{-1}\nabla J(\bftheta_t).
\tag{5}
\]
Hasta ahora hemos visto que \[
\nabla J(\bftheta)=\sum_{i=1}^n(T(y_i)-\media[T(Y)\mid\XX=\xx_i])\xx_i,
\]
donde \(\media[T(Y)\mid\XX=\xx_i]=A'(\eta_i)\). Además, por propiedad, \(A''(\eta_i)=\cov\left[T(Y_i)\mid \XX=\xx_i\right]\). Por lo tanto, resulta \[
\nabla^2 J(\bftheta)=-\sum_{i=1}^n\cov\left[T(Y)\mid\XX=\xx_i\right]\xx_i\xx_i^\top=-X^\top W X.
\]
donde \[
W:=\text{diag}\,(w_1,\ldots,w_n),\qquad w_i:=\cov\left[T(Y)\mid \XX=\xx_i\right]
\]
Por lo tanto, la iteración del método de Newton resulta \[
\bftheta_{t+1}:=\bftheta_t+(X^\top W X)^{-1}\left(\sum_{i=1}^n(T(y_i)-\media[T(Y)\mid\XX=\xx_i])\xx_i\right)
\tag{6}
\]
Si logramos proponer \(\zz\in\RR^n\) tal que \[
\sum_{i=1}^n(T(y_i)-\media[T(Y)\mid\XX=\xx_i])\xx_i=X^\top W\zz,
\] entonces la actualización en (6) del vector de parámetros \(\bftheta\) puede interpretarse a partir de la solución de un WLS con pesos \(W\) y pseudo-respuesta\(\zz\). Dado que \(X^\top W \zz=\sum_{i=1}^nw_iz_i\xx_i\), podemos elegir \[
z_i:=\frac{T(y_i)-\media[T(Y)\mid\XX=\xx_i]}{w_i}.
\]
Por lo tanto, (6) resulta en \[
\bftheta_{t+1}:=\bftheta_t+(X^\top W X)^{-1}X^\top W\zz.
\]
Esto se conoce como algoritmo IRLS (Iterative Reweighted Least Squares), ya que cada iteración consiste en resolver un problema de mínimos cuadrados ponderados con pesos \(W\) y pseudo-respuesta \(\zz\), para recién actualizar el vector de parámetros \(\bftheta\). Repitiendo este procedimiento se aproxima el estimador de máxima verosimilitud del GLM, aprovechando la estructura del Hessiano para generar pasos más precisos que un simple descenso por gradiente.
A continuación, formulamos IRLS con una modificación: para que el WLS quede en forma estándar se suele tomar \(z_i:=\eta_i+(T(y_i)-\media[T(Y)\mid\XX=\xx_i])/w_i\); esto permite absorber \(\bftheta_t\) en la fórmula iterativa, ya que \(\tilde{\bfeta}=X\bftheta_t\) implica \((X^\top W X)^{-1}X^T W\tilde{\bfeta}=\bftheta_t\).
Algoritmo IRLS
Dado un parámetro inicial\(\bftheta_0\in\RR^{p}\):
repetir para \(t = 0,1,2,\ldots\):
1° Calcular \(\tilde{\bfeta} = X \bftheta_t\) y \(\bfmu\in\RR^n\) tal que \(\mu_i=\media\left[T(Y)|\XX=\xx_i\right]\).
2° Construir la matriz diagonal de pesos \(W\in\RR^{n\times n}\) tal que \(w_{i} = \cov\left[T(Y)\mid X=\xx_i\right]\).
4° Actualizar \(\bftheta_{t+1}:=(X^\top W X)^{-1}X^\top W\zz\).
hasta que el criterio de parada se satisfaga.
Regresión Softmax
Consideremos un problema de clasificación en el que \(Y \in \{1, 2, \ldots, k\}\). Por ejemplo, en lugar de clasificar correos electrónicos en dos clases (spam o no spam), podríamos clasificarlos en tres clases: spam, correo personal y correo relacionado con el trabajo. Para modelar esta situación, usaremos la distribución multinomial con un único ensayo (también llamada distribución categórica) como una distribución de la familia exponencial.
Ejemplo 3 Distribución categórica
Sea \(Y\in\RR\) tal que \(Y\sim\text{Multinomial}(1,\boldsymbol{\phi})\), con \(\bfphi\in\RR^k\) tal que \(\sum_{i=1}^k\phi_i=1\). Entonces \[
p(y\mid\bfphi)=\phi_1^{\bfone\{y=1\}}\phi_2^{\bfone\{y=2\}}\ldots \phi_k^{\bfone\{y=k\}},
\]
donde \(\bfone\{\cdot\}\) es la función indicatriz, que vale 1 si su argumento es cierto y 0 en caso contrario. Definimos el estadístico suficiente \(T(y)\in\RR^{k-1}\) mediante \[
(T(y))_i=\left\{
\begin{array}{cl}
1 & \mbox{si }y=i, \\
0 & \mbox{si }y\neq i.
\end{array}
\right.
\]
El resto de los elementos son: \[
h(y)=1,\qquad A(\bfeta)=-\log\phi_k=-\log\left(1-\sum_{i=1}^{k-1}\phi_i\right).
\]
La expresión (7) es la función de enlace, conocida como logit generalizado. Su inversa, que permite el mapeo \(\bfeta\mapsto\bfphi\), es la función softmax determinada por \[
\phi_i=\frac{\exp\{\eta_i\}}{1+\sum_{j=1}^{k-1}\exp\{\eta_j\}},\qquad i=1,\ldots,k-1,
\]
El GLM asociado parte de la suposición \[
\bfeta=\bfTheta^\top\xx,\qquad\bfTheta=\begin{pmatrix}\bftheta_1,\ldots,\bftheta_{k-1}\end{pmatrix}\in\RR^{p\times(k-1)}.
\]
Así, la distribución condicional de \(Y|\XX\) es \[
p(Y=i\mid\xx;\bfTheta):=\phi_i=\frac{\exp\{\bftheta_i^\top\xx\}}{1+\sum_{j=1}^{k-1}\exp\{\bftheta_j^\top\xx\}},\qquad i=1,\ldots,k-1,
\]\[
p(Y=k\mid\xx;\bfTheta):=\phi_k=\frac{1}{1+\sum_{j=1}^{k-1}\exp\{\bftheta_j^\top\xx\}}.
\]
En el problema de optimización asociado \[
\begin{array}{ll}
\text{minimizar } & \sum_{i=1}^n\left[A(\bfeta_i)-\bfeta_i^\top T(y_i)\right],
\end{array}
\]
la función objetivo resulta \[
J(\bfTheta)=\sum_{i=1}^n\left[\log\left(1+\sum_{j=1}^{k-1}\exp\{\bftheta_j^\top\xx_i\}\right)-\xx_i^\top\bfTheta\, T(y_i)\right].
\]
Obtengamos el vector gradiente respecto a un \(\bftheta_l\) (\(l=1,\ldots,k-1\)): \[
\begin{align}
\nabla_{\bftheta_l}J(\bfTheta)&=\sum_{i=1}^n\left[\frac{\exp\{\bftheta_l^\top\xx_i\}}{1+\sum_{j=1}^{k-1}\exp\{\bftheta_j^\top\xx_i\}}-T(y_i)_l\right]\xx_i \\
&=\sum_{i=1}^n\left[\phi_l(\xx_i)-T(y_i)_l\right]\xx_i
\end{align}
\]
Si bien este gradiente de la función objetivo se puede usar para diseñar el método de descenso por gradiente, es más común emplear IRLS, que aprovecha la información del Hessiano para generar pasos más precisos y acelerar la convergencia.
