Métodos Estadísticos

• La modelización mediante métodos estadísticos permite

  • Encontrar patrones

  • Interpretar datos

tibble(x = runif(n=1e7, min = 0, max = 20)) %>% pull(x) %>% mean()

set.seed(9915)
tibble(x = runif(n=25, min = 0, max = 20)) %>% pull(x) %>% mean()

Procedimiento Bootstrap (cont.)

distribDF %>% ggplot(aes(x=medias)) + geom_density()

# Error estándar de la media
mediaSE <- distribDF %>% pull(medias) %>% sd()

# Intervalos de confianza al 95% para la media
mediaIC <- distribDF %>% pull(medias) %>% quantile(c(0.025, 0.975))

Aprendizaje supervisado (cont.)

• Dos formas de especificar la relación \(f\)

  • Modelo paramétrico: supone un forma de \(f\) que depende de parámetros desconocidos, p.e., lineal \(f(x) =\beta_0 + \beta_1 x_1 + \dots + \beta_k x_k\)

  • Modelo no paramétrico: ajustar \(f\) a los datos sin supuestos funcionales

    • más flexibilidad y mejor ajuste, pero más difícil de estimar e interpretar

• La capacidad predictiva de un modelo mejora si incluimos más variables explicativas (modelo más flexible)

• Dependiendo de la naturaleza de la variable dependiente, tenemos:

  1.- Problemas de Regresión: variable de respuesta cuantitativa (toma valores numéricos)

  2.- Problemas de Clasificación: variable de respuesta cualitativa (toma valores en una de \(C\) categorías o clases)

Regresión Lineal: Estimación

• Objetivo: estimar los coeficientes poblacionales desconocidos a partir de una muestra

• El error de estimación o residuo es \(\small \hat{e}_i = y_i - \hat{y}_i\), donde la predicción a partir del modelo estimado es \(\small \hat{y}_i=\hat{\beta}_0+\hat{\beta}_1 X_1+ \dots + \hat{\beta}_k X_p\)

• Los coeficientes estimados son los que minimizan la Suma Cuadrática de Residuos: la suma total de distancias entre los datos observados y predichos

\[ \small SCR=\sum_{i=1}^{n} \hat{e}_i^2= \sum_{i=1}^{n} ( y_i - \hat{y}_i)^2 \]

• Esto equivale a las condiciones derivadas de suponer \(\small E(\varepsilon|X)=0\)

Regresión Lineal: Estimación (cont.)

descuento <- read.csv("data/descuento.csv") %>% 
        mutate(zona = parse_factor(zona), mujer = as.factor(mujer))
modelo <- lm(data = descuento, ventas ~ renta + descuento)
(sum.modelo <- summary(modelo))
nobs(modelo)

pred <- predict(modelo, type = "response", 
                se.fit = T, interval = "confidence")
cbind(descuento$ventas, pred$fit, pred$se.fit) %>% head()

Ampliando el Modelo de Regresión Lineal

• Para incorporar información cualitativa en el modelo, se debe incluir una variable binaria para cada categoría, excepto una (grupo de referencia)

  • Los coeficientes de cada variable binaria son el efecto diferencial de ese grupo respecto al grupo de referencia

  • Al incluir un factor, se crean automáticamente las variables binarias

lm(data = descuento, ventas ~ zona) %>% summary()
lm(data = descuento, ventas ~ renta + descuento + zona) %>% summary()

Regresión Lineal: superando la linealidad

lm(data=descuento, log(ventas) ~ log(renta) + descuento) %>% summary()
lm(data=descuento, log(ventas) ~ renta + descuento) %>% summary()

lm(data=descuento, ventas ~ edad + I(edad^2) ) %>% summary()
ggplot(descuento, aes(x=edad, y=ventas)) + geom_point() + 
  geom_smooth(method = "lm", formula = y ~ x + I(x^2))

lm(data=descuento, ventas ~ 
     cut(edad, seq(20, 60, 5), include.lowest = T)) %>% summary()
ggplot(descuento, aes(x=edad, y=ventas)) + geom_point() + 
  geom_smooth(method = "lm", 
              formula = y ~ cut(x, seq(20, 60, 5), include.lowest = T))

Regresión Lineal: superando la linealidad (cont.)

• Incluyendo interacciones entre variables, el efecto de un regresor dependerá de otro regresor: la “pendiente” cambia con el valor de la otra variable

  • Al incluir una interacción siempre deben incluirse los factores principales (NO sólo edad:renta)

lm(data=descuento, ventas ~ edad*renta)  %>% summary()
lm(data=descuento, ventas ~ renta*mujer) %>% summary()

