lunes, 12 de marzo de 2018

Regresión Lineal Simple

Regresión Lineal Simple

1. Introducción

Un modelo de regresión lineal simple trata de analizar cuál es el comportamiento de una variable llamada dependiente como función lineal de una variable explicativa o independiente.

Para el conjunto de la población, el modelo de regresión se representaría así:

\(y_{i} = \beta_{0} + \beta_{1}x_{i} + \epsilon_{i}\)

donde y es la variable dependiente y x es la variable independiente. A su vez, \(\epsilon\) es el error o perturbación aleatoria.

El objetivo es estimar los valores de \(\beta_{0}\) y \(\beta_{1}\) a partir de los datos de una muestra determinada de forma que el modelo obtenido sería de la forma:

\(\hat{y}_{i} = \hat{\beta}_{0} + \hat{\beta}_{1}x_{i}\)

La diferencia entre \(y_{i}\) y \(\hat{y}_{i}\) para cada valor de \(x_{i}\) es el residuo \(\epsilon_{i}\), que se supone una variable aleatoria independiente de x y de y.

2. Estimación de un modelo de Regresión Lineal Simple

2.1. Obtención de los datos

Utilizaremos el data.set “Prestige”, incluido en el paquete {car}

# Cargamos el paquete correspondiente:
library (car)
# Incorporamos el data.set "Prestige":
data ("Prestige")

Este data.set recoge datos referidos al Prestige of Canadian Occupations. Con la función help repasamos la descripción, variables y origen de los datos. Las variables son:

  • education: Average education of occupational incumbents, years, in 1971.
  • income: Average income of incumbents, dollars, in 1971.
  • women: Percentage of incumbents who are women.
  • prestige: Pineo-Porter prestige score for occupation, from a social survey conducted in the mid-1960s.
  • census: Canadian Census occupational code.
  • type: Type of occupation. A factor with levels (note: out of order): bc, Blue Collar; prof, Professional, Managerial, and Technical; wc, White Collar.

La estructura del fichero puede comprobarse con la función str:

str (Prestige)
## 'data.frame':    102 obs. of  6 variables:
##  $ education: num  13.1 12.3 12.8 11.4 14.6 ...
##  $ income   : int  12351 25879 9271 8865 8403 11030 8258 14163 11377 11023 ...
##  $ women    : num  11.16 4.02 15.7 9.11 11.68 ...
##  $ prestige : num  68.8 69.1 63.4 56.8 73.5 77.6 72.6 78.1 73.1 68.8 ...
##  $ census   : int  1113 1130 1171 1175 2111 2113 2133 2141 2143 2153 ...
##  $ type     : Factor w/ 3 levels "bc","prof","wc": 2 2 2 2 2 2 2 2 2 2 ...

2.2. Exploración de la relación entre variables

Creamos un “scatterplot” (gráfico de puntos) para explorar visualmente la relación entre prestige y education

Primero cargamos los paquetes necesarios:

library (magrittr)  # Para el uso de la función "pipe" (%>%)
library (ggplot2)   # Para la creación de gráficos

Dibujamos el gráfico:

Prestige %>%
  ggplot(aes(x = education, y = prestige)) +
  geom_point(colour = "red")

La relación puede ser cuantificada usando el coeficiente de correlación de Pearson:

cor (Prestige$education, Prestige$prestige)
## [1] 0.8501769

Una información mucho más completa con:

cor.test (Prestige$education, Prestige$prestige)
## 
##  Pearson's product-moment correlation
## 
## data:  Prestige$education and Prestige$prestige
## t = 16.148, df = 100, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7855899 0.8964367
## sample estimates:
##       cor 
## 0.8501769

Volvemos a dibujar el scatterplot, añadiendo ahora la recta de regresión:

Prestige %>%
  ggplot (aes (x = education, y = prestige)) +
  geom_point (colour = "red") +
  geom_smooth (method = lm)

Como puede observarse a lo largo del análisis anterior, la correlación entre ambas variables es muy alta (cor = 0,85) y con un alto nivel de confianza (p-value < 2.2e-16). La última gráfica señala además el intervalo de confianza (95%) en cada tramo de la recta.

2.3. Creación del modelo de Regresión Lineal Simple

Creamos el modelo de regresión lineal simple núm. 1:

rls1 <- lm (prestige ~ education, data = Prestige)
# La pendiente de la recta y el punto de corte con el eje "y" son los coeficientes:
rls1$coefficients
## (Intercept)   education 
##  -10.731982    5.360878
# Con la función summary obtenemos la información más importante:
summary (rls1)
## 
## Call:
## lm(formula = prestige ~ education, data = Prestige)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.0397  -6.5228   0.6611   6.7430  18.1636 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -10.732      3.677  -2.919  0.00434 ** 
## education      5.361      0.332  16.148  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared:  0.7228, Adjusted R-squared:   0.72 
## F-statistic: 260.8 on 1 and 100 DF,  p-value: < 2.2e-16

