1. Modelo explicativo

La pregunta es ¿Cuál es el efecto de X sobre Y? No me interesa tanto acertar qué paciente tendrá el evento. Lo que quiero saber es si una variable produce, aumenta o disminuye el riesgo del desenlace. Elección de variables: Pongo la variable explicativa y los confundidores Las variables entran se seleccionan porque pueden modificar el estimador del efecto de la variable principal. En un modelo explicativo el coeficiente es el protagonista. Selección de variables: cuáles son las variables que confunden la relación X-Y? Utilizar DAGs, conocimiento clínico, confundidores conocidos.En un modelo explicativo. Solo uno (o muy pocos) coeficientes tienen interpretación causal. Los demás son variables de ajuste.

  1. Modelo predictivo: La pregunta cambia completamente. Ahora es

¿Qué combinación de variables me permite predecir mejor quién tendrá el evento? Ya no me interesa el efecto aislado de ninguna variable. Me interesa que el modelo haga la mejor predicción posible. Lo único que quiero es minimizar el error de predicción. Entonces las variables se eligen porque agregan información predictiva. No porque sean confundidores. Mejoran el AUC. En un modelo predictivo los OR pasan a segundo plano.

Lo importante es AUC calibración sensibilidad especificidad Brier Score validación

Selección de variables: Entonces pueden utilizarse

selección hacia atrás (backward elimination) penalización (como LASSO) árboles random forest gradient boosting

La selección depende del rendimiento predictivo, no de la significación estadística.

El cribado univariado (“univariate screening”) es una estrategia de selección de variables orientada principalmente a modelos predictivos, y no es la estrategia recomendada para modelos explicativos.Por qué se hace? Porque reduce el número de variables candidatas.

Supongamos que tenés 40 predictores. No querés empezar el modelo con los 40. Entonces primero hacés un filtro.No debería utilizarse para decidir qué variables ajustar en un modelo explicativo. En un modelo explicativo, las variables se eligen por su papel causal o como potenciales confundidores, no por su p-valor.

  1. Modelo asociativo o explicativo exploratorio.(qué variables se asocian con la mortalidad) Supongamos que recolectaron:

Edad, Sexo, Diabetes, Hipertensión, Insuficiencia renal, Albúmina, Lactato, APACHE II Ventilación mecánica, etc.

Primero hacen análisis univariado. Luego dicen:

“Variables con p < 0.10 fueron incluidas en el análisis multivariable.” y concluyen:

“Factores independientemente asociados con la mortalidad…” La pregunta que intenta responder es:

¿Qué variables se asocian independientemente con la mortalidad? No están intentando demostrar causalidad. No están intentando construir un score de riesgo.

Simplemente quieren identificar asociaciones ajustadas. Cuál es el problema? 1. Puede eliminar confundidores importantes, Puede eliminar variables que cobran importancia al ajustarse, El p del univariado depende del tamaño muestral. Un análisis univariado seguido de un modelo multivariable puede ser aceptable como una estrategia exploratoria, siempre que el autor sea consciente de sus limitaciones y no extraiga conclusiones causales. El error aparece cuando el trabajo concluye: “La hipoalbuminemia es un factor independiente de mortalidad.” En realidad, lo que el modelo mostró es: “En esta cohorte, la hipoalbuminemia se asoció con la mortalidad tras ajustar por las variables incluidas en el modelo.” Si el objetivo desde el inicio fue: “Identificar variables asociadas con la mortalidad.” Entonces sí es esperable interpretar varios coeficientes. Pero aun así hay que ser prudente. La conclusión debería ser “Variables independientemente asociadas…” No “Factores de riesgo.” Ni mucho menos “Predictores causales.”

Porque sigue siendo un estudio observacional.

Hay algo que suele confundir mucho, incluso a investigadores con experiencia. La falacia de la Tabla 2 no dice que esté prohibido interpretar más de un coeficiente. Lo que dice es que la interpretación debe ser coherente con la pregunta de investigación y con el papel que cada variable desempeña en el modelo.

Si el objetivo era estimar el efecto de la isquemia, los demás coeficientes son principalmente parámetros de ajuste y no deberían convertirse en resultados principales. Si el objetivo era identificar asociaciones ajustadas de varias variables con la mortalidad, entonces es razonable presentar e interpretar varios coeficientes, siempre dejando claro que se trata de asociaciones ajustadas y no de efectos causales.

FALACIA DE LA TABLA 2:

El artículo dice “Los factores independientes asociados con la mortalidad fueron…” El investigador armó el modelo porque quería estudiar, por ejemplo, el efecto de la diabetes sobre la mortalidad. Para controlar confusión incluyó edad, sexo, tabaquismo, hipertensión Hasta ahí, perfecto. Pero después hace algo que metodológicamente no corresponde.

Empieza a interpretar TODOS los coeficientes del modelo como si fueran el objetivo del estudio. Y escribe “La edad fue un factor independiente.”¿Por qué es un error? Porque esas variables no fueron incluidas para responder una pregunta sobre ellas. Fueron incluidas para ajustar la estimación del efecto de la diabetes.

##Base simulada

# =========================================================
# Base simulada para aprender regresión logística
# Outcome: muerte hospitalaria sí/no
# n = 500 pacientes
# Mortalidad aproximada: 20%
# =========================================================

set.seed(123)

n <- 500

