¿Cómo incluir un término lineal y cuadrático cuando también se incluye la interacción con esas variables?

Al agregar un predictor numérico con predictores categóricos y sus interacciones, generalmente se considera necesario centrar las variables en 0 de antemano. El razonamiento es que los efectos principales son difíciles de interpretar ya que se evalúan con el predictor numérico en 0.

Mi pregunta ahora es cómo centrar si uno no solo incluye la variable numérica original (como un término lineal) sino también el término cuadrático de esta variable. Aquí, dos enfoques diferentes son necesarios:

1. Centrando ambas variables en su media individual. Esto tiene el inconveniente de que el 0 ahora está en una posición diferente para ambas variables considerando la variable original.
2. Centrar ambas variables en la media de la variable original (es decir, restar la media de la variable original para el término lineal y restar el cuadrado de la media de la variable original del término cuadrático). Con este enfoque, el 0 representaría el mismo valor de la variable original, pero la variable cuadrática no estaría centrada en 0 (es decir, la media de la variable no sería 0).

Creo que el enfoque 2 parece razonable dada la razón para centrarse después de todo. Sin embargo, no puedo encontrar nada al respecto (tampoco en las preguntas relacionadas: una , y b ).

¿O es generalmente una mala idea incluir términos lineales y cuadráticos y sus interacciones con otras variables en un modelo?

Cuando se incluyen polinomios e interacciones entre ellos, la multicolinealidad puede ser un gran problema; Un enfoque es observar los polinomios ortogonales.

Generalmente, los polinomios ortogonales son una familia de polinomios que son ortogonales con respecto a algún producto interno.

Así, por ejemplo, en el caso de polinomios sobre alguna región con función de peso $w$, el producto interno es $\int_a^bw(x)p_m(x)p_n(x)dx$ - la ortogonalidad hace que ese producto interno $0$ a no ser que $m=n$.

El ejemplo más simple para polinomios continuos son los polinomios de Legendre, que tienen una función de peso constante durante un intervalo real finito (comúnmente sobre $[-1,1]$)

En nuestro caso, el espacio (las observaciones mismas) es discreto, y nuestra función de peso también es constante (por lo general), por lo que los polinomios ortogonales son una especie de equivalente discreto de los polinomios de Legendre. Con la constante incluida en nuestros predictores, el producto interno es simplemente$p_m(x)^Tp_n(x) = \sum_i p_m(x_i)p_n(x_i)$.

Por ejemplo, considere $x = 1,2,3,4,5$

Comience con la columna constante, $p_0(x) = x^0 = 1$. El siguiente polinomio es de la forma$ax-b$, pero no estamos preocupados por la escala en este momento, así que $p_1(x) = x-\bar x = x-3$. El próximo polinomio sería de la forma$ax^2+bx+c$; Resulta que$p_2(x)=(x-3)^2-2 = x^2-6x+7$ es ortogonal a los dos anteriores:

x         p0  p1  p2   
1          1  -2   2   
2          1  -1  -1
3          1   0  -2
4          1   1  -1
5          1   2   2

Con frecuencia, la base también se normaliza (produciendo una familia ortonormal), es decir, las sumas de cuadrados de cada término se configuran como constantes (por ejemplo, para $n$, o para $n-1$, de modo que la desviación estándar es 1, o quizás con mayor frecuencia, a $1$)

Las formas de ortogonalizar un conjunto de predictores polinomiales incluyen la ortogonalización de Gram-Schmidt y la descomposición de Cholesky, aunque existen muchos otros enfoques.

Algunas de las ventajas de los polinomios ortogonales:

1) la multicolinealidad no es un problema: estos predictores son todos ortogonales.

2) Los coeficientes de bajo orden no cambian a medida que agrega términos . Si te queda un grado$k$ polinomio a través de polinomios ortogonales, conoce los coeficientes de un ajuste de todos los polinomios de orden inferior sin volver a ajustar.

Ejemplo en R ( carsdatos, distancias de frenado contra velocidad):

Aquí consideramos la posibilidad de que un modelo cuadrático sea adecuado:

R usa la polyfunción para configurar predictores polinomiales ortogonales:

> p <- model.matrix(dist~poly(speed,2),cars)
> cbind(head(cars),head(p))
  speed dist (Intercept) poly(speed, 2)1 poly(speed, 2)2
1     4    2           1      -0.3079956      0.41625480
2     4   10           1      -0.3079956      0.41625480
3     7    4           1      -0.2269442      0.16583013
4     7   22           1      -0.2269442      0.16583013
5     8   16           1      -0.1999270      0.09974267
6     9   10           1      -0.1729098      0.04234892

Son ortogonales:

> round(crossprod(p),9)
                (Intercept) poly(speed, 2)1 poly(speed, 2)2
(Intercept)              50               0               0
poly(speed, 2)1           0               1               0
poly(speed, 2)2           0               0               1

Aquí hay una gráfica de los polinomios:

Aquí está la salida del modelo lineal:

> summary(carsp)

Call:
lm(formula = dist ~ poly(speed, 2), data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-28.720  -9.184  -3.188   4.628  45.152 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)       42.980      2.146  20.026  < 2e-16 ***
poly(speed, 2)1  145.552     15.176   9.591 1.21e-12 ***
poly(speed, 2)2   22.996     15.176   1.515    0.136    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.18 on 47 degrees of freedom
Multiple R-squared:  0.6673,    Adjusted R-squared:  0.6532 
F-statistic: 47.14 on 2 and 47 DF,  p-value: 5.852e-12

Aquí hay una gráfica del ajuste cuadrático:

No creo que el centrado valga la pena, y el centrado hace que la interpretación de las estimaciones de parámetros sea más compleja. Si utiliza un software moderno de álgebra matricial, la colinealidad algebraica no es un problema. Su motivación original de centrarse para poder interpretar los efectos principales en presencia de interacción no es fuerte. Los efectos principales cuando se estiman en cualquier valor elegido automáticamente de un factor de interacción continua son algo arbitrarios, y es mejor pensar en esto como un simple problema de estimación al comparar los valores pronosticados. En el rmspaquete Rcontrast.rmsfunción, por ejemplo, puede obtener cualquier contraste de interés independiente de las codificaciones variables. Aquí hay un ejemplo de una variable categórica x1 con niveles "a" "b" "c" y una variable continua x2, ajustada usando una spline cúbica restringida con 4 nudos predeterminados. Se permiten diferentes relaciones entre x2 e y para diferentes x1. Dos de los niveles de x1 se comparan en x2 = 10.

require(rms)
dd <- datadist(x1, x2); options(datadist='dd')
f <- ols(y ~ x1 * rcs(x2,4))
contrast(f, list(x1='b', x2=10), list(x1='c', x2=10))
# Now get all comparisons with c:
contrast(f, list(x1=c('a','b'), x2=10), list(x1='c', x2=10))
# add type ='joint' to get a 2 d.f. test, or conf.type='simultaneous'
# to get simultaneous individual confidence intervals

Con este enfoque, también puede estimar fácilmente los contrastes en varios valores de los factores que interactúan, p. Ej.

contrast(f, list(x1='b', x2=10:20), list(x1='c', x2=10:20))