ggplot(descuento, aes(x=renta, y=ventas, color = mujer)) + 
  geom_point() + geom_smooth(method = "lm", formula = y ~ x,
                               se = FALSE)

Regresión Lineal: comentarios finales

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

Regresión Logística

• Un modelo de regresión lineal para respuestas binarias no es adecuado porque genera predicciones fuera del rango \(\small [0,1]\)

• De manera que \(\small \Pr(Y=1|X)= f(x)= \color{blue}{\Lambda( \beta_0 + \beta_1 x_1 + \ldots + \beta_k x_k)} \in [0,1]\)

Regresión Logística (cont.)

• La regresión logística se puede estimar en R con la función glm(), de forma similar a lm(); podemos incluir variables explicativas tanto cuantitativas como cualitativas, interacciones, etc.

censo <- read_csv("data/census.csv") %>%
  mutate(across(where(is.character), ~parse_factor(.x)))

logit <- glm(data = censo, income ~  race + age + I(age^2) + log(hours_per_week)*sex + 
      education_1  + capital_gain + capital_loss, 
         family = "binomial") 
logit %>%  summary()

• Como en el modelo lineal, podemos calcular predicciones de la probabilidad tanto en los mismos datos como en unos nuevos

predict(logit, type="response")

logit2 <- logit <- glm(data = censo, income ~  age + education_1, 
         family = "binomial")
summary(logit2)
predict(logit2, type="response",
        newdata = tibble(age=c(20,60), education_1=c(9,16)))

Regresión Logística con más de dos clases

• La regresión logística se puede generalizar a situaciones con múltiples clases (modelos multinomiales) con un índice lineal para cada clase \[ \small \Pr(Y=c|X)=\frac{e^{\beta_{0c}+\beta_{1c}X_1+\dots+\beta_{kc}X_k}}{\sum_{l=1}^{C}e^{\beta_{0l}+\beta_{1l}X_1+\dots+\beta_{kl}X_k}} \]

• La librería glmnet() permite la estimación de estos modelos

health <- read_csv("data/health_insurance.csv") %>% 
  mutate(across(where(is.character), ~parse_factor(.x)))

library(glmnet)
x <- model.matrix(product ~ age + gender, data = health)
mod.glmnet <- glmnet(x = x, y = health$product, family = "multinomial", 
                     lambda = 0, type.multinomial = "grouped")
coef(mod.glmnet) 
predict(mod.glmnet, newx=x, type = "response")  # probabilidad de cada clase
predict(mod.glmnet, newx=x, type = "class")     # clase

Error de predicción

• Un modelo es mejor si sus predicciones se ajusten mejor a las observaciones

• El error de predicción es \(y - \widehat{y} = f(X) - \widehat{f}(X) + \varepsilon\)

  • \(f - \widehat{f}\) = error reducible (eligiendo modelo)

  • \(\varepsilon\) = error irreducible (variables no observadas)

• La función de pérdida (o coste) evalúa cómo valoramos las desviaciones

Métricas de error de predicción (cuantitativa)

• Mean Square Error (Error Cuadrático Medio): \(\small MSE(y,\widehat{y})={\frac{1}{n}\sum_{i=1}^{n}\left(y-\widehat{y}\right)^2}\)

  • penaliza grandes desviaciones

  • \(\small R^2\) y \(\small R^2\)-ajustado son variantes del MSE, pero solo sirven para comparar modelos con la misma variable dependiente.

• Root Mean Square Error: \(\small RMSE(y,\widehat{y})=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(y-\widehat{y}\right)^2}\)

  • mismas unidades que \(\small y\)

• Mean Absolute Error: \(\small MAE(y,\widehat{y})=\frac{1}{n}\sum_{i=1}^{n}\left|y-\widehat{y}\right|\)

• Otras medidas basadas en distintas funciones de pérdida, la verosimilitud del modelo (\(\small AIC\), \(\small BIC\)), etc

Seleccionar el mejor modelo

• ¿Podemos predecir el número de visitantes en función de la temperatura?

library(mosaicData)
data(RailTrail)
RailTrail %>% ggplot(aes(x = hightemp, y = volume)) + 
  geom_point() +  geom_smooth()

• ¿Cuál es el mejor modelo?

  • \(\small volume = \beta_0 + \beta_1 hightemp + \varepsilon\)

  • \(\small volume = \beta_0 + \beta_1 hightemp + \beta_2 hightemp^2 + \varepsilon\)

  • …

  • \(\small volume = \beta_0 + \beta_1 hightemp + \dots + \beta_{22} hightemp^{22} + \varepsilon\)

