1. Introducción

El problema de regresión consiste en describir la relación entre una variable aleatoria de respuesta $Y\in\R$ y un vector de variables aleatorias explicativas $\vX\in\R^k,k\in\Z^+$ bajo la consideración de un número $n\in\Z^+$ de observaciones dadas.

La primera solución sistemática suele atribuirse a Francis Galton en sus estudios sobre herencia y regresión hacia la media [Gal94]. Formalmente, el llamado modelo lineal supone la existencia de un predictor lineal $\beta\in\R^{k+1}$ tal que $$ Y|\{\vX = \vx\} \approx \beta_0 + \beta_1 \vx_1 + \cdots + \beta_k \vx_k, $$ donde la aproximación se da de acuerdo a un error normal de media $0$.

El desarrollo posterior de la teoría, particularmente mediante el uso de transformaciones sobre las variables, permitió ampliar su alcance a phenomena cuya relación subyacente no parecía ser estrictamente lineal. Sin embargo, estas modificaciones continuaban imponiendo restricciones importantes sobre la forma funcional de la relación entre las variables.

Hacia finales del siglo XX, y a la par de la creciente capacidad computacional, surgió el interés por plantear modelos más flexibles. Una primer propuesta a destacar es el modelo lineal generalizado introducido en 1972 [NW72], el cual ahora supone que el predictor lineal solo determina la media de la variable respuesta $\mu := \E[Y]$ y permite que esta siga cualquier distribución en la familia exponencial de distribuciones $\EF$.

Una segunda propuesta es el modelo aditivo introducido una década después en 1981 [FS81], que es más cercano a métodos computacionales y propone cambiar el paradigma del predictor lineal a uno de un predictor aditivo. De nuestro interés será presentar el modelo aditivo generalizado introducido en 1986 [HT86], que surge de incorporar las dos perspectivas anteriores y se sirve de la capacidad computacional actual (Figura 1).

2. Definición

Sean una variable aleatoria de respuesta $Y\in\R$ y un vector de variables aleatorias explicativas $\vX := \m{X_1 & \ldots & X_k}^t\in\R^k$, ambos en un espacio de probabilidad común $(\Omega, {\cal F}, \P)$.

En la literatura [Woo17][MPV06][RR26] se enuncia la definición de un modelo de regresión, que es el caso de un GAM, desde dos componentes. La suposición que $Y\mid\vX=\vx \sim \EF(\theta_\vx,\phi)$ se conoce como la componente aleatoria del modelo, mientras que la suposición referente al predictor aditivo se conoce como la componente sistemática del modelo. Es usual escribir $\mu_\vx := \E\big[Y\mid\vX=\vx\big]$ y $\eta_\vx := \alpha+f_1(\vx_1)+\cdots+f_k(\vx_k)$, de modo que la relación entre las componentes se escribe $g(\mu_\vx) = \eta_\vx$ o, simplemente, $g(\mu) = \eta$.

La Figura 2 ubica la propuesta de un GAM entre las alternativas al problema de explicar una relación no lineal. En comparación con el enfoque de transformaciones, que impone una restricción paramétrica global sobre todo el soporte de la covariable, los GAMs ofrecen una flexibilidad superior al exigir únicamente condiciones de regularidad local (funciones suaves a trozos). Es decir, un GAM realiza ajustes locales adaptativos directamente de los datos sin distorsionar la estimación en otras regiones del dominio.

Por otro lado, conforme la dimensión $k = \dim\vX$ aumenta, el volumen del espacio muestral en $\R^k$ crece exponencialmente y provoca una dispersión extrema de los datos (sparsity) que dificulta una estimación de un predictor no necesariamente aditivo. Este fenómeno se conoce como la maldición de la dimensión y decimos que limitarnos a un predictor aditivo lo resuelve porque reduce el problema a estimar $k$ funciones unidimensionales sobre $\R$.

3. Estimación del modelo

