Intervalos de confianza con confint()

• Bajo normalidad de \(u\) (o con \(n\) grande), \(\hat\beta \sim N(\beta, Var(\hat\beta))\), y el intervalo al 95% es \(\hat\beta \pm 1{,}96 \cdot se(\hat\beta)\). R lo da hecho:

modelo <- lm(ventas ~ renta + descuento + zona, data = tienda)
confint(modelo)                      # todos los coeficientes
confint(modelo, "renta", level = 0.95)

• A mano, con el valor crítico de la normal (qnorm) o de la t de Student (qt, que es lo que usa confint()):

b  <- coef(modelo)["renta"]
se <- summary(modelo)$coefficients["renta", "Std. Error"]
b + qnorm(c(0.025, 0.975)) * se

• Para la varianza del error, \((n - k)\, s^2 / \sigma^2 \sim \chi^2_{n-k}\), con \(s^2\) = summary(modelo)$sigma^2, \(n\) = nobs(modelo) y \(k\) el número de coeficientes:

n <- nobs(modelo); k <- length(coef(modelo)); s2 <- summary(modelo)$sigma^2
(n - k) * s2 / qchisq(c(0.975, 0.025), df = n - k)

• Con los datos del tema: renta 0,380 a 0,606 por fórmula, 0,370 a 0,604 por bootstrap. Las dos cosas miden lo mismo; el bootstrap no necesita el supuesto

Intervalos de confianza y de predicción

• predict() puede dar, además del valor esperado, dos intervalos distintos:

nuevos <- tibble(renta = 2000, descuento = 1,
                 zona = factor("ciudad",
                               levels = levels(tienda$zona)))
predict(modelo, newdata = nuevos, interval = "confidence")
predict(modelo, newdata = nuevos, interval = "prediction")

• Confianza ("confidence"): dónde está el gasto medio de todos los clientes con esas características; solo recoge la incertidumbre de los coeficientes

• Predicción ("prediction"): dónde caerá el gasto de un cliente concreto; añade la varianza de \(u\), y por eso es mucho más ancho. Es el que responde a “cuánto gastará este cliente”

rsample::bootstraps() y boot::boot()

• El bucle del tema es la forma de entenderlo. Para usarlo a diario hay funciones que lo hacen

library(rsample)
remuestras <- bootstraps(tienda,
                         times = 1000)
coefs <- map_dbl(remuestras$splits,
  function(s) {
    m <- lm(ventas ~ renta + descuento +
              zona, data = as_tibble(s))
    coef(m)["renta"]
  })
sd(coefs)
quantile(coefs, c(0.025, 0.975))

library(boot)
coef_renta <- function(datos, filas) {
  m <- lm(ventas ~ renta + descuento +
            zona, data = datos[filas, ])
  coef(m)["renta"]
}
boot(data = tienda,
     statistic = coef_renta, R = 1000)

• En el Tema 07, rsample también parte la muestra (initial_split()) y crea los bloques de la validación cruzada (vfold_cv())

Variantes de la validación cruzada

• K bloques (la del tema): K = 5 o 10. Cada observación hace de prueba una vez

• Dejar uno fuera (leave-one-out, LOOCV): K = \(n\). Se estima \(n\) veces con \(n - 1\) observaciones y se evalúa en la que falta. Sin azar, pero caro, y con mucha varianza entre iteraciones

• Particiones aleatorias repetidas: se repite muchas veces la partición entrenamiento/prueba con reemplazamiento del sorteo; una observación puede caer varias veces en prueba

• Validación cruzada repetida: K bloques, varias veces con distintas semillas, y se promedia todo. Es lo más estable

• Con rsample:

vfold_cv(tienda, v = 5)                    # 5 bloques
vfold_cv(tienda, v = 5, repeats = 10)      # 5 bloques, diez veces
loo_cv(tienda)                             # dejar uno fuera

• Regla: la validación cruzada elige el modelo o el hiperparámetro; el error que se publica se mide en una muestra de prueba apartada antes de todo

Validación cruzada de K bloques, paso a paso

set.seed(12345678)
bloque <- sample(rep(1:5,
  length.out = nrow(tienda)))
