Regresión lineal múltiple

Ejemplo Zarigüeyas

Author

Edimer David Jaramillo

Published

March 28, 2025

Bibliotecas

Code

library(tidyverse)
library(corrplot)
library(qqplotr)

library(yardstick) # métricas de desempeño

library(ggeffects)

Datos

Code

datos <- read_csv("../datos/zarigueyas.csv") |> 
  mutate(sexo = as.factor(sexo),
         sitio = as.factor(sitio),
         poblacion = as.factor(poblacion)) |> 
  select(-caso)
datos

1 variable categórica + 1 variable numérica

Modelo sin interacción

Code

modelo1 <-
  lm(edad ~ longitud_cabeza + sexo, data = datos)

modelo1 |> summary()


Call:
lm(formula = edad ~ longitud_cabeza + sexo, data = datos)

Residuals:
   Min     1Q Median     3Q    Max 
-3.589 -1.377 -0.268  1.206  4.850 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -13.1759     4.8046  -2.742 0.007252 ** 
longitud_cabeza   0.1857     0.0520   3.571 0.000553 ***
sexom            -0.3802     0.3694  -1.029 0.305955    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.817 on 98 degrees of freedom
Multiple R-squared:  0.1181,    Adjusted R-squared:  0.1001 
F-statistic: 6.562 on 2 and 98 DF,  p-value: 0.002116

Podemos obtener las predicciones de este modelo:

Code

pred_modelo1 <-
  predict_response(modelo1, terms = c("longitud_cabeza", "sexo"))

pred_modelo1

Ahora las graficamos:

Code

pred_modelo1 |>   
  plot()

Modelo con interacción

Code

modelo2 <-
  lm(edad ~ longitud_cabeza * sexo, data = datos)

modelo2 |> summary()


Call:
lm(formula = edad ~ longitud_cabeza * sexo, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5348 -1.3029 -0.3029  1.2508  4.7964 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)           -17.86944   10.35500  -1.726   0.0876 .
longitud_cabeza         0.23658    0.11222   2.108   0.0376 *
sexom                   5.62087   11.72196   0.480   0.6327  
longitud_cabeza:sexom  -0.06493    0.12677  -0.512   0.6097  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.824 on 97 degrees of freedom
Multiple R-squared:  0.1205,    Adjusted R-squared:  0.09328 
F-statistic: 4.429 on 3 and 97 DF,  p-value: 0.005808

Podemos obtener las predicciones de este modelo:

Code

pred_modelo2 <-
  predict_response(modelo2, terms = c("longitud_cabeza", "sexo"))

pred_modelo2

Ahora las graficamos:

Code

pred_modelo2 |>   
  plot()

Comparación de modelos

Obtener las predicciones de cada modelo:

Code

pred_mod1 <-
  predict(object = modelo1, newdata = datos)

pred_mod2 <-
  predict(object = modelo2, newdata = datos)

¿Cuál de los dos modelos es mejor?

Code

rmse_mod1 <- rmse_vec(truth = datos$edad, estimate = pred_mod1)
rmse_mod2 <- rmse_vec(truth = datos$edad, estimate = pred_mod2)

rmse_mod1

[1] 1.789608

Code

rmse_mod2

[1] 1.787193

Residuales del modelo

Primero creamos una tabla con los valores reales, predichos y residuales.

Code

tabla_residuales1 <- 
  data.frame(
    real = datos$edad,
    predicho = pred_mod2
  ) |> 
  mutate(residual = real - predicho)

tabla_residuales1

Validamos la normalidad:

Code

tabla_residuales1 |> 
  ggplot(aes(x = residual)) +
  geom_histogram(color = "black") +
  geom_vline(xintercept = 0, lty = 2, color = "red")

Podríamos hacer el cálculo de qué porcentaje de los datos está por encima o por debajo de 1.787193

Code

tabla_residuales1 |>
  mutate(resultado = if_else(abs(residual) > 1.787193, true = "si", false = "no")) |> 
  count(resultado) |> 
  mutate(porcentaje = n / sum(n))

Validamos la homocedasticidad:

Code

