Suavizados: P-splines y no linealidad

1. La vida no es lineal

En Cementos del Cauca, el equipo de control de calidad estudia la relación entre la temperatura de curado y la dureza del cemento. Un analista junior propone una regresión lineal: dureza ~ temperatura. La línea recta muestra una pendiente negativa marginal, pero el scatter plot revela algo que la recta oculta: la dureza es máxima alrededor de 25°C y decae en ambos extremos. La relación es cuadrática invertida, no lineal.

Peor aún, la variabilidad de la dureza aumenta en los extremos temperatura: a 5°C y a 60°C, las mediciones son mucho más dispersas que a 25°C. Un modelo lineal con varianza constante no puede capturar ni la media no lineal ni la heterocedasticidad.

2. P-splines: la idea fundamental

Los P-splines (penalized B-splines) combinan dos ideas:

1. B-splines: funciones base polinómicas a trozos definidas por un conjunto de nodos (knots) equiespaciados. Cada B-spline es no nula solo en un intervalo local, lo que los hace computacionalmente eficientes.

2. Penalización de diferencias: se imponen una penalización sobre las diferencias entre coeficientes adyacentes para controlar la rugosidad.

El predictor suavizado para un parámetro $\theta_k$ se escribe:

$\eta_k(x) = \sum_{j=1}^{M} \beta_{kj} B_j(x)$

donde $B_j(x)$ son las funciones base B-spline y la penalización es:

$\lambda_k \sum_{j=d+1}^{M} (\Delta^d \beta_{kj})^2$

con $\Delta^d$ el operador de diferencias de orden $d$ (típicamente $d=2$) y $\lambda_k$ el parámetro de suavizado que controla el balance entre fidelidad a los datos y suavidad.

1. 1

   Definir los nodos

   Por defecto, pb() coloca 20 nodos equiespaciados en el rango de x. Para datos con pocos puntos, se reduce automáticamente.

2. 2

   Construir las bases B-spline

   Se generan M funciones base polinómicas de grado 3 (cúbicas) centradas en los nodos. Cada base es local: no nula solo en un vecindario.

3. 3

   Penalizar las diferencias

   La penalización de 2ª diferencia obliga a los coeficientes a cambiar suavemente. λ grande → línea recta; λ pequeño → curva muy flexible.

4. 4

   Estimar λ por GAIC o CV

   gamlss() busca el λ óptimo minimizando el GAIC o por validación cruzada generalizada (GCV). El usuario puede fijar λ con el argumento spar.

3. Sintaxis en gamlss: pb() y cy()

La función pb() (penalized B-spline) es la interfaz principal:

gamlss(y ~ pb(x), family = NO, data = df)

Esto suaviza la relación entre x e y en el predictor de $\mu$. Para suavizar otros parámetros:

gamlss(y ~ pb(x),
       sigma.formula = ~ pb(x),   # suavizado para sigma
       nu.formula = ~ pb(x),      # suavizado para nu
       family = BCT, data = df)

4. Caso de estudio: Cementos del Cauca

Ajustemos tres modelos crecientes en complejidad:

El modelo 3 (P-splines para $\mu$ y $\sigma$) es claramente superior: $\Delta AIC \approx 240$ respecto al lineal. El gráfico de term.plot() muestra:

• Para $\mu$: una curva invertida con máximo en ~25°C, confirmando la intuición física de que existe una temperatura óptima de curado.
• Para $\sigma$: una forma de V con mínimo en ~25°C: la dispersión es mínima en la temperatura óptima y máxima en los extremos.

5. Control del suavizado

El parámetro $\lambda$ controla cuánto se suaviza la curva. En gamlss(), se estima automáticamente minimizando el GAIC. Pero el usuario puede intervenir:

6. Splines cíclicos para datos estacionales

Para datos con periodicidad (ej. demanda mensual con estacionalidad anual), cy() crea un spline cíclico que fuerza continuidad en los extremos:

gamlss(demanda ~ cy(mes, cyclic = TRUE), family = GA, data = df)

Esto es esencial para modelar la estacionalidad en la demanda de inventarios (Módulo 13) o la frecuencia de siniestros por mes (Módulo 17).

7. Cierre

¿Qué modelarías tú?

Toma una relación no lineal de tu dominio (temperatura vs rendimiento, experiencia vs salario, dosis vs respuesta). Ajusta un modelo con pb() para $\mu$ y otro con pb() también para $\sigma$. ¿La dispersión cambia con la covariable? Si sí, ¿qué implicancias tiene para tus decisiones?

Referencias: Eilers & Marx (1996), Stat. Science 11(2); Stasinopoulos et al. (2017), Cap. 5; Rigby et al. (2019), Cap. 3.