#1.- INTRODUCCIÓN

La regresión lineal es un método útil para predecir una respuesta cuantitativa Y partiendo de una sola variable predictora X, asumiendo que hay una relación aproximadamente lineal entre X e Y.

donde β0 (ordenada en el origen, valor esperado de Y cuando X=0 y β1 (pendiente, incremento medio de Y asociado con el aumento de X en una unidad) son las dos constantes o parámetros desconocidos en el modelo. Asumimos que el residuo o error ϵ (diferencia entre lo observado y estimado por el modelo) es independiente de X.

#2.- ESTIMACIÓN DE LOS COEFICIENTES

La verdadera recta de regresión poblacional suele ser desconocida, pero teniendo acceso a un conjunto de observaciones, podemos calcular un modelo aproximado, teniendo en cuenta que distintos conjuntos de datos pueden tender a generar rectas de regresión ligeramente distintas. Por tanto, en la práctica, β0 y β1 son desconocidos, por lo que para poder obtener una predicción de la variable respuesta, tenemos que obtener una estimación de los mismos utilizando los datos de entrenamiento.

El objetivo es obtener unos estimadores insesgados con los que el modelo lineal se ajuste bien a los datos disponibles. Para esto, la estrategia más comúnmente utilizada se basa en minimizar la suma de residuos al cuadrado (RSS), método conocido como mínimos cuadrado.

Los residuos no son más que la diferencia entre cada valor de la variable respuesta observada y la predicha por el modelo. Algunos residuos serán positivos (para observaciones por encima de la recta) y otros negativos (observaciones por debajo de la recta), siendo su promedio de 0. La recta que se ajuste bien a los datos tendrá residuos pequeños.

#3.- CONTRASTE DE HIPÓTESIS

El error estándar también puede usarse para llevar a cabo un test de hipótesis sobre los parámetros del modelo. El más común establece que

H0: β1 = 0 (no existe relación entre X e Y)

Ha: β1 ≠ 0 (existe alguna relación entre X e Y)

Para comprobar la hipótesis nula, necesitamos determinar si β1^ se aleja lo suficientemente de 0. La precisión con la que podemos determinar esto dependerá del SE(β1). Para ello llevamos a cabo un t-test, calculando el estadístico t, el cual mide el número de desviaciones estándar que el estimador β1^ y β1^ están del valor 0, y por último obtenemos el p-value.

Si el p-value es menor que el nivel de significancia establecido, podemos deducir que hay una relación entre el predictor y la variable respuesta.

#4.- BONDAD DE AJUSTE

En el caso de que la hipótesis nula sea rechazada, podemos cuantificar el grado con el que el modelo se ajusta a los datos. La calidad de un ajuste de regresión lineal es normalmente estimada usando dos valores relacionados: RSE y el estadístico o coeficiente de determinación R2, como fracción de la varianza explicada.

Cuanto más próximo sea a 1, mayor será la proporción de variabilidad en la variable respuesta explicada por el modelo. Determinar si su valor es lo suficientemente bueno dependerá de la aplicación en cuestión. Un valor bajo podría indicar que el modelo lineal no es adecuado o podría deberse a errores residuales, debido a variables no tenidas en cuenta.

A diferencia del RSE, R2 es independiente de la escala de medida de Y, es decir, adimensional, por lo que presenta la ventaja de ser más fácil de interpretar.

#5.- CONDICIONES REGRESIÓN LINEAL

  1. Linealidad: la relación entre el predictor y la variable respuesta ha de ser lineal. ¿Cómo comprobarlo? Diagrama de dispersión de los datos, graficar los residuos.

Distribución normal de los residuos: los residuos deben distribuirse de forma normal, en torno a 0. Esta condición podría no cumplirse en presencia de observaciones atípicas o outliers que no siguen el patrón del conjunto de datos. ¿Cómo comprobarlo? Histograma de los datos, distribución de quantiles (normal Q-Q plot), *test de hipótesis de normalidad

En la página [https://gallery.shinyapps.io/slr_diag/] se muestran ejemplos visuales intuitivos.

Variabilidad de los residuos constante (homocedasticidad): la variabilidad de los datos en torno a la recta de regresión ha de ser aproximadamente constante, lo cual implica que la variabilidad de los residuos debería ser constante también en torno a 0.

¿Cómo comprobarlo? Graficar los residuos, test de Breusch-Pagan

Independencia: las observaciones han de ser independientes unas de otras. Tener en cuenta en el caso de mediciones temporales. ¿Cómo comprobarlo? Graficar los residuos y estudiar si siguen un patrón o tendencia.

#6.- EJERCICIO

library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
head(trees)
##   Girth Height Volume
## 1   8.3     70   10.3
## 2   8.6     65   10.3
## 3   8.8     63   10.2
## 4  10.5     72   16.4
## 5  10.7     81   18.8
## 6  10.8     83   19.7
glimpse(trees)
## Rows: 31
## Columns: 3
## $ Girth  <dbl> 8.3, 8.6, 8.8, 10.5, 10.7, 10.8, 11.0, 11.0, 11.1, 11.2, 11.3, …
## $ Height <dbl> 70, 65, 63, 72, 81, 83, 66, 75, 80, 75, 79, 76, 76, 69, 75, 74,…
## $ Volume <dbl> 10.3, 10.3, 10.2, 16.4, 18.8, 19.7, 15.6, 18.2, 22.6, 19.9, 24.…
pairs(trees)

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
ggplot(data = trees, mapping = aes(x= Girth, y= Volume))+
  geom_point(color= "blue", size=2)+
  labs(title = "Scatter", x = "Girth", y = "Volume")+
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5))

