¿Resultados de la función lm() con una variable categórica ordenada dependiente?

Question

Intento comprender La respuesta de Ben Bolker a esta pregunta.

En primer lugar, creamos un marco de datos:

set.seed(101)
d <- data.frame(x=sample(1:4,size=30,replace=TRUE))
d$y <- rnorm(30,1+2*d$x,sd=0.01)

Entonces el Sr. Bolker dice:

x como factor ordenado

coef(lm(y~ordered(x),d))
##  (Intercept) ordered(x).L ordered(x).Q ordered(x).C 
##  5.998121421  4.472505514  0.006109021 -0.003125958 

Ahora el intercepto especifica el valor de y en el media nivel de factor (a medio camino entre 2 y 3); el L (lineal) da una medida de la tendencia lineal (no muy seguro que puedo explicar el valor particular ...), Q y C especifique los términos cuadrático y cúbico (que son próximos a cero en este caso porque el patrón es lineal); si hubiera más niveles, los contrastes de orden superior se numerarían 5, 6, ...

---

Mi pregunta es, ¿cómo es explícitamente la fórmula de regresión?

Pensaba lm() hace un modelo como este:

y = 5,9981 + 4,4725 (x_1) + 0,0061 (x_2) - 0,00312 (x_3)

donde, dado que el x_i son categorías, sólo pueden ser 0 ó 1.

No entiendo qué tienen que ver los términos cuadrático y cúbico con un modelo lineal. Aun así, elevar al cuadrado/cubo cualquiera de las variables no supondría ninguna diferencia, ya que 0 ^ 3= 0 y 1^3 = 1.

Answer 1

En este caso, ¿qué lm() está haciendo es convertir su variable "categórica" en una secuencia numérica en orden.

Para que esto quede más claro, adaptaré un poco el código de Bolker para que la variable X sea más obviamente categórica:

set.seed(101)
d <- data.frame(x=sample(1:4,size=30,replace=TRUE))
d$y <- rnorm(30,1+2*d$x,sd=0.01)
d$x = factor(d$x, labels=c("none", "some", "more", "a lot"))
coef(lm(y~x, d))
#  (Intercept)       xsome       xmore      xa lot 
#     3.001627    1.991260    3.995619    5.999098

En este caso, la media de x=None está en el intercepto, y se indica la desviación con respecto a éste para cada categoría.

coef(lm(y~ordered(x), d))
#  (Intercept) ordered(x).L ordered(x).Q ordered(x).C 
#  5.998121421  4.472505514  0.006109021 -0.003125958

Conceptualmente, lo que ha ocurrido aquí es que el ordered() función convertida x en newx utilizando (algo similar a):

if (x=="None") newx=-.67
if (x=="some") newx=-.22
if (x=="more") newx=.22
if (x=="a lot") newx=.67

y luego ajustó (algo así como) el modelo: $$y = a + b_0 \times newx + b_1 \times newx^2 + b_2 \times newx^3$$

donde tiene lineal $newx$ cuadrática $newx^2$ y cúbico $newx^3$ componentes.

Nota, he dicho que es algo así como que, porque el problema con el modelo descrito allí es que $newx$ , $newx^2$ y $newx^3$ no son en absoluto independientes. Lo que lm() utiliza un conjunto de contrastes generados por contr.poly(4) . Estos contrastes garantizan la ortogonalidad, de modo que los componentes lineal, cuadrático y cúbico son independientes. Pero el principio es similar - cuando se ajustan factores ordenados, lm() se ajusta a un componente lineal, cuadrático, cúbico, etc...

Puede comprobarlo comparando:

coef(lm(y~ordered(x), d))
#  (Intercept) ordered(x).L ordered(x).Q ordered(x).C 
#  5.998121421  4.472505514  0.006109021 -0.003125958

con

contrasts(d$x) <- contr.poly(4)
coef(lm(y~x, d))
#  (Intercept) ordered(x).L ordered(x).Q ordered(x).C 
#  5.998121421  4.472505514  0.006109021 -0.003125958

Exactamente idéntico. Así que si quieres una comprensión más completa de lo que pasó, echa un vistazo más de cerca a contr.poly() y contrastes polinómicos ortogonales en general.

Hay que tener en cuenta que aquí hay una suposición implícita: se supone que la diferencia entre cada dos niveles es igual. Así, "ninguno" está tan lejos de "algunos" como "algunos" de "más", y "más" de "muchos".