tabla_residuales1 |> 
  ggplot(aes(x = predicho, y = residual)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", lty = 2) +
  geom_hline(yintercept = -1.787193, color = "green", lty = 2) +
  geom_hline(yintercept = 1.787193, color = "green", lty = 2)

Modelo con 2 variables numéricas

Vamos a explorar la relación de la edad con la longitud de la cabeza y el tamaño del vientre.

Code

datos |> 
  ggplot(aes(x = longitud_cabeza, y = vientre, color = edad)) +
  geom_point() +
  scale_color_viridis_c()

Modelo sin interacción

Interpretación: para la longitud de la cabeza podemos decir: por cada unidad (centímetro) que aumenta la longitud de la cabeza esperamos que en promedio se incremente la edad en 0.10282, manteniendo constantes las demás variables.

Code

modelo3 <- lm(edad ~ longitud_cabeza + vientre, data = datos)
modelo3 |> summary()


Call:
lm(formula = edad ~ longitud_cabeza + vientre, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4963 -1.2630 -0.2295  1.1222  5.0440 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)     -11.61970    4.71707  -2.463   0.0155 *
longitud_cabeza   0.10282    0.06029   1.705   0.0913 .
vientre           0.18099    0.07777   2.327   0.0220 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.778 on 98 degrees of freedom
Multiple R-squared:  0.1553,    Adjusted R-squared:  0.138 
F-statistic: 9.006 on 2 and 98 DF,  p-value: 0.0002568

Modelo con interacción

¿Cómo puedo calcular la interacción en mis datos?

Code

datos |> 
  mutate(inter_lcab_vientre = longitud_cabeza * vientre) |> 
  ggplot(aes(x = inter_lcab_vientre, y = edad)) +
  geom_point() +
  geom_smooth(method = "lm", color = "dodgerblue", se = FALSE) +
  geom_smooth(method = "loess", color = "firebrick", se = FALSE)

Ajustamos el modelo con interacción:

Code

modelo4 <- lm(edad ~ longitud_cabeza * vientre, data = datos)
modelo4 |> summary()


Call:
lm(formula = edad ~ longitud_cabeza * vientre, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1177 -1.3225 -0.3367  0.9417  4.8105 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)             -163.73320   47.06153  -3.479 0.000755 ***
longitud_cabeza            1.74824    0.50999   3.428 0.000894 ***
vientre                    4.92291    1.46223   3.367 0.001092 ** 
longitud_cabeza:vientre   -0.05120    0.01577  -3.247 0.001601 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.697 on 97 degrees of freedom
Multiple R-squared:  0.2381,    Adjusted R-squared:  0.2145 
F-statistic:  10.1 on 3 and 97 DF,  p-value: 7.467e-06

Comparación de modelos

Obtener las predicciones de cada modelo:

Code

pred_mod3 <-
  predict(object = modelo3, newdata = datos)

pred_mod4 <-
  predict(object = modelo4, newdata = datos)

¿Cuál de los dos modelos es mejor?

Code

rmse_mod3 <- rmse_vec(truth = datos$edad, estimate = pred_mod3)
rmse_mod4 <- rmse_vec(truth = datos$edad, estimate = pred_mod4)

rmse_mod3

[1] 1.751507

Code

rmse_mod4

[1] 1.663432

Residuales del modelo

Primero creamos una tabla con los valores reales, predichos y residuales.

Code

tabla_residuales2 <- 
  data.frame(
    real = datos$edad,
    predicho = pred_mod4
  ) |> 
  mutate(residual = real - predicho)

tabla_residuales2

Validamos la normalidad:

Code

tabla_residuales2 |> 
  ggplot(aes(x = residual)) +
  geom_histogram(color = "black") +
  geom_vline(xintercept = 0, lty = 2, color = "red")

Podríamos hacer el cálculo de qué porcentaje de los datos está por encima o por debajo de 1.663432

Code

tabla_residuales2 |>
  mutate(resultado = if_else(abs(residual) > 1.663432, true = "si", false = "no")) |> 
  count(resultado) |> 
  mutate(porcentaje = n / sum(n))

Validamos la homocedasticidad:

Code

tabla_residuales2 |> 
  ggplot(aes(x = predicho, y = residual)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", lty = 2) +
  geom_hline(yintercept = -1.663432, color = "green", lty = 2) +
  geom_hline(yintercept = 1.663432, color = "green", lty = 2)

Modelo con muchas variables - Multicolinealidad

Ajuste del modelo

Code

modelo5 <- lm(edad ~ ., data = datos |> select(-poblacion))
modelo5 |> summary()


Call:
lm(formula = edad ~ ., data = select(datos, -poblacion))

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1817 -1.1100 -0.1268  0.8232  4.6205 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)     -16.54617    9.59997  -1.724   0.0885 .
sitio2            0.18269    0.79667   0.229   0.8192  
sitio3           -1.93996    1.49047  -1.302   0.1966  
sitio4           -3.73311    1.53219  -2.436   0.0169 *
sitio5           -1.18410    1.63469  -0.724   0.4709  
sitio6           -1.61635    1.72724  -0.936   0.3521  
sitio7           -1.71210    1.45234  -1.179   0.2418  
sexom            -0.15027    0.40835  -0.368   0.7138  
longitud_cabeza   0.13470    0.09796   1.375   0.1727  
anchura_craneo    0.09364    0.09223   1.015   0.3128  
longitud_total   -0.02198    0.09813  -0.224   0.8233  
cola              0.19131    0.17276   1.107   0.2713  
longitud_pata    -0.17645    0.09692  -1.820   0.0722 .
longitud_oreja   -0.01567    0.13406  -0.117   0.9073  
ojo               0.20058    0.21226   0.945   0.3474  
pecho             0.16045    0.14601   1.099   0.2749  
vientre           0.12142    0.09095   1.335   0.1855  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.762 on 84 degrees of freedom
Multiple R-squared:  0.2889,    Adjusted R-squared:  0.1534 
F-statistic: 2.132 on 16 and 84 DF,  p-value: 0.0138

¿Por qué no es posible estimar los coeficientes de la variable “Poblacion”?

Code

modelo_prueba <- lm(edad ~ poblacion + sitio, data = datos)

modelo_prueba |> summary()


Call:
lm(formula = edad ~ poblacion + sitio, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3077 -1.3077 -0.3077  0.8182  4.8182 

Coefficients: (1 not defined because of singularities)
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    3.4444     0.4536   7.593 2.26e-11 ***
poblacionVic   0.7374     0.5639   1.308    0.194    
sitio2        -0.7818     0.6947  -1.125    0.263    
sitio3         0.8413     0.8572   0.981    0.329    
sitio4        -0.1587     0.8572  -0.185    0.853    
sitio5         0.8632     0.7005   1.232    0.221    
sitio6        -0.1368     0.7005  -0.195    0.846    
sitio7             NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.925 on 94 degrees of freedom
Multiple R-squared:  0.05083,   Adjusted R-squared:  -0.009759 
F-statistic: 0.8389 on 6 and 94 DF,  p-value: 0.543

Code

datos |> 
  count(sitio, poblacion)

Diagnóstico de multicolinealidad

Matriz de correlación

Code

library(corrplot)

datos |> 
  select(where(is.numeric)) |> 
  cor() |> 
  corrplot(
    type = "lower",
    tl.srt = 45,
    method = "pie",
    diag = FALSE
  )

VIF

Code

library(car)
vif(modelo5)

                     GVIF Df GVIF^(1/(2*Df))
sitio           69.445661  6        1.423870
sexo             1.317471  1        1.147811
longitud_cabeza  3.825972  1        1.956009
anchura_craneo   2.636806  1        1.623824
longitud_total   5.462362  1        2.337170
cola             3.736692  1        1.933052
longitud_pata    5.892840  1        2.427517
longitud_oreja   9.541088  1        3.088865
ojo              1.626067  1        1.275173
pecho            2.803272  1        1.674298
vientre          1.982053  1        1.407854

Quitamos las variables con alto VIF:

Code

modelo6 <- lm(edad ~ ., data = datos |> select(-c(poblacion, sitio, longitud_oreja)))
modelo6 |> summary()


Call:
lm(formula = edad ~ ., data = select(datos, -c(poblacion, sitio, 
    longitud_oreja)))

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5282 -1.3212 -0.1776  0.9089  4.7681 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)     -10.94418    5.92918  -1.846   0.0682 .
sexom            -0.29534    0.40594  -0.728   0.4688  
longitud_cabeza   0.05930    0.09291   0.638   0.5249  
anchura_craneo    0.01378    0.08868   0.155   0.8769  
longitud_total    0.01797    0.08512   0.211   0.8333  
cola             -0.05298    0.13486  -0.393   0.6954  
longitud_pata    -0.03958    0.05720  -0.692   0.4907  
ojo               0.25224    0.19182   1.315   0.1918  
pecho             0.13398    0.14315   0.936   0.3518  
vientre           0.13280    0.09012   1.474   0.1440  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.81 on 91 degrees of freedom
Multiple R-squared:  0.1871,    Adjusted R-squared:  0.1067 
F-statistic: 2.328 on 9 and 91 DF,  p-value: 0.02084

Volvemos a calcular los VIF:

Code

vif(modelo6)

           sexo longitud_cabeza  anchura_craneo  longitud_total            cola 
       1.233962        3.262163        2.310531        3.894638        2.157994 
  longitud_pata             ojo           pecho         vientre 
       1.944891        1.258639        2.553996        1.844170

Calculemos los RMSE de los modelos 5 y 6:

Code

pred_mod5 <-
  predict(object = modelo5, newdata = datos)

pred_mod6 <-
  predict(object = modelo6, newdata = datos)

rmse_mod5 <- rmse_vec(truth = datos$edad, estimate = pred_mod5)
rmse_mod6 <- rmse_vec(truth = datos$edad, estimate = pred_mod6)

rmse_mod5

[1] 1.607042

Code

rmse_mod6

[1] 1.718152

Selección de predictoras

Método forward

Code

modelo_nulo <- lm(edad ~ 1, data = datos)

horizonte <- formula(
  edad ~ sexo + longitud_cabeza + anchura_craneo +
    longitud_total + cola + longitud_pata + ojo + pecho + vientre
)

modelo_forward <- step(object = modelo_nulo, scope = horizonte, direction = "forward")

Start:  AIC=132.26
edad ~ 1

                  Df Sum of Sq    RSS    AIC
+ vientre          1    47.752 319.04 120.17
+ pecho            1    41.171 325.62 122.23
+ longitud_cabeza  1    39.824 326.97 122.65
+ anchura_craneo   1    29.911 336.88 125.67
+ longitud_total   1    26.403 340.39 126.71
+ ojo              1    19.718 347.07 128.68
<none>                         366.79 132.26
+ longitud_pata    1     5.841 360.95 132.64
+ cola             1     5.300 361.49 132.79
+ sexo             1     1.226 365.57 133.92

Step:  AIC=120.17
edad ~ vientre

                  Df Sum of Sq    RSS    AIC
+ longitud_cabeza  1    9.1945 309.85 119.22
+ ojo              1    8.1068 310.93 119.57
+ pecho            1    7.7245 311.32 119.69
+ anchura_craneo   1    7.1724 311.87 119.87
<none>                         319.04 120.17
+ longitud_total   1    3.7632 315.28 120.97
+ sexo             1    0.2125 318.83 122.10
+ longitud_pata    1    0.0754 318.96 122.15
+ cola             1    0.0714 318.97 122.15

Step:  AIC=119.22
edad ~ vientre + longitud_cabeza

                 Df Sum of Sq    RSS    AIC
<none>                        309.85 119.22
+ ojo             1    4.4034 305.44 119.77
+ pecho           1    2.6361 307.21 120.35
+ sexo            1    1.2131 308.63 120.82
+ anchura_craneo  1    1.0499 308.80 120.87
+ longitud_pata   1    0.4833 309.36 121.06
+ longitud_total  1    0.1114 309.73 121.18
+ cola            1    0.0271 309.82 121.21

Code

modelo_forward


Call:
lm(formula = edad ~ vientre + longitud_cabeza, data = datos)

Coefficients:
    (Intercept)          vientre  longitud_cabeza  
       -11.6197           0.1810           0.1028

Método backward

Code

modelo_backward <- step(modelo6, direction = "backward")

Start:  AIC=129.33
edad ~ sexo + longitud_cabeza + anchura_craneo + longitud_total + 
    cola + longitud_pata + ojo + pecho + vientre

                  Df Sum of Sq    RSS    AIC
- anchura_craneo   1    0.0791 298.24 127.36
- longitud_total   1    0.1460 298.30 127.38
- cola             1    0.5056 298.66 127.50
- longitud_cabeza  1    1.3347 299.49 127.78
- longitud_pata    1    1.5693 299.73 127.86
- sexo             1    1.7343 299.89 127.92
- pecho            1    2.8700 301.03 128.30
- ojo              1    5.6655 303.82 129.23
<none>                         298.16 129.33
- vientre          1    7.1157 305.27 129.71

Step:  AIC=127.36
edad ~ sexo + longitud_cabeza + longitud_total + cola + longitud_pata + 
    ojo + pecho + vientre

                  Df Sum of Sq    RSS    AIC
- longitud_total   1    0.1425 298.38 125.41
- cola             1    0.4832 298.72 125.52
- longitud_pata    1    1.6040 299.84 125.90
- sexo             1    1.6863 299.92 125.93
- longitud_cabeza  1    1.9399 300.18 126.01
- pecho            1    3.6321 301.87 126.58
- ojo              1    5.9082 304.14 127.34
<none>                         298.24 127.36
- vientre          1    7.0533 305.29 127.72

Step:  AIC=125.41
edad ~ sexo + longitud_cabeza + cola + longitud_pata + ojo + 
    pecho + vientre

                  Df Sum of Sq    RSS    AIC
- cola             1    0.3504 298.73 123.53
- longitud_pata    1    1.5262 299.90 123.92
- sexo             1    2.0469 300.43 124.10
- longitud_cabeza  1    2.9276 301.31 124.39
- pecho            1    3.8553 302.23 124.70
<none>                         298.38 125.41
- ojo              1    6.0959 304.47 125.45
- vientre          1    7.0077 305.39 125.75

Step:  AIC=123.53
edad ~ sexo + longitud_cabeza + longitud_pata + ojo + pecho + 
    vientre

                  Df Sum of Sq    RSS    AIC
- longitud_pata    1    1.2453 299.97 121.95
- sexo             1    1.8924 300.62 122.16
- longitud_cabeza  1    2.6295 301.36 122.41
- pecho            1    3.9346 302.66 122.85
- ojo              1    5.9327 304.66 123.51
<none>                         298.73 123.53
- vientre          1    6.6647 305.39 123.75

Step:  AIC=121.95
edad ~ sexo + longitud_cabeza + ojo + pecho + vientre

                  Df Sum of Sq    RSS    AIC
- sexo             1    1.5567 301.53 120.47
- longitud_cabeza  1    1.9939 301.97 120.61
- pecho            1    3.0644 303.04 120.97
<none>                         299.97 121.95
- ojo              1    6.5768 306.55 122.14
- vientre          1    6.8038 306.78 122.21

Step:  AIC=120.47
edad ~ longitud_cabeza + ojo + pecho + vientre

                  Df Sum of Sq    RSS    AIC
- longitud_cabeza  1    1.3976 302.93 118.94
- pecho            1    3.9113 305.44 119.77
- ojo              1    5.6787 307.21 120.35
<none>                         301.53 120.47
- vientre          1    7.7453 309.28 121.03

Step:  AIC=118.94
edad ~ ojo + pecho + vientre

          Df Sum of Sq    RSS    AIC
<none>                 302.93 118.94
- pecho    1    8.0049 310.93 119.57
- ojo      1    8.3871 311.31 119.69
- vientre  1    9.6688 312.60 120.11

Code

modelo_backward


Call:
lm(formula = edad ~ ojo + pecho + vientre, data = select(datos, 
    -c(poblacion, sitio, longitud_oreja)))

Coefficients:
(Intercept)          ojo        pecho      vientre  
   -10.0007       0.2821       0.1767       0.1469

Método stepwise

Code

modelo_stepwise <- step(modelo6, direction = "both")

Start:  AIC=129.33
edad ~ sexo + longitud_cabeza + anchura_craneo + longitud_total + 
    cola + longitud_pata + ojo + pecho + vientre

                  Df Sum of Sq    RSS    AIC
- anchura_craneo   1    0.0791 298.24 127.36
- longitud_total   1    0.1460 298.30 127.38
- cola             1    0.5056 298.66 127.50
- longitud_cabeza  1    1.3347 299.49 127.78
- longitud_pata    1    1.5693 299.73 127.86
- sexo             1    1.7343 299.89 127.92
- pecho            1    2.8700 301.03 128.30
- ojo              1    5.6655 303.82 129.23
<none>                         298.16 129.33
- vientre          1    7.1157 305.27 129.71

Step:  AIC=127.36
edad ~ sexo + longitud_cabeza + longitud_total + cola + longitud_pata + 
    ojo + pecho + vientre

                  Df Sum of Sq    RSS    AIC
- longitud_total   1    0.1425 298.38 125.41
- cola             1    0.4832 298.72 125.52
- longitud_pata    1    1.6040 299.84 125.90
- sexo             1    1.6863 299.92 125.93
- longitud_cabeza  1    1.9399 300.18 126.01
- pecho            1    3.6321 301.87 126.58
- ojo              1    5.9082 304.14 127.34
<none>                         298.24 127.36
- vientre          1    7.0533 305.29 127.72
+ anchura_craneo   1    0.0791 298.16 129.33

Step:  AIC=125.41
edad ~ sexo + longitud_cabeza + cola + longitud_pata + ojo + 
    pecho + vientre

                  Df Sum of Sq    RSS    AIC
- cola             1    0.3504 298.73 123.53
- longitud_pata    1    1.5262 299.90 123.92
- sexo             1    2.0469 300.43 124.10
- longitud_cabeza  1    2.9276 301.31 124.39
- pecho            1    3.8553 302.23 124.70
<none>                         298.38 125.41
- ojo              1    6.0959 304.47 125.45
- vientre          1    7.0077 305.39 125.75
+ longitud_total   1    0.1425 298.24 127.36
+ anchura_craneo   1    0.0756 298.30 127.38

Step:  AIC=123.53
edad ~ sexo + longitud_cabeza + longitud_pata + ojo + pecho + 
    vientre

                  Df Sum of Sq    RSS    AIC
- longitud_pata    1    1.2453 299.97 121.95
- sexo             1    1.8924 300.62 122.16
- longitud_cabeza  1    2.6295 301.36 122.41
- pecho            1    3.9346 302.66 122.85
- ojo              1    5.9327 304.66 123.51
<none>                         298.73 123.53
- vientre          1    6.6647 305.39 123.75
+ cola             1    0.3504 298.38 125.41
+ anchura_craneo   1    0.0552 298.67 125.51
+ longitud_total   1    0.0098 298.72 125.52

Step:  AIC=121.95
edad ~ sexo + longitud_cabeza + ojo + pecho + vientre

                  Df Sum of Sq    RSS    AIC
- sexo             1    1.5567 301.53 120.47
- longitud_cabeza  1    1.9939 301.97 120.61
- pecho            1    3.0644 303.04 120.97
<none>                         299.97 121.95
- ojo              1    6.5768 306.55 122.14
- vientre          1    6.8038 306.78 122.21
+ longitud_pata    1    1.2453 298.73 123.53
+ longitud_total   1    0.1275 299.85 123.90
+ anchura_craneo   1    0.1048 299.87 123.91
+ cola             1    0.0696 299.90 123.92

Step:  AIC=120.47
edad ~ longitud_cabeza + ojo + pecho + vientre

                  Df Sum of Sq    RSS    AIC
- longitud_cabeza  1    1.3976 302.93 118.94
- pecho            1    3.9113 305.44 119.77
- ojo              1    5.6787 307.21 120.35
<none>                         301.53 120.47
- vientre          1    7.7453 309.28 121.03
+ sexo             1    1.5567 299.97 121.95
+ longitud_pata    1    0.9096 300.62 122.16
+ anchura_craneo   1    0.0441 301.49 122.45
+ cola             1    0.0312 301.50 122.46
+ longitud_total   1    0.0000 301.53 122.47

Step:  AIC=118.94
edad ~ ojo + pecho + vientre

                  Df Sum of Sq    RSS    AIC
<none>                         302.93 118.94
- pecho            1    8.0049 310.93 119.57
- ojo              1    8.3871 311.31 119.69
- vientre          1    9.6688 312.60 120.11
+ longitud_cabeza  1    1.3976 301.53 120.47
+ sexo             1    0.9604 301.97 120.61
+ anchura_craneo   1    0.5503 302.38 120.75
+ longitud_pata    1    0.4928 302.44 120.77
+ longitud_total   1    0.2648 302.66 120.85
+ cola             1    0.0000 302.93 120.94

Code

modelo_stepwise


Call:
lm(formula = edad ~ ojo + pecho + vientre, data = select(datos, 
    -c(poblacion, sitio, longitud_oreja)))

Coefficients:
(Intercept)          ojo        pecho      vientre  
   -10.0007       0.2821       0.1767       0.1469

Comparación de modelos

Calculemos los RMSE de los modelos forward, backward y stepwise:

Code

pred_mod_for <-
  predict(object = modelo_forward, newdata = datos)

pred_mod_back <-
  predict(object = modelo_backward, newdata = datos)

pred_mod_both <-
  predict(object = modelo_stepwise, newdata = datos)

rmse_mod_for <- rmse_vec(truth = datos$edad, estimate = pred_mod_for)
rmse_mod_back <- rmse_vec(truth = datos$edad, estimate = pred_mod_back)
rmse_mod_step <- rmse_vec(truth = datos$edad, estimate = pred_mod_both)

rmse_mod_for

[1] 1.751507

Code

rmse_mod_back

[1] 1.731846

Code

rmse_mod_step

[1] 1.731846

Configurando otros modelos

Code

modelo7 <- lm(edad ~ ojo + pecho + vientre + longitud_cabeza, data = datos)

modelo8 <- lm(edad ~ ojo*pecho + ojo*vientre + ojo*longitud_cabeza, data = datos)

modelo9 <- lm(edad ~ pecho*ojo + pecho*vientre + pecho*longitud_cabeza, data = datos)

modelo10 <- lm(edad ~ vientre*ojo + vientre*pecho + vientre*longitud_cabeza, data = datos)

modelo11 <- lm(edad ~ ojo*pecho*vientre + ojo*vientre*longitud_cabeza, data = datos)

modelo12 <- lm(edad ~ ojo*pecho*vientre*longitud_cabeza, data = datos)

pred_mod7 <-
  predict(object = modelo7, newdata = datos)

pred_mod8 <-
  predict(object = modelo8, newdata = datos)

pred_mod9 <-
  predict(object = modelo9, newdata = datos)

pred_mod10 <-
  predict(object = modelo10, newdata = datos)

pred_mod11 <-
  predict(object = modelo11, newdata = datos)

pred_mod12 <-
  predict(object = modelo12, newdata = datos)

rmse_mod7 <- rmse_vec(truth = datos$edad, estimate = pred_mod7)
rmse_mod8 <- rmse_vec(truth = datos$edad, estimate = pred_mod8)
rmse_mod9 <- rmse_vec(truth = datos$edad, estimate = pred_mod9)
rmse_mod10 <- rmse_vec(truth = datos$edad, estimate = pred_mod10)
rmse_mod11 <- rmse_vec(truth = datos$edad, estimate = pred_mod11)
rmse_mod12 <- rmse_vec(truth = datos$edad, estimate = pred_mod12)

rmse_mod7

[1] 1.727846

Code

rmse_mod8

[1] 1.697876

Code

rmse_mod9

[1] 1.609058

Code

rmse_mod10

[1] 1.635343

Code

rmse_mod11

[1] 1.614102

Code

rmse_mod12

[1] 1.556558

Agregando polinomios

Code

modelo13 <- lm(
  edad ~ poly(ojo, degree = 2) + pecho + vientre + longitud_cabeza,
  data = datos
)

modelo14 <- lm(
  edad ~ ojo + poly(pecho, degree = 2) + vientre + longitud_cabeza,
  data = datos
)

modelo15 <- lm(
  edad ~ ojo + pecho + poly(vientre, degree = 2) + longitud_cabeza,
  data = datos
)

modelo16 <- lm(
  edad ~ ojo + pecho + vientre + poly(longitud_cabeza, degree = 2),
  data = datos
)

modelo17 <- lm(
  edad ~ ojo + poly(pecho, degree = 2) + vientre + poly(longitud_cabeza, degree = 2),
  data = datos
)

modelo18 <- lm(
  edad ~ ojo + poly(pecho, degree = 2) + vientre + poly(longitud_cabeza, degree = 3),
  data = datos
)

pred_mod13 <-
  predict(object = modelo13, newdata = datos)

pred_mod14 <-
  predict(object = modelo14, newdata = datos)

pred_mod15 <-
  predict(object = modelo15, newdata = datos)

pred_mod16 <-
  predict(object = modelo16, newdata = datos)

pred_mod17 <-
  predict(object = modelo17, newdata = datos)

pred_mod18 <-
  predict(object = modelo18, newdata = datos)

rmse_mod13 <- rmse_vec(truth = datos$edad, estimate = pred_mod13)
rmse_mod14 <- rmse_vec(truth = datos$edad, estimate = pred_mod14)
rmse_mod15 <- rmse_vec(truth = datos$edad, estimate = pred_mod15)
rmse_mod16 <- rmse_vec(truth = datos$edad, estimate = pred_mod16)
rmse_mod17 <- rmse_vec(truth = datos$edad, estimate = pred_mod17)
rmse_mod18 <- rmse_vec(truth = datos$edad, estimate = pred_mod18)

rmse_mod13

[1] 1.727022

Code

rmse_mod14

[1] 1.664788

Code

rmse_mod15

[1] 1.716965

Code

rmse_mod16

[1] 1.6117

Code

rmse_mod17

[1] 1.59737

Code

rmse_mod18

[1] 1.578271

Tabla comparativa

Code

paste0("Modelo", 1:6)

[1] "Modelo1" "Modelo2" "Modelo3" "Modelo4" "Modelo5" "Modelo6"

Code

df_comparativa <-
  data.frame(
  modelo = c(
    paste0("Modelo", 1:6),
    "Forward",
    "Backward",
    "Stepwise",
    paste0("Modelo", 7:18)
  ),
  RMSE = c(
    rmse_mod1,
    rmse_mod2,
    rmse_mod3,
    rmse_mod4,
    rmse_mod5,
    rmse_mod6,
    rmse_mod_for,
    rmse_mod_back,
    rmse_mod_step,
    rmse_mod7,
    rmse_mod8,
    rmse_mod9,
    rmse_mod10,
    rmse_mod11,
    rmse_mod12,
    rmse_mod13,
    rmse_mod14,
    rmse_mod15,
    rmse_mod16,
    rmse_mod17,
    rmse_mod18
  )
)

df_comparativa |> 
  arrange(RMSE)

Diagnóstico de residuales en mejores modelos

Code

tabla_residuales_mod12 <- 
  data.frame(
    real = datos$edad,
    predicho = pred_mod12
  ) |> 
  mutate(residual = real - predicho)

tabla_residuales_mod12

Normalidad:

Code

tabla_residuales_mod12 |> 
  ggplot(aes(x = residual)) +
  geom_histogram(color = "black") +
  geom_vline(xintercept = 0, lty = 2, color = "red")

Homocedasticidad:

Code

tabla_residuales_mod12 |> 
  ggplot(aes(x = predicho, y = residual)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", lty = 2) +
  geom_hline(yintercept = -1.787193, color = "green", lty = 2) +
  geom_hline(yintercept = 1.787193, color = "green", lty = 2)

Code

tabla_residuales_mod18 <- 
  data.frame(
    real = datos$edad,
    predicho = pred_mod18
  ) |> 
  mutate(residual = real - predicho)

tabla_residuales_mod18

Normalidad:

Code

tabla_residuales_mod18 |> 
  ggplot(aes(x = residual)) +
  geom_histogram(color = "black") +
  geom_vline(xintercept = 0, lty = 2, color = "red")

Homocedasticidad:

Code

tabla_residuales_mod18 |> 
  ggplot(aes(x = predicho, y = residual)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", lty = 2) +
  geom_hline(yintercept = -1.787193, color = "green", lty = 2) +
  geom_hline(yintercept = 1.787193, color = "green", lty = 2)

Reales vs predichos

Code

tabla_residuales_mod12 |> 
  ggplot(aes(x = predicho, y = real)) +
  geom_point()

Code

tabla_residuales_mod18 |> 
  ggplot(aes(x = predicho, y = real)) +
  geom_point()

Resumen de modelos

Code

modelo12 |> summary()


Call:
lm(formula = edad ~ ojo * pecho * vientre * longitud_cabeza, 
    data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2162 -1.1196 -0.0589  0.8234  4.3252 

Coefficients:
                                    Estimate Std. Error t value Pr(>|t|)
(Intercept)                        4.483e+03  1.058e+04   0.424    0.673
ojo                               -2.610e+02  7.078e+02  -0.369    0.713
pecho                             -1.882e+02  4.085e+02  -0.461    0.646
vientre                           -1.429e+02  3.166e+02  -0.451    0.653
longitud_cabeza                   -4.626e+01  1.168e+02  -0.396    0.693
ojo:pecho                          1.094e+01  2.735e+01   0.400    0.690
ojo:vientre                        8.029e+00  2.116e+01   0.379    0.705
pecho:vientre                      5.854e+00  1.213e+01   0.482    0.631
ojo:longitud_cabeza                2.683e+00  7.807e+00   0.344    0.732
pecho:longitud_cabeza              1.941e+00  4.475e+00   0.434    0.666
vientre:longitud_cabeza            1.508e+00  3.469e+00   0.435    0.665
ojo:pecho:vientre                 -3.286e-01  8.112e-01  -0.405    0.686
ojo:pecho:longitud_cabeza         -1.124e-01  2.995e-01  -0.375    0.708
ojo:vientre:longitud_cabeza       -8.465e-02  2.318e-01  -0.365    0.716
pecho:vientre:longitud_cabeza     -6.155e-02  1.321e-01  -0.466    0.642
ojo:pecho:vientre:longitud_cabeza  3.452e-03  8.827e-03   0.391    0.697

Residual standard error: 1.697 on 85 degrees of freedom
Multiple R-squared:  0.3328,    Adjusted R-squared:  0.2151 
F-statistic: 2.827 on 15 and 85 DF,  p-value: 0.001308

Code

modelo18 |> summary()


Call:
lm(formula = edad ~ ojo + poly(pecho, degree = 2) + vientre + 
    poly(longitud_cabeza, degree = 3), data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3568 -1.1405 -0.0002  0.7443  4.6945 

Coefficients:
                                   Estimate Std. Error t value Pr(>|t|)   
(Intercept)                        -1.96436    3.50001  -0.561  0.57598   
ojo                                 0.17006    0.17592   0.967  0.33621   
poly(pecho, degree = 2)1            2.59430    2.33443   1.111  0.26929   
poly(pecho, degree = 2)2           -2.25295    1.86803  -1.206  0.23085   
vientre                             0.09886    0.08017   1.233  0.22063   
poly(longitud_cabeza, degree = 3)1  2.52755    2.30743   1.095  0.27617   
poly(longitud_cabeza, degree = 3)2 -5.35995    1.84018  -2.913  0.00449 **
poly(longitud_cabeza, degree = 3)3 -2.53235    1.68288  -1.505  0.13577   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.645 on 93 degrees of freedom
Multiple R-squared:  0.3141,    Adjusted R-squared:  0.2625 
F-statistic: 6.084 on 7 and 93 DF,  p-value: 7.351e-06

Estimaciones de modelos

Code

pred_modelo12 <-
  predict_response(modelo12, terms = c("ojo", "pecho"))

pred_modelo12 |> 
  plot()

Code

pred_modelo18 <-
  predict_response(modelo18, terms = c("longitud_cabeza", "ojo"))

pred_modelo18 |> 
  plot()