cor.test(x= trees$Girth, y= trees$Volume, method = "pearson", digits=2)
## 
##  Pearson's product-moment correlation
## 
## data:  trees$Girth and trees$Volume
## t = 20.478, df = 29, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9322519 0.9841887
## sample estimates:
##       cor 
## 0.9671194

La estimación de la correlación esde un 0.9671 (directa y muy fuerte). El IC al 95% asociado al coeficiente de correlación poblaicona (rho) toma valores entre 0.93 y 0.98 (la correlación es mayor al 90%, puesto que los extremos son mayores a 0.9).

El p valor asociado es prácticamente cero, indicando que la correlación es estadísticamente significativa.

library(GGally)
## Warning: package 'GGally' was built under R version 4.2.3
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggpairs(trees, lower = list (continuous = "smooth"), diag = list(continuous = "bar"), axisLabels = "none")
## Warning in check_and_set_ggpairs_defaults("diag", diag, continuous =
## "densityDiag", : Changing diag$continuous from 'bar' to 'barDiag'
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

De lo hasta ahora analizado podemos concluir que:

I. Observando los gráficos de dispersión podemos observar que la variable Girth (diámetro) está más linealmente asociada con la variable respuesta Volume, por lo que utilizaremos ésta para el modelo.

El coeficiente de correlación de Pearson es bastante alto (r = 0,967) y significativo (p-value = 2,2x10−16). Ello indica una correlación entre ambas variables bastante intensa.

Tiene sentido generar el modelo de regresión lineal cumpliéndose estos primeros requisitos.

REGRESIÓN

modelo <- lm(Volume ~ Girth, data = trees)
summary(modelo)
## 
## Call:
## lm(formula = Volume ~ Girth, data = trees)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.065 -3.107  0.152  3.495  9.587 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -36.9435     3.3651  -10.98 7.62e-12 ***
## Girth         5.0659     0.2474   20.48  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.252 on 29 degrees of freedom
## Multiple R-squared:  0.9353, Adjusted R-squared:  0.9331 
## F-statistic: 419.4 on 1 and 29 DF,  p-value: < 2.2e-16

Constante es -36.94; se estima que el volumen medio independientemente del diámetro sería de -36.94 unidades. El coeficiente es significativo al 5% (P=0.000<0.05)

Girth (5.9659). Se estima que el aumento del diámetro en una unidad con el resto de factores constantes provocaría un incremento del volumen de 5.06 uds. El coeficiente es significativo al 5% (P = 0.000 < 0.05).

El coeficiente R2 93.53%. Indica que las variaciones del diámetro (Girth) permiten explicar el 93.53% de las variaciones del volumen alrededor de la media. (RF2 > 0.7 –> las predicciones son fiables)

El estadístico del contraste global es una F=419.4 con un p valor de 0.000 > 0.05. Tenemos evidencia empírica suficiente para rechazar H0 al 5%. las pendientes son globalmente significativas.

names(modelo)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"
residuals <- modelo$residuals
plot(residuals)

mean(residuals)
## [1] 3.044719e-16

La media de los errores es prácticamente cero (exogeneidad estricta y el estimador OLS es INSESGADO)

confint(modelo)
##                  2.5 %     97.5 %
## (Intercept) -43.825953 -30.060965
## Girth         4.559914   5.571799

El término independiente de la ecuación toma valores entre -30.06 y -43.82 para una confianza del 95%. Ambos extremos son negativos, indicando que el término independiente es negativo.

El efecto parcial que acompaña al diámetro toma valores entre 4.559 y 5.571 para una confianza del 95%. Ambos extremos son positivos, indicando que el incremento del diámetro provoca aumentos en el volumen.

REPRESENTACIÓN DE LA REGRESIÓN

ggplot(data= trees, mapping = aes(x = Girth, y = Volume))+
  geom_point(color="blue", size =2)+
  geom_smooth(method = "lm", se= TRUE, color = "red")+
  labs(title = "Volume = a + bxGirth + e", x = "Girth", y = "Volume")+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula = 'y ~ x'

CONDICIONES

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

