Intuición detrás de las interacciones del producto tensorial en GAM (paquete MGCV en R)

Los modelos aditivos generalizados son aquellos donde

$$y = \alpha + f_1(x_1) + f_2(x_2) + e_i$$ por ejemplo. Las funciones son suaves y se deben estimar. Generalmente por splines penalizadas. MGCV es un paquete en R que lo hace, y el autor (Simon Wood) escribe un libro sobre su paquete con ejemplos de R. Ruppert y col. (2003) escriben un libro mucho más accesible sobre versiones más simples de lo mismo.

Mi pregunta es sobre las interacciones dentro de este tipo de modelos. ¿Qué pasa si quiero hacer algo como lo siguiente:

$$y = \alpha + f_1(x_1) + f_2(x_2) + f_3(x_1\times x_2) + e_i$$ si estábamos en OLS tierra (donde la $f$ es sólo una beta), no tendría ningún problema con la interpretación f 3 . Si estimamos a través de splines penalizadas, tampoco tengo problemas con la interpretación en el contexto aditivo. $\hat{f}_3$

Pero el paquete MGCV en GAM tiene estas cosas llamadas "suavizadores de productos tensoriales". Busco en Google "producto tensorial" y mis ojos se ven de inmediato tratando de leer las explicaciones que encuentro. O no soy lo suficientemente inteligente o las matemáticas no se explican muy bien, o ambas.

En lugar de codificar

normal = gam(y~s(x1)+s(x2)+s(x1*x2))

un producto tensor haría lo mismo (?) por

what = gam(y~te(x1,x2))

Cuando lo hago

plot(what)

o

vis.gam(what)

Tengo una salida realmente genial. Pero no tengo idea de lo que está sucediendo dentro de la caja negra que es te(), ni cómo interpretar la salida genial mencionada anteriormente. La otra noche tuve la pesadilla de dar un seminario. Les mostré a todos un gráfico genial, me preguntaron qué significaba y no lo sabía. Entonces descubrí que no tenía ropa puesta.

¿Podría alguien ayudarme tanto a mí como a la posteridad, dando un poco de mecánica e intuición sobre lo que está sucediendo debajo del capó aquí? ¿Idealmente al decir un poco sobre la diferencia entre el caso normal de interacción aditiva y el caso tensorial? Puntos de bonificación por decir todo en inglés simple antes de pasar a las matemáticas.

Voy a (intentar) responder esto en tres pasos: primero, identifiquemos exactamente lo que queremos decir con un suave univariante. A continuación, describiremos un suavizado multivariado (específicamente, un suavizado de dos variables). Finalmente, haré mi mejor intento de describir un producto tensorial sin problemas.

1) Univariante suave

Digamos que tenemos algunos datos de respuesta que conjeturamos que es una función desconocida f de una variable predictora x más algún error ε . El modelo sería:$y$$f$$x$$ε$

$$y=f(x)+ε$$

Ahora, para ajustar este modelo, tenemos que identificar la forma funcional de . La forma en que hacemos esto es identificando las funciones básicas, que se superponen para representar la función f en su totalidad. Un ejemplo muy simple es una regresión lineal, en la cual las funciones básicas son solo β 2 x y β 1 , la intersección. Aplicando la base de expansión, tenemos$f$$f$$β_2x$$β_1$

$$y=β_1+β_2x+ε$$

En forma de matriz, tendríamos:

$$Y=Xβ+ε$$

Donde es un vector de columna n-por-1, X es una matriz modelo n-por-2, β es un vector de columna 2-por-1 de coeficientes del modelo, y ε es un vector de errores de columna n-por-1 . X tiene dos columnas porque hay dos términos en nuestra expansión de base: el término lineal y la intersección.$Y$$X$$β$$ε$$X$

El mismo principio se aplica para la expansión de la base en MGCV, aunque las funciones básicas son mucho más sofisticadas. Específicamente, las funciones de base individual no necesitan definirse sobre el dominio completo de la variable independiente . Este suele ser el caso cuando se utilizan bases basadas en nudos (ver "ejemplo basado en nudos"$x$) El modelo se representa como la suma de las funciones básicas, cada una de las cuales se evalúa en cada valor de la variable independiente. Sin embargo, como mencioné, algunas de estas funciones básicas toman un valor de cero fuera de un intervalo dado y, por lo tanto, no contribuyen a la expansión de la base fuera de ese intervalo. Como ejemplo, considere una base spline cúbica en la que cada función base es simétrica respecto a un valor diferente (nudo) de la variable independiente; en otras palabras, cada función base se ve igual pero simplemente se desplaza a lo largo del eje de la variable independiente (esto es una simplificación excesiva, ya que cualquier base práctica también incluirá una intercepción y un término lineal, pero espero que entiendas la idea).

