Fundamentos de estadística y probabilidad para Data Science

4.2 Regresión lineal simple: Evaluando el modelo

El modelo de regresión lineal simple

Características del modelo lineal

Code

swedish_motor_insurance %>%
  ggplot(aes(x = n_claims, y = total_payment_sek)) +
  geom_point(aes(x = n_claims, y = total_payment_sek), size = 2, alpha = 0.5) +
  labs(
    title = "Precio del seguro vs partes",
    x = "Nº de partes",
    y = "Precio póliza"
  ) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", size = 1) +
  theme_bw()

¿Cuánto vale la ordenada en el origen?
¿Correlación positiva o negativa?
¿Correlación fuerte o débil?
¿Cómo es la pendiente de la recta, positiva, negativa o cercana a cero?

Ejercicio:

Dados los puntos \(p_1 = (40, 150)\) y \(p_2 = (110, 400)\), ¿Cómo calcularías la pendiente de la recta que los une?

\(m = \frac{y_2 - y_1}{x_2 - x_1} = \frac{400 - 150}{110 - 40} = \frac{250}{70}\)

cor(swedish_motor_insurance$n_claims, swedish_motor_insurance$total_payment_sek)

[1] 0.9128782

Generando el modelo de regresión lineal simple con la función `lm`

En su forma más básica, debemos especificar los parámetros: fórmula y data. Si miramos la documentación, vemos que el parámetro fórmula, debe ser un objeto de la clase fórmula. A modo simplificado, requiere una notación del tipo y = f(x) —> formula = y ~ x

¿Cuál es la variable objetivo, y la variable explicativa?

Si, a partir del número de reclamaciones queremos estimar el precio del seguro, nuestra variable objetivo será el precio y nuestra predictora el número de reclamaciones.

objeto_modelo <- lm(formula = total_payment_sek ~ n_claims, data = swedish_motor_insurance)
print(objeto_modelo)


Call:
lm(formula = total_payment_sek ~ n_claims, data = swedish_motor_insurance)

Coefficients:
(Intercept)     n_claims  
     19.994        3.414

Predicciones con `predict`:

Requiere el objeto modelo creado previamente y un dataframe que contenga un conjunto de datos (reales o sintéticos) de las variables predictoras.

str(swedish_motor_insurance)

'data.frame':   63 obs. of  2 variables:
 $ n_claims         : int  108 19 13 124 40 57 23 14 45 10 ...
 $ total_payment_sek: num  392.5 46.2 15.7 422.2 119.4 ...

Code

swedish_motor_insurance %>%
  mutate(prediction_payment_sek = predict(objeto_modelo, 
                                    newdata = select(., n_claims))) %>%
  ggplot() +
    geom_point(aes(x = n_claims, y = total_payment_sek), color = "blue") + 
    geom_point(aes(x = n_claims, y = prediction_payment_sek), color = "red", alpha = 0.5) +
    theme_bw()

Extrapolando datos fuera del rango observado:

Aunque la fórmula del modelo es la ecuación de una recta y su dominio e imagen son teóricamente infinitos, hay que ser cauteloso a la hora de extrapolar datos fuera de rango, usando un modelo creado con un rango determinado.

explanatory_data <- as.integer(round(seq(-10, 150, length.out = 63))) # Datos sintéticos
prediction_data <- predict(object = objeto_modelo, newdata = data.frame(n_claims = explanatory_data))
swedish_motor_insurance$prediction_data = prediction_data
swedish_motor_insurance$explanatory_data = explanatory_data
swedish_motor_insurance %>%
  ggplot() +
    geom_point(aes(x = n_claims, y = total_payment_sek), color = "blue") +
    geom_point(aes(x = explanatory_data, y = prediction_data), color = "red") +
    theme_bw()

Las métricas principales del modelo

mdl_lqsa_conv_vs_locura <- lm(formula = Convivencia ~ Locura, data = lqsa)
mdl_lqsa_conv_vs_locura


Call:
lm(formula = Convivencia ~ Locura, data = lqsa)

Coefficients:
(Intercept)       Locura  
    16.9515      -0.1096

typeof(mdl_lqsa_conv_vs_locura) ## Tipo interno en R

[1] "list"

class(mdl_lqsa_conv_vs_locura) ## Clase del objeto

[1] "lm"

Coeficiente de determinación R cuadrado (r-squared)

\(R^2 = 1 - \frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{\sum_{i=1}^{n}(y_i - \bar{y})^2}\)

\(y_i\) es el valor observado de la variable dependiente para la observación \(i\).
\(\hat{y}\) es el valor predicho por el modelo para la observación \(i\).
\(\bar{y}\) es la media de los valores observados \(y_i\)
\(n\) es el número total de observaciones

