Submuestrear sigue siendo pecado.Ejemplo

estadística
muestreo
2025
Author

José Luis Cañadas Reche

Published

May 16, 2025

Important

El post anterior a este no vale para mucho, puesto que Harrell tenía razón (as usual). Gracias a Carlos Gil que me avisaba de que había algo raro en mi post.

El problema es que el glm inicial no converge. El tema era movidas de configuración en mi slimbook con Linux. Tenía puesto para BLAS y LAPACK cosas incongruentes. Usaba intel mkl para LAPACK y en BLAS la de por defecto

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/libmkl_rt.so;

En problemas como el creditcard dónde hay una separación muy buena entre los 1’s y los 0’s, la estabilidad numérica en el algoritmoi de iteratively reweighted least squares (IWLS) que usa glm es importante. Y esta incoherencia en mi sistema hizo que glm no convergiera al usar un conjunto de train grande, y se dieran coeficientes muy grandes en valor absoluto.

Lo he arreglado y ahora tengo openblas tanto para BLAS como para LAPACK.

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

Y ya no hay problemas. De todas formas la librería rms de Harrell es muy buena y hace un ajuste penalizado si es necesario.

Hago lo mismo que en post anterior, pero para los glms uso la fantástica librería de Harrell

Credit card

El típico ejemplo de creditcard de kagle.

Show the code
# pruebo esto de la funcińo use nueva en r base.

use(package = "pROC", c("roc", "auc"))
use(package = "yardstick", c("mn_log_loss_vec"))

creditcard <- data.table::fread(here::here("data/creditcard.csv"))  |>  as.data.frame()
creditcard$Class  <-  as.factor(creditcard$Class)

table(creditcard$Class)
#> 
#>      0      1 
#> 284315    492

skimr::skim(creditcard)
Data summary
Name creditcard
Number of rows 284807
Number of columns 31
_______________________
Column type frequency:
factor 1
numeric 30
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Class 0 1 FALSE 2 0: 284315, 1: 492

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Time 0 1 94813.86 47488.15 0.00 54201.50 84692.00 139320.50 172792.00 ▃▇▅▆▇
V1 0 1 0.00 1.96 -56.41 -0.92 0.02 1.32 2.45 ▁▁▁▁▇
V2 0 1 0.00 1.65 -72.72 -0.60 0.07 0.80 22.06 ▁▁▁▇▁
V3 0 1 0.00 1.52 -48.33 -0.89 0.18 1.03 9.38 ▁▁▁▁▇
V4 0 1 0.00 1.42 -5.68 -0.85 -0.02 0.74 16.88 ▂▇▁▁▁
V5 0 1 0.00 1.38 -113.74 -0.69 -0.05 0.61 34.80 ▁▁▁▇▁
V6 0 1 0.00 1.33 -26.16 -0.77 -0.27 0.40 73.30 ▁▇▁▁▁
V7 0 1 0.00 1.24 -43.56 -0.55 0.04 0.57 120.59 ▁▇▁▁▁
V8 0 1 0.00 1.19 -73.22 -0.21 0.02 0.33 20.01 ▁▁▁▇▁
V9 0 1 0.00 1.10 -13.43 -0.64 -0.05 0.60 15.59 ▁▁▇▁▁
V10 0 1 0.00 1.09 -24.59 -0.54 -0.09 0.45 23.75 ▁▁▇▁▁
V11 0 1 0.00 1.02 -4.80 -0.76 -0.03 0.74 12.02 ▁▇▁▁▁
V12 0 1 0.00 1.00 -18.68 -0.41 0.14 0.62 7.85 ▁▁▁▇▁
V13 0 1 0.00 1.00 -5.79 -0.65 -0.01 0.66 7.13 ▁▃▇▁▁
V14 0 1 0.00 0.96 -19.21 -0.43 0.05 0.49 10.53 ▁▁▁▇▁
V15 0 1 0.00 0.92 -4.50 -0.58 0.05 0.65 8.88 ▁▇▂▁▁
V16 0 1 0.00 0.88 -14.13 -0.47 0.07 0.52 17.32 ▁▁▇▁▁
V17 0 1 0.00 0.85 -25.16 -0.48 -0.07 0.40 9.25 ▁▁▁▇▁
V18 0 1 0.00 0.84 -9.50 -0.50 0.00 0.50 5.04 ▁▁▂▇▁
V19 0 1 0.00 0.81 -7.21 -0.46 0.00 0.46 5.59 ▁▁▇▂▁
V20 0 1 0.00 0.77 -54.50 -0.21 -0.06 0.13 39.42 ▁▁▇▁▁
V21 0 1 0.00 0.73 -34.83 -0.23 -0.03 0.19 27.20 ▁▁▇▁▁
V22 0 1 0.00 0.73 -10.93 -0.54 0.01 0.53 10.50 ▁▁▇▁▁
V23 0 1 0.00 0.62 -44.81 -0.16 -0.01 0.15 22.53 ▁▁▁▇▁
V24 0 1 0.00 0.61 -2.84 -0.35 0.04 0.44 4.58 ▁▇▆▁▁
V25 0 1 0.00 0.52 -10.30 -0.32 0.02 0.35 7.52 ▁▁▇▂▁
V26 0 1 0.00 0.48 -2.60 -0.33 -0.05 0.24 3.52 ▁▆▇▁▁
V27 0 1 0.00 0.40 -22.57 -0.07 0.00 0.09 31.61 ▁▁▇▁▁
V28 0 1 0.00 0.33 -15.43 -0.05 0.01 0.08 33.85 ▁▇▁▁▁
Amount 0 1 88.35 250.12 0.00 5.60 22.00 77.16 25691.16 ▇▁▁▁▁