errores <- numeric(5)
for (k in 1:5) {
  entr <- tienda[bloque != k, ]
  prue <- tienda[bloque == k, ]
  m <- lm(ventas ~ poly(edad, 3),
          data = entr)
  pred <- predict(m, prue)
  e <- prue$ventas - pred
  errores[k] <- mean(e^2)
}
mean(errores)

Validación cruzada aleatoria y dejar uno fuera

Muestras de entrenamiento y de prueba

set.seed(12345678)
n <- nrow(tienda)
filas <- sample(n, size = round(0.7 * n))
entrena <- tienda[filas, ]
prueba  <- tienda[-filas, ]

Ajuste insuficiente, correcto y sobreajuste

• Figuras del libro ISL (An Introduction to Statistical Learning, figuras 2.9 a 2.11). Los grados de libertad, el número de valores del modelo que pueden variar, resumen la flexibilidad de una curva

Trade-off varianza-sesgo

• El MSE esperado en la muestra de prueba se descompone en:

\[ \small E\left[\left(y-\hat{f}(x)\right)^2\right] = \underbrace{\left[E\left(\hat{f}(x)\right)-f(x)\right]^2}_{\text{sesgo}^2} + \underbrace{E\left(\left[\hat{f}(x)-E\left(\hat{f}(x)\right)\right]^2\right)}_{\text{varianza}} + \text{Var}(u) \]

Sesgo y varianza con simulaciones

• 50 muestras de la relación conocida del tema y un ajuste de cada grado en cada una:

• Grado 1: las 50 rectas caen casi en el mismo sitio y todas lejos de la curva: sesgo. Grado 2: abrazan la curva. Grado 12: cada muestra da una curva distinta: varianza

• Con la tienda, 200 remuestras bootstrap y la predicción para un cliente de 65 años:

Tablas de resultados en Quarto

• tidy() (paquete broom) convierte la salida de lm() en una tabla, que se puede manipular y presentar con kable(). La celda necesita #| results: markup

library(broom); library(kableExtra)
modelo |> tidy() |>
  select(term, estimate, std.error) |>
  kbl(digits = 3) |> kable_classic()

• modelsummary() pone varios modelos en columnas, con estrellas de significatividad y las medidas de ajuste que se pidan:

library(modelsummary)
modelos <- list(
  "Lineal"      = lm(ventas ~ renta + descuento + zona, data = tienda),
  "Logaritmo"   = lm(ventas ~ log(renta) + descuento + zona, data = tienda),
  "Parábola"    = lm(ventas ~ renta + descuento + zona +
                       edad + I(edad^2), data = tienda),
  "Interacción" = lm(ventas ~ renta * mujer + descuento + zona, data = tienda)
)
modelsummary(modelos, stars = TRUE,
             gof_map = c("nobs", "r.squared", "adj.r.squared", "rmse"))

• Es la forma habitual de presentar la tabla de la entrega 2 y las comparaciones del trabajo final

Errores estándar robustos y las “pruebas de supuestos”

• La fórmula del error estándar de summary() supone varianza constante de \(u\). Si varía con las \(X\) (heterocedasticidad, lo habitual con gasto o renta), los coeficientes siguen bien y los errores estándar no. Se corrigen con errores estándar robustos:

library(sandwich); library(lmtest)
coeftest(modelo, vcov = vcovHC(modelo, type = "HC1"))
modelsummary(modelo, vcov = "HC1")         # lo mismo dentro de la tabla

• El bootstrap de filas del tema también es robusto a esto

• A menudo se proponen “pruebas” de los supuestos del modelo lineal (normalidad, linealidad, colinealidad, heterocedasticidad). En este curso se tratan así:

  • No linealidad: se modeliza (Parte III del tema), no se contrasta

  • No normalidad: innecesaria con \(n\) grande o con bootstrap

  • Colinealidad: dos variables que miden casi lo mismo; se quita una o se acepta el error estándar grande

  • Heterocedasticidad: errores estándar robustos

Métricas de la matriz de confusión

• Logit del tema, umbral 0,5: precisión 7/8 = 0,875 (de 8 contactados, 7 tenían descuento) y sensibilidad 0,117. Las dos miran los positivos predichos desde lados distintos. En yardstick (Tema 07): precision(), f_meas(), kap(), bal_accuracy()