La función práctica de cualquier propuesta de modelación estadística consiste en tomar un conjunto de $n\in\Z^+$ observaciones dadas $(y_1,\vx_{1 \bullet}),\ldots,(y_n,\vx_{n \bullet})$ y utilizarlas para estimar con la mayor precisión posible las componentes aleatoria y sistemática. Es aquí donde la teoría abstracta debe materializarse en algoritmos concretos.

Hasta este punto, en el caso de un GAM nos topamos con un enorme abismo metodológico para cumplir dicha tarea. A diferencia de lo que ocurre en un modelo lineal, el problema es en principio no paramétrico porque debemos evaluar e identificar funciones candidatas a ser el predictor aditivo y estas habitan en espacios de dimensión infinita.

Si permitiéramos que el algoritmo elija cualquier trayectoria libremente, el modelo arrojaría resultados objetivamente malos. Para cruzar este abismo y concretar la teoría, en lo siguiente construiremos una estrategia detallada de estimación dividida en tres etapas.

3.1 RKHS

Para bajar el problema de una estimación arbitraria no paramétrica a un terreno que se pueda trabajar de manera computacional, nos vemos obligados a sacrificar un poco de generalidad matemática en favor de la viabilidad práctica. La estrategia consiste en renunciar a la búsqueda ciega sobre el espacio de todas las funciones suaves a trozos y, en su lugar, limitar nuestro radio de acción a un espacio mucho más pequeño y controlado.

En particular, buscamos subespacios vectoriales de $\R^\R$ que tengan un sentido de geometría como es el que da tener un producto interno. Además, buscamos que la evaluación puntual sea compatible con dicha geometría para facilitar el proceso de estimación.

En un RKHS ${\cal H} \subseteq \R^\R$, como consecuencia del teorema de representación de Riesz, para cada $x\in\R$ existe una única función $K_x\in{\cal H}$ que satisface la identidad $$ f_j(x) = \delta_x(f_j) = \inp{f_j,K_x}, \qquad f_j\in{\cal H}. $$ El tener una familia $\{K_x\}_{x\in\R}$ de funciones representantes es precisamente la razón por la que nos interesa traducir la teoría de los GAMs a los RKHS.

Dado que $K_y(x) = \delta_x(K_y) = \inp{K_y, K_x}$, vemos que el núcleo de un RKHS es una función simétrica y definida positiva. En general, es complicado hacer el cálculo de $K$ para un ${\cal H}$ abstracto y en su lugar procedemos en orden inverso mediante el siguiente teorema.

Más aún, en muchos problemas de optimización sobre un RKHS ocurre que las soluciones dependen únicamente de las evaluaciones de las funciones en los datos observados aunque el espacio sea infinito-dimensional. Esto queda formalizado en el llamado teorema del representador y más adelante será nuestra principal herramienta para proceder con la estimación.

3.2 Verosimilitud para funciones

Una vez redujimos el espacio de candidatos a ser nuestro predictor aditivo, aún es necesario desarrollar una forma de evaluar e identificar al mejor de ellos. Para esto, haremos uso de la teoría de verosimilitud asintótica adaptada a funciones [Woo17].