En el caso de la regresión lineal simple, coincide con el coeficiente de correlación al cuadrado.

Code

r_squared <- lqsa %>%
  mutate(Conv_pred = predict(object = mdl_lqsa_conv_vs_locura, 
                        newdata = select(., Locura))) %>%
  summarise(r_squared = 1-sum((Convivencia - Conv_pred)**2) /
            sum((Convivencia - mean(Convivencia))**2)) %>% pull()
r_squared

[1] 0.2879345

La regresión lineal calcula una ecuación que minimiza la distancia entre la línea del modelo y los puntos de los datos reales. Técnicamente, se minimiza la distancia al cuadrado entre el dato y el valor ajustado al modelo.

El parámetro estadístico r-cuadrado es una medida que nos indica el nivel de cercanía entre los datos y la línea de regresión ajustada. Se conoce también como coeficiente de determinación. Es frecuente encontrar en la bibliografía el término r-squared en minúscula cuando hablamos de regresión lineal simple y R-squared cuando tenemos más de una variable predictora.

El coeficiente nos dice qué porcentaje de la varianza explica el modelo.

Puede tomar valores entre 0 y 1
1 significa ajuste perfecto
0 significa que el modelo no explica la variabilidad de los datos de respuesta para los predictores la muestra.
En general, cuanto mayor es r-cuadrado, mejor se ajusta el modelo a los datos.

El coeficiente de determinación debemos relativizarlo al contexto de los datos. En ocasiones, 0.5, puede ser suficiente, mientras que en otras circunstancias es pura aleatoriedad.

Cálculo: correlación elevada al cuadrado. (Para regresión lineal simple)

Coeficiente de determinación ajustado R cuadrado ajustado

En regresión simple, si la muestra es suficientemente grande, es un valor muy parecido, aunque varía ligeramente, porque se usa n-1. El coeficiente tiene sentido cuando la fórmula utiliza más variables predictoras (k), ya que el coeficiente de determinación aumenta al incrementar k, aunque algunas de esas otras variables tengan poca influencia sobre la variable de respuesta. R cuadrado ajustado ya no es entre 0 y 1, sino que podría ser negativo si el nº de variables es grande respecto al nº de filas, especialmente, si el coeficiente de determinación es bajo.

\(1 - \frac{{(n-1)}}{{(n-k-1)}} \cdot (1-R^2)\)

Code

n <- nrow(lqsa)
k <- 1 # Grados de libertad (nº de var. ind.)
adj_r_squared <- 1 - ( (n-1)/(n-k-1) )*(1 - r_squared)
adj_r_squared

[1] 0.2625036

Error residual estándar (RSE / sigma): Grados de libertad

Otra métrica del modelo es el error residual estándard. Es la diferencia “estándar” entre el valor real y el ajustado por el modelo. Para calcularlo, se eleva al cuadrado el vector de residuos y se realiza la suma total y se calcula la raíz cuadrada y se divide por los grados de libertad (nº de observaciones - nº de coeficientes del modelo).

Lo que representa es por cuánto se equivoca el modelo.

\(RSE = \sqrt{\frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{n - p}}\)

Code

residuos <- residuals(mdl_lqsa_conv_vs_locura)
n <- length(residuos)
p <- length(coef(mdl_lqsa_conv_vs_locura))

# Calcula el RSE
RSE <- sqrt(sum(residuos^2) / (n - p))
RSE

[1] 5.307924

Error cuadrático medio (RMSE)

Se calcula con la misma fórmula que RSE pero dividiendo por el número de observaciones, no por los grados de libertad. Es más frecuente usar RSE.

Code

RMSE <- sqrt(sum(residuos^2) / (n))
RMSE

[1] 5.127942

Funciones para analizar modelos y métricas detalladas

Generales

Resumen: summary()

Muestra las métricas detalladas del modelo. El formato de salida no es tabular.

resumen_mdl <- summary(mdl_lqsa_conv_vs_locura)
resumen_mdl


Call:
lm(formula = Convivencia ~ Locura, data = lqsa)

Residuals:
   Min     1Q Median     3Q    Max 
-8.403 -3.645 -1.250  3.838 12.024 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 16.95153    2.48096   6.833 2.01e-07 ***
Locura      -0.10963    0.03258  -3.365  0.00224 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.308 on 28 degrees of freedom
Multiple R-squared:  0.2879,    Adjusted R-squared:  0.2625 
F-statistic: 11.32 on 1 and 28 DF,  p-value: 0.002236

broom::tidy()

Muestra información algo más resumida que summary y en formato tabular