El dataset está muy desbalanceado.

Show the code
id_train <- sample(1:nrow(creditcard), 140000)
id_test <- setdiff(1:nrow(creditcard), id_train)


train <- creditcard[id_train, ]
test <-  creditcard[id_test, ]

(t1 <- table(train$Class))
#> 
#>      0      1 
#> 139759    241
(t2 <- table(test$Class))
#> 
#>      0      1 
#> 144556    251

prop.table(t1)
#> 
#>           0           1 
#> 0.998278571 0.001721429
prop.table(t2)
#> 
#>           0           1 
#> 0.998266658 0.001733342

Modelos glm normal, con submuestreo y con pesos

Modelo normal sin submuestrear

No uso la variable Time, ni me caliento la cabez en buscar interacciones (por el momento)

Show the code

f1  <-  as.formula("Class ~ V1 + V2 + V3 + V4 + V5 +V6 + V7 + V8 + V9 + V10 + V11 + V12 +V13 +V14 + Amount")

glm1 <- glm(f1, family = binomial, data = train)

## no ha ajustado y coeficentes muy altos y bajos, inestabilidad
# señal de que algo ha ido muy muy mal
summary(glm1)
#> 
#> Call:
#> glm(formula = f1, family = binomial, data = train)
#> 
#> Coefficients:
#>               Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -8.371e+00  1.836e-01 -45.601  < 2e-16 ***
#> V1          -2.118e-02  4.764e-02  -0.445  0.65657    
#> V2           2.052e-02  6.121e-02   0.335  0.73740    
#> V3           1.651e-01  6.234e-02   2.648  0.00809 ** 
#> V4           5.896e-01  8.560e-02   6.887 5.68e-12 ***
#> V5          -5.366e-02  7.350e-02  -0.730  0.46531    
#> V6          -1.772e-01  1.085e-01  -1.633  0.10246    
#> V7          -1.104e-02  5.857e-02  -0.189  0.85044    
#> V8          -1.576e-01  3.898e-02  -4.043 5.28e-05 ***
#> V9          -8.942e-03  1.184e-01  -0.075  0.93982    
#> V10         -4.673e-01  9.690e-02  -4.822 1.42e-06 ***
#> V11          1.389e-01  1.012e-01   1.372  0.16994    
#> V12         -2.502e-01  8.240e-02  -3.036  0.00240 ** 
#> V13         -2.786e-01  1.041e-01  -2.675  0.00747 ** 
#> V14         -6.445e-01  7.674e-02  -8.399  < 2e-16 ***
#> Amount      -6.062e-05  3.888e-04  -0.156  0.87610    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 3549.3  on 139999  degrees of freedom
#> Residual deviance: 1132.8  on 139984  degrees of freedom
#> AIC: 1164.8
#> 
#> Number of Fisher Scoring iterations: 11