base_logistica <- data.frame(
  id = 1:n,
  edad = round(rnorm(n, mean = 65, sd = 12)),
  sexo = factor(rbinom(n, 1, 0.55),
                levels = c(0, 1),
                labels = c("Femenino", "Masculino")),
  diabetes = factor(rbinom(n, 1, 0.30),
                    levels = c(0, 1),
                    labels = c("No", "Si")),
  hta = factor(rbinom(n, 1, 0.55),
               levels = c(0, 1),
               labels = c("No", "Si")),
  fumador = factor(rbinom(n, 1, 0.25),
                   levels = c(0, 1),
                   labels = c("No", "Si")),
  colesterol = round(rnorm(n, mean = 210, sd = 40)),
  creatinina = round(rlnorm(n, meanlog = log(1.1), sdlog = 0.35), 2)
)

# Limitar valores poco realistas
base_logistica$edad <- pmin(pmax(base_logistica$edad, 18), 95)
base_logistica$colesterol <- pmin(pmax(base_logistica$colesterol, 100), 380)

# Crear variables numéricas auxiliares para simular el evento
sexo_num <- ifelse(base_logistica$sexo == "Masculino", 1, 0)
dm_num <- ifelse(base_logistica$diabetes == "Si", 1, 0)
hta_num <- ifelse(base_logistica$hta == "Si", 1, 0)
fumador_num <- ifelse(base_logistica$fumador == "Si", 1, 0)

# Modelo "verdadero" que genera la probabilidad de muerte
logit_p <- -7 +
  0.055 * base_logistica$edad +
  0.55  * sexo_num +
  0.80  * dm_num +
  0.25  * hta_num +
  0.65  * fumador_num +
  0.006 * base_logistica$colesterol +
  0.70  * base_logistica$creatinina

p_muerte <- exp(logit_p) / (1 + exp(logit_p))

# Simular muerte según la probabilidad individual
base_logistica$muerte <- rbinom(n, 1, p_muerte)

base_logistica$muerte <- factor(base_logistica$muerte,
                                levels = c(0, 1),
                                labels = c("No", "Si"))

# Ver mortalidad observada
prop.table(table(base_logistica$muerte))
## 
##    No    Si 
## 0.602 0.398
table(base_logistica$muerte)
## 
##  No  Si 
## 301 199

#Etapas de la confección del modelo

##Parte I: Construcción de modelos explicativos.

Etapa 1. Antes de ajustar el modelo

Definir el objetivo (explicativo, exploratorio o predictivo). Seleccionar las variables según ese objetivo. Verificar el número de eventos por predictor. Explorar las variables continuas y evaluar la linealidad del logit (si vas a modelarlas como lineales).

Etapa 2. Ajustar el modelo

Estimar la regresión logística. Interpretar los coeficientes (OR, IC95%, p).

Etapa 3. Diagnóstico del modelo

Multicolinealidad. Observaciones influyentes. Bondad de ajuste. Discriminación y calibración (según el objetivo). Validación (especialmente en modelos predictivos).

##Selección de variables Modelo explicativo: El objetivo es estimar el efecto del tabaquismo sobre la mortalidad. Por lo tanto, la variable de exposición (tabaquismo) se incorpora obligatoriamente al modelo. Además, se incluyen variables potencialmente confundidoras, seleccionadas sobre la base del conocimiento clínico y la evidencia científica, independientemente de su significación estadística. En este ejemplo, la edad se incorpora porque es un potencial confundidor: está asociada tanto con el tabaquismo como con la mortalidad y podría distorsionar la estimación del efecto del tabaquismo si no se ajusta por ella.Las variables no se incorporan porque sean estadísticamente significativas, sino porque son necesarias para estimar de manera no sesgada el efecto de la exposición sobre el desenlace.Ojo! Necesito 10 eventos por variable continuar y 10 eventos por cantidad de categorías. Si pruebo edad y fumador (si/no) voy a necesitar al menos 30 eventos Antes de confeccionar el modelo: 2.Evaluar la Linealidad del logit a. Box-Tidwell. b. Restricted Cubic Splines. c. LOWESS sobre el logit. Si no es lineal, categorizar o splines, funciones binomiales, transformación

#Modelo explcativo
# Evaluar supuesto de linealidad entre variable continua y logit del evento
#Ajustar modelo de splines

#---------------
#Gráfico: RAZONABLEMENTE LINEAL
#---------------
library(rms)
library(splines)
dd <- datadist(base_logistica)
options(datadist = "dd")

fit <- lrm(muerte ~ rcs(edad,4) + fumador,
           data = base_logistica)

plot(Predict(fit, edad))

#-------------------------
#Comparación de modelos con y sin splines: NO HAY DIFERENCIA P 0.38
#-------------------------
modelo_lineal <- glm(
  muerte ~ edad + fumador,
  family = binomial,
  data = base_logistica
)
modelo_spline <- glm(
  muerte ~ ns(edad, df = 4) + fumador,
  family = binomial,
  data = base_logistica
)
anova(modelo_lineal, modelo_spline, test = "LRT")
## Analysis of Deviance Table
## 
## Model 1: muerte ~ edad + fumador
## Model 2: muerte ~ ns(edad, df = 4) + fumador
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       497        629                     
## 2       494        626  3     3.04     0.38
#------------------------
#BoxTidwell
#------------------------
# Box-Tidwell para edad : P>0.05
base_logistica$edad_logedad <- base_logistica$edad * log(base_logistica$edad)

modelo_bt <- glm(
  muerte ~ edad + edad_logedad + fumador,
  data = base_logistica,
  family = binomial
)

