Librerías

Revisamos si tenemos todas la librerías necesarias con estos comandos en la consola:

library(performance)
library(skimr)
library(dplyr)

En caso de no tenerlos instalados, utilizamos estas líneas de código en la consola:

install.packages("performance", dependencies=TRUE)
install.packages("see", dependencies=TRUE)
install.packages("skimr", dependencies=TRUE)
install.packages("dplyr", dependencies=TRUE)

Caso 1: BMI

Lectura de datos

databmi = read.csv("https://egutierreza.netlify.app/uploads/databmi.csv")
names(databmi)
## [1] "bmi"       "bmi.papa"  "bmi.mama"  "edad.papa" "edad.mama" "edad"

En estos datos, observaremos que BMI es el Indice de masa corporal (IMC). En el dataset tenemos datos sobre el IMC y edad de una persona además de los IMC y edades de sus padres.

Comenzamos realizando una exploración de cada variable.

library(skimr)
skim(databmi)
Data summary
Name databmi
Number of rows 311
Number of columns 6
_______________________
Column type frequency:
numeric 6
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
bmi 0 1 25.31 4.27 15.6 22.10 25.0 27.75 43.5 ▃▇▅▁▁
bmi.papa 0 1 26.47 4.24 16.7 23.35 26.0 29.40 40.6 ▂▇▇▂▁
bmi.mama 0 1 29.56 5.07 16.3 26.05 29.7 32.65 45.6 ▁▆▇▂▁
edad.papa 0 1 59.62 10.60 38.7 51.60 59.0 66.85 92.1 ▅▇▇▃▁
edad.mama 0 1 54.45 10.30 21.5 47.55 53.9 61.00 85.9 ▁▃▇▃▁
edad 0 1 27.42 8.11 18.0 21.10 25.5 31.05 56.6 ▇▅▂▁▁

Exploración bivariada

Obtenemos las correlaciones por pares de variables:

cor(databmi)
##                  bmi    bmi.papa    bmi.mama   edad.papa  edad.mama        edad
## bmi       1.00000000  0.22829796  0.25430493  0.08375241  0.1533664  0.20557922
## bmi.papa  0.22829796  1.00000000  0.08127217 -0.35633830 -0.2374230 -0.23220622
## bmi.mama  0.25430493  0.08127217  1.00000000 -0.12954055 -0.1579986 -0.08774887
## edad.papa 0.08375241 -0.35633830 -0.12954055  1.00000000  0.7910062  0.68414693
## edad.mama 0.15336643 -0.23742304 -0.15799864  0.79100622  1.0000000  0.68763701
## edad      0.20557922 -0.23220622 -0.08774887  0.68414693  0.6876370  1.00000000

En específico, observamos la relación entre el IMC de una persona y su padre:

plot(databmi$bmi.papa,databmi$bmi,
     xlab="Indice de masa corporal (padre)",
     ylab="Indice de masa corporal (hijo)",
     main="Diagrama de dispersión")
abline(c(18,0.1),lty=2, col="blue")
abline(c(22,0.3),lty=3, col="green")
abline(c(20,0.12),lty=4, col="orange")

Y ahora, obtenemos los gráficos de correlación entre todos los pares de variables.

pairs(databmi)

Se puede observar que entre algunas variables hay mayor fuerza en la correlación.

  • ¿Entre qué variables se evidencia de alta correlación?

Evaluación de supuestos

Iniciamos revisando la normalidad de la variable de interés que se esta modelando como variable dependiente.

En este caso, la variable de interés es la variable ‘bmi’, el IMC de una persona.

Utilizamos una prubea de hipótesis para revisar la normalidad de estos datos.

  • $H_0: $ Los datos tienen distribución normal.
  • $H_1: $ Los datos NO tienen distribución normal.
shapiro.test(databmi$bmi)
## 
##  Shapiro-Wilk normality test
## 
## data:  databmi$bmi
## W = 0.97206, p-value = 9.802e-06

Observamos que el p-valor es menor a la significancia 0.05 (o 5%). Por lo tanto, se rechaza la hipótesis nula. NO hay normalidad en los datos.

En este caso, la variable dependiente debe transformarse. Utilizaremos la siguiente transformación:

databmi$bmi_2  = databmi$bmi^0.2

La cual cumple con el supuesto de normalidad:

shapiro.test(databmi$bmi_2)
## 
##  Shapiro-Wilk normality test
## 
## data:  databmi$bmi_2
## W = 0.9927, p-value = 0.1328

En los siguientes modelos utilizaremos la variable original para realizar la interpretación pero nótese que en la práctica se require utilizar la variable transformada.

Modelos

Modelo 1: BMI del padre

m1_papa  = lm(bmi ~ bmi.papa,databmi)
m1_papa
## 
## Call:
## lm(formula = bmi ~ bmi.papa, data = databmi)
## 
## Coefficients:
## (Intercept)     bmi.papa  
##     19.2272       0.2298

Estos resultados me dicen que el modelo de regresión es el siguiente:

\[ bmi = 19.23+0.23*bmi_{papa} \]

Interpretamos los coeficientes:

  • El intercepto (19.23) no tiene interpretación en este caso ya que no existe bmi=0.
  • Por cada unidad de aumento en bmi del padre, el hijo tendra 0.23 de bmi adicional.
summary(m1_papa)
## 
## Call:
## lm(formula = bmi ~ bmi.papa, data = databmi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.5626 -3.2692 -0.2422  2.5950 16.9424 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 19.22724    1.49446  12.866  < 2e-16 ***
## bmi.papa     0.22979    0.05575   4.122 4.83e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.165 on 309 degrees of freedom
## Multiple R-squared:  0.05212,    Adjusted R-squared:  0.04905 
## F-statistic: 16.99 on 1 and 309 DF,  p-value: 4.831e-05

Evaluar el performance del modelo

library(performance)
performance(m1_papa)
AIC BIC R2 R2_adjusted RMSE Sigma
1773.924 1785.143 0.05212 0.0490524 4.151141 4.164553

De estos resultados, se observa un R2 bajo.

Modelo 2: BMI de la madre

Analicemos ahora evaluamos un modelo con el bmi de la madre:

m2_mama  = lm(bmi ~ bmi.mama,databmi)
m2_mama
## 
## Call:
## lm(formula = bmi ~ bmi.mama, data = databmi)
## 
## Coefficients:
## (Intercept)     bmi.mama  
##     18.9736       0.2143

En este caso, los resultados me dicen que la relación entre el bmi de una persona y su madre es el siguiente:

\[ bmi = 18.97+0.21*bmi_{mama} \]

¿Cómo saber cual es mejor modelo? ¿Cuál debo utilizar?

Evaluo el performance de cada modelo y elijo el mejor.

compare_performance(m1_papa,m2_mama)
Name Model AIC BIC R2 R2_adjusted RMSE Sigma
m1_papa lm 1773.924 1785.143 0.052120 0.0490524 4.151141 4.164553
m2_mama lm 1769.779 1780.998 0.064671 0.0616440 4.123566 4.136889 
plot(compare_performance(m1_papa, m2_mama, rank = TRUE))

Modelo 3: Ambos BMI

Veamos ahora un nuevo modelo utilizando el bmi de la madre y padre juntos.

m3_ambos  = lm(bmi ~ bmi.mama+bmi.papa,databmi)
m3_ambos
## 
## Call:
## lm(formula = bmi ~ bmi.mama + bmi.papa, data = databmi)
## 
## Coefficients:
## (Intercept)     bmi.mama     bmi.papa  
##     13.8280       0.2000       0.2104
compare_performance(m1_papa,m2_mama, m3_ambos)
Name Model AIC BIC R2 R2_adjusted RMSE Sigma
m1_papa lm 1773.924 1785.143 0.0521200 0.0490524 4.151141 4.164553
m2_mama lm 1769.779 1780.998 0.0646710 0.0616440 4.123566 4.136889
m3_ambos lm 1757.003 1771.963 0.1080679 0.1022761 4.026768 4.046332

Gráficamente:

plot(compare_performance(m1_papa, m2_mama,m3_ambos, rank = TRUE))

¿Cuál modelo es mejor?

El modelo resulto ser el mejor de todos los probados.

Uso de la transformación

Ajuste de modelos

AL usar la transformación, observemos que la interpretación de los coeficientes resulta complicada. Sin embargo, aún se pueden evaluar estos modelos utilizando los indicadores de R cuadrado, AIC, y BIC.