Para ser explícito, una expansión básica de la dimensión podría verse así:$i-2$

$$y=β_1+β_2x+β_3f_1(x)+β_4f_2(x)+...+β_if_{i-2} (x)+ε$$

donde cada función es, quizás, una función cúbica de la variable independiente x .$f$$x$

La ecuación matricial todavía se puede utilizar para representar nuestro modelo. La única diferencia es que X es ahora una matriz n-por-i; es decir, tiene una columna para cada término en la expansión de base (incluido el término de intercepción y lineal). Dado que el proceso de expansión de bases nos ha permitido representar el modelo en forma de una ecuación matricial, podemos usar mínimos cuadrados lineales para ajustar el modelo y encontrar los coeficientes β .$Y=Xβ+ε$$X$$β$

Este es un ejemplo de regresión sin potencializar, y una de las principales fortalezas de MGCV es su estimación de suavidad a través de una matriz de penalización y un parámetro de suavizado. En otras palabras, en lugar de:

$$β=(X^TX)^{-1}X^TY$$

tenemos:

$$β=(X^TX+λS)^{-1}X^TY$$

donde es una matriz de penalización cuadrática i -by- i y λ es un parámetro de suavizado escalar. No entraré en la especificación de la matriz de penalización aquí, pero debería ser suficiente decir que para cualquier expansión de base dada de alguna variable independiente y definición de una penalización de "meneo" cuadrática (por ejemplo, una penalización de segunda derivada), uno puede calcular la matriz de penalización S .$S$$i$$i$$λ$$S$

MGCV puede usar varios medios para estimar el parámetro de suavizado óptimo . No voy a entrar en ese tema ya que mi objetivo aquí era dar una visión general de cómo se construye una suavidad univariante, lo que creo que he hecho.$λ$

2) Multivariante suave

La explicación anterior se puede generalizar a múltiples dimensiones. Volvamos a nuestro modelo que da la respuesta como una función f de predictores x y z . La restricción a dos variables independientes evitará abarrotar la explicación con notación arcana. El modelo es entonces:$y$$f$$x$$z$

$$y=f(x,z)+ε$$

Ahora, debería ser intuitivamente obvio que vamos a representar con una expansión básica (es decir, una superposición de funciones básicas) tal como lo hicimos en el caso univariante de f ( x ) anterior. También debería ser obvio que al menos una, y casi seguramente muchas más, de estas funciones básicas deben ser funciones de x y z (si este no fuera el caso, entonces implícitamente f sería separable de tal manera que f ( x , z ) = f x ( x ) + $f(x,z)$$f(x)$$x$$z$$f$ ). Una ilustración visual de una base spline multidimensional se puede encontraraquí. Una expansión de base bidimensional completa de la dimensión i - 3 podría verse así:$f(x,z)=f_x(x)+f_z(z)$$i-3$

$$y=β_1+β_2x+β_3z+β_4f_1(x,z)+...+β_if_{i-3} (x,z)+ε$$

Creo que está bastante claro que todavía podemos representar esto en forma de matriz con:

$$Y=Xβ+ε$$

simplemente evaluando cada función base en cada combinación única de y z . La solución sigue siendo:$x$$z$

$$β=(X^TX)^{-1}X^TY$$

Calcular la matriz de penalización de la segunda derivada es muy similar a la del caso univariante, excepto que en lugar de integrar la segunda derivada de cada función base con respecto a una sola variable, integramos la suma de todas las segundas derivadas (incluidas las parciales) con respecto a todas las variables independientes. Los detalles de lo anterior no son especialmente importantes: el punto es que aún podemos construir la matriz de penalización y usar el mismo método para obtener el valor óptimo del parámetro de suavizado λ , y dado ese parámetro de suavizado, el vector de coeficientes sigue siendo:$S$$λ$

$$β=(X^TX+λS)^{-1}X^TY$$

