Repaso Regresión

Importar datos

x_enfermos <- seq(0, 8, 2)
y_hacinamiento <- c(9,11,11,14,15)
reg <- lm(y_hacinamiento ~ x_enfermos)
summary(reg)

Call:
lm(formula = y_hacinamiento ~ x_enfermos)

Residuals:
         1          2          3          4          5 
-1.069e-15  5.000e-01 -1.000e+00  5.000e-01 -6.080e-16 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9.0000     0.5477  16.432 0.000491 ***
x_enfermos    0.7500     0.1118   6.708 0.006760 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7071 on 3 degrees of freedom
Multiple R-squared:  0.9375,    Adjusted R-squared:  0.9167 
F-statistic:    45 on 1 and 3 DF,  p-value: 0.00676

\[ Y = \beta_1 + \beta_2 X + \epsilon \]

\[ Y = 9 + 0,75 \]

\[ \epsilon \sim N(0, \sigma^2) \]

Cómo verificamos el supuesto de normalidad?

1. QQPlot
2. Residuos vs predichos

Gráfico de QQplot para verificar normalidad:

plot(reg, 2)

Gráfico de residuos vs predichos para verificar homocedasticidad

plot(reg, 1)

Prueba F

Verificamos si el modelo es significativo, es decir, si la pendiente de nuestro modelo no es 0:

\[ H_0: \beta_2 = 0 \quad (\forall \beta_i = 0) \]

\[ H_1 : \beta_2 \neq 0 \quad (\exists \beta_i \neq 0) \]

Estadístico de prueba:

\[ \frac{CM_{regresión}}{CM_{error}} \sim F(gl_R, gl_E) \]

summary.aov(reg)

            Df Sum Sq Mean Sq F value  Pr(>F)   
x_enfermos   1   22.5    22.5      45 0.00676 **
Residuals    3    1.5     0.5                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nota

En este caso, existe un único \(\beta_i\), en particular \(\beta_2\).

Prueba T

summary(reg)

Call:
lm(formula = y_hacinamiento ~ x_enfermos)

Residuals:
         1          2          3          4          5 
-1.069e-15  5.000e-01 -1.000e+00  5.000e-01 -6.080e-16 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9.0000     0.5477  16.432 0.000491 ***
x_enfermos    0.7500     0.1118   6.708 0.006760 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7071 on 3 degrees of freedom
Multiple R-squared:  0.9375,    Adjusted R-squared:  0.9167 
F-statistic:    45 on 1 and 3 DF,  p-value: 0.00676

Tenemos una prueba t para el intercepto y otra para la pendiente:

Pendiente \(\beta_2\)

\[ H_0: \beta_2 = 0 \]

\[ H_1: \beta_2 \neq 0 \]

Estadístico de prueba:

\[ \frac{b_2-\beta_2}{S_{b_2}} \sim t_{n-2} \]

\[ \frac{0,75 - 0}{0,1118} = 6,35 \]

Prueba F vs Prueba T

Esta prueba es más específica que la prueba F, ya que esta última mide significancia del modelo, mientras que la prueba T mide la contribución de cada uno de los parámetros (incluso el intercepto). En el caso de regresión lineal simple, las pruebas T y F son equivalentes.

Nota

El P valor en esta prueba figura como Pr(>|t|)

Intervalos de confianza

Para hacer estimaciones puntuales, primero verificamos que el valor \(X\) pertenezca al intervalo de predicción

\[ Y_{X=5} = 9 + 0,75 \times 5 = 12,75 \]

No puedo conocer la recta poblacional, sino que utilizo la recta muestral. A partir de la recta muestral se puede crear intervalos de confianza

\[ \hat{Y}_i \sim N (\beta_1 + \beta_2 X, \sigma_{\hat{Y}}) \]

\[ IC_{(\hat{Y}/x=5)} = \hat{Y}_{x=5} \pm t_{n-2} \times \sqrt{\frac{s_e^2}{n} + s^2_e (\frac{(x_1-\bar{x})^2}{sc_x})} \]

\[ IC_{(\hat{Y}/x=5)} = \hat{Y}_{x=5} \pm t_{n-2} \times \underbrace{S_e \times \sqrt{\frac{1}{n} + \frac{(x_i - \bar{x})^2}{Sc_x}}}_{\text{factor de inflación}} \]

Reemplazando con los datos:

\[ IC_{(\hat{Y}/x=5)} = 12,75 \pm 3,18 \times 0,7071 \times \sqrt{\frac{1}{5} + \frac{(5-4)^2}{40}} \]

12.75 - 3.18 * 0.70711 * sqrt(1/5 + 1/40)

[1] 11.68339

12.75 + 3.18 * 0.70711 * sqrt(1/5 + 1/40)

[1] 13.81661

\[ IC_{(\hat{Y}/x=5)} = [11,68;13,81] \]

conf <- predict(reg, newdata = data.frame(x_enfermo=5), interval="confidence", level=0.95)

Warning: 'newdata' had 1 row but variables found have 5 rows

conf

   fit       lwr      upr
1  9.0  7.256902 10.74310
2 10.5  9.267444 11.73256
3 12.0 10.993622 13.00638
4 13.5 12.267444 14.73256
5 15.0 13.256902 16.74310

Interpretación del intervalo de confianza

Con un 95% de confianza puedo decir que el verdadero valor medio de hacinamiento para una cantidad de animales igual a 5 se encuentra entre 11,68 y 13,81.

Intervalo de predicción

Además de la variación de la recta, tenemos una variación propia de cada individuo (variación de la población de hacinamiento) \(S_e^2\):

\[ IP_{(Y/x=5)} = \hat{Y}_{x=5} \pm t_{n-2} \times \sqrt{ \underbrace{S_e^2}_{\text{variación propia}} + \underbrace{\frac{s_e^2}{n} + s^2_e (\frac{(x_1-\bar{x})^2}{sc_x})}_{\text{variación de la recta}}} \]

prediccion <- predict(reg, newdata = data.frame(x_enfermos = 5), interval = "predict", level = 0.95)
prediccion

    fit      lwr      upr
1 12.75 10.25934 15.24066

Interpretación del intervalo de predicción

Con una confianza del 95% puedo decir que el valor de hacinamiento para un individuo que tiene 5 animales, está entre 10,26 y 15,23.

Gráfico de los intervalos

valores <- seq(min(x_enfermos),max(x_enfermos), 1)
predictores <- 9 + 0.75*valores #Ecuacion de la recta
plot(x_enfermos, y_hacinamiento, pch=16, cex=1.25,
     main = "Hacinados = 9 + 0.75*Enfermos", sub = "Rojo = Int.Conf y Verde = Int.Predi") + lines(valores, predictores, col="blue", lwd=2)

integer(0)

# Calcular los valores de confianza Y prediccion
confianza <- predict(reg, newdata = data.frame(x_enfermos=valores), interval = "confidence", level=0.95)
prediccion <- predict(reg, newdata = data.frame(x_enfermos=valores), interval = "prediction", level=0.95)
# Agregar las bandas de confianza
lines(valores, confianza[, "lwr"], col="red", lty=2)  # limite inferior
lines(valores, confianza[, "upr"], col="red", lty=2)  # limite superior
# Agregar las bandas de prediccion
lines(valores, prediccion[, "lwr"], col="darkgreen", lty=2)  # limite inferior
lines(valores, prediccion[, "upr"], col="darkgreen", lty=2)  # limite superior