Ejemplo de regresión simple

En este ejemplo se trabajará con la base precargada, “state.x77”, la cual contiene 50 renglones correspondientes a información económica de cada uno estados que conforman la Unión Americana, con 8 columnas que hacen referencia a:

Estos datos pueden ser visualizados mediante el siguiente código

library(DT)
## Warning: package 'DT' was built under R version 4.2.2
datos = state.x77
colnames(datos) = c("Population", 'Income', 'Illiteracy', 'Life_Exp', 'Murder', 'HS_Grad', 'Frost', 'Area')
datos = data.frame(datos)
datatable(datos)

Se desea analizar la relación lineal que existe entre la tasa de asesinatos y la esperanza de vida, es decir, se ajustará el comportamiento de la tasa de homicidios, la cual se denotará con la variable \(y\), mediante la variable relacionada a la esperanza de vida en un determinado estado, la cual se denominará \(x\). Para ello se conformará un conjunto de vectores \((x_i, y_i)\) para cada pareja de observaciones, los cuales se grafican en el eje cartesiano, observando que la muestra sigue una tendencia inversamente proporcional, es decir, entre mayor sea la esperanza de vida en un estado entonces la tasa de homicidios tiende a disminuir de forma lineal.

\[\mbox{Murder} = \beta_0 + \beta_1 \mbox{Life_Exp} + \varepsilon_i\]

plot(datos$Life_Exp, datos$Murder, main='Relación vida vs homicidios')
library(plotly)

library(ggplot2)
grafica = ggplot(datos, aes(x=Life_Exp, y=Murder)) +
  geom_point(color="royalblue", size=5) +
  geom_smooth(method = 'lm', color="darkred")

ggplotly(grafica)

Calculando el modelo de regresión en R, en donde, \(Y= \mbox{Murder}\) y \(X=\mbox{Life_exp}\), así se tiene el siguiente código:

modelo = lm(Murder ~ Life_Exp, data = datos)
modelo
## 
## Call:
## lm(formula = Murder ~ Life_Exp, data = datos)
## 
## Coefficients:
## (Intercept)     Life_Exp  
##     159.576       -2.147

\[\mbox{Murder} = 159.576 -2.147 \mbox{Life_Exp} + \varepsilon_i\]

Analizando de nuevo, el modelo de regresión, se puede obtener varias relaciones y estadística por medio de la función summary

summary(modelo)
## 
## Call:
## lm(formula = Murder ~ Life_Exp, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7272 -1.6733 -0.1734  1.4909  4.8680 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  159.576     17.579   9.078 5.45e-12 ***
## Life_Exp      -2.147      0.248  -8.660 2.26e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.33 on 48 degrees of freedom
## Multiple R-squared:  0.6097, Adjusted R-squared:  0.6016 
## F-statistic: 74.99 on 1 and 48 DF,  p-value: 2.26e-11

La prueba de hipótesis de F-fisher es: \[ H_0 : \beta_0 = \beta_1 = \beta_2 = 0 \qquad \qquad vs \qquad \qquad H_a: \beta_0 \ne 0 \ o \ \beta_1 \ne 0 \ o \ \beta_2\ne 0 \] Consultando el anova

anova(modelo)
## Analysis of Variance Table
## 
## Response: Murder
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## Life_Exp   1 407.14  407.14  74.989 2.26e-11 ***
## Residuals 48 260.61    5.43                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Intervalos de confianza

confint(modelo, level = 0.95)
##                  2.5 %     97.5 %
## (Intercept) 124.231359 194.920026
## Life_Exp     -2.645874  -1.648729

La importancia de del modelo de regresión es que se puede dar estimaciones, por ejemplo, suponiendo que, en México hay una esperanza de vida de 70 años. ¿Cuál sería la tasa de homicidios? \[\hat{\mbox{Murder}} = 159.576 -2.147 * 70\] \[\hat{\mbox{Murder}} = 9.28\] La tasa de homicidios es de 9.28 asesinatos por cada 100,000 personas.

informacion = data.frame(Life_Exp = c(70) )
regresion = predict(modelo, newdata = informacion)
regresion
##        1 
## 9.264619

Gráficamente

plot(datos$Life_Exp, datos$Murder, main='Relación vida vs homocidios', pch=16, lwd=2)
points(x=70, y = as.numeric(regresion), col='red', pch=16, lwd=10)

Intervalos de regresión.

Los intervalos de regresión son intervalos de confianza sobre los valores sobre los valores de la regresión \[y = \hat{\beta_0} + \hat{\beta_1} x_i \]

regresion_intervalo = predict(modelo, newdata = informacion, interval = 'confidence', level = 0.95)
regresion_intervalo
##        fit      lwr      upr
## 1 9.264619 8.470351 10.05889

Intervalos de predicción

Los intervalos de predicción son intervalos de confianza sobre los valores de \(y = \beta_0 + \beta_1x_i + \varepsilon\) (son estimaciones de valores reales de y)

regresion_intervalo = predict(modelo, newdata = informacion, interval = 'prediction', level = 0.95)
regresion_intervalo
##        fit      lwr      upr
## 1 9.264619 4.512803 14.01643

Comparación entre intervalos

Algo interesante es comparar el intervalo de confianza y el de predicción, ya que se puede visualizar, que al estimar una observación nuestra variabilidad aumenta con ello las bandas son más grandes que estimar sólo un punto sobre la recta de regresión

extremos = data.frame(Life_Exp = c(60, 80))
regresion_limit = predict(modelo, newdata = extremos)

informacion_regresion = data.frame(Life_Exp = seq(60, 80, 0.5) )
regresion_intervalo = predict(modelo, newdata = informacion_regresion, interval = 'confidence', level=0.95)
predicion_intervalo = predict(modelo, newdata = informacion_regresion, interval = 'prediction', level=0.95)

plot(datos$Life_Exp, datos$Murder, main='Relación vida vs homocidios', pch=16, lwd=2)
points(x=70, y = as.numeric(regresion), col='red', pch=16, lwd=10)
lines(x=extremos$Life_Exp, y = regresion_limit, col='blue')
lines(informacion_regresion$Life_Exp, regresion_intervalo[,3], col='blue', lty = 3)
lines(informacion_regresion$Life_Exp, regresion_intervalo[,2], col='blue', lty = 3)
lines(informacion_regresion$Life_Exp, predicion_intervalo[,3], col='orange', lty = 2)
lines(informacion_regresion$Life_Exp, predicion_intervalo[,2], col='orange', lty = 2)