Ahora, este suavizado bidimensional tiene una penalización isotrópica : esto significa que se aplica un solo valor de en ambas direcciones. Esto funciona bien cuando tanto x como z están en aproximadamente la misma escala, como una aplicación espacial. Pero, ¿qué sucede si reemplazamos la variable espacial z con la variable temporal t ? Las unidades de t pueden ser mucho más grandes o más pequeñas que las unidades de x , y esto puede alterar la integración de nuestras segundas derivadas porque algunos de esos derivados contribuirán de manera desproporcionada a la integración general (por ejemplo, si medimos t en nanosegundos y $λ$$x$$z$$z$$t$$t$$x$$t$$x$en años luz, la integral de la segunda derivada con respecto a puede ser mucho más grande que la integral de la segunda derivada con respecto a x , y por lo tanto la "ondulación" a lo largo de la dirección x puede quedar en gran parte sin convertir). La diapositiva 15 de la "caja de herramientas suave" que vinculé tiene más detalles sobre este tema.$t$$x$$x$

Vale la pena señalar que no descomponemos las funciones básicas en bases marginales de y z . La implicación aquí es que los suavizados multivariados deben construirse a partir de bases que admitan múltiples variables. El producto tensor suaviza la construcción de bases multivariadas a partir de bases marginales univariadas, como explico a continuación.$x$$z$

3) El producto tensor suaviza

Los productos de tensor suavizan el problema de modelar respuestas a interacciones de múltiples entradas con diferentes unidades. Supongamos que tenemos una respuesta que es una función f de la variable espacial xy la variable temporal t . Nuestro modelo es entonces:$y$$f$$x$$t$

$$y=f(x,t)+ε$$

Lo que nos gustaría hacer es construir una base bidimensional para las variables y t . Esto será mucho más fácil si podemos representar f como:$x$$t$$f$

$$f(x,t)=f_x(x)f_t(t)$$

En un sentido algebraico / analítico, esto no es necesariamente posible. Pero recuerde, estamos discretizar los dominios de y t (imaginar una de dos dimensiones "celosía" definida por las ubicaciones de nudos en la x y t ejes) de tal manera que la "verdadera" función f está representado por la superposición de funciones de base . Así como asumimos que una función univariada muy compleja puede ser aproximada por una función cúbica simple en un intervalo específico de su dominio, podemos suponer que la función no separable f ( x , t ) puede ser aproximada por el producto de funciones más simples f x ( $x$$t$$x$$t$$f$$f(x,t)$ y f t ( t ) en un intervalo, siempre que nuestra elección de dimensiones base haga que esos intervalos sean lo suficientemente pequeños$f_x(x)$$f_t(t)$

Nuestra expansión de base, dada una base -dimensional en x y j -dimensional en t , se vería así:$i$$x$$j$$t$

$$\begin{align} y = &β_{1} + β_{2}x + β_{3}f_{x1}(x)+β_{4}f_{x2}(x)+...+ \\ &β_{i}f_{x(i-3)}(x)+ β_{i+1}t + β_{i+2}tx + β_{i+3}tf_{x1}(x)+β_{i+4}tf_{x2}(x)+...+ \\ &β_{2i}tf_{x(i-3)}(x)+ β_{2i+1}f_{t1}(t) + β_{2i+2}f_{t1}(t)x + β_{2i+3}f_{t1}(t)f_{x1}(x)+β_{i+4}f_{t1}(t)f_{x2}(x){\small +...+} \\ &β_{2i}f_{t1}(t)f_{x(i-3)}(x)+\ldots+ \\ &β_{ij}f_{t(j-3)}(t)f_{x(i-3)}(x) + ε \end{align}$$

$x$$t$$X$$T$$n^2$$ij$ $X \otimes T$$ij$$i$$j$$ij$

NOTA: Me gustaría señalar que, dado que construimos explícitamente las funciones de base de producto tensorial tomando productos de funciones de base marginal, las bases de producto tensorial pueden construirse a partir de bases marginales de cualquier tipo. No necesitan admitir más de una variable, a diferencia del suave multivariado discutido anteriormente.

$ij-i-j+1$$t$$β_{x1}$$j$$x$$β_{t1}$$i$

Entonces podemos representar esto como:

$$y=β_1+β_2x+β_3t+β_4f_1(x,t)+β_5f_2(x,t)+...+β_{ij-i-j+1}f_{ij-i-j-2}(x,t)+ε$$

$f$$x$$t$

$$Y=Xβ+ε$$

Cuál (todavía) tiene la solución:

$$β=(X^TX)^{-1}X^TY$$

$X$$ij-i-j+1$$J_x$$J_t$

$$J_x=β^T I_j \otimes S_x β$$

y,

$$J_t=β^T S_t \otimes I_i β$$

$x$$t$$λ_x$$λ_t$

Recomiendo leer todas las viñetas en el sitio web de MGCV, así como " Modelos aditivos generalizados: e introducción con R ". Larga vida a Simon Wood.