Algunos resultados obtenidos

  • El p-value de 2.2e-16 indica que la relación es estadísticamente significativa a nivel de p < 0.001
  • El multiple R-squared value (R-squared) de 0.7278 indica la varianza explicada por el modelo y puede usarse como medida de la fuerza predictiva (en ausencia de overfitting)
  • Por cada aumento de una unidad en “education” hay un incremento de 5.361 en “prestige”. Por tanto, el prestigio se incrementa con la educación

También podemos examinar de forma gráfica la relación entre education y prestige en función del tipo de ocupación (“type”):

ggplot (Prestige) + 
  aes (x = education, y = prestige, color = type) + 
  geom_point() + 
  scale_colour_hue (l = 50) + 
  geom_smooth (method = lm, se = FALSE)

3. Evaluación del modelo (bondad del ajuste)

3.1. Representación de los residuos

Pueden dibujarse los residuos respecto a la recta de regresión. Volvemos a dibujar el gráfico de dispersión con la recta de regresión y añadimos las líneas verticales que unen los valores registrados con los valores ajustados:

plot (Prestige$education, Prestige$prestige,
  pch = 16,
  cex = 0.8,
  col = "red",
  main = "Relación prestigio y educación",
  xlab = "Educación (años)",
  ylab = "Prestigio (escala 0-100)")
abline (rls1, lwd = 2)
# Creamos un objeto con los valores que predice el modelo:
pred <- predict (rls1)
segments (Prestige$education,
  Prestige$prestige,
  Prestige$education,
  pred,
  col="blue")

Las líneas de color azul representan los residuos del modelo, es decir la diferencia entre \(y_{i}\) y \(\hat{y}_{i}\) para cada valor de \(x_{i}\).

3.2. Examen gráfico de los residuos

Ahora examinamos la distribución de los residuos con objeto de comprobar si siguen alguna pauta concreta o se pueden considerar aleatorios y si su distribución estadística es normla o sesgada.

# Dibujamos los valores ajustados y los residuos, para observar su distribución, añadiendo una línea horizontal en y = 0:
ggplot (rls1) + aes (x = fitted (rls1), y = resid (rls1)) +
  geom_point (colour = "blue") + 
  geom_hline (yintercept = 0)

# Dibujamos el gráfico Q-Q de los residuos para observar su distribución normal:
qqnorm (resid (rls1))
qqline (resid (rls1))

Tal y como se observa por los gráficos trazados, la distribución de los residuos no presenta ninguna pauta apreciable, por lo que cabe suponerlos totalmente aleatorios y con una distribución que se aproxima a la de la curva normal.

3.3. Examen estadístico de los residuos

Aunque la inspección gráfica suele darse como suficiente cuando las pautas (o no pautas) aparecen claras, nunca está de más aplicar las pruebas de hipótesis para comprobar que el modelo cumple con las especificaciones de validez del modelo de regresión lineal.

3.3.1. Test de normalidad

La primera comprobación es que los residuos presentan una distribución normal. Como el número de observaciones es menor de 2.000 el test a aplicar sería el de Shapiro-Wilk:

shapiro.test (resid (rls1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(rls1)
## W = 0.98065, p-value = 0.1406

Al obtener un p-value > 0.05 no se puede rechazar la hipótesis nula (\(H_{0}\)) de que los residuos presentan una distribución normal.

3.3.2. Test de varianza constante u homocedasticidad

Mediante este test se comprueba que la varianza es constante a lo largo de todo el rango de valores de las observaciones. El test es el de contraste de heterocedasticidad de Breusch-Pagan, utilizando la función “bptest” del paquete {lmtest} aplicada al modelo de regresión:

lmtest::bptest (rls1)
## 
##  studentized Breusch-Pagan test
## 
## data:  rls1
## BP = 0.70577, df = 1, p-value = 0.4009

El estadístico de Breusch-Pagan arroja un valor de 0.71 con un p-valor de 0.40, lo que supone aceptar la hipótesis nula de varianza constante.

3.3.3. Test de ausencia de autocorrelación en las perturbaciones

El test de Durbin-Watson construye un estadístico con un rango de valores entre 0 y 4. Toma el valor 2 cuando no hay autocorrelación en los residuos y los valores extremos 0 y 4 denotan una correlación total positiva o negativa respectivamente. Utilizamos la función “durbinWatsonTest” del paquete {car} aplicada al modelo de regresión:

car::durbinWatsonTest (rls1, max.lag = 4) # "max.lag" compara la autocorrelaciones de orden 1 a 4
##  lag Autocorrelation D-W Statistic p-value
##    1       0.2752512      1.439170   0.006
##    2       0.3005009      1.364598   0.002
##    3       0.3488403      1.256238   0.000
##    4       0.1593670      1.624994   0.086
##  Alternative hypothesis: rho[lag] != 0

Por lo general valores de D-W Statistic por encima de 1 no suponen motivos de alarma, aunque en este caso los p-values tan bajos dejan una cierta duda en el aire.

No hay comentarios:

Publicar un comentario