Curva precisión-sensibilidad y otras pérdidas

• Con clases muy desequilibradas (fraude al 1%), la curva ROC puede parecer buena porque la tasa de falsos positivos se calcula sobre muchísimos negativos. La curva precisión-sensibilidad (PR) y su área (PR-AUC) miran solo lo que pasa con los positivos, y son más exigentes

library(yardstick)
logit <- glm(descuento ~ renta + edad + zona,
             data = tienda, family = "binomial")
datos <- tienda |>
  mutate(prob = predict(logit, type = "response"),
         descuento = factor(descuento, levels = c(1, 0)))
pr_curve(datos, truth = descuento, prob)
pr_auc(datos, truth = descuento, prob)
roc_curve(datos, truth = descuento, prob)
roc_auc(datos, truth = descuento, prob)

• Otras funciones de pérdida que no fijan umbral: la log-loss (menos la log-verosimilitud media), que penaliza dar probabilidad alta a la clase equivocada, y la curva de ganancia, que ordena por probabilidad y dice qué fracción de positivos se captura contactando al x% con mayor probabilidad

• Con más de dos clases se calcula una ROC por clase (esa clase frente a todas las demás) y se promedian, con o sin ponderar por el tamaño de cada clase (roc_aunu(), roc_aunp())

Curva ROC: cada punto es un umbral

Verosimilitud y log-odds

• El logit no minimiza la SCR: elige los \(\beta\) que hacen más probable observar los unos y ceros de la muestra (máxima verosimilitud):

\[ \small L(\beta) = \prod_{i:\, y_i = 1} p(x_i) \prod_{i:\, y_i = 0} \left(1 - p(x_i)\right), \qquad p(x) = \Lambda(\beta_0 + \beta_1 x_1 + \ldots) \]

• La probabilidad relativa (odds) es \(\dfrac{p}{1 - p} = e^{\beta_0 + \beta_1 x_1 + \ldots}\), así que el log-odds es lineal: \(\beta_j\) es el cambio en el logaritmo de la probabilidad relativa por una unidad de \(x_j\), y \(e^{\beta_j}\) el factor por el que se multiplica la odds

• El efecto sobre la probabilidad depende del punto: \(\partial p / \partial x_j = \beta_j\, p\, (1 - p)\). Máximo en \(p = 0{,}5\) y casi nulo en los extremos. Un resumen habitual es el efecto medio sobre los clientes de la muestra:

p <- predict(logit, type = "response")
mean(coef(logit)["renta"] * p * (1 - p))

Más de dos clases

• Con una variable de \(C\) categorías (producto elegido, tipo de contrato), el logit multinomial tiene un índice lineal por clase:

\[ \small \Pr(Y = c \mid X) = \frac{e^{\beta_{0c} + \beta_{1c} x_1 + \ldots}}{\sum_{l = 1}^{C} e^{\beta_{0l} + \beta_{1l} x_1 + \ldots}} \]

• Dos implementaciones. nnet::multinom() con fórmula; glmnet con matrices (es lo que usa el Tema 06 para regularizar):

library(nnet)
multinom(zona ~ renta + edad, data = tienda)

library(glmnet)
x <- model.matrix(zona ~ renta + edad, data = tienda)[, -1]
mod <- glmnet(x, tienda$zona, family = "multinomial", lambda = 0)
predict(mod, newx = x, type = "response")   # una probabilidad por clase
predict(mod, newx = x, type = "class")      # la clase más probable

• La clase predicha es la de mayor probabilidad; las métricas se calculan clase a clase o promediadas

Regresión lineal: comentarios finales

• Estos modelos informan de la relación entre los regresores y la variable dependiente. Correlación no implica causalidad, salvo en ensayos aleatorios controlados o pruebas A/B

• La selección de variables es importante, porque evita sesgos y afecta a la variabilidad:

lm(ventas ~ descuento, data = tienda) |> summary()
lm(ventas ~ descuento + renta, data = tienda) |> summary()
lm(ventas ~ renta, data = tienda) |> summary()
lm(ventas ~ renta + educ, data = tienda) |> summary()

