SML - Regresión

Regresión lineal simple

Cuando buscamos correlaciones entre variables es útil examinar todas las variables de form simultánea. La función pairs() sirve para visualizar variables de par en par en busca de relaciones entre ellas:

bw <- MASS::birthwt
pairs(~ bw$age + bw$bwt + bw$race)

El coeficiente de correlación de Pearson se puede obtener con cor():

# Una correlación simple:
#           cor(bw$age, bw$bwt)
#
# Una matriz de correlacion:
cor(bw[, c("age", "bwt", "race")])

             age         bwt       race
age   1.00000000  0.09031781 -0.1728180
bwt   0.09031781  1.00000000 -0.1947135
race -0.17281795 -0.19471349  1.0000000

El modelo lineal de regresión hace una estimación por el método de estimación por mínimos cuadrados:

modelo <- lm(age ~ bwt, data = bw)
summary(modelo)

Call:
lm(formula = age ~ bwt, data = bw)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.8920 -3.9614 -0.7409  3.1685 20.4196 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.131e+01  1.605e+00   13.27   <2e-16 ***
bwt         6.563e-04  5.292e-04    1.24    0.216    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.291 on 187 degrees of freedom
Multiple R-squared:  0.008157,  Adjusted R-squared:  0.002853 
F-statistic: 1.538 on 1 and 187 DF,  p-value: 0.2165

lm() devuelve un objeto devuelve toda la información sobre el análisis, con un coeficiente de correlación que habla sobre la bondad del ajuste (\(R^{2}\)) - RMSE (root mean squared error). Es una medida de cuanto en qué medida impacta la variable independiente sobre la variabilidad de los datos observada.

Cómo reportar R²

Cuando \(R^{2}\) se obtiene entre dos variables, \(\sqrt{R^{2}} = r\) (\(r\) siendo el coeficiente de correlación de Pearson).

Esto devuelve una información para una recta, donde \(y\) es la variable dependiente, \(x\) es la variable independiente.

\[ y = \alpha x + \beta \]

\(\alpha\) es la pendiente y \(\beta\) es el punto en el que la línea se cruza con el eje Y (intercepción).

plot(bw$bwt, bw$age, xlab = "bwt", ylab = "age")
abline(modelo)

Regresión lineal múltiple

Se pueden explorar todas las variables de un dataframe al mismo tiempo

cor(x = bw, method = "pearson")

              low         age         lwt         race       smoke          ptl
low    1.00000000 -0.11893928 -0.16962694  0.137792751  0.16140431  0.196087267
age   -0.11893928  1.00000000  0.18007315 -0.172817953 -0.04434618  0.071606386
lwt   -0.16962694  0.18007315  1.00000000 -0.165048544 -0.04417908 -0.140029003
race   0.13779275 -0.17281795 -0.16504854  1.000000000 -0.33903074  0.007951293
smoke  0.16140431 -0.04434618 -0.04417908 -0.339030745  1.00000000  0.187557063
ptl    0.19608727  0.07160639 -0.14002900  0.007951293  0.18755706  1.000000000
ht     0.15237025 -0.01583700  0.23636040  0.019929917  0.01340704 -0.015399579
ui     0.16904283 -0.07515558 -0.15276317  0.053602088  0.06215900  0.227585340
ftv   -0.06296026  0.21539394  0.14052746 -0.098336254 -0.02801314 -0.044429660
bwt   -0.78480516  0.09031781  0.18573328 -0.194713487 -0.19044806 -0.154653390
               ht          ui         ftv         bwt
low    0.15237025  0.16904283 -0.06296026 -0.78480516
age   -0.01583700 -0.07515558  0.21539394  0.09031781
lwt    0.23636040 -0.15276317  0.14052746  0.18573328
race   0.01992992  0.05360209 -0.09833625 -0.19471349
smoke  0.01340704  0.06215900 -0.02801314 -0.19044806
ptl   -0.01539958  0.22758534 -0.04442966 -0.15465339
ht     1.00000000 -0.10858506 -0.07237255 -0.14598189
ui    -0.10858506  1.00000000 -0.05952341 -0.28392741
ftv   -0.07237255 -0.05952341  1.00000000  0.05831777
bwt   -0.14598189 -0.28392741  0.05831777  1.00000000

multimodel <- lm(bwt ~ low + 
                 age + 
                 lwt + 
                 race + 
                 smoke + 
                 ptl + 
                 ht + 
                 ui + 
                 ftv, data = bw)
summary(multimodel)

Call:
lm(formula = bwt ~ low + age + lwt + race + smoke + ptl + ht + 
    ui + ftv, data = bw)

Residuals:
    Min      1Q  Median      3Q     Max 