summary(modelo_bt)
## 
## Call:
## glm(formula = muerte ~ edad + edad_logedad + fumador, family = binomial, 
##     data = base_logistica)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -11.9528     5.4457   -2.19   0.0282 * 
## edad           0.6901     0.4252    1.62   0.1046   
## edad_logedad  -0.1232     0.0817   -1.51   0.1314   
## fumadorSi      0.6168     0.2161    2.85   0.0043 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 672.19  on 499  degrees of freedom
## Residual deviance: 626.34  on 496  degrees of freedom
## AIC: 634.3
## 
## Number of Fisher Scoring iterations: 4
#------------------------------
#otra evaluación gráfica: RAZONABLEMENTE LINEAL
#-----------------------------
base_logistica <- base_logistica %>%
  mutate(
    muerte_num = ifelse(muerte == "Si", 1, 0)
  )
base_logistica<- base_logistica %>% mutate(edad_q = ntile(edad, 5))


base_logistica %>% dplyr::group_by(edad_q) %>% dplyr::summarise(tar=mean(muerte_num)) %>% ggplot(aes(x=edad_q, y=tar))+geom_point()+geom_smooth(method = "lm", col = "blue")
## `geom_smooth()` using formula = 'y ~ x'

Si la curva es aproximadamente recta la edad puede permanecer lineal. Si hace una S, una U, o cambia mucho de pendiente, la linealidad no se cumple. Si p < 0,05 en la comparación de modelos con y sin solines→ el spline mejora significativamente el ajuste. Hay evidencia de que la relación entre edad y el logit no es lineal. En el modelo BoxTidwell Si la p de edad_logedad es <0.05 hay evidencia de no linealidad del logit para edad

#Modelo
modelo_ex <- glm(
  muerte ~ edad+ fumador,
  data = base_logistica,
  family = binomial
)

summary(modelo_ex)
## 
## Call:
## glm(formula = muerte ~ edad + fumador, family = binomial, data = base_logistica)
## 
## Coefficients:
##             Estimate Std. Error z value      Pr(>|z|)    
## (Intercept) -3.86459    0.59953   -6.45 0.00000000011 ***
## edad         0.04979    0.00879    5.66 0.00000001473 ***
## fumadorSi    0.61987    0.21578    2.87        0.0041 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 672.19  on 499  degrees of freedom
## Residual deviance: 628.77  on 497  degrees of freedom
## AIC: 634.8
## 
## Number of Fisher Scoring iterations: 4
tab_model(modelo_ex)
  muerte
Predictors Odds Ratios CI p
(Intercept) 0.02 0.01 – 0.07 <0.001
edad 1.05 1.03 – 1.07 <0.001
fumador [Si] 1.86 1.22 – 2.84 0.004
Observations 500
R2 Tjur 0.083

Por cada aumento de un año de edad, los odds de morir aumentan aproximadamente un 5%, manteniendo constante el tabaquismo. Los fumadores presentan 1.86 veces los odds de morir respecto de los no fumadores, ajustando por la edad. O los fumadores tienen un 86% mayor odds de muerte que los no fumadores. R2=0.083 significa que, en promedio,las probabilidades predichas de los fallecidos son solamente 8.3 puntos porcentuales mayores que las de los sobrevivientes.Eso indica una capacidad discriminativa modesta. AIC: Solo sirve para comparar modelos.A menor valor, mejor. Mientras más variables tenga,mayor es el AIC por penalización.

Qué hacer? 1. Evaluar multicolinealidad ante variables como PAS, PAD o Col total, LDL a través del VIF (VIF <5 tranquiliza bastante. Aunque algunos usan 10.) la evaluación de la multicolinealidad es recomendable tanto en modelos explicativos como predictivos. Si están muy correlacionados, el modelo tiene dificultades para distinguir cuánto del efecto corresponde a cada uno. La consecuencia es: errores estándar más grandes, intervalos de confianza más amplios, p-valores más altos, coeficientes más inestables.

#Evaluación de multicolinealidad. VIF debe ser menor a 5
#ojo. Se busca entre los predictores que finalmente quedaron en el modelo.
vif(modelo_ex)
##      edad fumadorSi 
##         1         1
# Selección de variables. Solo numéricas
vars <- base_logistica[, c("edad", "colesterol", "creatinina")]
# Matriz de dispersión
cor(vars, method = "spearman")
##               edad colesterol creatinina
## edad        1.0000    -0.0593    -0.0192
## colesterol -0.0593     1.0000     0.0023
## creatinina -0.0192     0.0023     1.0000
pairs(vars, main = "Matriz de graficos de dispersion")

library(GGally)
## Warning: package 'GGally' was built under R version 4.4.3
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
# Selección de variables
vars <- base_logistica[, c("edad", "colesterol","fumador")]
# Matriz de dispersión con GGally
ggpairs(vars) 
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Qué hacer con las variables si hay multicolinealidad Modelo explicativo: ¿Las dos son confundidores importantes? Si la respuesta es sí, quizá prefieras mantenerlas y aceptar que los errores estándar sean más grandes. El objetivo principal es controlar la confusión, no obtener el menor VIF posible. Solo eliminaría una si realmente ambas representan casi la misma información y existe una justificación clínica para conservar una sola. Modelo predictivo: Si dos variables aportan prácticamente la misma información, no necesito las dos. Podría: eliminar una,usar LASSO,usar Ridge (que justamente nació para manejar colinealidad).

