Regresión simple en R

Pasos para ajustar un modelo de regresión simple

Variable de respuesta numérica vs. Variable explicativa numérica

En esta entrada me tomaré el tiempo para explicar cómo realizar un análisis de regresión lineal simple usando R.

1. Para n pares de observaciones \( (x_i,y_i) \), con \( i=1,2,.,n \), se construye un diagrama de dispersión entre X vs. Y. Este gráfico dará idea del tipo de relación entre ambas variables: si es directamente proporcional o inversamente proporcional, y si es fuerte (nube de puntos muy estrecha) o débil (nube de puntos muy dispersa).

2. Hallar la recta que atraviesa la nube de puntos del gráfico anterior. Esto consiste estimar los parámetros intercepto y pendiente. La ecuación resultante es lo que se conoce como “el modelo” y tiene la forma \( E[Y|x]=\beta_0+\beta_1 x \). La manera de estimar los parámetros del modelo es a partir de las \( n \) observaciones \( (x_i,y_i) \) utilizando las ecuaciones de mínimos cuadrados dadas por:

   Con esto se obtiene el modelo ajustado o estimado \( \bar y = \bar\beta_0+ \bar \beta_1 x \) .

3. Antes de utilizar este modelo (es decir, antes de realizar inferencias con él) se debe probar su validez. Esto se logra probando la 1) normalidad, 2) homocedasticidad e 3) independencia de los residuales (denotados con la letra e). Todo modelo cuenta con n residuales; cada uno de ellos se calcula como la diferencia entre el valor observado de la variable respuesta y el valor ajustado: \( e_i=y_i- \bar y_i \). Así, para saber si el modelo es válido se deben probar, estadísticamente, los tres requisitos mencionados.

   3.1 Normalidad de los residuales. La hipótesis nula de esta prueba estadística es Ho: los residuales se distribuyen normal. Se prueba con el test de Shapiro-Wilk, o el de Kolmogorov-Smirnov, o el de Kramer, o el de Chi-cuadrado, etc. La idea es no rechazar \( H_o \).

   3.2. Homocedasticidad de los residuales. La hipótesis nula de esta prueba estadística es \( H_o \): la varianza de los residuales es constante a lo largo de los valores de x. También existen varias formas de probar esta hipótesis, una de ellas es utilizando la prueba de Harrison-McCabe o Bartlett. La idea es no rechazar \( H_o \).

   3.3. Independencia de los residuales. La hipótesis nula de esta prueba estadística es \( H_o \): los residuales son independientes entre sí. También existen varias pruebas y una de ellas es la Durbin-Watson. La idea es no rechazar \( H_o \).

   En los tres casos es muy útil realizar gráficos que vislumbren el comportamiento de los residuales.

4. De incumplirse los requisitos de normalidad u homocedasticidad, se recurren a transformaciones y se comienza desde un principio. Una vez se pruebe que el modelo es válido, se procede a hacer uso de éste. Ello consiste en probar la hipótesis nula de que \( \beta_1=0 \).

5. De aquí, se procede a concluir sobre la relación entre la variable respuesta y la explicativa y realizar estimaciones de la variable respuesta.

6. Para concluir sobre qué tan fuerte o débil es la relación, utilizamos el coeficiente de correlación de Pearson, cuyo valor absoluto se obtiene extrayendo la raíz cuadrada del coeficiente de determinación.

Pasos en R

A continuación se describen las funciones utilizadas en los pasos para ajustar un modelo de regresión simple en R.

1.Datos y diagrama de dispersión

library(ggplot2)
temperatura <- c(88.6, 71.6, 93.3, 84.3, 80.6, 75.2, 69.7, 71.6, 69.4)
chirridos <- c(20, 16, 19.8, 18.4, 17.1, 15.5, 14.7, 15.7, 15.5)
qplot(temperatura, chirridos, xlab = "Temperatura", ylab = "Chirridos") + geom_smooth(method = "loess") + 
    theme_bw()

2.Ajuste del modelo

modelo <- lm(chirridos ~ temperatura)

3.Chequeo de supuestos

residuales <- modelo$residuals

3.1 Normalidad

shapiro.test(residuales)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.9472, p-value = 0.6591

qplot(sample = residuales, stat = "qq") + theme_bw()

Basado en el valor P arrojado por la prueba de Shapiro - WilK P=0.6591, no se rechaza \( H_o \), y se concluye que los datos sigue una distribución normal.

3.2 Homocedasticidad

library(lmtest)

## Warning: package 'lmtest' was built under R version 3.0.2

## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following object is masked from 'package:base':
## 
##     as.Date, as.Date.numeric

hmctest(modelo)

## 
##  Harrison-McCabe test
## 
## data:  modelo
## HMC = 0.4679, p-value = 0.556

Usando la prueba Harrison-McCabe obtenemos un P=0.52, con cual concluimos que la varianza de los residuales es constante a lo largo de los valores de x. Una vez asegurado el cumplimiento de los supuestos, el modelo ajustado es válido y puede ser usado para hacer inferencias.

4. ANOVA

anova(modelo)

## Analysis of Variance Table
## 
## Response: chirridos
##             Df Sum Sq Mean Sq F value Pr(>F)    
## temperatura  1  29.04   29.04    90.6  3e-05 ***
## Residuals    7   2.24    0.32                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.1 Cálculo de coeficientes

coef(modelo)

## (Intercept) temperatura 
##      0.0981      0.2156

qplot(temperatura, chirridos, xlab = "Temperatura", ylab = "Chirridos") + theme_bw() + 
    geom_abline(intercept = 0.09809541, slope = 0.21555749)

4.2 Cálculo coeficiente de correlación

cor(temperatura, chirridos, use = "everything", method = "pearson")

## [1] 0.9635

qplot(temperatura, chirridos, xlab = "Temperatura", ylab = "Chirridos") + theme_bw() + 
    geom_abline(intercept = 0.09809541, slope = 0.21555749) + geom_text(aes(75, 
    19, label = "Pearson-R = 0.96"))

5. Conclusiones

Basado en los resultados podemos ver claramente, que las dos varialbes muestra una correlación creciente y que con el aumento de la temperatura se aumenta los chirridos significativamente. Este modelo además de darmos la correlación, es válido para predecir los chirridos en función de la temperatura.