library(tidyverse)
library(tidymodels)
library(ggfortify) # autoplot para diagnósticos
library(rstanarm) # para bayesiana
library(poissonreg) # modelo Poisson
library(car) # para el VIF (multicolinealidad)
library(forecast) # series de tiempo, la forma clásica
library(modeltime) # series de tiempo, la forma tidymodels
library(timetk) # utilerías de series de tiempo
# Siempre que uses tidymodels
tidymodels_prefer()11 Regresiones (parte 1)
11.1 Paquetes y datos a utilizar
A lo largo de esta sección usaremos los siguientes paquetes:
En general usaremos la filosofía tidymodels para combinar los resultados con el tidyverse. Si quieres saber más de tidymodels te recomiendo checar su libro o su página web.
ImportanteEn este capítulo todo lo vas a ver de dos formas
En R casi siempre hay más de una manera de hacer lo mismo, y las regresiones son el caso más claro. Existen dos “dialectos”:
| La forma clásica | La forma tidymodels |
|
|---|---|---|
| Regresión lineal | lm() |
linear_reg() + fit() |
| Regresión logística | glm() |
logistic_reg() + fit() |
| Series de tiempo | forecast (auto.arima) |
modeltime (arima_reg) |
| ¿Cuántos años tiene? | Desde los inicios de R |
Desde ~2020 |
Vamos a mostrarte las dos para cada cosa, con el mismo resultado, porque las dos te vas a topar en la vida real: el código viejo (y la mayoría de los libros y respuestas de internet) está en la forma clásica, mientras que el código nuevo tiende a usar tidymodels.
Ninguna es mejor que la otra. La clásica es más corta y directa para un modelo suelto; tidymodels es más verbosa pero te da la misma sintaxis para cientos de modelos distintos y te resuelve la validación cruzada y el preprocesamiento, que es donde de verdad gana cuando el proyecto crece.
Tú aprende a leer las dos y escribe en la que te guste. Lo único que no debes hacer es mezclarlas a media línea.
11.1.1 Embarazo adolescente y pobreza
Para los datos usaremos la información de Utts y Heckard para determinar si hay una relación entre embarazo adolescente y pobreza.
emb_pob <- read_delim("https://online.stat.psu.edu/stat462/sites/onlinecourses.science.psu.edu.stat462/files/data/poverty/index.txt")Las variables Brth15to17 y Brth18to19 son las tasas brutas de natalidad por cada 1000 mujeres (en el año 2002) en adolescentes de 15 a 17 años y de 18 a 19 respectivamente. La variable PovPct representa la proporción (%) de la población que vive bajo la línea de pobreza en cada una de las entidades de EEUU (Location). Las variables ViolCrime y TeenBrth no se explican por lo que no las usaremos.
12 Regresiones lineales
Usaremos Brth15to17 y PovPct para estudiar si hay una relación entre la tasa de natalidad en adolescentes y el porcentaje de la población en pobreza. Para ello comenzaremos con graficar:
ggplot(emb_pob) +
geom_point(aes(x = PovPct, y = Brth15to17), color = "#bc5090", size = 3) +
labs(
x = "Porcentaje en pobreza",
y = "Tasa bruta de natalidad (por cada 1,000 mujeres adolescentes)",
title = "Relación entre tasa bruta de natalidad en adolescentes\nde 15 a 19 años y porcentaje en pobreza"
) +
theme_bw()
Parece que a mayor porcentaje de pobreza, mayor tasa bruta de natalidad en adolescentes.
12.1 Planteamiento clásico
Si la relación fuera perfecta todos los casos caerían exactamente en una línea recta como sigue:
donde es necesario especificar dos parámetros: el intercepto (el valor que toma cuando la \(x\) en este caso PovPct vale cero) y la pendiente (el valor que relaciona por cada unidad de aumento en PovPct cuánto aumenta Brth15to17).
La ecuación de la línea está dada por:
\[ y = \beta_0 + \beta_1 x \]
donde \(\beta_0\) es el intercepto y \(\beta_1\) la pendiente. Usando la terminología de arriba:
\[ \text{Brth15to17} = \text{Intercepto} + \text{Pendiente}\times \text{PovPct} \] en particular en ese ejemplo:
\[ \text{Brth15to17} = 5 + 3\cdot \text{PovPct} \]
La idea es que el intercepto (\(\beta_0\) ó \(5\)) te indica dónde comienza tu línea cuando no tienes \(x\)’s (es decir cuando \(\text{PovPct} = 0\)). El intercepto controla la altura de la línea como puedes ver en la siguiente gráfica donde puse varios interceptos distintos:
por otro lado la idea de la pendiente (\(\beta_1\) ó \(3\)) es retratar cómo cambia la \(y\) (en este caso \(\text{Brth15to17}\)) por cada unidad que cambia la \(x\) (en este caso \(\text{PovPct}\)). El valor de \(3\) por ejemplo indica que por cada aumento en 1 en \(\text{PovPct}\) la variable \(\text{Brth15to17}\) aumenta en \(3\). Este cambio es proporcional; es decir si ahora \(\text{PovPct}\) aumenta 4 (por decir algo) \(\text{Brth15to17}\) aumenta \(3\times 4 = 12\) unidades. La siguiente gráfica muestra varias líneas todas comenzando en el mismo intercepto de \(5\):
Como ya vimos en la primer figura el mundo no es tan perfecto que todo sea una línea recta. Hay un poco de aleatoriedad involucrada (sea por variables no medidas, por errores de medición o porque el mundo no sea determinista). Por lo cual se plantea que en lugar de que la \(y\) sea exactamente \(\beta_0 + \beta_1 x\) planteamos que la \(y\) proviene de una variable aleatoria normal donde la media de esa normal es \(\beta_0 + \beta_1 x\); es decir:
\[ y \sim \textrm{Normal}( \beta_0 + \beta_1 x, \sigma^2) \]
o dicho de otra manera:
\[ \text{Brth15to17} \sim \textrm{Normal}(\text{Intercepto} + \text{Pendiente}\times \text{PovPct}, \sigma^2) \]
donde la \(\sigma^2\) es la varianza de dicha normal. Puesto gráficamente lo que esto quiere decir es que si, por ejemplo, el intercepto es \(5\) y la pendiente \(3\) entonces cada medición de \(y\) viene de una normal ligeramente distinta:
\[ \text{Brth15to17} \sim \textrm{Normal}(5 + 3\times \text{PovPct}, \sigma^2) \]
Dicho de otra forma, suponemos que en un mundo perfecto los valores de \(y\) (Brth15to17) estarían completamente determinados por los de \(x\) (PovPct) mediante la ecuación de la recta. Pero como el mundo no es perfecto entonces la \(y\) proviene de una normal con promedio dado por la recta. Los puntos (como puedes ver an la siguiente gráfica) se centran más en torno a los promedios de las normales sin embargo están colocados aleatoriamente pues corresponden a distintas realizaciones de \(y\).
Por ejemplo si el porcentaje de pobreza (PovPct) es \(10\) entonces la normal de la que provienen estos datos es:
\[ \text{Brth15to17} \sim \textrm{Normal}(\underbrace{5 + 3\times 10}_{35}, \sigma^2) \]
mientras que si el porcentaje en pobreza es \(20\) entonces la normal es:
\[ \text{Brth15to17} \sim \textrm{Normal}(\underbrace{5 + 3\times 20}_{65}, \sigma^2) \]
Por supuesto que no hay nada de especial con el modelo normal y alguien podría elegir otra distribución (por ejemplo una Gamma) y establecer que:
\[ \text{Brth15to17} \sim \textrm{Gamma}(\text{Intercepto} + \text{Pendiente}\times \text{PovPct}, \beta) \]
Estos modelos son algunos de los lineales generalizados y los discutiremos más adelante. Por ahora nos quedaremos con la idea del modelo dado por:
\[ \text{Brth15to17} \sim \textrm{Normal}(\text{Intercepto} + \text{Pendiente}\times \text{PovPct}, \sigma^2) \]
Nota quizá conoces la regresión lineal bajo la idea clásica de que \[ y = \beta_0 + \beta_1 x + \epsilon \] donde \(\epsilon\sim\text{Normal}(0,\sigma^2)\) son los errores normales. Esta definición es equivalente a la que damos aquí pues por propiedades aditivas de la normal \(\epsilon + \beta_0 + \beta_1 x\) se sigue distribuyendo normal pero con la media ahora dada por lo agregado (\(\beta_0 + \beta_1 x\)). Las ventajas de esta notación es que un modelo para regresión Poisson es simplemente: \[ y \sim \textrm{Poisson}(\beta_0 + \beta_1 x) \] y un modelo para regresión logística es: \[ y \sim \textrm{Bernoulli}\big(\textrm{logit}(\beta_0 + \beta_1 x)\big) \]
12.2 Planteamiento en R
Lo que nos toca ahora es programar nuestro modelo para ello seguiremos la filosofía de tidymodels que pretende unificar bajo la misma notación todos los modelos. La notación básica es como sigue:
modelo %>%
set_engine("tipo de ajuste") %>%
fit("formula a ajustar", data = tus_datos)A partir del ajuste se pueden predecir cosas con predict:
modelo %>%
set_engine("tipo de ajuste") %>%
fit("formula a ajustar", data = tus_datos) %>%
predict()o extraer nuevos datos con extract_fit_engine y tidy:
modelo %>%
set_engine("tipo de ajuste") %>%
fit("formula a ajustar", data = tus_datos) %>%
extract_fit_engine() %>%
tidy()Comencemos con nuestro primer modelo: una regresión lineal clásica dada por:
\[ \text{Brth15to17} \sim \textrm{Normal}(\text{Intercepto} + \text{Pendiente}\times \text{PovPct}, \sigma^2) \]
En R el engine que necesitamos es “lm” que es el clásico:
# Ajusta un modelo lineal
# Brth15to17 = intercepto + pendiente*PovPct
modelo_ajustado <- linear_reg() %>%
set_engine("lm") %>%
fit(Brth15to17 ~ PovPct, data = emb_pob) # Notación y ~ xNota que a diferencia de Stata, R no arroja demasiados resultados. Podemos usar extract_fit_engine combinado con summary para obtenerlos:
# Ajusta un modelo lineal
# Brth15to17 = intercepto + pendiente*PovPct
modelo_ajustado %>%
extract_fit_engine() %>%
summary()
Call:
stats::lm(formula = Brth15to17 ~ PovPct, data = data)
Residuals:
Min 1Q Median 3Q Max
-11.2275 -3.6554 -0.0407 2.4972 10.5152
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.2673 2.5297 1.687 0.098 .
PovPct 1.3733 0.1835 7.483 1.19e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.551 on 49 degrees of freedom
Multiple R-squared: 0.5333, Adjusted R-squared: 0.5238
F-statistic: 56 on 1 and 49 DF, p-value: 1.188e-09o bien con tidy si deseamos nos devuelva una tabla de resultados:
# Ajusta un modelo lineal
# Brth15to17 = intercepto + pendiente*PovPct
modelo_ajustado %>%
extract_fit_engine() %>%
tidy()# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 4.27 2.53 1.69 0.0980
2 PovPct 1.37 0.184 7.48 0.0000000011912.2.1 Lo mismo con lm: la forma clásica
Todo lo anterior se puede hacer sin tidymodels, con la función lm que viene en R desde siempre. La fórmula (y ~ x) es idéntica; lo que cambia es que no hay que declarar el motor ni llamar a fit:
# La misma regresión, en un solo renglón
modelo_clasico <- lm(Brth15to17 ~ PovPct, data = emb_pob)Para ver los resultados, summary directo (sin el extract_fit_engine):
summary(modelo_clasico)
Call:
lm(formula = Brth15to17 ~ PovPct, data = emb_pob)
Residuals:
Min 1Q Median 3Q Max
-11.2275 -3.6554 -0.0407 2.4972 10.5152
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.2673 2.5297 1.687 0.098 .
PovPct 1.3733 0.1835 7.483 1.19e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.551 on 49 degrees of freedom
Multiple R-squared: 0.5333, Adjusted R-squared: 0.5238
F-statistic: 56 on 1 and 49 DF, p-value: 1.188e-09Y si quieres la tabla ordenada, tidy también funciona sobre un lm, porque en realidad tidy nunca fue de tidymodels sino del paquete broom:
tidy(modelo_clasico, conf.int = TRUE, conf.level = 0.90)# A tibble: 2 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 4.27 2.53 1.69 0.0980 0.0260 8.51
2 PovPct 1.37 0.184 7.48 0.00000000119 1.07 1.68Comprobemos que de verdad son el mismo modelo comparando los coeficientes:
# Los coeficientes de las dos rutas
coef(modelo_clasico)(Intercept) PovPct
4.267293 1.373345 coef(extract_fit_engine(modelo_ajustado))(Intercept) PovPct
4.267293 1.373345 # ¿Son iguales?
all.equal(coef(modelo_clasico), coef(extract_fit_engine(modelo_ajustado)))[1] TRUE
TipTabla de traducción entre los dos dialectos
Guarda esta tabla: te va a servir cada vez que encuentres código en el dialecto que no usas.
| Qué quieres | Forma clásica | Forma tidymodels |
|---|---|---|
| Ajustar el modelo | lm(y ~ x, data = d) |
linear_reg() %>% set_engine("lm") %>% fit(y ~ x, data = d) |
| Ver el resumen completo | summary(m) |
m %>% extract_fit_engine() %>% summary() |
| Tabla de coeficientes | tidy(m) |
tidy(m) |
| Coeficientes solos | coef(m) |
coef(extract_fit_engine(m)) |
| Predecir | predict(m, newdata = d) |
predict(m, new_data = d) |
| Intervalos de predicción | predict(m, newdata = d, interval = "prediction") |
predict(m, new_data = d, type = "pred_int") |
| Métricas del ajuste | glance(m) |
glance(m) |
| Diagnósticos | autoplot(m, which = 1:5) |
autoplot(m, which = 1:5) |
Ojo con dos detalles que muerden: en la forma clásica el argumento se llama newdata (junto) y en tidymodels se llama new_data (con guion bajo). Y los intervalos se piden distinto. Son las dos confusiones más comunes al pasar de un dialecto al otro.
Nota¿Y entonces para qué sirve
tidymodels si lm es más corto?
Buena pregunta, porque con un modelo lm gana por mucho. tidymodels empieza a valer la pena cuando:
- Quieres probar varios modelos distintos con la misma sintaxis (una lineal, un bosque aleatorio, una red elástica) y comparar. Con la ruta clásica cada uno tiene su función, sus argumentos y su forma de predecir; con
tidymodelstodos se ajustan igual. - Necesitas validación cruzada o partir en entrenamiento y prueba sin escribirlo a mano.
- Tienes preprocesamiento (estandarizar, crear variables indicadoras, imputar) y quieres que se aplique exactamente igual a los datos nuevos. Ésta es la fuente de errores silenciosos más común en modelos predictivos, y
tidymodelsla resuelve con lasrecipes.
Para un modelo suelto de inferencia —que es la mayoría de lo que se hace en salud pública— lm está perfecto.
Según el tipo de regresión que estemos haciendo es el tipo de tabla que regresa tidy (ver ?tidy). En particular, por ejemplo, podemos modificar para que devuelva intervalos de confianza al 90%:
# Ajusta un modelo lineal
# Brth15to17 = intercepto + pendiente*PovPct
modelo_ajustado %>%
tidy(conf.int = T, conf.level = 0.90)# A tibble: 2 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 4.27 2.53 1.69 0.0980 0.0260 8.51
2 PovPct 1.37 0.184 7.48 0.00000000119 1.07 1.68En este caso, el modelo estima que el intercepto (\(\beta_0\)) es \(4.27\) y la pendiente (\(\beta_1\)) es \(1.37\). Como son estimadores del verdadero valor se denotan con gorrito: \(\hat\beta_0 = 4.27\) y \(\hat\beta_1 = 1.37\).
Podemos utilizar la función de predict para que el modelo nos muestre cómo cree que son los verdaderos valores en relación a los ajustados:
# Ajusta un modelo lineal
# Brth15to17 = intercepto + pendiente*PovPct
predichos <- modelo_ajustado %>%
predict(new_data = emb_pob)
intervalo_predichos <- modelo_ajustado %>%
predict(new_data = emb_pob, type = "pred_int", level = 0.95)
# Juntamos los predichos con los observados
obs_y_modelo <- emb_pob %>%
cbind(predichos) %>%
cbind(intervalo_predichos)
# Graficamos
ggplot(obs_y_modelo) +
geom_ribbon(aes(x = PovPct, ymin = .pred_lower, ymax = .pred_upper),
fill = "#003f5c", linewidth = 1, alpha = 0.5
) +
geom_point(aes(x = PovPct, y = Brth15to17), color = "#bc5090", size = 3) +
geom_line(aes(x = PovPct, y = .pred),
color = "#003f5c", linewidth = 1
) +
labs(
x = "Porcentaje en pobreza",
y = "Tasa bruta de natalidad (por cada 1,000 mujeres adolescentes)",
title = "Relación entre tasa bruta de natalidad en adolescentes\nde 15 a 19 años y porcentaje en pobreza"
) +
theme_bw()
Nada más a ojo no parece que el modelo sea el mejor pues puedes ver que no explica bien la variabilidad (los observados varían mucho respecto al intervalo). La \(R^2\), una métrica que explica cuánto de la varianza captura el modelo tampoco es muy buena:
resumen_ajuste <- modelo_ajustado %>%
extract_fit_engine() %>%
summary()
# R^2 clásica
resumen_ajuste$r.squared[1] 0.533328# R^2 ajustada
resumen_ajuste$adj.r.squared[1] 0.523804Podemos checar las diferentes gráficas de diagnóstico:
library(ggfortify)
autoplot(modelo_ajustado, which = 1:5)
Veamos qué significa cada una de ellas y juguemos un poco con R para irlas modificando.
12.2.2 Residuales contra ajustados
Los residuales son la diferencia entre el modelo (\(\hat{y}\)) y lo real \(y\). En el caso que estábamos trabajando tenemos que con nuestro modelo podemos predecir los valores de Brth15to17 a partir del porcentaje en pobreza PovPct. A los valores predichos por el modelo de Brth15to17 les ponemos un gorro encima y los llamamos: \(\widehat{\text{Brth15to17}}\). Estos están dados por la siguiente función:
\[ \widehat{\text{Brth15to17}} = 4.26 + 1.37 \cdot \text{PovPct} \]
por ejemplo para el porcentaje en pobreza de \(20.1\) obtendríamos:
\[ \widehat{\text{Brth15to17}} = 4.26 + 1.37 \cdot \text{PovPct} = 31.797 \]
por otro lado el verdadero valor de cuando el PovPct es \(20.1\) (estado de Alabama) es \(\text{Brth15to17} = 31.5\). La diferencia entre el verdadero valor (\(\text{Brth15to17} = 31.5\)) y el predicho por el modelo (\(\widehat{\text{Brth15to17}} = 31.797\)) se conoce como el residual. La idea es que en un modelo bueno no debe haber patrones en los residuales (todos deben de flotar en torno al cero pero no mostrar un patrón).
Veamoslo en nuestra base:
obs_y_modelo <- obs_y_modelo %>%
mutate(residuales = Brth15to17 - .pred)
ggplot(obs_y_modelo) +
geom_point(aes(x = .pred, y = residuales), color = "#ff6361") +
labs(
x = "Valores ajustados (predichos)",
y = "Residuales",
title = "Residuales vs ajustados"
) +
theme_bw() +
geom_hline(aes(yintercept = 0), linetype = "dashed")
En esta gráfica el modelo predice mejor rumbo al final que en medio y esto parece estar corroborado por la gráfica del modelo (previa). Nada más para darnos una idea veamos una gráfica de malos residuales y una de buenos
12.2.3 Escala locación
Representa la escala locación contra los residuales estandarizados. La idea de la gráfica es ver que la varianza \(\sigma^2\) del modelo no cambie conforme cambia la \(x\) (propiedad de homoscedasticidad). Para ello graficamos los residuales estandarizados dados por los residuales mismos dividos entre su desviación estándar:
\[ r_{\text{Std}} = \frac{\hat{y} - y}{\text{sd}(\hat{y} - y)} = \frac{\text{Residuales}}{\text{sd}\big(\text{Residuales}\big)} \]
estos residuales estandarizados los podemos calcular en R como sigue:
obs_y_modelo <- obs_y_modelo %>%
mutate(residuales_std = residuales / sd(residuales))Si los graficamos contra los valores ajustados no deberíamos de ver ningún patrón:
ggplot(obs_y_modelo) +
geom_point(aes(x = .pred, y = residuales_std), color = "#ff6361") +
labs(
x = "Valores ajustados (predichos)",
y = "Residuales estandarizados",
title = "Escala Locación"
) +
theme_bw() +
geom_hline(aes(yintercept = 0), linetype = "dashed")
Podemos ver cómo se ven estos puntos en el modelo ideal vs en un modelo donde la \(\sigma^2\) depende de la \(x\):
12.2.4 Normal cuantil cuantil
La segunda gráfica corresponde a una gráfica cuantil cuantil. Esta la utilizamos para verificar la hipótesis de normalidad. En una gráfica cuantil cuantil se grafican los cuantiles de los residuales contra los cuantiles teóricos de la normal. Por ejemplo si la hacemos con sólo 4 puntos se vería algo así:
Podemos armar una gráfica cuantil cuantil con ggplot2:
ggplot(obs_y_modelo, aes(sample = residuales)) +
stat_qq(color = "#ff6361") +
stat_qq_line(color = "#58508d") +
theme_bw() +
labs(
x = "Cuantiles teóricos de la normal",
y = "Cuantiles observados de los residuales",
title = "Gráfica qq"
)
La idea de la gráfica cuantil cuantil es que los puntos sigan la línea lo más posible. Veamos cómo se ve con los datos bien (y los mal)
12.2.5 Residuales contra apalancamiento
El apalancamiento representa qué tanto cambia el modelo al quitar una sola observación. Para poner un ejemplo considera los siguientes datos donde hay un valor atípico y selecciono dos puntos de interés en dos colores:
Veamos cómo cambia la regresión si dejo todos los puntos, si quito el normal y si quito el influyente:
Nota que el modelo no cambia prácticamente nada cuando hago la regresión sin el dato que marqué como normal pero cambia mucho cuando quito el que marqué como influyente. La gráfica de residuales contra apalancamiento muestra también el valor extraño:
El apalancamiento mide la influencia de un dato y en R se puede calcular con hatvalues. Los datos con mayor apalancamiento siempre valen la pena checarlos para verificar que todo opera en orden.
apalancamiento <- modelo_ajustado %>%
extract_fit_engine() %>%
hatvalues()
obs_y_modelo <- obs_y_modelo %>%
cbind(apalancamiento)
# Graficamos residuales contra apalancamiento
ggplot(obs_y_modelo) +
geom_point(aes(x = apalancamiento, y = residuales), color = "#ff6361") +
labs(
x = "Apalancamiento",
y = "Residuales",
title = "Residuales vs ajustados"
) +
theme_bw() +
geom_hline(aes(yintercept = 0), linetype = "dashed")
12.2.6 Distancia de Cook
La distancia de Cook es un concepto similar al apalancamiento que identifica observaciones influyentes. Aquellos valores con distancia de Cook alta vale la pena revisar. En R podemos usar cooks.distance para calcular la distancia de Cook.
distanciaCook <- modelo_ajustado %>%
extract_fit_engine() %>%
cooks.distance()
obs_y_modelo <- obs_y_modelo %>%
cbind(distanciaCook)
# Graficamos residuales contra apalancamiento
ggplot(obs_y_modelo) +
geom_col(aes(x = 1:nrow(obs_y_modelo), y = distanciaCook), fill = "#ff6361") +
labs(
x = "Observación (número de entrada en la base)",
y = "Distancia de Cook",
title = "Distancia de Cook"
) +
theme_bw()
podemos ver que en el ejemplo anterior (el de la observación influyente) la distancia de Cook es exagerada, tan exagerada que ni se alcanzan a ver los otros:
12.3 Regresión múltiple: más de una variable
Hasta ahora la recta depende de una sola variable. En la vida real casi nunca es así: la tasa de natalidad adolescente no depende sólo de la pobreza. Agregar variables se hace sumándolas en la fórmula con un +.
La base también trae ViolCrime (una tasa de crimen violento). Agreguémosla:
modelo_multiple <- lm(Brth15to17 ~ PovPct + ViolCrime, data = emb_pob)
tidy(modelo_multiple, conf.int = TRUE)# A tibble: 3 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 5.98 2.27 2.64 0.0112 1.43 10.5
2 PovPct 1.04 0.183 5.67 0.000000790 0.669 1.40
3 ViolCrime 0.344 0.0877 3.93 0.000275 0.168 0.520modelo_multiple_tm <- linear_reg() %>%
set_engine("lm") %>%
fit(Brth15to17 ~ PovPct + ViolCrime, data = emb_pob)
tidy(modelo_multiple_tm, conf.int = TRUE)# A tibble: 3 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 5.98 2.27 2.64 0.0112 1.43 10.5
2 PovPct 1.04 0.183 5.67 0.000000790 0.669 1.40
3 ViolCrime 0.344 0.0877 3.93 0.000275 0.168 0.520Fíjate en algo importante: al agregar ViolCrime, el coeficiente de PovPct cambió. Pasó de \(1.37\) a \(1.04\). Esto no es un error ni un problema: es la naturaleza de la regresión múltiple.
ImportanteCómo se lee un coeficiente en regresión múltiple
En el modelo de una sola variable, el \(1.37\) significaba: “por cada punto más de pobreza, la tasa sube 1.37”.
En el modelo múltiple, el \(1.04\) significa: “por cada punto más de pobreza, la tasa sube 1.04, manteniendo la tasa de crimen violento constante”.
Ese “manteniendo lo demás constante” (en latín, ceteris paribus) es la frase más importante de toda la regresión múltiple, y es la que más se olvida al escribir los resultados. Los dos coeficientes son correctos: contestan preguntas distintas.
Advertencia
La documentación de esta base no explica qué es exactamente ViolCrime ni TeenBrth. La usamos aquí para mostrar la mecánica del código, pero fíjate en lo que eso implica: si no sabes qué mide una variable, no puedes interpretar su coeficiente ni justificar por qué la incluiste. En un análisis de verdad, eso solo descalifica la variable para un modelo de inferencia.
12.4 Comparación de modelos: AIC, BIC y pruebas anidadas
Ya tenemos dos modelos. ¿Cuál es mejor? No podemos usar la \(R^2\) para decidir, por una razón incómoda:
Importante
La \(R^2\) siempre sube (o se queda igual) cuando agregas variables, aunque las variables sean basura pura. Si le agregas a tu modelo una columna de números aleatorios, la \(R^2\) sube. Por eso la \(R^2\) no sirve para comparar modelos con distinto número de variables.
Para eso existen los criterios de información. Los dos más comunes:
- El AIC (Akaike Information Criterion)
- El BIC (Bayesian Information Criterion)
Los dos funcionan igual de raro: entre más chico, mejor. Los dos miden qué tan bien ajusta el modelo, pero le cobran una multa por cada parámetro que uses. Así, una variable sólo “vale la pena” si mejora el ajuste más de lo que cuesta la multa. La diferencia entre ellos es que el BIC cobra una multa más cara, así que tiende a elegir modelos más pequeños.
# Entre más chico, mejor
AIC(modelo_clasico, modelo_multiple) df AIC
modelo_clasico 3 323.5094
modelo_multiple 4 311.3072BIC(modelo_clasico, modelo_multiple) df BIC
modelo_clasico 3 329.3049
modelo_multiple 4 319.0345glance te da el AIC, el BIC y otras métricas de un jalón, y funciona igual en los dos dialectos:
bind_rows(
glance(modelo_ajustado) %>% mutate(modelo = "PovPct"),
glance(modelo_multiple_tm) %>% mutate(modelo = "PovPct + ViolCrime")
) %>%
select(modelo, r.squared, adj.r.squared, AIC, BIC, df)# A tibble: 2 × 6
modelo r.squared adj.r.squared AIC BIC df
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 PovPct 0.533 0.524 324. 329. 1
2 PovPct + ViolCrime 0.647 0.632 311. 319. 2El modelo con las dos variables tiene AIC y BIC más chicos, así que gana.
12.4.1 Prueba anidada con anova
Cuando un modelo está contenido en el otro (le llamamos anidado: son las mismas variables más otras), hay una prueba formal para preguntar si las variables extra aportan algo:
anova(modelo_clasico, modelo_multiple)Analysis of Variance Table
Model 1: Brth15to17 ~ PovPct
Model 2: Brth15to17 ~ PovPct + ViolCrime
Res.Df RSS Df Sum of Sq F Pr(>F)
1 49 1509.6
2 48 1142.7 1 366.94 15.413 0.0002753 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El valor p chiquito nos dice que agregar ViolCrime sí mejora el ajuste más de lo que se esperaría por azar.
AdvertenciaTres cosas que hay que tener claras del AIC
- El AIC no dice si tu modelo es bueno, sólo cuál de los que le diste es menos malo. Si todos tus modelos están mal, el AIC te va a elegir uno con toda confianza.
- Sólo se pueden comparar modelos ajustados a exactamente los mismos datos. Si un modelo usa una variable con datos faltantes,
Rva a tirar esas filas y los dos modelos ya no serán comparables. Es un error silencioso y frecuente. - Los números del AIC no significan nada por sí solos. Un AIC de
311no es bueno ni malo; sólo sirve comparado con el323del otro modelo. No los reportes sueltos.
12.5 Multicolinealidad y el VIF
Aquí viene el problema que más confunde en regresión múltiple, y lo vamos a ver con un ejemplo donde se nota muchísimo.
La base tiene tres variables de natalidad: Brth15to17 (tasa en adolescentes de 15 a 17 años), Brth18to19 (de 18 a 19) y TeenBrth. Por cómo están definidas, la tercera contiene información de las dos primeras. Metámoslas todas al mismo modelo y veamos qué pasa:
modelo_colineal <- lm(Brth15to17 ~ PovPct + Brth18to19 + TeenBrth, data = emb_pob)
tidy(modelo_colineal)# A tibble: 4 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -0.583 0.644 -0.906 3.69e- 1
2 PovPct -0.0534 0.0464 -1.15 2.56e- 1
3 Brth18to19 -0.520 0.0502 -10.4 1.04e-13
4 TeenBrth 1.44 0.0827 17.5 3.56e-22Y ahora fíjate en las métricas de ajuste:
glance(modelo_colineal) %>% select(r.squared, adj.r.squared, AIC, BIC)# A tibble: 1 × 4
r.squared adj.r.squared AIC BIC
<dbl> <dbl> <dbl> <dbl>
1 0.988 0.988 140. 149.
ImportanteUn modelo espectacular y completamente inservible
Compara con los modelos anteriores: la \(R^2\) pasó de \(0.53\) a \(0.99\) y el AIC se derrumbó de \(323\) a \(140\). Por cualquier criterio de los que acabamos de ver, éste es por mucho el mejor modelo.
Y sin embargo mira los coeficientes:
- El de
PovPctcambió de signo: pasó de \(+1.37\) a \(-0.05\), y dejó de ser distinguible de cero. ¿La pobreza ahora reduce la natalidad adolescente? - El de
Brth18to19es negativo (\(-0.52\)): ¿más nacimientos en las de 18 a 19 significan menos nacimientos en las de 15 a 17?
Ninguna de las dos cosas tiene el menor sentido. Éste es el mejor ejemplo posible de lo que vimos en el capítulo de modelos: un modelo puede ser excelente para predecir y una fuente de conclusiones falsas si lo lees como si fuera de inferencia.
¿Por qué pasó esto? Porque las variables se explican unas a otras. Cuando dos variables cargan casi la misma información, el modelo no puede decidir a cuál darle el crédito, así que reparte los coeficientes de forma inestable: uno gigante positivo, otro gigante negativo, y las predicciones salen bien de milagro.
12.5.1 Cómo se detecta: el VIF
El VIF (Variance Inflation Factor, factor de inflación de la varianza) mide exactamente esto: para cada variable, cuánto se infla su varianza por culpa de estar correlacionada con las demás. Se calcula con vif del paquete car:
# Nuestro modelo bueno
vif(modelo_multiple) PovPct ViolCrime
1.282857 1.282857 # El modelo con el problema
vif(modelo_colineal) PovPct Brth18to19 TeenBrth
2.442371 56.289092 64.369636
TipCómo leer el VIF
El VIF vale como mínimo 1 (esa variable no comparte nada con las otras) y no tiene máximo.
| VIF | Interpretación |
|---|---|
| \(\approx 1\) | Sin problema. |
| \(< 5\) | Tolerable, es lo normal en datos reales. |
| \(5\) a \(10\) | Sospechoso, vale la pena revisar. |
| \(> 10\) | Problema serio de multicolinealidad. |
En nuestro modelo bueno los dos VIF son \(1.28\): perfecto. En el modelo colineal son \(56\) y \(64\), o sea cinco a seis veces arriba del umbral de alarma.
Hay una regla mnemotécnica útil: \(\sqrt{\text{VIF}}\) te dice por cuánto se multiplicó el error estándar de ese coeficiente. Un VIF de \(64\) significa que el error estándar es 8 veces más grande de lo que sería sin colinealidad.
12.5.2 ¿Y qué se hace al respecto?
La multicolinealidad no se “arregla” con un comando; se resuelve pensando:
- Quita variables redundantes. Si
TeenBrthya contiene a las otras dos, escoge una y ya. Es la solución correcta el 90% de las veces. - Combínalas en un solo indicador si tiene sentido conceptual.
- No hagas nada — si tu modelo es puramente de predicción y no vas a interpretar ningún coeficiente, la colinealidad no te estorba. Las predicciones del modelo colineal son de hecho buenísimas.
Advertencia
La decisión de qué variable quitar es conceptual, no estadística. No dejes que el VIF elija por ti: elige tú con base en qué significa cada variable y qué pregunta estás contestando.
12.5.3 Ejercicio (VIF y AIC)
Ajusta el modelo
Brth15to17 ~ PovPct + TeenBrth(sólo dos predictores) y calcula suVIFy suAIC. ¿Sigue habiendo problema de colinealidad? ¿Cómo se compara elAICcon el del modelo de tres variables?Calcula la matriz de correlaciones de todas las variables numéricas con
cor(emb_pob[,-1]). Antes de verla, predice cuáles dos variables van a tener la correlación más alta. ¿Le atinaste?Si tuvieras que reportar en un informe un solo modelo para contestar “¿la pobreza está asociada a la natalidad adolescente?”, ¿cuál de los tres elegirías y por qué? Nota que no es el del mejor AIC. Escribe tu justificación en tres renglones.
NotaRespuesta del punto 3
El de Brth15to17 ~ PovPct, o cuando mucho el que agrega ViolCrime si pudieras justificar conceptualmente por qué ajustas por crimen violento.
El modelo de tres variables tiene el mejor AIC y es inutilizable para esta pregunta: ajustar por Brth18to19 y TeenBrth es ajustar por variables que son prácticamente el mismo desenlace que estás midiendo. El coeficiente de pobreza que sale de ahí no contesta la pregunta que te hicieron.
Es el punto central del capítulo de modelos: primero decides qué pregunta tienes, y sólo entonces eliges el criterio. El AIC es una herramienta para elegir entre modelos que ya tienen sentido, no un sustituto de pensar.
12.6 Series de tiempo: forecast clásico y modeltime
Las regresiones que hemos visto no tienen tiempo: cada fila es un estado y el orden no importa. Cuando los datos son una serie de tiempo (una medición por día, por semana, por año) el orden sí importa y hacen falta otras herramientas.
Aquí también hay dos dialectos: forecast (el clásico, de Rob Hyndman) y modeltime (que es la versión tidymodels). Usaremos los casos diarios de enfermedades respiratorias del SINAVE:
serie_covid <- read_csv("datasets/casos_covid_agosto_2022.csv") %>%
rename(fecha = FECHA_SINTOMAS, casos = n) %>%
arrange(fecha)
# Así se ven los datos
head(serie_covid)# A tibble: 6 × 2
fecha casos
<date> <dbl>
1 2020-01-01 287
2 2020-01-02 233
3 2020-01-03 247
4 2020-01-04 245
5 2020-01-05 349
6 2020-01-06 345ggplot(serie_covid) +
geom_line(aes(x = fecha, y = casos), color = "#003f5c") +
labs(x = "", y = "Casos") +
scale_y_continuous(labels = scales::label_comma()) +
theme_bw()
12.6.1 La forma clásica: forecast
En el mundo clásico, una serie de tiempo se guarda en un objeto especial (ts) donde le dices cada cuántas observaciones se repite el patrón. Como son datos diarios con patrón semanal, la frecuencia es 7:
# frequency = 7 porque el patrón se repite cada semana
serie_ts <- ts(serie_covid$casos, frequency = 7)
# auto.arima elige solo el tipo de modelo ARIMA
modelo_arima <- auto.arima(serie_ts)
modelo_arimaSeries: serie_ts
ARIMA(2,0,2)(2,1,0)[7]
Coefficients:
ar1 ar2 ma1 ma2 sar1 sar2
1.9146 -0.9276 -1.4958 0.6152 -0.5220 -0.2266
s.e. 0.0178 0.0174 0.0318 0.0320 0.0329 0.0325
sigma^2 = 11627635: log likelihood = -9074.54
AIC=18163.09 AICc=18163.21 BIC=18197.08Para pronosticar las siguientes dos semanas:
pronostico <- forecast(modelo_arima, h = 14)
autoplot(pronostico) +
labs(x = "Semanas", y = "Casos") +
theme_bw()
AdvertenciaRevisa el paso 3 del flujo de trabajo
Mira los valores que pronosticó el modelo:
round(head(as.numeric(pronostico$mean), 5), 1)[1] 7785.5 2058.5 -10.9 -2781.3 -3821.2Hay un pronóstico negativo. El modelo está prediciendo un número negativo de casos de una enfermedad, que es imposible.
Esto no es un error de programación: es que un ARIMA supone que los datos son continuos y pueden tomar cualquier valor, y nosotros le dimos conteos que nunca pueden bajar de cero. Es exactamente el problema del paso 3 del flujo de trabajo: si hubieras simulado del modelo antes de ajustarlo, lo habrías visto.
La solución honesta es modelar los conteos con algo apropiado (por ejemplo un modelo Poisson o binomial negativo, o modelar el logaritmo de los casos), no truncar los negativos a cero y seguir como si nada.
12.6.2 La forma tidymodels: modeltime
modeltime hace lo mismo con la sintaxis de tidymodels. La diferencia grande es que te obliga (para bien) a partir los datos en entrenamiento y prueba, para que evalúes el pronóstico contra datos que el modelo no vio:
# Apartamos los últimos 14 días para evaluar
particion <- time_series_split(
serie_covid,
date_var = fecha,
assess = "14 days",
cumulative = TRUE
)
# Mismo auto.arima, sintaxis tidymodels
modelo_arima_tm <- arima_reg() %>%
set_engine("auto_arima") %>%
fit(casos ~ fecha, data = training(particion))Ahora lo calibramos contra los datos que apartamos y medimos qué tan bien le fue:
tabla_modelos <- modeltime_table(modelo_arima_tm) %>%
modeltime_calibrate(new_data = testing(particion))
tabla_modelos %>%
modeltime_accuracy() %>%
select(.model_desc, mae, rmse, rsq)# A tibble: 1 × 4
.model_desc mae rmse rsq
<chr> <dbl> <dbl> <dbl>
1 ARIMA(2,0,2)(2,1,0)[7] 6570. 8599. 0.295
TipAquí sí gana
tidymodels
Nota lo que acaba de pasar: la ruta clásica te dio un pronóstico bonito y nunca te dijo si era bueno. La ruta modeltime te obligó a apartar datos y te dio un error medio absoluto (mae) medido contra datos reales que el modelo no vio.
Ésa es la diferencia entre el paso 4 y el paso 7 del flujo de trabajo. Y en predicción, saltarse el paso 7 es la falla más grave que puedes cometer: un pronóstico sin evaluar no es un pronóstico, es un dibujo.
Se puede hacer perfectamente con forecast a mano (partiendo la serie tú), pero modeltime no te deja olvidarlo.
12.6.3 Ejercicio (series de tiempo)
Cambia el
assess = "14 days"porassess = "60 days"y vuelve a calcular la exactitud. ¿Mejoró o empeoró elmae? ¿Por qué tiene sentido?Ajusta la serie con
frequency = 365en lugar de7en la ruta clásica. ¿Cambió el modelo que eligióauto.arima? ¿Qué patrón estás suponiendo con cada frecuencia?A mano. Mira la gráfica de la serie completa y dibuja en tu cuaderno cómo crees que debería verse un pronóstico razonable a 14 días. Después compáralo con lo que produjo el modelo. ¿El modelo captó la tendencia que tú viste?
Ajusta un segundo modelo con
exp_smoothing()(suavizamiento exponencial), agrégalo a lamodeltime_tablejunto alARIMAy compara los dos conmodeltime_accuracy(). ¿Cuál pronostica mejor?
12.7 Ejercicio
- Corre el siguiente código para generar una base de datos de nombre
datLongque contiene \(4\) grupos. Para cada grupo genere una regresión lineal de la forma:
\[ y = \beta_0 + \beta_1 x \]
Identifica cuáles regresiones sí ajustan bien y cuáles no mediante los gráficos de diagnóstico. Finalmente grafica tus datos \(x\) contra \(y\) y la regresión para ver que lo hayas hecho bien. ¿Hay alguna forma de corregir alguna de las que no ajusta bien?
dat <- datasets::anscombe
datLong <- data.frame(
grupo = rep(1:4, each = 11),
x = unlist(dat[, c(1:4)]),
y = unlist(dat[, c(5:8)])
)
rownames(datLong) <- NULLLee la base de datos de \(n = 345\) niños entre \(6\) y \(10\) años de Kahn, Michael (2005). Las variables de interés son \(y = \text{FEV}\) el volumen de expiración forzada y \(x = \text{edad}\) en años. Realiza una regresión lineal. Justifica que no se cumple la homocedasticidad mediante una gráfica de escala locación.
Plantee una regresión lineal usando los datos de esta liga para determinar si el sexo influye en el salario. Ojo en la regresión es necesario incluir otras covariables.
12.8 ¿Y si lo hacemos bayesiano?
Para hacer la misma regresión lineal pero con estadística bayesiana podemos nada más cambiar el engine:
# Ajusta un modelo lineal
# Brth15to17 = intercepto + pendiente*PovPct
modelo_bayesiano <- linear_reg() %>%
set_engine("stan") %>%
fit(Brth15to17 ~ PovPct, data = emb_pob) # Notación y ~ xY podemos ver el ajuste:
modelo_bayesiano %>%
extract_fit_engine() %>%
summary()
Model Info:
function: stan_glm
family: gaussian [identity]
formula: Brth15to17 ~ PovPct
algorithm: sampling
sample: 4000 (posterior sample size)
priors: see help('prior_summary')
observations: 51
predictors: 2
Estimates:
mean sd 10% 50% 90%
(Intercept) 4.2 2.6 0.9 4.2 7.5
PovPct 1.4 0.2 1.1 1.4 1.6
sigma 5.6 0.6 4.9 5.6 6.4
Fit Diagnostics:
mean sd 10% 50% 90%
mean_PPD 22.3 1.1 20.9 22.3 23.7
The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
MCMC diagnostics
mcse Rhat n_eff
(Intercept) 0.0 1.0 3652
PovPct 0.0 1.0 3604
sigma 0.0 1.0 3679
mean_PPD 0.0 1.0 4063
log-posterior 0.0 1.0 1558
For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).o realizar predicciones:
predichos <- modelo_bayesiano %>%
extract_fit_engine() %>%
predict(new_data = emb_pob)
ic <- modelo_bayesiano %>%
extract_fit_engine() %>%
predictive_interval(newdata = emb_pob, prob = 0.95) # Para bayesiana
# Juntamos los predichos con los observados
obs_y_modelo <- emb_pob %>%
cbind(predichos) %>%
cbind(ic)
# Graficamos
ggplot(obs_y_modelo) +
geom_ribbon(aes(x = PovPct, ymin = `2.5%`, ymax = `97.5%`),
fill = "#003f5c", linewidth = 1, alpha = 0.5
) +
geom_point(aes(x = PovPct, y = Brth15to17), color = "#ffa600", size = 3) +
geom_line(aes(x = PovPct, y = predichos),
color = "#003f5c", linewidth = 1
) +
labs(
x = "Porcentaje en pobreza",
y = "Tasa bruta de natalidad (por cada 1,000 mujeres adolescentes)",
title = "Relación entre tasa bruta de natalidad en adolescentes\nde 15 a 19 años y porcentaje en pobreza"
) +
theme_bw()
Para validación del modelo puedes checar esta página que explica loo (ver paper:
modelo_bayesiano %>%
extract_fit_engine() %>%
loo()
Computed from 4000 by 51 log-likelihood matrix.
Estimate SE
elpd_loo -161.6 4.2
p_loo 2.4 0.5
looic 323.2 8.4
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume independent draws (r_eff=1).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.asi como su visualización (si el modelo no fuera bueno)
modelo_bayesiano %>%
extract_fit_engine() %>%
loo() %>%
plot(label_points = TRUE)
12.9 Ejercicio
- Utilice las opciones de
enginepara cambiar elprior_interceptdel intercepto a unatde Student y elpriorde los coeficientes a unaLaplace. ¿Cambia mucho el resultado?
13 Continuación
14 Sistema
sessioninfo::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.5.3 (2026-03-11)
os macOS Sequoia 15.7.3
system x86_64, darwin20
ui X11
language (EN)
collate es_ES.UTF-8
ctype es_ES.UTF-8
tz America/Mexico_City
date 2026-08-05
pandoc 3.10.1 @ /usr/local/bin/ (via rmarkdown)
quarto 1.10.18 @ /usr/local/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
abind 1.4-8 2024-09-12 [1] CRAN (R 4.5.0)
backports 1.5.1 2026-04-03 [1] CRAN (R 4.5.2)
base64enc 0.1-6 2026-02-02 [1] CRAN (R 4.5.1)
bayesplot 1.15.0.9000 2026-04-17 [1] https://stan-dev.r-universe.dev (R 4.5.3)
bit 4.6.0 2025-03-06 [1] CRAN (R 4.5.0)
bit64 4.8.2 2026-05-19 [1] CRAN (R 4.5.2)
boot 1.3-32 2025-08-29 [1] CRAN (R 4.5.3)
broom * 1.0.12 2026-01-27 [1] CRAN (R 4.5.1)
cachem 1.1.0 2024-05-16 [1] CRAN (R 4.5.0)
car * 3.1-5 2026-02-03 [1] CRAN (R 4.5.1)
carData * 3.0-6 2026-01-30 [1] CRAN (R 4.5.1)
checkmate 2.3.4 2026-02-03 [1] CRAN (R 4.5.1)
class 7.3-23 2025-01-01 [1] CRAN (R 4.5.3)
cli 3.6.6 2026-04-09 [1] CRAN (R 4.5.2)
codetools 0.2-20 2024-03-31 [1] CRAN (R 4.5.3)
colorspace 2.1-2 2025-09-22 [1] CRAN (R 4.5.1)
colourpicker 1.3.0 2023-08-21 [1] CRAN (R 4.5.0)
conflicted 1.2.0 2023-02-01 [1] CRAN (R 4.5.0)
cowplot * 1.2.0 2025-07-07 [1] CRAN (R 4.5.1)
crayon 1.5.3 2024-06-20 [1] CRAN (R 4.5.0)
crosstalk 1.2.2 2025-08-26 [1] CRAN (R 4.5.1)
curl 7.1.0 2026-04-22 [1] CRAN (R 4.5.2)
data.table 1.18.4 2026-05-06 [1] CRAN (R 4.5.2)
dials * 1.4.2 2025-09-04 [1] CRAN (R 4.5.1)
DiceDesign 1.10 2023-12-07 [1] CRAN (R 4.5.0)
digest 0.6.39 2025-11-19 [1] CRAN (R 4.5.1)
distributional 0.8.1 2026-06-27 [1] CRAN (R 4.5.2)
dplyr * 1.2.1 2026-04-03 [1] CRAN (R 4.5.2)
DT 0.34.0 2025-09-02 [1] CRAN (R 4.5.1)
dygraphs 1.1.1.6 2018-07-11 [1] CRAN (R 4.5.0)
evaluate 1.0.5 2025-08-27 [1] CRAN (R 4.5.1)
farver 2.1.2 2024-05-13 [1] CRAN (R 4.5.0)
fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.5.0)
forcats * 1.0.1 2025-09-25 [1] CRAN (R 4.5.1)
forecast * 9.0.2 2026-03-18 [1] CRAN (R 4.5.2)
Formula 1.2-5 2023-02-24 [1] CRAN (R 4.5.0)
fracdiff 1.5-3 2024-02-01 [1] CRAN (R 4.5.0)
furrr 0.3.1 2022-08-15 [1] CRAN (R 4.5.0)
future 1.70.0 2026-03-14 [1] CRAN (R 4.5.2)
future.apply 1.20.2 2026-02-20 [1] CRAN (R 4.5.2)
generics 0.1.4 2025-05-09 [1] CRAN (R 4.5.0)
ggfortify * 0.4.19 2025-07-27 [1] CRAN (R 4.5.1)
ggplot2 * 4.0.3 2026-04-22 [1] CRAN (R 4.5.2)
ggridges * 0.5.7 2025-08-27 [1] CRAN (R 4.5.1)
globals 0.19.1 2026-03-13 [1] CRAN (R 4.5.2)
glue 1.8.1 2026-04-17 [1] CRAN (R 4.5.2)
gower 1.0.2 2024-12-17 [1] CRAN (R 4.5.0)
GPfit 1.0-9 2025-04-12 [1] CRAN (R 4.5.0)
gridExtra 2.3 2017-09-09 [1] CRAN (R 4.5.0)
gtable 0.3.6 2024-10-25 [1] CRAN (R 4.5.0)
gtools 3.9.5 2023-11-20 [1] CRAN (R 4.5.0)
hardhat 1.4.2 2025-08-20 [1] CRAN (R 4.5.1)
hms 1.1.4 2025-10-17 [1] CRAN (R 4.5.1)
htmltools 0.5.9 2025-12-04 [1] CRAN (R 4.5.1)
htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.5.0)
httpuv 1.6.17 2026-03-18 [1] CRAN (R 4.5.2)
igraph 2.2.2 2026-02-12 [1] CRAN (R 4.5.2)
infer * 1.1.0 2025-12-18 [1] CRAN (R 4.5.1)
inline 0.3.21 2025-01-09 [1] CRAN (R 4.5.0)
ipred 0.9-15 2024-07-18 [1] CRAN (R 4.5.0)
janitor 2.2.1 2024-12-22 [1] CRAN (R 4.5.0)
jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.5.0)
knitr 1.51 2025-12-20 [1] CRAN (R 4.5.1)
labeling 0.4.3 2023-08-29 [1] CRAN (R 4.5.0)
later 1.4.8 2026-03-05 [1] CRAN (R 4.5.2)
lattice 0.22-9 2026-02-09 [1] CRAN (R 4.5.3)
lava 1.8.2 2025-10-30 [1] CRAN (R 4.5.1)
lhs 1.2.1 2026-03-01 [1] CRAN (R 4.5.2)
lifecycle 1.0.5 2026-01-08 [1] CRAN (R 4.5.1)
listenv 0.10.1 2026-03-10 [1] CRAN (R 4.5.2)
lme4 2.0-6 2026-07-16 [1] CRAN (R 4.5.2)
loo 2.10.0.9000 2026-07-10 [1] https://stan-dev.r-universe.dev (R 4.5.3)
lubridate * 1.9.5 2026-02-04 [1] CRAN (R 4.5.1)
magrittr 2.0.5 2026-04-04 [1] CRAN (R 4.5.2)
markdown 2.0 2025-03-23 [1] CRAN (R 4.5.0)
MASS 7.3-65 2025-02-28 [1] CRAN (R 4.5.3)
Matrix 1.7-5 2026-03-21 [1] CRAN (R 4.5.2)
matrixStats 1.5.0 2025-01-07 [1] CRAN (R 4.5.0)
memoise 2.0.1 2021-11-26 [1] CRAN (R 4.5.0)
mgcv 1.9-4 2025-11-07 [1] CRAN (R 4.5.3)
mime 0.13 2025-03-17 [1] CRAN (R 4.5.0)
miniUI 0.1.2 2025-04-17 [1] CRAN (R 4.5.0)
minqa 1.2.8 2024-08-17 [1] CRAN (R 4.5.0)
modeldata * 1.5.1 2025-08-22 [1] CRAN (R 4.5.1)
modeltime * 1.3.5 2026-01-31 [1] CRAN (R 4.5.1)
nlme 3.1-168 2025-03-31 [1] CRAN (R 4.5.3)
nloptr 2.2.1 2025-03-17 [1] CRAN (R 4.5.0)
nnet 7.3-20 2025-01-01 [1] CRAN (R 4.5.3)
otel 0.2.0 2025-08-29 [1] CRAN (R 4.5.1)
parallelly 1.47.0 2026-04-17 [1] CRAN (R 4.5.2)
parsnip * 1.4.1 2026-01-11 [1] CRAN (R 4.5.1)
pillar 1.11.1 2025-09-17 [1] CRAN (R 4.5.1)
pkgbuild 1.4.8 2025-05-26 [1] CRAN (R 4.5.0)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.5.0)
plyr 1.8.9 2023-10-02 [1] CRAN (R 4.5.0)
poissonreg * 1.0.2 2026-04-20 [1] CRAN (R 4.5.2)
posterior 1.7.1 2026-07-15 [1] https://stan-dev.r-universe.dev (R 4.5.3)
prodlim 2026.03.11 2026-03-11 [1] CRAN (R 4.5.2)
promises 1.5.0 2025-11-01 [1] CRAN (R 4.5.1)
purrr * 1.2.2 2026-04-10 [1] CRAN (R 4.5.2)
QuickJSR 1.10.0 2026-05-17 [1] CRAN (R 4.5.2)
R6 2.6.1 2025-02-15 [1] CRAN (R 4.5.0)
rbibutils 2.4.1 2026-01-21 [1] CRAN (R 4.5.1)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.5.0)
Rcpp * 1.1.2 2026-07-05 [1] CRAN (R 4.5.2)
RcppParallel 5.1.11-2 2026-03-05 [1] CRAN (R 4.5.2)
Rdpack 2.6.6 2026-02-08 [1] CRAN (R 4.5.1)
readr * 2.2.0 2026-02-19 [1] CRAN (R 4.5.2)
recipes * 1.3.1 2025-05-21 [1] CRAN (R 4.5.0)
reformulas 0.4.4 2026-02-02 [1] CRAN (R 4.5.1)
reshape2 1.4.5 2025-11-12 [1] CRAN (R 4.5.1)
rlang 1.3.0 2026-07-05 [1] CRAN (R 4.5.2)
rmarkdown 2.31 2026-03-26 [1] CRAN (R 4.5.2)
rpart 4.1.24 2025-01-07 [1] CRAN (R 4.5.3)
rsample * 1.3.2 2026-01-30 [1] CRAN (R 4.5.1)
rstan 2.36.0.9000 2026-05-22 [1] https://stan-dev.r-universe.dev (R 4.5.3)
rstanarm * 2.32.2 2025-09-30 [1] CRAN (R 4.5.1)
rstantools 2.6.0.9000 2026-01-10 [1] https://stan-dev.r-universe.dev (R 4.5.3)
rstudioapi 0.19.0 2026-06-11 [1] CRAN (R 4.5.2)
S7 0.2.2 2026-04-22 [1] CRAN (R 4.5.2)
scales * 1.4.0 2025-04-24 [1] CRAN (R 4.5.0)
sessioninfo 1.2.3 2025-02-05 [1] CRAN (R 4.5.0)
shiny 1.13.0 2026-02-20 [1] CRAN (R 4.5.2)
shinyjs 2.1.1 2026-01-15 [1] CRAN (R 4.5.1)
shinystan 2.7.0 2025-12-12 [1] CRAN (R 4.5.1)
shinythemes 1.2.0 2021-01-25 [1] CRAN (R 4.5.0)
snakecase 0.11.1 2023-08-27 [1] CRAN (R 4.5.0)
sparsevctrs 0.3.6 2026-01-27 [1] CRAN (R 4.5.1)
StanHeaders 2.36.0.9000 2026-05-22 [1] https://stan-dev.r-universe.dev (R 4.5.3)
stringi 1.8.7 2025-03-27 [1] CRAN (R 4.5.0)
stringr * 1.6.0 2025-11-04 [1] CRAN (R 4.5.1)
survival 3.8-6 2026-01-16 [1] CRAN (R 4.5.3)
tailor * 0.1.0 2025-08-25 [1] CRAN (R 4.5.1)
tensorA 0.36.2.1 2023-12-13 [1] CRAN (R 4.5.0)
threejs 0.3.4 2025-04-21 [1] CRAN (R 4.5.0)
tibble * 3.3.1 2026-01-11 [1] CRAN (R 4.5.1)
tidymodels * 1.4.1 2025-09-08 [1] CRAN (R 4.5.1)
tidyr * 1.3.2 2025-12-19 [1] CRAN (R 4.5.1)
tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.5.0)
tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.5.0)
timechange 0.4.0 2026-01-29 [1] CRAN (R 4.5.1)
timeDate 4052.112 2026-01-28 [1] CRAN (R 4.5.1)
timetk * 2.9.1 2025-08-29 [1] CRAN (R 4.5.1)
tune * 2.0.1 2025-10-17 [1] CRAN (R 4.5.1)
tzdb 0.5.0 2025-03-15 [1] CRAN (R 4.5.0)
urca 1.3-4 2024-05-27 [1] CRAN (R 4.5.0)
utf8 1.2.6 2025-06-08 [1] CRAN (R 4.5.0)
V8 8.2.0 2026-04-21 [1] CRAN (R 4.5.2)
vctrs 0.7.3 2026-04-11 [1] CRAN (R 4.5.2)
vroom 1.7.1 2026-03-31 [1] CRAN (R 4.5.2)
withr 3.0.3 2026-06-19 [1] CRAN (R 4.5.2)
workflows * 1.3.0 2025-08-27 [1] CRAN (R 4.5.1)
workflowsets * 1.1.1 2025-05-27 [1] CRAN (R 4.5.0)
xfun 0.60 2026-07-09 [1] CRAN (R 4.5.2)
xtable 1.8-8 2026-02-22 [1] CRAN (R 4.5.2)
xts 0.14.2 2026-02-28 [1] CRAN (R 4.5.2)
yaml 2.3.12 2025-12-10 [1] CRAN (R 4.5.1)
yardstick * 1.3.2 2025-01-22 [1] CRAN (R 4.5.0)
zoo 1.8-15 2025-12-15 [1] CRAN (R 4.5.1)
[1] /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/library
* ── Packages attached to the search path.
──────────────────────────────────────────────────────────────────────────────