(roc_glm <- roc(test$Class, predict(glm1, newdata = test, type = "response")))
#> 
#> Call:
#> roc.default(response = test$Class, predictor = predict(glm1,     newdata = test, type = "response"))
#> 
#> Data: predict(glm1, newdata = test, type = "response") in 144556 controls (test$Class 0) < 251 cases (test$Class 1).
#> Area under the curve: 0.9772



(logloss_glm <- mn_log_loss_vec(as.factor(test$Class), predict(glm1, newdata = test, type = "response"), event_level = "second"))
#> [1] 0.004255575

# logloss(ytrue, predict(glm1, newdata = test, type = "response"))

Usando la libreríarms que da mucha más información del modelo ajustado, como ínidices de discriminación o C-index.

Show the code
library(rms)

dd <- datadist(train)
options(datadist = "dd")

modelo_rms <- lrm(Class ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 +
                    V9 + V10 + V11 + V12 + V13 + V14 + Amount,
                  data = train,
                  x = TRUE, y = TRUE)

modelo_rms
#> Logistic Regression Model
#> 
#> lrm(formula = Class ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + 
#>     V9 + V10 + V11 + V12 + V13 + V14 + Amount, data = train, 
#>     x = TRUE, y = TRUE)
#> 
#>                         Model Likelihood        Discrimination    Rank Discrim.    
#>                               Ratio Test               Indexes          Indexes    
#> Obs         140000    LR chi2    2416.53        R2       0.684    C       0.978    
#>  0          139759    d.f.            15    R2(15,140000)0.017    Dxy     0.956    
#>  1             241    Pr(> chi2) <0.0001     R2(15,721.8)0.964    gamma   0.956    
#> max |deriv| 0.0005                              Brier    0.001    tau-a   0.003    
#> 
#>           Coef    S.E.   Wald Z Pr(>|Z|)
#> Intercept -8.3712 0.1836 -45.60 <0.0001 
#> V1        -0.0212 0.0476  -0.44 0.6566  
#> V2         0.0205 0.0612   0.34 0.7374  
#> V3         0.1651 0.0623   2.65 0.0081  
#> V4         0.5896 0.0856   6.89 <0.0001 
#> V5        -0.0537 0.0735  -0.73 0.4653  
#> V6        -0.1772 0.1085  -1.63 0.1025  
#> V7        -0.0110 0.0586  -0.19 0.8504  
#> V8        -0.1576 0.0390  -4.04 <0.0001 
#> V9        -0.0089 0.1184  -0.08 0.9398  
#> V10       -0.4673 0.0969  -4.82 <0.0001 
#> V11        0.1389 0.1012   1.37 0.1699  
#> V12       -0.2502 0.0824  -3.04 0.0024  
#> V13       -0.2786 0.1041  -2.68 0.0075  
#> V14       -0.6445 0.0767  -8.40 <0.0001 
#> Amount    -0.0001 0.0004  -0.16 0.8761


(roc_glm_harrel <- roc(test$Class, predict(modelo_rms, newdata = test, type = "fitted")))
#> 
#> Call:
#> roc.default(response = test$Class, predictor = predict(modelo_rms,     newdata = test, type = "fitted"))
#> 
#> Data: predict(modelo_rms, newdata = test, type = "fitted") in 144556 controls (test$Class 0) < 251 cases (test$Class 1).
#> Area under the curve: 0.9772