Recordemos que, bajo la suposición $Y\mid\vX=\vx \sim \EF(\theta_\vx,\phi)$ con funciones de la distribución $a,b,c$ dadas, la media está determinada por $\mu := \E[Y\mid\vX=\vx] = b'(\theta)$. Esto induce el mapa $$ \theta_\vx \longrightarrow_{b'} \mu_\vx = \E[Y|X = \vx] \longrightarrow_{g} \eta_\vx = \alpha+f_1(x_1)+\cdots+f_k(x_k), $$ para cada $\vx = \m{x_1 & \ldots & x_k}^t\in\R^k$. Asumiendo inyectividad, tenemos el mapeo inverso $$ f \longrightarrow \m{\eta_{\vx_{1 \bullet}} \\ \vdots \\ \eta_{\vx_{n \bullet}}} \longrightarrow_{g^{-1}} \m{\mu_{\vx_{1 \bullet}} \\ \vdots \\ \mu_{\vx_{1 \bullet}}} \longrightarrow_{(b')^{-1}} \m{\theta_{\vx_{1 \bullet}} \\ \vdots \\ \theta_{\vx_{1 \bullet}}} =: \theta_f. $$ De modo que a cada predictor aditivo $f$ le podemos asociar un parámetro $\theta_f$, no necesariamente de forma inyectiva. Usaremos esto para definir una verosimilitud en nuestros candidatos.

Observemos que si probamos dos funciones $f,f^*$ que evalúan lo mismo en las observaciones $\vx_{1\bullet},\ldots,\vx_{n\bullet}$, entonces arrojarán la misma log-verosimilitud y no podremos diferenciar una de la otra o distinguir cuál es una mejor estimación. Apelando al principio de que lo más simple es mejor en caso de empate, enunciaremos las siguientes dos definiciones.

Aunque la elección del vector $\lambda\in\R^k$ ha sido tradicionalmente considerada como parte del diseño del modelo, en la actualidad es común dejar que la computadora las calcule algorítmicamente a través de métodos de validación cruzada o parecidos [Woo17]. Una observación es que $\lambda_j,j\in[k]$ determina el peso que le damos a la penalización en $\vx_{\bullet j}$, ya que si $\lambda \to\infty$ el mínimo de $l_P$ solo percibe líneas y si $\lambda\to0^+$ la penalización es despreciable.

Como mencionamos antes, existe un teorema de análisis funcional conocido como teorema del representador que podemos aplicar en el problema de maximizar el efecto de cada $f_j$ en $l_P$. Este nos asegura que dicho máximo vive en el subespacio de ${\cal H}_j$ definido por $${\cal H}_j|\vx_{1j},\ldots,\vx_{nj} := \Span\Big\{K(\vx_{1j},\bullet),\ldots,K(\vx_{nj},\bullet)\Big\}. $$ Esto es, la solución es una combinación lineal de las funciones determinadas por el núcleo generador de ${\cal H}_j$ evaluado en cada una de las observaciones. De este modo, nuestro problema se reduce a estimar parámetros $\beta_{1j},\ldots,\beta_{kj}\in\R,j\in[k]$ para los que el máximo predictor aditivo es de la forma $$ \alpha+\sum_{j=1}^k\sum_{i=1}^n \beta_{ij}K(\vx_{ij},x_j), \qquad \m{x_1 & \ldots & x_k}^t \in\R^k. $$ Este problema es paramétrico y se resuelve del mismo modo en que se resuelve un GLM, mediante métodos numéricos como es el llamado Penalized Iterative Least Squares.

3.3 Elección de los parámetros de suavidad

Ya hemos reducido el problema original de estimar funciones a un problema paramétrico gobernado por la log-verosimilitud penalizada. Sin embargo, aún queda pendiente determinar el valor adecuado de los parámetros de suavidad $\lambda_1,\ldots,\lambda_k$.

En particular, buscamos obtener un modelo con buena capacidad de predicción sobre nuevos datos y para ello se utilizan métodos de validación estadística. La idea general es separar una parte de los datos, ajustar el modelo sobre el conjunto restante y posteriormente comparar la predicción obtenida con los valores reales observados.

Si denotamos por $\widehat f_\lambda^{(-i)}$ al predictor ajustado excluyendo la observación $(y_i,\vx_{i\bullet})$, entonces el error de predicción asociado a dicha observación puede medirse mediante la discrepancia $$ L\Big( y_i, \widehat f_\lambda^{(-i)}(\vx_{i\bullet}) \Big), $$ donde $L$ es una función de pérdida apropiada, como es la devianza. Este mismo principio se puede aplicar a la exclusión de subconjuntos de observaciones.

Promediando estos errores sobre múltiples particiones de los datos se obtiene una estimación del error de generalización del modelo para cada $\lambda\in\R^k$ dado. Los parámetros de suavidad se eligen entonces minimizando dicho error.

En la actualidad, los GAMs suelen determinar automáticamente a $\lambda\in\R^k$ mediante algoritmos más modernos como validación cruzada generalizada (GCV) o máxima verosimilitud restringida (REML).

3.4 Validación del modelo

Una vez ajustado el modelo y determinados los parámetros de suavidad, es necesario evaluar qué tan bien el GAM describe los datos observados y qué tan adecuada es su capacidad de predecir nuevos datos.

Due to los GAMs se construyen sobre la familia exponencial de distribuciones, la principal medida de ajuste utilizada es la devianza [MPV06], como mencionamos antes. Recordemos que esta cuantifica la discrepancia entre el modelo ajustado y un modelo saturado que reproduce exactamente las observaciones.

En particular, si $D_{\mathrm{res}}$ denota la devianza residual del modelo ajustado y $D_{\mathrm{null}}$ la devianza del modelo nulo, es común considerar la proporción de devianza explicada $$ D_{\mathrm{exp}} := 1-\frac{D_{\mathrm{res}}}{D_{\mathrm{null}}}, $$ como un análogo no lineal del coeficiente de determinación $R^2$ de los modelos lineales.

4. Ejemplo: GAM para la abundancia de camarón rosado

En el artículo "Generalized Additive Models Used to Predict Species Abundance in the Gulf of Mexico: An Ecosystem Modeling Tool" [DA13] se estudia la abundancia del camarón rosado (Farfantepenaeus duorarum) mediante el ajuste de un GAM.

Para el ajuste del modelo se consideran las variables ambientales profundidad (depth), clorofila a (chlorophyll a), temperatura, oxígeno disuelto, tipo de sedimento y una variable de corrección por esfuerzo de muestreo (offset(effort)).

Una de las principales razones para utilizar un GAM es que las especies responden de manera no lineal a las condiciones ambientales. Por ejemplo, no es razonable asumir que un incremento continuo en la temperatura produzca necesariamente un aumento proporcional en la abundancia de una especie. En muchos casos existen rangos óptimos de temperatura, profundidad u oxígeno donde la especie alcanza mayor abundancia.

Bajo este enfoque, la variable respuesta debe suponerse perteneciente a la familia exponencial de distribuciones. En este caso, la abundancia fue modelada mediante una distribución Binomial Negativa debido a que los datos presentan una gran cantidad de observaciones iguales a cero, como puede apreciarse en el mapa de captura.

En comparación con una distribución Poisson, cuya restricción teórica establece que la media y la varianza deben coincidir, $$ \mathrm{Var}(Y)=\mathbb{E}[Y], $$ la distribución Binomial Negativa permite modelar situaciones donde la variabilidad es considerablemente mayor que la media, $$ \mathrm{Var}(Y)>\mathbb{E}[Y], $$ lo cual resulta adecuado para datos ecológicos y pesqueros. Esto ocurre porque el camarón rosado no se distribuye uniformemente a lo largo del océano; existen regiones con muy baja abundancia y otras donde se concentran grandes cantidades de individuos.

Considerando lo anterior, el modelo ajustado tiene la forma $$ g(\mu)= s(\text{depth})+ s(\text{chla})+ s(\text{temp})+ s(\text{oxygen})+ \text{factor}(\text{sedimenttype})+ \text{offset}(\log(\text{effort})). $$

En este modelo se ajustan funciones suaves para las variables profundidad, clorofila a, temperatura y oxígeno disuelto. La variable tipo de sedimento es categórica y el término offset corrige por diferencias en el esfuerzo de muestreo.

Por ejemplo, supóngase que en dos estaciones de muestreo se realizan arrastres con distinta duración:

• Estación A: arrastre de 1 hora.
• Estación B: arrastre de 3 horas.

Naturalmente, la estación B capturará una mayor cantidad de camarones simplemente porque el tiempo de pesca fue mayor. El término offset evita interpretar esto como una mayor abundancia real de la especie.

Las distintas mediciones ambientales utilizadas en el estudio fueron obtenidas por diferentes instituciones, como se muestra en la siguiente tabla.

Due to dichas instituciones aplican distintos procesos de validación y filtrado de datos, no es posible reproducir exactamente los resultados originales del artículo.

De acuerdo con los autores, el ajuste del modelo puede realizarse mediante la librería mgcv de R. Un posible código para ajustar el modelo es el siguiente:

library(mgcv)

modelo <- gam(
  abundance ~ s(depth) +
    s(chla) +
    s(temp) +
    s(oxygen) +
    factor(sedimenttype) +
    offset(log(effort)),
  family = negative.binomial(theta = theta_est),
  method = "GCV.Cp",
  data = datos
)

La función s() indica que el modelo debe ajustar funciones suaves para las variables continuas. Además, el argumento

method = "GCV.Cp"

utiliza el criterio Generalized Cross Validation (GCV), el cual selecciona automáticamente el nivel óptimo de suavizamiento para evitar tanto el sobreajuste como modelos excesivamente rígidos.

Con este ajuste se obtuvo información acerca de la contribución de cada variable regresora sobre la abundancia estimada del camarón rosado.

Como puede observarse en la gráfica correspondiente a la profundidad, el efecto positivo de esta variable sobre la abundancia del camarón rosado se presenta aproximadamente entre los 25 y 40 metros, lo cual coincide con los valores reportados en la literatura para el hábitat natural de la especie, ubicado principalmente entre 30 y 40 metros de profundidad.

Por otra parte, la temperatura óptima de hábitat para el camarón rosado se encuentra aproximadamente entre 18 y 32 $^\circ$C. En concordancia con ello, la gráfica de temperatura muestra un efecto positivo de esta variable en un intervalo cercano a los 15 y 30 $^\circ$C.

De manera similar, los resultados asociados a la clorofila a indican que concentraciones menores a 15 mg/m$^3$ favorecen la abundancia del camarón rosado.

Finalmente, para el oxígeno disuelto se sabe que valores entre 2.0 y 4.5 mL/L son adecuados para el desarrollo de la especie. Los resultados obtenidos en el ajuste del modelo muestran efectos positivos de esta variable principalmente entre 2.5 y 4.2 mL/L, en concordancia con los rangos reportados para el hábitat del camarón rosado.

Referencias

1. [Gal94] Francis Galton. Natural Inheritance. Macmillan and Company, New York, 5th edition, 1894.
2. [NW72] J. A. Nelder and R. W. M. Wedderburn. Generalized linear models. Journal of the Royal Statistical Society: Series A (General), 135(3):370–384, 1972.
3. [FS81] Jerome H. Friedman and Werner Stuetzle. Projection pursuit regression. Journal of the American Statistical Association, 76(376):817–823, 1981.
4. [HT86] T. Hastie and R. Tibshirani. Generalized additive models. Statistical Science, 1(3):297–318, 1986.
5. [Woo17] Woo17, S. N. (2017). Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC.
6. [MPV06] Simon N. Woo17. Generalized Additive Models. Texts in Statistical Science. Chapman and Hall/CRC, 2nd edition, 2017.
7. [RR26] L. Ramírez Ramírez. Modelos estadísticos i. http://personal.cimat.mx:8181/~RR26cia.ramirez/Modelos_Estadisticos_I/, 2026. Accesado: 27-06-2026.
8. [Aro50] N. Aronszajn. Theory of reproducing kernels. Transactions of the American Mathematical Society, 68(3):337–505, 1950.
9. [Cla19] M. Cla19. Generalized additive models. https://m-Cla19.github.io/generalized-additive-models/, 2019. Accesado: 06-05-2026.
10. [VV23] Vijaya Veparala and Kalpana Vattikunta. Big data y diferentes enfoques de clustering subespacial: De la promoción en redes sociales al mapeo genómico. Salud, Ciencia y Tecnología, 3:413, 06 2023.
11. [DA13]Michael Drexler and Cameron H. Ainsworth. Generalized additive models used to predict species abundance in the gulf of mexico: An ecosystem modeling tool. PLOS ONE, 8(5):1–7, 05 2013.