Regresión Logistica

Independencia

Si el hecho de que aparezca o se presente una categoría junto con otra no es ni más ni menos probable de que se presenten las dos categorías por separado, se dice que las variables son independientes y, en general, se dice que la tabla es homogénea.

Dadas dos variables aleatorias \(X\) e \(Y\) , son independientes si

\[P\left(X=x_{i}, Y=y_{j}\right)=P\left(X=x_{i}\right) \cdot P\left(Y=y_{j}\right)\]

para todo (\(i,j\))

\[ \begin{aligned} p_{i j} &=\frac{n_{i j}}{n . .} \\ p_{i .} &=\frac{n_{i .}}{n . .} \\ p_{. j} &=\frac{n_{. j}}{n . .} \end{aligned} \] Así, si \[ P\left(X=x_{i}, Y=y_{j}\right)=p_{i j}=p_{i} \cdot \times p_{\cdot j} \]

para todo \(i\), \(j\), las variables \(X\) e \(Y\) son independientes y la tabla es homogénea

En el caso de de ser cierta la hipótesis de independencia esperaremos encontrar \(E_{ij}\) objetos dentro de la casilla \((i, j)- ésima\), donde

\[ E_{i j}=n.. p_{i j}=n.. p_{i\cdot} p_{\cdot j}=\frac{n_{i\cdot} n_{\cdot j}}{n_{. .}} \]

library("factoextra")
library("FactoMineR")
library("gplots")
library("dplyr")

library(sjPlot)
library(sjmisc)
library(sjlabelled)

O=matrix(c(11,3,8,2,9,14,12,13,28), nrow = 3, byrow = T)

colnames(O)=c("Mano der", "Mano izq", "Manos iguales")
rownames(O)=c("pie der", "pie izq", "pies iguales")

library(kableExtra)

kable(addmargins(O))

O_chisq = chisq.test(O)

kable(addmargins(O_chisq$expected))

par(mfrow = c(1, 2))

mosaicplot(O,  col = 3:5, main= "Valores observados",)
mosaicplot(O_chisq$expected, col = 3:5, main= "Valores esperados")

par(mfrow = c(1, 1))

O_chisq = chisq.test(O)
O_chisq

## 
##  Pearson's Chi-squared test
## 
## data:  O
## X-squared = 11.942, df = 4, p-value = 0.01779

Regresión Logistica

El objetivo de este procedimiento es establecer el mejor modelo de explicación de una variable dicotómica cualitativa \(Y [0,1]\), denominada variable de respuesta o dependiente, en base a otra serie. de variables denominadas covariables, variables predictoras o simplemente variables independientes \(X_1\) , \(X_2\) , \(X_3\) , …, \(X_m\) que indican las características del sujeto, y que pueden ser: - Discretas - Continuas.

Cuando la variable respuesta \(Y\) es de tipo dicotómico nos referimos a regresión logística, y cuando la variable dependiente \(Y\) es cualitativa con más de dos categorías (policotómica), nos referimos a regresión logística multinomial.

Además de evitar las limitaciones de la regresión lineal cuando la variable de resultado es dicotómica, esta técnica permite interpretar los parámetros de forma sencilla en términos de odds ratio (OR).

Definiciones previas

Regresión logística binaria

La definición del mejor modelo depende del tipo y objetivo del estudio.

El modelo suele tener dos tipos de objetivo:

• predictivo o

• explicativo.

  • En un modelo con objetivos predictivos pretendemos establecer un modelo parsimonioso, es decir, un modelo que involucre el menor número de variables que mejor expliquen la variable dependiente.

  • En el caso de un modelo con objetivos explicativos, se pretende estudiar la relación causal entre una variable “causa” y una variable “efecto”, con el debido control de las posibles variables de confusión (definidas en el Apartado ‘Interacción y factores de confusión’) o variables modificadoras del efecto (interacción) en esta relación causal.

Construcción del modelo:

sigmoid <- function(x) {
  1 / (1 + exp(-x))
}