3.Evaluar observaciones influyentes Cook’s distance leverage DFBETAs Qué hago si hay observ influyentes? En un modelo explicativo: revisar que no haya errores, realizar un análisis de sensibilidad. En un modelo predictivo ¿El modelo va a predecir correctamente en pacientes futuros? Si un único paciente cambia muchísimo los coeficientes, eso puede indicar sobreajuste (overfitting) o un modelo poco estable.LAs obs inf se analizan más desde la estabilidad del modelo, evaluando el auc con y sin la observación

#observaciones influyentes
#distancia de cook
cooks <- cooks.distance(modelo_ex)

plot(cooks,
     type = "h",
     main = "Cook's distance",
     ylab = "Cook's distance")

abline(h = 4 / nrow(base_logistica),
       col = "red",
       lty = 2)

which(cooks > 4 / nrow(base_logistica))
##  44  98 113 143 198 201 246 307 318 360 386 471 473 
##  44  98 113 143 198 201 246 307 318 360 386 471 473
#leverage: El leverage alto indica pacientes con combinaciones raras de predictores.
lev <- hatvalues(modelo_ex)

plot(lev,
     type = "h",
     main = "Leverage",
     ylab = "Hat values")

abline(h = 2 * length(coef(modelo_ex)) / nrow(base_logistica),
       col = "red",
       lty = 2)

which(lev > 2 * length(coef(modelo_ex)) / nrow(base_logistica))
##  26  44  70  97 125 137 145 147 149 164 174 196 201 209 216 246 265 295 310 343 
##  26  44  70  97 125 137 145 147 149 164 174 196 201 209 216 246 265 295 310 343 
## 360 371 380 392 427 456 471 477 
## 360 371 380 392 427 456 471 477
#DFbetas
#Los DFBETAs muestran cuánto cambia cada coeficiente si se elimina una observación.
dfb <- dfbetas(modelo_ex)

head(dfb)
##   (Intercept)       edad fumadorSi
## 1     0.03726 -0.0379275    0.0898
## 2     0.04039 -0.0277630   -0.0383
## 3     0.07334 -0.0869811    0.0277
## 4     0.01809 -0.0054446   -0.0358
## 5    -0.00866  0.0000413    0.0240
## 6    -0.06326  0.0737381   -0.0205
plot(dfb[, "edad"],
     type = "h",
     main = "DFBETAs para edad",
     ylab = "DFBETA")

abline(h = c(-2 / sqrt(nrow(base_logistica)),
              2 / sqrt(nrow(base_logistica))),
       col = "red",
       lty = 2)

plot(dfb[, "fumadorSi"],
     type = "h",
     main = "DFBETAs para fumador",
     ylab = "DFBETA")

abline(h = c(-2 / sqrt(nrow(base_logistica)),
              2 / sqrt(nrow(base_logistica))),
       col = "red",
       lty = 2)

2 / sqrt(500)
## [1] 0.0894
#observaciones por arriba deese valor pueden ser influyentes

4.Bondad del ajuste ¿El modelo representa adecuadamente los datos observados? Las probabilidades que predice el modelo se parecen a las que realmente ocurrieron? - Hosmer y Lemershow: Divide a los pacientes según su probabilidad predicha.Generalmente en 10 grupos (deciles). Luego compara Esperado vs Observado mediante un Chi-cuadrado.H₀=El modelo ajusta correctamente. No existen diferencias entre lo esperado y lo observado. Ojo.Depende mucho del tamaño de la muestra. Con muestras pequeñas puede dar 0.7 y con muestra muy grandes <0.05 - Gráfico de calibración (mejor en modelos predictivos)Es preferible. muestra la magnitud del error, muestra dónde ocurre, permite detectar sobreestimación o subestimación, no depende tanto de un único p-valor. Ojo. En un modelo explicativo no es obligatorio hacerlo porque el objetivo del modelo no es predecir probabilidades. Es estimar OR ajustados. En un modelo explicativo, la calibración es una propiedad interesante pero secundaria. En un modelo predictivo, la calibración es una propiedad esencial, porque el modelo se utilizará para estimar riesgos individuales.

#bondad del ajuste
#HyL
library(ResourceSelection)
## Warning: package 'ResourceSelection' was built under R version 4.4.3
## ResourceSelection 0.3-6   2023-06-27
## 
## Adjuntando el paquete: 'ResourceSelection'
## The following object is masked from 'package:prodlim':
## 
##     sindex
hoslem.test(
  x = base_logistica$muerte_num,
  y = fitted(modelo_ex),
  g = 10
)
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  base_logistica$muerte_num, fitted(modelo_ex)
## X-squared = 11, df = 8, p-value = 0.2
# gráfico de calibración de Harrell
library(rms)

modelo <- lrm(
  muerte ~ edad + fumador,
  data = base_logistica,
  x = TRUE,
  y = TRUE
)

cal <- calibrate(
  modelo,
  method = "boot",
  B = 200
)

plot(cal)

## 
## n=500   Mean absolute error=0.033   Mean squared error=0.00146
## 0.9 Quantile of absolute error=0.058

La línea continua es la corregida por bootstrap En valores bajos, La curva está ligeramente por debajo de la línea ideal. Eso significa que el modelo sobreestima ligeramente el riesgo. Error 0.033 Significa que en promedio las probabilidades predichas difieren 3.3 % de las observadas.(más del 10% es malo el error)

  1. Discriminación ROC. AUC (es más para predictivos pero puede informarse en explicativos) Un modelo puede tener una excelente discriminación (AUC alta) y, sin embargo, estar mal calibrado si sistemáticamente sobreestima o subestima el riesgo. Por eso, en modelos predictivos modernos, la calibración y la discriminación son propiedades complementarias y ambas deben evaluarse. ¿El modelo es capaz de distinguir correctamente a los pacientes que presentarán el evento de aquellos que no? Es decir, si separa correctamente los pacientes que van a tener el evento de los que no AUC: AUC = 0.5 El modelo discrimina igual que tirar una moneda. AUC=1 Discriminación perfecta. El AUC de 0.80 significa que en el 80% de las veces el modelo asignará un riesgo mayor al paciente que realmente tuvo el evento.Un AUC >0.7 es aceptable