-991.22 -300.96   -5.39  277.74 1637.80 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3612.508    229.457  15.744  < 2e-16 ***
low         -1131.217     73.957 -15.296  < 2e-16 ***
age            -6.245      6.347  -0.984 0.326416    
lwt             1.051      1.133   0.927 0.355085    
race         -100.905     38.544  -2.618 0.009605 ** 
smoke        -174.116     72.000  -2.418 0.016597 *  
ptl            81.340     68.552   1.187 0.236980    
ht           -181.955    137.661  -1.322 0.187934    
ui           -336.776     93.314  -3.609 0.000399 ***
ftv            -7.578     30.992  -0.245 0.807118    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 433.7 on 179 degrees of freedom
Multiple R-squared:  0.6632,    Adjusted R-squared:  0.6462 
F-statistic: 39.16 on 9 and 179 DF,  p-value: < 2.2e-16

Hecho esto, se pueden seleccionar las mejores correlaciones utilizando el criterio de información de Aikaike o AIC utilizando step

step(object = multimodel, 
     directio = "both",
     trace = 1)

Start:  AIC=2305.09
bwt ~ low + age + lwt + race + smoke + ptl + ht + ui + ftv

        Df Sum of Sq      RSS    AIC
- ftv    1     11246 33682899 2303.2
- lwt    1    161709 33833361 2304.0
- age    1    182159 33853812 2304.1
- ptl    1    264837 33936490 2304.6
- ht     1    328638 34000291 2304.9
<none>               33671653 2305.1
- smoke  1   1100074 34771727 2309.2
- race   1   1289179 34960832 2310.2
- ui     1   2450185 36121838 2316.4
- low    1  44009294 77680946 2461.1

Step:  AIC=2303.15
bwt ~ low + age + lwt + race + smoke + ptl + ht + ui

        Df Sum of Sq      RSS    AIC
- lwt    1    154400 33837299 2302.0
- age    1    205543 33888442 2302.3
- ptl    1    268799 33951698 2302.7
- ht     1    319845 34002744 2302.9
<none>               33682899 2303.2
+ ftv    1     11246 33671653 2305.1
- smoke  1   1095243 34778142 2307.2
- race   1   1280252 34963151 2308.2
- ui     1   2442561 36125460 2314.4
- low    1  44045331 77728230 2459.2

Step:  AIC=2302.02
bwt ~ low + age + race + smoke + ptl + ht + ui

        Df Sum of Sq      RSS    AIC
- age    1    158960 33996259 2300.9
- ht     1    229750 34067049 2301.3
- ptl    1    232303 34069602 2301.3
<none>               33837299 2302.0
+ lwt    1    154400 33682899 2303.2
+ ftv    1      3937 33833361 2304.0
- smoke  1   1132887 34970186 2306.2
- race   1   1427486 35264785 2307.8
- ui     1   2527757 36365056 2313.6
- low    1  45611646 79448944 2461.3

Step:  AIC=2300.9
bwt ~ low + race + smoke + ptl + ht + ui

        Df Sum of Sq      RSS    AIC
- ptl    1    189792 34186051 2299.9
- ht     1    228505 34224763 2300.2
<none>               33996259 2300.9
+ age    1    158960 33837299 2302.0
+ lwt    1    107817 33888442 2302.3
+ ftv    1     19622 33976637 2302.8
- smoke  1   1054854 35051113 2304.7
- race   1   1301334 35297593 2306.0
- ui     1   2450537 36446796 2312.1
- low    1  45476560 79472818 2459.4

Step:  AIC=2299.95
bwt ~ low + race + smoke + ht + ui

        Df Sum of Sq      RSS    AIC
- ht     1    236483 34422534 2299.3
<none>               34186051 2299.9
+ ptl    1    189792 33996259 2300.9
+ age    1    116450 34069602 2301.3
+ lwt    1     85589 34100462 2301.5
+ ftv    1     21990 34164061 2301.8
- smoke  1    941319 35127371 2303.1
- race   1   1268492 35454543 2304.8
- ui     1   2279585 36465636 2310.2
- low    1  45501243 79687294 2457.9

Step:  AIC=2299.26
bwt ~ low + race + smoke + ui

        Df Sum of Sq      RSS    AIC
<none>               34422534 2299.3
+ ht     1    236483 34186051 2299.9
+ ptl    1    197771 34224763 2300.2
+ age    1    114599 34307935 2300.6
+ lwt    1     24448 34398086 2301.1
+ ftv    1     12780 34409754 2301.2
- smoke  1    935878 35358412 2302.3
- race   1   1269413 35691947 2304.1
- ui     1   2122280 36544815 2308.6
- low    1  48004723 82427257 2462.3

Call:
lm(formula = bwt ~ low + race + smoke + ui, data = bw)

Coefficients:
(Intercept)          low         race        smoke           ui  
    3586.50     -1139.20       -97.34      -157.42      -303.19  

El modelo óptimo es el que tiene el AIC más bajo.

Akaike information criterion

Este algoritmo intenta minimizar la pérdida de información de un modelo (bondad del ajuste o goodness of fit), intentando al mismo tiempo que sea lo más simple posible.

Utiliza la función de verosimilitud para la bondad de ajuste, penalizando el número de variables.

Wikipedia

Podemos tomar la información del análisis previo y generar un nuevo modelo

multimodel2 <- lm(bwt ~ low +
                  race + 
                  smoke + 
                  ui, data = bw)
summary(multimodel2)