broom::tidy(mdl_lqsa_conv_vs_locura)

# A tibble: 2 × 5
  term        estimate std.error statistic     p.value
  <chr>          <dbl>     <dbl>     <dbl>       <dbl>
1 (Intercept)   17.0      2.48        6.83 0.000000201
2 Locura        -0.110    0.0326     -3.36 0.00224

Glance:

glance(mdl_lqsa_conv_vs_locura)

# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.288         0.263  5.31      11.3 0.00224     1  -91.6  189.  193.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

sigma: estimación de la desviación estándar de los errores: n-k-1
statistic: F-statistic. Es un parámetro que se obtiene comparando dos modelos, uno con todos los coeficientes y otro sin incluir predictores. Si el modelo que incluye los predictores es significativamente explica mejor la variabilidad de la variable dependiente en función de la independiente, entonces se acepta el modelo, dado un nivel de significancia previamente establecido. \(f(x; d_1, d_2) = \frac{\sqrt{\frac{(d_1 x)^{d_1} \cdot d_2^{d_2}}{(d_1 x + d_2)^{d_1 + d_2}}}}{x \cdot \text{Beta}\left(\frac{d_1}{2}, \frac{d_2}{2}\right)}\) \(F = \frac{\text{SSR} / k}{\text{SSE} / (n - k - 1)}\)
p.value: A partir del parámetro anterior, y usando la distribución F, se obtiene el p-valor como la probabilidad de encontrar un valor tan extremo o más que F_statictic según la distribución F. Usamos la función de distrubución de probabilidad acumulada pf
df: Grados de libertad (nº de variables predictoras)
logLik: El logaritmo de la verosimilitud del modelo. Es una medida de cuán bien se ajusta el modelo a los datos, y se utiliza en la comparación de modelos. en la práctica, se usa junto con otros criterios, como AIC y BIC. Un valor más alto indica un mejor ajuste.
AIC (Criterio de Información de Akaike): Un criterio de selección de modelos que tiene en cuenta tanto la bondad de ajuste del modelo como la complejidad del modelo. Se utiliza para comparar modelos, siendo preferible un valor más bajo.
BIC (Criterio de Información Bayesiana): Similar al AIC, es otro criterio de selección de modelos que penaliza la complejidad del modelo. Al igual que el AIC, se utiliza para comparar modelos, y un valor más bajo es preferible.
deviance: Una medida de la falta de ajuste del modelo en comparación con un modelo nulo (modelo con solo el intercepto). Cuanto menor sea la deviance, mejor se ajusta el modelo a los datos.
df.residual: Los grados de libertad residuales, que representan el número de observaciones menos el número de coeficientes en el modelo.
nobs: El número total de observaciones en el modelo

Code

#sigma
sqrt(sum((mdl_lqsa_conv_vs_locura$residuals**2 - (mean(mdl_lqsa_conv_vs_locura$residuals))**2)) / 28)

[1] 5.307924

Code

# statistic (F)
n <- 30  # Número total de observaciones
k <- 1   # Número de variables predictoras

# Suma de Cuadrados Total (SST)
y_mean <- mean(mdl_lqsa_conv_vs_locura$model$Convivencia)

# Suma de Cuadrados de la Regresión (SSR)
y_pred <- predict(mdl_lqsa_conv_vs_locura)
SSR <- sum((y_pred - y_mean)**2)

# Suma de Cuadrados del Error (SSE)
SSE <- sum((mdl_lqsa_conv_vs_locura$model$Convivencia - y_pred)**2)

# Cálculo del Estadístico F
F_statistic <- (SSR / k) / (SSE / (n - k - 1))
F_statistic

[1] 11.32222

Code

# p.value
p_v <- 1- pf(F_statistic, k, n-k-1)
p_v

[1] 0.002236027

Code

# Parámetros
F_value <- 1.32
#F_value <- F_statistic
k <- 1
n <- 30
df1 <- k
df2 <- n - k - 1

# Crear un data frame con valores de F y sus densidades
F_values <- seq(0, 15, length.out = 1000)
densities <- df(F_values, df1, df2)
data_f <- data.frame(F_values, densities)

# Gráfico de la distribución F con ggplot2
ggplot(data_f, aes(x = F_values, y = densities)) +
    geom_line(color = "blue") +
    geom_area(data = subset(data_f, F_values >= F_value), fill = "lightblue", aes(y = densities)) +
    geom_vline(xintercept = F_value, linetype = "dashed", color = "red") +
    annotate("text", x = F_value + 1, y = max(densities), label = paste("F =", F_value), hjust = 1, vjust = 2) +
    labs(title = "Distribución F con df1 = 1 y df2 = 28",
         x = "Valor F", y = "Densidad") +
    theme_minimal()