shapiro.test(residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals
## W = 0.97889, p-value = 0.7811

El estadistico Shapiro-Wilks asciende a 0.97 con un p-valor de 0.78 > 0.05. Tenemos evidencia empírica suficiente para aceptar la hipótesis de normalidad.

library(lmtest)
## Warning: package 'lmtest' was built under R version 4.2.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.2.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 5.6197, df = 1, p-value = 0.01776

El p-valor es 0.01 < 0.05. Tenemos evidencia empírica suficiente para rechazar H0 homocedasticidad. El error tiene varianza heterocedástica.

ggplot(data = trees, aes(x = seq_along(modelo$residuals),
                         y = modelo$residuals)) + 
  geom_point(aes(color = modelo$residuals, size = 2)) + 
  scale_color_gradient2(low = "darkgreen", mid = "green", high = "lightgreen")+
  geom_line(size = 0.75)+
  labs(title = "Residuals", x = "obs", y = "residual") +
  geom_hline(yintercept = 0) + 
  theme(plot.title = element_text(hjust = 0.5), legend.position = "bottom")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

En este caso, la normalidad de los residuos podemos aceptarla, y tampoco parecen seguir una clara tendencia según el orden de registro de las observaciones, pero la condición de homocedasticidad parece no cumplirse. Al observar los gráficos podríamos sospechar que la observación 31 podría estar influyendo al modelo. Para analizar en qué medida pueda estar influyendo esta u otras observaciones, reajustaremos el modelo excluyendo posibles observaciones sospechosas. Dependiendo de la finalidad del modelo, la exclusión de posibles outliers debe analizarse con detalles, ya que estas observaciones podrían ser errores de medida, pero también podrían representar casos interesantes.

library(car)
## Warning: package 'car' was built under R version 4.2.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.2.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
summary(influence.measures(model = modelo))
## Potentially influential observations of
##   lm(formula = Volume ~ Girth, data = trees) :
## 
##    dfb.1_  dfb.Grth dffit   cov.r cook.d  hat    
## 31 -1.20_*  1.37_*   1.49_*  0.82  0.89_*  0.22_*
influencePlot(model = modelo)

##      StudRes        Hat      CookD
## 1   1.315573 0.11514037 0.10983618
## 20 -2.030440 0.03328798 0.06408063
## 31  2.837732 0.21519431 0.88805814

Tenemos un dato atípico influyente en la observación 31.

ggplot(data = trees, mapping = aes(x = Girth, y = Volume))+
  geom_point(color = "blue", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "red") + 
  geom_point(data = trees[31, ], color = "purple", size = 3, pch = 5)+
  geom_smooth(data = trees[-31, ], method = "lm", se = TRUE, color = "green") + 
  labs(x = "Girth", y= "Volume") +
  theme_bw()+theme(plot.title = element_text(hjust = 0.5), legend.position = "bottom")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

MODELO CORREGIDO

modelo2 <- lm(Volume ~Girth, data = trees[-31, ])
summary(modelo2)
## 
## Call:
## lm(formula = Volume ~ Girth, data = trees[-31, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5036 -2.3834 -0.0489  2.3951  6.3726 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -33.3104     3.2784  -10.16 6.76e-11 ***
## Girth         4.7619     0.2464   19.33  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.813 on 28 degrees of freedom
## Multiple R-squared:  0.9303, Adjusted R-squared:  0.9278 
## F-statistic: 373.6 on 1 and 28 DF,  p-value: < 2.2e-16

PREDICCIÓN

predict(modelo, data.frame(Girth = 15), interval = "confidence")
##        fit      lwr      upr
## 1 39.04439 37.24858 40.84019

La predicción en base al modelo 1 asciende a 39.04 unds. El intervalo de predicción al 95% de confianza indica que el valor predicho estaría entre 37.24 y 40.84 uds.

predict(modelo2, data.frame(Girth = 15), interval = "confidence")
##       fit      lwr      upr
## 1 38.1179 36.37165 39.86415

La predicción en base al modelo 1 asciende a 38.11 unds. El intervalo de predicción al 95% de confianza indica que el valor predicho estaría entre 36.37 y 39.86 uds.

Como es de esperar, el intervalo de predicción es mayor que el de confianza. Podemos interpretar la predicción de la siguiente manera: se espera que en promedio un árbol de diámetro 15 tenga un volumen de 39,04. Podemos afirmar con un 95% de confianza que el verdadero valor promedio se encuentra entre (37,24 – 40.84), mientras que el intervalo de predicción para un solo árbol de este diámetro será de (36.37 - 39.86).

#7.- INFO DE LA SESIÓN

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)
## # A tibble: 79 × 3
##    package    loadedversion source        
##    <chr>      <chr>         <chr>         
##  1 abind      1.4-5         CRAN (R 4.2.0)
##  2 bslib      0.4.2         CRAN (R 4.2.2)
##  3 cachem     1.0.6         CRAN (R 4.2.2)
##  4 callr      3.7.3         CRAN (R 4.2.3)
##  5 car        3.1-2         CRAN (R 4.2.3)
##  6 carData    3.0-5         CRAN (R 4.2.3)
##  7 cli        3.6.0         CRAN (R 4.2.2)
##  8 colorspace 2.1-0         CRAN (R 4.2.3)
##  9 crayon     1.5.2         CRAN (R 4.2.3)
## 10 devtools   2.4.5         CRAN (R 4.2.3)
## # ℹ 69 more rows