#Discriminación
library(pROC)

roc1 <- roc(base_logistica$muerte, fitted(modelo_ex))
## Setting levels: control = No, case = Si
## Setting direction: controls < cases
plot(roc1,
     col = "blue",
     lwd = 2)

auc(roc1)
## Area under the curve: 0.664
ci.auc(roc1)
## 95% CI: 0.616-0.711 (DeLong)
  1. Interacciones
#Evaluación de interacciones. Las mismas debe estar registradas en el protocolo por plausibilidad clínica antes de comenzar con el estudio
#modelo ejemplo
modelo_sin_int <- glm(
  muerte ~ edad + fumador + diabetes,
  family = binomial,
  data = base_logistica
)
summary(modelo_sin_int)
## 
## Call:
## glm(formula = muerte ~ edad + fumador + diabetes, family = binomial, 
##     data = base_logistica)
## 
## Coefficients:
##             Estimate Std. Error z value       Pr(>|z|)    
## (Intercept) -4.07251    0.61086   -6.67 0.000000000026 ***
## edad         0.05056    0.00886    5.71 0.000000011365 ***
## fumadorSi    0.60719    0.21737    2.79         0.0052 ** 
## diabetesSi   0.54386    0.20987    2.59         0.0096 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 672.19  on 499  degrees of freedom
## Residual deviance: 622.05  on 496  degrees of freedom
## AIC: 630
## 
## Number of Fisher Scoring iterations: 4
tab_model(modelo_sin_int)
  muerte
Predictors Odds Ratios CI p
(Intercept) 0.02 0.00 – 0.05 <0.001
edad 1.05 1.03 – 1.07 <0.001
fumador [Si] 1.84 1.20 – 2.82 0.005
diabetes [Si] 1.72 1.14 – 2.60 0.010
Observations 500
R2 Tjur 0.096 
modelo_int <- glm(
  muerte ~ edad + fumador + diabetes +
    fumador:diabetes,
  family = binomial,
  data = base_logistica
)
summary(modelo_int)
## 
## Call:
## glm(formula = muerte ~ edad + fumador + diabetes + fumador:diabetes, 
##     family = binomial, data = base_logistica)
## 
## Coefficients:
##                      Estimate Std. Error z value       Pr(>|z|)    
## (Intercept)          -4.02027    0.61242   -6.56 0.000000000052 ***
## edad                  0.05066    0.00888    5.71 0.000000011582 ***
## fumadorSi             0.38222    0.26192    1.46           0.14    
## diabetesSi            0.34415    0.24671    1.39           0.16    
## fumadorSi:diabetesSi  0.74650    0.48248    1.55           0.12    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 672.19  on 499  degrees of freedom
## Residual deviance: 619.60  on 495  degrees of freedom
## AIC: 629.6
## 
## Number of Fisher Scoring iterations: 4
tab_model(modelo_int)
  muerte
Predictors Odds Ratios CI p
(Intercept) 0.02 0.01 – 0.06 <0.001
edad 1.05 1.03 – 1.07 <0.001
fumador [Si] 1.47 0.87 – 2.45 0.144
diabetes [Si] 1.41 0.87 – 2.29 0.163
fumador [Si] × diabetes
[Si] 		2.11 0.83 – 5.53 0.122
Observations 500
R2 Tjur 0.101 
#modelo_int <-glm(
#  muerte ~ edad +  fumador*diabetes,
#  family = binomial,
#  data = base_logistica
#)
anova(modelo_sin_int,
      modelo_int,
      test = "LRT")
## Analysis of Deviance Table
## 
## Model 1: muerte ~ edad + fumador + diabetes
## Model 2: muerte ~ edad + fumador + diabetes + fumador:diabetes
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       496        622                     
## 2       495        620  1     2.44     0.12

Si la p de ANOVA es menor a 0.05, quiere decir que la interacción mejora el modelo. Lo mismo para el término de interacción. Interpretación Supongamos

OR fumador = 1.47

OR interacción = 2.11 OR fumador no diabéticos: 1.47 OR fumador en diabéticos: 1.47 x 2.11 =3.1

  1. Confusión residual. Evaluar por investigadores: La confusión residual aparece por dos grandes motivos.

1.Un confundidor importante no fue medido.

2.Fue medido pero mal ajustado.(por la escala utilizada por ejemplo). No categorizar variables continuas innecesariamente. Usar splines, usar DAGS, preguntas de sensibilidad.

APARTADO DE SPLINES Es una función del paquete splines que transforma una variable continua (edad) en varias funciones matemáticas para permitir que la relación con el desenlace sea curva en lugar de obligatoriamente recta. El spline Divide el eje X (edad) en varios segmentos. df = degrees of freedom (grados de libertad). Significa cuánta flexibilidad le das al spline.DF=10 probablemente termine sobreajustando.entre 3 y 5 grados de libertad funciona muy bien para la mayoría de las variables clínicas.

##Parte II: Construcción de modelos predictivos

  1. Selección de predictores (sin caer en el error del univariado p < 0,10).Me interesa encontrar la combinación de variables que mejor prediga el desenlace en pacientes futuros. Selección univariado p<0.10 Stepwise, forward, backward Lasso Elastic net *Ridge