Funciones y datos más específicos

Coeficientes: coefficients()

coefficients(mdl_lqsa_conv_vs_locura)

(Intercept)      Locura 
 16.9515284  -0.1096271

# mdl_lqsa_conv_vs_locura$coefficients
# 
# mdl_lqsa_conv_vs_locura %>%
#   coefficients()
# 
# mdl_lqsa_conv_vs_locura[["coefficients"]]

Ajuste - Fitted values: fitted()

fitted(mdl_lqsa_conv_vs_locura) %>% head() # Lo mismo que llamar a predict, con los datos de la variable explicativa

        1         2         3         4         5         6 
13.114579  5.988816  6.536951  5.988816  6.098443 14.978240

Residuos residuals():

Nos genera un vector con las diferencias entre los valores originales del dataset y los valores ajustados por el modelo

residuals(mdl_lqsa_conv_vs_locura) %>% head()

        1         2         3         4         5         6 
 3.885421 -3.988816 -2.536951 -2.988816 -1.098443  9.021760

Equivale a:

(lqsa$Convivencia - fitted(mdl_lqsa_conv_vs_locura)) %>% head()

        1         2         3         4         5         6 
 3.885421 -3.988816 -2.536951 -2.988816 -1.098443  9.021760

Varianza residual:

resumen_mdl$sigma**2 # Para la fórmula de la varianza, se usa n-p

[1] 28.17406

media_residuos <- mean(residuals(mdl_lqsa_conv_vs_locura)**2)
media_residuos

[1] 26.29579

media_residuos * 30 / 28

[1] 28.17406

augment: Obtiene métricas específicas para cada dato.

Code

augment(mdl_lqsa_conv_vs_locura) %>% head()

# A tibble: 6 × 8
  Convivencia Locura .fitted .resid   .hat .sigma .cooksd .std.resid
        <dbl>  <dbl>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl>      <dbl>
1          17     35   13.1    3.89 0.0797   5.35 0.0252       0.763
2           2    100    5.99  -3.99 0.0670   5.35 0.0217      -0.778
3           4     95    6.54  -2.54 0.0567   5.38 0.00728     -0.492
4           3    100    5.99  -2.99 0.0670   5.37 0.0122      -0.583
5           5     99    6.10  -1.10 0.0648   5.40 0.00159     -0.214
6          24     18   15.0    9.02 0.136    5.07 0.262        1.83

Convivencia: variable objetivo
Locura: variable explicativa
.fitted: Los valores ajustados (o predicciones) del modelo.
.resid: Los residuos del modelo. Un residuo es la diferencia entre el valor observado de la variable dependiente y el valor ajustado (predicho) por el modelo.
.hat: Los valores de apalancamiento son una medida de la influencia de cada observación en el ajuste del modelo. Cuanto mayor es el valor de .hat para una observación, mayor es su influencia. Estos valores son importantes en la identificación de valores atípicos.
.sigma: Una estimación de la desviación estándar de los residuos del modelo. En cada observación, .sigma excluye la observación en cuestión de la estimación. Es útil para entender la variabilidad de los residuos.
.cooksd: Las distancias de Cook son una medida de la influencia de cada observación en el conjunto de parámetros estimados del modelo. Una observación con una distancia de Cook alta puede ser un punto influyente, lo que significa que su eliminación del conjunto de datos podría cambiar significativamente la estimación de los coeficientes.
.std.resid: Los residuos estandarizados. Son los residuos divididos por su desviación estándar estimada. Son útiles para identificar valores atípicos y para verificar si los residuos se distribuyen uniformemente y cumplen con las suposiciones del modelo.

Visualizando el ajuste del modelo

Code

mdl_fz_vs_tam <- lm(fuerza ~ tamano, data = lotr)
mdl_fz_vs_tam_sin_an <- lm(fuerza ~ tamano, data = lotr_sin_anillo)

Code

glance(mdl_fz_vs_tam)

# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.292         0.267  21.7      11.6 0.00204     1  -134.  274.  278.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Code

glance(mdl_fz_vs_tam_sin_an)

# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic     p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>       <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.619         0.605  15.6      43.9 0.000000417     1  -120.  246.  250.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Code

augment(mdl_fz_vs_tam) %>% head(9)

# A tibble: 9 × 8
  fuerza tamano .fitted .resid   .hat .sigma  .cooksd .std.resid
   <dbl>  <dbl>   <dbl>  <dbl>  <dbl>  <dbl>    <dbl>      <dbl>