Call:
lm(formula = bwt ~ low + race + smoke + ui, data = bw)

Residuals:
    Min      1Q  Median      3Q     Max 
-1025.8  -351.0    30.8   285.8  1500.8 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3586.50      86.68  41.379  < 2e-16 ***
low         -1139.20      71.12 -16.019  < 2e-16 ***
race          -97.34      37.37  -2.605 0.009942 ** 
smoke        -157.42      70.38  -2.237 0.026510 *  
ui           -303.19      90.02  -3.368 0.000922 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 432.5 on 184 degrees of freedom
Multiple R-squared:  0.6557,    Adjusted R-squared:  0.6482 
F-statistic: 87.59 on 4 and 184 DF,  p-value: < 2.2e-16

Evaluación del modelo

La forma de evaluar la bondad del modelo generado es utilizando la raíz del error cuadrático medio (RMSE; root mean squared edrror) en regresión; y precisión predictiva del modelo para la clasificación.

Para calcular estas medidas no se deben utilizar muestras que contengan elementos incluidos en la muestra de entrenamiento porque sobreestiman la bondad del modelo.

Rendimiento en regresión

El paquete caret (Classification And REgression Training) contiene muchas herramientas útiles en SML, como caret::predict.

in-sample y out-of-sample error

Cuando el modelo se valida contra sí mismo decimos que la raíz del error cuadrático medio es in-sample error - pertenece a la misma muestra con la que se ha entrando el modelo y es mucho menos precisa:

library(caret)

# Generamos un modelo con una regresión lineal
diam <- lm(price ~ ., diamonds)

# Generamos la predicción
diam_pred <- predict(diam, diamonds)

# Obtenemos RMSE 
error <- diam_pred - diamonds$price
rmse_in <- sqrt(mean(error^2))
rmse_in

[1] 1129.843

Para poner a prueba el modelo con una muestra ajena a la generación del mismo(out-of-sample error), entrenamos al modelo con una fracción de la muestra y la validamos con la fracción restante:

# 80% de la muestra para entrenar y 20% para validar:
# Crea números secuenciales del primer al último elemento
all <- 1:length(diamonds$price)

# Selecciona el 80% para entrenar el modelo
trn <- sample(length(diamonds$price), 
              length(diamonds$price)*0.8)

# Entrena el modelo con la muestra de entrenamiento...
diam2 <- lm(price ~ ., diamonds[trn,])

# ... lo valida con la reserva
diam2_pred <- predict(diam2, diamonds[-trn,])

# ... y obtiene el RMSE
error2 <- diam2_pred - diamonds$price[-trn]
rmse_out <- sqrt(mean(error2^2))
rmse_out

[1] 1136.894

Muestras pequeñas

Estos métodos no son fiables en muestras pequeñas. La presencia de un solo outlier en estas circunstancias puede influir de forma desmesurada en el modelo.

La validación cruzada es un método más robusto para demostrar la bondad de un modelo.

Validación cruzada

La validación cruzada realiza el mismo proceso anterior un número \(n\) de veces; una vez se demuestra la bondad de los \(n\) modelos parciales, se genera un modelo definitivo con toda la muestra.

Para este proceso se utiliza caret::train.

Note

caret también se encarga de que las proporciones de las clases se mantengan constantes en las muestras que se utilizan para generar los sub-modelos (folds).

caret::train

model <- train(price ~ ., diamonds,
  # algoritmo del modelo
  method = "lm",       
       # 'cv' - cross-validation
  trControl = trainControl(
       method = "cv",
       # número de modelos parciales que se crean
       number = 10,
       # esta opción puede estar en TRUE en modo interactivo
       verboseIter = FALSE))
model

Linear Regression 

53940 samples
    9 predictor

No pre-processing
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 48547, 48546, 48545, 48546, 48546, 48546, ... 
Resampling results:

  RMSE      Rsquared   MAE     
  1140.036  0.9183058  741.3911

Tuning parameter 'intercept' was held constant at a value of TRUE

diam_cv <- predict(model, diamonds)
error <- diam_cv - diamonds$price
rmse_cv <- sqrt(mean(error^2))
rmse_cv

[1] 1129.843

1.

Chollet F. Deep learning with R. Second edition. Shelter Island: Manning; 2022.

2.

Gatto L. An Introduction to Machine Learning with R [Internet]. 2020 [cited 2024 Dec 6]. Available from: https://lgatto.github.io/IntroMachineLearningWithR/index.html

3.

Song Y, Millidge B, Salvatori T, Lukasiewicz T, Xu Z, Bogacz R. Inferring neural activity before plasticity as a foundation for learning beyond backpropagation. Nature Neuroscience [Internet]. 2024 Feb [cited 2024 Dec 20];27(2):348–58. Available from: https://www.nature.com/articles/s41593-023-01514-1

4.

Jones E, Harden S, Crawley MJ. The R book. Third edition. Hoboken, NJ: Wiley; 2022.

5.

Field A, Miles J, Field Z. Discovering statistics using R. Repr. Los Angeles, CA, USA: Sage; 2014.