#métodos de selección de variables
#1. Métodos basados en conocimiento previo (preferidos): Aquí los predictores se seleccionan porque existe evidencia de que mejoran la predicción.(diabetes, edad, hta, etc)
#2. Métodos automáticos

#a. Selección por univariado (p < 0,10)

library(dplyr)
library(purrr)
library(broom)

variables <- c(
  "edad",
  "sexo",
  "fumador",
  "diabetes",
  "hta",
  "colesterol",
  "creatinina"
)

tabla_uni <- map_dfr(variables, function(x){

  form <- as.formula(
    paste("muerte ~", x)
  )

  mod <- glm(form,
             family = binomial,
             data = base_logistica)

  tidy(mod) %>%
    filter(term != "(Intercept)") %>%
    mutate(variable = x)

})

tabla_uni %>%
  dplyr::select(variable, term, estimate, std.error,
                statistic, p.value) %>%
  dplyr::arrange(p.value)
## # A tibble: 7 × 6
##   variable   term          estimate std.error statistic      p.value
##   <chr>      <chr>            <dbl>     <dbl>     <dbl>        <dbl>
## 1 edad       edad          0.0489     0.00866    5.65   0.0000000163
## 2 sexo       sexoMasculino 0.810      0.190      4.26   0.0000205   
## 3 fumador    fumadorSi     0.583      0.208      2.81   0.00500     
## 4 creatinina creatinina    0.674      0.241      2.80   0.00518     
## 5 diabetes   diabetesSi    0.510      0.200      2.55   0.0108      
## 6 hta        htaSi         0.201      0.184      1.09   0.275       
## 7 colesterol colesterol    0.000118   0.00223    0.0530 0.958

Método clásico: las variables seleccionadas se incorporan todas juntas al modelo multivariado.

#forward, backward y stepwise
# Modelo nulo
modelo_nulo <- glm(
  muerte ~ 1,
  data = base_logistica,
  family = binomial
)

# Modelo completo
modelo_completo <- glm(
  muerte ~ edad + sexo + fumador + diabetes + hta +
    colesterol + creatinina,
  data = base_logistica,
  family = binomial
)

# Stepwise en ambas direcciones
modelo_step <- step(
  modelo_nulo,
  scope = list(
    lower = modelo_nulo,
    upper = modelo_completo
  ),
  direction = "both",
  trace = TRUE
)
## Start:  AIC=674
## muerte ~ 1
## 
##              Df Deviance AIC
## + edad        1      637 641
## + sexo        1      653 657
## + creatinina  1      664 668
## + fumador     1      664 668
## + diabetes    1      666 670
## <none>               672 674
## + hta         1      671 675
## + colesterol  1      672 676
## 
## Step:  AIC=641
## muerte ~ edad
## 
##              Df Deviance AIC
## + sexo        1      616 622
## + creatinina  1      628 634
## + fumador     1      629 635
## + diabetes    1      630 636
## + hta         1      635 641
## <none>               637 641
## + colesterol  1      637 643
## - edad        1      672 674
## 
## Step:  AIC=622
## muerte ~ edad + sexo
## 
##              Df Deviance AIC
## + fumador     1      606 614
## + diabetes    1      607 615
## + creatinina  1      608 616
## <none>               616 622
## + hta         1      614 622
## + colesterol  1      616 624
## - sexo        1      637 641
## - edad        1      653 657
## 
## Step:  AIC=614
## muerte ~ edad + sexo + fumador
## 
##              Df Deviance AIC
## + diabetes    1      597 607
## + creatinina  1      597 607
## + hta         1      604 614
## <none>               606 614
## + colesterol  1      606 616
## - fumador     1      616 622
## - sexo        1      629 635
## - edad        1      644 650
## 
## Step:  AIC=607
## muerte ~ edad + sexo + fumador + diabetes
## 
##              Df Deviance AIC
## + creatinina  1      589 601
## + hta         1      595 607
## <none>               597 607
## + colesterol  1      597 609
## - diabetes    1      606 614
## - fumador     1      607 615
## - sexo        1      622 630
## - edad        1      636 644
## 
## Step:  AIC=601
## muerte ~ edad + sexo + fumador + diabetes + creatinina
## 
##              Df Deviance AIC
## <none>               589 601
## + hta         1      588 602
## + colesterol  1      589 603
## - creatinina  1      597 607
## - diabetes    1      597 607
## - fumador     1      599 609
## - sexo        1      613 623
## - edad        1      629 639
summary(modelo_step)
## 
## Call:
## glm(formula = muerte ~ edad + sexo + fumador + diabetes + creatinina, 
##     family = binomial, data = base_logistica)
## 
## Coefficients:
##               Estimate Std. Error z value          Pr(>|z|)    
## (Intercept)   -5.83315    0.75749   -7.70 0.000000000000014 ***
## edad           0.05542    0.00932    5.95 0.000000002712921 ***
## sexoMasculino  0.98438    0.20700    4.76 0.000001981054436 ***
## fumadorSi      0.69990    0.22627    3.09            0.0020 ** 
## diabetesSi     0.61890    0.21905    2.83            0.0047 ** 
## creatinina     0.72082    0.26053    2.77            0.0057 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 672.19  on 499  degrees of freedom
## Residual deviance: 589.47  on 494  degrees of freedom
## AIC: 601.5
## 
## Number of Fisher Scoring iterations: 3
# OR e IC95%
exp(cbind(
  OR = coef(modelo_step),
  confint(modelo_step)
))
## Waiting for profiling to be done...
##                    OR    2.5 % 97.5 %
## (Intercept)   0.00293 0.000633 0.0124
## edad          1.05699 1.038276 1.0770
## sexoMasculino 2.67614 1.792315 4.0392
## fumadorSi     2.01355 1.294148 3.1465
## diabetesSi    1.85688 1.210064 2.8596
## creatinina    2.05612 1.241798 3.4572

