Aprendizaje Supervisado I

Métodos de predicción lineal

Grado en Ciencia de Datos Aplicada • curso 2023-2024

Empieza lo bueno…empieza la modelización

Vamos a juntar las piezas del puzzle para hacer «magia»

¡Buenas!

Correo: . Despacho: 722 (3ª planta). Tutorías: lunes (14:30-16:00), martes (13:00-15:00) y viernes (13:00-14:00).

  • Javier Álvarez Liébana, de Carabanchel (Bajo).

  • Licenciado en Matemáticas (UCM). Doctorado en estadística (UGR).

  • Encargado de la visualización y análisis de datos covid del Principado de Asturias (2021-2022).

  • Miembro de la Sociedad Española de Estadística e IO y la Real Sociedad Matemática Española.

Actualmente, investigador y docente en la Facultad de Estadística de la UCM. Divulgando por Twitter e Instagram

Objetivos

  • Empezar a relacionar asignaturas como matemáticas, inferencia y R.

  • Aprender los fundamentos del aprendizaje estadístico (ahora llamado Machine Learning o data science)

  • Pasar de los descriptivo a lo predictivo: construir nuestros primeros modelos

  • Entender en profundidad el contexto de la predicción lineal.

Evaluación

  • Asistencia. Se valorará muy positivamente la participación. Si se restarán puntos si eres expulsado de clase: -0.3 la 1ª vez, -0.6 la 2ª, -1.2 la 3ª…
  • Evaluación continua: 1 examen a papel (25%), 2 entregas R en clase (20%-30%) y entrega grupal último día (3-6p obligatoriamente, 25%)
  • Examen final:
    • Menos de un 6 (sin la grupal) -> examen de R en lugar de la grupal.
    • Más de un 6.5 (con la grupal) -> podrás decidir peso del final entre un 20% y un 100%.
    • Entre 5 y 6.5 (con la grupal) -> podrás decidir peso del final entre un 40% y un 100%.
    • Por debajo de 5 (con la grupal) -> podrás decidir peso del final entre un 70% y un 100%.

Planificación

  • Entrega I (20%): semana 27-29 de febrero de 2024

  • Entrega II (30%): semana 2-4 de abril de 2024

  • Parcial con papel y boli (25%): 16-18 de abril de 2024

  • Entrega grupal (25%): 9 de mayo de 2024 (incluyendo lo visto hasta el 30 de abril)

 

(se admiten propuestas)

Clase a recuperar otro día: clase del 21 de marzo.

Materiales

En el menú de las diapositivas (abajo a la izquierda) tienes una opción para descargarlas en pdf en Tools (consejo: no lo hagas hasta el final del curso ya que irán modificándose)

 

Requisitos

Para el curso los únicos requisitos serán:

  1. Conexión a internet (para la descarga de algunos datos y paquetes).

  2. Instalar R y RStudio: la descarga la haremos (gratuitamente) desde https://cran.r-project.org/ y https://posit.co/download/rstudio-desktop/

  3. Se darán por asumido conocimientos aprendidos de R base, tidyverse y ggplot

  4. Se darán por asumido conocimientos aprendidos de Quarto, diapositivas en Quarto y Github. Para las entregas SOLO SE VALORARÁ la salida html correspodiente.

  5. Recomendable: saber usar la calculadora en modo estadístico.

Listado de clases

  • Clase 1: repaso de descriptiva
  • Clase 2: importancia de la visualización
  • Clase 3: intro al aprendizaje estadístico. Sesgo vs varianza. Supervisado vs no supervisado. Correlación vs dependencia vs causalidad.

Clase 1: repaso

Objetivo de la predicción lineal. Concepto de linealidad. Repaso de estadística descriptiva

¿Qué es predecir?

Como veremos más adelante, en esta asignatura vamos a tratar principalmente lo que se conoce el aprendizaje estadístico como predicción (continua)

Dada una variable objetivo (variable dependiente), y con la información aportada por un conjunto de variables predictoras (covariables), el objetivo será construir un modelo que consiga dar una estimación/predicción lo «mejor posible»

Es importante que - de momento - distingamos dos conceptos:

  • Estimación: el modelo aprende de unos datos e intenta estimar dichos valores que ha usado.
  • Predicción: el modelo aprende de unos datos e intenta estimar valores que el modelo no conoce.

Más adelante los llamaremos «predicción en train» y «predicción en test»

¿Qué es la linealidad?

En matemáticas decimos que una función \(f(x)\) es lineal cuando se cumple:

  • Propiedad aditiva: \(f(x + y) = f(x) + f(y)\)

  • Propiedad homogénea: \(f(k*x) = k*f(x)\) (donde \(k\) es una constante en \(\mathbb{R}\)).

Ambas se pueden resumir en \(f(a*x + b*y) = a*f(x) + b*f(y)\)

En estadística llamamos modelo de predicción lineal a un modelo que usa la información de covariables \(X_1, X_2, \ldots, X_p\), de manera que su información siempre se relacionen entre sí con sumas y restas.

  • Ejemplos lineales: \(y = 2*x_1 - 3\) o \(y = 4 - \frac{x_1}{2} + 3*x_2\)

  • Ejemplos no lineales: \(y = 2*\frac{1}{x_1}\) o \(y = 4 - x_{1}^{2} - x_2\) o \(y = ln(x_1) + cos(x_2)\)

En general son modelos que se pueden representar con rectas (una sola predictora), planos (dos predictoras), etc.

Repaso: continua vs discreta

  • Caracteres: cada una de las características o cualidades que se podrían medir o analizar para cada individuo de la población (y de los que disponemos el valor para cada individuo de la muestra).

  • Modalidades: diferentes valores que puede adoptar una característica o variable.

 

Como veremos más adelante, en el ámbito del aprendizaje estadístico va a ser fundamental tener clara la diferencia entre una variable cualitativa y cuantitativa y entre variable continua y otra discreta. ¿Cuál es la diferencia?

Repaso: continua vs discreta

Imagina las siguientes variables:

  • ¿Tienes hermanos? (si/no)
  • Color de zapatillas
  • Nivel de estudios
  • Número de hermanos
  • Número de pelos en la cabeza
  • Resultado de un dado
  • Temperatura ºC
  • Estatura o peso

 

¿Cuál es la diferencia entre ellas

Repaso: continua vs discreta

  • Cualitativas: cualidades o categorías. Ejemplos: sexo, estado civil, estudios, etc.
    • Ordinales: admiten relación jerarquica (suspenso-aprobado-bien, sano-leve-grave).
    • Nominales: no tienen asociada una jerarquía (sexo, religión, color de zapatillas, etc).
  • Cuantitativas: cuantificación numérica.
    • Discretas (barras): se pueden contar y enumerar (aunque sean infinitos) (granos de arena, nº hermanos, etc) → se suman
    • Continuas: toman infinitos valores y entre dos valores cualesquiera hay a su vez infinitas opciones (estatura, peso, etc) → se integran

Repaso: continua vs discreta

Repaso: población vs muestra

En estadística llamaremos población al universo teórico, al colectivo de individuos a estudiar, o de posibles elementos/eventos de los podríamos tener observaciones (ejemplo: 47 millones de españoles).

Cada uno de esos elementos o eventos se llaman individuos.

Aunque nuestro objetivo será conocer algunas de las propiedades de la población, la población suele ser inaccesible en su totalidad → SELECCIÓN de un conjunto de individuos

Repaso: población vs muestra

Para ello en estadística usamos lo que se conoce como muestra: un subconjunto finito, de tamaño \(n\), «representativo» de la población (en estudio estadístico realizado sobre la totalidad de una población se denomina censo).

Para ello existe una rama conocida como muestreo

  • Probabilístico: todos los individuos tienen oportunidad de ser seleccionados, y la probabilidad de que suceda puede ser modelizada (al azar, por estratos, etc).

  • No probabilístico: algunos elementos de la población no tienen posibilidad de selección (sesgo de exclusión), o su probabilidad no puede ser conocida.

🤔 ¿Sería adecuado hacer una encuesta sobre el streamer favorito de los jóvenes a través de una encuesta realizada por teléfono fijo?

Sesgos en el muestreo

Sesgo de selección: aparece cuando no se tiene en cuenta la forma en la que se han recogido los datos.

Sesgos en el muestreo

El ejemplo más famoso es el caso «Dewey defeats Truman» (Dewer derrota a Truman), el titular con el que abrió el Chicago Tribune en 1948, el mismo día en el que Truman ganó al repúblicano Dewer en las elecciones de 1948: sin esperar a los resultados, se basaron en una encuesta telefónica (sin contar con el sesgo que, en aquella época, solo la clase alta tenía teléfono).

Sesgos en el muestreo

¿Dónde reforzarías los aviones?

Sesgos en el muestreo

El sesgo del superviviente (un tipo de sesgo de selección) aparece cuando se toma una muestra de un fenómeno ignorando si los individuos elegidos tienen las mismas opciones respecto al mismo.

Repaso: medidas de centralización

  • Media: dada una muestra \(\boldsymbol{x} =\left(x_1, \ldots, x_n \right)\), la media muestral \(\overline{x}\) se define como la suma de todos los valores dividida por el tamaño muestral

\[\overline{x} = \frac{1}{n} \sum_{i=1}^{n} x_i\]

Geométricamente: es el valor «más cercano» de todos los datos a la vez (minimiza las distancias al cuadrado)

Media muestral

VENTAJAS

  • Fácil de calcular y entender
  • Fácil y eficiente de programar
  • Siempre existe (para cuantitativas)

DESVENTAJAS

  • No es un valor de los datos (la media de {1, 2, 3, 4} es 2.5)
  • Poco robusta (valores atípicos le afectan mucho)
  • Solo se puede definir para variables cuantitativas

Repaso: medidas de centralización

  • Mediana: dada una muestra \(\boldsymbol{x} =\left(x_1, \ldots, x_n \right)\), la mediana muestral se define como el valor que es mayor o igual que al menos el 50%, y menor igual que al menos el 50% de los datos

\[Me_{x} = \arg \min_{x_i} \left\lbrace F_i > 0.5 \right\rbrace, \quad Me_x = e_{i-1} + \frac{0.5 - F_{i-1}}{F_i - F_{i-1} }a_i\]

La mediana es el valor de en medio si ordenamos los datos (y si se pueden ordenar…)

Mediana muestral

VENTAJAS

  • Suele ser un valor de la muestra
  • Un poco más robusta que la media

DESVENTAJAS

  • Muy ineficiente (requiere un algoritmo de ordenación)
  • Solo definida para cuantitativas o cualitativas ordinales

Repaso: medidas de centralización

  • Moda: dada una muestra \(\boldsymbol{x} =\left(x_1, \ldots, x_n \right)\), la moda muestral se define como el valor o valores más repetidos (en caso de que existan).

\[Mo_x = \arg \max_{x_i} f_i, \quad Mo_x = e_{i-1} + \frac{d_i - d_{i-1}}{\left(d_i - d_{i-1} \right) + \left(d_i - d_{i+1} \right)}a_i\]

Gráficamente: representa el «pico» de un diagrama de barras o un histograma

Moda muestral

VENTAJAS

  • Es un valor de la muestra
  • Muy robusta
  • Se puede calcular para cualquier cuanti o cuali

DESVENTAJAS

  • No siempre existe (amodal) y pueden existir varias (bimodal, trimodal, etc)
  • Poco usada en inferencia

Repaso: medidas de centralización

¿Cuál es la mediana, la media y la moda?

Repaso: medidas de centralización

Repaso: medidas de dispersión

¿Qué tiene que ver la imagen con la dispersión?

Repaso: medidas de dispersión

El cambio climático no solo es porque aumente la temperatura media (centralización) sino por la aparición cada vez más frecuente de fenómenos extremos

Aumento de la variabilidad → aumento de la DISPERSIÓN

Repaso: medidas de dispersión

¿Cómo medir lo que se alejan los datos de la media?

Una primera idea podría ser medir la distancia de cada dato al centro, es decir, restar cada dato de la media, y después realizar su promedio.

Repaso: medidas de dispersión

Imagina que tenemos la siguiente muestra \(X = \left\lbrace -5, -3, -1, 0, 1, 3, 5 \right\rbrace\).

¿Cuánto vale la media?

La media vale 0 y la distancia a ella es…la propia muestra \(\left\lbrace -5, -3, -1, 0, 1, 3, 5 \right\rbrace\). ¿Cuál es el promedio de dichas distancias?

Pues…de nuevo vale 0.

Si la dispersión es 0…¿no hay dispersión? ¿No debería de dar 0 solo cuando los datos sean constantes?

Repaso: medidas de dispersión

Para evitar que se cancelen los signos lo que haremos será calcular el promedio PERO de las distancias al cuadrado, la conocida como varianza

\[s_{x}^{2} = \frac{1}{n} \sum_{i=1}^{n} \left(x_i - \overline{x} \right)^2 = \overline{x^2} - \overline{x}^2 \]

Cuidado

Tomar el valor absoluto (para evitar que se cancelen los signos) suele ser una mala idea en matemáticas (no es derivable como función).

Repaso: medidas de dispersión

Problema: si los datos están en metros, la varianza estará en…metros cuadrados

¿Tiene sentido medir la dispersión de nuestra estatura en baldosas?

Repaso: medidas de dispersión

Para tener una medida de dispersión en las unidades de los datos calcularemos la desviación típica, como la raíz cuadrada de la varianza

\[s_{x} = \sqrt{s_{x}^{2}} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} \left(x_i - \overline{x} \right)^2} = \sqrt{\overline{x^2} - \overline{x}^2}\]

Repaso: medidas de dispersión

Todavía tenemos un pequeño problema.

Imagina que queremos comparar la dispersión de dos conjuntos de datos, estaturas de personas y diámetros de núcleos de células. Y Supongamos que las medias son 170 cm y 5 micrómetros, y la desviación típica de 1 cm y 1.5 micrómetros.

¿Qué conjunto de datos es más disperso?

Para tener una medida de dispersión adimensional definiremos el coeficiente de variación

\[CV_{x} = \frac{s_{x}}{\left| \overline{x} \right|}\]

Repaso: medidas de localización

Las medidas de posición o localización nos localizan los datos, siendo valores que nos dividen un conjunto ordenado en subconjuntos del mismo tamaño (ejemplo: mediana es percentil 50).

  • Percentil: valores \(P_{\alpha}\) del conjunto ordenado que dejan por debajo, al menos, el \(\alpha\)% de datos y \(\alpha\)% por encima.

  • Decil: valores \(D_{\alpha}\) que dividen los datos en 10 partes iguales.

  • Cuartil: valores \(C_{\alpha}\) o \(q_{\alpha}\) que dividen los datos en 4 partes iguales.

Repaso: covarianza y correlación

¿Qué es en realidad la varianza?

La varianza es el promedio de las desviaciones al cuadrado (respecto a la media), apareciendo dos veces dicha desviación: puede ser entendida como una medida que cuantifica la relación de una variable CONSIGO MISMA

¿Y si quiésemos medir la relación de una variable X respecto a otra variable Y (en lugar de consigo misma)?

Repaso: covarianza y correlación

\[s_{x}^{2} = \frac{1}{n} \sum_{i=1}^{n} \left(x_i - \overline{x} \right)^2 = \overline{x^2} - \overline{x}^2 \quad \text{(varianza)}\]

La idea detrás de la covarianza es justo esa: sustituir una de esas desviaciones de la X por la desviación de la Y.

\[s_{xy} = \frac{1}{n} \sum_{i=1}^{n} \left(x_i - \overline{x} \right)\left(y_i - \overline{y} \right) = \overline{x*y} - \overline{x}*\overline{y}\]

Repaso: covarianza y correlación

Es importante entender algunas propiedades de la covarianza

  • Signo: la covarianza puede ser tanto positiva como negativa como 0: al eliminar el cuadrado de la varianza, ya no es necesario que sea positiva
  • ¿Qué cuantifica? La covarianza mide la relación LINEAL (en torno a una recta) entre dos variables
  • ¿Qué dice su signo? El signo de la covarianza nos indicará la dirección de la dependencia lineal: si es positiva, la relación será creciente (cuando X crece, Y crece); si es negativa, la relación será decreciente (cuando X crece, Y decrece)

Repaso: covarianza y correlación

Al igual que pasaba con la varianza, la covarianza depende de las unidades y magnitudes de los datos, así que lo que haremos será estandarizar la covarianza. Definiremos la coeficiente correlación lineal (de Pearson) como la covarianza dividida entre el producto de las desviaciones típicas (adimensional)

\[r_{xy} = \rho_{xy} = \frac{s_{xy}}{s_x s_y}\]

Tiene el mismo signo que la covarianza (el denominador es siempre positivo) y sus valores siempre están entre -1 y 1

  • más cerca de -1 o 1 → relación lineal más fuerte
  • más cerca de 0 → ausencia de relación LINEAL

Repaso: covarianza y correlación

Clase 2: importancia del dataviz

¿Basta con calcular la correlación para cuantificar si un ajuste es adecuado?

Conjunto anscombe

En el paquete {datasets} se encuentra el dataset conocido como cuarteto de Anscombe, un dataset que cuenta con 4 conjuntos de datos.

anscombe_tb <- as_tibble(datasets::anscombe)
anscombe_tb
# A tibble: 11 × 8
      x1    x2    x3    x4    y1    y2    y3    y4
   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1    10    10    10     8  8.04  9.14  7.46  6.58
 2     8     8     8     8  6.95  8.14  6.77  5.76
 3    13    13    13     8  7.58  8.74 12.7   7.71
 4     9     9     9     8  8.81  8.77  7.11  8.84
 5    11    11    11     8  8.33  9.26  7.81  8.47
 6    14    14    14     8  9.96  8.1   8.84  7.04
 7     6     6     6     8  7.24  6.13  6.08  5.25
 8     4     4     4    19  4.26  3.1   5.39 12.5 
 9    12    12    12     8 10.8   9.13  8.15  5.56
10     7     7     7     8  4.82  7.26  6.42  7.91
11     5     5     5     8  5.68  4.74  5.73  6.89

Primera tarea: convierte a tidydata

Conjunto anscombe

anscombe_x <- 
  anscombe_tb |>
  pivot_longer(cols = x1:x4, names_to = "dataset", values_to = "x") |> 
  select(-contains("y")) |> 
  mutate(dataset = str_remove_all(dataset, "x"))

anscombe_y <- 
  anscombe_tb |>
  pivot_longer(cols = y1:y4, names_to = "dataset", values_to = "y") |> 
  select(-contains("x"), -dataset)

anscombe_tidy <-
  anscombe_x |> mutate("y" = anscombe_y$y)
anscombe_tidy
# A tibble: 44 × 3
   dataset     x     y
   <chr>   <dbl> <dbl>
 1 1          10  8.04
 2 2          10  9.14
 3 3          10  7.46
 4 4           8  6.58
 5 1           8  6.95
 6 2           8  8.14
 7 3           8  6.77
 8 4           8  5.76
 9 1          13  7.58
10 2          13  8.74
# ℹ 34 more rows

Conjunto anscombe

¿Qué características muestrales tienen? Calcula la media, varianza, desv. típica, covarianza y correlación en cada dataset

anscombe_tidy |> 
  summarise(mean_x = mean(x), mean_y = mean(y),
            var_x = var(x), var_y = var(y),
            sd_x = sd(x), sd_y = sd(y),
            cov = cov(x, y), cor = cor(x, y), .by = dataset)
# A tibble: 4 × 9
  dataset mean_x mean_y var_x var_y  sd_x  sd_y   cov   cor
  <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1            9   7.50    11  4.13  3.32  2.03  5.50 0.816
2 2            9   7.50    11  4.13  3.32  2.03  5.5  0.816
3 3            9   7.5     11  4.12  3.32  2.03  5.50 0.816
4 4            9   7.50    11  4.12  3.32  2.03  5.50 0.817

Si te fijas todos los datasets tienen los mismos momentos muestrales, incluso tendrían el mismo ajuste de regresión…¿serán el mismo dataset desordenado?

Conjunto anscombe

Visualiza los 4 datasets

Código
ggplot(anscombe_tidy) +
  geom_point(aes(x = x, y = y, color = dataset), size = 3, alpha = 0.7) +
  ggthemes::scale_color_colorblind() +
  facet_wrap(~dataset) + theme_minimal()

Conjunto anscombe

Por suerte o por desgracia no todo son matemáticas: antes de pensar que modelo es mejor para nuestros datos, es importantísimo realizar un análisis exploratorio de los datos (incluyendo visualización)

Datasaurus

Podemos visualizarlo de manera aún más extrema con el dataset datasaurus_dozen del paquete {datasauRus} (ver más en https://www.research.autodesk.com/publications/same-stats-different-graphs/)

library(datasauRus)
datasaurus_dozen |>
  summarise(mean_x = mean(x), mean_y = mean(y),
            var_x = var(x), var_y = var(y),
            sd_x = sd(x), sd_y = sd(y),
            cov = cov(x, y), cor = cor(x, y), .by = dataset)
# A tibble: 13 × 9
   dataset    mean_x mean_y var_x var_y  sd_x  sd_y   cov     cor
   <chr>       <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
 1 dino         54.3   47.8  281.  726.  16.8  26.9 -29.1 -0.0645
 2 away         54.3   47.8  281.  726.  16.8  26.9 -29.0 -0.0641
 3 h_lines      54.3   47.8  281.  726.  16.8  26.9 -27.9 -0.0617
 4 v_lines      54.3   47.8  281.  726.  16.8  26.9 -31.4 -0.0694
 5 x_shape      54.3   47.8  281.  725.  16.8  26.9 -29.6 -0.0656
 6 star         54.3   47.8  281.  725.  16.8  26.9 -28.4 -0.0630
 7 high_lines   54.3   47.8  281.  726.  16.8  26.9 -30.9 -0.0685
 8 dots         54.3   47.8  281.  725.  16.8  26.9 -27.2 -0.0603
 9 circle       54.3   47.8  281.  725.  16.8  26.9 -30.8 -0.0683
10 bullseye     54.3   47.8  281.  726.  16.8  26.9 -31.0 -0.0686
11 slant_up     54.3   47.8  281.  726.  16.8  26.9 -31.0 -0.0686
12 slant_down   54.3   47.8  281.  726.  16.8  26.9 -31.2 -0.0690
13 wide_lines   54.3   47.8  281.  726.  16.8  26.9 -30.1 -0.0666

Datasaurus

Código
ggplot(datasaurus_dozen |> filter(dataset != "wide_lines"),
       aes(x = x, y = y, color = dataset)) +
  geom_point(size = 2.5, alpha = 0.75) +
  geom_smooth(method = "lm", se = FALSE)  +
  facet_wrap(~dataset, ncol = 3) +
  theme_minimal()

Clase 3: profundizando en la correlación

Matrices de correlación y covarianza. Correlación vs causalidad vs dependencia

Correlación lineal: sin agrupar

Como decíamos, la idea detrás de la covarianza es una “varianza” entre dos variales (la varianza es una covarianza de una variable consigo misma), midiendo el promedio de lo que se desvía cada una respecto a su media

\[s_{xy} = \frac{1}{n} \sum_{i=1}^{n} \left(x_i - \overline{x} \right)\left(y_i - \overline{y} \right) = \overline{x*y} - \overline{x}*\overline{y}\]

Para calcularla en R basta con usar la función cov()

starwars |> 
  drop_na(mass, height) |> 
  summarise(cov(mass, height))
# A tibble: 1 × 1
  `cov(mass, height)`
                <dbl>
1                806.

Correlación lineal: sin agrupar

Vamos a practicar una vez más como hacerlo a mano con el siguiente ejercicio.

Cor lineal: sin agrupar

En la tabla inferior se han recopilado (del 2013 al 2022) la temperatura media en el mes de abril en Madrid (variable X, en ºC) y el número de días (variable Y) en el que el nivel de ozono superó las 0.20 ppm (partes por millón)

  • ¿Cuál fue media de días en los que se superó umbral de ozono de 0.20 ppm?
  • ¿Cuál fue media de días en los que se superó umbral de ozono en los años que la temperatura media en marzo superó los 17.4ºC?
  • ¿Cuál es su covarianza?

Correlación lineal: sin agrupar

Repite el ejercicio con pocas líneas de código R

  • ¿Cuál fue la media de días en los que se superó el umbral de ozono de 0.20 ppm?
  • ¿Cuál fue la media de días en los que se superó el umbral de ozono en los años que la temperatura media en marzo superó los 17.4ºC?
  • ¿Cuál es su covarianza?

Correlación lineal: sin agrupar

Realiza lo que consideres tanto a mano como en R

  • ¿Existe alguna relación de dependencia entre las variables? ¿De qué tipo? ¿Cómo de fuerte o débil es dicha relación? ¿En qué dirección es dicha relación?

\[s_{xy} = \frac{1}{n} \sum_{i=1}^{n} \left(x_i - \overline{x} \right)\left(y_i - \overline{y} \right) = \overline{x*y} - \overline{x}*\overline{y}\]

\[r_{xy} = \rho_{xy} = \frac{s_{xy}}{s_x s_y}\]

Correlación lineal: sin agrupar

No sé si te has fijado qué sucede cuando intentamos calcular la covarianza/correlación de varias variables, por ejemplo vamos a calcular la (cuasi)covarianza de todas las variables numéricas de starwars.

starwars |> 
  select(where(is.numeric)) |> 
  drop_na() |>
  cov()
               height       mass birth_year
height       955.5071   715.3057  -2138.334
mass         715.3057 46313.2034  17402.547
birth_year -2138.3336 17402.5466  28603.045

Podemos usar la función cov() sin más, fuera de un resumen, obteniendo lo que se conoce como matriz de (cuasi)covarianzas y que tendrá un papel fundamental en estadística ya que contiene la información (= varianza) del dataset.

Matriz de covarianzas

starwars |> 
  select(where(is.numeric)) |> 
  drop_na() |>
  cov()
               height       mass birth_year
height       955.5071   715.3057  -2138.334
mass         715.3057 46313.2034  17402.547
birth_year -2138.3336 17402.5466  28603.045

Además de ser simétrica…¿qué tenemos en la diagonal?

La matriz de (cuasi)covarianzas se denota como \(\Sigma\) y sus elementos se define como \(\Sigma_{ii} = s_{x_i}^{2}\) para la diagonal y \(\Sigma_{ij} = \Sigma_{ji} = s_{x_i x_j}\) fuera de ella.

Importante

Recuerda que los softwares estadísticos nos devuelven siempre la cuasi covarianza, dividido entre \(n-1\) y no entre \(n\). La cuasivarianza y la cuasicovarianza son los mejores estimadores muestrales (insesgados) de los respectivos parámetros poblaciones

Matriz de correlaciones

De la misma manera con cor() podemos calcular la matriz de correlaciones (en este caso sin el cuasi ya que se cancelan denominadores)

starwars |> 
  select(where(is.numeric)) |> 
  drop_na() |>
  cor()
               height      mass birth_year
height      1.0000000 0.1075282 -0.4090274
mass        0.1075282 1.0000000  0.4781391
birth_year -0.4090274 0.4781391  1.0000000

La matriz de correlaciones se denota como \(R\) y sus elementos se define como \(r_{ii} = 1\) para la diagonal y \(r_{ij} = r_{x_ix_j}\) fuera de ella, y nos proporciona la dependencia lineal entre variables ya de manera estandarizada.

Matriz de correlaciones

¿Se te ocurre alguna manera de calcular la matriz de correlaciones a partir de la de covarianzas?

Correlación lineal: datos agrupados

Imagina ahora que queremos calcular la covarianza/correlación de los siguientes datos que representan 24 tiradas de dados

¿Cómo calcular la covarianza/correlación agrupando los datos?

Correlación lineal: datos agrupados

Correlación lineal: datos agrupados

Correlación vs dependencia

Podemos tener variables incorreladas, con correlación nula, pero que exista dependencia entre ellas: la covarianza/correlación SOLO CAPTURA relaciones lineales, nada más.

Veamos un ejemplo sencillo con \(X = \left\lbrace -1, 0, 1 \right\rbrace\) y \(Y = X^2 = \left\lbrace 1, 0, 1 \right\rbrace\).

  • La media de ambas es nula
  • La media del producto es la media de \(XY = \left\lbrace -1, 0, 1 \right\rbrace\), que es de nuevo nula
  • Así la covarianza \(\overline{x*y} - \overline{x}*\overline{y}\) es nula a pesar de tener la mayor dependencia posible (dependencia funcional)

Correlación vs dependencia

En relaciones no lineales como la de la imagen, la correlación estará cercana a cero (ya que no hay relación lineal) pero existe una dependencia. Diremos que dos variables son dependientes entre sí cuando existe un patrón numérico que las relaciona

  • Independencia implica incorrelación
  • Incorrelación NO implica independencia

Correlación vs dependencia

Dependencia vs causalidad

Si \(X\) e \(Y\) fuesen variables dependientes (lineales o no)…¿implicaría que el ascenso/descenso de una provoca el ascenso/descenso de la otra?

Imagina por ejemplo dos variables: nivel de bronceado (X) y consumo de helados (Y). ¿Son dependientes? Aparentemente sí ya que su comportamiento es similar. ¿Una causa la otra?

Dependencia vs causalidad

Diremos que dos variables tienen una relación causal o de causalidad cuando una variable es la CAUSA del comportamiento de la otra, algo que no sabremos solo con estadística (sino con conocimiento experto, en este caso de nutricionistas y médicos)

Correlación NO IMPLICA causalidad (al menos no tiene por qué)

Dependencia vs causalidad

Este fenómeno es conocido como correlaciones espúreas, y aparece cuando dos variables presentan una relación matemática pero sin ningún tipo de relación causal o lógica. Puedes ver más en https://www.tylervigen.com/spurious-correlations

Dependencia vs causalidad

Dicho patrón matemático puede deberse al mero azar o la existencia de lo que se conoce como variables de confusión (en el caso del helado, el verano es el que causa ambas).

La rama que se dedica al estudio de la causalidad, combinando las ramas de la estadística, filosofía, sociología y psicología se le conoce como inferencia causal y es fundamental en un análisis más profunde de las relaciones entre las variables (sobre todo en campos como la economía o la sociología)

Clase 4: primeras regresiones

Historia regresión. Aprendizaje supervisado. Regresión lineal univariante

Historia de la reg

Hay una regla universal: cualquier pariente tuyo es probablemente más mediocre que tú»

La historia de regresión se remonta a Francis Galton, primo de Charles Darwin, que además de estadístico fue psicólogo, geógrafo y, por desgracia, el primer eugenésico (de hecho acuñó el termino)

También fue el primero en proponer métodos de clasificación de huellas en medicina forense e incluso se le atribuye el primer mapa meteorológico de la historia

Regresión y Darwin

Galton mostró fascinación por «El origen de la especies» de su primo. Sin embargo, Galton no se centraba en los mejor adaptados sino en los que él llama mediocres

Según Galton, las sociedades estaban fomentando la mediocridad, interfiriendo en la selección natural, así que empezó a estudiar si el talento era o no hereditario.

¿Su conclusión? El talento se disipaba a lo largo de las generaciones

Regresión a la mediocridad

En 1886 publicó «Regression towards mediocrity in hereditary stature», un artículo que cambiaría la estadística: fue el primer uso documentado de lo que hoy conocemos como recta de regresión

Galton analizó la estatura de 205 hijos y padres, observando que, de nuevo, los valores extremos se disipaban: a lo largo de las generaciones había una regresión (un retroceso) a la mediocridad (entendida como la media)

Hijos de altos eran algo más bajitos, e hijos de bajitos eran algo más altos.

Regresión a la mediocridad

Galton no solo observó que las estaturas «regresaban» a un valor medio sino que lo hacían con un patrón, con un factor constante de \(2/3\): si los padres se desviaban \(+3\) por encima de la media, los hijos se desviaban solo \((2/3)*3 = +2\) por encima de la media.

Aprendizaje estadístico

La regresión lineal es el modelo más simple de lo que se conoce como aprendizaje estadístico (Machine Learning), en concreto del conocido como aprendizaje supervisado

Ciencia de Datos

La ciencia de datos es precisamente la rama que integra las matemáticas, la estadística, la probabilidad, el Machine Learning e incluso el Big Data

Aprendizaje ¿supervisado?

Aprendizaje ¿supervisado?

En el campo del Machine Learning hay principalmente dos tipos de metodologías:

  • Aprendizaje supervisado: tendremos dos tipos de variables, la variable dependiente (output/target) que se quiere predecir/clasificar, normalmente denotada como \(Y\), y las variables independientes (inputs) o explicativas o predictoras, que contienen la información disponible. Ejemplos: regresión, knn, árboles, etc.

Aprendizaje ¿supervisado?

En el campo del Machine Learning hay principalmente dos tipos de metodologías:

  • Aprendizaje no supervisado: no existe la distinción entre target y variables explicativas ya que no tenemos etiquetados los datos, no sabemos a priori la respuesta correcta. El aprendizaje no supervisado buscará patrones basados en similitudes/diferencias. Ejemplos: PCA, clustering, redes neuronales, etc.

Clasificación vs predicción

Como hemos comentado, la regresión lineal se enmarca dentro del predicción supervisada

  • Predicción: la variable objetivo es una variable cuantitativa continua (por ejemplo, precio, glucosa, peso, etc).

  • Clasificación: la variable objetivo es una variable cualitativa (por ejemplo, especie de flor, ausencia/presencia de enfermedad, si/no, etc) o cuantitativa discreta (por ejemplo, número de accidentes). La etiqueta tomará un valor dentro del conjunto de modalidades permitidas, pudiendo ser binaria (si/no) o multiclase (A, B, C, D).

 

📚 Ver «The elements of Statistical Learning» (Hastie et al., 2008)

Clase 5: ajuste de regresión

Interpretación de coeficientes. Método mínimos cuadrados

Modelo predictivo

Dentro del marco de la predicción supervisada un modelo tendrá siempre la siguiente forma:

\[Y = f(\mathbf{X}) + \varepsilon = f\left(X_1, \ldots, X_p \right) + \varepsilon, \quad E \left[Y | \boldsymbol{X} = x \right] =f\left(X_1, \ldots, X_p \right) \]

  • \(X\) serán los datos

  • \(f(\cdot)\) será nuestro modelo, es decir, el valor esperado de \(Y\) con la información que tenemos.

  • \(\mathbf{X} = \left(X_1, \ldots, X_p \right)\) serán nuestras predictoras o variables independientes

  • \(\varepsilon\) será el error o ruido, una variable aleatoria de media 0 \(E \left[\varepsilon | \boldsymbol{X} = x \right] = 0\) (el error debería ser reducido a algo aleatorio (irreducible), aunque en estadística SIEMPRE nos vamos a equivocar).

Modelo predictivo

\[Y = f(\mathbf{X}) + \varepsilon = f\left(X_1, \ldots, X_p \right) + \varepsilon, \quad E \left[Y | \boldsymbol{X} = x \right] =f\left(X_1, \ldots, X_p \right) \]

El objetivo es intentar estimar dicha variable objetivo minimizando al máximo el error cometido (la parte que podemos evitar).

El modelo estimado se definirá como

\[\widehat{Y} := \widehat{E \left[Y | \boldsymbol{X} = x \right]} = \widehat{f}\left(X_1, \ldots, X_p \right), \quad \widehat{\varepsilon} = \widehat{Y} - Y\]

  • \(\widehat{Y}\) serán las estimaciones, definidas como la estimación del valor esperado de \(Y\) con la información que tenemos.

  • \(\widehat{f}\) será el modelo estimado

  • \(\widehat{\varepsilon}\) los errores estimados (con media muestral igual a cero).

Reg. lineal univariante

En el caso concreto de la regresión lineal nuestro modelo será un hiperplano lineal (en el caso de una variable, una simple recta):

\[E \left[Y | \boldsymbol{X} = x \right] = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p\]

Su estimación será por tanto

\[\widehat{Y} := \widehat{E \left[Y | \boldsymbol{X} = x \right]} = \widehat{f}(\mathbf{X}) = \widehat{\beta_0} + \widehat{\beta_1} X_1 + \ldots + \widehat{\beta_p} X_p \]

  • El modelo \(\widehat{f}(\cdot)\) será una recta/hiperplano lineal

  • En el caso de la regresión lineal univariante tendremos por tanto \(E \left[Y | \boldsymbol{X} = x \right] = \beta_0 + \beta_1X\) (\(p = 1\))

  • El objetivo será obtener la estimación de los \(\widehat{\beta}\) tal que minimicemos el error

Varianza residual

Como ya hemos comentado, una manera de cuantificar lo que nuestro modelo se equivoca es lo que llamamos varianza residual o error cuadrático medio (varianza del error), que podemos estimar como sigue (ya que la media muestral de los errores será cero)

\[s_{r}^{2} = \frac{1}{n} \displaystyle \sum_{i=1}^{n} \left( \widehat{\varepsilon}_{i} - \overline{\widehat{\varepsilon}}_{i}\right)^2 = \frac{1}{n} \displaystyle \sum_{i=1}^{n} \widehat{\varepsilon}_{i}^{2} \] ¿Cómo quedaría la fórmula si desarrollo lo que vale los errores estimados \(\widehat{\varepsilon}\)?

\[s_{r}^{2} := s_{r}^{2}\left(\widehat{\beta}_0, \widehat{\beta}_1 \right)= \frac{1}{n} \sum_{i=1}^{n} \left(Y_i - \widehat{Y}_i \right)^2 = \frac{1}{n} \sum_{i=1}^{n} \left[Y_i - \left(\widehat{\beta}_0 + \widehat{\beta}_1X_i \right)\right]^2\]

Método mínimos cuadrados

El método de los mínimos cuadrados, originariamente planteado a la vez por Legendre y Gauss, consiste en establecer que los parámetros óptimos serán aquellos \(\left(\widehat{\beta}_0, \widehat{\beta}_1 \right)\) que minimicen dicha suma de cuadrados (que minimicen la varianza residual)

\[\widehat{\boldsymbol{\beta}} = \left(\widehat{\beta}_0, \widehat{\beta}_1 \right) = \arg \min_{\left(\beta_0, \beta_1 \right) \in \mathbb{R}^2} s_{r}^{2}\left(\beta_0, \beta_1 \right) = \arg \min_{\left(\beta_0, \beta_1 \right) \in \mathbb{R}^2} n*s_{r}^{2}\left(\beta_0, \beta_1 \right)\]

Vamos a desarrollar ese cuadrado…

Método mínimos cuadrados

\[s_{r}^{2}\left(\beta_0, \beta_1 \right) = \frac{1}{n} \sum_{i=1}^{n} \left[Y_i - \left(\beta_0 + \beta_1 X_i \right)\right]^2\] Vamos a desarrollar ese cuadrado…

\[\begin{eqnarray} s_{r}^{2}\left(\beta_0, \beta_1 \right) &=& \frac{1}{n} \sum_{i=1}^{n} \left[Y_{i}^{2} + \left(\beta_0 + \beta_1 X_i \right)^2 - 2Y_{i}\left(\beta_0 + \beta_1 X_i \right)\right] \nonumber \\ &=& \frac{1}{n} \sum_{i=1}^{n} \left[Y_{i}^{2} + \left(\beta_{0}^{2} + \beta_{1}^{2}X_{i}^{2} + 2\beta_{0}\beta_{1} X_i \right) - 2Y_{i}\beta_{0} - 2Y_{i}\beta_{1} X_i \right] \end{eqnarray}\]

Método mínimos cuadrados

\[\begin{eqnarray} s_{r}^{2}\left(\beta_0, \beta_1 \right) &=& \frac{1}{n} \sum_{i=1}^{n} \left[Y_{i}^{2} + \left(\beta_0 + \beta_1 X_i \right)^2 - 2Y_{i}\left(\beta_0 + \beta_1 X_i \right)\right] \nonumber \\ &=& \frac{1}{n} \sum_{i=1}^{n} \left[Y_{i}^{2} + \left(\beta_{0}^{2} + \beta_{1}^{2}X_{i}^{2} + 2\beta_{0}\beta_{1} X_i \right) - 2Y_{i}\beta_{0} - 2Y_{i}\beta_{1} X_i \right] \end{eqnarray}\]

¿Cómo encontrar el mínimo de una función?

Efectivamente, derivando.

\[\frac{\partial s_{r}^{2}}{\partial \beta_0} = \frac{1}{n} \sum_{i=1}^{n} \left[ 2\beta_{0} + 2 \beta_1 X_i - 2Y_i \right] = 2 \left( \beta_{0} + \beta_{1} \overline{x} - \overline{y} \right)\]

\[\frac{\partial s_{r}^{2}}{\partial \beta_1} = \frac{1}{n} \sum_{i=1}^{n} \left[ 2\beta_{1}X_{i}^{2} + 2 \beta_0 X_i - 2Y_iX_i \right] = 2 \left( \beta_{1} \overline{x^2}+\beta_{0}\overline{x} - \overline{xy} \right)\]

Método mínimos cuadrados

Si lo agrupamos en un sistema y lo igualamos a cero para encontrar el óptimo tenemos

\[\left\{\beta_{0} + \beta_{1} \overline{x} = \overline{y} \atop \overline{x} \beta_{0} + \overline{x^2} \beta_{1} = \overline{xy} \right.\]

¿Cómo resolver un sistema?

Haciendo la regla de Cramer tenemos que

\[\widehat{\beta}_1 = \frac{\overline{xy} - \overline{x}*\overline{y}}{\overline{x^2} - \overline{x}^2} = \frac{s_{xy}}{s_{x}^{2}}\]

\[\widehat{\beta}_0 = \frac{\overline{y}*\overline{x^2} - \overline{x}*\overline{xy}}{\overline{x^2} - \overline{x}^2} = \overline{y} + \frac{\overline{y}*\overline{x}^2- \overline{x}*\overline{xy}}{\overline{x^2} - \overline{x}^2} = \overline{y} - \widehat{\beta}_1 \overline{x}\]

Estimación reg. univariante

Los estimadores de los parámetros de una regresión lineal univariante serán por tanto \(\widehat{\beta}_1 = \frac{s_{xy}}{s_{x}^{2}}\) (pendiente de la recta) y \(\widehat{\beta}_0 = \overline{y} - \widehat{\beta}_1 \overline{x}\) (ordenada en el origen)

Es importante distinguir lo población de lo muestral

  • Los parámetros \(\left( \beta_0, \beta_1 \right)\) son desconocidos, los parámetros poblacionales

  • Los parámetros \(\left( \widehat{\beta}_0, \widehat{\beta}_1 \right)\) son estimados a partir de los datos, son variables aleatorias ya que han sido calculados en función de una muestra aleatoria \(\left\lbrace \left(x_i, y_i \right) \right\rbrace_{i=1}^{n}\) de la población $(X, Y ) $

Estimación reg. univariante

Los estimadores de los parámetros de una regresión lineal univariante serán por tanto \(\widehat{\beta}_1 = \frac{s_{xy}}{s_{x}^{2}}\) (pendiente de la recta) y \(\widehat{\beta}_0 = \overline{y} - \widehat{\beta}_1 \overline{x}\) (ordenada en el origen)

¿Cuál es su interpretación?

  • Ordenada en el origen: también llamado intercepto, y denotado como \(\beta_0\) su valor real, es el valor de \(Y\) cuando \(X=0\). Es decir, \(\widehat{\beta}_0\) se puede interpretar como la estimación \(\widehat{Y}\) cuando \(X = 0\) (cuando dicha estimación tenga sentido).
  • Pendiente: denotado com \(\beta_1\) su valor real, cuantifica el incremento de \(Y\) cuando \(X\) aumenta una unidad. Es decir, \(\widehat{\beta}_1\) se puede interpretar como la variación de la estimación \(\widehat{Y}\) cuando \(X\) tiene un incremento unitario.

Regresión en R

Para hacer un ajuste de regresión lineal univarainte en R solo necesitamos utilizar la función lm() (de linear model) con dos argumentos:

  • data = ...: los datos de los que queremos sacar el ajuste (sin ausentes de momento)
  • formula = ...: la fórmula del ajuste (y frente a x, expresado como y ~ x)
datos <- 
  starwars |>
  drop_na(mass, height) 
ajuste_lineal <- lm(data = datos, formula = mass ~ height)

Regresión en R

Una vez realizado el ajuste podemos resumirlo con summary()

ajuste_lineal |> summary()

Call:
lm(formula = mass ~ height, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
 -61.43  -30.03  -21.13  -17.73 1260.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.8103   111.1545  -0.124    0.902
height        0.6386     0.6261   1.020    0.312

Residual standard error: 169.4 on 57 degrees of freedom
Multiple R-squared:  0.01792,   Adjusted R-squared:  0.0006956 
F-statistic:  1.04 on 1 and 57 DF,  p-value: 0.312

¿Qué representa cada bloque de la salida?

Regresión en R

ajuste_lineal |> summary()

Call:
lm(formula = mass ~ height, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
 -61.43  -30.03  -21.13  -17.73 1260.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.8103   111.1545  -0.124    0.902
height        0.6386     0.6261   1.020    0.312

Residual standard error: 169.4 on 57 degrees of freedom
Multiple R-squared:  0.01792,   Adjusted R-squared:  0.0006956 
F-statistic:  1.04 on 1 and 57 DF,  p-value: 0.312
  • Call: ...: la orden que hemos ejecutado

  • Residuals: ...: un resumen en forma de cuartiles de los residuales, es decir, de los errores (fíjate que aunque la media de los errores siempre es 0, no tiene porque serlo la mediana, de hecho que no lo sea ya nos indica que normales no van a ser seguramente).

Regresión en R

ajuste_lineal |> summary()

Call:
lm(formula = mass ~ height, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
 -61.43  -30.03  -21.13  -17.73 1260.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.8103   111.1545  -0.124    0.902
height        0.6386     0.6261   1.020    0.312

Residual standard error: 169.4 on 57 degrees of freedom
Multiple R-squared:  0.01792,   Adjusted R-squared:  0.0006956 
F-statistic:  1.04 on 1 and 57 DF,  p-value: 0.312
  • Coefficients: ...: de momento solo nos interesa la columna Estimate que nos da las estimaciones de los parámetros. En la fila Intercept siempre irá \(\widehat{\beta}_0\), y el resto de filas tendrá el nombre de la variable predictora a la que multiplica el parámetro (en este caso la fila height corresponde a la estimación \(\widehat{\beta}_1\)).

Regresión en R

ajuste_lineal |> summary()

Call:
lm(formula = mass ~ height, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
 -61.43  -30.03  -21.13  -17.73 1260.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.8103   111.1545  -0.124    0.902
height        0.6386     0.6261   1.020    0.312

Residual standard error: 169.4 on 57 degrees of freedom
Multiple R-squared:  0.01792,   Adjusted R-squared:  0.0006956 
F-statistic:  1.04 on 1 and 57 DF,  p-value: 0.312

En nuestro caso particular, el ajuste de regresión viene dado por \(\widehat{Y} = \widehat{\beta}_0 + \widehat{\beta}_1 X = -11.4868 + 0.6240*X\), donde \(X\) es la estatura: por cada cm extra que mida un personaje, el modelo estima que el peso aumenta 0.6240 kg, y la estimación de peso para un persona que mida 0 cm es de -11.4868 kg (obviamente esta estimación mucho sentido no tiene)

Regresión en R

ajuste_lineal |> summary()

Call:
lm(formula = mass ~ height, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
 -61.43  -30.03  -21.13  -17.73 1260.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.8103   111.1545  -0.124    0.902
height        0.6386     0.6261   1.020    0.312

Residual standard error: 169.4 on 57 degrees of freedom
Multiple R-squared:  0.01792,   Adjusted R-squared:  0.0006956 
F-statistic:  1.04 on 1 and 57 DF,  p-value: 0.312
  • Multiple R-squared: ... (de momento lo demás lo ignoramos): es el conocido como \(R^2\) o bondad de ajuste. Hablaremos de él en las próximas clases, pero de momento, nos basta saber que es una métrica de calidad del ajuste que va de 0 (peor) a 1 (mejor).

Predicción en R

Una vez estimado el modelo (sus coeficientes) podemos usarlo, bien para estimar valores que el modelo conoce (de los que aprendió) o para predecir nuevos valores de los que el modelo no ha aprendido

Simplemente necesitamos la función predict(modelo, nuevos_datos), usando como argumentos el modelo entrenado y un nuevo dataset que contenga una/s columna/s llamadas igual que nuestras variables predictores con los valores de la predictora para los que queremos predecir

predicciones <- 
  tibble("height" = c(15, 81.3, 153.1, 189.3, 210.3),
         "mass_predict" = predict(ajuste_lineal, tibble("height" = height)))
predicciones
# A tibble: 5 × 2
  height mass_predict
   <dbl>        <dbl>
1   15          -4.23
2   81.3        38.1 
3  153.         84.0 
4  189.        107.  
5  210.        120.  

Rango de fiabilidad

Aunque hablaremos en profundidad de la bondad de ajuste, supongamos que tenemos un modelo extremadamente bueno, con una bondad de ajuste \(R^2 = 0.99\). ¿Sería fiable nuestro modelo anterior para predecir el peso de un personaje con \(X = 15\) cm?

predicciones
# A tibble: 5 × 2
  height mass_predict
   <dbl>        <dbl>
1   15          -4.23
2   81.3        38.1 
3  153.         84.0 
4  189.        107.  
5  210.        120.  

Si te fijas la ¡predicción del peso es negativa!: por muy bebé que sea, algo pesará. ¿Por qué sucede esto?

Rango de fiabilidad

Nuestro modelo solo puede aprender de los datos que conoce: solo podremos predecir para valores nuevos de las predictoras dentro del rango de las predictoras con las que entrenamos el modelo

Da igual lo bueno que sea el modelo: las predicciones fuera del rango de las predictoras de entrenamiento no serán fiables ya que el modelo no sabe lo que sucede fuera.

Clase 6: diagnosis

Hipótesis e inferencia del modelo

Repaso

  • Modelo lineal: descrito como \(Y = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p + \varepsilon\), buscamos estimar el valor esperado \(E \left[Y | \boldsymbol{X} = x \right] = \beta_0 + \beta_1 X_1 + \ldots + \beta_p x_p\), donde \(\boldsymbol{X}\) es la población y \(x\) la muestra.
  • Caso univariante: \(p = 1\) tal que \(E \left[Y | \boldsymbol{X} = x \right] = \beta_0 + \beta_1 x\).
  • Estimación: \(\widehat{Y}\) en función de \(\left(\widehat{\beta}_0, \widehat{\beta}_1 \right)\). De todas las rectas posibles, queremos la recta de regresión mínimos-cuadrados: minimiza la varianza residual (calculamos sus derivadas parciales, igualamos a cero y resolvemos el sistema). Esos estimadores serán \(\widehat{\beta}_1 = \frac{s_{xy}}{s_{x}^{2}}\) y \(\widehat{\beta}_0 = \overline{y} - \widehat{\beta}_1 \overline{x}\).
  • Predicción: podemos estimar o predecir valores con la recta \(\widehat{Y} = \widehat{\beta}_0 + \widehat{\beta}_0 X\), siempre y cuando \(X = x\) esté dentro del rango de valores de entrenamiento del modelo.

Diagnosis

Hasta ahora no hemos pedido nada a la muestra \(\left\lbrace \left( X_i, Y_i \right) \right\rbrace_{i=1}^{n}\), ¿para qué necesitaríamos hipótesis entonces?

La razón es que, hasta ahora, lo único que hemos podido realizar es una estimación puntual de los parámetros, pero dado que dichos estimadores serán variables aleatorias, necesitaremos realizar inferencia estadística sobre ellos (recuerda: los parámetros son simpleme estimaciones para esa muestra de la población, de forma que dada otra muestra, la recta será distinta).

Para poder cuantificar la variabilidad y precisión de nuestras estimaciones necesitaremos que los datos cumplan ciertas hipótesis probabilísticas: lo interesante no es la estimación puntual de los parámetros a partir de la muestra sino lo que podamos inferir de ellos a la población

Diagnosis

En el caso de la regresión lineal univariante pediremos 4 hipótesis

  1. Linealidad: el valor esperado de \(Y\) es \(E \left[Y | \boldsymbol{X} = x \right] = \beta_0 + \beta_1 x\)
  1. Homocedasticidad: el modelo no podrá explicar toda la información (a veces se equivocará por arriba y otras por abajo), pero necesitamos que la varianza del error sea finita y constante, tal que \(\sigma_{r}^2 = \sigma_{\varepsilon}^2 = {\rm Var} \left[\varepsilon | \boldsymbol{X} = x \right] = cte < \infty\), que no varíe según aumenta o disminuye la \(X\).
  1. Normalidad: pediremos que \(\varepsilon \sim N \left(0, \sigma_{r}^2 \right)\)
  1. Independencia: los errores \(\left\lbrace \varepsilon_i \right\rbrace_{i=1}^{n}\) deben ser independientes entre sí (el error en una observación no depende de otras). En particular, serán incorrelados

\[{\rm Cor}_{\varepsilon_i \varepsilon_j} = E \left[\varepsilon_i \varepsilon_j \right] - E \left[\varepsilon_i \right] E \left[\varepsilon_j \right] = E \left[\varepsilon_i \varepsilon_j \right] = 0, \quad i \neq j\]

Diagnosis

Las 4 hipótesis se pueden resumir de manera teórica en

\[Y | \boldsymbol{X} = x \sim N \left(\beta_0 + \beta_1x, \sigma_{\varepsilon}^2 \right)\]

Su versión muestral sería simplemente

\[y_i | \boldsymbol{X} = x_i \sim N \left(\beta_0 + \beta_1x_i, \sigma_{\varepsilon}^2 \right), \quad i=1,\ldots,n\]

Inferencia de los parámetros

Las hipótesis nos permiten decir (lo demostraremos más adelante) que los parámetros estimados siguen una distribución (condicionada) normal de media el parámetro a estimar y de varianza el conocido como standard error

\[\widehat{\beta}_j | \left( X=x \right) \sim N \left(\beta_j, SE \left(\widehat{\beta}_j\right)^2 \right), \quad j=0,1\]

Por las propiedades de la normal y el teorema central del límite

\[\frac{\widehat{\beta}_j - \beta_j}{SE \left(\widehat{\beta}_j\right)}| \left( X=x \right) \sim N \left(0, 1 \right), \quad j=0,1\]

Inferencia de los parámetros

\[\widehat{\beta}_j | \left( X=x \right) \sim N \left(\beta_j, SE \left(\widehat{\beta}_j\right)^2 \right), \quad j=0,1\]

Veremos más adelante porqué pero de momento…estos son los \(SE\)

\[SE \left(\widehat{\beta}_0\right)^2 = \frac{\sigma_{\varepsilon}^{2}}{n}\left[ 1+ \frac{\overline{x}^2}{s_{x}^2} \right], \quad SE \left(\widehat{\beta}_1\right)^2 = \frac{\sigma_{\varepsilon}^{2}}{n s_{x}^2}\]

¿Qué propiedades tienen estos estimadores?

Inferencia de los parámetros

\[\widehat{\beta}_j | \left( X=x \right) \sim N \left(\beta_j, SE \left(\widehat{\beta}_j\right)^2 \right), \quad j=0,1\]

\[SE \left(\widehat{\beta}_0\right)^2 = \frac{\sigma_{\varepsilon}^{2}}{n}\left[ 1+ \frac{\overline{x}^2}{s_{x}^2} \right], \quad SE \left(\widehat{\beta}_1\right)^2 = \frac{\sigma_{\varepsilon}^{2}}{n s_{x}^2}\]

  • Insesgados: ambos son estimadores insesgados ya que la esperanza de la estimación es el valor a estimar. \(E \left[ \widehat{\beta}_j \right] = \beta_j\)
  • Precisión vs tamaño muestral: cuando \(n \to \infty\), sus varianzas tienden a cero, es decir, \(SE \left(\widehat{\beta}_j\right)^2 \to^{n \to \infty} 0\) ya que \(n\) está dividiendo. Traducción: a más datos, mayor precision (menos varianza tendrán los estimadores si repetimos la toma de muestras)

Inferencia de los parámetros

\[SE \left(\widehat{\beta}_0\right)^2 = \frac{\sigma_{\varepsilon}^{2}}{n}\left[ 1+ \frac{\overline{x}^2}{s_{x}^2} \right], \quad SE \left(\widehat{\beta}_1\right)^2 = \frac{\sigma_{\varepsilon}^{2}}{n s_{x}^2}\]

  • Precisión vs var residual: cuanto más grande sea la varianza del error \(\sigma_{\varepsilon}^{2}\), \(SE \left(\widehat{\beta}_j\right)^2\) crecerá (es decir, más ruido implicará más imprecisión)
  • Precisión vs varianza de X: cuanto más grande sea la varianza de la predictora \(s_{x}^2\), \(SE \left(\widehat{\beta}_j\right)^2\) decrecerá, tal que \(SE \left(\widehat{\beta}_j\right)^2 \to^{s_{x}^2 \to \infty} 0\). Dicho de otra forma: cuánta más información (varianza) contenga nuestra tabla, mayor precisión.
  • Precisión vs media X: solo afecta a la estimación de \(\beta_0\), cuya precisión decrece cuando la media aumenta. Esto tiene sentido ya que cuánto más grande en media sean los datos, menos fiable será la predicción para \(X=0\).

Inferencia de los parámetros

\[\widehat{\beta}_j | \left( X=x \right) \sim N \left(\beta_j, SE \left(\widehat{\beta}_j\right)^2 \right), \quad j=0,1\]

\[SE \left(\widehat{\beta}_0\right)^2 = \frac{\sigma_{\varepsilon}^{2}}{n}\left[ 1+ \frac{\overline{x}^2}{s_{x}^2} \right], \quad SE \left(\widehat{\beta}_1\right)^2 = \frac{\sigma_{\varepsilon}^{2}}{n s_{x}^2}\]

Pero tenemos un problema: la varianza residual poblacional (varianza población del error) \(\sigma_{\varepsilon}^{2}\) es desconocida. Como suele pasar, la sustituiremos por un estimador (insesgado) de la misma (ya demostaremos…), donde \(p\) es el número de variables (lo que R llama std error en la salida)

\[\widehat{\sigma_{\varepsilon}^{2}} = \frac{n}{n-p-1} s_{r}^{2} = \frac{1}{n-p-1} \displaystyle \sum_{i=1}^{n} \widehat{\varepsilon}_{i}^2 =_{p=1} \frac{1}{n-2} \displaystyle \sum_{i=1}^{n} \widehat{\varepsilon}_{i}^2\]

Inferencia de los parámetros

\[\frac{\widehat{\beta}_j - \beta_j}{SE \left(\widehat{\beta}_j\right)} | \left( X=x \right) \sim N \left(0, 1 \right), \quad j=0,1\]

¿Qué suele suceder cuando en una normal desconocemos la varianza real y la sustituimos por su estimador insesgado?

Si sustituimos el error standard por su estimador (sustituyendo la varianza residual por su estimador insesgado), obtenemos que sigue una t-Student

\[\frac{\widehat{\beta}_j - \beta_j}{\widehat{SE} \left(\widehat{\beta}_j\right)} | \left( X=x \right) \sim t_{n-p-1} , \quad j=0,1\] donde \(\widehat{SE} \left(\widehat{\beta}_0\right)^2 = \frac{\widehat{\sigma_{\varepsilon}^{2}}}{n}\left[ 1+ \frac{\overline{x}^2}{s_{x}^2} \right]\) y \(\widehat{SE} \left(\widehat{\beta}_1\right)^2 = \frac{\widehat{\sigma_{\varepsilon}^{2}}}{n s_{x}^2}\)

Inferencia de los parámetros

ajuste_lineal |> summary()

Call:
lm(formula = mass ~ height, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
 -61.43  -30.03  -21.13  -17.73 1260.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.8103   111.1545  -0.124    0.902
height        0.6386     0.6261   1.020    0.312

Residual standard error: 169.4 on 57 degrees of freedom
Multiple R-squared:  0.01792,   Adjusted R-squared:  0.0006956 
F-statistic:  1.04 on 1 and 57 DF,  p-value: 0.312
  • std error: es lo que hemos denotado como \(\widehat{SE} \left(\widehat{\beta}_j\right)^2\)
  • t value: es el estadístico \(\frac{\widehat{\beta}_j - 0}{\widehat{SE} \left(\widehat{\beta}_j\right)}\) que hemos dicho que sigue una t-Student (fíjate que hemos puesto \(\beta_j = 0\))

Inferencia de los parámetros

ajuste_lineal |> summary()

Call:
lm(formula = mass ~ height, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
 -61.43  -30.03  -21.13  -17.73 1260.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.8103   111.1545  -0.124    0.902
height        0.6386     0.6261   1.020    0.312

Residual standard error: 169.4 on 57 degrees of freedom
Multiple R-squared:  0.01792,   Adjusted R-squared:  0.0006956 
F-statistic:  1.04 on 1 and 57 DF,  p-value: 0.312
  • Pr(>|t|): representa un p-valor. ¿De qué contraste se te ocurre?
  • Residual standard error: es el estimador insesgado de la varianza residual, que hemos denotado como \(\widehat{\sigma_{\varepsilon}^{2}}\)

Inferencia de los parámetros

\[\frac{\widehat{\beta}_j - \beta_j}{\widehat{SE} \left(\widehat{\beta}_j\right)} | \left( X=x \right) \sim t_{n-p-1} , \quad j=0,1\]

Si te acuerdas de inferencia, esto implica que los intervalos de confianza de nuestros parámetros son

\[IC \left(\alpha, \widehat{\beta}_j \right) = \left[\widehat{\beta}_j - t_{n-p-1; \alpha/2} \widehat{SE} \left(\widehat{\beta}_j\right), \widehat{\beta}_j + t_{n-p-1; \alpha/2} \widehat{SE} \left(\widehat{\beta}_j\right) \right]\] donde \(t_{n-p-1; \alpha/2}\) es el cuantil \(1-\alpha/2\) (deja a la izquierda una probabilidad de \(1-\alpha/2\)).

Inferencia de los parámetros

Podemos calcular los intervalos de confianza de los parámetros en R con confint()

confint(ajuste_lineal, level = 0.99)
                  0.5 %     99.5 %
(Intercept) -310.022729 282.402102
height        -1.029793   2.306935
confint(ajuste_lineal, level = 0.95)
                   2.5 %     97.5 %
(Intercept) -236.3934126 208.772785
height        -0.6150891   1.892231
confint(ajuste_lineal, level = 0.9)
                     5 %       95 %
(Intercept) -199.6638922 172.043265
height        -0.4082166   1.685359

Inferencia de los parámetros

¿Qué significa realmente un intervalo de confianza?

La probabilidad de que el parámetro real caiga dentro NO EXISTE: el parámetro es desconocido pero fijo, no aleatorio, así que no sentido calcular su probabilidad.

Si obtenemos para un parámetro un intervalo \([-1, 1]\) al 95%, no significa que \(P(parametro \in [-1, 1]) = 0.95\): significa que si tomamos 1000 muestras distintas de la población y calculamos para cada una el intervalo de confianza, aproximadamente 950 intervalos de confianza contendrán dentro el parámetro real

Un intervalo al 95% implica que habrá una frecuencia esperada de 0.95 de que intervalos que no conocemos (porque se derivan de muestras que no hemos tomado) contengan al parámetro real, pero no es la probabilidad de que tu intervalo calculado contenga a dicho parámetro: nos habla la precisión de nuestra metodología de estimación, no del parámetro.

Inferencia de los parámetros

Deberes. Dada una población normal \(\mu = 3\) y \(\sigma = 1.2\), crea un código que genere 500 muestras distintas (tamaño \(n = 100\) cada una), de manera que para cada una apliques un t.test() para calcular un IC para \(\mu\). Tras ello, gráfica los intervalos como se muestra en la imagen (haz uso de geom_segment())

Inferencia de los parámetros

Y si tenemos inferencia, tenemos contrastes: ¿te acuerdas de los p-valores que devuelve la tabla para cada parámetro?

Para cada parámetro se realiza un contraste de significancia: ¿cuánta evidencia hay en mis datos para poder decir que el valor estimado de mi parámetro es distinto de 0?

\[H_0:~\widehat{\beta}_j = 0 \quad vs \quad H_1:~\widehat{\beta}_j \neq 0\]

El estadístico usado es el t value que devuelve R definido como \(\frac{\widehat{\beta}_j - 0}{\widehat{SE} \left(\widehat{\beta}_j\right)}\) y que, BAJO EL SUPUESTO DE QUE LA HIPÓTESIS NULA SEA CIERTA, sigue una t-Student \(t_{n-p-1}\)

Inferencia de los parámetros

\[H_0:~\widehat{\beta}_j = 0 \quad vs \quad H_1:~\widehat{\beta}_j \neq 0\]

Este contraste significancia en realidad contesta a la siguiente pregunta

¿Tiene mi variable \(X_j\) un efecto lineal SIGNIFICATIVO en mi variable \(Y\)?

 

Esa pregunta puede ser respondida desde un punto de vista frecuentista haciendo uso de los p-valores: si \(p-value < \alpha\), se suele rechazar la nula (rechazamos que la variable no tenga efecto lineal significativo, es decir, aceptamos que sí lo tiene); en caso contrario no se suele rechazar (que no es lo mismo que aceptarla). Pero…

Paréntesis: p-valor

¿Qué representa un p-valor?

Hablaremos más adelante de ellos pero este es un pequeño resumen respecto al consenso recogido en 2016 por la ASA (American Statistical Association) en su trabajo https://www.tandfonline.com/doi/epdf/10.1080/00031305.2016.1154108?needAccess=true

  • Habla sobre los datos. Los p-valores nos pueden servir como indicadores de cómo de incompatible son los datos respecto a un modelo/hipotesis/explicación asumida: habla sobre los datos, no sobre la veracidad de la hipotesis nula per se o la probabilidad de que los datos hayan salido tan extremos por azar.

Paréntesis: p-valor

¿Qué representa un p-valor?

Hablaremos más adelante de ellos pero este es un pequeño resumen respecto al consenso recogido en 2016 por la ASA (American Statistical Association) en su trabajo https://www.tandfonline.com/doi/epdf/10.1080/00031305.2016.1154108?needAccess=true

  • Realidades complejas, decisiones complejas. Cuidado con reducir a decisiones binarias realidades complejas. La regla del p-valor es una herramienta más, pero no debe ser la única en la que nos basemos para decidir. Otros aspectos a considerar: calidad de las medidas, diseño del estudio, evidencia externa en la literatura respecto a la causalidad subyacente, validez de las hipótesis planteadas, etc

Paréntesis: p-valor

¿Qué representa un p-valor?

Hablaremos más adelante de ellos pero este es un pequeño resumen respecto al consenso recogido en 2016 por la ASA (American Statistical Association) en su trabajo https://www.tandfonline.com/doi/epdf/10.1080/00031305.2016.1154108?needAccess=true

  • No hagas cherry picking: muestra de manera transparentes que has probado, que ha salido y que no, y no te quedes solo con lo que sale bien o los p-valores que te convenga (p-hacking).

Paréntesis: p-valor

¿Qué representa un p-valor?

Hablaremos más adelante de ellos pero este es un pequeño resumen respecto al consenso recogido en 2016 por la ASA (American Statistical Association) en su trabajo https://www.tandfonline.com/doi/epdf/10.1080/00031305.2016.1154108?needAccess=true

  • Significancia estadística no implica significacia real (científica, humana, económica, etc). Obtener p-valores muy pequeños no implica un mayor efecto que otros p-valores no tan pequeños (y viceversa). Cualquier efecto, por pequeño que sea, puede derivar en p-valores pequeños si el tamaño uestral o la precisión de las medidas es suficiente alto, y al contrario (efectos evidentes pueden derivar en p-valores altos si \(n\) es pequeño)

Paréntesis: p-valor

¿Qué representa un p-valor?

Hablaremos más adelante de ellos pero este es un pequeño resumen respecto al consenso recogido en 2016 por la ASA (American Statistical Association) en su trabajo https://www.tandfonline.com/doi/epdf/10.1080/00031305.2016.1154108?needAccess=true

  • Alternativas: intervalos de confianza, de credibilidad, métodos bayesianos, false discovery rates.

 

📚 Ver más en https://www.tandfonline.com/doi/epdf/10.1080/00031305.2016.1154108?needAccess=true, https://anabelforte.com/2020/11/15/contraste/ y https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4877414/

Inferencia de los parámetros

ajuste_lineal |> summary()

Call:
lm(formula = mass ~ height, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
 -61.43  -30.03  -21.13  -17.73 1260.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.8103   111.1545  -0.124    0.902
height        0.6386     0.6261   1.020    0.312

Residual standard error: 169.4 on 57 degrees of freedom
Multiple R-squared:  0.01792,   Adjusted R-squared:  0.0006956 
F-statistic:  1.04 on 1 and 57 DF,  p-value: 0.312

Si te fijas ambos p-valores están por encima de \(\alpha = 0.05\) (umbral adoptado habitualmente) los que nos dice que no hay evidencias de los datos sean compatibles con afirmar que haya efectos lineales significativos

¿Y si probamos a quitar \(\beta_0\) (es decir, la respuesta está centrada)?

Inferencia de los parámetros

Para ello basta añadir un -1 antes de las variables predictoras.

lm(data = datos, formula = mass ~ -1 + height) |> summary()

Call:
lm(formula = mass ~ -1 + height, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
 -60.53  -29.45  -21.98  -18.50 1259.59 

Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
height   0.5623     0.1232   4.566 2.64e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 168 on 58 degrees of freedom
Multiple R-squared:  0.2644,    Adjusted R-squared:  0.2517 
F-statistic: 20.85 on 1 and 58 DF,  p-value: 2.638e-05

Fíjate que ahora, amén que hay un efecto lineal significativo de la altura y de tener \(R^2\) mayor, \(\beta_1\) tiene una precisión mayor. En este caso solo podía quitar uno (perderíamos \(X\)), pero veremos más adelante cómo decidir cuál quitar si tuviésemos varias variables.

Clases 7 y 8: caso práctico

¿De qué depende el precio del vino?

Caso práctico

Vamos a poner en práctica lo aprendido con el dataset wine.csv disponible en el campus.

datos <- read_csv(file = "./datos/wine.csv")
datos
# A tibble: 27 × 7
    Year Price WinterRain  AGST HarvestRain   Age FrancePop
   <dbl> <dbl>      <dbl> <dbl>       <dbl> <dbl>     <dbl>
 1  1952  7.50        600  17.1         160    31    43184.
 2  1953  8.04        690  16.7          80    30    43495.
 3  1955  7.69        502  17.2         130    28    44218.
 4  1957  6.98        420  16.1         110    26    45152.
 5  1958  6.78        582  16.4         187    25    45654.
 6  1959  8.08        485  17.5         187    24    46129.
 7  1960  6.52        763  16.4         290    23    46584.
 8  1961  8.49        830  17.3          38    22    47128.
 9  1962  7.39        697  16.3          52    21    48089.
10  1963  6.71        608  15.7         155    20    48799.
# ℹ 17 more rows

Caso práctico

El conjunto de datos está formado por 27 observaciones (cosechas de vino rojo Burdeos) y 7 variables

  • Year, Age: año de la cosecha y número de años en barrica.
  • Price: precio en 1990-1991 de cada cosecha (en escala log)
  • WinterRain: lluvia (en mm) que cayó ese año en invierno.
  • AGST: crecimiento medio de la temperatura (en grados Celsius) durante la temporada.
  • HarvestRain: lluvia (en mm) que cayó ese año durante la cosecha.
  • FrancePop: población (miles de habitantes) de Francia.

Ver más en detalles en https://doi.org/10.1080/09332480.1995.10542468

El objetivo: predecir el precio

Caso práctico

datos
# A tibble: 27 × 7
    Year Price WinterRain  AGST HarvestRain   Age FrancePop
   <dbl> <dbl>      <dbl> <dbl>       <dbl> <dbl>     <dbl>
 1  1952  7.50        600  17.1         160    31    43184.
 2  1953  8.04        690  16.7          80    30    43495.
 3  1955  7.69        502  17.2         130    28    44218.
 4  1957  6.98        420  16.1         110    26    45152.
 5  1958  6.78        582  16.4         187    25    45654.
 6  1959  8.08        485  17.5         187    24    46129.
 7  1960  6.52        763  16.4         290    23    46584.
 8  1961  8.49        830  17.3          38    22    47128.
 9  1962  7.39        697  16.3          52    21    48089.
10  1963  6.71        608  15.7         155    20    48799.
# ℹ 17 more rows

Para predecir el precio vamos a usar (de momento) una regresión lineal univariante, donde \(Y = precio\) y deberemos elegir la predictora \(X\) más apropiada.

Pasos a seguir

  1. Análisis exploratorio inicial:
  • ¿Las variables son numéricas (continuas)?
  • ¿Tienen problemas de rango (por ejemplo, pesos negativos)? ¿Tienen datos ausentes?
  • ¿Cómo se distribuyen las variables? Ideas: resúmen numérico, histogramas/densidades, boxplots, gráficos de violín, etc
  • ¿Hay datos atípicos?

Pasos a seguir

  1. Análisis de dependencia:
  • ¿Qué predictora está más correlacionada (linealmente) con la variable objetivo a predecir? ¿Existe otro tipo de dependencia (pendiente implementar en R)?
  • ¿Cómo se relacionan las predictoras entre sí? ¿Están correlacionadas? Ideas: matriz de correlaciones, diagramas de dispersión vs Y, corrplots, etc
  1. Formulación del modelo
  1. Fase de estimación:
  • ¿Cuánto valen los parámetros estimados? ¿Cómo queda el ajuste?
  • ¿Qué interpretación tienen?

Pasos a seguir

  1. Fase de diagnosis (paquetes {performance} y {olsrr}):
  • ¿Cumplen los datos las hipótesis parámetricas requeridas para poder hacer inferencia?
  • ¿Cómo modificar los datos para que se cumplan?
  • Análisis de residuales
  1. Fase de inferencia:
  • ¿Qué variabilidad tienen las estimaciones de nuestro parámetros?
  • ¿Las predictoras/intercepto tienen un efecto lineal significativo?
  • ¿Debemos re-entrenar el modelo sin alguno de ellos?

Pasos a seguir

  1. Fase de evaluación:
  • ¿Es significativo el modelo? ANOVA: análisis de la varianza
  • ¿Qué información de la predictora explica el modelo? Parámetros de bondad de ajuste (\(R^2\) por ejemplo)
  1. Fase de predicción

Análisis exploratorio inicial

  1. Análisis exploratorio inicial:
  • ¿Las variables son numéricas (continuas)?
  • ¿Tienen problemas de rango (por ejemplo, pesos negativos)? ¿Tienen datos ausentes?

 

Para responder a dichas preguntas lo más sencillos es hacer uso de la función skim() del paquete {skimr}

library(skimr)
datos |> skim()

En este caso no tenemos ausentes ni problemas de codificación

Análisis exploratorio inicial

  • ¿Cómo se distribuyen las variables?
  • ¿Hay datos atípicos?

 

Para ello podemos visualizar la distribución de cada variable haciendo uso de densidades, histogramas y/o boxplots (por ejemplo), pero si queremos ahorrar tiempo en visualizar convertimos nuestro dataset en tidydata

datos_tidy <-
  datos |> 
  pivot_longer(cols = everything(), names_to = "variable", values_to = "values")
datos_tidy
# A tibble: 189 × 2
   variable      values
   <chr>          <dbl>
 1 Year         1952   
 2 Price           7.50
 3 WinterRain    600   
 4 AGST           17.1 
 5 HarvestRain   160   
 6 Age            31   
 7 FrancePop   43184.  
 8 Year         1953   
 9 Price           8.04
10 WinterRain    690   
# ℹ 179 more rows

Análisis exploratorio inicial

Código
ggplot(datos_tidy, aes(x = values)) +
  geom_density(aes(color = variable, fill = variable),
               alpha = 0.75) +
  ggthemes::scale_color_colorblind() +
  ggthemes::scale_fill_colorblind() +
  facet_wrap(~variable, scales = "free") +
  theme_minimal()

Análisis exploratorio inicial

Código
ggplot(datos_tidy, aes(x = values)) +
  geom_histogram(aes(color = variable, fill = variable),
               alpha = 0.75, bins = 10) +
  ggthemes::scale_color_colorblind() +
  ggthemes::scale_fill_colorblind() +
  facet_wrap(~variable, scales = "free") +
  theme_minimal()

Análisis exploratorio inicial

Código
ggplot(datos_tidy, aes(x = values)) +
  geom_histogram(aes(color = variable, fill = variable),
               alpha = 0.75, bins = 10) +
  ggthemes::scale_color_colorblind() +
  ggthemes::scale_fill_colorblind() +
  facet_wrap(~variable, scales = "free") +
  theme_minimal()

Análisis exploratorio inicial

Código
ggplot(datos_tidy, aes(y = values, x = variable)) +
  geom_boxplot(aes(color = variable, fill = variable),
               alpha = 0.75) +
  ggthemes::scale_color_colorblind() +
  ggthemes::scale_fill_colorblind() +
  facet_wrap(~variable, scales = "free") +
  theme_minimal()

No hay valores atípicos (respecto a los percentiles)

Análisis exploratorio inicial

Podemos incluso gráficar todos los scatter plot sin transformar los datos con facet_matrix() del paquete {ggforce}

Código
library(ggforce)
ggplot(datos) +
  geom_point(aes(x = .panel_x, y = .panel_y)) +
  facet_matrix(vars(everything())) +
  ggthemes::scale_color_colorblind() +
  ggthemes::scale_fill_colorblind() +
  theme_minimal()

Análisis exploratorio inicial

Podemos también visualizar con un scatter plot todas las variables y además sus correlaciones, con el paquete {GGally} y la función `ggpairs(), sin necesidad de convertir a tidydata.

Código
library(GGally)
ggpairs(datos) +
  theme_minimal()

Análisis de dependencia

  1. Análisis de dependencia:
  • ¿Qué predictora está más correlacionada (linealmente) con la variable objetivo a predecir? ¿Existe otro tipo de dependencia (pendiente implementar en R)?
  • ¿Cómo se relacionan las predictoras entre sí? ¿Están correlacionadas? Ideas: matriz de correlaciones, diagramas de dispersión vs Y, corrplots, etc

Este paso será crucial en el contexto multivariante pero en este caso simplemente vamos a ver como se relacionan linealmente las predictoras entre sí, y cuál de ellas es la más adecuada para predecir linealmente precio

Análisis de dependencia

El primer paso es la matriz de correlaciones con la función cor() o con la función correlate() del paquete {corrr} (importa en tibble más visual)

library(corrr)
datos |> correlate()
# A tibble: 7 × 8
  term           Year  Price WinterRain    AGST HarvestRain     Age FrancePop
  <chr>         <dbl>  <dbl>      <dbl>   <dbl>       <dbl>   <dbl>     <dbl>
1 Year        NA      -0.460     0.0512 -0.295      -0.0588 -1         0.992 
2 Price       -0.460  NA         0.135   0.668      -0.507   0.460    -0.481 
3 WinterRain   0.0512  0.135    NA      -0.321      -0.268  -0.0512    0.0295
4 AGST        -0.295   0.668    -0.321  NA          -0.0271  0.295    -0.301 
5 HarvestRain -0.0588 -0.507    -0.268  -0.0271     NA       0.0588   -0.0320
6 Age         -1       0.460    -0.0512  0.295       0.0588 NA        -0.992 
7 FrancePop    0.992  -0.481     0.0295 -0.301      -0.0320 -0.992    NA     
  • Respecto a Y: predictoras con mayor cor lineal son AGST (más calor, menos cosechas, sube el precio) y HarvestRain (más lluvias, más cosechas, baja el precio, ¡el signo importa!)

  • Dependencia entre predictoras: las variables Age, Year y FrancePop presentan la misma información.

Análisis de dependencia

También podemos usar corrplot() del paquete {corrplot}, al que le pasamos una matriz de correlaciones clásica y nos la visualiza.

library(corrplot)
datos |> cor() |> corrplot(method = "ellipse")

Puedes ver distintas opciones de visualización en https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html

Análisis de dependencia

Otra opción es visualizar con un scatter plots todas las predictoras vs Y, pivotando antes nuestro dataset (solo pivotamos las predictoras)

Código
datos_tidy <-
  datos |> pivot_longer(cols = -Price, names_to = "variable",
                        values_to = "values")
ggplot(datos_tidy, aes(x = values, y = Price)) + 
  geom_point(aes(color = variable), alpha = 0.7) +
  geom_smooth(method = "lm") +
  ggthemes::scale_color_colorblind() +
  facet_wrap(~variable, scales = "free_x") +
  theme_minimal()

No solo comprobamos que las rectas con más pendiente son AGST y HarvestRain, además los puntos parecen poder ajustarse a una recta sin otro patrón identificable.

Esto es importante hacerlo ya que debemos descartar posibles correlaciones espúreas (ver ejemplo datasaurus)

Formulación del modelo

Una vez que hemos decidido que dos predictoras usaemos, vamos por tanto a plantear dos modelos univariantes

\[Price = \beta_0 + \beta_1*AGST + \varepsilon\] \[Price = \beta_0 + \beta_1*HarvestRain + \varepsilon\]

Fase de estimación

  1. Fase de estimación:
  • ¿Cuánto valen los parámetros estimados? ¿Cómo queda el ajuste?
  • ¿Qué interpretación tienen?

 

Para ello ejecutaremos ambos modelos con lm()

ajuste_AGST <- lm(data = datos, formula = Price ~ AGST)
ajuste_harvest <- lm(data = datos, formula = Price ~ HarvestRain)

Fase de estimación

Ajuste con AGST

ajuste_AGST |> summary()

Call:
lm(formula = Price ~ AGST, data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.78370 -0.23827 -0.03421  0.29973  0.90198 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.5469     2.3641  -1.500 0.146052    
AGST          0.6426     0.1434   4.483 0.000143 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4819 on 25 degrees of freedom
Multiple R-squared:  0.4456,    Adjusted R-squared:  0.4234 
F-statistic: 20.09 on 1 and 25 DF,  p-value: 0.0001425
  • \(\beta_0=\) -3.547: predicción del precio cuando \(AGST = 0\) es de -3 (recuerda que está en escala logartímica)

  • \(\beta_1=\) 0.643: por cada grado de aumento, el precio sube 0.643 (escala log).

Fase de estimación

Ajuste con AGST


Call:
lm(formula = Price ~ AGST, data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.78370 -0.23827 -0.03421  0.29973  0.90198 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.5469     2.3641  -1.500 0.146052    
AGST          0.6426     0.1434   4.483 0.000143 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4819 on 25 degrees of freedom
Multiple R-squared:  0.4456,    Adjusted R-squared:  0.4234 
F-statistic: 20.09 on 1 and 25 DF,  p-value: 0.0001425
  • Residuales: además de media cero, parecen presentar una distribución simétrica con la mediana en torno al cero. Además se tiene que \(\widehat{\sigma}_{\varepsilon}^{2}=\frac{1}{n-2}\sum_{i=1}^{n} \widehat{\varepsilon}_{i}^{2} = 0.4819\) (estimador insesgado de la varianza residual) y \(R^2 = 0.4456\) (bondad de ajuste)

Fase de estimación

Ajuste con harvestRain

ajuste_harvest |> summary()

Call:
lm(formula = Price ~ HarvestRain, data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.03792 -0.27679 -0.07892  0.40434  1.21958 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.679856   0.241911  31.747  < 2e-16 ***
HarvestRain -0.004405   0.001497  -2.942  0.00693 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5577 on 25 degrees of freedom
Multiple R-squared:  0.2572,    Adjusted R-squared:  0.2275 
F-statistic: 8.658 on 1 and 25 DF,  p-value: 0.00693
  • \(\beta_0=\) 7.68: predicción del precio cuando la lluvia fue nula es de 7.679 (recuerda que está en escala logartímica)

  • \(\beta_1=\) -0.004: por cada litro de lluvia, precio baja 642.614 (escala log).

Fase de estimación

Ajuste con harvestRain


Call:
lm(formula = Price ~ HarvestRain, data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.03792 -0.27679 -0.07892  0.40434  1.21958 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.679856   0.241911  31.747  < 2e-16 ***
HarvestRain -0.004405   0.001497  -2.942  0.00693 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5577 on 25 degrees of freedom
Multiple R-squared:  0.2572,    Adjusted R-squared:  0.2275 
F-statistic: 8.658 on 1 and 25 DF,  p-value: 0.00693
  • Residuales: además de media cero, parecen presentar una distribución simétrica con la mediana en torno al cero. Además se tiene que \(\widehat{\sigma}_{\varepsilon}^{2}=\frac{1}{n-2}\sum_{i=1}^{n} \widehat{\varepsilon}_{i}^{2} = 0.5577\) (algo más grande que el otro ajuste) y \(R^2 = 0.2572\) (algo más pequeño que el otro ajuste) -> de momento es mejor el primer modelo.

Fase de diagnosis

  1. Fase de diagnosis (paquetes {performance} y {olsrr}):
  • ¿Cumplen los datos las hipótesis parámetricas requeridas para poder hacer inferencia?
  • ¿Cómo modificar los datos para que se cumplan?
  • Análisis de residuales

 

Recuerda que necesitamos verificar antes las hipótesis para poder hacer inferencia con los parámetros, así que vamos a ello con el paquete {performance} y el paquete {olsrr}

Fase de diagnosis

library(performance)
check_model(ajuste_AGST)

Diagnosis: linealidad

  1. Linealidad: el valor esperado de \(Y\) es \(E \left[Y | \boldsymbol{X} = x \right] = \beta_0 + \beta_1 x\)
Código
check_model(ajuste_AGST)

Si te fijas el gráfico que se refiere a ello está visualizando residuales vs valores estimados: está volviendo a plantear un segundo modelo de regresión donde ahora \(\widehat{\varepsilon}_i = \gamma_0 + \gamma_1 \widehat{y}_i\)

Diagnosis: linealidad

linealidad <- lm(data = tibble("fitted" = ajuste_AGST$fitted.values,
                               "residuals" = ajuste_AGST$residuals),
                 formula = residuals ~ fitted)
linealidad |> summary()

Call:
lm(formula = residuals ~ fitted, data = tibble(fitted = ajuste_AGST$fitted.values, 
    residuals = ajuste_AGST$residuals))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.78370 -0.23827 -0.03421  0.29973  0.90198 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.148e-15  1.574e+00       0        1
fitted      -3.084e-16  2.231e-01       0        1

Residual standard error: 0.4819 on 25 degrees of freedom
Multiple R-squared:  3.694e-32, Adjusted R-squared:  -0.04 
F-statistic: 9.234e-31 on 1 and 25 DF,  p-value: 1

Si te fijas ambos parámetros no son significativamente distintos de 0: no presentan una tendencia (lineal al menos)

Diagnosis: linealidad

Más adelante probaremos alguna otra cosa pero de momento nos basta con eso. También podemos visualizar nosotros ese scatter plot residuales vs estimaciones

Código
ggplot(tibble("fitted" = ajuste_AGST$fitted.values,
              "residuals" = ajuste_AGST$residuals),
       aes(x = fitted, y = residuals)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm") +
  theme_minimal()

Diagnosis: homocedasticidad

  1. Homocedasticidad: necesitamos que la varianza del error sea finita y constante, tal que \(\sigma_{r}^2 = \sigma_{\varepsilon}^2 = {\rm Var} \left[\varepsilon | \boldsymbol{X} = x \right] = cte < \infty\).
check_heteroscedasticity(ajuste_AGST)
OK: Error variance appears to be homoscedastic (p = 0.137).

El gráfico titulado Homogeneity of variance nos visualiza la raíz cuadrada del valor absoluto de los residuos estandarizados frente a las predicciones (se conoce como gráfico de escala-localización)

Diagnosis: homocedasticidad

Si visualizamos los residuales deberían estar en torno a 0, dentro de una banda constante (varianza constante)

Código
ggplot(tibble("id" = 1:length(ajuste_AGST$residuals),
              "residuals" = ajuste_AGST$residuals),
       aes(x = id, y = residuals)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal()

Diagnosis: homocedasticidad

Si visualizamos el gráfico de escala-localización deberíamos obtener un diagrama de dispersión cuya recta de regresión saliese casi plana en torno al 1.

Código
ggplot(tibble("fitted" = ajuste_AGST$fitted.values,
              "sqrt_std_residuals" = sqrt(abs((ajuste_AGST$residuals - mean(ajuste_AGST$residuals)) / sd(ajuste_AGST$residuals)))),
       aes(x = fitted, y = sqrt_std_residuals)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal()

Según el gráfico no deberíamos asumir homocedasticidad. ¿Por qué el contraste no la rechaza?

Diagnosis: homocedasticidad

Código
ggplot(tibble("fitted" = ajuste_AGST$fitted.values,
              "sqrt_std_residuals" = sqrt(abs((ajuste_AGST$residuals - mean(ajuste_AGST$residuals)) / sd(ajuste_AGST$residuals)))),
       aes(x = fitted, y = sqrt_std_residuals)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal()

Con el poco tamaño muestral que tenemos, es complicado tener evidencias que refuten la hipótesis nula (y el gráfico puede estar parcialmente diseñado). Por eso es la hipótesis más difícil de cumplir. Lo importante es que en la recta de regresión al dibujar los residuos no se aprecia una banda cuya anchura se modifique groseramente, más o menos constante

Diagnosis: normalidad

  1. Normalidad: pediremos que \(\varepsilon \sim N \left(0, \sigma_{r}^2 \right)\)

Con la función ols_test_normality() del paquete {olsrr} podemos obtener diferentes contrastes de normalidad

library(olsrr)
ols_test_normality(ajuste_AGST)
-----------------------------------------------
       Test             Statistic       pvalue  
-----------------------------------------------
Shapiro-Wilk              0.9685         0.5633 
Kolmogorov-Smirnov        0.0904         0.9657 
Cramer-von Mises          3.3055         0.0000 
Anderson-Darling          0.2422         0.7457 
-----------------------------------------------

Nos centraremos en los contrastes de Shapiro-Wilk, Kolmogorov-Smirnov y Anderson-Darling: no se rechaza normalidad

Diagnosis: normalidad

Además del contraste podemos visualizar con stat_qq() y stat_qq_line() el conocido como Q-Q plot: enfrenta los cuantiles de una muestra con los cuantiles de una normal teórica, teniendo que obtener los puntos en torno a una recta (especilamente en el centro).

Código
ggplot(tibble("residuals" = ajuste_AGST$residuals)) +
  stat_qq(aes(sample = residuals)) +
  stat_qq_line(aes(sample = residuals))

Diagnosis: independencia

  1. Independencia: los errores \(\left\lbrace \varepsilon_i \right\rbrace_{i=1}^{n}\) deben ser independientes entre sí (el error en una observación no depende de otras). En particular, serán incorrelados

\[{\rm Cor}_{\varepsilon_i \varepsilon_j} = E \left[\varepsilon_i \varepsilon_j \right] - E \left[\varepsilon_i \right] E \left[\varepsilon_j \right] = E \left[\varepsilon_i \varepsilon_j \right] = 0, \quad i \neq j\]

check_autocorrelation(ajuste_AGST)
OK: Residuals appear to be independent and not autocorrelated (p = 0.598).

Por último, check_autocorrelation() comprueba como efectivamente los residuales/errores son independientes, haciendo un test de autocorrelación (nos tiene que salir lo contrario a una serie temporal, que el error i no depende del i-1).

Diagnosis: independencia

Otra forma de verlo es visualizando los residuos respecto a su versión con retardo (por ejemplo, \(\left(\widehat{\varepsilon}_1, \widehat{\varepsilon}_2, \widehat{\varepsilon}_3, \ldots, \widehat{\varepsilon}_{n-1} \right)\) vs \(\left(\widehat{\varepsilon}_2, \widehat{\varepsilon}_3, \ldots, \widehat{\varepsilon}_n \right)\)

Código
ggplot(tibble("lag1" = ajuste_AGST$residuals[-length(ajuste_AGST$residuals)],
              "residuals" = ajuste_AGST$residuals[-1]),
       aes(x = residuals, y = lag1)) +
  geom_point(size = 3, alpha = 0.7) + 
  geom_smooth(method = "lm") + 
  theme_minimal()

Fase de diagnosis

check_model(ajuste_AGST)

En nuestro caso se cumplen todas las hipótesis (algunas más fuertemente que otras).

Repite el proceso con el otro modelo

Fase de diagnosis

check_predictions(ajuste_AGST)

Nos faltan dos gráficas por comentar:

  • Posterior Predictive Checks: simula distintas variables respuesta suponiendo que el modelo fuese cierto (añadiendo ruido aleatorio) y lo compara con la muestra. Si lo observado se distancia mucho de las simulaciones es que el modelo planteado no ajusta bien a la muestra.

Fase de diagnosis

check_outliers(ajuste_AGST)
OK: No outliers detected.
- Based on the following method and threshold: cook (0.713).
- For variable: (Whole model)
  • Influential Observations: nos permite identificar observaciones influyentes, marcando aquellas (con su id de fila) que se salgan fuera de la banda definida por la conocida como distancia de Cook denotada como \(D_i\) (realiza, para cada observación, la suma de todos los cambios de la regresión cuando la observación \(i\) es retirada: si hay muchos cambios al cambiar una observación, es que era muy influyente)

Diferencia dos tipos:

  • outliers: valor atípico de la respuesta pudiendo perturbar la varianza residual
  • high-leverage points: valor atípico en alguna de las predictoras

Fase de inferencia

  1. Fase de inferencia:
  • ¿Qué variabilidad tienen las estimaciones de nuestro parámetros?
  • ¿Las predictoras/intercepto tienen un efecto lineal significativo?
  • ¿Debemos re-entrenar el modelo sin alguno de ellos?

 

Una vez verificadas las hipótesis lo que haremos será inferir conclusiones de la población en función de la muestra

Fase de inferencia

ajuste_AGST |> summary()

Call:
lm(formula = Price ~ AGST, data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.78370 -0.23827 -0.03421  0.29973  0.90198 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.5469     2.3641  -1.500 0.146052    
AGST          0.6426     0.1434   4.483 0.000143 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4819 on 25 degrees of freedom
Multiple R-squared:  0.4456,    Adjusted R-squared:  0.4234 
F-statistic: 20.09 on 1 and 25 DF,  p-value: 0.0001425
  • Variabilidad de las estimaciones de nuestro parámetros

    • \(\widehat{SE} \left( \widehat{\beta}_0 \right)\) igual a 2.344 por lo que (aprox) \(\widehat{\beta_0} \sim N(-3.547, 2.344)\)
    • \(\widehat{SE} \left( \widehat{\beta}_1 \right)\) igual a 0.143 por lo que (aprox) que \(\widehat{\beta_1} \sim N(0.643, 0.143)\)

Fase de inferencia

ajuste_AGST |> summary()

Call:
lm(formula = Price ~ AGST, data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.78370 -0.23827 -0.03421  0.29973  0.90198 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.5469     2.3641  -1.500 0.146052    
AGST          0.6426     0.1434   4.483 0.000143 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4819 on 25 degrees of freedom
Multiple R-squared:  0.4456,    Adjusted R-squared:  0.4234 
F-statistic: 20.09 on 1 and 25 DF,  p-value: 0.0001425
  • Estadístico del contraste

    • \(\frac{\widehat{\beta}_0 - 0}{\widehat{SE} \left( \widehat{\beta}_0 \right)}\) igual a -1.5 (valor que tendrías que buscar en las tablas a mano)
    • \(\frac{\widehat{\beta}_1 - 0}{\widehat{SE} \left( \widehat{\beta}_1 \right)}\) igual a 4.483 (valor que tendrías que buscar en las tablas a mano)

Fase de inferencia

ajuste_AGST |> summary()

Call:
lm(formula = Price ~ AGST, data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.78370 -0.23827 -0.03421  0.29973  0.90198 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.5469     2.3641  -1.500 0.146052    
AGST          0.6426     0.1434   4.483 0.000143 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4819 on 25 degrees of freedom
Multiple R-squared:  0.4456,    Adjusted R-squared:  0.4234 
F-statistic: 20.09 on 1 and 25 DF,  p-value: 0.0001425
  • Efecto (lineal): si nos fijamos en la tabla, el p-valor de \(\beta_0\) es 0.146052. Si adoptamos \(\alpha = 0.05\) como suele ser habitual, el contraste \(H_0:~\beta_0 = 0\) vs \(H_1:~\beta_0 \neq 0\) nos dice que no podemos rechazar de forma significativa la hipótesis nula (no sucede con \(\beta_1\), si sucediese no habría modelo)

Fase de inferencia

ajuste_AGST |> summary()

Call:
lm(formula = Price ~ AGST, data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.78370 -0.23827 -0.03421  0.29973  0.90198 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.5469     2.3641  -1.500 0.146052    
AGST          0.6426     0.1434   4.483 0.000143 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4819 on 25 degrees of freedom
Multiple R-squared:  0.4456,    Adjusted R-squared:  0.4234 
F-statistic: 20.09 on 1 and 25 DF,  p-value: 0.0001425

¿Y si quitamos dicho parámetro?

Para quitarlo añadimos un -1 al modelo

ajuste_AGST_sin_beta0 <- lm(data = datos, formula = Price ~ -1 + AGST)

Re-aprendiendo

ajuste_AGST_sin_beta0 |> summary()

Call:
lm(formula = Price ~ -1 + AGST, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.7324 -0.3756 -0.0092  0.3838  1.0804 

Coefficients:
     Estimate Std. Error t value Pr(>|t|)    
AGST 0.427691   0.005757   74.29   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4933 on 26 degrees of freedom
Multiple R-squared:  0.9953,    Adjusted R-squared:  0.9951 
F-statistic:  5519 on 1 and 26 DF,  p-value: < 2.2e-16
  • La bondad de ajuste ha pasado de \(R^2 = 0.446\) a \(R^2 = 0.9953\)

  • La variabilidad de la estimación \(\widehat{SE} \left( \widehat{\beta}_1 \right)\) ha pasado de 0.143 a 0.005757.

Comparar modelos

Aunque no hemos hablado en profundidad de las métricas de evaluación podemos comparar los modelos con compare_performance() del paquete {performance}

compare_performance(ajuste_AGST, ajuste_AGST_sin_beta0,
                    ajuste_harvest, ajuste_harvest_sin_beta0)
# Comparison of Model Performance Indices

Name                     | Model | AIC (weights) | AICc (weights) | BIC (weights) |    R2 | R2 (adj.) |  RMSE | Sigma
---------------------------------------------------------------------------------------------------------------------
ajuste_AGST              |    lm |  41.1 (0.535) |   42.2 (0.469) |  45.0 (0.379) | 0.446 |     0.423 | 0.464 | 0.482
ajuste_AGST_sin_beta0    |    lm |  41.4 (0.454) |   41.9 (0.522) |  44.0 (0.614) | 0.995 |     0.995 | 0.484 | 0.493
ajuste_harvest           |    lm |  49.0 (0.010) |   50.1 (0.009) |  52.9 (0.007) | 0.257 |     0.228 | 0.537 | 0.558
ajuste_harvest_sin_beta0 |    lm | 147.5 (<.001) |  148.0 (<.001) | 150.1 (<.001) | 0.762 |     0.753 | 3.450 | 3.515

Clases 9: evaluación y predicción

¿Cómo se puede evaluar un modelo? ¿Qué métricas existen? ¿Cómo predecir?

Repitamos el proceso

Para interiorizar lo aprendido, vamos a repetir todo el proceso con el conjunto datos_linearreg.csv (las variables predictoras representan diferentes variables meteorológicas y la variable objetivo y la temperatura media en primavera, para distintas ciudades).

datos <- read_csv(file = "./datos/datos_linearreg.csv")
datos
# A tibble: 1,000 × 4
      x1     y     x2       x3
   <dbl> <dbl>  <dbl>    <dbl>
 1  5.64  17.3 -2.43   2.19   
 2  4.00  14.5  0.881  6.06   
 3  5.02  15.7  0.457  6.91   
 4  5.86  18.4 -4.03   1.27   
 5  4.65  13.3 -2.37   3.68   
 6  4.97  14.7 -5.48  -0.00953
 7  4.05  13.0 -1.24   3.79   
 8  5.38  17.8 -3.36   2.21   
 9  5.16  17.2 -3.52   1.42   
10  4.43  13.3  1.42   6.84   
# ℹ 990 more rows

Regresión lineal

Debes seguir los siguientes pasos de la manera más detallada posible

  1. Análisis exploratorio inicial tanto numérico como visualizando. ¿Son numéricas sin problemas de codificación? ¿Cómo se distribuyen? ¿Hay datos atípicos?

  2. Análisis de dependencia. ¿Qué predictora correlaciona más con la objetivo? ¿Cómo se relacionan las predictoras entre sí?

  3. Formulación del modelo

  4. Fase de estimación. ¿Cuáles son los parámetros? ¿Cómo se interpretan?

  5. Fase de diagnosis

  6. Fase de inferencia. ¿Qué variabilidad tiene la estimación? ¿Hay efecto significativo?

Fase de evaluación

ajuste_lineal <- lm(data = datos, formula = y ~ x1)
ajuste_lineal |> summary()

Call:
lm(formula = y ~ x1, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2967 -1.0477  0.0599  1.0080  4.8827 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.77574    0.23885  -3.248   0.0012 ** 
x1           3.10868    0.04687  66.324   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.498 on 998 degrees of freedom
Multiple R-squared:  0.8151,    Adjusted R-squared:  0.8149 
F-statistic:  4399 on 1 and 998 DF,  p-value: < 2.2e-16

Teníamos pendiente la fase final: fase de evaluación

Fase de evaluación

  1. Fase de evaluación:
  • ¿Es significativo el modelo? ANOVA: análisis de la varianza
  • ¿Qué información de la predictora explica el modelo? Parámetros de bondad de ajuste (\(R^2\) por ejemplo)
  • ¿Qué otras métricas o herramientas podemos usar para cuantificar la calidad predictora de nuestro ajuste

Una de las herramientas más útiles para evaluar nuestro modelo es enfrentar los valores ajustados con los valores reales (dado que los conocemos al ser aprendizaje supervisado)

Fase de evaluación

Código
ggplot(tibble("y" = datos$y, "y_est" = ajuste_lineal$fitted.values),
       aes(x = y, y = y_est)) +
  geom_point(size = 1.2, alpha = 0.75) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal() +
  labs(x = "Valores reales", y = "Valores estimados")

En el gráfica podemos ver como los valores reales vs estimados están muy cercanos a la diagonal: el error cometido es muy pequeño.

Fase de evaluación

Podemos considerar algunas métricas para cuantificar el acierto del modelo

  • SSE (sum of squared errors): definido como la suma de errores al cuadrado

\[SSE = \sum_{i=1}^{n} \widehat{\varepsilon}_{i}^2 = \sum_{i=1}^{n} \left( Y_i - \widehat{Y}_i \right)^2\] Fíjate que dado que \(\widehat{\sigma}_{\varepsilon} = \widehat{\sigma}_{r} = \frac{n}{n-p-1} s_{r}^{2} = \frac{1}{n-p-1} \sum_{i=1}^{n} \widehat{\varepsilon}_{i}^2\) (el que R llama Residual standard error), tenemos que \(SSE = (n-p-1)\widehat{\sigma}_{\varepsilon}\)

  • MSE: media de lo anterior \(MSE = s_{r}^{2} = \frac{1}{n} SSE = \frac{1}{n} \sum_{i=1}^{n} \widehat{\varepsilon}_{i}^2 = \frac{n-p-1}{n}\widehat{\sigma}_{\varepsilon}\)

Ambas métricas nos hablan de la varianza del error, es decir, de lo que el modelo no es capaz de explicar

Fase de evaluación

  • SSR (regressions sum of squares): definido como la suma de las desviaciones de cada predicción a su media (al tener estimadores insesgados, la media de las estimaciones \(\overline{\widehat{Y}}\) es la misma que la de la variable a estimar \(\overline{Y}\))

\[SSR = \sum_{i=1}^{n} \left( \overline{Y} - \widehat{Y}_i \right)^2 = \sum_{i=1}^{n} \left( \overline{\widehat{Y}}_i - \widehat{Y}_i \right)^2\]

  • MSR: media de lo anterior \(MSR = s_{\widehat{y}}^2 = \frac{1}{n} SSR = \frac{1}{n} \sum_{i=1}^{n} \left( \overline{Y} - \widehat{Y}_i \right)^2\)

Ambas métricas nos hablan de la variación de Y en torno a la regresión, es decir, la variación de \(\overline{Y}\) que es explicada por la media condicional estimada \((Y_i|X=x_i) \sim \widehat{\beta}_0 + \widehat{\beta}_1 X_i\), cuantifica la información de Y explicada por el modelo

Fase de evaluación

  • SST (total sum of squares): definido como la suma de las desviaciones de la variable objetivo \(Y\) a su media.

\[SST = \sum_{i=1}^{n} \left( Y_i - \overline{Y} \right)^2 \]

  • Varianza muestral de Y: la media de lo anterior (la varianza de \(Y\))

\[s_{y}^2= \frac{1}{n} SST = \frac{1}{n}\sum_{i=1}^{n} \left( Y_i - \overline{Y} \right)^2\]

Ambas métricas nos hablan de la variación total de Y, es decir, la información total de nuestra variable objetivo

ANOVA

Así tenemos 3 tipos de métricas:

  • SST y \(s_{y}^2\): el total de info a explicar

  • SSR y \(MSR\): el total de info explicada por el modelo

  • SSE y \(MSE\): el total de info NO explicada por el modelo (a veces se usa la raíz cuadrada del MSE, conocido como \(RMSE\), o el \(MAE\), tomando valor absoluto en los errores).

Se pueden demostrar que, SOLO EN EL CASO LINEAL, desarrollando el sumatorio \(SST = \sum_{i=1}^{n} \left( Y_i - \overline{Y} \right)^2 = \sum_{i=1}^{n} \left( \left(Y_i - \widehat{Y}_i \right) + \left( \widehat{Y}_i - \overline{Y} \right) \right)^2\) se llega a que

SST = SSR + SSE

ANOVA

SST = SSR + SSE

\(s_{y^2}\) = \(s_{\widehat{y}}^2\) + \(s_{r}^2\) (equivalente, \(s_{y^2}\) = \(s_{\widehat{y}}^2\) + \(\frac{n-p-1}{n}\widehat{\sigma}_{\varepsilon}\))

Esta decomposición (solo se cumple en el caso lineal) se conoce como ANOVA o análisis de la varianza y podemos hacerlo en R con aov() o anova()

ajuste_lineal |> aov()
Call:
   aov(formula = ajuste_lineal)

Terms:
                      x1 Residuals
Sum of Squares  9870.948  2239.494
Deg. of Freedom        1       998

Residual standard error: 1.497993
Estimated effects may be unbalanced
ajuste_lineal |> anova()
Analysis of Variance Table

Response: y
           Df Sum Sq Mean Sq F value    Pr(>F)    
x1          1 9870.9  9870.9  4398.9 < 2.2e-16 ***
Residuals 998 2239.5     2.2                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA

ajuste_lineal |> aov()
Call:
   aov(formula = ajuste_lineal)

Terms:
                      x1 Residuals
Sum of Squares  9870.948  2239.494
Deg. of Freedom        1       998

Residual standard error: 1.497993
Estimated effects may be unbalanced
Terms x1 (predictora) Residuals
Sum of Squares SSR SSE
Deg. of Freedom (grados libertad) p n - p - 1

 

Residual standard error: \(\widehat{\sigma}_{\varepsilon}^{2}= \frac{n}{n-p-1} s_{r}^{2} = \frac{1}{n-p-1} \sum_{i=1}^{n} \widehat{\varepsilon}_{i}^{2}\)

ANOVA

ajuste_lineal |> anova()
Analysis of Variance Table

Response: y
           Df Sum Sq Mean Sq F value    Pr(>F)    
x1          1 9870.9  9870.9  4398.9 < 2.2e-16 ***
Residuals 998 2239.5     2.2                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Df (grados) Sum Sq Mean Sq F value Pr(>F)
x1 p SSR \(\frac{SSR}{p}\) \(F-value = \frac{\frac{SSR}{p}}{\frac{SSE}{n-p-1}}\) p-valor F test
Residuals n - p - 1 SSE \(\widehat{\sigma}_{\varepsilon}^{2} = \frac{SSE}{n-p-1}\)

El F-value es el estadítico \(F = \frac{\frac{SSR}{p}}{\frac{SSE}{n-p-1}} \sim F_{p, n-p-1}\) asociado al contraste de significación global

\[H_0:~\beta_1 = \ldots = \beta_p = 0 \quad vs \quad H_1:~\text{existe al menos un} \quad \beta_j \neq 0\]

ANOVA

ajuste_lineal |> anova()
Analysis of Variance Table

Response: y
           Df Sum Sq Mean Sq F value    Pr(>F)    
x1          1 9870.9  9870.9  4398.9 < 2.2e-16 ***
Residuals 998 2239.5     2.2                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\[H_0:~\beta_1 = \ldots = \beta_p = 0 \quad vs \quad H_1:~\text{existe al menos un} \quad \beta_j \neq 0\]

El contraste pretende responder a: ¿existe una dependencia lineal entre \(Y\) y el CONJUNTO de predictoras? (global, no parámetro a parámetro).

Si se rechaza significa que existe al menos un predictor cuyo efecto LINEAL sobre Y es significativo.

Importante: en el caso de la reg. lineal univariante, \(F-value\) y \(p-value\) del ANOVA es equivalente al \(t-value\) y \(p-value\) del contraste de significación para \(\beta_1\) (ya que…no hay más).

Bondad de ajuste

El conocido como \(R^2\) o bondad de ajuste o el coeficiente de determinación es una forma sencilla, y relacionada con el ANOVA, para cuantificar la cantidad de información explicada por el modelo

De hecho es, literal, un ratio de información explicada (lo que R llama Multiple R-squared, ya hablaremos de su versión ajustada en regresión multivariante)

\[R^2 = \text{Ratio info explicada} = \frac{SSR}{SST}\]

En el caso lineal, por lo visto en el ANOVA

\[R^2= \frac{SSR}{SST} = \frac{SST - SSE}{SST} = 1 - \frac{SSE}{SST} = 1 - \frac{s_{r}^{2}}{s_{y}^{2}} = 1 - \frac{\frac{n-p-1}{n}\widehat{\sigma}_{\varepsilon}^{2}}{s_{y}^{2}}\]

Por definición, la bondad de ajuste está entre 0 y 1.

Bondad de ajuste

Deberes: demuestra que \(R^2 = r_{y \widehat{y}}^2\) y que en el caso de \(p=1\) coincide además con \(R^2 = r_{x,y}^2\)

\[R^2 = r_{y \widehat{y}}^2 =_{p=1} r_{x,y}^2\]

La bondad de ajuste tiene un problema importante: no solo depende del modelo sino que también de los datos. ¿De qué depende? ¿Por qué es peligroso usar ciegamente \(R^2\) para decidir?

Bondad de ajuste

\[R^2= \frac{SSR}{SST} = \frac{SST - SSE}{SST} = 1 - \frac{SSE}{SST} = 1 - \frac{s_{r}^{2}}{s_{y}^{2}} = 1 - \frac{\frac{n-p-1}{n}\widehat{\sigma}_{\varepsilon}^{2}}{s_{y}^{2}}\]

¿De qué depende?

  • Más predictoras implica que \(R^2\) crece, ¡incluso aunque dichas predictoras no sean significativas! (esto lo arreglaremos con el \(R^2\) ajustado en el futuro)

  • Más ruido, menor \(R^2\)

  • Ignora si el modelo cumple las hipótesis: un modelo con un alto \(R^2\) puede dar predicciones nefastas.

Predictoras vs R2

Vamos a repetir pero añadiendo 5 y 22 variables más, sin efecto lineal con \(Y\).

datos_extras <-
  tibble("y" = datos$y, "x1" = datos$x1, "x2" = rnorm(1e3), "x3" = rnorm(1e3),
         "x4" = rnorm(1e3), "x5" = rnorm(1e3), "x6" = rnorm(1e3), "x7" = rnorm(1e3),
         "x8" = rnorm(1e3), "x9" = rnorm(1e3), "x10" = rnorm(1e3), "x11" = rnorm(1e3),
         "x12" = rnorm(1e3), "x13" = rnorm(1e3), "x14" = rnorm(1e3), "x15" = rnorm(1e3),
         "x16" = rnorm(1e3), "x17" = rnorm(1e3), "x18" = rnorm(1e3), "x19" = rnorm(1e3),
         "x20" = rnorm(1e3), "x21" = rnorm(1e3), "x22" = rnorm(1e3), "x23" = rnorm(1e3))
ajuste_6pred <- lm(data = datos_extras, formula = y ~ x1 + x2 + x3 + x4 + x5 + x6)
ajuste_23pred <- lm(data = datos_extras, formula = y ~ .)
compare_performance(ajuste_lineal, ajuste_6pred, ajuste_23pred)
# Comparison of Model Performance Indices

Name          | Model |  AIC (weights) | AICc (weights) |  BIC (weights) |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------------------------------------------------
ajuste_lineal |    lm | 3650.1 (0.906) | 3650.2 (0.912) | 3664.9 (>.999) | 0.815 |     0.815 | 1.496 | 1.498
ajuste_6pred  |    lm | 3654.7 (0.093) | 3654.8 (0.088) | 3693.9 (<.001) | 0.816 |     0.815 | 1.492 | 1.498
ajuste_23pred |    lm | 3663.5 (0.001) | 3664.8 (<.001) | 3786.2 (<.001) | 0.821 |     0.816 | 1.474 | 1.492

Con compare_performance() podemos comparar métricas de modelos.

Ruido vs R2

Hemos dicho que el ruido afecta…¿qué crees que pasará si fijamos la parte determinística y solo modificamos el ruido? Vamos a simular 6 modelos, con exactamente los mismos \(\beta_0\) y \(\beta_1\) (es decir, mismo ajuste) pero con diferente varianza en el ruido (supongamos que \(X \sim N(3, 1.5)\)).

\[\text{Modelo 1: } Y = -1.2 + 3.2X + N(0, 0.25)\] \[\text{Modelo 2: } Y = -1.2 + 3.2X + N(0, 1)\]

\[\text{Modelo 3: } Y = -1.2 + 3.2X + N(0, 1.5)\]

\[\text{Modelo 4: } Y = -1.2 + 3.2X + N(0, 2)\]

\[\text{Modelo 5: } Y = -1.2 + 3.2X + N(0, 4)\]

\[\text{Modelo 6: } Y = -1.2 + 3.2X + N(0, 8)\]

Ruido vs R2

Código
x <- rnorm(n = 1000, mean = 3, sd = 1.5)
eps1 <- rnorm(n = 1000, mean = 0, sd = 0.25)
eps2 <- rnorm(n = 1000, mean = 0, sd = 1)
eps3 <- rnorm(n = 1000, mean = 0, sd = 1.5)
eps4 <- rnorm(n = 1000, mean = 0, sd = 2)
eps5 <- rnorm(n = 1000, mean = 0, sd = 4)
eps6 <- rnorm(n = 1000, mean = 0, sd = 8)
datos <- tibble("x" = x, "y1" = -1.2 + 3.2*x + eps1,
                "y2" = -1.2 + 3.2*x + eps2, "y3" = -1.2 + 3.2*x + eps3,
                "y4" = -1.2 + 3.2*x + eps4, "y5" = -1.2 + 3.2*x + eps5,
                "y6" = -1.2 + 3.2*x + eps6)
ajuste_1 <- lm(data = datos, formula = y1 ~ x)
ajuste_2 <- lm(data = datos, formula = y2 ~ x)
ajuste_3 <- lm(data = datos, formula = y3 ~ x)
ajuste_4 <- lm(data = datos, formula = y4 ~ x)
ajuste_5 <- lm(data = datos, formula = y5 ~ x)
ajuste_6 <- lm(data = datos, formula = y6 ~ x)
compare_performance(ajuste_1, ajuste_2, ajuste_3, ajuste_4, ajuste_5, ajuste_6)
# Comparison of Model Performance Indices

Name     | Model |  AIC (weights) | AICc (weights) |  BIC (weights) |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------------------------------------------------------
ajuste_1 |    lm |  123.3 (>.999) |  123.3 (>.999) |  138.0 (>.999) | 0.997 |     0.997 | 0.257 | 0.257
ajuste_2 |    lm | 2868.3 (<.001) | 2868.3 (<.001) | 2883.0 (<.001) | 0.955 |     0.955 | 1.012 | 1.013
ajuste_3 |    lm | 3649.1 (<.001) | 3649.1 (<.001) | 3663.8 (<.001) | 0.904 |     0.904 | 1.496 | 1.497
ajuste_4 |    lm | 4294.8 (<.001) | 4294.8 (<.001) | 4309.5 (<.001) | 0.842 |     0.842 | 2.066 | 2.068
ajuste_5 |    lm | 5620.5 (<.001) | 5620.6 (<.001) | 5635.3 (<.001) | 0.571 |     0.570 | 4.008 | 4.012
ajuste_6 |    lm | 7087.1 (<.001) | 7087.1 (<.001) | 7101.8 (<.001) | 0.241 |     0.240 | 8.344 | 8.353

Ruido vs R2

En todos los casos el ajuste debe ser (aprox) el mismo ya que el modelo de regresión busca modelizar los efectos lineales no aleatorios entre la variable objetivo y los predictores.

Moraleja: tener un \(R^2\) no implica que el modelo sea malo, ya que la cantidad de información no modelizable puede deberse a una cantidad alta de ruido (algo aleatorio no predecible). Por ello es importante usar más herramientas que un mero coeficiente para valorar un ajuste (por ejemplo, en campos como la sociología o la economía la bondad de ajuste será generalmente bajo)

 

Deberes: ¿cómo ilustrar gráficamente que a mayor varianza del ruido, menor es \(R^2\)? Diseña un estudio de simulación para ello con distintos modelos y gráfica la caída de \(R^2\).

Diagnosis vs R2

Entonces, si tenemos un modelo con un alto \(R^2\), ¿no hace falta que cumpla las hipótesis?

Vamos a simular un modelo que incumple

  • Linealidad

  • Homocedasticidad

La variable predictora \(x_i = 0.01 + 0.01*(i-1)\) (\(i = 1,\ldots, n = 200\)) será la siguiente

x <- seq(0.01, 2, l = 200)

¿Cómo podríamos crear una \(Y\) cuya relación con \(X\) sea no lineal?

Diagnosis vs R2

¿Cómo podríamos crear una \(Y\) cuya relación con \(X\) sea no lineal?

Tenemos muchas maneras, por ejemplo:

\[Y = X + X^2 + \varepsilon, \quad \varepsilon \sim N(0, \sigma_{\varepsilon})\]

\[Y = \log(X^2) - cos(X) + \varepsilon, \quad \varepsilon \sim N(0, \sigma_{\varepsilon})\]

\[Y = \frac{1}{X + 1} + \varepsilon, \quad \varepsilon \sim N(0, \sigma_{\varepsilon})\]

\[Y = 1 - 2X(1 + 0.25 \sin(4 \pi X)) + \varepsilon, \quad \varepsilon \sim N(0, \sigma_{\varepsilon})\]

Simula este último modelo (con \(\sigma_{\varepsilon} = 0.5\)) y realiza el ajuste

Diagnosis vs R2

x <- seq(0.15, 2, l = 200)
eps <- rnorm(n = 200, sd = 0.5)
y <- 1 - 2 * x * (1 + 0.25*sin(4 * pi * x)) + eps
datos <- tibble("y" = y, "x" = x)
ajuste_lineal <- lm(data = datos, formula = y ~ x)
ajuste_lineal |> summary()

Call:
lm(formula = y ~ x, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5514 -0.3966  0.0531  0.4308  1.4051 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.93752    0.09947   9.425   <2e-16 ***
x           -1.89214    0.08279 -22.855   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6284 on 198 degrees of freedom
Multiple R-squared:  0.7251,    Adjusted R-squared:  0.7237 
F-statistic: 522.3 on 1 and 198 DF,  p-value: < 2.2e-16

¿Cómo podríamos hacer que la hipótesis de homocedasticidad no se cumpla

Diagnosis vs R2

Vamos a considerar que la varianza del error no es cte, crece según aumenta x

\[y_i = f(x_i) + \varepsilon_i, \quad \varepsilon_i \sim N(0, \sigma_{\varepsilon,i}), \quad \sigma_{\varepsilon,i} = 0.25 * x_{i}^2 \quad i=1, \ldots, 200\]

Código
x <- seq(0.01, 2, l = 200)
eps <- rnorm(n = 200, sd = 0.25 * x^2)
y <- 1 - 2 * x * (1 + 0.25*sin(4 * pi * x)) + eps
datos <- tibble("y" = y, "x" = x, "sigma" =  0.25 * x^2)
ajuste_lineal <- lm(data = datos, formula = y ~ x)
ajuste_lineal |> summary()

Call:
lm(formula = y ~ x, data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.88732 -0.31063  0.03454  0.20884  2.52492 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.88798    0.08718   10.19   <2e-16 ***
x           -1.82274    0.07522  -24.23   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6141 on 198 degrees of freedom
Multiple R-squared:  0.7478,    Adjusted R-squared:  0.7466 
F-statistic: 587.2 on 1 and 198 DF,  p-value: < 2.2e-16

Diagnosis vs R2

Código
ggplot(datos, aes(x = x, y = y)) +
  geom_point(aes(color = sigma), size = 3, alpha = 0.75) +
  geom_smooth(method = "lm") +
  theme_minimal()

Aunque el \(R^2\) es bastante alto, el modelo no tiene sentido, dando estimaciones cada vez más erradas según nos movemos hacia la derecha en el eje X.

Cosicas a examen

Algunas cosas que podrían caer en la primera entrega:

  • Realizar un análisis exploratorio inicial lo más completo posible. Elegir de manera justificada la mejor predictora
  • Saber ajustar e INTERPRETAR un modelo de regresión.
  • Realizar e interpretar la fase de diagnosis lo más completa posible para luego interpretar (si procede) la inferencia del modelo.
  • Saber hacer una correcta fase de evaluación dle modelo.
  • Saber simular modelos de regresión (tanto que cumplan las hipótesis como que no las cumplan).

Clase 10: predicción y resumen

Fase de predicción. Resumen

Resumen

Un breve resumen de lo aprendido sobre reg. lineal univariante

  • Modelo supervisado de predicción: hay una variable objetivo \(Y\) continua (numérica) cuyo valor real conocemos
  • La visualización importa: no fies tu análisis solo a los parámetros matemáticos, la visualización ayuda a entender los datos.
  • Relación entre variables: buscamos predictoras muy correladas (linealmente) con \(Y\) y lo más incorreladas/independientes entre ellas.
  • Modelo: el modelo parámetrico se resume en \(Y = \beta_0 + \beta_1 X + \varepsilon\), donde \(\varepsilon\) es una variable aleatoria (ruido).
  • Estimación: la estimación viene modelizada bajo la hipótesis de linealidad

\[E[Y|X=x] = \beta_0 + \beta_1 x\]

Clase 11: entrega I

Entrega I

Clase 12: incumpliendo hipótesis

¿Cómo simular datos que incumplan las hipótesis?

Incumpliendo hipótesis

Los contenidos vistos en esta clase se subirán resumidos en formato notebook.

Clase 13: reg. multivariante

Introducción a la regresión multivariante. Formulación del modelo y estimación

Formulación multivariante

De aquí en adelante llamaremos modelo multivariante a todo modelo en el que \(p > 1\) (es decir, tenemos más de una variable predictora).

\[Y = f\left(X_1, \ldots, X_p \right) + \varepsilon, \quad E \left[Y | \boldsymbol{X} = x \right] =f\left(X_1, \ldots, X_p \right)\] tal que \(E \left[ \varepsilon | \left( X_1 = x_1, \ldots, X_p = x_p \right) \right] = 0\).

En el caso del modelo lineal multivariante se traducirá por tanto en

\[E \left[Y | \boldsymbol{X} = x \right] = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p, \quad \widehat{Y} = \widehat{\beta_0} + \displaystyle \sum_{j=1}^{p} \widehat{\beta_j} X_j \]

El objetivo seguirá siendo obtener la estimación de los \(\widehat{\beta}\) tal que minimicemos el error

Formulación matricial

Su formulación muestral la podemos expresar mediante la matriz de diseño

\[\mathbf{X} =\begin{pmatrix} 1 & X_{11} & \ldots & X_{1p} \\ 1 & X_{21} & \ldots & X_{2p} \\ \vdots & \vdots & \ddots & \vdots\\ 1 & X_{n1} & \ldots & X_{np} \end{pmatrix}_{n\times(p+1)}\]

tal que \(Y_i = \beta_0 + \sum_{j=1}^{p} \beta_j X_j + \varepsilon_i\), para todo elemento de la muestra \(i=1,\ldots, n\)

\[ \mathbf{Y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} = \begin{pmatrix} 1 & X_{11} & \ldots & X_{1p} \\ 1 & X_{21} & \ldots & X_{2p} \\ \vdots & \vdots & \ddots & \vdots\\ 1 & X_{n1} & \ldots & X_{np} \end{pmatrix}_{n\times(p+1)}\begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{pmatrix}_{(p+1)\times1} + \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{pmatrix}_{n\times1}\]

Formulación matricial

La estimación será por tanto

\[\widehat{\mathbf{Y}} = \begin{pmatrix} \widehat{Y}_1 \\ \widehat{Y}_2 \\ \vdots \\ \widehat{Y}_n \end{pmatrix}_{n\times1} = \mathbf{X} \widehat{\boldsymbol{\beta}} = \begin{pmatrix} 1 & X_{11} & \ldots & X_{1p} \\ 1 & X_{21} & \ldots & X_{2p} \\ \vdots & \vdots & \ddots & \vdots\\ 1 & X_{n1} & \ldots & X_{np} \end{pmatrix}_{n\times(p+1)}\begin{pmatrix} \widehat{\beta}_1 \\ \widehat{\beta}_2 \\ \vdots \\ \widehat{\beta}_p \end{pmatrix}_{(p+1)\times1}\]

Varianza residual

Como sucedía antes, el objetivo será minimizar la varianza residual o error cuadrático medio (varianza del error), o lo que hemos llamado \(SSE\)

\[SSE \left(\boldsymbol{\beta} \right) = \displaystyle \sum_{i=1}^{n} \widehat{\varepsilon}_{i}^{2} = \displaystyle \sum_{i=1}^{n} \left(Y_i - \beta_0 - \beta_1 X_{i1} - \ldots - \beta_p X_{ip} \right)^2 \] ¿Cómo quedaría matricialmente?

\[SSE \left(\boldsymbol{\beta} \right) = \displaystyle \sum_{i=1}^{n} \widehat{\varepsilon}_{i}^{2} = \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)^{T}\left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right) = \varepsilon^{T} \varepsilon\]

Así, aplicando el método de los mínimos cuadrados, los parámetros óptimos serán aquellos \(\widehat{\boldsymbol{\beta}} = \left(\widehat{\beta}_0, \widehat{\beta}_1, \ldots, \widehat{\beta}_p\right)\) que minimicen dicha suma

\[\widehat{\boldsymbol{\beta}} = \left(\widehat{\beta}_0, \widehat{\beta}_1, \ldots, \widehat{\beta}_p \right) = \arg \min_{\boldsymbol{\beta} \in \mathbb{R}^{p+1}} \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)^{T}\left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)\]

Y repetimos la idea: calcularemos la derivada respecto a \(\widehat{\boldsymbol{\beta}}\) e igualamos a cero

Método mínimos cuadrados

\[\widehat{\boldsymbol{\beta}} = \left(\widehat{\beta}_0, \widehat{\beta}_1, \ldots, \widehat{\beta}_p \right) = \arg \min_{\boldsymbol{\beta} \in \mathbb{R}^{p+1}} \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)^{T}\left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)\]

Derivando matricialmente tenemos (recuerda: \(\left(A B \right)^{T} = B^{T} A^{T}\))

\[\begin{eqnarray}\frac{\partial SSE \left( \boldsymbol{\beta} \right)}{\partial \boldsymbol{\beta}} &=& \frac{\partial \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)^{T}\left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)}{\partial \boldsymbol{\beta}} \nonumber \\ &=& \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)^{T} \frac{\partial \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)}{\partial \boldsymbol{\beta}} + \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right) \frac{\partial \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)^{T}}{\partial \boldsymbol{\beta}} \nonumber \\ &=& -\left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)^{T} \mathbf{X} - \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right) \mathbf{X}^{T} \nonumber \\ &=& -\mathbf{X}^{T}\left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right) - \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right) \mathbf{X}^{T} \nonumber \\ &=& -2\mathbf{X}^{T} \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right) = 0 \end{eqnarray}\]

Método mínimos cuadrados

\[\begin{eqnarray}\frac{\partial SSE \left( \boldsymbol{\beta} \right)}{\partial \boldsymbol{\beta}} &=& -2\mathbf{X}^{T} \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right) = 0 \end{eqnarray}\]

Si descartamos el \(-2\) como constante, tenemos que

\[\mathbf{X}^{T} \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right) = 0 \longrightarrow \mathbf{X}^{T} \mathbf{Y} = \mathbf{X}^{T} \mathbf{X} \boldsymbol{\beta}\]

Si ahora despejamos \(\beta\) multiplicando en ambos lados por \(\left(\mathbf{X}^{T}\mathbf{X} \right)^{-1}\) tenemos finalmente

\[\widehat{\boldsymbol{\beta}}=\left(\mathbf{X}^{T}\mathbf{X} \right)^{-1}\mathbf{X}^{T}\mathbf{Y}\]

Estimación

Fíjate que cuando \(p = 1\) tenemos que

\[\widehat{\mathbf{Y}} = \begin{pmatrix} \widehat{Y}_1 \\ \vdots \\ \widehat{Y}_n \end{pmatrix}_{n\times1} = \mathbf{X} \widehat{\boldsymbol{\beta}} = \begin{pmatrix} 1 & X_{11} \\\vdots & \vdots \\ 1 & X_{n1} \end{pmatrix}_{n \times 2}\begin{pmatrix} \widehat{\beta}_0 \\ \widehat{\beta}_1 \end{pmatrix}_{2\times1}\]

y que por tanto \(\widehat{\boldsymbol{\beta}}= \left(\begin{pmatrix} 1 & \ldots & 1 \\ X_{11} & \ldots & X_{n1} \end{pmatrix} \begin{pmatrix} 1 & X_{11} \\ 1 & X_{21} \\ \vdots & \vdots \\ 1 & X_{n1} \end{pmatrix} \right)^{-1}\mathbf{X}^{T}\mathbf{Y}\)

Deberes: sería interesante que para el examen supieses demostrar que cuando \(p=1\) esa expresión acaba en

\[\widehat{\beta}_1 = \frac{s_{xy}}{s_{x}^{2}}, \quad \widehat{\beta}_0 = \overline{y} - \hat{\beta}_1 \overline{x}\]

Estimación

\[\widehat{\boldsymbol{\beta}}=\left(\mathbf{X}^{T}\mathbf{X} \right)^{-1}\mathbf{X}^{T}\mathbf{Y}\]

  • La estimación será por tanto

\[\widehat{\mathbf{Y}} = \mathbf{X}\hat{\boldsymbol{\beta}}=\mathbf{X}\left(\mathbf{X}^{T}\mathbf{X} \right)^{-1}\mathbf{X}^{T}\mathbf{Y} = \mathbf{H}\mathbf{Y}\]

  • La matriz \(\mathbf{H} = \mathbf{X}\left(\mathbf{X}^{T}\mathbf{X} \right)^{-1}\mathbf{X}^{T}\) se conoce como hat matrix o matriz de proyección (ya que hace que las estimaciones \(\hat{y}\) sean en realidad los valores \(y\) proyectados verticalmente sobre el plano de regresión ajustado).

Interpretación

Al igual que pasaba antes será importante distinguir lo población de lo muestral

  • Los parámetros \(\left( \beta_0, \beta_1, \ldots, \beta_p \right)\) son desconocidos, los parámetros poblacionales

  • Los parámetros \(\left( \widehat{\beta}_0, \widehat{\beta}_1 , \ldots, \widehat{\beta}_p\right)\) son estimados a partir de los datos, son variables aleatorias ya que han sido calculados en función de una muestra aleatoria.

 

¿Cómo se interpretan ahora los parámetros?

Interpretación

\[\widehat{\mathbf{Y}} = \mathbf{X}\widehat{\boldsymbol{\beta}}= \widehat{\beta}_0 + \sum_{j=1}^{p} \widehat{\beta}_j X_j\]

  • Ordenada en el origen: también llamado intercepto, y denotado como \(\beta_0\) su valor real, es el valor de \(Y\) cuando \(X_1 = \ldots = X_p = 0\). Es decir, \(\widehat{\beta}_0\) se puede interpretar como la estimación \(\widehat{Y}\) cuando TODAS las predictoras son nulas
  • Pendientes: denotadas como \(\beta_j\), para \(j=1,\ldots, p\), su valor real, cuantifica el incremento de \(Y\) cuando solo \(X_j\) aumenta una unidad. Es decir, \(\widehat{\beta}_j\) se puede interpretar como la variación de la estimación \(\widehat{Y}\) cuando \(X_j\) tiene un incremento unitario y el RESTO DE PREDICTORAS permanecen fijas.

Diagnosis

En el caso multivariante las 4 hipótesis se convierten en

  1. Linealidad: el valor esperado de \(Y\) es

\[E \left[Y | \left( \boldsymbol{X}_1 = x_1, \ldots, \boldsymbol{X}_p = x_p \right) \right] = \beta_0 + \beta_1 x_1 + \ldots + \beta_p x_p\]

  1. Homocedasticidad: necesitamos que la varianza del error sea finita y constante \(\sigma_{r}^2 = \sigma_{\varepsilon}^2 = {\rm Var} \left[\varepsilon | \left( \boldsymbol{X}_1 = x_1, \ldots, \boldsymbol{X}_p = x_p \right) \right] = cte < \infty\)

  2. Normalidad: pediremos que \(\varepsilon \sim N \left(0, \sigma_{r}^2 \right)\)

  3. Independencia: los errores \(\left\lbrace \varepsilon_i \right\rbrace_{i=1}^{n}\) deben ser en independientes, y en particular incorrelados \({\rm Cor}_{\varepsilon_i \varepsilon_j} = E \left[\varepsilon_i \varepsilon_j \right] = 0, \quad i \neq j\)

\[Y | \left( \boldsymbol{X}_1 = x_1, \ldots, \boldsymbol{X}_p = x_p \right) \sim N \left(\beta_0 + \beta_1x_1 + \ldots + \beta_p x_p, \sigma_{\varepsilon}^2 \right)\]

Inferencia

Ahora las hipótesis nos permiten decir que los parámetros estimados siguen una distribución normal multivariante de media el vector de parámetros a estimar y de matriz de covarianzas la varianza residual por \(\left(\mathbf{X}^{T}\mathbf{X} \right)^{-1}\)

\[\widehat{\boldsymbol{\beta}} \sim N_{p+1}\left(\boldsymbol{\beta}, \sigma_{\varepsilon}^2 \left( \mathbf{X}^{T}\mathbf{X} \right)^{-1}\right)\]

Esa normal multivariante, componente a componente, deriva en

\[\widehat{\beta}_j\sim N \left(\beta_j, SE \left(\widehat{\beta}_j\right)^2 \right), \quad \frac{\widehat{\beta}_j - \beta_j}{SE \left(\widehat{\beta}_j\right)} \sim N \left(0, 1 \right), \quad \frac{\widehat{\beta}_j - \beta_j}{\widehat{SE} \left(\widehat{\beta}_j\right)} \sim t_{n-p-1}\]

donde \(SE \left(\widehat{\beta}_j\right)^2 = \sigma_{\varepsilon}^2 v_j\) y \(\widehat{SE} \left(\widehat{\beta}_j\right)^2 = \widehat{\sigma}_{\varepsilon}^2 v_j\), con \(v_j\) el elemento \(j\)-ésimo de la diagonal de \(\left( \mathbf{X}^{T}\mathbf{X} \right)^{-1}\) (matriz que representa la variabilidad de predictores)

Caso práctico: wine.csv

Vamos a volver a usar nuestros datos wine.csv

datos <- read_csv(file = "./datos/wine.csv")
datos
# A tibble: 27 × 7
    Year Price WinterRain  AGST HarvestRain   Age FrancePop
   <dbl> <dbl>      <dbl> <dbl>       <dbl> <dbl>     <dbl>
 1  1952  7.50        600  17.1         160    31    43184.
 2  1953  8.04        690  16.7          80    30    43495.
 3  1955  7.69        502  17.2         130    28    44218.
 4  1957  6.98        420  16.1         110    26    45152.
 5  1958  6.78        582  16.4         187    25    45654.
 6  1959  8.08        485  17.5         187    24    46129.
 7  1960  6.52        763  16.4         290    23    46584.
 8  1961  8.49        830  17.3          38    22    47128.
 9  1962  7.39        697  16.3          52    21    48089.
10  1963  6.71        608  15.7         155    20    48799.
# ℹ 17 more rows

Pero esta vez no vamos a seleccionar ninguna variable previamente. Para ajustar un modelo multivariante basta con añadir variables +

ajuste_uni <- lm(data = datos, formula = Price ~ AGST)
ajuste_2 <- lm(data = datos, formula = Price ~ AGST + HarvestRain)
ajuste_full <- lm(data = datos, formula = Price ~ .)

Caso práctico: wine.csv

Fíjate que ahora la tabla coefficients tiene una línea por covariable (más \(\beta_0\)) y además en este caso dice (1 not defined because of singularities): la matriz de covarianzas no es invertible ya que el determinante es 0 (en este caso Year es “igual” que Age)

ajuste_full |> summary()

Call:
lm(formula = Price ~ ., data = datos)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46541 -0.24133  0.00413  0.18974  0.52495 

Coefficients: (1 not defined because of singularities)
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.496e+01  1.081e+02   0.231  0.81956    
Year        -1.377e-02  5.821e-02  -0.237  0.81531    
WinterRain   1.153e-03  4.991e-04   2.311  0.03109 *  
AGST         6.144e-01  9.799e-02   6.270 3.22e-06 ***
HarvestRain -3.837e-03  8.366e-04  -4.587  0.00016 ***
Age                 NA         NA      NA       NA    
FrancePop   -2.213e-05  1.268e-04  -0.175  0.86313    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.293 on 21 degrees of freedom
Multiple R-squared:  0.8278,    Adjusted R-squared:  0.7868 
F-statistic: 20.19 on 5 and 21 DF,  p-value: 2.232e-07

Caso práctico: wine.csv

Una vez eliminada, en la última línea F-statistic: ... p-value: 2.232e-07: se está rechazando la hipótesis nula del contraste de significación global (existe alguna predictora cuyo efecto lineal es significativo). Recuerda que la nula es \(H_0:~\beta_1 = \ldots = \beta_5 = 0\)

ajuste_full <- lm(data = datos |> select(-Age), formula = Price ~ .)
ajuste_full |> summary()

Call:
lm(formula = Price ~ ., data = select(datos, -Age))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46541 -0.24133  0.00413  0.18974  0.52495 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.496e+01  1.081e+02   0.231  0.81956    
Year        -1.377e-02  5.821e-02  -0.237  0.81531    
WinterRain   1.153e-03  4.991e-04   2.311  0.03109 *  
AGST         6.144e-01  9.799e-02   6.270 3.22e-06 ***
HarvestRain -3.837e-03  8.366e-04  -4.587  0.00016 ***
FrancePop   -2.213e-05  1.268e-04  -0.175  0.86313    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.293 on 21 degrees of freedom
Multiple R-squared:  0.8278,    Adjusted R-squared:  0.7868 
F-statistic: 20.19 on 5 and 21 DF,  p-value: 2.232e-07

Caso práctico: wine.csv

Si nos fijamos en la tabla de coeficientes el modelo ajustado es

\[\begin{eqnarray}\widehat{Precio} &=& 2.496 - 0.0137 * Year + 0.0012 * WinterRain + 0.614 * AGST \nonumber \\ &-& 0.0038 * HarvestRain - 0.00002 * FrancePop \nonumber \end{eqnarray}\]

ajuste_full |> summary()

Call:
lm(formula = Price ~ ., data = select(datos, -Age))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46541 -0.24133  0.00413  0.18974  0.52495 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.496e+01  1.081e+02   0.231  0.81956    
Year        -1.377e-02  5.821e-02  -0.237  0.81531    
WinterRain   1.153e-03  4.991e-04   2.311  0.03109 *  
AGST         6.144e-01  9.799e-02   6.270 3.22e-06 ***
HarvestRain -3.837e-03  8.366e-04  -4.587  0.00016 ***
FrancePop   -2.213e-05  1.268e-04  -0.175  0.86313    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.293 on 21 degrees of freedom
Multiple R-squared:  0.8278,    Adjusted R-squared:  0.7868 
F-statistic: 20.19 on 5 and 21 DF,  p-value: 2.232e-07

Caso práctico: wine.csv

\[\begin{eqnarray}\widehat{Precio} &=& 2.496 - 0.0137 * Year + 0.0012 * WinterRain + 0.614 * AGST \nonumber \\ &-& 0.0038 * HarvestRain - 0.00002 * FrancePop \nonumber \end{eqnarray}\]

La interpretación es

  • La estimación del modelo cuando todas las covariables son 0 (sin población, sin lluvia ni temperatura, año 0) es de 2.496 (escala log)

  • Por ejemplo, si el resto de variables permaneciesen fijas, por cada litro de lluvia que caiga de más durante la cosecha (agosto-septiembre), el modelo estima que el precio baja 3.8 (escala log)

Caso práctico: wine.csv

Si nos fijamos en la tabla de coeficientes tenemos 2 predictoras cuyo efecto lineal no se acepta que sea significativo: Year, FrancePop (además de \(\beta_0\)) ya que el contraste de significación \(H_0:~\beta_j = 0\) vs \(H_1:~\beta_j \neq 0\) nos devuelve un \(p-valor > \alpha\)

ajuste_full |> summary()

Call:
lm(formula = Price ~ ., data = select(datos, -Age))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46541 -0.24133  0.00413  0.18974  0.52495 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.496e+01  1.081e+02   0.231  0.81956    
Year        -1.377e-02  5.821e-02  -0.237  0.81531    
WinterRain   1.153e-03  4.991e-04   2.311  0.03109 *  
AGST         6.144e-01  9.799e-02   6.270 3.22e-06 ***
HarvestRain -3.837e-03  8.366e-04  -4.587  0.00016 ***
FrancePop   -2.213e-05  1.268e-04  -0.175  0.86313    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.293 on 21 degrees of freedom
Multiple R-squared:  0.8278,    Adjusted R-squared:  0.7868 
F-statistic: 20.19 on 5 and 21 DF,  p-value: 2.232e-07

Caso práctico: wine.csv

La pregunta que intentaremos resolver en futuras clases es:

¿Cómo decidir que predictoras seleccionamos? El problema si quitamos todas las no significativas de manera simultánea es que no sabemos qué efectos puede haber entre las propias predictoras

ajuste_full |> summary()

Call:
lm(formula = Price ~ ., data = select(datos, -Age))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46541 -0.24133  0.00413  0.18974  0.52495 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.496e+01  1.081e+02   0.231  0.81956    
Year        -1.377e-02  5.821e-02  -0.237  0.81531    
WinterRain   1.153e-03  4.991e-04   2.311  0.03109 *  
AGST         6.144e-01  9.799e-02   6.270 3.22e-06 ***
HarvestRain -3.837e-03  8.366e-04  -4.587  0.00016 ***
FrancePop   -2.213e-05  1.268e-04  -0.175  0.86313    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.293 on 21 degrees of freedom
Multiple R-squared:  0.8278,    Adjusted R-squared:  0.7868 
F-statistic: 20.19 on 5 and 21 DF,  p-value: 2.232e-07

Caso práctico: wine.csv

Fíjate que ahora en la fase de diagnosis tenemos una sexta gráfica a chequear: una gráfica que nos calcula la conocida como VIF (Variance Inflaction), que nos cuantifica la colinealidad (efectos lineales) entre las predictoras

check_model(ajuste_full)

ANOVA

Ahora anova() en realidad nos hace una ANOVA secuencial (variable a variable)

ajuste_full |> anova()
Analysis of Variance Table

Response: Price
            Df Sum Sq Mean Sq F value    Pr(>F)    
Year         1 2.2195  2.2195 25.8480 4.916e-05 ***
WinterRain   1 0.2635  0.2635  3.0693 0.0943762 .  
AGST         1 4.3055  4.3055 50.1419 5.495e-07 ***
HarvestRain  1 1.8760  1.8760 21.8485 0.0001298 ***
FrancePop    1 0.0026  0.0026  0.0305 0.8631279    
Residuals   21 1.8032  0.0859                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA

Ahora anova() en realidad nos hace una ANOVA secuencial (variable a variable)

Df (grados) Sum Sq Mean Sq F value Pr(>F)
\(x_1\) 1 \(SSR_1\) \(\frac{SSR_1}{1}\) \(F-value = \frac{\frac{SSR_1}{1}}{\frac{SSE}{n-p-1}}\) \(p_1\)
\(F-value = \frac{\frac{SSR}{p}}{\frac{SSE}{n-p-1}}\) p-valor F test
\(x_p\) 1 \(SSR_p\) \(\frac{SSR_p}{1}\) \(F-value = \frac{\frac{SSR_p}{1}}{\frac{SSE}{n-p-1}}\) \(p_p\)
Residuals n - p - 1 SSE \(\widehat{\sigma}_{\varepsilon}^{2} = \frac{SSE}{n-p-1}\)

ANOVA

ajuste_full |> anova()
Analysis of Variance Table

Response: Price
            Df Sum Sq Mean Sq F value    Pr(>F)    
Year         1 2.2195  2.2195 25.8480 4.916e-05 ***
WinterRain   1 0.2635  0.2635  3.0693 0.0943762 .  
AGST         1 4.3055  4.3055 50.1419 5.495e-07 ***
HarvestRain  1 1.8760  1.8760 21.8485 0.0001298 ***
FrancePop    1 0.0026  0.0026  0.0305 0.8631279    
Residuals   21 1.8032  0.0859                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ahora \(SSR_j\), con \(j=1,\ldots,p\), representa la suma de residuos al cuadrado \(SSR\) cuando incluimos la predictora j-ésima \(X_j\), frente a cuando no lo hacemos, tal que

\[SSR_j = SSR \left(X_1, \ldots, X_j \right) - SSR \left(X_1, \ldots, X_{j-1} \right)\]

Así los p-valores \(p_j\) son individuales (los mismos que los de la tabla de coeficientes).

Clase 14: selección de modelos

Selección secuencial de modelos

Bondad de ajuste

Como ya hemos hablado, el conocido como \(R^2\) o bondad de ajuste o el coeficiente de determinación es una forma sencilla, y relacionada con el ANOVA, para cuantificar la cantidad de información explicada por el modelo definida como

\[R^2 = \text{Ratio info explicada} = \frac{SSR}{SST} =_{lineal} 1 - \frac{SSE}{SST} = 1 - \frac{s_{r}^{2}}{s_{y}^{2}} = 1 - \frac{\frac{n-p-1}{n}\widehat{\sigma}_{\varepsilon}^{2}}{s_{y}^{2}}\]

Por definición, la bondad de ajuste está entre 0 y 1 pero además presenta un problema en el caso multivariante

  • Más predictoras implica que \(R^2\) crece (al aumentar \(p\)), ¡incluso aunque dichas predictoras no sean significativas!

R2 ajustado

Para evitar dicho problema vamos a definir la conocida como bondad de ajuste ajustada

\[R_{adj}^2 = 1 - \frac{\frac{SSE}{n-p-1}}{\frac{SST}{n-1}} = 1 - \frac{SSE}{SST}\frac{n-1}{n-p-1} = 1 - \frac{\widehat{\sigma}_{\varepsilon}^{2}}{SST}(n-1)\]

Así el \(R_{adj}^2\) no depende del número de predictoras de manera directa (si sigue dependiendo del ruido, de manera que descenderá solo cuando incrementar \(p\) implica reducir el error, es decir, variables con un efecto significativo). Se cumple además que que si \(p=1\)

\[\lim_{n\to \infty} R_{adj}^2 = \lim_{n\to \infty} \left(1 - \frac{\widehat{\sigma}_{\varepsilon}^{2}}{SST}(n-1) \right) = \lim_{n\to \infty} \left(1 - \frac{\widehat{\sigma}_{\varepsilon}^{2}}{SST}(n-p-1) \right) = R^2\]

R2 ajustado

Para ver mejor su diferencia debes realizar el siguiente estudio de simulación (si quieres fija semilla set.seed(12345)):

  1. Considera dos predictoras \(X_1 \sim N(0, 1)\) y \(X_2 \sim N(0, 1)\) de tamaño muestral \(n = 200\). Considera el ruido como \(\varepsilon \sim N(0, \sigma = 12)\)
  1. Simula el modelo \(Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \varepsilon_i\) con \(\beta_0 = 0.01\), \(\beta_1 = 1.5\) y \(\beta_2 = -1.5\)
  1. Realiza el ajuste con todas las predictoras y calcula el \(R^2\) y el \(R_{adj}^2\) (lo necesitamos guardar, así usa la fórmula no la salida del lm())
  1. Repite este proceso \(M = 300\) veces (vuelve a simular las variables, vuelve a construir el modelo, y obtén de nuevo las bondades de ajuste, de manera que tengamos \(300\) de cada uno)

R2 ajustado

Lo que debería salirte si fijas set.seed(12345)

R2 ajustado

Repite el estudio de simulación añadiendo

  1. Una nueva predictora “basura” \(X_3 \sim N(0, 2)\) tal que \(Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + 10^{-10} X_{i3} + \varepsilon_i\).

  2. Dos nuevas predictora “basura” \(X_3 \sim N(0, 2)\) y \(X_4 \sim N(0, 2)\) tal que \(Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + 10^{-10} X_{i3} + 10^{-10} X_{i4} + \varepsilon_i\).

  1. 195 nuevas predictora “basura” \(X_j \sim N(0, 2)\) tal que \(Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \displaystyle \sum_{j=3}^{195} 10^{-10} X_{ij} + \varepsilon_i\).

Ahora debes obtener \(300*195\) valores de \(R^2\) y \(R_{adj}^2\) (300 simulaciones por cada nueva predictora que añadimos)

R2 ajustado

Multicolinealidad

Uno de los mayores problemas de los modelos lineales multivariantes es el conocido como problema de colinealidad

Imagina un modelo \(Y = \beta_0 + \beta_1 X_1 + \beta_2X_2 + \beta_3 X_3\), donde \(X_3\) es definida como (literal) la suma de \(X_1\) y \(X_2\).

n <- 200
x_1 <- rnorm(n, sd = 2)
x_2 <- rnorm(n, sd = 2)
x_3 <- x_1 + x_2
eps <- rnorm(n)

y <- 0.5 - 2*x_1 + 3*x_2 -2.5*x_3
datos <- tibble(y, x_1, x_2, x_3)
ajuste <- lm(data = datos, formula = y ~ .)

Multicolinealidad

Como ya vimos en el ejemplo de wine.csv, cuando tenemos dos o más predictoras que dependen linealmente entre sí la matriz \(\left(X^{T} X \right)^{-1}\) no se puede invertir ya que es singular (determinante igual a 0), así que debe eliminar una de las ecuaciones para que el problema sea de rango completo

Esto se conoce como colinealidad exacta

ajuste |> summary()

Call:
lm(formula = y ~ ., data = datos)

Residuals:
       Min         1Q     Median         3Q        Max 
-5.595e-15 -1.073e-15 -4.250e-16  1.710e-16  8.243e-14 

Coefficients: (1 not defined because of singularities)
              Estimate Std. Error    t value Pr(>|t|)    
(Intercept)  5.000e-01  4.264e-16  1.172e+15   <2e-16 ***
x_1         -4.500e+00  2.098e-16 -2.145e+16   <2e-16 ***
x_2          5.000e-01  2.169e-16  2.305e+15   <2e-16 ***
x_3                 NA         NA         NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.03e-15 on 197 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 2.342e+32 on 2 and 197 DF,  p-value: < 2.2e-16

Multicolinealidad

La colinealidad exacta es solo el ejemplo más extremo de lo que se conoce como problema de colinealidad: un problema que aparece cuando varias predictoras están altamente correladas. Un problema de colinealidad tiene principalmente dos consecuencias:

  • Reduce la precisión de los estimadores ya que la matriz \(X^{T} X\), según aumenta la dependencia, tiene un determinante cada vez más cercano a 0, por lo que \(\left(\mathbf{X}^{T}\mathbf{X} \right)^{-1}\), que siempre depende de \(1/det\left( \mathbf{X}^{T}\mathbf{X} \right)\), tendrá valores cada vez más grandes. Dado que \(\widehat{\beta}_j\sim N \left(\beta_j, SE \left(\widehat{\beta}_j\right)^2 \right)\), con \(SE \left(\widehat{\beta}_j\right)^2 = \sigma_{\varepsilon}^2 v_j\) siendo \(v_j\) el elemento \(j\)-ésimo de la diagonal de \(\left( \mathbf{X}^{T}\mathbf{X} \right)^{-1}\), implica que a mayor colinealidad, mayor varianza de las estimaciones
  • Distorsiona los efectos de los predictores sobre la respuesta

Multicolinealidad

¿Pero qué sucede cuando esa colinealidad no es tan obvia? Veamos un ejemplo sencillo: simula el siguiente modelo (para \(n = 100\))

\[Y = 1 + 0.5 X_1 + 2 X_2 - 3 X_3 - X_4 + \varepsilon, \quad \varepsilon \sim N(0, 1)\]

  • \(X_1 \sim N(0, 1)\), \(X_2 = 0.5*X_1 + N(0, 1)\) y \(X_3 = 0.5*X_2 + N(0, 1)\)

  • \(X_4 = -X_1 + X_2 + N(0, 0.5)\)

Código
set.seed(12345)
n <- 100

x1 <- rnorm(n)
x2 <- 0.5 * x1 + rnorm(n)
x3 <- 0.5 * x2 + rnorm(n)
x4 <- -x1 + x2 + rnorm(n, sd = 0.5)
eps <- rnorm(100)

y <- 1 + 0.5 * x1 + 2 * x2 - 3 * x3 - x4 + eps
datos <- tibble(y, x1, x2, x3, x4)

¿Cómo serán sus correlaciones (lineales)?

Multicolinealidad

A priori no se observa un problema obvio de colinealidad ya que no hay dos predictoras altamente correladas.

datos |> 
  correlate()
# A tibble: 5 × 6
  term       y     x1     x2     x3     x4
  <chr>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1 y     NA      0.441  0.147 -0.653 -0.213
2 x1     0.441 NA      0.550  0.279 -0.311
3 x2     0.147  0.550 NA      0.464  0.554
4 x3    -0.653  0.279  0.464 NA      0.142
5 x4    -0.213 -0.311  0.554  0.142 NA    

Multicolinealidad

Por la definición realizada es obvio que hay una dependencia lineal entre todas las predictoras, pero al ser una relación lineal más compleja (con ruido de por medio y distintas predictoras interactuando a la vez), una simple revisión de las correlaciones es necesaria pero no suficiente ya que puede ocultar problemas de colinealidad que sí existen.

datos |> 
  cor() |> 
  corrplot(method = "color", addCoef.col = "#121212")

Multicolinealidad

Necesitamos por tanto una forma de cuantificar dicha multicolinealidad, y para ello usaremos el conocido como factor de inflación de la varianza (VIF)

Para cada coeficiente \(\beta_j\), y su estimador \(\widehat{\beta}_j\), se define el VIF como

\[VIF\left(\widehat {\beta_j} \right) = \frac{1}{1 - R^2_{X_j | X_{-j}}}, \quad R_{j}^2:= R^2_{X_j | X_{-j}}\]

donde \(R_{j}^2\) representa la bondad de ajuste \(R^2\) cuando predecimos con una regresión lineal multivariante la predictora \(X_j\) en función del resto de predictoras (sin la variable objetivo)

\[X_j = \gamma_0 + \gamma_1 X_1 + \ldots + \gamma_{j-1} X_{j-1} + \gamma_{j+1} X_{j+1} + \ldots + \gamma_j X_j\]

Multicolinealidad

\[VIF\left(\widehat {\beta_j} \right) = \frac{1}{1 - R^2_{X_j | X_{-j}}}, \quad R_{j}^2:= R^2_{X_j | X_{-j}}\]

La idea es cuantificar como de relacionado está cada predictor respecto al resto (de manera conjunta, no dos a dos). Nota: siempre es mayor que 1

  • Si \(VIF\left(\widehat {\beta_j} \right)\) cercano a 1, entonces \(R_{j}^{2}\) cercano a 0 –> variabilidad del predictor \(X_j\) no se puede explicar con la combinación lineal de otras
  • Si \(VIF\left(\widehat {\beta_j} \right)\) mayor que 5, entonces \(1 - R_{j}^{2} < 0.2\), ergo \(R_{j}^{2} > 0.8\) –> problema moderado de colinealidad
  • Si además \(VIF\left(\widehat {\beta_j} \right)\) mayor que 10, entonces \(R_{j}^{2} > 0.9\) –> problema grave de colinealidad

Multicolinealidad

Para calcular el \(VIF\) podemos o bien hacer regresión con cada parámetro, de manera manual, o bien hacer uso de check_collinearity() del paquete {performance}

ajuste <- lm(data = datos, y ~ .)
performance::check_collinearity(ajuste)
# Check for Multicollinearity

Low Correlation

 Term  VIF    VIF 95% CI Increased SE Tolerance Tolerance 95% CI
   x3 1.38 [1.16,  1.88]         1.18      0.72     [0.53, 0.86]

Moderate Correlation

 Term  VIF    VIF 95% CI Increased SE Tolerance Tolerance 95% CI
   x1 7.07 [5.10,  9.99]         2.66      0.14     [0.10, 0.20]
   x4 7.29 [5.25, 10.30]         2.70      0.14     [0.10, 0.19]

High Correlation

 Term   VIF    VIF 95% CI Increased SE Tolerance Tolerance 95% CI
   x2 10.41 [7.42, 14.79]         3.23      0.10     [0.07, 0.13]

Multicolinealidad

Si te fijas \(VIF \left(\widehat{\beta}_2 \right) = 10.41\), que es justamente el \(1/(1-R^2)\) si hacemos la regresión \(X_2\) vs el resto de predictoras

check_collinearity(ajuste)$VIF
[1]  7.068885 10.408006  1.381455  7.285520
ajuste_X2 <- lm(data = datos |> select(-y), x2 ~ .)
ajuste_X2 |> summary()

Call:
lm(formula = x2 ~ ., data = select(datos, -y))

Residuals:
     Min       1Q   Median       3Q      Max 
-1.14613 -0.20489  0.04263  0.21392  1.09615 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.05919    0.03906  -1.516    0.133    
x1           0.80757    0.03829  21.090  < 2e-16 ***
x3           0.16945    0.03925   4.317 3.85e-05 ***
x4           0.74732    0.03355  22.274  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3791 on 96 degrees of freedom
Multiple R-squared:  0.9039,    Adjusted R-squared:  0.9009 
F-statistic: 301.1 on 3 and 96 DF,  p-value: < 2.2e-16

Multicolinealidad

El VIF podemos además visualizarlo en una sexta gráfica cuando hacemos check_model()

check_model(ajuste)

Multicolinealidad

El VIF podemos además visualizarlo de manera manual

Código
VIF <- tibble("variable" = c("x1", "x2", "x3", "x4"),
              "VIF"= check_collinearity(ajuste)$VIF)
ggplot(VIF) +
  geom_col(aes(x = variable, y = VIF, fill = VIF),
           alpha = 0.75) +
  scale_fill_gradient2(low = "#1B8A54", mid = "#F7E865",
                       high = "#CE2424", midpoint = 5) +
  theme_minimal()

Selección de modelos

Por lo tanto parece que obvio que en un modelo de regresión multivariante vamos a tener que que seleccionar distintas variables para solventar esos problemas de colinealidad

Llamaremos modelo saturado al modelo con todas las \(p\) predictoras, sin seleccionar.

Y como hemos visto, no podremos simplemente comparar 2 a 2 correlaciones o p-valores ya que el efecto de una variable nos puede afectar en otras. ¿Qué hacer?

Selección de modelos

Un protocolo de actuación habitual podría ser el siguiente:

  1. Hacemos el ajuste del modelo saturado

  2. Calculamos el VIF de cada estimador

  3. Adoptamos un umbral (por ejemplo, \(VIF > 10\)), de manera que toda predictora que lo supere se elimina ya de antemano

  4. Volvemos a chequear el VIF. Si hay que eliminar, volver al paso 3.

  5. Del resto de predictoras se hace una selección más fina. Algunas opciones: BIC/AIC, regresión penalizada (LASSO, ridge, elastic net), PCA.

Selección de modelos

Lo que nos dice el VIF y lo visto hasta ahora es que, incluso aunque tengan efecto, añadir predictoras complejizando el modelo no es gratis, ya que lo hacemos a costa de sobreajustar el modelo (la varianza de los estimadores se incrementa considerablemente)

¿Cuál es el número mínimo/máximo de predictoras que podremos incluir?

  • El número mínimo de predictoras será \(p = 1\)
  • El número máximo será \(p = n-2\), o dicho de otra forma, necesito al menos \(n \geq p+2\) observaciones. Piensa en \(p=1\): si quiero una recta, necesito al menos 2 puntos para que la recta existe y al menos 3 para poder calcular la varianza residual estimada (acuérdate que se divide entre \(n-p-1\))

De manera resumida: dado un modelo con \(p\) predictores, necesitaremos estimar \(p+2\) incógnitas (\(p+1\) coeficientes y la varianza residual).

BIC y AIC

¿Cómo seleccionar los predictores más adecuados para el CONJUNTO del ajuste?

Los métodos más conocidos son los conocidos como selección de modelos stepwise (paso-a-paso), que de manera iterativa, va incluyendo y descartando distintos predictores, y comparando su calidad, para decidir que combinación de parámetros es la más óptima.

La idea es seleccionar el modelo más óptimo en función de un criterio de información que combina la calidad del modelo con el número de predictoras empleadas: vamos a penalizar el uso de variables que no mejore suficiente el modelo

BIC y AIC

Los dos criterios de información más famosos son el Bayesian Information Criterion (BIC) y el Akaike Information Criterion (AIC) definidos ambos como

\[AIC/BIC_{modelo} = -Calidad + \underbrace{\text{npar(modelo)} * \text{penalización}}_{\text{Complejidad}}\]

  • \(\text{npar(modelo)}\) es el número de parámetros a estimar del modelo, que en el caso que nos ocupa es \(\text{npar(modelo)} = p+2\) (coeficientes + var residual)

  • \(\text{Calidad}\) una medida que nos diga como de bueno es nuestro ajuste: al tener signo negativo, buscamos el menor valor de AIC/BIC

  • \(\text{penalización}\): cuando sube \(p\), entonces suba el \(AIC/BIC\)

El objetivo será probar distintos modelos y quedarnos con el que tenga AIC/BIC más pequeño

BIC y AIC

Es habitual que la calidad del modelo venga cuantificada por \(2\ell(\text{modelo})\), donde \(\ell(\text{modelo})\) es lo que se conoce como log-verosimilitud del modelo

\[\ell(\text{modelo}) = \log \left(P \left(\text{datos} | \left(\text{parámetros}, \text{modelo} \right) \right) \right)\]

es decir, suponiendo que el modelo fuese correcto y los parámetros valiesen la estimación obtenida, ¿cómo de probable es que mis datos hayan sido los que han sido?

Así nuestros criterios quedan como

\[AIC/BIC_{modelo} = -2\ell(\text{modelo}) + \underbrace{(p+2) * \text{penalización}}_{\text{Complejidad}}\]

BIC y AIC

La diferencia entre ambos está en la penalización usada

\[BIC(modelo_p) = -2\ell(\text{modelo}) + \underbrace{(p+2) * \log(n)}_{\text{Complejidad}}\]

\[AIC (modelo_p) = -2\ell(\text{modelo}) + \underbrace{(p+2) * 2}_{\text{Complejidad}}\]

El criterio BIC es más agresivo seleccionando variables ya que, si \(n\) crece, \(\log(n) >> 2\) (penaliza más el sobreajuste) y, además, depende del tamaño muestral

BIC y AIC

Vamos a ver un pequeño ejemplo en R con wine.csv, aplicando un modelo saturado (habiendo hecho un filtro inicial por VIF), un modelo solo con AGST, otro con AGST + HarvestRain y otro con AGST + HarvestRain + Year, y calcular su BIC/AIC con las funciones BIC() y AIC()

datos <- read_csv(file = "./datos/wine.csv")
ajuste_saturado <- lm(data = datos, formula = Price ~ .)
ajuste_1 <- lm(data = datos, formula = Price ~ AGST)
ajuste_2 <- lm(data = datos, formula = Price ~ AGST + HarvestRain)
ajuste_3 <- lm(data = datos, formula = Price ~ AGST + HarvestRain + Age)
BIC(ajuste_saturado)
[1] 26.62387
BIC(ajuste_1)
[1] 45.00735
BIC(ajuste_2)
[1] 33.04011
BIC(ajuste_3)
[1] 26.44583

BIC y AIC

Vamos a ver un pequeño ejemplo en R con wine.csv, aplicando un modelo saturado (habiendo hecho un filtro inicial por VIF), un modelo solo con AGST, otro con AGST + HarvestRain y otro con AGST + HarvestRain + Year, y calcular su BIC/AIC con las funciones BIC() y AIC()

AIC(ajuste_saturado)
[1] 17.55302
AIC(ajuste_1)
[1] 41.11983
AIC(ajuste_2)
[1] 27.85676
AIC(ajuste_3)
[1] 19.96665

BIC y AIC

Podemos usar compare_performance() (\(AICc = AIC +\frac{2(p+2)^{2} + 2(p+2)}{n-(p+2)-1}\) es una versión corregida para tamaños muestrales pequeños)

performance::compare_performance(ajuste_saturado, ajuste_1, ajuste_2, ajuste_3)
# Comparison of Model Performance Indices

Name            | Model | AIC (weights) | AICc (weights) | BIC (weights) |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------------------------------------------------
ajuste_saturado |    lm |  17.6 (0.766) |   23.4 (0.415) |  26.6 (0.469) | 0.828 |     0.787 | 0.258 | 0.293
ajuste_1        |    lm |  41.1 (<.001) |   42.2 (<.001) |  45.0 (<.001) | 0.446 |     0.423 | 0.464 | 0.482
ajuste_2        |    lm |  27.9 (0.004) |   29.7 (0.018) |  33.0 (0.019) | 0.685 |     0.659 | 0.350 | 0.371
ajuste_3        |    lm |  20.0 (0.229) |   22.8 (0.567) |  26.4 (0.512) | 0.782 |     0.753 | 0.291 | 0.315

En ambos criterios la conclusión es que de los 4 modelos, el “mejor” es el 3º ya que AIC/BIC más bajos aunque tenga un \(R_{adj}^2\) menor: la mejora que produce meter todas las variables no es suficiente para lo que se complica el modelo.

BIC y AIC

En realidad lo correcto sería antes chequear si podemos eliminar alguna variable ya de manera preliminar usando el VIF

check_collinearity(ajuste_saturado)
# Check for Multicollinearity

Low Correlation

        Term  VIF      VIF 95% CI Increased SE Tolerance Tolerance 95% CI
  WinterRain 1.26 [ 1.06,   2.13]         1.12      0.80     [0.47, 0.95]
        AGST 1.26 [ 1.06,   2.13]         1.12      0.79     [0.47, 0.94]
 HarvestRain 1.13 [ 1.01,   2.48]         1.06      0.88     [0.40, 0.99]

High Correlation

      Term   VIF      VIF 95% CI Increased SE Tolerance Tolerance 95% CI
      Year 69.78 [42.40, 115.27]         8.35      0.01     [0.01, 0.02]
 FrancePop 70.07 [42.58, 115.75]         8.37      0.01     [0.01, 0.02]

BIC y AIC

Según el VIF debemos ya eliminar de antemano Year y FrancePop así que lo hacemos y volvemos calcularlo para las restantes

datos_VIF <- datos |> select(-Year, -FrancePop)
ajuste_saturado <- lm(data = datos_VIF, formula = Price ~ .)
check_collinearity(ajuste_saturado)
# Check for Multicollinearity

Low Correlation

        Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
  WinterRain 1.23 [1.04, 2.35]         1.11      0.82     [0.42, 0.96]
        AGST 1.25 [1.05, 2.33]         1.12      0.80     [0.43, 0.96]
 HarvestRain 1.10 [1.00, 3.70]         1.05      0.91     [0.27, 1.00]
         Age 1.11 [1.00, 3.55]         1.05      0.90     [0.28, 1.00]

BIC y AIC

Tras este primer filtro grosero volvemos a implementar los 3 modelos

ajuste_saturado <- lm(data = datos_VIF, formula = Price ~ .)
ajuste_1 <- lm(data = datos_VIF, formula = Price ~ AGST)
ajuste_2 <- lm(data = datos_VIF, formula = Price ~ AGST + HarvestRain)
ajuste_3 <- lm(data = datos_VIF, formula = Price ~ AGST + HarvestRain + Age)
compare_performance(ajuste_saturado, ajuste_1, ajuste_2, ajuste_3)
# Comparison of Model Performance Indices

Name            | Model | AIC (weights) | AICc (weights) | BIC (weights) |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------------------------------------------------
ajuste_saturado |    lm |  15.6 (0.897) |   19.8 (0.815) |  23.4 (0.818) | 0.828 |     0.796 | 0.259 | 0.286
ajuste_1        |    lm |  41.1 (<.001) |   42.2 (<.001) |  45.0 (<.001) | 0.446 |     0.423 | 0.464 | 0.482
ajuste_2        |    lm |  27.9 (0.002) |   29.7 (0.006) |  33.0 (0.006) | 0.685 |     0.659 | 0.350 | 0.371
ajuste_3        |    lm |  20.0 (0.101) |   22.8 (0.179) |  26.4 (0.175) | 0.782 |     0.753 | 0.291 | 0.315

Ahora nos indica que el mejor modelo es el saturado (pero que ya no tiene 6 variables sino 4): la mejora de bondad de ajuste compensa por incrementar una sola variable (entre ajuste3 y saturado)

BIC y AIC

Fíjate que además al haber quitado dos variables con alta dependenica lineal del resto las hipótesis se cumplen (antes no)

check_model(ajuste_saturado)

BIC y AIC

Vamos a repetirlo con un ejemplo simulado \(Y = 0.01 + 1.5*X_1 -1.5X_2 + \varepsilon\), donde \(X_1,X_2 \sim N(0,1)\) y \(\varepsilon \sim N(0, 2)\), y al que vamos añadiendo predictoras basura \(X_{2+j} \sim N(0, 2)\) con \(\beta_{2+j} = 10^{-8}\), probando cuando \(j=1\), \(j = 5\), \(j = 10\), \(j=25\), \(j=50\) y \(j = 100\). El mejor modelo es el que solo tiene las dos predictoras con un efecto real sobre y.

Código
n <- 200
p <- c(1, 5, 10, 25, 50, 100)

set.seed(12345)
x_1 <- rnorm(n)
x_2 <- rnorm(n)
eps <- rnorm(n, mean = 0, sd = 2)
y <- 0.01 + 1.5*x_1 - 1.5*x_2 + eps
ajuste_0 <- lm(data = tibble("y" = y, "x_1" = x_1, "x_2" = x_2), y ~ .)

X <- matrix(rnorm(n * max(p), mean = 0, sd = 2), nrow = n)
ajuste <- list()
for (j in 1:length(p)) {
  
  ynew <- y + 1e-8 * sum(X[, 1:p[j]])
  datos <- tibble("y" = ynew, "x_1" = x_1, "x_2" = x_2, X[, 1:p[j]])
  ajuste[[j]] <- lm(data = datos, y ~ .)
}
compare_performance(ajuste_0, ajuste[[1]], ajuste[[2]],
                    ajuste[[3]], ajuste[[4]], ajuste[[5]])
# Comparison of Model Performance Indices

Name                         | Model | AIC (weights) | AICc (weights) |  BIC (weights) |    R2 | R2 (adj.) |  RMSE | Sigma
--------------------------------------------------------------------------------------------------------------------------
ajuste_0                     |    lm | 847.1 (0.052) |  847.3 (0.066) |  860.3 (0.407) | 0.531 |     0.527 | 1.972 | 1.987
c(\[[\, \ajuste\, \1\) |    lm | 843.1 (0.393) |  843.4 (0.474) |  859.6 (0.592) | 0.545 |     0.538 | 1.942 | 1.962
c(\[[\, \ajuste\, \2\) |    lm | 842.6 (0.509) |  843.5 (0.445) |  872.3 (0.001) | 0.564 |     0.548 | 1.901 | 1.940
c(\[[\, \ajuste\, \3\) |    lm | 848.0 (0.034) |  850.3 (0.015) |  894.2 (<.001) | 0.574 |     0.547 | 1.880 | 1.944
c(\[[\, \ajuste\, \4\) |    lm | 850.0 (0.013) |  860.2 (<.001) |  945.6 (<.001) | 0.630 |     0.572 | 1.752 | 1.890
c(\[[\, \ajuste\, \5\) |    lm | 865.4 (<.001) |  906.4 (<.001) | 1043.5 (<.001) | 0.689 |     0.578 | 1.607 | 1.875

Sobreajuste

Algo importante a tener en cuenta es que, aunque ambos criterios nos ayudan a seleccionar modelos, ambos funcionan de manera aceptable bajo la hipótesis de que \(n >> p +2\): en caso contrario, si \(n\) se acerca a \(p +2\), el sobreajuste seguirá produciéndose

Dado que la penalización es más grande en el BIC, el criterio BIC nos garantiza una más temprana detección del sobreajuste

Sobreajuste

Lo anterior se puede ilustrar calculando BIC e AIC del anterior estudio de simulación de \(R^2\) vs \(R_{adj}^2\): fíjate como el BIC tarda más en bajar, tal que \(BIC(p = 195) > BIC(p=2)\) (nos haría quedarnos con el modelo sin sobreajustar) mientras que \(AIC(p = 156) < AIC(p=2)\) (y para todos los que van detrás), por lo que acabaríamos eligiendo el modelo más sobreajustado.

Consistencia

Otra enorme ventaja del BIC es que se ha demostrado matemáticamente que es consistente: si el tamaño muestral fuese lo suficientemente grande, el BIC garantiza elegir el modelo correcto que genera los datos. Matemáticamente se cumple que, si tenemos una serie de modelos \(M_1, \ldots, M_m\), y el modelo real que genera los datos \(M_0\) (que pretendemos estimar), entonces

\[P\left[\arg\min_{k=0,\ldots,m}\text{BIC}(\widehat{M}_k)=0\right]\to 1, \quad n \to \infty\]

Esto solo suponiendo que el modelo subyacente sea lineal claro… Esto no sucede con el AIC (para vuestro yo del futuro: en https://doi.org/10.2307/2290328. se prueba que el AIC es equivalente a usar validación cruzada leave-one-out, el cual es inconsistente)

stepAIC

Y aquí nos puede nacer una duda: si tengo muchos predictores, ¿tengo que calcular el BIC/AIC de todas las combinaciones posibles?

Sí, pero lo hará por nosotros MASS::stepAIC(), donde en el parámetro k = ... le indicamos la penalización (si k = 2 es AIC y si k = log(n) es el BIC)

datos <- read_csv(file = "./datos/wine.csv")
ajuste_saturado <- lm(data = datos, Price ~ .)
MASS::stepAIC(ajuste_saturado, k = log(nrow(datos)))

stepAIC

datos <- read_csv(file = "./datos/wine.csv")
ajuste_saturado <- lm(data = datos, Price ~ .)
MASS::stepAIC(ajuste_saturado, k = log(nrow(datos)))
Start:  AIC=-53.29
Price ~ Year + WinterRain + AGST + HarvestRain + Age + FrancePop


Step:  AIC=-53.29
Price ~ Year + WinterRain + AGST + HarvestRain + FrancePop

              Df Sum of Sq    RSS     AIC
- FrancePop    1    0.0026 1.8058 -56.551
- Year         1    0.0048 1.8080 -56.519
<none>                     1.8032 -53.295
- WinterRain   1    0.4585 2.2617 -50.473
- HarvestRain  1    1.8063 3.6095 -37.852
- AGST         1    3.3756 5.1788 -28.105

Step:  AIC=-56.55
Price ~ Year + WinterRain + AGST + HarvestRain

              Df Sum of Sq    RSS     AIC
<none>                     1.8058 -56.551
- WinterRain   1    0.4809 2.2867 -53.473
- Year         1    0.9089 2.7147 -48.840
- HarvestRain  1    1.8760 3.6818 -40.612
- AGST         1    3.4428 5.2486 -31.039

Call:
lm(formula = Price ~ Year + WinterRain + AGST + HarvestRain, 
    data = datos)

Coefficients:
(Intercept)         Year   WinterRain         AGST  HarvestRain  
  43.639042    -0.023848     0.001167     0.616392    -0.003861  

stepAIC

Lo que obtendremos será una selección secuencial de modelos, de manera que el irá probando las combinaciones, nos muestra el ajuste y el valor del BIC/AIC y se parará cuando el AIC/BIC no mejore.

datos <- read_csv(file = "./datos/wine.csv")
ajuste_saturado <- lm(data = datos, Price ~ .)
MASS::stepAIC(ajuste_saturado, k = log(nrow(datos)))

stepAIC

MASS::stepAIC(ajuste_saturado, k = log(nrow(datos)), direction = "forward")
Start:  AIC=-53.29
Price ~ Year + WinterRain + AGST + HarvestRain + Age + FrancePop

Call:
lm(formula = Price ~ Year + WinterRain + AGST + HarvestRain + 
    Age + FrancePop, data = datos)

Coefficients:
(Intercept)         Year   WinterRain         AGST  HarvestRain          Age  
  2.496e+01   -1.377e-02    1.153e-03    6.144e-01   -3.837e-03           NA  
  FrancePop  
 -2.213e-05  

El argumento direction = ... puede tomar los valores "both", "forward" y "backward" (por defecto) que nos determina la dirección de búsqueda:

  • direction = "forward": empieza con el modelo proporcionado y va añadiendo predictoras haciendo modelos cada vez más complejos.

stepAIC

MASS::stepAIC(ajuste_saturado, k = log(nrow(datos)), direction = "forward")
Start:  AIC=-53.29
Price ~ Year + WinterRain + AGST + HarvestRain + Age + FrancePop

Call:
lm(formula = Price ~ Year + WinterRain + AGST + HarvestRain + 
    Age + FrancePop, data = datos)

Coefficients:
(Intercept)         Year   WinterRain         AGST  HarvestRain          Age  
  2.496e+01   -1.377e-02    1.153e-03    6.144e-01   -3.837e-03           NA  
  FrancePop  
 -2.213e-05  

Fíjate que ahora como le hemos dado el modelo saturado y direction = "forward", entonces no hace nada. Para controlar esto podemos incluir la variable scope = list(lower = mod_easy, upper = mod_complex), donde podemos pasarle dos modelosp: los modelos más complejos y más sencillo pentre los que queremos que se mueva en la búsqueda.

stepAIC

MASS::stepAIC(ajuste_saturado, k = log(nrow(datos)), direction = "backward")
  • direction = "backward": empieza con el modelo proporcionado y va eliminando predictoras haciendo modelos cada vez más sencillos.
MASS::stepAIC(ajuste_saturado, k = log(nrow(datos)))
  • direction = "both": empieza con el modelo proporcionado y va añadiendo y eliminando predictoras según sea más conveniente (pudiendo añadir/eliminar una variable que previamente e había eliminado/añadido)

stepAIC

MASS::stepAIC(ajuste_saturado, k = log(nrow(datos)), direction = "forward")
Start:  AIC=-53.29
Price ~ Year + WinterRain + AGST + HarvestRain + Age + FrancePop

Call:
lm(formula = Price ~ Year + WinterRain + AGST + HarvestRain + 
    Age + FrancePop, data = datos)

Coefficients:
(Intercept)         Year   WinterRain         AGST  HarvestRain          Age  
  2.496e+01   -1.377e-02    1.153e-03    6.144e-01   -3.837e-03           NA  
  FrancePop  
 -2.213e-05  

Aviso: en realidad MASS::stepAIC() no calcula exactamente el mismo BIC/AIC que las funciones BIC() e AIC(), sino que les suma una constante \(n(\log(2\pi) + 1) + \log(n)\) para BIC y \(n(\log(2\pi) + 1) + 2\) para AIC

Dado que es una constante nos da igual para comparar modelos, pero hace que NO podamos comparar salidas de BIC() y AIC() con salidas de MASS::stepAIC().

Clase 15: casos prácticos

Casos prácticos: datos de viviendas de Boston y seatpos dataset

Casos prácticos

Realiza todo el ajuste completo multivariante con los datasets:

  • seatpos del paquete {faraway}: datos de 38 conductores donde el objetivo es predecir hipcenter, la posición del asiento del conductor, en función de distintas variables (ver ? faraway::seatpos)

  • Boston del paquete {MASS}: datos de 560 suburbios de Boston, en el que se han medido 14 variables en cada uno, con el objetivo de predecir medv el precio mediano de inmuebles (en millones de dolares), en función de variables estructurales (rm y age), variables de vecindario (crim, zn, indus, chas, tax, ptratio, black y lstat), variables de accesibilidad (dis y rad) y variables de calidad del aire (nox)

Clase 16: variables cualitativas

¿Cómo introducir predictoras cualitativas?

Cosas que faltan por añadir

  • ejemplo seatpos

  • ejemplo one-hot-encoding

  • ejemplo boston

Clase n:

Regresión penalizada

Hagamos un resumen: ¿cuál era el objetivo de la regresión ordinaria?

Asumiendo que

\[E \left[Y | \boldsymbol{X} = x \right] = \beta_0 + \displaystyle \sum_{j=1}^{p} \beta_j X_j, \quad \widehat{Y} = \widehat{\beta_0} + \displaystyle \sum_{j=1}^{p} \widehat{\beta_j} X_j, \quad E \left[ \varepsilon | \left( X_1 = x_1, \ldots, X_p = x_p \right) \right] = 0 \]

El objetivo ha sido obtener la estimación de los \(\widehat{\beta}\) tal que minimicemos la varianza residual o suma de errores al cuadrado \(SSE\)

\[SSE \left(\boldsymbol{\beta} \right) = \displaystyle \sum_{i=1}^{n} \widehat{\varepsilon}_{i}^{2} = \displaystyle \sum_{i=1}^{n} \left(Y_i - \beta_0 - \beta_1 X_{i1} - \ldots - \beta_p X_{ip} \right)^2 = \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)^{T}\left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right) \]

Tras derivar e igualar a cero obteníamos \(\widehat{\boldsymbol{\beta}}=\left(\mathbf{X}^{T}\mathbf{X} \right)^{-1}\mathbf{X}^{T}\mathbf{Y}\)

Regresión penalizada

El estimador \(\widehat{\boldsymbol{\beta}}\) es aquel que nos garantiza

\[\widehat{\boldsymbol{\beta}} = \arg \min_{\beta \in \mathbb{R}^{p+1}} SSE \left(\boldsymbol{\beta} \right)\]

También tenemos garantizado que si se cumple las hipótesis tenemos que

\[\widehat{\boldsymbol{\beta}} \sim N_{p+1}\left(\boldsymbol{\beta}, \sigma_{\varepsilon}^2 \left( \mathbf{X}^{T}\mathbf{X} \right)^{-1}\right)\]

lo que implica que \(E \left[ \widehat{\boldsymbol{\beta}}\right] = \boldsymbol{\beta}\), es decir, tenemos un estimador insesgado.

El sesgo en una estimación es un error sistemático en la misma. Pero hay otra componente importante que hasta ahora solo hemos mencionado pero que no hemos tenido en cuenta: la varianza de la estimación

Balance sesgo-varianza

Si descomponemos el error cuadrático medio de una estimación obtenemos que

\[\begin{eqnarray}MSE \left(\widehat{Y} \right) &=& E \left[ \left( \widehat{Y} - Y \right)^2 \right] = E \left[ \left( \left(\widehat{Y} - E \left[ \widehat{Y}\right] \right) - \left(E \left[ \widehat{Y}\right] - Y \right)\right)^2 \right] \nonumber \\ &=& E \left[ \left(\widehat{Y} - E \left[ \widehat{Y}\right] \right)^2 \right] +\left( E \left[ \widehat{Y}\right] - Y \right)^2 + 2*cosas \nonumber \\ &=& \underbrace{\left( E \left[ \widehat{Y}\right] - Y \right)^2}_{Sesgo^2}+\underbrace{Var \left[ \widehat{Y} \right]}_{Varianza} + ruido\end{eqnarray}\]

El eror cuadrático medio podemos descomponerlo en una suma de sesgo (al cuadrado), varianza y ruido. El último es intrínseco a los datos, por lo que la metodología de estimación puede bajar el sesgo (subiendo la varianza) o bajar la varianza (subiendo el sesgo).

Balance sesgo-varianza

  • Bajoajuste (underfitting): modelos muy simples proporcionan un sesgo muy grande y poca varianza ya que la predicción siempre será muy parecida (errores altos en train).

  • Sobreajuste (overfitting): modelos muy complicados proporcionan un sesgo bajo pero al ser tan complejas proporcionarán una mayor varianza para cada intento (errores altos en test).

Balance sesgo-varianza

Si lo visto en las estimaciones lo trasladamos a los coeficientes estimados tenemos

\[\begin{eqnarray}MSE \left(\widehat{\beta}_j \right) = \underbrace{\left( E \left[ \widehat{\beta}_j \right] - \beta_j \right)^2}_{Sesgo^2}+\underbrace{Var \left[ \widehat{\beta}_j \right]}_{Varianza} + ruido \end{eqnarray} \]

En este caso además sabemos que \(Var \left[ \widehat{\beta}_j \right] = SE \left(\widehat{\beta}_j\right)^2 = \sigma_{\varepsilon}^2 v_j\), con \(v_j\) el elemento \(j\)-ésimo de la diagonal de \(\left( \mathbf{X}^{T}\mathbf{X} \right)^{-1}\).

La forma más sencilla de reducir drásticamente la varianza residual es aumentando el número de parámetros, es decir sobreajustando. Dicho de otra forma: para conseguir simplificar modelos vamos a reducir la varianza de las estimaciones, y para ello lo que haremos será introducir un pequeño sesgo

Regresión penalizada

\[\begin{eqnarray}MSE \left(\widehat{\beta}_j \right) &=& \underbrace{\left( E \left[ \widehat{\beta}_j \right] - \beta_j \right)^2}_{Sesgo^2}+\underbrace{Var \left[ \widehat{\beta}_j \right]}_{Varianza} + ruido = cte\end{eqnarray} \]

La idea será sesgar los coeficientes de tal manera que se busque obtener coeficientes que, en caso de no ser tan importantes como otros, tiendan a cero, un sesgo hacia lo que se conoce como sparsity.

Para ello haremos uso de lo que se conoce como shrinkage methods, y en particular hablaremos de regresión penalizada

Regresión penalizada

Si en el caso de la regresión lineal ordinaria teníamos

\[\widehat{\boldsymbol{\beta}} = \arg \min_{\beta \in \mathbb{R}^{p+1}} SSE \left(\boldsymbol{\beta} \right)\]

Ahora incluiremos una penalización \(\lambda \geq 0\)

\[\widehat{\boldsymbol{\beta}}_{\lambda} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left\| \beta_{-1} \right\| \right), \quad \lambda \geq 0\]

donde \(\left\| \beta \right\|\) será alguna norma del vector de coeficientes (excluyendo el intercepto) (distancia del vector a 0) y \(\lambda\) deberemos determinar su valor óptimo.

Por tanto ahora el objetivo no será solo minimizar los errores al cuadrado sino también minimizar la penalización, y la única forma de hacerlo (dado que una norma siempre devuelve algo positiva) es que los coeficientes sean lo más pequeño posibles

Penalización

\[\widehat{\boldsymbol{\beta}}_{\lambda} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left\| \beta_{-1} \right\| \right), \quad \lambda \geq 0\]

¿Qué sucede cuando \(\lambda\) crece o decrece?

  • Si \(\lambda \to 0\) –> estamos en el caso de mínimos cuadrados ordinarios (no hay penalización)

  • Si \(\lambda \to \infty\) –> \(\beta_j \to 0\) para todo \(j \geq 0\) (el modelo tiende a desaparecer)

Familia de normas

\[\widehat{\boldsymbol{\beta}}_{\lambda} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left\| \beta_{-1} \right\| \right), \quad \lambda \geq 0\]

Una de las familias de normas más habituales son las normas \(\ell^p\)

\[\left\| x \right\|_{p} = \left\| \left(x_1, \ldots, x_k \right)\right\|_{p} = \left( \sum_{j=1}^{k} \left| x_{j} \right|^{p} \right)^{1/p}, \quad x \in \mathbb{R}^{k}\]

  • Si \(p=1\) –> \(\left\| x \right\|_{1} = \sum_{j=1}^{k} \left| x_{j} \right|\) se conoce como norma taxicab o de Manhattan (la métrica que usaríamos para recorrer calles en un mapa)

  • Si \(p=2\) –> \(\left\| x \right\|_{2} = \left(\sum_{j=1}^{k} x_{j}^{2} \right)^{1/2}\) se conoce como norma Euclídea (la habitual)

  • Si \(p=\infty\) –> \(\left\| x \right\|_{\infty} = \lim_{p \to \infty} \left(\sum_{j=1}^{k} \left| x_{j} \right|^{p} \right)^{1/p} = \displaystyle \max_{j} \left| x_{j} \right|\) se conoce como norma infinita o de Chebyshev (movimientos de rey en ajedrez)

Regresión elastic-net

Dado que las más usadas son la norma Manhattan y la Euclídea, podemos generalizar la regresión penalizada anterior combinando ambas penalizaciones

\[\widehat{\boldsymbol{\beta}}_{\lambda, \alpha} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left(\alpha \left\| \beta_{-1} \right\|_1 + (1-\alpha) \left\| \beta_{-1} \right\|_{2}^{2} \right)\right), \quad \lambda \geq 0,~0 \leq \alpha \leq 1\]

La regresión elastic-net es una regresión penalizada que combina una penalización \(\lambda\) (que nos cuantifica lo agresivo que queremos hacer tender los coeficientes a cero) y una proporción \(\alpha\) (que nos permite combinar ambas penalizaciones)

Importante

Dado que la penalización depende de la norma de los coeficientes, su magnitud será importante, por lo que deberemos estandarizar antes los predictores.

Ridge y lasso

\[\widehat{\boldsymbol{\beta}}_{\lambda, \alpha} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left(\alpha \left\| \beta_{-1} \right\|_1 + (1-\alpha) \left\| \beta_{-1} \right\|_{2}^{2} \right)\right), \quad \lambda \geq 0,~0 \leq \alpha \leq 1\]

Hay dos escenarios de penalización principalmente conocidos:

  • Regresión ridge: cuando \(\alpha = 0\) (penalización cuadrática)

\[\widehat{\boldsymbol{\beta}}_{\lambda, ridge} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left\| \beta_{-1} \right\|_{2}^{2} \right), \quad \lambda \geq 0\]

  • Regresión LASSO: cuando \(\alpha = 1\) (penalización absoluta)

\[\widehat{\boldsymbol{\beta}}_{\lambda, LASSO} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left\| \beta_{-1} \right\|_1 \right), \quad \lambda \geq 0\]

Ridge y lasso

  • Regresión ridge: cuando \(\alpha = 0\) (penalización cuadrática)

\[\widehat{\boldsymbol{\beta}}_{\lambda, ridge} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left\| \beta_{-1} \right\|_{2}^{2} \right), \quad \lambda \geq 0\]

Una de las ventajas de la regresión ridge es que podemos calcular su expresión explícita (algo que no sucede en la regresión lasso, solo se puede obtener mediante simulación numérica)

\[\widehat{\boldsymbol{\beta}}_{\lambda, ridge} = \left(\mathbf{X}^{T}\mathbf{X} + \lambda \mathbf{I}_{p \times p} \right)^{-1}\mathbf{X}^{T}\mathbf{Y}\] donde \(\mathbf{I}_{p \times p}\) es la matriz identidad.

Ridge y lasso

  • Regresión LASSO: cuando \(\alpha = 1\) (penalización absoluta)

\[\widehat{\boldsymbol{\beta}}_{\lambda, LASSO} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left\| \beta_{-1} \right\|_1 \right), \quad \lambda \geq 0\]

Una de las ventajas de la regresión LASSO es que nos realiza una selección de variables ya hace tender más rápido los coeficientes a cero.

 

¿Por qué?

Regresión elastic-net

\[\widehat{\boldsymbol{\beta}}_{\lambda, \alpha} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left(\alpha \left\| \beta_{-1} \right\|_1 + (1-\alpha) \left\| \beta_{-1} \right\|_{2}^{2} \right)\right), \quad \lambda \geq 0,~0 \leq \alpha \leq 1\]

La regresión elastic-net puede ser también vista como un problema de minimización de la suma de los residuos al cuadrado con una restricción en los argumentos

\[\widehat{\boldsymbol{\beta}}_{s_\lambda, \alpha} = \arg \min_{\beta \in \mathbb{R}^{p+1}:~\sum_{j=1}^{p} \left(\alpha \left| \beta_{j} \right| + (1-\alpha) \left| \beta_{j} \right|^{2} \right) \leq s_\lambda} SSE \left(\boldsymbol{\beta} \right) , \quad 0 \leq \alpha \leq 1\]

donde \(s_\lambda\) es un valor numérico que no depende de \(\beta\).

Regresión elastic-net

Por ejemplo, si tenemos solo dos predictoras…

Las elipses nos muestran los valores para los que \(SSE \left( \boldsymbol{\beta} \right)\) es constante (lo que queremos minimizar), y las regiones sombreadas representan el espacio de valores posibles debido a la penalización (ver más en https://link.springer.com/book/10.1007/978-1-4614-7138-7)

LASSO

Si nos fijamos en ambas regiones, dado que el LASSO usa una norma absoluta, el espacio de restricciones es más “picudo”, lo que hará más probable que cuando la elipse corte con la zona sombreada lo haga en los vértices (donde uno de los coeficientes es cero) en lugar de en los laterales.

Ver más en https://link.springer.com/book/10.1007/978-1-4614-7138-7

Selección de la penalización

En todo este proceso de regresión penalizada hay un “elefante en la habitación”: ¿cómo seleccionar el \(\lambda\) más óptimo?

Una de las más formas más habituales para ello es lo que se conoce como validación cruzada. La idea se basa en un ideal: poder disponer de muchas muestras tal que, probando en cada una distintos valores de \(\lambda\), consigamos obtener aquel que nos proporcione (en media) el menor error.

¿El problema? Casi nunca disponemos de un gran número de muestras de la misma población, así que una solución es “simularlas” nosotros mismos a partir de la muestra original

Validación cruzada

Supongamos que tenemos una muestra tal que contamos con una variable objetivo a predecir en función de \(p\) predictoras \(\left\lbrace Y_i, \left(X_{i1}, \ldots, X_{ip} \right) \right\rbrace_{i=1,\ldots,n}\), y que pretendemos usar un estimador \(\widehat{f}_{\lambda}\) que depende de un parámetro \(\lambda\), tal que

\[\widehat{Y} =\widehat{f}_{\lambda} \left(X_{1}, \ldots, X_{p} \right)\]

La opción más sencilla es lo que se conoce como validación simple:

Validación simple

\[\widehat{Y} =\widehat{f}_{\lambda} \left(X_{1}, \ldots, X_{p} \right)\] La opción más sencilla es lo que se conoce como validación simple:

  1. Hacemos una partición de train y test. De los datos de train, realizamos una segunda partición que llamamos conjunto de validación.
  1. Entrenamos solo con train el modelo para distintos valores de \(\lambda\)
  1. Evaluamos el error que produce cada \(\widehat{f}_{\lambda}\) en el conjunto de validación (un conjunto del que el modelo no ha aprendido) para elegir el mejor parámetro \(\lambda\).
  1. Proporcionamos una evaluación final en la tercera partición (test), que no ha sido usada ni para la estimación ni para la elección del mejor parámetro.

Validación cruzada

La validación simple tiene un problema: depende muchísimo de la partición aleatoria que se realice (si tenemos mala o buena suerta) y, además, puede provocar un problema de tamaño muestral si \(n\) no es muy grande (ya que necesitamos 3 particiones de los datos).

Para solventarlo existe una alternativa conocida como validación cruzada: en lugar de extraer un subconjunto de train, lo que se hace es que cada observación pueda jugar ambos roles: el rol de entrenamiento y el rol de validación.

Validación cruzada

La más famosa es la conocida como leave-one-out cross-validation (LOOCV):

  1. Hacemos una partición de train (con \(m\) elementos) y test (\(n - m\)).

  2. De los datos de train \(\left\lbrace X_{i} \right\rbrace_{i=1,\ldots,m}\), entrenamos el modelo sin una de las observaciones de entrenamiento \(\boldsymbol{X}_{-1} = \left\lbrace X_{i} \right\rbrace_{i=2,\ldots,m}\), y usamos esa observación \(X_{1}\) para obtener una métrica del error para cada \(\lambda\) (un conjunto de validación de tamaño muestral igual a 1).

  3. Repetimos el proceso para cada una de ellas, considerando como entrenamiento \(\boldsymbol{X}_{-j} = \left\lbrace X_{i} \right\rbrace_{i=1,\ldots, j-1, j+1, \ldots,m}\) y \(X_j\) como conjunto de validación.

  4. En total se entrena el modelo con \(n\) conjuntos train (de \(n-1\) observaciones cada uno) y \(n\) conjuntos de validación (de una observación cada uno), obteniendo al final una media del error en entrenamiento y otra en validación, usando esta última para elegir el mejor parámetro \(\lambda\).

Validación cruzada

La anterior idea puede ser generalizada a lo que se conoce como leave-k-out cross-validation (LOOCV):

  1. Hacemos una partición de train (con \(m\) elementos) y test (\(n - m\)).

  2. De los datos de train \(\left\lbrace X_{i} \right\rbrace_{i=1,\ldots,m}\), entrenamos el modelo sin \(k\) observaciones de entrenamiento \(\boldsymbol{X}_{-i_1, \ldots, -i_k}\), y usamos esa \(k\) observaciónes \(\left\lbrace X_{i_1}, \ldots, X_{i_k} \right\rbrace\) para obtener una métrica del error (un conjunto de validación de tamaño muestral igual a \(k\)).

  3. Repetimos el proceso para cada posible forma de combinar \(m\) elementos de \(k\) en \(k\)

  4. En total se entrena el modelo con \(C_{k}^{m}\) conjuntos train (de \(m-k\) observaciones cada uno) y \(C_{k}^{m}\) conjuntos de validación (de \(k\) observaciones cada uno), obteniendo estos últimos para elegir el mejor parámetro \(\lambda\).

 

Aunque esta opción nos da una métrica más robusta del error, si hacemos unas sencillas cuentas vemos que puede no ser factible si \(m\) y/o \(k\) toman valores elevados. Por ejemplo, si \(m = 100\) y \(k = 30\), tendríamos que entrenar \(C_{k}^{m} \simeq 10^{25}\) modelos

Validación cruzada k-folds

Dado que las dos opciones anteriores pueden no servirnos existe una tercera vía, la conocida como validación cruzada k-folds:

Validación cruzada k-folds

Dado que las dos opciones anteriores pueden no servirnos existe una tercera vía, la conocida como validación cruzada k-folds:

  1. Hacemos una partición de train (con \(m\) elementos) y test (\(n - m\)).

  2. De los datos de train \(\left\lbrace X_{i} \right\rbrace_{i=1,\ldots,m}\), dividimos dicho conjunto en \(k\) particiones (\(k\) folds) que denotaremos como \(\boldsymbol{X}_{fold_1}, \ldots, \boldsymbol{X}_{fold_k}\) (con \(m/k\) observaciones cada una aproximadamente).

  3. Entrenamos \(k\) modelos, de manera que en la iteración \(j\), usaremos las \(m-m/k\) observaciones de las particiones \(\boldsymbol{X}_{fold_1}, \ldots, \boldsymbol{X}_{fold_{j-1}}, \boldsymbol{X}_{fold_{j+1}}, \ldots, \boldsymbol{X}_{fold_{k}}\) para entrenar y las \(m/k\) observaciones de \(\boldsymbol{X}_{fold_{j}}\) para validar.

  4. De esta manera cada observación jugará \(k-1\) veces el rol de entrenamiento y una sola vez el rol de validación, teniendo \(k\) validaciones cuyos resultados promediaremos para elegir el mejor parámetro \(\lambda\).

Selección de la penalización

Así llamaremos \(\widehat{\lambda}_{k-CV}\) a la penalización óptima elegida a través de validación k-fold

\[\widehat{\lambda}_{k-CV} = \arg \min_{\lambda \geq 0} CV_k(\lambda),\quad CV_k(\lambda)=\sum_{j=1}^k SSE \left(\widehat{\boldsymbol{\beta}_{(j)} } \right)\] donde \(SSE \left(\widehat{\boldsymbol{\beta}}_{(k)} \right)\) es la suma de residuos al cuadrado del \(j\)-ésimo slot de validación (con los parámetros obtenidos en el entrenamiento sin ese slot).

Selección de la penalización

\[\widehat{\lambda}_{k-CV} = \arg \min_{\lambda \geq 0} CV_k(\lambda),\quad CV_k(\lambda)=\sum_{j=1}^k SSE \left(\widehat{\boldsymbol{\beta}_{(j)} } \right)\]

Con el objetivo de simplificar aún más el modelo existe una variante de penalización \(\widehat{\lambda}_{k-1SE}\) conocida como «one standard error rule»: la idea es elegir el modelo más sencillo de los modelos aceptanles, la penalización más grande posible (es decir, el modelo más sencillo posible) acercándose lo máximo que pueda a \(\widehat{\lambda}_{k-CV}\). La idea consiste en estimar \(\widehat{\lambda}_{k-CV}\) y calcular su correspondiente \(CV_k(\widehat{\lambda}_{k-CV})\) para cada partición de validación(que es en sí una variable aleatoria), para posteriormente estimar su variabilidad (\(SE\) igual a desviación típica dividida entre \(\sqrt{n}\)) entre todas las particiones de validación, tal que

\[\widehat{\lambda}_{k-1SE}= \max \left\lbrace \lambda \geq 0: CV_k(\lambda) \in \left(CV_k(\widehat{\lambda}_{k-CV})\pm\widehat{\mathrm{SE}}\left(CV_k(\widehat{\lambda}_{k-CV})\right)\right)\right\rbrace\]

Regresión penalizada en R

Vamos a realizar un pequeño ejemplo con el famoso dataset iris, al que primero lo preprocesaremos convenientemente para tener solo variables numéricas (la variable objetivo será Sepal.Length)

Código
datos <- as_tibble(iris)
datos_preproc <-
  datos |> 
  fastDummies::dummy_cols(select_columns = "Species", remove_first_dummy = TRUE) |> 
  select(-Species)

# Ajuste Saturado
ajuste_saturado <- lm(data = datos_preproc, Sepal.Length ~ .)
ajuste_saturado |> summary()

Call:
lm(formula = Sepal.Length ~ ., data = datos_preproc)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.79424 -0.21874  0.00899  0.20255  0.73103 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         2.17127    0.27979   7.760 1.43e-12 ***
Sepal.Width         0.49589    0.08607   5.761 4.87e-08 ***
Petal.Length        0.82924    0.06853  12.101  < 2e-16 ***
Petal.Width        -0.31516    0.15120  -2.084  0.03889 *  
Species_versicolor -0.72356    0.24017  -3.013  0.00306 ** 
Species_virginica  -1.02350    0.33373  -3.067  0.00258 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3068 on 144 degrees of freedom
Multiple R-squared:  0.8673,    Adjusted R-squared:  0.8627 
F-statistic: 188.3 on 5 and 144 DF,  p-value: < 2.2e-16

Regresión penalizada en R

Tras realizar el ajuste saturado (vamos a hacer una prueba sencilla, no hacemos colinealidad ni diagnosis, solo ajuste), seleccionamos modelos con AIC y BIC

ajuste_AIC <- MASS::stepAIC(ajuste_saturado, direction = "both", k = 2)
ajuste_BIC <- MASS::stepAIC(ajuste_saturado, direction = "both", k = log(nrow(datos_preproc)))
ajuste_AIC |> summary()
ajuste_BIC |> summary()

Regresión ridge en R

\[\widehat{\boldsymbol{\beta}}_{\lambda, \alpha} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left(\alpha \left\| \beta_{-1} \right\|_1 + (1-\alpha) \left\| \beta_{-1} \right\|_{2}^{2} \right)\right), \quad \lambda \geq 0,~0 \leq \alpha \leq 1\]

¿Cómo realizar una regresión penalizada en R?. Para ello vamos a usar el completísimo paquete {glmnet}, cuya documentación puedes consultar en https://glmnet.stanford.edu/articles/glmnet.html#linear-regression-family-gaussian-default

Para realizar una ridge regresión, es decir \(\alpha = 0\), basta con hacer uso de la función glmnet(), dándole en x matriz de variables (solo numéricas en formato matriz), en y la variable objetivo y alpha = 0

library(glmnet)

y <- datos_preproc |> select(Sepal.Length) |> as.matrix()
x <- datos_preproc |> select(-Sepal.Length) |> as.matrix()

ajuste_ridge <- glmnet(x = x, y = y, alpha = 0)

Regresión ridge en R

¿Qué tenemos guardado?

  • $lambda: tenemos los valores de \(\lambda\) probados (por defecto, hay un argumento nlambda en glmnet() que vale 100, y la secuencia de valores los decide en función de otro argumento lambda.min.ratio que por defecto vale 0.01)
ajuste_ridge$lambda
  [1] 719.45951736 655.54471675 597.30793087 544.24474054 495.89553779
  [6] 451.84154496 411.70118743 375.12678864 341.80155863 311.43685020
 [11] 283.76965879 258.56034440 235.59055603 214.66134035 195.59141850
 [16] 178.21561595 162.38343182 147.95773528 134.81357787 122.83711117
 [21] 111.92460077 101.98152771  92.92176986  84.66685593  77.14528580
 [26]  70.29191123  64.04737157  58.35757960  53.17325307  48.44948782
 [31]  44.14536885  40.22361596  36.65026077  33.39435262  30.42769037
 [36]  27.72457822  25.26160310  23.01743191  20.97262670  19.10947634
 [41]  17.41184313  15.86502298  14.45561806  13.17142080  12.00130807
 [46]  10.93514494   9.96369680   9.07854944   8.27203613   7.53717124
 [51]   6.86758972   6.25749198   5.70159364   5.19507977   4.73356320
 [56]   4.31304649   3.92988733   3.58076697   3.26266150   2.97281565
 [61]   2.70871890   2.46808378   2.24882602   2.04904651   1.86701487
 [66]   1.70115441   1.55002854   1.41232828   1.28686093   1.17253975
 [71]   1.06837455   0.97346310   0.88698332   0.80818617   0.73638914
 [76]   0.67097035   0.61136320   0.55705138   0.50756447   0.46247384
 [81]   0.42138894   0.38395391   0.34984450   0.31876528   0.29044705
 [86]   0.26464453   0.24113424   0.21971254   0.20019388   0.18240921
 [91]   0.16620448   0.15143933   0.13798588   0.12572759   0.11455830
 [96]   0.10438126   0.09510831   0.08665915   0.07896059   0.07194595

Regresión ridge en R

¿Qué tenemos guardado?

  • $dev.ratio: tenemos lo que se conoce como «deviance» para las penalizaciones probadas, que es una generalización de \(R^2\) para modelos más complejos (en nuestro caso es equivalente)
ajuste_ridge$dev.ratio
  [1] 3.749855e-36 4.649958e-03 5.100797e-03 5.595078e-03 6.136933e-03
  [6] 6.730875e-03 7.381832e-03 8.095183e-03 8.876794e-03 9.733060e-03
 [11] 1.067095e-02 1.169804e-02 1.282258e-02 1.405354e-02 1.540066e-02
 [16] 1.687446e-02 1.848640e-02 2.024882e-02 2.217508e-02 2.427958e-02
 [21] 2.657782e-02 2.908643e-02 3.182326e-02 3.480739e-02 3.805914e-02
 [26] 4.160012e-02 4.545324e-02 4.964263e-02 5.419367e-02 5.913286e-02
 [31] 6.448777e-02 7.028683e-02 7.655922e-02 8.333455e-02 9.064263e-02
 [36] 9.851306e-02 1.069748e-01 1.160558e-01 1.257821e-01 1.361777e-01
 [41] 1.472633e-01 1.590559e-01 1.715678e-01 1.848060e-01 1.987710e-01
 [46] 2.134560e-01 2.288464e-01 2.449191e-01 2.616415e-01 2.789718e-01
 [51] 2.968584e-01 3.152404e-01 3.340476e-01 3.532019e-01 3.726178e-01
 [56] 3.922041e-01 4.118659e-01 4.315057e-01 4.510264e-01 4.703327e-01
 [61] 4.893335e-01 5.079439e-01 5.260868e-01 5.436942e-01 5.607114e-01
 [66] 5.770852e-01 5.927825e-01 6.077783e-01 6.220574e-01 6.356142e-01
 [71] 6.484515e-01 6.605788e-01 6.720118e-01 6.827706e-01 6.928786e-01
 [76] 7.023514e-01 7.112339e-01 7.195478e-01 7.273222e-01 7.345867e-01
 [81] 7.413706e-01 7.477221e-01 7.536297e-01 7.591428e-01 7.642901e-01
 [86] 7.691003e-01 7.736015e-01 7.778445e-01 7.818415e-01 7.856024e-01
 [91] 7.891877e-01 7.926015e-01 7.958499e-01 7.989911e-01 8.020236e-01
 [96] 8.049648e-01 8.078121e-01 8.106137e-01 8.133596e-01 8.160575e-01

Regresión ridge en R

Fíjate que a mayor penalización (modelos más simples, más sparsity), el \(R^2\) disminuye (o lo que es lo mismo, aumenta el SSE)

ggplot(tibble("lambda" = log(ajuste_ridge$lambda),
              "R2" = ajuste_ridge$dev.ratio)) +
  geom_line(aes(x = lambda, y = R2)) +
  theme_minimal() + labs(x = "log-lambda", y = "R2")

Regresión ridge en R

  • $a0: tenemos guardados los interceptos \(\widehat{\beta}_0\) de las regresiones ajustadas para cada penalización (en este caso 100 valores distintos)

  • $beta: tenemos guardados en formato matricial los coeficientes \(\widehat{\beta}_j\) (para \(j \geq 1\)) de las regresiones ajustadas para cada penalización (en este caso una matriz de 5 filas (5 predictoras) y 100 columnas (100 penalizaciones)

dim(ajuste_ridge$beta)
[1]   5 100
length(ajuste_ridge$a0)
[1] 100

Regresión ridge en R

Fíjate que a mayor penalización (modelos más simples, más sparsity), la norma Euclídea \(\left\| \cdot \right\|_2\) de los coeficientes disminuye (ya que penaliza más)

ggplot(tibble("lambda" = log(ajuste_ridge$lambda),
              "norma" = sqrt(colSums(ajuste_ridge$beta^2)))) +
  geom_line(aes(x = lambda, y = norma)) +
  theme_minimal() + labs(x = "log-lambda", y = "norma")

Selección de lambda en R

Para no tener que elegir manualmente el mejor \(\lambda\), la función cv.glmnet() nos permite elegir el mejor \(\widehat{\lambda}_{k-CV}\) (penalización óptima elegida a través de validación k-fold)

  • nfolds: indica el número \(k\) de particiones en la validación cruzada.

  • type.measure: métrica de error usada (por defecto lo que hemos llamado \(SSE\) o \(MSE\))

set.seed(12345)
kfolds_lambda_ridge <- cv.glmnet(x = x, y = y, alpha = 0, nfolds = 20)

Selección de lambda en R

set.seed(12345)
kfolds_lambda_ridge <- cv.glmnet(x = x, y = y, alpha = 0, nfolds = 20)
  • $lambda: los valores \(\lambda\) probados.
  • $lambda.min: el valor que nos proporciona el promedio de errores en validación más pequeño
  • $cvm: valores \(CV_k(\lambda)\)
# Lambda mínimo
kfolds_lambda_ridge$lambda.min
[1] 0.07194595
# Lambda mínimo manualmente
cv_min <- which.min(kfolds_lambda_ridge$cvm)
kfolds_lambda_ridge$lambda[cv_min]
[1] 0.07194595

Selección de lambda en R

Es importante chequear que el mínimo no sea un extremo de la secuencia de \(\lambda\) (ya que sino no sabemos si el mínimo está fuera del rango).

min(kfolds_lambda_ridge$lambda)
[1] 0.07194595
kfolds_lambda_ridge$lambda.min
[1] 0.07194595

Para solventarlo podemos proporcionar manualmente un grid (normalmente logarítmico) de penalizaciones con las que debe probar la función glmnet(..., lambda = ...)

# grid logarítmico
lambda_seq <- 10^seq(log10(max(kfolds_lambda_ridge$lambda)), log10(0.0001), l = 200)

set.seed(12345)
kfolds_lambda_ridge <-
  cv.glmnet(x = x, y = y, alpha = 0, nfolds = 20, lambda = lambda_seq)
kfolds_lambda_ridge$lambda.min
[1] 0.000786847

Selección de lambda en R

  • $lambda.1se: mejor penalización \(\widehat{\lambda}_{k-1SE}\) según «one standard error rule»
kfolds_lambda_ridge$lambda.1se
[1] 0.01604236
kfolds_lambda_ridge

Call:  cv.glmnet(x = x, y = y, lambda = lambda_seq, nfolds = 20, alpha = 0) 

Measure: Mean-Squared Error 

      Lambda Index Measure      SE Nonzero
min 0.000787   174 0.09912 0.01013       5
1se 0.016042   136 0.10881 0.01192       5

Selección de lambda en R

Con la función plot() podemos dibujar de manera automática los distintos valores \(CV_k(\lambda)\), resultado de promediar en los \(k\)-folds de la validación (la incertidumbre se representa con las barras de error). Con abline() podemos pintar las lineas horizontales para \(\widehat{\lambda}_{k-CV}\) y \(\widehat{\lambda}_{k-1SE}\).

plot(kfolds_lambda_ridge)
abline(h = min(kfolds_lambda_ridge$cvm) + c(0, min(kfolds_lambda_ridge$cvsd)))

Selección de lambda en R

Podemos replicarlo en ggplot ya que en $cvlo y $cvup tenemos guardado los valores superiores e inferiores en las validaciones probadas (los extremos de las barras de error)

Código
ggplot() +
  geom_errorbar(data =
                  tibble("log_lambda" = log(kfolds_lambda_ridge$lambda),
                         "cvlo" = kfolds_lambda_ridge$cvlo,
                         "cvup" = kfolds_lambda_ridge$cvup),
                aes(x = log_lambda,
                    ymin = cvlo, ymax = cvup), color = "gray60") + 
        geom_point(data =
                     tibble("log_lambda" = log(kfolds_lambda_ridge$lambda),
                            "MSE" = kfolds_lambda_ridge$cvm),
                   aes(x = log_lambda, y = MSE), color = "red") +
        geom_vline(xintercept = log(kfolds_lambda_ridge$lambda.min),
                   linetype = "dashed") +
        geom_vline(xintercept = log(kfolds_lambda_ridge$lambda.1se),
                   linetype = "dashed") +
  theme_minimal()

Regresión ridge en R

Para aplicar el modelo de regresión ridge con la penalización óptima considerada se puede hacer uso de predict(kfolds_lambda_ridge, ...), indicándole que quieres «predecir» los coeficientes (type = "coefficients") o con coef()

predict(kfolds_lambda_ridge, type = "coefficients",
        s = kfolds_lambda_ridge$lambda.1se)
6 x 1 sparse Matrix of class "dgCMatrix"
                           s1
(Intercept)         2.3047587
Sepal.Width         0.5600840
Petal.Length        0.5496809
Petal.Width        -0.1088164
Species_versicolor -0.1158683
Species_virginica  -0.2110796
coef(kfolds_lambda_ridge)
6 x 1 sparse Matrix of class "dgCMatrix"
                           s1
(Intercept)         2.3047587
Sepal.Width         0.5600840
Petal.Length        0.5496809
Petal.Width        -0.1088164
Species_versicolor -0.1158683
Species_virginica  -0.2110796

Regresión ridge en R

Importante

Los coeficientes calculados son los promedios en validación. Si se quiere un ajuste “exacto” en el dataset original…

ajuste_ridge_CV <- glmnet(x = x, y = y, alpha = 0,
                          lambda = kfolds_lambda_ridge$lambda.1se)
ajuste_ridge_CV$dev.ratio
[1] 0.8519883
ajuste_ridge_CV$a0
      s0 
2.304107 
ajuste_ridge_CV$beta
5 x 1 sparse Matrix of class "dgCMatrix"
                           s0
Sepal.Width         0.5600173
Petal.Length        0.5503804
Petal.Width        -0.1095168
Species_versicolor -0.1170738
Species_virginica  -0.2126737
predict(kfolds_lambda_ridge, type = "coefficients",
        s = kfolds_lambda_ridge$lambda.1se, newx = x)
6 x 1 sparse Matrix of class "dgCMatrix"
                           s1
(Intercept)         2.3047587
Sepal.Width         0.5600840
Petal.Length        0.5496809
Petal.Width        -0.1088164
Species_versicolor -0.1158683
Species_virginica  -0.2110796

Regresión lasso en R

\[\widehat{\boldsymbol{\beta}}_{\lambda, \alpha} = \arg \min_{\beta \in \mathbb{R}^{p+1}} \left( SSE \left(\boldsymbol{\beta} \right) + \lambda \left(\alpha \left\| \beta_{-1} \right\|_1 + (1-\alpha) \left\| \beta_{-1} \right\|_{2}^{2} \right)\right), \quad \lambda \geq 0,~0 \leq \alpha \leq 1\]

La única diferencia para implementar la regresión lasso en R es \(\alpha = 1\) (con la misma función glmnet())

ajuste_lasso <- glmnet(x = x, y = y, alpha = 1)
dim(ajuste_lasso$beta)
[1]  5 83
length(ajuste_lasso$a0)
[1] 83
ajuste_lasso$lambda
 [1] 0.7194595174 0.6555447167 0.5973079309 0.5442447405 0.4958955378
 [6] 0.4518415450 0.4117011874 0.3751267886 0.3418015586 0.3114368502
[11] 0.2837696588 0.2585603444 0.2355905560 0.2146613403 0.1955914185
[16] 0.1782156160 0.1623834318 0.1479577353 0.1348135779 0.1228371112
[21] 0.1119246008 0.1019815277 0.0929217699 0.0846668559 0.0771452858
[26] 0.0702919112 0.0640473716 0.0583575796 0.0531732531 0.0484494878
[31] 0.0441453689 0.0402236160 0.0366502608 0.0333943526 0.0304276904
[36] 0.0277245782 0.0252616031 0.0230174319 0.0209726267 0.0191094763
[41] 0.0174118431 0.0158650230 0.0144556181 0.0131714208 0.0120013081
[46] 0.0109351449 0.0099636968 0.0090785494 0.0082720361 0.0075371712
[51] 0.0068675897 0.0062574920 0.0057015936 0.0051950798 0.0047335632
[56] 0.0043130465 0.0039298873 0.0035807670 0.0032626615 0.0029728156
[61] 0.0027087189 0.0024680838 0.0022488260 0.0020490465 0.0018670149
[66] 0.0017011544 0.0015500285 0.0014123283 0.0012868609 0.0011725397
[71] 0.0010683745 0.0009734631 0.0008869833 0.0008081862 0.0007363891
[76] 0.0006709704 0.0006113632 0.0005570514 0.0005075645 0.0004624738
[81] 0.0004213889 0.0003839539 0.0003498445

Regresión lasso en R

Fíjate que, de nuevo, a mayor penalización (modelos más simples, más sparsity), el \(R^2\) disminuye (o lo que es lo mismo, aumenta el SSE)

ggplot(tibble("lambda" = log(ajuste_lasso$lambda),
              "R2" = ajuste_lasso$dev.ratio)) +
  geom_line(aes(x = lambda, y = R2)) +
  theme_minimal() + labs(x = "log-lambda", y = "R2")

Regresión lasso en R

De nuevo a mayor penalización (modelos más simples, más sparsity), la norma absoluta \(\left\| \cdot \right\|_1\) de los coeficientes disminuye (ya que penaliza más)

ggplot(tibble("lambda" = log(ajuste_lasso$lambda),
              "norma" = colSums(abs(ajuste_lasso$beta)))) +
  geom_line(aes(x = lambda, y = norma)) +
  theme_minimal() + labs(x = "log-lambda", y = "norma")

Selección de lambda en R

De nuevo podemos hacer uso de cv.glmnet() para elegir la (penalización óptima elegida a través de validación k-fold)

set.seed(12345)
kfolds_lambda_lasso <- cv.glmnet(x = x, y = y, alpha = 1, nfolds = 20)
min(kfolds_lambda_lasso$lambda)
[1] 0.0003498445
kfolds_lambda_lasso$lambda.min
[1] 0.0003498445
kfolds_lambda_lasso$lambda.1se
[1] 0.009078549
lambda_seq <- 10^seq(log10(max(kfolds_lambda_ridge$lambda)), log10(0.0001), l = 200)

set.seed(12345)
kfolds_lambda_lasso <- cv.glmnet(x = x, y = y, alpha = 1, nfolds = 20, lambda = lambda_seq)
min(kfolds_lambda_lasso$lambda)
[1] 1e-04
kfolds_lambda_lasso$lambda.min
[1] 0.0002042283
kfolds_lambda_lasso$lambda.1se
[1] 0.009205919

Selección de lambda en R

Con la función plot() podemos dibujar de nuevo de manera automática los distintos valores \(CV_k(\lambda)\): los números que aparecen arriba indican el número de predictores distintos de 0

plot(kfolds_lambda_lasso)
abline(h = min(kfolds_lambda_lasso$cvm) + c(0, min(kfolds_lambda_lasso$cvsd)))

Regresión lasso en R

La predicción se realiza de la misma manera (depende de si queremos el promedio en validación o del conjunto de train original): ahora algunos coeficientes están a 0 lo cual nos puede servir no solo como modelo sino como método de selección de variables (ejecutando ahora una regresión sin penalizar pero sin las variables cuyos coeficientes son 0).

predict(kfolds_lambda_lasso, type = "coefficients",
        s = kfolds_lambda_lasso$lambda.1se, newx = x)
6 x 1 sparse Matrix of class "dgCMatrix"
                            s1
(Intercept)         2.21977156
Sepal.Width         0.58624860
Petal.Length        0.54669558
Petal.Width        -0.17323359
Species_versicolor  .         
Species_virginica  -0.04653838

Mínimo global

Uno de los problemas de la validación cruzada es que podemos encontrarnos con lo que se conoce como L-shaped problem: la curva de la cual pretendemos encontrar el mínimo no tiene un mínimo global. Veamos un ejemplo.

Código
set.seed(12345)
p <- 100
n <- 300
x <- matrix(rnorm(n * p), n, p)
y <- 1 + rnorm(n)
lambdaGrid <- exp(seq(-10, 3, l = 200))

# Validación leave-1-out
plot(cv.glmnet(x = x, y = y, alpha = 1, nfolds = n, lambda = lambdaGrid))

Esto puede ser un indicio de que el intercepto es significativo y no así el resto de predictores

Consistencia reg. penalizada

Como sucedía entre el BIC/AIC, en Zhao and Yu (2006) se prueba como la selección lasso es consistente bajo la one standard rule (bajo ciertas condiciones), pero no así cuando \(\lambda\) es seleccionado por validación cruzada (incluye demasiados más predictores de los necesarios, como le sucedía al AIC)

Imagen obtenida de https://bookdown.org/egarpor/PM-UC3M/lm-iii-shrink.html#lm-iii-shrink-varsel

Deberes

Diseña un estudio que, con el dataset iris y considerando el modelo general de elastic-net, nos determina el mejor par \(\left(\alpha, \lambda \right)\) (de manera conjunta).

Clase n+1: regresión logística

Introducción a los modelos lineales generalizados (glm)

Regresión logística

Hasta ahora nuestra variable objetivo era una variable continua. De hecho si echamos la vista atrás, siempre y cuando se cumplan las hipótesis, el modelo lineal cumple que

\[Y | \left( \boldsymbol{X} = \left( x_1, \ldots, x_p \right) \right) \sim N \left(\beta_0 + \beta_1x_1 + \ldots + \beta_p x_p, \sigma_{\varepsilon}^2 \right)\]

¿Pero qué sucede cuando nuestra variable objetivo es cualitativa (o cuantitativa discreta)?

Regresión logística

Vamos a empezar por el caso más sencillo: una variable objetivo binaria. Para ello vamos a cargar el dataset manif.rda (basta con usar load() ya que la extensión .rda es un fichero nativo de R)

load(file = "./datos/manif.rda")
data <- as_tibble(manif)

 

En este caso la variable objetivo será man (binaria, asistió una persona o no a un manifestación) y contamos solo con una predictora, la edad, que es una cuantitativa continua.

data |> count(man)
# A tibble: 2 × 2
    man     n
  <int> <int>
1     0    76
2     1    24

Regresión logística

Código
ggplot(manif |> mutate(man = factor(man, labels = c("Sí", "No"))),
       aes(x = edad, y = man)) +
  geom_point(aes(color = man), size = 3, alpha = 0.7) +
  theme_minimal() +
  labs(x = "Edad (años)", y = "¿Asistió a manifestación?",
       color = "Asistencia",
       title = "¿Cómo usar la regresión lineal para CLASIFICAR?",
       subtitle = "Reminder: la regresión lineal debe cumplir unas hipótesis")

Dado que man es una variable binaria podemos pensarla de la misma manera que pensamos en el lanzamiento de una moneda: tiene dos estados posibles (1 y 0), y cada uno con una probabilidad asignada

Regresión logística

Si lo equiparamos al lanzamiento de una moneda, podemos asumir que la variable objetivo man sigue una distribución binomial (con n = 1) o de Bernoulli con probabilidad de éxito (asistir) \(p\) (y \(1-p\) de no asistir) tal que

\[Y | \left( \boldsymbol{X} = \left( x_1, \ldots, x_p \right) \right) \sim Ber \left(p \right)\]

Como sucedía en el caso de la regresión lineal, nuestro objetivo por tanto será estimar la esperanza de la variable objetivo condicionada a la información de las predictoras

\[E \left[ Y | \left( \boldsymbol{X} = \left( x_1, \ldots, x_p \right) \right) \right] = p, \quad \widehat{E} \left[ Y | \left( \boldsymbol{X} = \left( x_1, \ldots, x_p \right) \right) \right] = \widehat{p}\] Lo importante: ¿cómo aprovechar lo que sabemos de modelización lineal para estimar \(\widehat{p}\)?

Regresión logística

La idea más inmediata es, precisamente, asignarle a esa esperanza estimada la misma estimación de la regresión, es decir,

\[\widehat{E} \left[ Y | \left( \boldsymbol{X} = \left( x_1, \ldots, x_p \right) \right) \right] = \widehat{p} = \widehat{\beta}_0 + \widehat{\beta}_1 X_1 + \ldots + \widehat{\beta}_p X_p\]

de manera que, si \(G\) son las modalidades que puede adoptar la variable objetivo (en nuestro caso \(G = \left\lbrace 0,1 \right\rbrace\))

\[\widehat{y_i} = j \quad \text{si} \quad \widehat{P}(Y = j | X = \left(x_{1}, \ldots, x_{p} \right)) = \arg \max_{g \in G} \widehat{P}(Y = g | X = \left(x_{1}, \ldots, x_{p} \right))\]

Regresión logística

¿Qué pasaría si aplicásemos un modelo lineal?

set.seed(12345)
library(rsample)
split <- initial_split(data, prop = 0.8, strata = man)
train <- training(split)
test <- testing(split)

ajuste_lineal <- lm(data, formula = man ~ .)
predict(ajuste_lineal, test)
          1           2           3           4           5           6 
 0.55209036  0.22791718  0.15424146  0.05109545  0.55209036  0.33106319 
          7           8           9          10          11          12 
 0.12477117  0.24265233  0.16897660  0.09530088  0.03636030  0.28685776 
         13          14          15          16          17          18 
 0.28685776  0.06583059  0.15424146 -0.03731542  0.44894435  0.40473892 
         19          20          21 
 0.39000377  0.31632805  0.08056574 

Tenemos un problema: si aplicamos una regresión lineal, no hay nada que nos garantice que la probabilidad predicha esté entre 0 y 1. ¡Puede salir incluso negativa!

Regresión logística

La idea será garantizar que la salida de la regresión lineal acabe dentro del rango [0,1], encapsulando esa salida con una función de enlace \(g^{-1}:\mathbb{R} \to [0,1]\) tal que

\[P(Y = 1 | X = \left(x_{1}, \ldots, x_{p} \right)) = p = g^{-1} \left(\beta_0 + \beta_1 x_1 + \ldots + \beta_p x_p\right) \]

\[g(p) = \beta_0 + \beta_1 x_1 + \ldots + \beta_p x_p\]

 

¿Qué condiciones debe cumplir dicha función? ¿Nos vale cualquiera?

Regresión logística

Necesitamos unas mínimas condiciones de regularidad

  • Invertible: la función \(g:~[0,1] \to \mathbb{R}\) debe ser invertible, tal que \(g^{-1}:~\mathbb{R} \to [0,1]\)

  • Soporte en [0, 1]: la función \(g:~[0,1] \to \mathbb{R}\) debe estar definida para cualquier valor en \([0, 1]\)

  • Codiminio real: la función \(g^{-1}:~\mathbb{R} \to [0,1]\) debe estar definida para todo valor real (para cualquier salida de una regresión)

  • Monótona creciente: dado que \(\beta_0 + \beta_1 x_1 + \ldots + \beta_p x_p\) cuantifica el efecto de los predictores en la probabilidad de éxito de la variable objetivo, la probabilidad debe aumentar según crezca dicha cantidad, así que \(g^{-1}\) nunca debe decrecer

Regresión logística

En función de las distintas funciones de enlace podemos obtener distintos modelos (siempre y cuando cumplan las 4 propiedades anteriores):

  • Enlace uniforme (unit): función \(g^{−1}(x) = x I_{0< x < 1} + I_{x \geq 1}\) (nos devuelve la salida de la regresión si está entre 0 y 1, y lo trunca a 0 en cualquier otro caso).
  • Enlace probit: función \(g^{−1}(x) = \Phi (x)\) donde \(\Phi(x)=\frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} e^{-\frac{u^2}{2}} du\) es la función de distribución acumulada de una normal
  • Enlace logit: función \(g^{−1}(x) = logistic(x):= \frac{e^{x}}{1 + e^{x}} = \frac{1}{1 + e^{-x}}\) basada en la distribución acumulada de una distribución logística.

 

Esta última es la abordaremos en este curso por ser la más común y fácil de interpretar

Regresión logística

Regresión logística

La función de enlace más usada es la función logit, que no es más que la inversa de la función \(logistic\) que hemos definido arriba.

 

Para ahorar notación llamaremos \(\eta = \beta_0 + \beta_1 x_1 + \ldots + \beta_p x_p\) a la salida de la regresión y \(\widehat{\eta} = \widehat{\beta}_0 + \widehat{\beta}_1 x_1 + \ldots + \widehat{\beta}_p x_p\) a su estimación

\[g(p) = logit(p) = \ln \left( \frac{p}{1-p} \right) = \eta\]

\[p = g^{−1}(\eta) = logistic(\eta):= \frac{1}{1 + e^{-\eta}}\]

Regresión logística

\[g(p) = logit(p) = \ln \left( \frac{p}{1-p} \right) = \eta, \quad p = g^{−1}(\eta) = logistic(\eta):= \frac{1}{1 + e^{-\eta}}\]

Así la estimación buscada será

\[\widehat{p} = g^{−1}(\widehat{\eta} ) = \frac{1}{1 + e^{-\widehat{\eta} }} = \frac{1}{1 + e^{-\left( \widehat{\beta}_0 + \widehat{\beta}_1 x_1 + \ldots + \widehat{\beta}_p x_p\right) }}\]

¿Cómo interpretar la salida de la regresión? ¿Qué sucederá con \(\widehat{p}\) cuando \(\widehat{\eta} > 0\)? ¿Y \(\widehat{\eta}< 0\)?

Regresión logística

\[\widehat{p} = g^{−1}(\widehat{\eta} ) = \frac{1}{1 + e^{-\widehat{\eta} }} = \frac{1}{1 + e^{-\left( \widehat{\beta}_0 + \widehat{\beta}_1 x_1 + \ldots + \widehat{\beta}_p x_p\right) }}\]

  • Si \(\widehat{\eta} = \widehat{\beta}_0 + \widehat{\beta}_1x_1 +\ldots + \widehat{\beta}_p x_p = 0\), \(\widehat{p} = g^{−1}(\widehat{\eta} ) = \frac{1}{1 + e^{-\widehat{\eta} }} = 0.5\) –> la probabilidad de éxito es la misma que la de fracaso

  • Si \(\widehat{\eta} = \widehat{\beta}_0 + \widehat{\beta}_1x_1 +\ldots + \widehat{\beta}_p x_p > 0\), \(\widehat{p} = g^{−1}(\widehat{\eta} ) = \frac{1}{1 + e^{-\widehat{\eta} }} > \frac{1}{1 + e^{0}} > 0.5\) –> la probabilidad de éxito es mayor que la de fracaso –> \(\widehat{y} = 1\)

  • Si \(\widehat{\eta} = \widehat{\beta}_0 + \widehat{\beta}_1x_1 +\ldots + \widehat{\beta}_p x_p < 0\), \(\widehat{p} = g^{−1}(\widehat{\eta} ) = \frac{1}{1 + e^{-\widehat{\eta} }} > \frac{1}{1 + e^{0}} < 0.5\) –> la probabilidad de éxito es menor que la de fracaso –> \(\widehat{y} = 0\)

Cuotas u odds

\[g(p) = logit(p) = \ln \left( \frac{p}{1-p} \right) = \eta, \quad p = g^{−1}(\eta) = logistic(\eta):= \frac{1}{1 + e^{-\eta}}\]

Como hemos visto antes, la inversa es lo que se conoce como función logística tal que \(g(p) = logit(p) = \ln \left( \frac{p}{1-p} \right) = \eta\) es precisamente la salida de la regresión (lo que nos facilita su interpretación)

De hecho al cociente entre \(p\) y \(1-p\) es lo que se conoce como odds o cuotas: probabilidad de éxito vs probabilidad de fracaso (como las cuotas de las casas de apuestas), cuantificando cuántas veces es más o menos probable el éxito que el fracaso

Cuotas u odds

Así podemos definir los log-odds como el logaritmo de las cuotas

\[log-odds(Y)= \ln \left(\frac{p}{1-p} \right) = \eta = \beta_0 + \beta_1 x_1 + \ldots + \beta_p x_p\] \[odds (Y) = \frac{p}{1-p} = e^{\eta} = e^{\beta_0 + \beta_1 x_1 + \ldots + \beta_p x_p} = e^{\beta_0} e^{\beta_1x_1} \ldots e^{\beta_px_p} \]

  • \(e^{\beta_0}\): estimación (media) de odds cuando todas las predictoras son cero –> el ratio esperado p:1−p (éxito vs fracaso) cuando las predictoras se anulan

  • \(e^{\beta_j}\) para \(j\geq 1\): incremento MULTIPLICATVO medio de odds para un incremento ADITIVO unitario de la predictora \(X_j\) (siempre y cuando el resto permanezcan fijas) –> por cada unidad que se incremente \(X_j\), el ratio éxito/fracaso se verá MULTIPLICADO por \(e^{\widehat{\beta}_j}\) (si \(\beta_j\) es positivo, aumenta; en caso contrario, disminuye).

Ajuste en R

Vamos a realizar el ajuste en R a nuestro dataset haciendo uso de la función glm(), indicándole además de los datos y la variable objetivo, la distribución de la variable objetivo (binomial) y la función de enlace (logit)

ajuste_logit <- glm(data = train, formula = man ~ ., family = binomial(link = "logit"))
ajuste_logit

Call:  glm(formula = man ~ ., family = binomial(link = "logit"), data = train)

Coefficients:
(Intercept)         edad  
    2.59211     -0.08926  

Degrees of Freedom: 78 Total (i.e. Null);  77 Residual
Null Deviance:      87.16 
Residual Deviance: 74.57    AIC: 78.57

El ajuste obtenido es \(\widehat{\eta} = 2.592 - 0.089 X\), es decir,

\[\widehat{p} = \widehat{P} \left(Y = 1 | X = x\right) = \frac{1}{1 + e^{-\left(2.592 - 0.089 X \right)}}, \quad \widehat{odds} (Y) = \frac{\widehat{p} }{1-\widehat{p} } = e^{\widehat{\eta}} = e^{2.592} e^{- 0.089x}\]

Ajuste en R

\[\widehat{p} = \widehat{P} \left(Y = 1 | X = x\right) = \frac{1}{1 + e^{-\left(2.592 - 0.089 X \right)}}, \quad \widehat{odds} (Y) = \frac{\widehat{p} }{1-\widehat{p} } = e^{\widehat{\eta}} = e^{2.592} e^{- 0.089x}\]

Si calculamos la exponencial de los coeficientes

exp(coef(ajuste_logit))
(Intercept)        edad 
 13.3579179   0.9146036 
  • Es 13 veces más probable que una persona de 0 años asista a la manifestación a que no lo haga (estimación poco útil ya que \(X=0\) no está dentro del rango)

  • Por cada año que cumpla la persona, la probabilidad de asistir frente a no asistir se MULTIPLICA por 0.9146, es decir, se reduce un 8.54%.

Evaluación: deviance

ajuste_logit

Call:  glm(formula = man ~ ., family = binomial(link = "logit"), data = train)

Coefficients:
(Intercept)         edad  
    2.59211     -0.08926  

Degrees of Freedom: 78 Total (i.e. Null);  77 Residual
Null Deviance:      87.16 
Residual Deviance: 74.57    AIC: 78.57

En los modelos de clasificación no podemos calcular métricas como el error cuadrático medio, var residual o \(R^2\) ya que no tenemos una objetivo continua.

Como vemos en la salida se nos proporciona dos medidas: Null Deviance y Residual Deviance

Evaluación: deviance

ajuste_logit

Call:  glm(formula = man ~ ., family = binomial(link = "logit"), data = train)

Coefficients:
(Intercept)         edad  
    2.59211     -0.08926  

Degrees of Freedom: 78 Total (i.e. Null);  77 Residual
Null Deviance:      87.16 
Residual Deviance: 74.57    AIC: 78.57
  • Null Deviance: al igual que en la bondad de ajuste de un modelo lineal se usa como referencia la SST, la anomalía nula o null deviance es la «diferencia» al comparar la log-verosimilitud del modelo perfectamente sobreajustado vs un modelo sin parámetros (solo intercepto)

\[D_0 = -2(\mathcal{L}_{\hat{\beta}_0} - \mathcal{L}_{saturado}) \geq 0, \quad \mathcal{L}_{modelo} = ln(P(observado | modelo))\]

Evaluación: deviance

ajuste_logit

Call:  glm(formula = man ~ ., family = binomial(link = "logit"), data = train)

Coefficients:
(Intercept)         edad  
    2.59211     -0.08926  

Degrees of Freedom: 78 Total (i.e. Null);  77 Residual
Null Deviance:      87.16 
Residual Deviance: 74.57    AIC: 78.57
  • Residual Deviance: la anomalía residual o anomalía del modelo o residual deviance es la «diferencia» al comparar la log-verosimilitud del modelo perfectamente sobreajustado vs nuestro modelo

\[D = -2(\mathcal{L}_{\widehat{\beta}} - \mathcal{L}_{saturado}) \geq 0\]

Evaluación: deviance

Dado que \(D_0\) mide la máxima distancia posible al modelo «perfecto» (sobreajustado) y \(D\) mide la distancia de nuestro modelo, una forma de evaluar una regresión logística es con el conocido como pseudo-\(R^2\) o coeficiente de McFadden (de todo lo mal que lo podía hacer, ¿cuánto lo he hecho?)

\[pseudo-R^2 = \frac{D_0 - D}{D_0} = 1 - \frac{D}{D_0}\]

El coeficiente de McFadden podemos extraerlo con la función pR2() del paquete {pscl}.

library(pscl)
pR2(ajuste_logit)
fitting null model for pseudo-r2
        llh     llhNull          G2    McFadden        r2ML        r2CU 
-37.2838586 -43.5813660  12.5950149   0.1445000   0.1473708   0.2205377 
1 - ajuste_logit$deviance/ajuste_logit$null.deviance
[1] 0.1445

Predicción logística

Una vez realizado el ajuste podemos usar como es habitual predict(..., type = "response") para obtener las predicciones de un conjunto de datos (en este caso lo vamos a hacer en train), teniendo claro que nos devuelve la probabilidad de ser 1

train |>
  mutate(prob_1 = predict(ajuste_logit, train, type = "response"),
         prob_0 = 1 - prob_1)
# A tibble: 79 × 4
     man  edad prob_1 prob_0
   <int> <int>  <dbl>  <dbl>
 1     0    40 0.273   0.727
 2     0    41 0.256   0.744
 3     0    40 0.273   0.727
 4     0    51 0.123   0.877
 5     0    56 0.0827  0.917
 6     0    54 0.0972  0.903
 7     0    50 0.133   0.867
 8     0    50 0.133   0.867
 9     0    50 0.133   0.867
10     0    48 0.155   0.845
# ℹ 69 more rows

Predicción logística

Con dichas probabilidades podemos construir un dataset donde proporcionemos las cuotas, los log-odds y la clasificación de la clase propiamente dicha (fijando un umbral a partir del cual una observación es clasificada como 1)

umbral <- 0.5
pred_train <-
  train |> 
  mutate(prob_1 = predict(ajuste_logit, train, type = "response"), prob_0 = 1 - prob_1,
         odds = prob_1 / prob_0, log.odds = log(odds),
         pred_class = if_else(prob_1 > umbral, 1, 0))
pred_train
# A tibble: 79 × 7
     man  edad prob_1 prob_0   odds log.odds pred_class
   <int> <int>  <dbl>  <dbl>  <dbl>    <dbl>      <dbl>
 1     0    40 0.273   0.727 0.376    -0.978          0
 2     0    41 0.256   0.744 0.344    -1.07           0
 3     0    40 0.273   0.727 0.376    -0.978          0
 4     0    51 0.123   0.877 0.141    -1.96           0
 5     0    56 0.0827  0.917 0.0901   -2.41           0
 6     0    54 0.0972  0.903 0.108    -2.23           0
 7     0    50 0.133   0.867 0.154    -1.87           0
 8     0    50 0.133   0.867 0.154    -1.87           0
 9     0    50 0.133   0.867 0.154    -1.87           0
10     0    48 0.155   0.845 0.184    -1.69           0
# ℹ 69 more rows

Matriz de confusión

Una de las formas más habituales de evaluar un modelo de clasificación es con la matriz de confusión: ¿cuántos 1’s reales son clasificados como 1? ¿cuántos 0’s? ¿cuántos 1’s reales son erroneamente clasificados como 0?

Para ello usamos conf_mat() del paquete {yardstick} (paquete de {tidymodels}), indicándole la columna con la clase real (truth = ...) y la columna con la clase predicha (estimate = ...), ambas como factor.

library(yardstick)
conf_mat_train <-
  pred_train |>
  mutate(man = as_factor(man), pred_class = as_factor(pred_class)) |> 
  conf_mat(truth = man, estimate = pred_class)
conf_mat_train
          Truth
Prediction  0  1
         0 59 10
         1  1  9

El modelo ha acertado 59 de las 60 observaciones etiquetadas con 0 y ha acertado solo 9 de las 19 observaciones etiquetadas con 1

Métricas de evaluación

conf_mat_train
          Truth
Prediction  0  1
         0 59 10
         1  1  9

Esa matriz de confusión puede ser analizada de diversas formas en función de la pregunta que queramos contestar.

 

De ahora en adelante llamaremos verdadero positivo/negativo (VP/VN) a las observaciones 1’s/0’s que han sido clasificadas correctamente como 1’s/0’s, respectivamente. En nuestro ejemplo tenemos 59 VN y 9 VP.

Llamaremos falsos positivo/negativo (FP/FN) a las observaciones 0’s/1’s que han sido erroneamente clasificadas como 1’s/0’s, respectivamente. En nuestro caso tenemos 10 FN y 1 FP.

Métricas de evaluación

conf_mat_train
          Truth
Prediction  0  1
         0 59 10
         1  1  9
  • ¿Qué proporción de observaciones han sido bien clasificadas

 

Dicha métrica se conoce como accuracy (la suma de la diagonal de la matriz de confusión entre el total de observaciones, en nuestro caso un \(0.86076\))

\[accuracy = \frac{VP + VN}{total}\]

Métricas de evaluación

conf_mat_train
          Truth
Prediction  0  1
         0 59 10
         1  1  9
  • ¿Qué proporción de positivos reales han sido bien clasificados?

 

Dicha métrica se conoce como sensibilidad (sensitivity) y aproxima la probabilidad de que una observación positiva sea clasificada como tal por el modelo (la probabilidad de que una PCR dé positiva en una persona que tiene realmente el COVID), en nuestro caso un \(0.47368\)

\[sensitivity = \frac{VP}{VP + FN}\]

Métricas de evaluación

conf_mat_train
          Truth
Prediction  0  1
         0 59 10
         1  1  9
  • ¿Qué proporción de negativos reales han sido bien clasificados?

 

Dicha métrica se conoce como especificidad (specificity) y aproxima la probabilidad de que una observación negativa sea clasificada como tal por el modelo (la probabilidad de que una PCR dé negativa en una persona que no tiene COVID), en nuestro caso un \(0.98333\)

\[specificity = \frac{VN}{VN + FP}\]

Métricas de evaluación

conf_mat_train
          Truth
Prediction  0  1
         0 59 10
         1  1  9
  • ¿Qué concordancia (no aleatoria) ahí entre la realidad y la predicción?

Dicha métrica se conoce como coeficiente kappa (de Cohen) y aproxima la probabilidad de que la clase real y la clase predicha coincidan más allá de por motivos de azar (asumiendo independencia)

\[\kappa = \frac{P(a) - P(e)}{1-P(e)} = \frac{acc - P(e)}{1-P(e)}\]

donde \(P(a)\) es la probabilidad empírica de acuerdo (es decir, el accuracy) y \(P(e)\) la probabilidad empírica de acuerdo espúreo (acuerdo por azar)

Métricas de evaluación

conf_mat_train
          Truth
Prediction  0  1
         0 59 10
         1  1  9

\[\kappa = \frac{P(a) - P(e)}{1-P(e)} = \frac{acc - P(e)}{1-P(e)}\]

\[\begin{eqnarray}P(e) &=& P(real = 0, pred = 0) + P(real = 1, pred = 1) \nonumber \\ &=& \frac{VN + FN}{total}* \frac{VN + FP}{total} + \frac{FP + VP}{total}* \frac{FN + VP}{total}\end{eqnarray}\]

En este caso \(P(e) = 0.6937991\) y \(\kappa = 0.54605\).

Métricas de evaluación

conf_mat_train
          Truth
Prediction  0  1
         0 59 10
         1  1  9
  • ¿Qué proporción de los clasificados como positivos son realmente verdaderos positivos

 

Dicha métrica se conoce como precision y aproxima la probabilidad de que una observación clasificada como positiva fuese realmente positiva (la probabilidad de que estés con covid si una PCR te da positiva), en nuestro caso un \(0.9\)

\[precision = \frac{VP}{VP + FP}\]

Métricas de evaluación

conf_mat_train
          Truth
Prediction  0  1
         0 59 10
         1  1  9
  • ¿Qué proporción de los clasificados como positivos son mal clasificados como negativos

 

Dicha métrica se conoce como tasa de descubrimiento falso o false discovery rate (FDR) (tasa de error de tipo I) y aproxima la probabilidad de que una observación clasificada como positiva no fuese realmente positiva (la probabilidad de que estés sano a pesar de que una PCR te da positiva), en nuestro caso un \(0.1\)

\[FDR = \frac{FP}{VP + FP} = 1 - precision\]

Métricas de evaluación

Las métricas mencionadas, amén de otras muchas, pueden ser obtenidas haciendo summary() de la matriz de confusión (indicando que el evento de interés es el segundo nivel del factor, los 1’s)

conf_mat_train |> summary(event_level = "second")
# A tibble: 13 × 3
   .metric              .estimator .estimate
   <chr>                <chr>          <dbl>
 1 accuracy             binary         0.861
 2 kap                  binary         0.545
 3 sens                 binary         0.474
 4 spec                 binary         0.983
 5 ppv                  binary         0.9  
 6 npv                  binary         0.855
 7 mcc                  binary         0.587
 8 j_index              binary         0.457
 9 bal_accuracy         binary         0.729
10 detection_prevalence binary         0.127
11 precision            binary         0.9  
12 recall               binary         0.474
13 f_meas               binary         0.621

Métricas de evaluación

conf_mat_train |>
  summary(event_level = "second") |> 
  slice(1:4)
# A tibble: 4 × 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy binary         0.861
2 kap      binary         0.545
3 sens     binary         0.474
4 spec     binary         0.983

Si nos fijamos en las primeras

  • a priori parece que el modelo funciona bastante bien ya que acierta el 86% de las veces
  • de hecho la probabilidad de acierto en el caso de personas que no asisten a la manifestación (0’s) es del 98.3% pero …
  • la probabilidad de acierto en el caso de los asistentes (1’s) es solo del 47.4%

Métricas de evaluación

Código
conf_mat_train |>
  summary(event_level = "second") |> 
  slice(1:4)
# A tibble: 4 × 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy binary         0.861
2 kap      binary         0.545
3 sens     binary         0.474
4 spec     binary         0.983

Si pensamos nuestro modelo como una PCR

  • la prueba acierta el 86% de las veces

  • de hecho la probabilidad de que la PCR dé negativa estando sano es del 98.3% (no hay apenas falsos positivos) pero …

  • la probabilidad de que la PCR dé positiva estando enfermo (1’s) es solo del 47.4% (hay muchísimos falsos negativos)

Si el objetivo fuese, por ejemplo evitar contagios…no parece muy buen método

Métricas de evaluación

Esto puede suceder por muchos motivos pero el principal es que, en nuestro caso, la muestra está desbalanceada: el modelo ha podido aprender muy bien de los 0’s pero mucho menos de los 1’s (ya que solo hay un 25%).

train |> count(man)
# A tibble: 2 × 2
    man     n
  <int> <int>
1     0    60
2     1    19

Por eso es importante no quedarse solo con el accuracy en la evaluación del modelo. De hecho, según \(\kappa\), la probabilidad de concordancia entre realidad y predicción, excluyendo el azar, es solo de \(0.545\).

Curva ROC

Todas las métricas anteriores se basan en una decisión: todo depende del umbral (threshold) que hayamos usado para decidir que una observación es clasificada como 1 o como 0 (ya que, recuerda, el ajuste nos devuelve solo la probabilidad estimada de serlo \(\widehat{p}\))

¿Qué pasaría si ese umbral fuese variando, desde 0 hasta 1? ¿Cómo de bueno sería nuestro método si movemos ese umbral?

  • Si el umbral tiende a cero –> casi todas las observaciones serán predichas como 1’s –> no habrá falsos negativos (todo positivo será clasificado como positivo) pero ningún negativo real será clasificado como tal –> sensibilidad tiende a 1, especificidad tiende a 0

  • Si el umbral tiende a 1 –> casi todas las observaciones serán predichas como 0’s –> no habrá falsos positivos pero ningún positivo real será clasificado como tal –> sensibilidad tiende a 0, especificidad tiende a 1

Curva ROC

Podemos calcular cada valor de sensibilidad y especificidad con roc_curve() indicando la clase real a predecir (como factor), indicando el evento de interés (los 1’s) y, en este caso, necesitamos proporcionar la probabilidad estimada de ser 1 (en lugar de la clase predicha en sí como antes)

Lo que nos devuelve es la sensibilidad y especificidad para cada umbral

roc_data <- 
  pred_train |>
  mutate(man = as_factor(man)) |> 
  roc_curve(man, prob_1, event_level = "second")
roc_data
# A tibble: 33 × 3
   .threshold specificity sensitivity
        <dbl>       <dbl>       <dbl>
 1  -Inf           0            1    
 2     0.0545      0            1    
 3     0.0593      0            0.947
 4     0.0645      0.0167       0.947
 5     0.0761      0.0333       0.895
 6     0.0827      0.0500       0.895
 7     0.0897      0.117        0.789
 8     0.0972      0.133        0.789
 9     0.105       0.183        0.789
10     0.114       0.233        0.789
# ℹ 23 more rows

Curva ROC

Podemos visualizar dicho dataset añadiendo autoplot() al dataset, lo que nos proporciona la conocida como curva ROC: una curva que nos permite visualizar de manera global la calidad de nuestro clasificador para cada posible decisión

roc_data |> autoplot()

Curva ROC

Fíjate que en el eje X se visualiza 1 - especificidad, es decir, la probabilidad empírica de que el método clasifique mal los 0’s: la curva ROC visualiza la probabilidad de falso positivo vs verdadero positivo

Cuando el umbral es muy bajo, dado que todo será clasificado como positivo, clasifica perfecto los verdaderos positivos (sensibilidad = 1) pero todos los negativos son falsos positivos (1-especificidad = 1). Y al contrario.

Curva ROC

Lo importante de esta curva es lo que se conoce como AUC o área debajo de la curva, que podemos obtener con roc_auc(): buscamos un clasificador cuya área sea lo más próxima a 1.

pred_train |>
  mutate(man = as_factor(man)) |> 
  roc_auc(man, prob_1, event_level = "second")
# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.693

Curva ROC

Sobre/bajomuestreo

Como hemos visto, tener la muestra desbalanceada puede ser un problema. Para resolverlo, podemos (antes de realizar el ajuste) lo que se conoce como sobre/bajo muestreo

La única información que vamos a usar para ello (para no hacer «trampas») es el conjunto de entrenamiento

train |>
  count(man) |> 
  mutate(porc = 100*n/sum(n))
# A tibble: 2 × 3
    man     n  porc
  <int> <int> <dbl>
1     0    60  75.9
2     1    19  24.1

En él observamos que contamos con un 75.9% de 0’s y 19% de 1’s. ¿Cómo equilibrar los porcentajes? Parece obvio que hay dos formas: subir el número de 1’s o bajar el número de 0’s.

Sobre/bajomuestreo

Las dos ideas buscan equilibrar dichos porcentajes

  • Bajomuestreo: dejamos fija la cantidad de la clase minoritaria (los 1’s en este caso) y reducimos la cantidad de la clase mayoritaria (los 0’s) para que haya algo parecido a un equilibrio.

Para ello filtramos los individuos de la clase mayoritaria, seleccionamos un % aleatorio de ellas (en este caso un 33% ya que basta con tener 20 de las 60 para ese equilibrio) y unimos la clase minoritaria.

set.seed(12345)
downsampling_data <-
  train |> filter(man == 0) |> 
  slice_sample(prop = 0.33, replace = FALSE) |> 
  bind_rows(train |> filter(man == 1))

downsampling_data |>
  count(man) |> mutate(porc = 100*n/sum(n))
# A tibble: 2 × 3
    man     n  porc
  <int> <int> <dbl>
1     0    19    50
2     1    19    50

Sobre/bajomuestreo

  • Sobremuestreo: dejamos fija la cantidad de la clase mayoritaria (los 0’s en este caso) y aumentamos la cantidad de la clase minoritaria (los 1’s) para que haya algo parecido a un equilibrio.

Para ello filtramos los individuos de la clase minoritaria, seleccionamos un % aleatorio de ellas (en este caso un 300% ya pasar de 20 a 60, x3, para ese equilibrio) y unimos la clase mayoritaria.

set.seed(12345)
oversampling_data <-
  train |> filter(man == 1) |> 
  slice_sample(prop = 3, replace = TRUE) |> 
  bind_rows(train |> filter(man == 0))

oversampling_data|>
  count(man) |>  mutate(porc = 100*n/sum(n))
# A tibble: 2 × 3
    man     n  porc
  <int> <int> <dbl>
1     0    60  51.3
2     1    57  48.7

Sobre/bajomuestreo

set.seed(12345)
downsampling_data <-
  train |> filter(man == 0) |> 
  slice_sample(prop = 0.33, replace = FALSE) |> 
  bind_rows(train |> filter(man == 1))
set.seed(12345)
oversampling_data <-
  train |> filter(man == 1) |> 
  slice_sample(prop = 3, replace = TRUE) |> 
  bind_rows(train |> filter(man == 0))
  • Bajomuestreo: la principal ventaja es que los nuevos datos son una submuestra real de los datos originales; la desventaja es que se reduce considerablemente el tamaño muestral.

  • Sobremuestreo: la principal ventaja es que no reducimos el tamaño muestral; la desventaja es que los datos nuevos son generados artificialmente (en este caso, son observaciones repetidas, de ahí el replace = TRUE)

Existen otras técnicas (SMOTE, ROSE, etc) que lo que hacen es simular datos sintéticos, bien usando la distribución de cada variable o bien tomando los datos originales y simulando nuevos perturbando con ruido.

Ajuste con bajomuestreo

Código
ajuste_logit_down <-
  glm(data = downsampling_data, formula = man ~ ., family = binomial(link = "logit"))
pred_train_down <-
  downsampling_data |> 
  mutate(prob_1 = predict(ajuste_logit_down, downsampling_data, type = "response"),
         prob_0 = 1 - prob_1, odds = prob_1 / prob_0, log.odds = log(odds),
         pred_class = if_else(prob_1 > 0.5, 1, 0))

conf_mat_train_down <-
  pred_train_down |>
  mutate(man = as_factor(man), pred_class = as_factor(pred_class)) |> 
  conf_mat(truth = man, estimate = pred_class)

pred_train_down |>
  mutate(man = as_factor(man)) |> roc_auc(man, prob_1, event_level = "second")
# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.661
Código
conf_mat_train_down |> summary(event_level = "second")
# A tibble: 13 × 3
   .metric              .estimator .estimate
   <chr>                <chr>          <dbl>
 1 accuracy             binary         0.711
 2 kap                  binary         0.421
 3 sens                 binary         0.579
 4 spec                 binary         0.842
 5 ppv                  binary         0.786
 6 npv                  binary         0.667
 7 mcc                  binary         0.436
 8 j_index              binary         0.421
 9 bal_accuracy         binary         0.711
10 detection_prevalence binary         0.368
11 precision            binary         0.786
12 recall               binary         0.579
13 f_meas               binary         0.667

Empeora el accuracy y un poco la especificidad pero mejora la sensibilidad

Ajuste con sobremuestreo

Código
ajuste_logit_over <- 
  glm(data = oversampling_data, formula = man ~ ., family = binomial(link = "logit"))
pred_train_over <-
  oversampling_data |> 
  mutate(prob_1 = predict(ajuste_logit_over, oversampling_data, type = "response"),
         prob_0 = 1 - prob_1, odds = prob_1 / prob_0, log.odds = log(odds),
         pred_class = if_else(prob_1 > 0.5, 1, 0))

conf_mat_train_over <-
  pred_train_over |>
  mutate(man = as_factor(man), pred_class = as_factor(pred_class)) |> 
  conf_mat(truth = man, estimate = pred_class)

pred_train_over |>
  mutate(man = as_factor(man)) |> roc_auc(man, prob_1, event_level = "second")
# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.664
Código
conf_mat_train_over |> summary(event_level = "second")
# A tibble: 13 × 3
   .metric              .estimator .estimate
   <chr>                <chr>          <dbl>
 1 accuracy             binary         0.624
 2 kap                  binary         0.243
 3 sens                 binary         0.509
 4 spec                 binary         0.733
 5 ppv                  binary         0.644
 6 npv                  binary         0.611
 7 mcc                  binary         0.249
 8 j_index              binary         0.242
 9 bal_accuracy         binary         0.621
10 detection_prevalence binary         0.385
11 precision            binary         0.644
12 recall               binary         0.509
13 f_meas               binary         0.569

Empeora el accuracy y mucho la especificidad pero mejora la sensibilidad

Caso real: challenger

Práctica con el dataset challenger.txt subido al campus. El objetivo será obtener la probabilidad de accidente en distintos lanzamientos espaciales, donde fail.field es nuestra variable objetivo (si hubo o no accidente debido a un fallo en las juntas) en función de si hubo o no fallo en los inyectores, la temperatura exterior, la presión de las juntas y la presión de los inyectores

datos <-
  read_delim(file = "./datos/challenger.txt")
datos <- 
  datos |> select(-contains("nfails"))
datos
# A tibble: 23 × 7
   flight date     fail.field fail.nozzle  temp pres.field pres.nozzle
   <chr>  <chr>         <dbl>       <dbl> <dbl>      <dbl>       <dbl>
 1 1      12/04/81          0           0  18.9         50          50
 2 2      12/11/81          1           0  21.1         50          50
 3 3      22/03/82          0           0  20.6         50          50
 4 5      11/11/82          0           0  20           50          50
 5 6      04/04/83          0           1  19.4         50          50
 6 7      18/06/83          0           0  22.2         50          50
 7 8      30/08/83          0           0  22.8        100          50
 8 9      28/11/83          0           0  21.1        100         100
 9 41-B   03/02/84          1           1  13.9        100         100
10 41-C   06/04/84          1           1  17.2        200         200
# ℹ 13 more rows

Pues…ha sido un placer

Espero que hayáis aprendido algo <3