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)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.
El objetivo de esta clase es 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
En este caso estamos modelando la probabilidad de la siguiente manera:
\(P(X)= \beta_0 + \sum\limits_{j=1}^p \beta_j X_j\)
Veamos que tan bueno es el modelo lineal para esto, usando el peso como predictor.
mrl <- train %>%
lm(formula = cardio ~ weight)
tdy = mrl %>% tidy()
tdy# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -0.131 0.0880 -1.49 1.38e- 1
2 weight 0.00845 0.00117 7.19 1.48e-12mrl %>% glance()# A tibble: 1 × 12
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.0609 0.0597 0.485 51.7 1.48e-12 1 -555. 1117. 1131.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>Los estimadores son significativos y el test de significatividad global del modelo también es significativo.
Veamos un gráfico de nuestro modelo.
Parece tener bastantes problemas para estimar la probabilidad de supervivencia de los individuo: no existe un punto de corte claro, la predicción podría ser mayor a 1 o menor a cero llegado el caso.
Cociente entre la probabilidad de éxito de un evento y la probabilidad de que no ocurra.
\(Odds = \frac{P(x)}{1-P(x)}\)
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, 200 no la tienen.
\(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
Para evitar los problemas vistos con regresión lineal, usamos la función logística.
\(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\)
La funcíón glm() nos permite crear un modelo lineal generalizado (Generalized Linear Model). Al igual que la función lm() toma como argumentos una formula y los datos pero también se debe especificar el argumento family: indicamos la distribución del error y la función link que vamos a utilizar en el modelo.
Algunas familias son:
Binomial: link=logit
Poisson: link=log
Gaussiana: link=identidad
Como estamos trabajando con un fenómeno que suponemos tiene una distribución binomial, así lo especificamos en el parámetro family.
Realizamos un modelo de regresión logística para predecir la 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 age + \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.02}\)
El odds ratio para la variable género es 0.98, lo que indica que los varones tienen un 2 % 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: 0x000002490961be78>
[[2]]
cardio ~ weight
<environment: 0x000002490960f5f0>
[[3]]
cardio ~ cholesterol
<environment: 0x0000024909015008>
[[4]]
cardio ~ gender + weight
<environment: 0x0000024908dbe730>
[[5]]
cardio ~ gender + cholesterol
<environment: 0x0000024908d8f568>
[[6]]
cardio ~ weight + cholesterol
<environment: 0x0000024908d695e0>
[[7]]
cardio ~ gender + weight + cholesterol
<environment: 0x0000024908d50760>#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),
deviance = map_dbl(modelos, ~ .$deviance),
null_deviance = map_dbl(modelos, ~ .$null.deviance),
AIC = map_dbl(modelos, AIC)
) %>%
mutate(dev_explicada = 1 - deviance / null_deviance)%>%
arrange(desc(dev_explicada))# se ordena segñun devianza explicada
resultados# A tibble: 7 × 5
formula deviance null_deviance AIC dev_explicada
<chr> <dbl> <dbl> <dbl> <dbl>
1 cardio ~ gender + weight + cholest… 1037. 1109. 1047. 0.0644
2 cardio ~ weight + cholesterol 1038. 1109. 1046. 0.0637
3 cardio ~ gender + weight 1057. 1109. 1063. 0.0470
4 cardio ~ weight 1058. 1109. 1062. 0.0455
5 cardio ~ gender + cholesterol 1080. 1109. 1088. 0.0262
6 cardio ~ cholesterol 1080. 1109. 1086. 0.0262
7 cardio ~ gender 1109. 1109. 1113. 0.0000850prediction_full <- map2(modelos, formulas, ~ {
tibble(
.fitted = predict(.x, newdata = train, type = "response"),
cardio = train$cardio
)
})roc_df <- map2_dfr(modelos, resultados$formula, ~ {
roc_obj <- roc(train$cardio, predict(.x, type = "response"))
tibble(
tpr = roc_obj$sensitivities,
fpr = 1 - roc_obj$specificities,
formula = .y
)
})
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") +
theme_minimal()
auc_values <- map_dbl(modelos, ~ auc(roc(train$cardio, predict(.x, type = "response"))))
resultados <- resultados %>%
mutate(AUC = auc_values)
resultados# 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.505
2 cardio ~ weight + cholesterol 1038. 1109. 1046. 0.0637 0.644
3 cardio ~ gender + weight 1057. 1109. 1063. 0.0470 0.581
4 cardio ~ weight 1058. 1109. 1062. 0.0455 0.648
5 cardio ~ gender + cholesterol 1080. 1109. 1088. 0.0262 0.576
6 cardio ~ cholesterol 1080. 1109. 1086. 0.0262 0.665
7 cardio ~ gender 1109. 1109. 1113. 0.0000850 0.667Hasta 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 = test,
cardio ~ weight + gender+ cholesterol,
family = 'binomial')
# generamos las predicciones para el mejor modelo
test_model <- tibble(
.fitted = predict(best_model, newdata = test, type = "response"),
cardio = test$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 = test_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 261 151
1 147 241
Accuracy : 0.6275
95% CI : (0.593, 0.6611)
No Information Rate : 0.51
P-Value [Acc > NIR] : 1.434e-11
Kappa : 0.2546
Mcnemar's Test P-Value : 0.862
Sensitivity : 0.6148
Specificity : 0.6397
Pos Pred Value : 0.6211
Neg Pred Value : 0.6335
Prevalence : 0.4900
Detection Rate : 0.3013
Detection Prevalence : 0.4850
Balanced Accuracy : 0.6273
'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 69 29
1 35 67
Accuracy : 0.68
95% CI : (0.6105, 0.744)
No Information Rate : 0.52
P-Value [Acc > NIR] : 3.184e-06
Kappa : 0.3605
Mcnemar's Test P-Value : 0.532
Sensitivity : 0.6979
Specificity : 0.6635
Pos Pred Value : 0.6569
Neg Pred Value : 0.7041
Prevalence : 0.4800
Detection Rate : 0.3350
Detection Prevalence : 0.5100
Balanced Accuracy : 0.6807
'Positive' Class : 1
El 69.8 % de los casos positivos reales fueron correctamente identificados. (Probabilidad de detectar correctamente una ECV si la persona realmente la tiene.)
El 66.4 % de los casos negativos fueron correctamente identificados como tales.
De todos los que el modelo predijo como negativos, el 70.4 % realmente no tienen ECV.
best_model
Call: glm(formula = cardio ~ weight + gender + cholesterol, family = "binomial",
data = test)
Coefficients:
(Intercept) weight gender2 cholesterol2 cholesterol3
-3.38811 0.04432 -0.51426 0.53028 1.27021
Degrees of Freedom: 199 Total (i.e. Null); 195 Residual
Null Deviance: 276.9
Residual Deviance: 249.3 AIC: 259.3Los 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(test_model, '.fitted', 'cardio', 10, 1) +
labs(subtitle="Modelo final")
La prueba de Hosmer Lemeshow permite evaluar el buen ajuste del modelo
library(ResourceSelection)
h1 <- hoslem.test (test$cardio ,
fitted(best_model), g=8)
h1
Hosmer and Lemeshow goodness of fit (GOF) test
data: test$cardio, fitted(best_model)
X-squared = 2.0033, df = 6, p-value = 0.9194results_ <- 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.05 [1.02, 1.07] 0.000249
2 gender2 0.6 [0.31, 1.13] 0.118
3 cholesterol2 1.7 [0.69, 4.32] 0.254
4 cholesterol3 3.56 [1.35, 10.57] 0.0139 weight
Por cada aumento de una unidad de peso, los odds de presentar una ECV aumentan un 5%, 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 3.56 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 256% 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)