Stepwise, forward, backward actualmente se desaconseja como estrategia principal para modelos predictivos por su inestabilidad y riesgo de sobreajuste. Tenemos un modelo con muchas variables. Queremos: evitar el sobreajuste, reducir la inestabilidad de los coeficientes, si es posible, eliminar variables poco útiles.

Entonces aparecen tres métodos:Son herramientas para armar modelos predictivos, similar al stepwise. Su objetivo principal es mejorar la capacidad de generalización del modelo, es decir, que prediga bien en pacientes nuevos reduciendo el sobreajuste. Para evitar el sobreajuste aparece λ, que es un parámetro que indica en cuánto debemos reducir los coeficientes para evitar el sobreajuste. λ = 0 es igual a una reg log común (elige el programa) λ = 1 queda solo el intercepto alfa: decide que método usar 1=lasso 0=ridge (elige el investigador) lambda decide cuánto penalizar lambda.min es el que produce el menor error

           Penalización

                │
    ┌───────────┴───────────┐
    │                       │
  Ridge                  LASSO

(reduce β) (reduce β y elimina) │ │ └───────────┬───────────┘ │ Elastic Net (combina ambas estrategias)

Ridge: Reduce todos los coeficientes. No elimina ninguna variable.Muy útil cuando hay multicolinealidad.

LASSO: cuando hay muchas variables y sospecho que muchas son innecesarias. Reduce los coeficientes.Algunos llegan exactamente a cero.Esas variables desaparecen del modelo

Elastic Net: cuando tengo muchas variables que tienen mucha correlación entre sí y busca un modelo parsimonioso. Reduce todos los coeficientes.Puede eliminar algunos.

LASSO: es probablemente el concepto que más cambia la forma de pensar los modelos predictivos.agrega un castigo (penalización).Empuja los coeficientes hacia 0. Hace selección de variables de forma continua, no con decisiones abruptas basadas en un p-valor. La penalización la controla lambda. λ = 0 No hay penalización. Es una regresión logística común.A diferencia de los métodos stepwise, que seleccionan variables mediante decisiones binarias (entra/no entra), LASSO realiza una selección continua penalizando los coeficientes de regresión. Los predictores con menor contribución son progresivamente reducidos hasta que algunos alcanzan exactamente un coeficiente igual a cero, siendo excluidos del modelo. Esta estrategia suele generar modelos más estables y con menor riesgo de sobreajuste.

#LASSO 
library(glmnet)
#este es el comando que se requiere para lasso, elastic net y ridge
##################################################################
x <- model.matrix(
  muerte ~ edad + sexo + fumador +
    diabetes + hta +
    colesterol + creatinina,
  data = base_logistica
)[,-1]
y <- ifelse(base_logistica$muerte=="Si",1,0)
#################################################################
lasso <- cv.glmnet(
  x,
  y,
  family="binomial",
  alpha=1
)
plot(lasso)

coef(lasso,
     s="lambda.min") # en cero las variables eliminadas
## 8 x 1 sparse Matrix of class "dgCMatrix"
##               lambda.min
## (Intercept)     -5.54670
## edad             0.05052
## sexoMasculino    0.86521
## fumadorSi        0.60425
## diabetesSi       0.52207
## htaSi            0.16041
## colesterol       0.00111
## creatinina       0.58173
coef_lasso <- coef(lasso, s = "lambda.min")

coef_lasso
## 8 x 1 sparse Matrix of class "dgCMatrix"
##               lambda.min
## (Intercept)     -5.54670
## edad             0.05052
## sexoMasculino    0.86521
## fumadorSi        0.60425
## diabetesSi       0.52207
## htaSi            0.16041
## colesterol       0.00111
## creatinina       0.58173
tabla_lasso <- data.frame(
  Variable = rownames(as.matrix(coef_lasso)),
  Beta = as.numeric(coef_lasso)
)

tabla_lasso$OR <- exp(tabla_lasso$Beta)

tabla_lasso #Coeficientes penalizados estimados mediante LASSO
##        Variable     Beta     OR
## 1   (Intercept) -5.54670 0.0039
## 2          edad  0.05052 1.0518
## 3 sexoMasculino  0.86521 2.3755
## 4     fumadorSi  0.60425 1.8299
## 5    diabetesSi  0.52207 1.6855
## 6         htaSi  0.16041 1.1740
## 7    colesterol  0.00111 1.0011
## 8    creatinina  0.58173 1.7891
#CALCULADORA DE RIESGO
calcular_riesgo_lasso <- function(edad, sexo, fumador, diabetes, hta,
                                  colesterol, creatinina) {
  
  beta <- as.numeric(coef(lasso, s = "lambda.min"))
  names(beta) <- rownames(as.matrix(coef(lasso, s = "lambda.min")))
  
  lp <- beta["(Intercept)"] +
    beta["edad"] * edad +
    beta["sexoMasculino"] * ifelse(sexo == "Masculino", 1, 0) +
    beta["fumadorSi"] * ifelse(fumador == "Si", 1, 0) +
    beta["diabetesSi"] * ifelse(diabetes == "Si", 1, 0) +
    beta["htaSi"] * ifelse(hta == "Si", 1, 0) +
    beta["colesterol"] * colesterol +
    beta["creatinina"] * creatinina
  
  riesgo <- 1 / (1 + exp(-lp))
  
  return(riesgo)
}
#prueba de la calculadora de riesgo
calcular_riesgo_lasso(
  edad = 35, #número
  sexo = "Masculino", #Masculino/Femenino
  fumador = "Si", #Si/No
  diabetes = "Si", #Si/No
  hta = "No", #Si/No
  colesterol = 220,
  creatinina = 1.3
) * 100
## (Intercept) 
##        31.3

