Ejemplo de análisis de regresión simple

Un analista de deportes quiere saber si existe una relación entre el número de bateos que realiza un equipo de béisbol y el número de runs que consigue. En caso de existir y de establecer un modelo, podría predecir el resultado del partido.

library(ggplot2)
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
datos<-read_excel("baseball_data.xlsx")

Representación gráfica de las observaciones

Representamos los datos para poder intuir si existe una relación y cuantificar dicha relación mediante un coeficiente de correlación. Si en este paso no se detecta la posible relación lineal, no tiene sentido seguir adelante generando un modelo lineal

ggplot(data = datos, mapping = aes(x = Numero_Bateos, y = Runs)) +
  geom_point(color = "firebrick", size = 2) +
  labs(title  =  'Diagrama de dispersión', x  =  'número  de bateos') +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

cor.test(x = datos$Numero_Bateos, y = datos$Runs, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  datos$Numero_Bateos and datos$Runs
## t = 4.0801, df = 28, p-value = 0.0003388
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3209675 0.7958231
## sample estimates:
##      cor 
## 0.610627

El gráfico y el test de correlación muestran una relación lineal, de intensidad considerable (r = 0.61) y significativa (p-value = 0.0003388). Tiene sentido intentar generar un modelo de regresión lineal que permita predecir el número de runs en función del número de bateos del equipo

2. Cálculo del modelo de regresión lineal simple

modelo_lineal <- lm(Runs ~ Numero_Bateos, datos)
# lm() devuelve el valor de la variable y para x=0 (intersección) 
# junto con la pendiente de la recta.
# Para ver la información del modelo se requiere summary().
summary(modelo_lineal)
## 
## Call:
## lm(formula = Runs ~ Numero_Bateos, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -125.58  -47.05  -16.59   54.40  176.87 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -2789.2429   853.6957  -3.267 0.002871 ** 
## Numero_Bateos     0.6305     0.1545   4.080 0.000339 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 66.47 on 28 degrees of freedom
## Multiple R-squared:  0.3729, Adjusted R-squared:  0.3505 
## F-statistic: 16.65 on 1 and 28 DF,  p-value: 0.0003388

La primera columna (Estimate) devuelve el valor estimado para los dos parámetros de la ecuación del modelo lineal (β0 y β1) que equivalen a la ordenada en el origen y la pendiente.

Se muestran los errores estándar, el valor del estadístico t y el p-value (dos colas) de cada uno de los dos parámetros. Esto permite determinar si los parámetros son significativamente distintos de 0, es decir, que tienen importancia en el modelo.

Para el modelo generado, tanto la ordenada en el origen como la pendiente son significativas (p-values < 0.05).

El R cuadrado es 0.3729, lo que significa que aproximadamente el 37.29% de la variabilidad en los Runs puede ser explicada por el número de bateos. Es una medida de la bondad del ajuste del modelo.

El R cuadrado ajustado es ligeramente más bajo (0.3505) y tiene en cuenta el número de predictores en el modelo (en este caso, solo uno). Es más preciso para comparar modelos con diferentes números de predictores.

El error estándar residual es 66.47, que proporciona una estimación de la desviación estándar de los términos de error (residuos).

El valor del Estadístico F es 16.65 con 1 y 28 grados de libertad, indicando que el modelo es significativamente diferente de un modelo sin predictores. El valor p asociado es muy bajo (0.000338), lo que confirma aún más la significancia estadística del modelo.

El modelo lineal generado sigue la ecuación Runs = -2789.2429 + 0.6305Bateos. Por cada unidad que se incrementa el número de bateos, el número de runs aumenta en promedio 0.6305 unidades.

3. Intervalos de confianza para los parámetros del modelo

confint(modelo_lineal)
##                       2.5 %        97.5 %
## (Intercept)   -4537.9592982 -1040.5264727
## Numero_Bateos     0.3139863     0.9471137

Tanto el intercepto como el número de bateos tienen intervalos de confianza que no incluyen el cero, lo que sugiere que ambos son significativos en el modelo

4.Representación gráfica del modelo

ggplot(data = datos, mapping = aes(x = Numero_Bateos, y = Runs)) +
  geom_point(color = "firebrick", size = 2) +
  labs(title  =  'Runs ~ número de bateos', x  =  'número  de bateos') +
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula = 'y ~ x'

Además de la línea de mínimos cuadrados es recomendable incluir los límites superior e inferior del intervalo de confianza. Esto permite identificar la región en la que, según el modelo generado y para un determinado nivel de confianza, se encuentra el valor promedio de la variable dependiente.

# Se genera una secuencia de valores x_i que abarquen todo el rango de las
# observaciones de la variable X
puntos <- seq(from = min(datos$Numero_Bateos),
              to = max(datos$Numero_Bateos),
              length.out = 100)
# Se predice el valor de la variable Y junto con su intervalo de confianza para
# cada uno de los puntos generados. En la función predict() hay que nombrar a 
# los nuevos puntos con el mismo nombre que la variable X del modelo.
# Devuelve una matriz.
limites_intervalo <- predict(object = modelo_lineal,
                             newdata = data.frame(Numero_Bateos = puntos),
                             interval = "confidence", level = 0.95)
head(limites_intervalo, 3) 
##        fit      lwr      upr
## 1 626.4464 584.5579 668.3350
## 2 628.3126 587.1743 669.4509
## 3 630.1788 589.7830 670.5745
# Finalmente se añaden al gráfico las líneas formadas por los límites
# superior e inferior.
plot(datos$Numero_Bateos, datos$Runs, col = "firebrick", pch = 19, ylab = "Runs",
     xlab = "Número de bateos", main = 'Runs ~ Número de bateos')
abline(modelo_lineal, col = 1)
lines(x = puntos, y = limites_intervalo[,2],type = "l", col = 2, lty = 3)
lines(x = puntos, y = limites_intervalo[,3],type = "l", col = 3, lty = 3)

La función geom_smooth() del paquete ggplot2 genera la regresión y su intervalo de forma directa.

# Por defecto incluye la región de 95% de confianza
ggplot(data = datos, mapping = aes(x = Numero_Bateos, y = Runs)) +
  geom_point(color = "firebrick", size = 2) +
  labs(title  =  'Diagrama de dispersión', x  =  'Número  de bateos') +
  geom_smooth(method = "lm", se = TRUE, color = "black") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula = 'y ~ x'

5.Verificar condiciones para poder aceptar un modelo lineal

Relación lineal entre variable dependiente e independiente:

Se calculan los residuos para cada observación y se representan (scatterplot). Si las observaciones siguen la línea del modelo, los residuos se deben distribuir aleatoriamente entorno al valor 0.

# La función lm() calcula y almacena los valores predichos por el modelo y los
# residuos.
datos$prediccion <- modelo_lineal$fitted.values
datos$residuos   <- modelo_lineal$residuals
head(datos)
## # A tibble: 6 × 5
##   Equipo  Numero_Bateos  Runs prediccion residuos
##   <chr>           <dbl> <dbl>      <dbl>    <dbl>
## 1 Texas            5659   855       779.     76.0
## 2 Boston           5710   875       811.     63.8
## 3 Detroit          5563   787       719.     68.5
## 4 Kansas           5672   730       787.    -57.2
## 5 St.              5532   762       699.     63.0
## 6 New_S.           5600   718       742.    -23.8
ggplot(data = datos, aes(x = prediccion, y = residuos)) +
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
  geom_hline(yintercept = 0) +
  geom_segment(aes(xend = prediccion, yend = 0), alpha = 0.2) +
  labs(title = "Distribución de los residuos", x = "predicción modelo",
       y = "residuo") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")

Los residuos se distribuyen de forma aleatoria entorno al 0 por lo que se acepta la linealidad.

Distribución normal de los residuos:

Los residuos se deben distribuir de forma normal con media 0. Para comprobarlo se recurre a histogramas, a los cuantiles normales o a un test de contraste de normalidad.

ggplot(data = datos, aes(x = residuos)) +
  geom_histogram(aes(y = ..density..)) +
  labs(title = "histograma de los residuos") +
  theme_light()
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

qqnorm(modelo_lineal$residuals)
qqline(modelo_lineal$residuals)

shapiro.test(modelo_lineal$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_lineal$residuals
## W = 0.96144, p-value = 0.337

Tanto la representación gráfica como el contraste de hipótesis confirman la distribución normal de los residuos.

Varianza constante de los residuos (Homocedasticidad):

La variabilidad de los residuos debe de ser constante a lo largo del eje X. Un patrón cónico es indicativo de falta de homogeneidad en la varianza.

ggplot(data = datos, aes(x = prediccion, y = residuos)) +
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
  geom_segment(aes(xend = prediccion, yend = 0), alpha = 0.2) +
  geom_smooth(se = FALSE, color = "firebrick") +
  labs(title = "Distribución de los residuos", x = "predicción modelo",
       y = "residuo") +
  geom_hline(yintercept = 0) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# Test de Breush-Pagan
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(modelo_lineal)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_lineal
## BP = 0.01269, df = 1, p-value = 0.9103

el valor de prueba de Breusch-Pagan es alto (p = 0.9103). Esto significa que no hay suficiente evidencia en los datos para rechazar la hipótesis nula de homocedasticidad. En otras palabras, no hay indicios estadísticamente significativos de que la varianza de los errores varíe con el nivel de las variables independientes en el modelo.

Ni la representación gráfica ni el contraste de hipótesis muestran evidencias que haga sospechar falta de homocedasticidad.

Autocorrelación de residuos:

Cuando se trabaja con intervalos de tiempo, es muy importante comprobar que no existe aoutocorrelación de los residuos, es decir que son independientes. Esto puede hacerse detectando visualmente patrones en la distribución de los residuos cuando se ordenan según se han registrado o con el test de Durbin-Watson dwt() del paquete Car.

ggplot(data = datos, aes(x = seq_along(residuos), y = residuos)) +
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
  geom_line(size = 0.3) +
  labs(title = "Distribución de los residuos", x = "index", y = "residuo") +
  geom_hline(yintercept = 0) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
## 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 representación de los residuos no muestra ninguna tendencia.