Interpretación de gráficos de residuos de regresión lineal con R

Question

Estoy investigando si existe una relación entre el día de la semana y un valor de resultado mediante regresión lineal en R y me gustaría saber cómo interpretar los gráficos de residuos.

Datos

Ejemplo de datos ficticios (la media y la DE se basan en datos reales que tengo):

set.seed(14)
mon <- data.frame(id=seq(6, 60*7, by=7), value = rnorm(60, 4372, 145))
tue <- data.frame(id=seq(7, 60*7, by=7), value = rnorm(60, 4433, 206))
wed <- data.frame(id=seq(1, 60*7, by=7), value = rnorm(60, 4671, 143))
thu <- data.frame(id=seq(2, 60*7, by=7), value = rnorm(60, 4555, 154))
fri <- data.frame(id=seq(3, 60*7, by=7), value = rnorm(60, 4268, 149))
sat <- data.frame(id=seq(4, 60*7, by=7), value = rnorm(60, 1579, 110))
sun <- data.frame(id=seq(5, 60*7, by=7), value = rnorm(60, 1136, 68))
startdate <- seq.Date(as.Date("2014-01-01"), by="day", length.out=(60*7) )
id <- seq(1, 60*7)
wd <- weekdays(startdate)
df <- data.frame(id, startdate, wd)
days <- rbind(mon, tue, wed, thu, fri, sat, sun)
df <- merge(df, days)

head(df)
  id  startdate        wd    value
1  1 2014-01-01 Wednesday 4593.117
2  2 2014-01-02  Thursday 4686.159
3  3 2014-01-03    Friday 4352.982
4  4 2014-01-04  Saturday 1825.172
5  5 2014-01-05    Sunday 1206.759
6  6 2014-01-06    Monday 4276.032

que se parece a

library(ggplot2)
ggplot(data=df, aes(x=startdate, y=value, colour=wd)) +
  geom_point() +
  geom_smooth( alpha=.3, size=1, aes(fill=wd)) +
  facet_wrap(~wd) 

Modelo

Modelización de los datos mediante fit <- lm(data=df, value ~ wd) produce los coeficientes:

summary(fit)
....
Coefficients:
             Estimate Std. Error  t value Pr(>|t|)    
(Intercept)  4242.60      18.83  225.319  < 2e-16 ***
wdMonday      148.93      26.63    5.593 4.07e-08 ***
wdSaturday  -2661.93      26.63  -99.965  < 2e-16 ***
wdSunday    -3113.78      26.63 -116.933  < 2e-16 ***
wdThursday    299.65      26.63   11.253  < 2e-16 ***
wdTuesday     189.04      26.63    7.099 5.51e-12 ***
wdWednesday   412.52      26.63   15.492  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 145.9 on 413 degrees of freedom
Multiple R-squared:  0.9896,    Adjusted R-squared:  0.9894 
F-statistic:  6539 on 6 and 413 DF,  p-value: < 2.2e-16

El gráfico de los datos y los coeficientes parecen indicar que existe una relación entre el día de la semana y el valor del resultado.

Sin embargo, sé que también debo tener en cuenta los gráficos de residuos a la hora de interpretar la validez de un modelo. Para este ejemplo los gráficos residuales son:

par(mfrow=c(2,2))
plot(fit)

Pregunta

A través de diversos cursos/universidades/investigaciones sobre estadística (por ejemplo esta pregunta ) Sé que para un buen modelo lineal se buscan residuos homocedásticos insesgados. Pero mis conocimientos sobre este tema están un poco oxidados. Por lo tanto, ¿sugieren mis residuos que un modelo lineal no se ajusta adecuadamente a los datos? ¿O hay algún otro aspecto que debería tener en cuenta, o es que no he entendido nada?

Answer 1

A juzgar por sus gráficos de diagnóstico, yo no me preocuparía por la heteroscedasticidad y utilizaría este modelo.

Sin embargo, puede probar si la modelización de la varianza de los residuos mejoraría el modelo:

library(nlme)
fit1 <- gls(value ~ wd, data=df)

fit2 <- gls(value ~ wd, data=df, weights = varIdent(form = ~ 1 | wd))

anova(fit1, fit2)

#     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
#fit1     1  8 5332.321 5364.509 -2658.160                        
#fit2     2 14 5286.911 5343.239 -2629.456 1 vs 2 57.41009  <.0001

Como ve, el modelo que incluye parámetros de varianza parece ligeramente mejor. Sin embargo, los errores estándar o los valores t sólo difieren ligeramente de su modelo original y no inferiría conclusiones diferentes:

summary(fit2)
#Generalized least squares fit by REML
#  Model: value ~ wd 
#  Data: df 
#       AIC      BIC    logLik
#  5286.911 5343.239 -2629.455
#
#Variance function:
# Structure: Different standard deviations per stratum
# Formula: ~1 | wd 
# Parameter estimates:
#Wednesday  Thursday    Friday  Saturday    Sunday    Monday   Tuesday 
#1.0000000 1.2789517 1.2876130 0.9901914 0.5975224 1.1128592 1.5904740 
#
#Coefficients:
#                Value Std.Error    t-value p-value
#(Intercept)  4242.602  20.92298  202.77237       0
#wdMonday      148.931  27.65462    5.38538       0
#wdSaturday  -2661.925  26.39433 -100.85215       0
#wdSunday    -3113.783  23.06607 -134.99409       0
#wdThursday    299.647  29.49021   10.16091       0
#wdTuesday     189.044  33.25205    5.68520       0
#wdWednesday   412.524  26.49180   15.57175       0

anova(fit)
#Analysis of Variance Table
#
#Response: value
#           Df    Sum Sq   Mean Sq F value    Pr(>F)    
#wd          6 834554455 139092409  6538.6 < 2.2e-16 ***
#Residuals 413   8785585     21273 

anova(fit2)
#Denom. DF: 413 
#            numDF  F-value p-value
#(Intercept)     1 210610.5  <.0001
#wd              6  11982.3  <.0001

Podría modelar los días de la semana como un factor ordenado, lo que haría que las funciones de regresión utilizaran contrastes polinómicos.