1 Introducción

Para analizar cómo realizar una regresión, vamos a usar una base de datos de pingüinos, que me parece bastante… intuitiva para este tipo de análisis.

#Cargamos la base de datos
library(palmerpenguins)
## Warning: package 'palmerpenguins' was built under R version 4.3.3
data("penguins")

#Cargo un par de librerías útiles
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.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
data <- penguins %>% 
  select(species, bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g, sex, island) %>%
  na.omit() #elimino valores faltantes.

2 ¿Qué queremos predecir o explicar?

La variable de interés va a ser el índice de masa corporal y lo vamos a querer explicar a partir de algunas otras variables físicas.

Antes de arrancar miremos un poco cómo está distribuida la variable de interés.

ggplot(data, aes(x = body_mass_g, fill = sex))+
  geom_histogram(stat = "count")+
  theme_light()+
  facet_wrap(~ island, scales = "free")
## Warning in geom_histogram(stat = "count"): Ignoring unknown parameters:
## `binwidth`, `bins`, and `pad`

Vemos que a simple vista el sexo parece una variable importante ya que se ve un patrón. Los individuos machos parecen tener mayor índice de masa corporal. La verdad es que en función de la isla hay cierta variación pero no demasiada. Parece más una tendencia general.

Veamos la relación entre el ala y el índice de masa corporal (dos columnas continuas).

ggplot(data, aes(x = body_mass_g, y = flipper_length_mm,
                 color = sex))+
  geom_point(position = position_jitter(width = 0.2,
                                        height = 0.2), alpha = 0.6,
             size = 2)+
  theme_light()

Bueno, las cosas se están poniendo extrañas. Así que vamos a ajustar un modelo simplón y veremos cosa succede.

3 Modelo Simple

ModeloSimple <- lm(body_mass_g ~ sex, data = data)
summary(ModeloSimple) 
## 
## Call:
## lm(formula = body_mass_g ~ sex, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1295.7  -595.7  -237.3   737.7  1754.3 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3862.27      56.83  67.963  < 2e-16 ***
## sexmale       683.41      80.01   8.542  4.9e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 730 on 331 degrees of freedom
## Multiple R-squared:  0.1806, Adjusted R-squared:  0.1781 
## F-statistic: 72.96 on 1 and 331 DF,  p-value: 4.897e-16
#No me gustan los residuos. Mirémoslos che.
residuos <- residuals(ModeloSimple)
plot(residuos)

Jajaja pero la PU’. ¿Qué es eso, hermano?

Le voy a agregar una dimensión más al análisis a ver si eso lo mejora un poco. Ya te digo que el ajuste es pésimo.

4 Modelo Múltiple

ModeloMultiple <- lm(body_mass_g ~ sex + flipper_length_mm,
                     data = data)
summary(ModeloMultiple)
## 
## Call:
## lm(formula = body_mass_g ~ sex + flipper_length_mm, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -910.28 -243.89   -2.94  238.85 1067.73 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -5410.300    285.798 -18.931  < 2e-16 ***
## sexmale             347.850     40.342   8.623 2.78e-16 ***
## flipper_length_mm    46.982      1.441  32.598  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 355.9 on 330 degrees of freedom
## Multiple R-squared:  0.8058, Adjusted R-squared:  0.8047 
## F-statistic: 684.8 on 2 and 330 DF,  p-value: < 2.2e-16
residuosMul <- residuals(ModeloMultiple)
plot(residuosMul)

5 Supuestos del Modelo

Bueno, amplia mejora. Pasamos de un r cuadrado de 18 a uno de 80. Una barbaridad. Esta variable es claramente un gran predictor.

Analicemos los supuestos con este modelo simplón.

#Normalidad
shapiro.test(residuosMul) #Residuos normales.
## 
##  Shapiro-Wilk normality test
## 
## data:  residuosMul
## W = 0.99793, p-value = 0.9569
#Homocedasticidad
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(ModeloMultiple) #auch, casi.
## 
##  studentized Breusch-Pagan test
## 
## data:  ModeloMultiple
## BP = 13.529, df = 2, p-value = 0.001154
#Multicolinealidad
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
vif(ModeloMultiple)
##               sex flipper_length_mm 
##          1.069646          1.069646