1     65      9    76.0 -11.0  0.0626   22.0 0.00905     -0.520 
2     76      7    64.3  11.7  0.0335   22.0 0.00522      0.549 
3     80     10    81.8  -1.80 0.0959   22.1 0.000400    -0.0869
4     49      9    76.0 -27.0  0.0626   21.5 0.0548      -1.28  
5     40      4    46.8  -6.75 0.0835   22.1 0.00479     -0.324 
6     49      7    64.3 -15.3  0.0335   21.9 0.00885     -0.715 
7     95      9    76.0  19.0  0.0626   21.8 0.0274       0.905 
8     56      8    70.1 -14.1  0.0418   22.0 0.00960     -0.663 
9    100      1    29.2  70.8  0.246    15.6 2.29         3.75

Code

lotr %>%
ggplot() +
  geom_point(aes(tamano, fuerza, size = magia), alpha = 0.33) +
  geom_smooth(aes(tamano, fuerza, color = "Con anillo"), method = "lm", se = FALSE) +
  geom_smooth(data = lotr_sin_anillo, aes(tamano, fuerza, color = "Sin anillo"), method = "lm", se=FALSE)+
  scale_color_manual(
    values = c("Sin anillo" = "red",
               "Con anillo" = "blue")
  ) +
  geom_text_repel(aes(x = tamano, y = fuerza, label = nombre), size = 3, vjust = -1) +
  theme_bw()

Code

library(ggfortify)

Warning: package 'ggfortify' was built under R version 4.3.2

Code

autoplot(mdl_fz_vs_tam,which = 1, ncol=1)

Code

autoplot(mdl_fz_vs_tam_sin_an,which = 1, ncol=1)

Code

autoplot(mdl_fz_vs_tam_sin_an,which = 2)

Code

autoplot(mdl_fz_vs_tam_sin_an)

Code

autoplot(mdl_fz_vs_tam)

¿Qué se calcula cuando se ajusta un modelo en R (Opcional):

Estimación de coeficientes: Coeficientes que minimizan la suma de los cuadrados de los residuos.
Cálculo del error Estándar de los Coeficientes: El error estándar de los coeficientes se calcula utilizando la matriz de varianza-covarianza de los coeficientes. Esta matriz se deriva de la varianza residual y de la matriz de diseño (X’X, donde X es la matriz de diseño que contiene las variables independientes). El error estándar de un coeficiente se calcula como la raíz cuadrada de la varianza de ese coeficiente. La fórmula general para el error estándar de un coeficiente (SE) es:

\(SE(\hat{\beta}_i) = \sqrt{Var(\hat{\beta}_i)}\)

Code

#Matriz de diseño: 
X_design <- model.matrix(mdl_lqsa_conv_vs_locura)
residual_variance <- summary(mdl_lqsa_conv_vs_locura)$sigma**2 

# Matriz de varianza-covarianza de los coeficientes
var_cov_matrix <- residual_variance * solve(t(X_design) %*% X_design)

# Mostrar la matriz de varianza-covarianza
var_cov_matrix

            (Intercept)       Locura
(Intercept)  6.15516858 -0.074408463
Locura      -0.07440846  0.001061462

Code

ee_Inter <- sqrt(var_cov_matrix[1,1])
ee_Inter

[1] 2.480961

Code

ee_Locur <- sqrt(var_cov_matrix[2,2])
ee_Locur

[1] 0.03258008

Cálculo del estadístico t para cada coeficiente:

\(t = \frac{\hat{\beta}{SE(\hat{\beta}_i)}_i}\)

Code

(mdl_lqsa_conv_vs_locura$coefficients["(Intercept)"]) / ee_Inter

(Intercept) 
   6.832645

Code

(mdl_lqsa_conv_vs_locura$coefficients["Locura"]) / ee_Locur

   Locura 
-3.364851

Cálculo del valor p para cada coeficiente

Se realiza usando una distribución t-student con n-p-1 grados de libertad

El modelo de regresión lineal simple

Características del modelo lineal

Generando el modelo de regresión lineal simple con la función lm

Predicciones con predict:

Extrapolando datos fuera del rango observado:

Las métricas principales del modelo

Coeficiente de determinación R cuadrado (r-squared)

Coeficiente de determinación ajustado R cuadrado ajustado

Error residual estándar (RSE / sigma): Grados de libertad

Error cuadrático medio (RMSE)

Funciones para analizar modelos y métricas detalladas

Generales

Funciones y datos más específicos

Visualizando el ajuste del modelo

¿Qué se calcula cuando se ajusta un modelo en R (Opcional):

Generando el modelo de regresión lineal simple con la función `lm`

Predicciones con `predict`: