knitr::opts_chunk$set(warning = FALSE, message = FALSE)
# Cargamos las librerías que vamos a utilizar
library(tidyverse)
library(tidymodels)
library(modelr)
library(GGally)
library(pROC)
library(cowplot)
library(OneR)
library(rlang)
library(caret)
set.seed(15)Un GLM es un modelo lineal para la esperanza transformada de una variable respuesta que tiene una distribución de probabilidad de la familia exponencial. Esto permite limitar los valores de Y al rango esperado para la variable (por ejemplo valores mayores o iguales a 0 para conteos o entre 0 y 1 para proporciones). La estimación de los parámetros se hace por máxima verosimilitud e iteraciones. Son modelos lineales ya que la variable respuesta es descripta a través de una función de enlace por una combinación lineal de los predictores X.
Los 3 componentes
Un componente aleatorio que identifica a la variable respuesta Y y sigue una de las distribuciones de la familia exponencial.
Un componente o predictor lineal que identifica a un conjunto de predictoras X en una función lineal.
Una función de enlace, que describe cómo se relacionan el valor esperado o medio de la variable respuesta con el componente lineal. Debe restringir los valores posibles de la variable respuesta a su rango teórico.
mieleros <- read.csv("nminer.csv")
glimpse(mieleros)Rows: 31
Columns: 8
$ Miners <int> 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0,…
$ Eucs <int> 2, 10, 16, 20, 19, 18, 12, 16, 3, 12, 32, 2, 16, 7, 10, 15, 30…
$ Area <int> 22, 11, 51, 22, 4, 61, 16, 14, 5, 6, 3, 13, 24, 32, 22, 27, 18…
$ Grazed <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1,…
$ Shrubs <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,…
$ Bulokes <int> 120, 67, 85, 45, 160, 75, 100, 321, 275, 227, 23, 277, 243, 11…
$ Timber <int> 16, 25, 13, 12, 14, 6, 12, 15, 8, 10, 61, 22, 25, 9, 35, 25, 1…
$ Minerab <int> 0, 0, 3, 2, 8, 1, 8, 5, 0, 4, 19, 0, 2, 0, 3, 1, 7, 1, 0, 7, 0…ggplot(mieleros, aes(y = Minerab, x = Eucs)) +
geom_point(colour=2)+
xlab("Cant. de Eucaliptos cada 2 Ha") +
ylab("Abundancia de mieleros")+
theme_minimal()
Propongo el modelo: \[Y \sim Poisson(\mu) \] con \[ \mu_i = e^{\beta_0 + \beta_1 x_i} \] o \[\log(\mu) = \beta_0 + \beta_1X \] \[E(Y_i) = Var(Y_i) = \lambda_i = \mu_i \]
Modelo en la escala del predictor lineal
\[log(E(Y_i)) = \beta_0 + \beta_1X_i \]
Modelo en la escala de la variable respuesta
\[E(Y_i) = \mu_i = e^{\beta_0 + \beta_1X_i} \]
La funcíón glm() nos permite crear un modelo lineal generalizado. Al igual que la función lm() toma como argumentos una fórmula y los datos pero también se debe especificar el argumento family donde indicamos la distribución de probabilidades y la función link que vamos a utilizar en el modelo.
Algunos ejemplos:
Binomial: link=logit
Poisson: link=log
Gaussiana: link=identidad
library(lme4)# glmmTMB si se quiere trabajar con modelos mixtos
poisson <- glm(Minerab~ Eucs,
data=mieleros,
family = "poisson")
summary(poisson)
Call:
glm(formula = Minerab ~ Eucs, family = "poisson", data = mieleros)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.87621 0.28279 -3.098 0.00195 **
Eucs 0.11398 0.01243 9.169 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 150.545 on 30 degrees of freedom
Residual deviance: 63.318 on 29 degrees of freedom
AIC: 121.47
Number of Fisher Scoring iterations: 5El modelo ajustado es:
\[ \log(\mu)=\beta_0+\beta_1 x \]
donde:
Con las estimaciones obtenidas:
\[ \log(\mu)= -0.876 + 0.114x \]
Interpretación de \(\beta_0\):
\[ \beta_0=-0.876 \]
representa el logaritmo de la cantidad media esperada de mieleros cuando \(x=0\), es decir, cuando no hay eucaliptos.
Aplicando exponencial:
\[ \mu(0)=e^{-0.876}\approx 0.42 \]
Por lo tanto, el modelo estima aproximadamente 0.42 mieleros cuando no hay eucaliptos.
Interpretación de \(\beta_1\):
\[ \beta_1=0.114 \]
Como el vínculo es logarítmico, el efecto es multiplicativo:
\[ e^{0.114}\approx 1.12 \]
Esto significa que por cada eucalipto adicional,la cantidad esperada de mieleros aumenta aproximadamente un 12%.
Independencia entre observaciones
Linealidad para las variables explicatorias cuantitativas a escala del predictor lineal.
Ausencia de sub o sobredispersión. La varianza de la variable es mayor o menor, respectivamente, a la que se espera según la distribucion de Poisson. Esto provoca estimaciones sesgadas de los errores estándar de los coeficientes de las variables explicatorias.
AER::dispersiontest(poisson)
Overdispersion test
data: poisson
z = 2.073, p-value = 0.01909
alternative hypothesis: true dispersion is greater than 1
sample estimates:
dispersion
2.337731 Ignorar la sobredispersión hace que haya mas chances de tener falsos positivos (encontrar diferencias que no son reales). Si hay subdispersión hay mas chances de tener falsos negativos.
¿Cómo tratar la sobredispersión?
Se debe ajustar a otras distribuciones de probabilidades para la variable respuesta: Binomial Negativa, Conway-Maxwell-Poisson, Tweedie.
¿ Cómo tratar la subdispersión?
Se debe ajustar a otras distribuciones de probabilidades para la variable respuesta: Conway-Maxwell-Poisson, Tweedie, Poisson generalizada (family=“genpois”)
Modelo estimado
ggplot(mieleros, aes(x = Eucs, y = Minerab)) +
geom_point(color = "red", alpha = 0.7, size = 2) + # Puntos observados
stat_smooth(
method = "glm",
method.args = list(family = "poisson"),
se = TRUE, # Muestra el intervalo de confianza (cambiar a FALSE si no lo deseas)
color = "violet", # Color de la curva estimada
size = 1
) +
labs(
title = "GLM Poisson",
x = "Cantidad de Eucaliptos cada 2 Ha",
y = "Abundancia de mieleros"
) +
theme_minimal()
Vamos a trabajar con una base de datos extraída de Kaggle. Es una base de datos que consiste de 69301 registros de pacientes a los que se les registró 12 caracteristicas asociadas a una enfermedad cardiovascular.
Vamos a utilizar diferentes modelos de regresión logística para predecir la presencia o ausencia de enfermedad cardiovascular (ECV) utilizando los resultados del examen del paciente.
Descripción de los datos
df <- read.csv("cardio_train.csv", sep = ";")
colSums(is.na(df))# chequeamos si hay datos faltantes id age gender height weight ap_hi
0 0 0 0 0 0
ap_lo cholesterol gluc smoke alco active
0 0 0 0 0 0
cardio
0 sum(duplicated(df))#chequeamos duplicados[1] 0glimpse(df)Rows: 70,000
Columns: 13
$ id <int> 0, 1, 2, 3, 4, 8, 9, 12, 13, 14, 15, 16, 18, 21, 23, 24, 2…
$ age <int> 18393, 20228, 18857, 17623, 17474, 21914, 22113, 22584, 17…
$ gender <int> 2, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 2, 1, 2, 2, 1, 1, 1, 2…
$ height <int> 168, 156, 165, 169, 156, 151, 157, 178, 158, 164, 169, 173…
$ weight <dbl> 62, 85, 64, 82, 56, 67, 93, 95, 71, 68, 80, 60, 60, 78, 95…
$ ap_hi <int> 110, 140, 130, 150, 100, 120, 130, 130, 110, 110, 120, 120…
$ ap_lo <int> 80, 90, 70, 100, 60, 80, 80, 90, 70, 60, 80, 80, 80, 70, 9…
$ cholesterol <int> 1, 3, 3, 1, 1, 2, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ gluc <int> 1, 1, 1, 1, 1, 2, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1…
$ smoke <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1…
$ alco <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0…
$ active <int> 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1…
$ cardio <int> 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0…Se genera una muestra de tamaño 1000 con la que se realizaran los análisis y se analiza si hay desbalance en la variable respuesta.
df <- df %>%
sample_n(1000)
df %>%
group_by(cardio) %>%
summarise(cnt = n()) %>%
mutate(freq = round(cnt / sum(cnt), 2))# A tibble: 2 × 3
cardio cnt freq
<int> <int> <dbl>
1 0 512 0.51
2 1 488 0.49df_split <- initial_split(df,
prop = 0.8)
train<- df_split %>%
training()
test <- df_split %>%
testing()
# Número de datos en train
paste0("Total del dataset de entrenamiento: ", nrow(train))[1] "Total del dataset de entrenamiento: 800"Analizamos cómo quedo el balance de clases para cardio en cada dataset.
# calculamos la distribución de clase en cada dataset
train_ <- train %>%
group_by(cardio) %>%
summarise(numero_casos=n()) %>%
mutate(prop = round(prop.table(numero_casos)*100,2))
test_ <- test %>%
group_by(cardio) %>%
summarise(numero_casos=n()) %>%
mutate(prop = round(prop.table(numero_casos)*100,2))
# armamos tabla conjunta para graficar
distrib = cbind(rbind(train_, test_), dataset = c("train", "train", "test", "test"))
# graficamos las distribuciones
ggplot(distrib, aes(x = cardio, y = prop, fill = factor(cardio), label = prop)) +
geom_bar(stat="identity", position = "dodge") + facet_wrap(~ dataset) +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
labs(x = "", y = "Proporción en %", title = "Proporción de cardio por dataset") +
guides(fill=guide_legend(title="cardio"))+
theme_bw() +
scale_fill_brewer(palette="Set2")
glimpse(train)Rows: 800
Columns: 13
$ id <int> 56139, 22387, 68567, 99348, 82747, 34471, 6559, 95883, 513…
$ age <int> 23106, 14799, 20578, 18994, 14463, 20587, 14849, 15050, 16…
$ gender <int> 1, 2, 1, 1, 1, 2, 2, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ height <int> 152, 171, 168, 163, 157, 168, 172, 158, 164, 169, 160, 155…
$ weight <dbl> 74, 72, 62, 104, 52, 70, 75, 50, 64, 65, 76, 90, 56, 89, 5…
$ ap_hi <int> 160, 120, 120, 130, 110, 150, 120, 120, 100, 110, 130, 130…
$ ap_lo <int> 100, 90, 80, 79, 70, 100, 90, 80, 70, 80, 70, 100, 60, 90,…
$ cholesterol <int> 2, 1, 1, 1, 2, 1, 2, 1, 1, 2, 2, 2, 1, 1, 2, 1, 1, 1, 2, 1…
$ gluc <int> 2, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1, 1, 2, 1…
$ smoke <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ alco <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ active <int> 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1…
$ cardio <int> 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1…Seleccioamos las variables con las que vamos a trabajar y les asignamos el tipo de dato correspondiente.
cols <- c("age",
"height",
"weight",
"gender",
"cholesterol",
"gluc",
"cardio")
train$gender <- as.factor(train$gender)
train$cholesterol <- as.factor(train$cholesterol)
train$gluc <- as.factor(train$gluc)
test$gender <- as.factor(test$gender)
test$cholesterol <- as.factor(test$cholesterol)
test$gluc <- as.factor(test$gluc)# graficamos con ggpairs coloreando por variable a predecir
g <- train %>%
mutate(cardio = factor(cardio)) %>%
select(all_of(cols)) %>%
ggpairs(title = "Correlograma de variables",
mapping = aes(colour= cardio),
progress = FALSE,
lower=list(combo=wrap("facethist", binwidth=0.8))) +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
theme_bw() +
scale_fill_brewer(palette="Set2") +
scale_color_brewer(palette="Set2")
g
Cociente entre la probabilidad de éxito de un evento y la probabilidad de que no ocurra.
\(Odds = \frac{P(x)}{1-P(x)}= \frac{\text{Nro Eventos a favor}}{\text{Nro Eventos en contra}}\)
Ej. si Odds es 3 indica que por cada 4 repeticiones del evento se espera que ocurran 3 y que una no ocurra.
Odds < 1 : es mas probable que no ocurra a que si
Odds = 1 : equiprobabilidad
Odds > 1: es mas probable que el evento ocurra a que no.
En nuestro caso:
\(Odds = \frac{391/800}{409/800}= 0.9558 \approx \frac{191}{200}\)
Por cada 191 pacientes con ECV, hay 200 pacientes sin ECV.
\(Logit = \log {\frac{P(x)}{1-P(x)}}= \log Odds\)
Es una medida de magnitud de efecto. Si tenemos el siguiente modelo de regresión logística
\(\log\left(\frac{p}{1 - p}\right) = \beta_0 + \beta_1 X\)
El odds ratio asociado a la variable \(\ X\) es:
\(\text{OR} = e^{\beta_1}\)
OR = 1 No hay efecto
OR > 1 asociación positiva
OR < 1 asociación negativa
1- Predictor lineal
2- Odds
3- Variable respuesta
\(P(Y=1|X)= \frac{e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X_j}}{1+e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X_j}}\)
Esta función acota el resultado entre 0 y 1, lo cual es mucho más adecuado para modelar una probabilidad.
Luego de hacer algunas operaciones, podemos llegar a la expresión:
\(\log {\frac{P(x)}{1-P(x)}}= \beta_0 + \sum\limits_{j=1}^p \beta_j X_j\)
Realizamos un modelo de regresión logística para predecir si el paciente tiene una ECV en función de weight y gender.
\(\log {\frac{P(x)}{1-P(x)}}= \beta_0 + \beta_1 weight + \beta_2 varón\)
# modelo de regresión logística
glm1 <- glm(data = train,
cardio ~ weight + gender,
family = 'binomial')
# veo los resultados
tidy(glm1)# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -2.71 0.402 -6.75 1.45e-11
2 weight 0.0373 0.00547 6.82 8.87e-12
3 gender2 -0.203 0.159 -1.27 2.02e- 1glance(glm1)# A tibble: 1 × 8
null.deviance df.null logLik AIC BIC deviance df.residual nobs
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 1109. 799 -528. 1063. 1077. 1057. 797 800Modelo estimado
\(\log {\frac{P(x)}{1-P(x)}}= -2.71 + 0.03* weight - 0.2*varón\)
Peso
Se espera un aumento de 0.03 unidades en el log de odds por cada aumento de 1kgr en el peso (manteniendo el resto de las variables constantes). La probabilidad de contraer una ECV aumenta con el aumento de peso del paciente.
\(\ OR = e^{\beta_{weight}}= e^{0.03}= 1.03\)
Por cada aumento de 1 kgr se estima que el odds de contraer una ECV aumente 3.05%
Género
\(\ Odds_{varón} = {\frac{Prob_{ECV} /varón}{Prob_{noECV}/varón}}\)
\(\ Odds_{mujer} = {\frac{Prob_{ECV} /mujer}{Prob_{noECV}/mujer}}\)
\(\ OR = {\frac{Odds_{varón}}{Odds_{mujer}}}= e^{\beta_{género}}= e^{-0.203}\)
El odds ratio para la variable género es 0.816, lo que indica que los varones tienen un 18,37 % menos de odds de presentar una ECV que las mujeres, manteniendo constante el peso. En consecuencia, ser mujer se asocia con un riesgo mayor riesgo relativo de ECV.
Es una medida de la falta de ajuste del modelo. Cuanto menor es la devianza residual, mejor se ajusta el modelo.
\(\ Devianza Explicada = {\frac{Dev_{mod null} - Dev_{Mod}}{Dev_{mod null}}}\)
\(\ Devianza Explicada = {\frac{1108.71 - 1056.63}{1108.71}}= 0.046\)
1 - (glm1$deviance / glm1$null.deviance)[1] 0.04698285Se generan la combinación de todos los modelos posibles a partir de las variables gender, weight, y cholesterol. Antes hacemos un boxplot con weight, cholesterol y cardio.
ggplot(train, aes(x = cholesterol, y = weight, fill = factor(cardio))) +
geom_boxplot() +
labs(x = "cholesterol", y = "Peso", fill = "cardio") +
theme_bw() +
scale_fill_brewer(palette="Set2")
vars <- c("gender",
"weight",
"cholesterol")
# Crear todas las combinaciones posibles de predictores
combs <- unlist(
lapply(1:length(vars), function(i) combn(vars, i, simplify = FALSE)),
recursive = FALSE
)
# Crear una lista de fórmulas
formulas <- map(combs, ~ reformulate(termlabels = .x, response = "cardio"))
formulas[[1]]
cardio ~ gender
<environment: 0x0000022c3d65ae78>
[[2]]
cardio ~ weight
<environment: 0x0000022c324a2d68>
[[3]]
cardio ~ cholesterol
<environment: 0x0000022c2b91e1d0>
[[4]]
cardio ~ gender + weight
<environment: 0x0000022c3fc33628>
[[5]]
cardio ~ gender + cholesterol
<environment: 0x0000022c3fb14518>
[[6]]
cardio ~ weight + cholesterol
<environment: 0x0000022c3fb76318>
[[7]]
cardio ~ gender + weight + cholesterol
<environment: 0x0000022c3fa495c0>#se realizan las regresiones logisticas para todos los modelos
modelos <- map(formulas, ~ glm(.x, data = train, family = binomial))Analizamos la performance de los modelos generados
resultados <- tibble(
formula = map_chr(formulas, format),
modelo = modelos
) %>%
mutate(
deviance = map_dbl(modelo, ~ .x$deviance),
null_deviance = map_dbl(modelo, ~ .x$null.deviance),
AIC = map_dbl(modelo, AIC),
dev_explicada = 1 - deviance / null_deviance,
roc_obj = map(modelo, ~ roc(train$cardio,
predict(.x, type = "response"),
quiet = TRUE)),
AUC = map_dbl(roc_obj, ~ as.numeric(auc(.x)))
) %>%
arrange(desc(dev_explicada))
# Tabla para mostrar (sin las list-columns)
resultados %>%
select(formula, deviance, null_deviance, AIC, dev_explicada, AUC)# A tibble: 7 × 6
formula deviance null_deviance AIC dev_explicada AUC
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 cardio ~ gender + weight + c… 1037. 1109. 1047. 0.0644 0.667
2 cardio ~ weight + cholesterol 1038. 1109. 1046. 0.0637 0.665
3 cardio ~ gender + weight 1057. 1109. 1063. 0.0470 0.648
4 cardio ~ weight 1058. 1109. 1062. 0.0455 0.644
5 cardio ~ gender + cholesterol 1080. 1109. 1088. 0.0262 0.576
6 cardio ~ cholesterol 1080. 1109. 1086. 0.0262 0.581
7 cardio ~ gender 1109. 1109. 1113. 0.0000850 0.505roc_df <- resultados %>%
mutate(curva = map(roc_obj, ~ tibble(tpr = .x$sensitivities,
fpr = 1 - .x$specificities))) %>%
select(formula, curva) %>%
unnest(curva)
ggplot(roc_df, aes(x = fpr, y = tpr, color = formula)) +
geom_line(size = 1.2) +
geom_abline(linetype = "dashed") +
labs(x = "1 - Especificidad", y = "Sensibilidad", title = "Curvas ROC de todos los modelos evaluado en TRAIN") +
theme_minimal()
Hasta ahora hemos evaluado el modelo de manera general, pero el resultado final del modelo debe consistir en asignar a la persona una clase predicha. En nuestro caso debemos establecer un punto de corte según el cual vamos a separar a los pacientes con ECV de los que no tienen ECV.
Probamos varios puntos de corte y graficamos el accuracy, la sensibilidad o recall, la especificidad y la precisión para cada uno de ellos.
| Clases predichas / Clases | Negativa | Positiva |
|---|---|---|
| Negativa | True Neg | False Neg |
| Positiva | False Pos | True Pos |
Recordemos que:
\(accuracy = \frac{TP+TN}{TP+FP+FN+TN}\)
\(sensitivity = recall = \frac{TP}{TP+FN}\)
\(specificity = \frac{TN}{TN+FP}\)
\(precision = \frac{TP}{TP+FP}\)
# Supongamos que 'mejor_modelo'
best_model <- glm(data = train,
cardio ~ weight + gender+ cholesterol,
family = 'binomial')
# generamos las predicciones para el mejor modelo
train_model <- tibble(
.fitted = predict(best_model, newdata = train, type = "response"),
cardio = train$cardio
)
# generamos una funcion que evalua diferentes umbrales y calcula las metricas
prediction_metrics <- function(cutoff, predictions) {
tab <- predictions %>%
mutate(predicted_class = factor(if_else(.fitted > cutoff, 1, 0), levels = c(0,1)),
cardio = factor(cardio, levels = c(0,1)))
cm <- confusionMatrix(tab$predicted_class, tab$cardio, positive = "1")
tidy(cm) %>%
select(term, estimate) %>%
filter(term %in% c('accuracy','sensitivity','specificity','precision')) %>%
mutate(cutoff = cutoff)
}# se calculan las metricas para diferentes valores de umbrales
cutoffs <- seq(0.05, 0.95, 0.05)
metricas_df <- map_df(cutoffs, ~ prediction_metrics(.x, predictions = train_model)) %>%
mutate(term = as.factor(term), estimate = round(estimate, 3))
#graficamos los resultados
ggplot(metricas_df, aes(x = cutoff, y = estimate, color = term, group = term)) +
geom_line(size = 1.2) +
geom_point(size = 2) +
scale_x_continuous(breaks = seq(0, 1, by = 0.1))+
theme_bw() +
labs(
title = "Accuracy, Sensitivity, Specificity y Precision vs Cutoff",
x = "Punto de corte",
y = "Valor de la métrica",
color = "Métrica"
)
¿Cuál seria el umbral adecuado?
Punto de cruce entre sensitivity y specificity dado que se minimizan tanto FP y FN. Útil cuando no hay una clara preferencia por minimizar FP o FN.
sel_cutoff = 0.45
# calculamos las predicciones sobre el dataset de train
table_train = augment(x = best_model,
newdata = train,
type.predict='response')
# Clasificamos utilizamos el punto de corte
table_train = table_train %>%
mutate(predicted_class = if_else(.fitted>sel_cutoff, 1, 0) %>% as.factor(),
cardio = factor(cardio))
# Creamos la matriz de confusión
confusionMatrix(table(table_train$predicted_class, table_train$cardio), positive = "1")Confusion Matrix and Statistics
0 1
0 250 137
1 158 255
Accuracy : 0.6312
95% CI : (0.5968, 0.6648)
No Information Rate : 0.51
P-Value [Acc > NIR] : 3.229e-12
Kappa : 0.263
Mcnemar's Test P-Value : 0.2442
Sensitivity : 0.6505
Specificity : 0.6127
Pos Pred Value : 0.6174
Neg Pred Value : 0.6460
Prevalence : 0.4900
Detection Rate : 0.3187
Detection Prevalence : 0.5162
Balanced Accuracy : 0.6316
'Positive' Class : 1
# Agregamos la predicciones al dataset de testeo
table_test = augment(x = best_model,
newdata=test,
type.predict='response')
# Clasificamos utilizando el punto de corte
table_test = table_test %>%
mutate(predicted_class = if_else(.fitted>sel_cutoff, 1, 0))
# Creamos la matriz de confusión
confusionMatrix(table(table_test$predicted_class, table_test$cardio), positive = "1")Confusion Matrix and Statistics
0 1
0 68 28
1 36 68
Accuracy : 0.68
95% CI : (0.6105, 0.744)
No Information Rate : 0.52
P-Value [Acc > NIR] : 3.184e-06
Kappa : 0.361
Mcnemar's Test P-Value : 0.3816
Sensitivity : 0.7083
Specificity : 0.6538
Pos Pred Value : 0.6538
Neg Pred Value : 0.7083
Prevalence : 0.4800
Detection Rate : 0.3400
Detection Prevalence : 0.5200
Balanced Accuracy : 0.6811
'Positive' Class : 1
Sensibilidad: El 70.8 % de los casos con ECV fueron correctamente identificados por el modelo. (Probabilidad de que el modelo prediga ECV si la persona realmente la tiene.)
Especificidad: El 65.4 % de los casos sin ECV fueron correctamente identificados como tales. (Probabilidad de que el modelo prediga ausencia de ECV si la persona realmente no la tiene.)
De todos los pacientes que el modelo predijo como negativos, el 70.8 % realmente no tiene ECV.
Los gráficos de residuos vs predichos no son informativos para evaluar el ajuste del modelo, ya que la variable solo toma dos valores.
En el eje X se observa la probabilidad predicha de contrarer una ECV.
En el eje Y se observa la frecuencia de clase, el cociente entre cantidad de individuos con ECV y el total de pacientes.
La línea punteada designa la igualdad entre probabilidad predicha y frecuencia de clase.
Los círculos, que se construyen de la siguiente manera:
Aquellos círculos que se ubiquen por encima de la línea punteada indican que el modelo está subestimando la probabilidad para dichos grupos. Mientras que si los círculos se ubican por debajo el modelo está sobreestimando la probabilidad para dichos grupos.
¿Para qué valores parece existir una sobreestimación de la probabilidad? ¿Para cuáles subestimación?
Hosmer_Lemeshow_plot <- function(dataset, predicted_column, class_column, bins = 10,
positive_value = 1, color = 'forestgreen', nudge_x = 0, nudge_y = 0.05){
dataset <- dataset %>%
mutate(group = cut(!!sym(predicted_column),
breaks = quantile(!!sym(predicted_column), probs = seq(0, 1, length.out = bins + 1)),
include.lowest = TRUE, labels = 1:bins))
positive_class <- dataset %>%
filter(!!sym(class_column) == positive_value) %>%
group_by(group) %>%
summarise(n_positive = n())
HL_df <- dataset %>%
group_by(group) %>%
summarise(pred = mean(!!sym(predicted_column)), n_total = n()) %>%
left_join(positive_class, by = "group") %>%
mutate(n_positive = replace_na(n_positive, 0),
freq = n_positive / n_total)
ggplot(HL_df, aes(x = pred, y = freq)) +
geom_point(aes(size = n_total), color = color) +
geom_text(aes(label = n_total), nudge_x = nudge_x, nudge_y = nudge_y, size = 3) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
labs(title = "Hosmer-Lemeshow",
subtitle = "Modelo final",
x = "Probabilidad Predicha",
y = "Frecuencia observada",
size = "Casos") +
theme_minimal()
}# modelo completo
Hosmer_Lemeshow_plot(train_model, '.fitted', 'cardio', 10, 1) +
labs(subtitle="Modelo final")
La prueba de Hosmer Lemeshow tiene la siguiente Ho
Ho buen ajuste del modelo
library(ResourceSelection)
p_test <- predict(best_model, newdata = test, type = "response")
h1 <- hoslem.test(test$cardio, p_test, g = 10)
h1
Hosmer and Lemeshow goodness of fit (GOF) test
data: test$cardio, p_test
X-squared = 7.3766, df = 8, p-value = 0.4966results_ <- tidy(best_model, exponentiate = TRUE, conf.int = TRUE)
results_ %>%
filter(term != "(Intercept)") %>% # elimina el intercepto
mutate(across(c(estimate, conf.low, conf.high), round, 2)) %>%
mutate(IC95 = paste0("[", conf.low, ", ", conf.high, "]")) %>%
select(term, OR = estimate, IC95, p.value)# A tibble: 4 × 4
term OR IC95 p.value
<chr> <dbl> <chr> <dbl>
1 weight 1.04 [1.02, 1.05] 4.96e-10
2 gender2 0.87 [0.63, 1.19] 3.81e- 1
3 cholesterol2 1.46 [0.95, 2.25] 8.83e- 2
4 cholesterol3 2.71 [1.7, 4.41] 3.68e- 5weight
Por cada aumento de una unidad de peso, los odds de presentar una ECV aumentan un 4%, manteniendo constantes el resto de las variables. Al no incluir el 1 se puede decir que el efecto es estadísticamente significativo
Conclusión: el peso se asocia positivamente con la probabilidad de contraer ECV.
gender
Aunque parece haber una tendencia a menor riesgo en varones, el intervalo incluye al 1, por lo tanto no hay evidencia estadística suficiente para afirmarlo.
cholesterol
Los pacientes con colesterol en nivel 3 tienen 2.71 veces más odds de presentar ECV que aquellos con colesterol normal (nivel base). Es decir, los odds de tener ECV se incrementan en un 171% en comparación con el grupo de referencia, manteniendo constantes las demás variables del modelo.
A continuación se muestra las curvas que representan la probabilidad predicha de presentar una ECV en función del peso, comparando varones y mujeres y manteniendo fijo el colesterol en el nivel “1”.
# Crear nueva grilla de datos
new_data <- expand.grid(
weight = seq(min(train$weight), max(train$weight), length.out = 150),
cholesterol = "1", # fijamos un valor de colesterol
gender = levels(train$gender)
)
# Obtener predicciones en escala del logit con error estándar
pred <- predict(best_model, newdata = new_data, type = "link", se.fit = TRUE)
# Calcular límites de confianza en la escala del logit
new_data <- new_data %>%
mutate(
fit = pred$fit,
se.fit = pred$se.fit,
lower = fit - 1.96 * se.fit,
upper = fit + 1.96 * se.fit,
predicted = plogis(fit), # transformar a probabilidad
lower_ci = plogis(lower),
upper_ci = plogis(upper)
)
# Graficar con bandas de confianza
ggplot(new_data, aes(x = weight, y = predicted, color = gender, fill = gender)) +
geom_line(size = 1.2) +
geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), alpha = 0.2, color = NA) +
labs(
x = "Peso",
y = "Probabilidad predicha de ECV",
color = "Género",
fill = "Género"
) +
theme_minimal(base_size = 14)