En el gráfico: a la izquierda, lambda grande, penaliza mucho. a la derecha, lambda más pequeño, penaliza menos. Es decir, hacia la derecha el modelo se parece cada vez más a una regresión logística clásica.En el eje de la y está el error de predicción, cuanto menor es, mejor predice Línea punteada izq: El λ que produce el menor error de predicción. Línea punteada der: el modelo más simple cuyo error todavía está dentro de un error estándar del mínimo. Los números arriba significan cuántas variables siguen vivas en el modelo. Si avanzás hacia la derecha, −log(λ) aumenta, lo que significa que λ disminuye. Y como λ es el parámetro de penalización, un λ más pequeño implica menos penalización.Dejó las 7 variables porque algo aportan. Ninguna vale 0

Ridge: Ridge no realiza selección de variables: conserva todos los predictores, pero reduce la magnitud de sus coeficientes, especialmente los menos informativos. El valor de alfa lo pone el investigadr, de acuerdo a si quiere parecerse más a lasso (0) o a Ridge (1). No importa cómo se llame el objeto (ridge, lasso, modelo, pepito…). Lo que define el método es el valor de alpha.

#ojo! antes debe realzar el comando común que define x e y
#Ridge
ridge <- cv.glmnet(
  x,
  y,
  family = "binomial",
  alpha = 1 #el valor de alfa lo elige el investigador. 
)
coef(ridge, s = "lambda.min")
## 8 x 1 sparse Matrix of class "dgCMatrix"
##               lambda.min
## (Intercept)     -6.46173
## edad             0.05596
## sexoMasculino    0.97301
## fumadorSi        0.71647
## diabetesSi       0.62426
## htaSi            0.24377
## colesterol       0.00239
## creatinina       0.67759
#elastic
elastic <- cv.glmnet(
  x,
  y,
  family = "binomial",
  alpha = 0.5
)
#lasso
elastic <- cv.glmnet(
  x,
  y,
  family = "binomial",
  alpha = 0
)
coef(lasso, s = "lambda.min")
## 8 x 1 sparse Matrix of class "dgCMatrix"
##               lambda.min
## (Intercept)     -5.54670
## edad             0.05052
## sexoMasculino    0.86521
## fumadorSi        0.60425
## diabetesSi       0.52207
## htaSi            0.16041
## colesterol       0.00111
## creatinina       0.58173
coef(elastic, s = "lambda.min")
## 8 x 1 sparse Matrix of class "dgCMatrix"
##               lambda.min
## (Intercept)     -5.66104
## edad             0.04858
## sexoMasculino    0.84760
## fumadorSi        0.62977
## diabetesSi       0.54827
## htaSi            0.21837
## colesterol       0.00195
## creatinina       0.60682
coef(ridge, s = "lambda.min")
## 8 x 1 sparse Matrix of class "dgCMatrix"
##               lambda.min
## (Intercept)     -6.46173
## edad             0.05596
## sexoMasculino    0.97301
## fumadorSi        0.71647
## diabetesSi       0.62426
## htaSi            0.24377
## colesterol       0.00239
## creatinina       0.67759

EPV (eventos por variable) y tamaño muestral en predicción: Los pacientes que no presentan el evento también aportan información, pero la cantidad de eventos es el principal determinante de la precisión de las estimaciones. 10 eventos por variable: Fue propuesta por Peduzzi y colaboradores en la década de 1990.La edad con 4 splines genera 3 grados de libertad. Lo que realmente contamos son grados de libertad (degrees of freedom) No solamente variables. 10 a 20 eventos por parámetro es aceptable. LASSO tolera EPV más bajos

#tamaño muestral
#install.packages("pmsampsize")
library(pmsampsize)


pmsampsize(
  type = "b",
  nagrsquared = 0.20,
  parameters = 8,
  prevalence = 0.20
)
## NB: Assuming 0.05 acceptable difference in apparent & adjusted R-squared 
## NB: Assuming 0.05 margin of error in estimation of intercept 
## NB: Events per Predictor Parameter (EPP) assumes prevalence = 0.2  
##  
##            Samp_size Shrinkage Parameter CS_Rsq Max_Rsq Nag_Rsq   EPP
## Criteria 1       531     0.900         8  0.126   0.632   0.199 13.28
## Criteria 2       233     0.799         8  0.126   0.632   0.199  5.82
## Criteria 3       246     0.900         8  0.126   0.632   0.199  6.15
## Final            531     0.900         8  0.126   0.632   0.199 13.28
##  
##  Minimum sample size required for new model development based on user inputs = 531, 
##  with 107 events (assuming an outcome prevalence = 0.2) and an EPP = 13.3 
##  
##

Sobreajuste (overfitting).

Validación interna (bootstrap vs. cross-validation). Discriminación. Calibración. Optimismo. Presentación del modelo (nomograma, calculadora de riesgo, etc.).