RailTrail %>% ggplot(aes(x = hightemp, y = volume)) + 
  geom_point() + 
  geom_smooth(method = 'lm', formula = y ~ poly(x,22) ) + 
  coord_cartesian(ylim = c(100,750))

Muestras de entrenamiento y de prueba

• Las métricas de error (ej., \(\small MSE\)) suelen calcularse para predicciones de los mismos datos usados para ajustar/estimar el modelo (in-sample prediction)

  • Esta muestra se denomina muestra de entrenamiento (training sample)

• PERO queremos saber qué tal se predicen casos nuevos (out-sample prediction)

• Usar las métricas en muestras de entrenamiento implica problemas de “overfitting”: sobreajuste a las características de la muestra concreta

  • Un modelo menos flexible podría tener menor error de predicción con casos nuevos

• Debemos calcular las métricas de error con observaciones que el modelo NO ha usado antes: muestra de prueba (test sample)

“Overfitting”

• Siempre que aumenta la flexibilidad, el MSE

  • disminuye en la muestra de entrenamiento

  • tiene forma de U en la muestra de prueba

• Nota: el MSE en entrenamiento es siempre menor que en prueba

MSE en la muestra de prueba

\[ \small E\left[\left(y-\widehat{f}(x)\right)^2\right] = E\left[\left(f(x)-\widehat{f}(x) + \varepsilon + E\left[\widehat{f}(x)\right]-E\left[\widehat{f}(x)\right] \right)^2\right] = \]

\[ \small =\underbrace{\left[E\left(\widehat{f}(x)\right)-f(x)\right]^2}_{(1)} + \underbrace{E\left(\left[\widehat{f}(x)-E\left(\widehat{f}(x)\right)\right]^2\right)}_{(2)}+Var(\varepsilon) \]

• \(\small (1)=\left[Sesgo\left(\widehat{f}(x)\right)\right]^2\): error por supuestos erróneos en \(f\)

  • ajuste insuficiente (underfit) al perder relaciones relevantes entre \(X\) e \(Y\)

• \(\small (2)=Var\left(\widehat{f}(x)\right)\): sensibilidad a fluctuaciones en el entrenamiento

  • si el algoritmo modela puro ruido en entrenamiento, ajustará bien allí, pero predecirá mal casos nuevos (overfit)

“Trade-off” Varianza–Sesgo

• El sesgo se reduce y la varianza aumenta con la complejidad del modelo \(\Rightarrow\) encontrar un método (ej., flexibilidad) para el que ambos sean bajos

• NO es posible minimizar simultáneamente ambas fuentes de error: memorización (en entrenamiento) vs. generalización de resultados

Medir el Error en la Clasificación

• Los modelos de clasificación NO predicen directamente la categoría, sino la probabilidad de que una observación pertenezca a cada categoría

• La clase predicha será aquella con mayor probabilidad. En el caso binario, implicar superar el umbral de 0.5 (se deben probar varios valores)

censo <- read_csv("data/census.csv") %>%
  mutate(income = parse_factor(income))
logit <- glm(data = censo, income ~ capital_gain,  
             family = "binomial")
prob.predict <- predict(logit, type = "response")

umbral <- 0.5
cat.predict  <- if_else(prob.predict > umbral, 1, 0) 
cbind(censo$income, cat.predict, prob.predict) %>% head(10)

• Como no tiene sentido diferencia de clases (variables categóricas), NO se pueden calcular medidas como el MSE y otros relacionados

Matriz de Confusión

• Tabular categorías observadas frente a las categorías predichas

table(cat.predict, censo$income)

• Existen varias medidas derivadas de la matriz de confusión

Métricas con la matriz de confusión

• Tasa de observaciones correctamente clasificadas (exactitud o accuracy)

\[ \scriptsize ACCUR=\frac{VP+VN}{VP+FP+VN+FN} \]

• No es informativo cuando algunas clases son infrecuentes (datos desequilibrados)

  • si hay poco fraude/enfermos (ej., 5%), predecir que nunca hay fraude implica \(\scriptsize ACCUR=95\%\), PERO NO detecta fraude/enfermedad

• El estadístico Kappa (\(\small \kappa\)) es una medida similar, pero que ajusta por lo se esperaría solo por azar (corrigiendo en parte el desequilibrio entre clases).

Métricas con la matriz de confusión (cont.)

• La tasa de verdaderos positivos o sensibilidad (recall) es el porcentaje de verdaderos positivos sobre el total de positivos observados \[ \scriptsize TVP=SENSIT=\frac{VP}{VP+FN} \]

  • ej., tasa de fraude/enfermos existentes que se detectan correctamente

• La tasa de verdaderos negativos o especificidad es el porcentaje de verdaderos negativos sobre el total de negativos observados \[ \scriptsize TVN=ESPECIF=\frac{VN}{VN+FP} \]

  • ej., tasa de “otras” opciones que se clasifican correctamente

  • Tasa de falsos positivos: \(\scriptsize TFP = 1 - TVN = 1 - ESPECIF\)

Métricas con la matriz de confusión (y 3)

• La exactitud equilibrada (Balanced Accuracy) es una media de la sensibilidad y de la especificidad

• La precisión o valor de predicción positivo es la cantidad de verdaderos positivos sobre el total de positivos predichos

\[ \scriptsize PREC=\frac{VP}{VP+FP} \]

• La familia de medidas \(\small F_{\beta}\) es una ratio de la importancia ponderada de la sensibilidad y de la precisión: \(\scriptsize F_{\beta}=\frac{(1+\beta)^2 \times SENSIT \times PREC}{\beta^2 \times SENSIT + PREC}\)

  • Para \(\scriptsize \beta<1\), se da menos importancia a la sensibilidad: los falsos positivos se consideran más costosos

  • Para \(\scriptsize \beta>1\), los falsos negativos son más costosos y para \(\scriptsize \beta=1\) son igualmente costosos

Curva ROC (“Receiver Operating Characteristic”)

• Representa TVP (eje y) frente a TFP (eje x) en diferentes umbrales : reducir el umbral clasifica más elementos como positivos (verdaderos y falsos)

• La curva ROC informa del grado de separabilidad: dado un nivel de TFP, el clasificador es mejor cuanto mayor sea TVP

AUC (“area under the curve”)

• La AUC es el área bajo la curva ROC: ofrece una medida agregada de rendimiento entre 0 (todas las clasificaciones incorrectas) y 1 (todas correctas)

• Resume la curva ROC y permite comparar curvas que se cruzan

Extensiones. Métricas adicionales

• Con más de dos clases, se realiza un análisis AUC-ROC para cada categoría (frente a las demás) y se promedian (ej., ponderando por casos en cada clase)

• Con clases desequilibradas, se puede preferir en lugar de la ROC un gráfico de precisión frente sensibilidad (precision-recall) y su correspondiente AUC (PR-AUC)

• Existen múltiples funciones de pérdida (o coste de clasificación) posibles.

  • Las relacionadas con la curva de ganancia consideran el coste de alcanzar un cierto nivel de sensibilidad

  • Otras se basan en la función de verosimilud o la entropía como medidas de pérdida (ej. mean log loss)

Evaluación de Modelos: entrenamiento y prueba

• Para minimizar problemas de underfit y, sobre todo, de overfit, DEBEMOS dividir aleatoriamente el conjunto de datos en dos partes:

• Entrenamiento (80-90%): datos sobre los que se construye/estima el modelo

• Prueba(20-10%): se usa el modelo construido para predecir y se evalúa con datos no vistos antes

• ¿Por qué renunciar a parte de los datos si sabemos que un tamaño muestral grande es importante? Evaluar correctamente un modelo lo es mucho más

• La estimación del error en prueba puede ser volátil dependiendo de las observaciones incluidas en cada grupo

Evaluación de Modelos: Validación cruzada

• Para evitar que los datos sean sensibles a una partición concreta, se usa validación cruzada (cross-validation o rotation estimation)

• Se repite varias veces y de forma ordenada el proceso de remuestreo para la partición en grupos de entrenamiento y prueba (similar a bootstrap)

• Permite utilizar todas las observaciones de la muestra, tanto para estimar como para evaluar el modelo (aunque no a la vez)

Validación cruzada de K bloques

• Se divide, aleatoriamente y ex-ante, la muestra en K subconjuntos (normalmente, 5 o 10)

• Un subconjunto se usa como prueba y el K-1 restantes como entrenamiento

• Se repite el proceso durante k iteraciones, con cada posible subconjunto de datos de prueba.

• Se obtiene una métrica de error en cada iteración; se promedian para obtener un único resultado de evaluación

• Es el tipo más habitual de validación cruzada

Validación cruzada aleatoria (RCV) y LOOCV

• RCV: en cada iteración se realiza la particion aleatoria (con reemplazamiento) entre entrenamiento y prueba

• Las observaciones pueden “repetir” como prueba

• LOOCV (leave one out CV): solo una observación se usa como prueba en cada iteración y el resto como entrenamiento

• Se realizan \(n\) iteraciones; se calcula una media sobre \(n\) resultados