Modelo glm submuestreado

Vamos a submuestrear la clase 0, pero no demasiado, nos quedamos con 3000

Show the code

samples_0  <-  train[sample(rownames(train[train$Class == "0",]), 3000),]

table(samples_0$Class)
#> 
#>    0    1 
#> 3000    0

train_subsample  <- rbind(samples_0, train[train$Class=="1", ])

dd_sub <- datadist(train_subsample)
options(datadist = "dd_sub")
    
table(train_subsample$Class)
#> 
#>    0    1 
#> 3000  241

glm2  <-  lrm(f1, data = train_subsample, x = TRUE,  y = TRUE)

(roc_glm_submuestreo <- pROC::roc(test$Class, predict(glm2, newdata = test, type = "fitted")))
#> 
#> Call:
#> roc.default(response = test$Class, predictor = predict(glm2,     newdata = test, type = "fitted"))
#> 
#> Data: predict(glm2, newdata = test, type = "fitted") in 144556 controls (test$Class 0) < 251 cases (test$Class 1).
#> Area under the curve: 0.9791

(logloss_glm_submuestreo <- mn_log_loss_vec(as.factor(test$Class), predict(glm2, newdata = test, type = "fitted"), event_level = "second"))
#> [1] 0.01903171

Y en este caso , submuestrear es básicamente igual que no hacerlo.

Modelo glm con pesos

Show the code

# Pesos inversamente proporcionales al tamaño de la clase
pesos <- ifelse(train$Class == "1", sum(train$Class == "0")/sum(train$Class == "1"), 1)

# Modelo con pesos 
# glm_pesos <- glm(f1, family = binomial, data = train, weights = pesos)
# 
# (roc_glm_pesos <- pROC::roc(test$Class, predict(glm_pesos, newdata = test, type = "response")))
# 
# (logloss_glm_pesos<- mn_log_loss_vec(as.factor(test$Class), predict(glm_pesos, newdata = test, type = "response"), event_level = "second"))

options(datadist = "dd")

modelo_rms_pesos <- lrm(Class ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 +
                            V9 + V10 + V11 + V12 + V13 + V14 + Amount,
                        data = train,
                        weights = pesos,
                        x = TRUE, y = TRUE)


(roc_rms_pesos <- pROC::roc(test$Class, predict(modelo_rms_pesos, newdata = test, type = "fitted")))
#> 
#> Call:
#> roc.default(response = test$Class, predictor = predict(modelo_rms_pesos,     newdata = test, type = "fitted"))
#> 
#> Data: predict(modelo_rms_pesos, newdata = test, type = "fitted") in 144556 controls (test$Class 0) < 251 cases (test$Class 1).
#> Area under the curve: 0.9811

(logloss_rms_pesos <- mn_log_loss_vec(as.factor(test$Class), predict(modelo_rms_pesos, newdata = test, type = "fitted"), event_level = "second"))
#> [1] 0.1139893

Pues amigos, tanto con submuestreo, pesos o sin hacer nada, la regresión logística ha funcionado estupendamente aún estando el dataset tan desbalanceado.

De hecho

Show the code
print(glue::glue("Auc sin submuestrear: {round(roc_glm_harrel$auc, 3)}"))
#> Auc sin submuestrear: 0.977
print(glue::glue("Auc muestreando: {round(roc_glm_submuestreo$auc, 3)}"))
#> Auc muestreando: 0.979
print(glue::glue("Auc glm con pesos: {round(roc_rms_pesos$auc, 3)}"))
#> Auc glm con pesos: 0.981

Coda

Tuve el atrevimiento de pensar que Harrell se equivocaba, iluso de mi. El error en la confi de mi sistema me la jugó. Cuando uno lee a Harrell hay que tomarse muy en serio todo lo que diga, aunque a veces parezca ser contraintuitivo.

No volveré a cometer semejante error. Buen finde