m1_v2_papa  = lm(bmi_2 ~ bmi.papa,databmi)
m1_v2_papa
## 
## Call:
## lm(formula = bmi_2 ~ bmi.papa, data = databmi)
## 
## Coefficients:
## (Intercept)     bmi.papa  
##    1.815451     0.003351
m2_v2_mama  = lm(bmi_2 ~ bmi.mama,databmi)
m2_v2_mama
## 
## Call:
## lm(formula = bmi_2 ~ bmi.mama, data = databmi)
## 
## Coefficients:
## (Intercept)     bmi.mama  
##    1.808976     0.003219
m3_v2_ambos  = lm(bmi_2 ~ bmi.papa+bmi.mama,databmi)
m3_v2_ambos
## 
## Call:
## lm(formula = bmi_2 ~ bmi.papa + bmi.mama, data = databmi)
## 
## Coefficients:
## (Intercept)     bmi.papa     bmi.mama  
##    1.734171     0.003058     0.003011

Evaluación

compare_performance(m1_v2_papa,m2_v2_mama, m3_v2_ambos)
Name Model AIC BIC R2 R2_adjusted RMSE Sigma
m1_v2_papa lm -845.9057 -834.6863 0.0505542 0.0474816 0.0615096 0.0617084
m2_v2_mama lm -851.1917 -839.9724 0.0665556 0.0635347 0.0609891 0.0611861
m3_v2_ambos lm -863.4544 -848.4953 0.1083973 0.1026076 0.0596065 0.0598961 
plot(compare_performance(m1_v2_papa,m2_v2_mama, m3_v2_ambos, rank = TRUE))

Se obtiene que el modelo 3 con ambas variables es el mejor modelo.

Ejercicio:

Se ha creado un nuevo modelo adicionando la edad del papa, edad de la mama y edad de la persona (edad.papa, edad.mama y edad)

nuevo_modelo  = lm(bmi_2 ~ bmi.papa + bmi.mama + 
                     edad.papa + edad.mama + edad,databmi)
nuevo_modelo
## 
## Call:
## lm(formula = bmi_2 ~ bmi.papa + bmi.mama + edad.papa + edad.mama + 
##     edad, data = databmi)
## 
## Coefficients:
## (Intercept)     bmi.papa     bmi.mama    edad.papa    edad.mama         edad  
##   1.6191812    0.0040226    0.0034008   -0.0004548    0.0009621    0.0019210

Evaluar si este modelo es mejor que los anteriores. Usar la función performance(nuevo_modelo) para ver los resultados del nuevo modelo, y comparelos con los modelos anteriores. Modificar el siguiente código cambiando mnuevo por modelo_nuevo:

performance(mnuevo)
compare_performance(m1_v2_papa,
                    m2_v2_mama, 
                    m3_v2_ambos,
                    mnuevo)

Caso 2: Arboles

Cargar datos

En esta data, tenemos los siguiente datos de arboles:

  • Diámetro: pulgadas
  • Altura en pies.
  • Volumen en pies cúbicos.
library(dplyr)
data_arboles <- trees
names(data_arboles) <- c("diametro","altura","volumen")
head(data_arboles)
diametro altura volumen
8.3 70 10.3
8.6 65 10.3
8.8 63 10.2
10.5 72 16.4
10.7 81 18.8
10.8 83 19.7

Realizamos un análisis exploratorio univariado para cada variable.

library(skimr)
skim(data_arboles)
Data summary
Name data_arboles
Number of rows 31
Number of columns 3
_______________________
Column type frequency:
numeric 3
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
diametro 0 1 13.25 3.14 8.3 11.05 12.9 15.25 20.6 ▃▇▃▅▁
altura 0 1 76.00 6.37 63.0 72.00 76.0 80.00 87.0 ▃▃▆▇▃
volumen 0 1 30.17 16.44 10.2 19.40 24.2 37.30 77.0 ▇▅▁▂▁

Exploración bivariada

Ahora, exploramos las relaciones entre pares de variables.

cor(data_arboles)
##           diametro    altura   volumen
## diametro 1.0000000 0.5192801 0.9671194
## altura   0.5192801 1.0000000 0.5982497
## volumen  0.9671194 0.5982497 1.0000000
pairs(data_arboles)

¿Qué tipo de relaciones observamos entre las variables? ¿Entre qué variables encontramos alta correlación?

  • Observamos relaciones lineales positivas entre las variables.

Otra forma de ver un gráfico de dispersión:

library(ggplot2)
ggplot(data = data_arboles, mapping = aes(x = diametro, y = volumen)) +
geom_point(color = "firebrick", size = 2) +
labs(title = "Diagrama de dispersión", x = "Diámetro", y = "Volumen") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))

Correlación:

cor(data_arboles$diametro,data_arboles$volumen)
## [1] 0.9671194

¿Existe una alta correlación entre estas variables?

Modelo estadístico

Modelo 1: Con el diámetro

Obtenemos los coeficientes del modelo con el código:

m1_diametro <- lm(volumen~diametro,data_arboles)
m1_diametro
## 
## Call:
## lm(formula = volumen ~ diametro, data = data_arboles)
## 
## Coefficients:
## (Intercept)     diametro  
##     -36.943        5.066

Modelo final: \(Volumen = -36.943 + 0.0564*Diametro\) Interpretación de coeficientes:

  • El intercepto (-36.943) no tiene interpretación en este caso.
  • Por cada pulgada de aumento en el diámetro, el volumen aumenta en 5.066 pies cúbicos.
summary(m1_diametro)
## 
## Call:
## lm(formula = volumen ~ diametro, data = data_arboles)
## 
## 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 ***
## diametro      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

Modelo 2: Con la altura:

Ahora, observemos el modelo con la altura:

m2_altura <- lm(volumen~altura,data_arboles)
m2_altura
## 
## Call:
## lm(formula = volumen ~ altura, data = data_arboles)
## 
## Coefficients:
## (Intercept)       altura  
##     -87.124        1.543

Modelo final: \(Volumen = -87.124 + 1.543*Altura\) Interpretación de coeficientes:

  • El intercepto (-87.124) no tiene interpretación en este caso.
  • Por un aumento en 1 pie en la altura, el volumen aumenta en 1.543 pies cúbicos.
summary(m2_altura)
## 
## Call:
## lm(formula = volumen ~ altura, data = data_arboles)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.274  -9.894  -2.894  12.068  29.852 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -87.1236    29.2731  -2.976 0.005835 ** 
## altura        1.5433     0.3839   4.021 0.000378 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.4 on 29 degrees of freedom
## Multiple R-squared:  0.3579, Adjusted R-squared:  0.3358 
## F-statistic: 16.16 on 1 and 29 DF,  p-value: 0.0003784

¿Cuál de los modelos es mejor?

Veamos sus estadísticos y gráficos:

library(performance)
compare_performance(m1_diametro, m2_altura)
Name Model AIC BIC R2 R2_adjusted RMSE Sigma
m1_diametro lm 181.6447 185.9467 0.9353199 0.9330895 4.11254 4.251988
m2_altura lm 252.7986 257.1005 0.3579026 0.3357614 12.95762 13.396982 
plot(compare_performance(m1_diametro, m2_altura, rank = TRUE))

Ahora, observemos el modelo con la altura y el diametro:

m3_diam_altura <- lm(volumen~altura+diametro,data_arboles)
m3_diam_altura
## 
## Call:
## lm(formula = volumen ~ altura + diametro, data = data_arboles)
## 
## Coefficients:
## (Intercept)       altura     diametro  
##    -57.9877       0.3393       4.7082

Observamos sus resultados:

summary(m3_diam_altura)
## 
## Call:
## lm(formula = volumen ~ altura + diametro, data = data_arboles)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4065 -2.6493 -0.2876  2.2003  8.4847 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -57.9877     8.6382  -6.713 2.75e-07 ***
## altura        0.3393     0.1302   2.607   0.0145 *  
## diametro      4.7082     0.2643  17.816  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.882 on 28 degrees of freedom
## Multiple R-squared:  0.948,  Adjusted R-squared:  0.9442 
## F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16

Veamos sus estadísticos y gráficos:

compare_performance(m1_diametro, m2_altura, m3_diam_altura)
Name Model AIC BIC R2 R2_adjusted RMSE Sigma
m1_diametro lm 181.6447 185.9467 0.9353199 0.9330895 4.112540 4.251988
m2_altura lm 252.7986 257.1005 0.3579026 0.3357614 12.957617 13.396982
m3_diam_altura lm 176.9100 182.6459 0.9479500 0.9442322 3.689223 3.881832 
plot(compare_performance(m1_diametro, 
                         m2_altura, 
                         m3_diam_altura,
                         rank = TRUE))