• Los valores atípicos pueden tratarse sin afectar a la validez del modelo

Lo que hizo la tienda: elegir a mano

ggplot(tienda, aes(renta, ventas,
  colour = factor(descuento))) +
  geom_point(alpha = 0.5) +
  scale_colour_manual(
    values = c("grey55", "red"),
    labels = c("sin descuento",
               "con descuento")) +
  labs(colour = NULL) + theme_minimal()

Lo que arregla la regresión, y lo que no

• Con la renta en el modelo, el coeficiente compara clientes de la misma renta: lm(ventas ~ descuento + renta) da +320 euros
• La regresión iguala los grupos en lo que está en la tabla. La tienda eligió a mano, y pudo mirar cosas que la tabla no tiene: fidelidad, quejas, lo que sabía el vendedor
  • si alguna de esas cosas también mueve el gasto, el coeficiente sigue mezclando
  • con datos observacionales nunca se sabe que no queda nada fuera
• Por eso un coeficiente predictivo, por preciso que sea, responde a la pregunta 1 y solo a veces a la 2

Repartir al azar: el experimento

set.seed(2027)
dif <- numeric(1000)
for (b in 1:1000) {
  azar <- sample(tienda$descuento)
  dif[b] <- mean(tienda$renta[azar == 1]) -
            mean(tienda$renta[azar == 0])
}
ggplot(tibble(dif), aes(dif)) +
  geom_histogram(bins = 30) +
  geom_vline(xintercept = -839,
             colour = "red") +
  theme_minimal()

El sesgo por variable omitida

• Modelo verdadero \(y = \beta_0 + \beta_1 x + \beta_2 z + u\), y se estima sin \(z\). El coeficiente de \(x\) que sale estima:

\[ \small \beta_1 + \beta_2 \cdot \delta, \qquad \text{donde } \delta \text{ es la pendiente de } z \text{ sobre } x \]

• El signo del sesgo es el producto de los signos: efecto de \(z\) sobre \(y\) (\(\beta_2\)) por relación entre \(z\) y \(x\) (\(\delta\))

• El caso del tema: \(x\) = descuento, \(z\) = renta. La renta sube el gasto (\(\beta_2 > 0\)) y el descuento fue a rentas bajas (\(\delta < 0\)): sesgo negativo, y el +500 se queda en -97 al omitir la renta

• Lo mismo con el modelo completo: si la tienda eligió además por algo que no está en la tabla y que se relaciona con el gasto, el coeficiente sigue sesgado. Con datos observacionales nunca se sabe que no queda nada fuera

Ensayos aleatorios y pruebas A/B

• La única forma de que \(\delta = 0\) para todas las variables, observadas o no, es que \(x\) se asigne al azar: ensayo aleatorio controlado. En empresa, prueba A/B: la mitad de los clientes (al azar) recibe el descuento, la otra mitad no, y se compara el gasto

• Con asignación al azar, lm(ventas ~ descuento) sí estima el efecto: los dos grupos son iguales en todo lo demás, en promedio

• Lo que puede fallar incluso así: que los clientes no usen el descuento (incumplimiento), que se vayan antes de medir (desgaste), que el descuento de unos afecte a otros (interferencia)

• Sin experimento hay métodos para aproximarse (diferencias en diferencias, variables instrumentales, discontinuidades) que son materia de Econometría. En este curso la regla es la del tema: un coeficiente predictivo no dice qué pasaría si cambiamos la variable

Un ejemplo con otros datos: descuentos y ventas

• discount.csv (página Datos): ventas y porcentaje medio de descuento de 404 clientes, con su renta

datos <- read_csv("data/discount.csv")
lm(sales ~ discount, data = datos) |> coef()
lm(sales ~ discount + income, data = datos) |> coef()
datos |> mutate(renta_baja = income < 7500) |>
  ggplot(aes(x = discount, y = sales)) + geom_point() +
  geom_smooth(method = "lm") + facet_wrap(~ renta_baja)

• La misma estructura que en la tienda: los descuentos mayores fueron a rentas bajas. Separando por renta, la pendiente cambia de signo

• Es un ejercicio de cinco minutos para comprobar que se ha entendido la regla del signo