predecir(lm) con dos interacciones ¿es posible o incorrecto?

Question

Para entender la salida de uno de mis modelos lme he producido un ejemplo un poco más simple usando lm (así que no hay factor aleatorio). Me he dado cuenta de que mi modelo ajustado no parece ajustarse correctamente a los datos, ya que los valores y predichos se desvían mucho de los valores dados cuando se utiliza predict(lm).

Mi conjunto de datos es:

a = c(rep(1:10, 4))
b = c(10,20,30,40,50,60,70,80,90,100, 5,8,10,14,17,22,27,35,42,50, 90,82,73,64,56,48,40,33,25,18, 5,6,8,10,12,14,17,20,23,26)
c = c(rep("male", 10), rep("female", 10), rep("male", 10), rep("female", 10))
d = c(rep("low", 20), rep("high", 20))
e = data.frame(yval = b, xval = a, sex = c, education = d)

Gráficamente se ve así:

library(car)
scatterplot(yval~xval | education, smooth = F, grid = T, spread = F, reg.line = T, data = e, xlab = "x", ylab = "y") 
scatterplot(yval~xval | sex, smooth = F, grid = T, spread = F, reg.line = T, data = e, xlab = "x", ylab = "y") 

y el modelo lineal es:

lm2 = lm(yval~xval+sex+xval:sex+education+xval:education, data = e)
summary(lm2)

Utilizando

e_pred = e
e_pred$pred = predict(lm2)

Obtengo los valores predichos que no coinciden en absoluto con los datos reales:

e = cbind(e,e_pred$pred)

¿Se debe esto a que hay más de una interacción significativa?

¡¡Muchas gracias (de antemano) por leer y quizás responder!!

Answer 1

Hay algo muy erróneo en ese modelo en cuanto al ajuste a los datos reales. Siempre es una buena idea mirar los gráficos de los aspectos del ajuste del modelo, en particular en este caso para los patrones en los residuos y las observaciones influyentes.

Para ver varios gráficos de diagnóstico, haga

layout(matrix(1:4, ncol = 2))
plot(lm2)
layout(1)

Esto produce el siguiente gráfico

Los patrones claros en los residuos sugieren que falta una covariable o un término en el modelo. Como los residuos no son independientes, eso pone en duda las inferencias de la tabla producida por summary() . Obsérvese también el gráfico inferior derecho, en el que la mayoría de sus observaciones se sitúan fuera de las líneas rojas (discontinuas), lo que indica una fuerte influencia y apalancamiento en los coeficientes del modelo.

Ahora tenemos que entender qué puede estar causando esto. Si reutilizamos el scatterplot pero con resid(lm2) como variable del eje Y, podemos examinar cómo varían los residuos con xval condicionada a los dos factores

scatterplot(resid(lm2) ~ xval | education, smooth = F, grid = T, spread = F, 
            reg.line = T, data = e)
scatterplot(resid(lm2) ~ xval | sex, smooth = F, grid = T, spread = F, 
            reg.line = T, data = e)

A continuación se muestran los dos gráficos, el primero condicionado por education

y el siguiente por sex

Esto demuestra que el efecto de xval en yval varía considerablemente (cambio de pendiente) en función de los conjuntos de observaciones que consideremos. Si se trazan los residuos frente a sex o contra education (como boxplots, no mostrados aquí) no se ven tan mal, así que el problema parece ser con xval .

Como menciona que está haciendo esto para ver los resultados de un modelo mixto, se supone que las observaciones están agrupadas por sujetos. Esto podría explicar las aparentes diferencias en la pendiente de la relación para xval en yval donde diferentes individuos tenían diferentes intercepciones y pendientes.

Para ello, debe dibujar gráficos de los residuos condicionados a Subject .

Si no hay Subject en los datos, entonces el patrón de los residuos sugiere que hay otro factor no medido (o no proporcionado aquí) que está modificando la relación entre yval y xval .

Pero para responder a la pregunta concreta, hay efectos en los datos que no se tienen en cuenta en el modelo lineal simple que se ha ajustado. Una opción que hay que seguir es si faltan interacciones de orden superior. Por ejemplo, si permitimos interacciones de tres vías de las variables en el modelo, se observa una fuerte interacción de tres vías y el xval término es ahora significativo.

lm4 <- lm(yval ~ (xval + sex + education)^3, data = e)

> summary(lm4)

Call:
lm(formula = yval ~ (xval + sex + education)^3, data = e)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5394 -0.7318  0.0000  0.3379  4.8545 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 1.0667     1.2218   0.873  0.38915    
xval                        2.3697     0.1969  12.034 2.03e-13 ***
sexmale                    96.0667     1.7279  55.598  < 2e-16 ***
educationlow               -5.1333     1.7279  -2.971  0.00559 ** 
xval:sexmale              -10.4121     0.2785 -37.390  < 2e-16 ***
xval:educationlow           2.5515     0.2785   9.162 1.84e-10 ***
sexmale:educationlow      -92.0000     2.4436 -37.649  < 2e-16 ***
xval:sexmale:educationlow  15.4909     0.3938  39.335  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.789 on 32 degrees of freedom
Multiple R-squared:  0.9965,    Adjusted R-squared:  0.9957 
F-statistic:  1295 on 7 and 32 DF,  p-value: < 2.2e-16

Los datos se ajustan ahora casi perfectamente, todos los términos del modelo son significativos y los gráficos de diagnóstico del modelo se ven mucho mejor, aunque todavía hay algunas observaciones influyentes y la distribución de los residuos es un poco más pesada en las colas:

layout(matrix(1:4, ncol = 2))
plot(lm4)
layout(1)

que produce

Así que no es perfecto, pero quizás ilustra bastante bien el problema de los predictores o términos que faltan en el modelo original.