curve(sigmoid, from = -10, to = 10, n = 100, xlab = "x", ylab = "sigmoid(x)", main = "Función Sigmoide")

Codificación e interpretación de los coeficientes del modelo

Coeficientes del modelo logístico como cuantificadores de riesgo

Es recomendable seguir las siguientes recomendaciones a la hora de codificar las variables de un modelo de regresión logística, ya que facilita la interpretación de los resultados:

• Variable dependiente : codificamos como 1 la ocurrencia del evento de interés y como 0 la ausencia del evento.

Ejemplo:

• 1: enfermedad sí,

• 0: enfermedad no.

• Variables independientes : estas pueden ser de diferentes tipos:

  • Variable numérica: para introducir la variable en el modelo, debe satisfacer la hipótesis de linealidad, es decir, por cada unidad de incremento en la variable numérica, el OR (\(exp(B_j\)) aumenta en un valor multiplicativo constante. En este caso podemos utilizar la variable tal y como está en el modelo. Si no se cumple la hipótesis de linealidad, podemos transformar la variable numérica, categorizándola. Cuando la variable es continua, la OR asociada se interpreta como el incremento del riesgo por incremento unitario en la covariable analizada, es decir, si el factor de riesgo fuera la edad, la variable resultado sería tener o no la enfermedad, y si la OR fuera 1.87, entonces este valor se interpretaría como que indica que por cada año adicional de edad del sujeto, el riesgo de padecer la enfermedad aumenta en 1,87.

  Es decir, la probabilidad de padecer la enfermedad aumenta un 87% por cada año adicional de edad, ajustando por el resto de variables del modelo en el caso multivariado.

  • Variable dicotómica : codificamos como 1 el caso que se cree que favorece la ocurrencia del evento, mientras que el caso contrario se codifica como 0, es decir, codificamos como 1 aquellos individuos expuestos a un supuesto factor de riesgo, y como 0 aquellos individuos no expuestos al mismo. factor mencionado. Así, si codificamos a las madres fumadoras como 1 y a las madres no fumadoras como 0, y se toma como evento de interés un niño asmático (codificado como 1) versus un niño no asmático (codificado como 0), y obtenemos OR=2.23, entonces se puede concluir que las madres fumadoras tienen un riesgo 2.23 veces mayor de tener un hijo asmático que las madres no fumadoras.

  Alternativamente, esto puede interpretarse como que el riesgo de tener un hijo con asma en el caso de una madre fumadora es algo más del doble que en el caso de una madre no fumadora, ajustando en todos los casos por el resto de covariables del modelo.

  • Variable categórica: este tipo de variable se divide en diferentes variables dicotómicas que representan las distintas categorías. Estas variables se conocen como indicadores, variables internas, variables de diseño o variables “ficticias”.

  La mayoría de las aplicaciones de software estadístico generan internamente estas variables al incluir una variable de tipo factorial en el modelo (es decir, con más de dos categorías). En este tipo de codificación solemos seleccionar una categoría de referencia, y el coeficiente de la ecuación de regresión para cada “dummy La variable” (siempre transformada con la función exponencial) corresponde al OR de esa categoría con respecto al nivel de referencia.

  Así, cuantificamos el cambio en el riesgo de sufrir el evento de cada una de las categorías con respecto a la categoría de referencia.

  

Como regla general, si una variable cualitativa tiene \(k\) niveles, se necesitarán de \(k-1\) variables indicadoras para resumir la variable cualitativa.

Conjunto de datos

Estos datos provienen de un estudio de pacientes con accidente cerebrovascular hospitalizados. El conjunto de datos original contiene 12 variables, pero nuestras principales variables de interés son:

• status : Estado del paciente durante la hospitalización (vivo o muerto)
• gcs: Escala de coma de Glasgow al ingreso (rango de 3 a 15)
• stroke_type : IS (ACV isquémico) o HS (ACV hemorrágico) sexo: femenino o masculino
• dm : Antecedentes de Diabetes (sí o no)
• sbp : Presión arterial sistólica (mmHg)
• edad : edad del paciente al ingreso

La variable objetivo es status . Se etiqueta como vivo o muerto, que es el resultado de cada paciente durante la hospitalización.

library(haven)
library(tidyverse)
library(gtsummary)
library(broom)
#library(LogisticDx)
library(here)
library(kableExtra)
library(rio)
library(DT)

fatal=rio::import("https://github.com/Wilsonsr/Metodos-Estadisticos/raw/main/BASES/basestronk.csv")

DT::datatable(fatal)

glimpse(fatal)

## Rows: 226
## Columns: 7
## $ sex         <chr> "Masculino", "Masculino", "Masculino", "Femenino", "Mascul…
## $ status      <chr> "Vivo", "Vivo", "Vivo", "Vivo", "Vivo", "Vivo", "Muerto", …
## $ gcs         <int> 13, 15, 15, 15, 15, 15, 13, 15, 15, 10, 15, 15, 15, 15, 15…
## $ sbp         <int> 143, 150, 152, 215, 162, 169, 178, 180, 186, 185, 122, 211…
## $ dm          <chr> "NO", "NO", "NO", "SI", "SI", "SI", "SI", "NO", "SI", "SI"…
## $ age         <int> 50, 58, 64, 50, 65, 78, 66, 72, 61, 64, 63, 59, 64, 62, 40…
## $ stroke_type <chr> "IS", "IS", "IS", "IS", "IS", "IS", "IS", "IS", "IS", "IS"…

fatal$status<-as.factor(fatal$status)
fatal$sex <-as.factor(fatal$sex)
fatal$dm <-as.factor(fatal$dm)
fatal$stroke_type <-as.factor(fatal$stroke_type)

fatal %>%
  tbl_summary() %>%
  as_gt()

fatal %>%
  tbl_summary(by = status) %>%
  as_gt()

`

#library(ExPanDaR)
## Create report
#ExPanD(fatal)

str(fatal)

## 'data.frame':    226 obs. of  7 variables:
##  $ sex        : Factor w/ 2 levels "Femenino","Masculino": 2 2 2 1 2 1 1 2 1 1 ...
##  $ status     : Factor w/ 2 levels "Muerto","Vivo": 2 2 2 2 2 2 1 2 2 1 ...
##  $ gcs        : int  13 15 15 15 15 15 13 15 15 10 ...
##  $ sbp        : int  143 150 152 215 162 169 178 180 186 185 ...
##  $ dm         : Factor w/ 2 levels "NO","SI": 1 1 1 2 2 2 2 1 2 2 ...
##  $ age        : int  50 58 64 50 65 78 66 72 61 64 ...
##  $ stroke_type: Factor w/ 2 levels "HS","IS": 2 2 2 2 2 2 2 2 2 2 ...

variables_categoricas <- sapply(fatal, is.factor) | sapply(fatal, is.character)
datos_categoricos <- fatal[, variables_categoricas]

variables <- names(datos_categoricos)
combinaciones <- combn(variables, m = 2, simplify = FALSE)

resultados_chi2 <- list()

for (i in 1:length(combinaciones)) {
  tabla_contingencia <- table(datos_categoricos[, combinaciones[[i]]])
  resultado_chi2 <- chisq.test(tabla_contingencia)
  resultados_chi2[[i]] <- resultado_chi2
}

for (i in 1:length(resultados_chi2)) {
  print(paste("Variables:", combinaciones[[i]]))
  print( resultados_chi2[[i]])
  print("-----------------------------------------------------")
  print("-----------------------------------------------------")
  
}

## [1] "Variables: sex"    "Variables: status"
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabla_contingencia
## X-squared = 4.9531, df = 1, p-value = 0.02604
## 
## [1] "-----------------------------------------------------"
## [1] "-----------------------------------------------------"
## [1] "Variables: sex" "Variables: dm" 
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabla_contingencia
## X-squared = 0.04049, df = 1, p-value = 0.8405
## 
## [1] "-----------------------------------------------------"
## [1] "-----------------------------------------------------"
## [1] "Variables: sex"         "Variables: stroke_type"
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabla_contingencia
## X-squared = 1.0126, df = 1, p-value = 0.3143
## 
## [1] "-----------------------------------------------------"
## [1] "-----------------------------------------------------"
## [1] "Variables: status" "Variables: dm"    
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabla_contingencia
## X-squared = 1.5498, df = 1, p-value = 0.2132
## 
## [1] "-----------------------------------------------------"
## [1] "-----------------------------------------------------"
## [1] "Variables: status"      "Variables: stroke_type"
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabla_contingencia
## X-squared = 37.653, df = 1, p-value = 8.45e-10
## 
## [1] "-----------------------------------------------------"
## [1] "-----------------------------------------------------"
## [1] "Variables: dm"          "Variables: stroke_type"
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabla_contingencia
## X-squared = 0.5105, df = 1, p-value = 0.4749
## 
## [1] "-----------------------------------------------------"
## [1] "-----------------------------------------------------"

fatal_glm_1 <-glm(status ~ gcs, data = fatal, family = binomial(link = 'logit'))
resumen=summary(fatal_glm_1)
tabla <- resumen$table
tabla_formateada <- kable(tabla)

# Imprimir la tabla formateada
tabla_formateada

summary(fatal_glm_1)

## 
## Call:
## glm(formula = status ~ gcs, family = binomial(link = "logit"), 
##     data = fatal)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.29479    0.60432  -5.452 4.98e-08 ***
## gcs          0.38811    0.05213   7.446 9.64e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 250.83  on 225  degrees of freedom
## Residual deviance: 170.92  on 224  degrees of freedom
## AIC: 174.92
## 
## Number of Fisher Scoring iterations: 5

library(sjPlot)
library(sjmisc)
library(sjlabelled)


fatal_glm_1 <-glm(status ~ gcs, data = fatal, family = binomial(link = 'logit'))
summary(fatal_glm_1)

## 
## Call:
## glm(formula = status ~ gcs, family = binomial(link = "logit"), 
##     data = fatal)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.29479    0.60432  -5.452 4.98e-08 ***
## gcs          0.38811    0.05213   7.446 9.64e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 250.83  on 225  degrees of freedom
## Residual deviance: 170.92  on 224  degrees of freedom
## AIC: 174.92
## 
## Number of Fisher Scoring iterations: 5

tab_model(fatal_glm_1, show.est = TRUE, show.se = TRUE, show.std = TRUE,  collapse.ci = FALSE)

tab_model(fatal_glm_1)

tidy(fatal_glm_1, conf.int = TRUE)

## # A tibble: 2 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)   -3.29     0.604      -5.45 4.98e- 8   -4.55     -2.17 
## 2 gcs            0.388    0.0521      7.45 9.64e-14    0.292     0.497

• Coeficiente del intercepto (Intercept): El coeficiente del intercepto es -3.29479. En un modelo de regresión logística, el intercepto representa el logaritmo de la razón de probabilidades (log-odds) de que la respuesta sea “vivo” en ausencia de cualquier efecto de las variables predictoras.

Para interpretar en términos de probabilidades,se aplicar la función exponencial al coeficiente del intercepto. En este caso, sería exp(-3.29479), lo cual es aproximadamente 0.036. Esto significa que, cuando todas las variables predictoras son cero, la probabilidad de que la respuesta sea “vivo” es de aproximadamente 0.036, o alrededor del 3.6%.

• Coeficiente de la variable gcs: El coeficiente para la variable gcs es 0.38811. Este coeficiente indica cómo cambia la razón de probabilidades (odds ratio) de la respuesta “vivo” por cada unidad de cambio en la variable “gcs”.

Para interpretarlo en términos de la dirección del efecto, un coeficiente positivo significa que un aumento en la variable “gcs” está asociado con un aumento en la probabilidad de que la respuesta sea “vivo”. Por ejemplo, si “gcs” aumenta en una unidad, la razón de probabilidades de estar “vivo” aumentará en aproximadamente exp(0.38811) = 1.474 veces.

Es importante tener en cuenta que los coeficientes en un modelo de regresión logística se interpretan en términos de log-odds o razones de probabilidades. Para obtener interpretaciones más intuitivas en términos de probabilidades, puedes aplicar funciones exponenciales a los coeficientes. Además, es recomendable considerar los intervalos de confianza y los p-valores asociados para evaluar la significancia estadística de los coeficientes y la incertidumbre en las estimaciones.

Ahora, usemos otra covariable, stroke_type. El tipo de trazo tiene 2 niveles o categorías; Accidente cerebrovascular hemorrágico (HS) y Accidente cerebrovascular isquémico (IS). Se sabe que la HS causa un mayor riesgo de muerte por accidente cerebrovascular.

fatal_glm_2 <-  glm(status ~ stroke_type,   data = fatal,   family = binomial(link = 'logit'))

summary(fatal_glm_2)

## 
## Call:
## glm(formula = status ~ stroke_type, family = binomial(link = "logit"), 
##     data = fatal)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.02598    0.22794   0.114    0.909    
## stroke_typeIS  2.02361    0.34403   5.882 4.05e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 250.83  on 225  degrees of freedom
## Residual deviance: 212.52  on 224  degrees of freedom
## AIC: 216.52
## 
## Number of Fisher Scoring iterations: 4

kable(tidy(fatal_glm_2, conf.int = TRUE)) %>% kable_styling("striped")

• Coeficiente del intercepto (Intercept): El coeficiente del intercepto es 0.02598. En este modelo, el intercepto representa el logaritmo de la razón de probabilidades (log-odds) de que la respuesta sea “vivo” (o la categoría de referencia) cuando “stroke_type” es el valor de referencia. Sin embargo, en este caso, el valor \(p\) para el coeficiente del intercepto es alto (0.909), lo que indica que no hay evidencia estadística significativa de que haya una diferencia en las probabilidades de estar “vivo” en función del tipo de accidente cerebrovascular.

• Coeficiente para la variable stroke_typeIS: El coeficiente para la categoría “IS” de la variable “stroke_type” es 2.02361. Este coeficiente indica cómo cambia la razón de probabilidades (odds ratio) de la respuesta “vivo” en comparación con el tipo de accidente cerebrovascular de referencia (en este caso, “HS”).

El coeficiente para la variable “stroke_typeIS” es 2.02361, lo que significa que la razón de probabilidades (odds ratio) de estar “vivo” en el caso de un accidente cerebrovascular isquémico (IS) en comparación con un accidente cerebrovascular hemorrágico (HS) es exp(2.02361) = 7.559.

Esto indica que los pacientes con un accidente cerebrovascular isquémico (IS) tienen una razón de probabilidades 7.559 veces mayor de estar “vivos” en comparación con aquellos con un accidente cerebrovascular hemorrágico (HS), manteniendo constantes las demás variables en el modelo.

En este caso, el valor p para el coeficiente de “stroke_typeIS” es muy bajo (< 0.001), lo que indica que hay evidencia estadística significativa de que el tipo de accidente cerebrovascular isquémico (IS) tiene un efecto significativo en la probabilidad de estar “vivo” en comparación con el tipo de accidente cerebrovascular hemorrágico (HS).

Regresión logística binaria múltiple

Hay múltiples factores que pueden contribuir al resultado de un accidente cerebrovascular. Por lo tanto, existe una fuerte motivación para incluir otras variables independientes o covariables en el modelo. Por ejemplo, en el caso de un accidente cerebrovascular:

• Es poco probable que solo una variable (gcs o tipo de ictus) esté relacionada con el ictus. El accidente cerebrovascular, al igual que otras enfermedades cardiovasculares, tiene muchos factores que afectan el resultado. Tiene más sentido considerar agregar otras variables independientes que creemos que son variables independientes importantes para el resultado del accidente cerebrovascular en el modelo.

al agregar más covariables en el modelo, podemos estimar las probabilidades logarítmicas ajustadas. Estas probabilidades logarítmicas indican la relación de una covariable particular independiente de otras covariables en el modelo. En epidemiología, siempre podemos esto como ajuste. Un ajuste es importante particularmente cuando tenemos efectos de confusión de otras variables independientes. Se puede generar un término de interacción (el producto de dos covariables) y agregarlo al modelo que se va a estimar. Agregar o no agregar variables es un gran tema en sí mismo. Por lo general, se rige por la experiencia clínica, los expertos en la materia y algunos análisis preliminares. Puede leer el último capítulo del libro para comprender más sobre la construcción de modelos y la selección de variables.

Ampliemos nuestro modelo e incluyamos gcs, tipo de brazada, sexo, dm, sbp y edad en el modelo. Llamaremos a este modelo como fatal_mv. Como tenemos más de una variable independiente en el modelo, lo llamaremos regresión logística binaria múltiple.

Para estimar el modelo de regresión logística múltiple en R:

fatal_mv1 <-   glm(status ~ gcs + stroke_type + sex + dm + sbp + age, data = fatal,     family = binomial(link = 'logit'))
summary(fatal_mv1)

## 
## Call:
## glm(formula = status ~ gcs + stroke_type + sex + dm + sbp + age, 
##     family = binomial(link = "logit"), data = fatal)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.5377396  1.4658253  -1.049  0.29415    
## gcs            0.3284640  0.0557574   5.891 3.84e-09 ***
## stroke_typeIS  1.2662764  0.4365882   2.900  0.00373 ** 
## sexMasculino   0.4302901  0.4362742   0.986  0.32399    
## dmSI          -0.4736670  0.4362309  -1.086  0.27756    
## sbp           -0.0008612  0.0060619  -0.142  0.88703    
## age           -0.0242321  0.0154010  -1.573  0.11562    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 250.83  on 225  degrees of freedom
## Residual deviance: 159.34  on 219  degrees of freedom
## AIC: 173.34
## 
## Number of Fisher Scoring iterations: 5

log_odds <- tidy(fatal_mv1, conf.int = TRUE)
log_odds

## # A tibble: 7 × 7
##   term           estimate std.error statistic       p.value conf.low conf.high
##   <chr>             <dbl>     <dbl>     <dbl>         <dbl>    <dbl>     <dbl>
## 1 (Intercept)   -1.54       1.47       -1.05  0.294          -4.46     1.33   
## 2 gcs            0.328      0.0558      5.89  0.00000000384   0.224    0.444  
## 3 stroke_typeIS  1.27       0.437       2.90  0.00373         0.411    2.13   
## 4 sexMasculino   0.430      0.436       0.986 0.324          -0.420    1.30   
## 5 dmSI          -0.474      0.436      -1.09  0.278          -1.35     0.368  
## 6 sbp           -0.000861   0.00606    -0.142 0.887          -0.0129   0.0110 
## 7 age           -0.0242     0.0154     -1.57  0.116          -0.0555   0.00520

1. Coeficiente del intercepto (Intercept): El coeficiente del intercepto es -1.5377396. Este coeficiente representa el logaritmo de la razón de probabilidades (log-odds) de estar “vivo” en ausencia de cualquier efecto de las variables predictoras. Sin embargo, en este caso, el valor p para el coeficiente del intercepto es alto (0.29415), lo que indica que no hay evidencia estadística significativa de que haya una diferencia en las probabilidades de estar “vivo” en función de las demás variables en el modelo.

2. Coeficiente para la variable “gcs”: El coeficiente para la variable “gcs” es 0.3284640. Este coeficiente indica cómo cambia la razón de probabilidades (odds ratio) de estar “vivo” por cada unidad de aumento en el valor de la variable “gcs”. Un coeficiente positivo indica que un mayor valor de “gcs” se asocia con una mayor probabilidad de estar “vivo”.

3. Coeficiente para la variable “stroke_typeIS”: El coeficiente para la variable “stroke_typeIS” es 1.2662764. Este coeficiente indica cómo cambia la razón de probabilidades (odds ratio) de estar “vivo” en comparación con el tipo de accidente cerebrovascular de referencia (en este caso, “HS”). Un coeficiente positivo indica que tener un accidente cerebrovascular isquémico (IS) en lugar de un accidente cerebrovascular hemorrágico (HS) se asocia con una mayor probabilidad de estar “vivo”.

4. Coeficiente para la variable “sexMasculino”: El coeficiente para la variable “sexMasculino” es 0.4302901. Este coeficiente indica cómo cambia la razón de probabilidades (odds ratio) de estar “vivo” para los pacientes masculinos en comparación con los pacientes femeninos. Un coeficiente positivo indica que los pacientes masculinos tienen una mayor probabilidad de estar “vivos” en comparación con los pacientes femeninos.

5. Coeficiente para la variable “dmSI”: El coeficiente para la variable “dmSI” es -0.4736670. Este coeficiente indica cómo cambia la razón de probabilidades (odds ratio) de estar “vivo” para los pacientes con diabetes en comparación con los pacientes sin diabetes. Un coeficiente negativo indica que los pacientes con diabetes tienen una menor probabilidad de estar “vivos” en comparación con los pacientes sin diabetes.

6. Coeficiente para la variable “sbp”: El coeficiente para la variable “sbp” es -0.0008612. Este coeficiente indica cómo cambia la razón de probabilidades (odds ratio) de estar “vivo” por cada unidad de cambio en la presión arterial sistólica (sbp). En este caso, el coeficiente es muy pequeño y no significativo (p > 0.05), lo que sugiere que la presión arterial sistólica no tiene un efecto significativo en la probabilidad de estar “vivo”.

7. Coeficiente para la variable “age”: El coeficiente para la variable “age” es -0.0242321. Este coeficiente indica cómo cambia la razón de probabilidades (odds ratio) de estar “vivo” por cada unidad de aumento en la edad. Un coeficiente negativo indica que a medida que aumenta la edad, la probabilidad de estar “vivo” disminuye. Sin embargo, en este caso, el valor p para el coeficiente de “age” no es significativo (p > 0.05), lo que sugiere que la edad no tiene un efecto estadísticamente significativo en la probabilidad de estar “vivo”.

odds_ratio <- tidy(fatal_mv1,
                   exponentiate = TRUE,  
                   conf.int = TRUE)
odds_ratio

## # A tibble: 7 × 7
##   term          estimate std.error statistic       p.value conf.low conf.high
##   <chr>            <dbl>     <dbl>     <dbl>         <dbl>    <dbl>     <dbl>
## 1 (Intercept)      0.215   1.47       -1.05  0.294           0.0116      3.78
## 2 gcs              1.39    0.0558      5.89  0.00000000384   1.25        1.56
## 3 stroke_typeIS    3.55    0.437       2.90  0.00373         1.51        8.45
## 4 sexMasculino     1.54    0.436       0.986 0.324           0.657       3.69
## 5 dmSI             0.623   0.436      -1.09  0.278           0.258       1.45
## 6 sbp              0.999   0.00606    -0.142 0.887           0.987       1.01
## 7 age              0.976   0.0154     -1.57  0.116           0.946       1.01

tab_logistic <- bind_cols(log_odds, odds_ratio) 

## New names:
## • `term` -> `term...1`
## • `estimate` -> `estimate...2`
## • `std.error` -> `std.error...3`
## • `statistic` -> `statistic...4`
## • `p.value` -> `p.value...5`
## • `conf.low` -> `conf.low...6`
## • `conf.high` -> `conf.high...7`
## • `term` -> `term...8`
## • `estimate` -> `estimate...9`
## • `std.error` -> `std.error...10`
## • `statistic` -> `statistic...11`
## • `p.value` -> `p.value...12`
## • `conf.low` -> `conf.low...13`
## • `conf.high` -> `conf.high...14`

tab_logistic %>% 
  select(term...1, estimate...2, std.error...3, 
         estimate...9, conf.low...13, conf.high...14 ,p.value...5) %>%
  rename(covariate = term...1, 
         log_odds = estimate...2,
         SE = std.error...3,
         odds_ratio = estimate...9,
         lower_OR = conf.low...13, 
         upper_OR = conf.high...14,
         p.val = p.value...5) 

## # A tibble: 7 × 7
##   covariate      log_odds      SE odds_ratio lower_OR upper_OR         p.val
##   <chr>             <dbl>   <dbl>      <dbl>    <dbl>    <dbl>         <dbl>
## 1 (Intercept)   -1.54     1.47         0.215   0.0116     3.78 0.294        
## 2 gcs            0.328    0.0558       1.39    1.25       1.56 0.00000000384
## 3 stroke_typeIS  1.27     0.437        3.55    1.51       8.45 0.00373      
## 4 sexMasculino   0.430    0.436        1.54    0.657      3.69 0.324        
## 5 dmSI          -0.474    0.436        0.623   0.258      1.45 0.278        
## 6 sbp           -0.000861 0.00606      0.999   0.987      1.01 0.887        
## 7 age           -0.0242   0.0154       0.976   0.946      1.01 0.116

Bondad de ajuste en R

anova( fatal_glm_1, fatal_mv1, test = 'Chisq')

## Analysis of Deviance Table
## 
## Model 1: status ~ gcs
## Model 2: status ~ gcs + stroke_type + sex + dm + sbp + age
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       224     170.92                       
## 2       219     159.34  5   11.582  0.04098 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fatal_mv2 <-   glm(status ~ gcs + stroke_type + age,       data = fatal,       family = binomial(link = 'logit'))

anova( fatal_mv1, fatal_mv2, test = 'Chisq')

## Analysis of Deviance Table
## 
## Model 1: status ~ gcs + stroke_type + sex + dm + sbp + age
## Model 2: status ~ gcs + stroke_type + age
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       219     159.34                     
## 2       222     161.51 -3  -2.1743    0.537

fatal_mv2_ia <-   glm(status ~ gcs + stroke_type + stroke_type:gcs + age, data = fatal,   family = binomial(link = 'logit'))
tidy(fatal_mv2_ia)

## # A tibble: 5 × 5
##   term              estimate std.error statistic    p.value
##   <chr>                <dbl>     <dbl>     <dbl>      <dbl>
## 1 (Intercept)        -2.12      1.15      -1.84  0.0662    
## 2 gcs                 0.355     0.0775     4.58  0.00000468
## 3 stroke_typeIS       1.61      1.30       1.24  0.217     
## 4 age                -0.0236    0.0147    -1.60  0.109     
## 5 gcs:stroke_typeIS  -0.0347    0.111     -0.312 0.755

prob_mv2 <-  augment(fatal_mv2, type.predict = "response")

DT::datatable(prob_mv2)

library(sjPlot)
plot_model(fatal_mv1)

library(stargazer)

## 
## Please cite as:

##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer

stargazer(fatal_glm_1, fatal_mv1, 
          type = 'text', header=FALSE,
          apply.coef = exp,
          apply.se   = exp)

## 
## ==============================================
##                       Dependent variable:     
##                   ----------------------------
##                              status           
##                        (1)            (2)     
## ----------------------------------------------
## gcs                   1.474          1.389    
##                      (1.054)        (1.057)   
##                                               
## stroke_typeIS                       3.548**   
##                                     (1.547)   
##                                               
## sexMasculino                         1.538    
##                                     (1.547)   
##                                               
## dmSI                                 0.623    
##                                     (1.547)   
##                                               
## sbp                                  0.999    
##                                     (1.006)   
##                                               
## age                                  0.976    
##                                     (1.016)   
##                                               
## Constant              0.037          0.215    
##                      (1.830)        (4.331)   
##                                               
## ----------------------------------------------
## Observations           226            226     
## Log Likelihood       -85.460        -79.669   
## Akaike Inf. Crit.    174.920        173.337   
## ==============================================
## Note:              *p<0.1; **p<0.05; ***p<0.01