Taller regresión logística y linear

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

SetWD

setwd("~/MAESTRIA EPIDEMIOLOGIA ICESI/recoleccion_d/taller_s6")

Llamar paquetes

library(readxl)

## Warning: package 'readxl' was built under R version 4.5.3

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.5.3

## 
## Adjuntando el paquete: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.5.3

library(moments)

## Warning: package 'moments' was built under R version 4.5.2

library(scales)

## Warning: package 'scales' was built under R version 4.5.3

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.5.3

## Warning: package 'tibble' was built under R version 4.5.2

## Warning: package 'tidyr' was built under R version 4.5.3

## Warning: package 'readr' was built under R version 4.5.2

## Warning: package 'purrr' was built under R version 4.5.2

## Warning: package 'stringr' was built under R version 4.5.2

## Warning: package 'forcats' was built under R version 4.5.2

## Warning: package 'lubridate' was built under R version 4.5.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ lubridate 1.9.5     ✔ tibble    3.3.1
## ✔ purrr     1.2.1     ✔ tidyr     1.3.2
## ✔ readr     2.2.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ readr::col_factor() masks scales::col_factor()
## ✖ purrr::discard()    masks scales::discard()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::lag()        masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(survival)

## Warning: package 'survival' was built under R version 4.5.3

library(dagitty)

## Warning: package 'dagitty' was built under R version 4.5.3

library(psych)

## Warning: package 'psych' was built under R version 4.5.3

## 
## Adjuntando el paquete: 'psych'
## 
## The following objects are masked from 'package:scales':
## 
##     alpha, rescale
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(epiR)

## Warning: package 'epiR' was built under R version 4.5.2

## Package epiR 2.0.91 is loaded
## Type help(epi.about) for summary information
## Type browseVignettes(package = 'epiR') to learn how to use epiR for applied epidemiological analyses

library(epitools)

## 
## Adjuntando el paquete: 'epitools'
## 
## The following object is masked from 'package:survival':
## 
##     ratetable

library(e1071)

## Warning: package 'e1071' was built under R version 4.5.3

## 
## Adjuntando el paquete: 'e1071'
## 
## The following objects are masked from 'package:moments':
## 
##     kurtosis, moment, skewness
## 
## The following object is masked from 'package:ggplot2':
## 
##     element

library(moments)
library(nortest)

## Warning: package 'nortest' was built under R version 4.5.2

library (ggplot2)
library(ggdist)

## Warning: package 'ggdist' was built under R version 4.5.3

library(ggrain)

## Warning: package 'ggrain' was built under R version 4.5.3

## Registered S3 methods overwritten by 'ggpp':
##   method                  from   
##   heightDetails.titleGrob ggplot2
##   widthDetails.titleGrob  ggplot2

library(tidyr)
library(broom)

## Warning: package 'broom' was built under R version 4.5.2

Cargar base de datos

library(readr)
data <- read_csv("data.csv")

## Rows: 1491 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Group
## dbl (7): pt.id, Sexe, age, hx_diabete, time, cTnT, creat
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

View(data)

REGRESIÓN LINEAL (Troponinas)

EJERCICIO 1

Rain cloud

ggplot(data, aes(x = Group, y = cTnT, fill = Group)) +
  geom_rain(
    alpha = 0.6,
    rain.side = "f",
    boxplot.args = list(width = 0.15, outlier.shape = NA),
    violin.args = list(alpha = 0.5)
  ) +
  geom_hline(yintercept = 14, linetype = "dashed") +
  coord_flip() +
  theme_minimal() +
  labs(
    title = "Raincloud plot de cTnT por Group ",
    x = "cTnT",
    y = "cTnT"
  )

## Warning: Option rain.side 'flanking' is being used with a side argument in violin.args.pos!!!
##       
## 
##       If you want the nudging position defaults for a flanking 1-by-1 raincloud use (rain.side = 'f1x1')
## 
##       If you want the nudging position defaults for a flanking 2-by-2 raincloud use (rain.side = 'f2x2')
## 
## 
##       Now defaulting to a 2-by-2

## Warning: Duplicated aesthetics after name standardisation: width

## Warning: Using the `size` aesthetic with geom_polygon was deprecated in ggplot2 3.4.0.
## ℹ Please use the `linewidth` aesthetic instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

En este gráfico presenta que la mayor densidad de observaciones de troponina en los 4 grupos se encuentra debajo del umbral de 14. Los boxplot, sugieren que la tendencia central se concentra en menores que 10 para cTnT. De igual forma, se observa mayor dispersión en aquellos con antecedentes vasculares, con valores extremos en el grupo de angina inestable y los pacientes con múltiples eventos.

DIABETES

data$hx_diabete <- trimws(as.character(data$hx_diabete))

data$hx_diabete <- factor(
  data$hx_diabete,
  levels = c("0", "1"),
  labels = c("Sin diabetes", "Con diabetes")
)

table(data$hx_diabete, useNA = "ifany")

## 
## Sin diabetes Con diabetes 
##         1210          281

p2 <- ggplot(data, aes(x = cTnT, y = hx_diabete, fill = hx_diabete)) +
  geom_rain(
    alpha = 0.6,
    rain.side = "f",
    boxplot.args = list(width = 0.15, outlier.shape = NA)
  ) +
  geom_vline(xintercept = 14, linetype = "dashed") +
  scale_fill_manual(values = c("Sin diabetes" = "grey", "Con diabetes" = "pink")) +
  labs(
    title = "Distribución de cTnT por Diabetes",
    x = "cTnT",
    y = "Diabetes",
    fill = "Diabetes"
  ) +
  theme_minimal()

## Warning: Option rain.side 'flanking' is being used with a side argument in violin.args.pos!!!
##       
## 
##       If you want the nudging position defaults for a flanking 1-by-1 raincloud use (rain.side = 'f1x1')
## 
##       If you want the nudging position defaults for a flanking 2-by-2 raincloud use (rain.side = 'f2x2')
## 
## 
##       Now defaulting to a 2-by-2

## Warning: Duplicated aesthetics after name standardisation: width

p2

## Warning in compute_layer(..., self = self): Argument 'x' longer than data: some
## values dropped!

El raincloud plot sugiere que la cTnT tiende a ser más alta y más dispersa en pacientes con diabetes que en aquellos sin diabetes. Se observa que la mediana de los datos en ambos grupos es menor a 10. Sin embargo, la diabetes parece asociarse con un desplazamiento hacia valores mayores de cTnT.

SEXO

data$Sexe <- factor(data$Sexe,
                    levels = c(1, 2),
                    labels = c("Masculino", "Femenino"))

p3 <- ggplot(data, aes(x = Sexe, y = cTnT, fill = Sexe)) +
  geom_rain(
    alpha = 0.6,
    color = "black"
  ) +
  geom_hline(yintercept = 14, linetype = "dashed") +
  coord_flip() +
  theme_minimal() +
  scale_fill_manual(values = c("Masculino" = "lightblue", 
                               "Femenino" = "purple")) +
  labs(
    title = "Distribución de cTnT por Sexo",
    x = "Sexo",
    y = "cTnT",
    fill = "Sexo"
  )

p3

El raincloud plot sugiere que la mayoría de las mediciones de cTnT en ambos sexos se encuentra por debajo del umbral clínico de 14. Sin embargo, el grupo femenino presenta una distribución desplazada hacia valores más altos, con mayor dispersión y más valores extremos que el grupo masculino.A pesar de que el grupo femenino parece desplazado hacia valores mayores, existe superposición importante entre ambos sexos en los rangos bajos e intermedios. Eso significa que el sexo no separa completamente las distribuciones de cTnT.

data$cTnT_log <- log(data$cTnT + 0.01)

pairs(data[, c("cTnT_log", "age", "creat")],
      main = "Relaciones (log cTnT, age, creat)")

pairs(data[, c("cTnT_log", "age", "creat")],
      main = "Pairs plot",
      pch = 19,
      col = rgb(0, 0, 0, 0.4))

El pairs plot sugiere que cTnT se relaciona positivamente tanto con la edad como con la creatinina, siendo la asociación visualmente más clara con esta última. La edad también muestra una relación positiva con cTnT, aunque más dispersa.

EJERCICIO 2

#A.

data <- data %>%
  mutate(
    Sexe = as.factor(Sexe),
    Group = as.factor(Group),
    hx_diabete = as.factor(hx_diabete)
  )

#modelo lineal

i. Realice una regresión lineal estándar con cTnT como variable dependiente y age, Sexe, Group, hx_diabete y creat como variables independientes. Almacene el modelo

modelo_ctnt <- lm(cTnT ~ age + Sexe + Group + hx_diabete + creat, data = data)
summary(modelo_ctnt)

## 
## Call:
## lm(formula = cTnT ~ age + Sexe + Group + hx_diabete + creat, 
##     data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.755 -2.888 -0.893  1.855 46.245 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -12.026988   1.235859  -9.732  < 2e-16 ***
## age                      0.064317   0.019878   3.236  0.00124 ** 
## SexeFemenino             1.612530   0.375731   4.292 1.89e-05 ***
## GroupMultiple Events     3.482445   0.419374   8.304 2.24e-16 ***
## GroupSingle MI          -0.223000   0.401394  -0.556  0.57859    
## GroupStable Angina       0.181713   0.419880   0.433  0.66524    
## hx_diabeteCon diabetes   4.940508   0.392056  12.602  < 2e-16 ***
## creat                    0.118525   0.007607  15.582  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.337 on 1483 degrees of freedom
## Multiple R-squared:  0.4092, Adjusted R-squared:  0.4064 
## F-statistic: 146.7 on 7 and 1483 DF,  p-value: < 2.2e-16

table(data$Sexe, useNA = "ifany")

## 
## Masculino  Femenino 
##       254      1237

table(data$Group, useNA = "ifany")

## 
##         Control Multiple Events       Single MI   Stable Angina 
##             375             370             375             371

table(data$hx_diabete, useNA = "ifany")

## 
## Sin diabetes Con diabetes 
##         1210          281

ii. Use tidy para mostrar los resultados del modelo

resultadosmodelo_tidy <- tidy(modelo_ctnt)

resultadosmodelo_tidy

## # A tibble: 8 × 5
##   term                   estimate std.error statistic  p.value
##   <chr>                     <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)            -12.0      1.24       -9.73  9.72e-22
## 2 age                      0.0643   0.0199      3.24  1.24e- 3
## 3 SexeFemenino             1.61     0.376       4.29  1.89e- 5
## 4 GroupMultiple Events     3.48     0.419       8.30  2.24e-16
## 5 GroupSingle MI          -0.223    0.401      -0.556 5.79e- 1
## 6 GroupStable Angina       0.182    0.420       0.433 6.65e- 1
## 7 hx_diabeteCon diabetes   4.94     0.392      12.6   1.15e-34
## 8 creat                    0.119    0.00761    15.6   8.23e-51

iii. Evalúe el modelo utilizando plot

par(mfrow = c(1, 1))
plot(modelo_ctnt, which = 1)

plot(modelo_ctnt, which = 2)

plot(modelo_ctnt, which = 3)

plot(modelo_ctnt, which = 5)

Se evaluó el ajuste del modelo lineal mediante los gráficos diagnósticos, que sugiere falta de linealidad. El gráfico Q-Q mostró una desviación marcada respecto a la línea teórica, especialmente en la cola superior, El gráfico Scale-Location evidenció un aumento de la dispersión de los residuos a medida que aumentan los valores ajustados. Finalmente, el gráfico de residuos versus leverage mostró algunas observaciones potencialmente influyentes. En conjunto, estos resultados indican que el modelo lineal estándar no cumple adecuadamente los supuestos del modelo y su ajuste no es satisfactorio.

iv. Cree un modelo en el que las variables continuas estén transformadas con logaritmo (Corrección de modelo)

modelo_ctnt_log <- lm(cTnT_log ~ age + Sexe + Group + hx_diabete + creat, data = data)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(modelo_ctnt_log)

par(mfrow = c(1, 1))

Modelo corregido presentado por tidy

tidy(modelo_ctnt_log, conf.int = TRUE)

## # A tibble: 8 × 7
##   term                  estimate std.error statistic  p.value conf.low conf.high
##   <chr>                    <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)           -2.45      0.236    -10.4    2.04e-24 -2.92      -1.99  
## 2 age                    0.0247    0.00380    6.50   1.07e-10  0.0173     0.0322
## 3 SexeFemenino           0.485     0.0719     6.75   2.11e-11  0.344      0.626 
## 4 GroupMultiple Events   0.806     0.0802    10.0    5.41e-23  0.648      0.963 
## 5 GroupSingle MI         0.00245   0.0768     0.0319 9.75e- 1 -0.148      0.153 
## 6 GroupStable Angina     0.148     0.0803     1.84   6.53e- 2 -0.00939    0.306 
## 7 hx_diabeteCon diabet…  0.641     0.0750     8.54   3.25e-17  0.494      0.788 
## 8 creat                  0.0133    0.00146    9.15   1.85e-19  0.0105     0.0162

glance(modelo_ctnt_log)

## # 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.333         0.330  1.02      106. 7.28e-126     7 -2143. 4304. 4352.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Se decidió transformar la variable dependiente cTnT usando logaritmo, debido a que los gráficos diagnósticos del modelo lineal estándar mostraron incumplimiento de supuestos, especialmente falta de normalidad de los residuos y heterocedasticidad. La transformación logarítmica se realizó para mejorar el ajuste del modelo y estabilizar la varianza.

En este modelo ajustado no hubo un patrón tan marcado como en el modelo original, lo que sugiere una mejora en la linealidad. El gráfico Scale-Location mostró una dispersión más estable, indicando mejor homocedasticidad. Sin embargo, en el gráfico Q-Q los residuos aún se apartan de la línea teórica, por lo que la normalidad no se cumple completamente. El gráfico de residuos versus leverage no mostró observaciones extremadamente influyentes, aunque sí algunos casos a revisar. En conjunto, la transformación logarítmica mejoró el ajuste del modelo

v. Aceptando este segundo modelo, ¿cuál es el valor predicho de log cTnT para un hombre de 50 años, sin diabetes, creatinina de 70 y que pertenece al grupo de múltiples eventos? ¿Cuál sería el valor predicho para un paciente idéntico con diabetes? ¿Qué observa sobre la diferencia?

nuevo_sin_dm <- data.frame(
  age = 50,
  Sexe = factor("Masculino", levels = levels(data$Sexe)),
  Group = factor("Multiple Events", levels = levels(data$Group)),
  hx_diabete = factor("Sin diabetes", levels = levels(data$hx_diabete)),
  creat = 70
)

nuevo_con_dm <- data.frame(
  age = 50,
  Sexe = factor("Masculino", levels = levels(data$Sexe)),
  Group = factor("Multiple Events", levels = levels(data$Group)),
  hx_diabete = factor("Con diabetes", levels = levels(data$hx_diabete)),
  creat = 70
)

predict(modelo_ctnt_log, newdata = nuevo_sin_dm)

##         1 
## 0.5198216

predict(modelo_ctnt_log, newdata = nuevo_con_dm)

##        1 
## 1.160531

Para un hombre de 50 años, sin diabetes, con creatinina de 70 y perteneciente al grupo de múltiples eventos sin diabetes mellitus, el valor predicho de log(cTnT) es aproximadamente 0.522. Para un paciente idéntico pero con diabetes, el valor predicho es aproximadamente 1.163. La diferencia entre ambos es 0.641, que corresponde exactamente al coeficiente estimado para diabetes en el modelo. Al realizar el exponencial de esta diferencia, volvemos a la escala original de la variable que da un resultado de 1.9 veces mayor que el paciente sin diabetes.

vi. Considerando los resultados de la regresión como medidas de asociación, ¿cómo interpretaría los coeficientes? (Pista: ¿son medidas absolutas o relativas, es decir, Razones o Diferencias?)

Debido a que la variable dependiente fue transformada usando logaritmo, los coeficientes del modelo no deben interpretarse como diferencias absolutas en cTnT, sino como asociaciones relativas.

¿Cómo interpretaría el coeficiente para hx_diabete y los IC 95%?

El coeficiente de hx_diabete fue 0.641. Esto quiere decir que, comparando pacientes iguales en edad, sexo, grupo y creatinina, los que tienen diabetes tienden a tener valores más altos de cTnT que los que no tienen diabetes. Como el modelo se hizo con el logaritmo de cTnT, ese número primero se interpreta en escala logarítmica. Al pasarlo a la escala original con la exponencial, se obtiene aproximadamente 1.9. Esto significa que los pacientes con diabetes tienen un valor esperado de cTnT cerca de 1.9 veces mayor que los pacientes sin diabetes.

El intervalo de confianza del 95% fue de 0.494 a 0.788 en la escala logarítmica. Al pasarlo a la escala original, corresponde aproximadamente a 1.64 a 2.20. Esto indica que el verdadero efecto de la diabetes probablemente está entre 1.64 y 2.20 veces el valor de los pacientes sin diabetes. Como este intervalo no incluye 1, la asociación parece ser estadísticamente significativa.

EJERCICIO 3

i Utilizando la información disponible en el conjunto de datos, las variables utilizadas en el modelo (ejercicio

#3) y conocimiento general, elabore un DAG que represente la asociación entre hx_diabete y cTnT ###

install.packages("dagitty")

## Warning: package 'dagitty' is in use and will not be installed

library(dagitty)

dag_modelo <- dagitty('dag {
Antecedentes_CV [pos="-0.634,1.546"]
Creatinina      [pos="0.045,1.227"]
Diabetes        [exposure,pos="-1.461,1.608"]
Edad            [pos="-0.493,1.369"]
Sexo            [pos="-0.766,1.196"]
cTnT            [outcome,pos="0.491,1.616"]

Antecedentes_CV -> Creatinina
Antecedentes_CV -> Diabetes
Antecedentes_CV -> cTnT
Creatinina -> cTnT
Diabetes -> cTnT
Edad -> Antecedentes_CV
Edad -> Creatinina
Edad -> Diabetes
Edad -> cTnT
Sexo -> Antecedentes_CV
Sexo -> Creatinina
Sexo -> Diabetes
Sexo -> cTnT
}')

plot(dag_modelo)

ii. Indique si considera que existe confusión potencial y cuáles serían los posibles factores de confusión de esta relación. Pista: ¿Tenemos todos los factores de confusión medidos y no medidos? [1 punto

Sí, considero que existe confusión potencial en la relación entre hx_diabete y cTnT. Según el DAG propuesto, los principales factores de confusión medidos serían Edad, Sexo y Antecedentes CV , ya que estas variables se relacionan tanto con la exposición como con el desenlace. En cambio, Creatinina, no actuaría como confusor, sino como mediador porque no es una causa de diabetes. Como se muestra en la figura derecha de colores creada en la plataforma Dagitty, se ajusta por sexo, edad y por el grupo de antecedentes Cardiovasculares. Adicionalmente, considero que otras variables posiblemente confusoras son la hipertensión, tabaquismo, obesidad que no son incluídas en el análisis.

iii. Si considera que hay confusión, indique el conjunto mínimo de variables que se deben controlar/ajustar en este modelo para evitar confusión. Nota: Su respuesta debe ser consistente con el DAG que dibujó en el punto i). [2 puntos]

Como lo nombrado en el anterior punto, considero que es pertinente ajustar por edad, sexo, y antecedentes cardiovasculares. Esto significa incluir esas variables en el mismo modelo de regresión lineal, esto permitirá estimar el efecto de la exposición controlando por otras variables “confusoras”

iv. Indique si la relación entre hx_diabete y cTnT es una asociación “causal” o no. Explique por qué lo es o por qué no.

REGRESIÓN LOGÍSTICA

Cargar Base de datos Diabetes

setwd("~/MAESTRIA EPIDEMIOLOGIA ICESI/recoleccion_d/taller_s6")

library(readr)
diabetes <- read_csv("diabetes.csv")

## Rows: 392 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (9): pregnancies, glucose, bloodpressure, skinthickness, insulin, bmi, d...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

View(diabetes)

i. ¿Cuántas observaciones hay en el conjunto de datos?

str(diabetes)

## spc_tbl_ [392 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ pregnancies  : num [1:392] -0.7165 -1.0279 -0.0937 -0.4051 -0.7165 ...
##  $ glucose      : num [1:392] -1.09 0.466 -1.446 2.41 2.151 ...
##  $ bloodpressure: num [1:392] -0.3732 -2.4538 -1.6536 -0.0531 -0.8533 ...
##  $ skinthickness: num [1:392] -0.584 0.557 0.271 1.508 -0.584 ...
##  $ insulin      : num [1:392] -0.522 0.101 -0.573 3.256 5.806 ...
##  $ bmi          : num [1:392] -0.71 1.425 -0.297 -0.368 -0.425 ...
##  $ dpf          : num [1:392] -1.031 5.109 -0.796 -1.057 -0.362 ...
##  $ age          : num [1:392] -0.967 0.209 -0.477 2.17 2.758 ...
##  $ outcome      : num [1:392] 0 1 1 1 1 1 1 0 1 0 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   pregnancies = col_double(),
##   ..   glucose = col_double(),
##   ..   bloodpressure = col_double(),
##   ..   skinthickness = col_double(),
##   ..   insulin = col_double(),
##   ..   bmi = col_double(),
##   ..   dpf = col_double(),
##   ..   age = col_double(),
##   ..   outcome = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

ii. ¿Cuántos parámetros hay en el conjunto de datos?

Hay 9 parámetros.

iii. ¿Cuántos eventos de outcome ocurrieron?

table(diabetes$outcome)

## 
##   0   1 
## 262 130

Hubo 130 outcomes.

iv. Graficar el histograma de las variables estandarizadas pregnancy y age

library(ggplot2)

ggplot(diabetes, aes(x = pregnancies)) +
  geom_histogram(bins = 30, fill = "pink", color = "white") +
  labs(title = "Histograma de Embarazos",
       x = "Embarazos estandarizada",
       y = "Frecuencia")

ggplot(diabetes, aes(x = age)) +
  geom_histogram(bins = 30, fill = "blue", color = "white") +
  labs(title = "Histograma de edad",
       x = "Edad estandarizada",
       y = "Frecuencia")

v. Graficar un diagrama de caja (boxplot) de bmi versus outcome y otro para age

library(ggplot2)
library(dplyr)

# Convertir outcome a factor
diabetes <- diabetes %>%
  mutate(
    outcome_f = factor(outcome,
                       levels = c(0, 1),
                       labels = c("No diabetes", "Diabetes"))
  )

# Boxplot de BMI vs desenlace
ggplot(diabetes, aes(x = outcome_f, y = bmi)) +
  geom_boxplot(fill = "orange", color = "black", width = 0.6) +
  labs(
    title = "Boxplot de BMI según desenlace",
    x = "Desenlace",
    y = "BMI estandarizado"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(hjust = 0.5),
    panel.grid.minor = element_blank()
  )

# Boxplot de edad vs desenlace
ggplot(diabetes, aes(x = outcome_f, y = age)) +
  geom_boxplot(fill = "blue", color = "black", width = 0.6) +
  labs(
    title = "Boxplot de edad según desenlace",
    x = "Desenlace",
    y = "Edad estandarizada"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(hjust = 0.5),
    panel.grid.minor = element_blank()
  )

*El boxplot de BMI según el desenlace sugiere que las personas con diabetes tienden a tener valores de BMI estandarizado más altos que las personas sin diabetes. Esto se observa porque la mediana del grupo con diabetes está por encima de la mediana del grupo sin diabetes. Además, el grupo con diabetes parece presentar algunos valores extremos altos.

*El boxplot de edad según el desenlace sugiere que las personas con diabetes tienden a tener edades más altas que las personas sin diabetes. Esto se observa porque la mediana del grupo con diabetes es mayor que la del grupo sin diabetes. Al contrario que en pacientes con tendencia a BMI más altos tienden a presnetar diabetes, en este caso el grupo sin diabetes presenta valores extremos altos de edad.

*En general la mediana del grupo con diabetes es mayor en ambos contextos, lo que sugiere una posible asociación positiva de estas variables con el desenlace.

vi. Crear un gráfico de correlación usando corrplot del paquete del mismo nombre. Coloca outcome en la primera fila

library(corrplot)

## Warning: package 'corrplot' was built under R version 4.5.3

## corrplot 0.95 loaded

library(dplyr)

# Si outcome quedó como factor en pasos anteriores, volverlo numérico
diabetes <- diabetes %>%
  mutate(outcome = as.numeric(as.character(outcome)))

# Reordenar columnas para que outcome quede primero
diabetes_corr <- diabetes %>%
  select(outcome, pregnancies, glucose, bloodpressure,
         skinthickness, insulin, bmi, dpf, age)

# Matriz de correlación
mat_cor <- cor(diabetes_corr, use = "complete.obs")

# Gráfico de correlación
corrplot(mat_cor,
         method = "color",
         type = "upper",
         addCoef.col = "black",
         tl.col = "black",
         tl.srt = 45)

vii. Usar ggpairs del paquete GGally para mostrar un Exploratory Data Analysis más completo que examine la correlación entre variables y con el outcome

library(GGally)

## Warning: package 'GGally' was built under R version 4.5.3

library(ggplot2)
library(dplyr)

# Preparar base
diabetes_eda <- diabetes %>%
  mutate(
    outcome_f = factor(outcome, levels = c(0, 1),
                       labels = c("No diabetes", "Diabetes"))
  ) %>%
  select(outcome_f, glucose, bmi, age, pregnancies, insulin)

# Panel de correlación limpio
panel_cor <- function(data, mapping, ...) {
  x <- GGally::eval_data_col(data, mapping$x)
  y <- GGally::eval_data_col(data, mapping$y)
  r <- cor(x, y, use = "complete.obs")
  
  ggally_text(
    label = paste0("r = ", round(r, 2)),
    xP = 0.5, yP = 0.5,
    size = 5
  ) +
    theme_void()
}

# ggpairs más limpio
p_eda <- ggpairs(
  diabetes_eda,
  mapping = aes(color = outcome_f, fill = outcome_f),
  upper = list(
    continuous = panel_cor
  ),
  lower = list(
    continuous = wrap("points", alpha = 0.35, size = 0.7)
  ),
  diag = list(
    continuous = wrap("densityDiag", alpha = 0.35),
    discrete = wrap("barDiag", alpha = 0.8)
  )
) +
  scale_color_manual(values = c("No diabetes" = "grey60",
                                "Diabetes" = "pink3")) +
  scale_fill_manual(values = c("No diabetes" = "grey80",
                               "Diabetes" = "pink")) +
  theme_bw() +
  theme(
    strip.background = element_rect(fill = "white"),
    panel.grid = element_blank(),
    legend.position = "bottom",
    legend.title = element_blank(),
    text = element_text(size = 10)
  )

p_eda

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

El gráfico ggpairs muestra un análisis exploratorio más completo de las variables. En la comparación por desenlace, glucosa y edad son las variables que presentan mayor separación entre pacientes con diabetes y sin diabetes, seguidas en menor medida por embarazos, BMI e insulina. En la diagonal se observa que varias variables tienen distribuciones asimétricas, especialmente embarazos, insulina y edad. En cuanto a la correlación entre predictores, las relaciones más altas se observan entre edad y embarazos (r = 0.68) y entre glucosa e insulina (r = 0.58). En conjunto, el gráfico sugiere que glucose y age podrían ser variables importantes para explicar el desenlace

EJERCICIO 2

a) Generar dos nuevas variables

library(dplyr)

diabetes <- diabetes %>%
  mutate(
    hgl = if_else(glucose >= 0.6, 1, 0),
    hbp = if_else(bloodpressure >= 0.5, 1, 0)
  )

b) Regresión logística estándar

library(dplyr)
library(ggplot2)

# 1. Ajustar modelo logístico
modelo_bmi <- glm(outcome ~ bmi,
                  data = diabetes,
                  family = binomial(link = "logit"))

summary(modelo_bmi)

## 
## Call:
## glm(formula = outcome ~ bmi, family = binomial(link = "logit"), 
##     data = diabetes)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.7498     0.1127  -6.651 2.92e-11 ***
## bmi           0.6067     0.1198   5.062 4.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 498.10  on 391  degrees of freedom
## Residual deviance: 469.03  on 390  degrees of freedom
## AIC: 473.03
## 
## Number of Fisher Scoring iterations: 4

# 2. Generar log-odds y probabilidades predichas
diabetes <- diabetes %>%
  mutate(
    log_odds_pred = predict(modelo_bmi, type = "link"),
    prob_pred = predict(modelo_bmi, type = "response")
  )

# Revisar primeras filas
head(diabetes[, c("bmi", "outcome", "log_odds_pred", "prob_pred")])

## # A tibble: 6 × 4
##      bmi outcome log_odds_pred prob_pred
##    <dbl>   <dbl>         <dbl>     <dbl>
## 1 -0.710       0        -1.18      0.235
## 2  1.42        1         0.115     0.529
## 3 -0.297       1        -0.930     0.283
## 4 -0.368       1        -0.973     0.274
## 5 -0.425       1        -1.01      0.267
## 6 -1.04        1        -1.38      0.201

# 3. Ordenar para que las líneas se vean bien en los gráficos
diabetes_plot <- diabetes %>%
  arrange(bmi)

# 4. Gráfico de log-odds predichos vs bmi
ggplot(diabetes_plot, aes(x = bmi, y = log_odds_pred)) +
  geom_point(alpha = 0.5, color = "orange") +
  geom_line(color = "black", linewidth = 1) +
  theme_bw() +
  labs(
    title = "Log-odds predichos de outcome = 1 según BMI",
    x = "BMI",
    y = "Log-odds predichos"
  )

# 5. Gráfico de probabilidad predicha vs bmi
ggplot(diabetes_plot, aes(x = bmi, y = prob_pred)) +
  geom_point(alpha = 0.5, color = "steelblue") +
  geom_line(color = "black", linewidth = 1) +
  theme_bw() +
  labs(
    title = "Probabilidad predicha de outcome = 1 según BMI",
    x = "BMI",
    y = "Probabilidad predicha"
  )

La línea recta es la de los log-odds (naranja) y la curva es la de la probabilidad (Azul) . Esto quiere decir que en regresión logística el efecto es constante en la escala de log-odds, pero no en la escala de probabilidad de riesgo. Entonces, la diferencia de riesgo no es la misma en todos los valores de BMI, mientras que la diferencia en lnOR sí sigue una relación más constante lineal.

c) Modelo Ajustado Ajuste un modelo de regresión logística frecuentista estándar (incondicional) para el resultado como función de hgl y hbp. Grafique la probabilidad predicha de outcome según los niveles de hgl y bmi de acuerdo con este modelo. Proporcione comentarios sobre las diferencias relativas y absolutas esperadas entre hgl = 0 y hgl = 1.

modelo_c <- glm(outcome ~ bmi + hgl + hbp,
                data = diabetes,
                family = binomial(link = "logit"))

library(dplyr)
library(ggplot2)

# Modelo ajustado incluyendo bmi
modelo_c <- glm(outcome ~ bmi + hgl + hbp,
                data = diabetes,
                family = binomial(link = "logit"))

summary(modelo_c)

## 
## Call:
## glm(formula = outcome ~ bmi + hgl + hbp, family = binomial(link = "logit"), 
##     data = diabetes)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.4391     0.1707  -8.432  < 2e-16 ***
## bmi           0.5043     0.1354   3.724 0.000196 ***
## hgl           1.9456     0.2618   7.431 1.08e-13 ***
## hbp           0.2858     0.2705   1.056 0.290765    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 498.10  on 391  degrees of freedom
## Residual deviance: 405.91  on 388  degrees of freedom
## AIC: 413.91
## 
## Number of Fisher Scoring iterations: 4

# Crear datos para predicción
pred_data <- expand.grid(
  bmi = seq(min(diabetes$bmi, na.rm = TRUE),
            max(diabetes$bmi, na.rm = TRUE),
            length.out = 200),
  hgl = c(0, 1),
  hbp = 0   # fijamos hbp en 0 para hacer un gráfico simple como tu ejemplo
)

# Probabilidad predicha
pred_data$prob_pred <- predict(modelo_c,
                               newdata = pred_data,
                               type = "response")

# Etiquetas
pred_data$hgl <- factor(pred_data$hgl,
                        levels = c(0, 1),
                        labels = c("0", "1"))

# Gráfico
ggplot(pred_data, aes(x = bmi, y = prob_pred, color = hgl)) +
  geom_line(linewidth = 1.2) +
  theme_minimal() +
  labs(
    title = "Predicted probability by BMI and hgl",
    x = "BMI",
    y = "Predicted probability of outcome",
    color = "hgl"
  )

pred_data <- expand.grid(
  bmi = seq(min(diabetes$bmi, na.rm = TRUE),
            max(diabetes$bmi, na.rm = TRUE),
            length.out = 200),
  hgl = c(0, 1),
  hbp = c(0, 1)
)

pred_data$prob_pred <- predict(modelo_c,
                               newdata = pred_data,
                               type = "response")

pred_data <- pred_data %>%
  mutate(
    hgl = factor(hgl, levels = c(0, 1), labels = c("0", "1")),
    hbp = factor(hbp, levels = c(0, 1), labels = c("0", "1"))
  )

ggplot(pred_data, aes(x = bmi, y = prob_pred, color = hgl)) +
  geom_line(linewidth = 1.2) +
  facet_wrap(~ hbp) +
  theme_minimal() +
  labs(
    title = "Predicted probability by BMI, hgl and hbp",
    x = "BMI",
    y = "Predicted probability of outcome",
    color = "hgl"
  )

En el gráfico se observa que, a medida que aumenta el BMI, también aumenta la probabilidad predicha del desenlace, tanto tener glucosa alta o glucosa alta. Además, en ambos niveles de presión alta, la curva de glucosa alta se encuentra por encima de la deglucosa baja, lo que sugiere que tener glucosa alta se asocia con una mayor probabilidad del desenlace. También se ve que, cuando hay presión alta, las curvas están un poco más arriba que cuando no la hay presión alta, lo que indica que la presión alta también podría aumentar la probabilidad del desenlace.

d) Riesgo predicho bajo exposición contrafactual: Genera una nueva variable en el conjunto de datos que sea el riesgo predicho de outcome si cada sujeto tiene glucosa BAJA (hgl). Para las personas no expuestas, esto es simplemente su probabilidad predicha. Para las personas expuestas, se aplica la fórmula expit a su combinación lineal específica con la variable hgl ajustada a 0 en lugar de 1.

library(dplyr)

# Logit observado
diabetes <- diabetes %>%
  mutate(
    lp_obs = predict(modelo_c, type = "link"),
    riesgo_predicho = predict(modelo_c, type = "response")
  )

# Coeficiente de hgl
beta_hgl <- coef(modelo_c)["hgl"]

# Riesgo contrafactual si hgl = 0
diabetes <- diabetes %>%
  mutate(
    lp_cf_hgl0 = if_else(hgl == 0,
                         lp_obs,
                         lp_obs - beta_hgl),
    riesgo_cf_hgl0 = plogis(lp_cf_hgl0)
  )

head(diabetes[, c("hgl", "lp_obs", "lp_cf_hgl0", "riesgo_predicho", "riesgo_cf_hgl0")])

## # A tibble: 6 × 5
##     hgl  lp_obs lp_cf_hgl0 riesgo_predicho riesgo_cf_hgl0
##   <dbl>   <dbl>      <dbl>           <dbl>          <dbl>
## 1     0 -1.80       -1.80            0.142          0.142
## 2     0 -0.720      -0.720           0.327          0.327
## 3     0 -1.59       -1.59            0.170          0.170
## 4     1  0.321      -1.62            0.580          0.165
## 5     1  0.292      -1.65            0.573          0.161
## 6     1 -0.0164     -1.96            0.496          0.123

A continuación, genera una nueva variable en el conjunto de datos que sea el riesgo predicho de outcome si cada sujeto tiene glucosa ALTA. Para las personas expuestas, esto es simplemente su probabilidad predicha. Para las personas no expuestas, se aplica la fórmula expit a su combinación lineal con hgl ajustada a 1 en lugar de 0. Ahora tendrás una columna con el riesgo bajo glucosa alta y otra columna con el riesgo bajo glucosa baja para cada persona. Crea dos nuevas variables:

# Predictor lineal observado
diabetes <- diabetes %>%
  mutate(
    lp_obs = predict(modelo_c, type = "link"),
    riesgo_predicho = predict(modelo_c, type = "response")
  )

# Coeficiente de hgl
beta_hgl <- coef(modelo_c)["hgl"]

# Riesgo contrafactual si hgl = 1
diabetes <- diabetes %>%
  mutate(
    lp_cf_hgl1 = if_else(hgl == 1,
                         lp_obs,
                         lp_obs + beta_hgl),
    riesgo_cf_hgl1 = plogis(lp_cf_hgl1)
  )

library(dplyr)

diabetes <- diabetes %>%
  mutate(
    RD = riesgo_cf_hgl1 - riesgo_cf_hgl0,
    RR = riesgo_cf_hgl1 / riesgo_cf_hgl0
  )

library(dplyr)

diabetes <- diabetes %>%
  mutate(
    riesgo_cf_hgl0 = predict(modelo_c,
                             newdata = mutate(., hgl = 0),
                             type = "response"),
    riesgo_cf_hgl1 = predict(modelo_c,
                             newdata = mutate(., hgl = 1),
                             type = "response"),
    RD = riesgo_cf_hgl1 - riesgo_cf_hgl0,
    RR = riesgo_cf_hgl1 / riesgo_cf_hgl0
  )

library(ggplot2)

ggplot(diabetes, aes(x = bmi, y = RD)) +
  geom_point(alpha = 0.5, color = "orange") +
  geom_smooth(se = FALSE, color = "black") +
  theme_bw() +
  labs(
    title = "Diferencia de riesgos (RD) según BMI",
    x = "BMI",
    y = "RD"
  )

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(diabetes, aes(x = bmi, y = RR)) +
  geom_point(alpha = 0.5, color = "steelblue") +
  geom_smooth(se = FALSE, color = "black") +
  theme_bw() +
  labs(
    title = "Razón de riesgos (RR) según BMI",
    x = "BMI",
    y = "RR"
  )

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

diabetes <- diabetes %>%
  mutate(
    hbp_f = factor(hbp, levels = c(0, 1),
                   labels = c("Sin presión alta", "Con presión alta"))
  )

ggplot(diabetes, aes(x = bmi, y = RD, color = hbp_f)) +
  geom_point(alpha = 0.5) +
  geom_smooth(se = FALSE) +
  theme_bw() +
  labs(
    title = "RD según BMI y presión alta",
    x = "BMI",
    y = "RD",
    color = "hbp"
  )

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(diabetes, aes(x = bmi, y = RR, color = hbp_f)) +
  geom_point(alpha = 0.5) +
  geom_smooth(se = FALSE) +
  theme_bw() +
  labs(
    title = "RR según BMI y presión alta",
    x = "BMI",
    y = "RR",
    color = "hbp"
  )

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

diabetes %>%
  select(bmi, hbp, RD, RR) %>%
  arrange(RD)

## # A tibble: 392 × 4
##      bmi   hbp    RD    RR
##    <dbl> <dbl> <dbl> <dbl>
##  1  4.84     1 0.178  1.23
##  2  3.74     1 0.260  1.38
##  3  3.45     1 0.284  1.44
##  4 -2.12     0 0.288  4.82
##  5 -1.93     0 0.303  4.69
##  6 -1.93     0 0.303  4.69
##  7 -1.92     0 0.304  4.68
##  8 -1.85     0 0.310  4.63
##  9 -1.81     0 0.313  4.60
## 10 -1.81     0 0.313  4.60
## # ℹ 382 more rows

diabetes %>%
  select(bmi, hbp, RD, RR) %>%
  arrange(desc(RD))

## # A tibble: 392 × 4
##      bmi   hbp    RD    RR
##    <dbl> <dbl> <dbl> <dbl>
##  1 0.372     1 0.451  2.64
##  2 0.343     1 0.451  2.65
##  3 0.941     0 0.451  2.64
##  4 0.898     0 0.451  2.66
##  5 0.898     0 0.451  2.66
##  6 0.898     0 0.451  2.66
##  7 0.386     1 0.451  2.63
##  8 0.400     1 0.451  2.62
##  9 0.970     0 0.451  2.62
## 10 0.301     1 0.451  2.68
## # ℹ 382 more rows

diabetes %>%
  select(bmi, hbp, RD, RR) %>%
  arrange(RR)

## # A tibble: 392 × 4
##      bmi   hbp    RD    RR
##    <dbl> <dbl> <dbl> <dbl>
##  1  4.84     1 0.178  1.23
##  2  3.74     1 0.260  1.38
##  3  3.45     1 0.284  1.44
##  4  2.73     1 0.341  1.61
##  5  3.12     0 0.356  1.67
##  6  2.36     1 0.369  1.72
##  7  2.86     0 0.374  1.75
##  8  1.95     1 0.397  1.87
##  9  1.94     1 0.398  1.87
## 10  1.91     1 0.400  1.88
## # ℹ 382 more rows

diabetes %>%
  select(bmi, hbp, RD, RR) %>%
  arrange(desc(RR))

## # A tibble: 392 × 4
##      bmi   hbp    RD    RR
##    <dbl> <dbl> <dbl> <dbl>
##  1 -2.12     0 0.288  4.82
##  2 -1.93     0 0.303  4.69
##  3 -1.93     0 0.303  4.69
##  4 -1.92     0 0.304  4.68
##  5 -1.85     0 0.310  4.63
##  6 -1.81     0 0.313  4.60
##  7 -1.81     0 0.313  4.60
##  8 -1.75     0 0.318  4.56
##  9 -1.75     0 0.318  4.56
## 10 -1.71     0 0.321  4.52
## # ℹ 382 more rows

#¿Dónde es mayor y dónde es menor la RD? La diferencia de riesgo es mayor en los sujetos con valores intermedios de bmi, donde la diferencia absoluta entre el riesgo bajo glucosa alta y el riesgo bajo glucosa baja es más grande. En cambio, la diferencia de riesgo es menor en los sujetos con bmi muy bajo o muy alto, donde esa diferencia absoluta es más pequeña.

#¿Dónde es mayor y dónde es menor la RR? La razón de riesgos es mayor al inicio del rango de bmi y va disminuyendo a medida que el bmi aumenta. Esto pasa porque, cuando el riesgo inicial es bajo, el cambio relativo entre glucosa alta y glucosa baja se ve más grande.

EJERCICIO 3

Ajusta un modelo de regresión logística frecuentista estándar (incondicional) para outcome como función de las variables: pregnancies, hgl, hbp, skinthickness, insulin, bmi, dpf y age.

library(dplyr)
library(broom)

# Asegúrate de que outcome esté como 0/1 numérico
table(diabetes$outcome, useNA = "ifany")

## 
##   0   1 
## 262 130

# Ajustar el modelo de regresión logística
modelo_p3 <- glm(
  outcome ~ pregnancies + hgl + hbp + skinthickness + insulin + bmi + dpf + age,
  data = diabetes,
  family = binomial(link = "logit")
)

# Ver resumen del modelo
summary(modelo_p3)

## 
## Call:
## glm(formula = outcome ~ pregnancies + hgl + hbp + skinthickness + 
##     insulin + bmi + dpf + age, family = binomial(link = "logit"), 
##     data = diabetes)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.3420     0.1809  -7.418 1.19e-13 ***
## pregnancies     0.1998     0.1703   1.173  0.24082    
## hgl             1.5064     0.3094   4.869 1.12e-06 ***
## hbp             0.1240     0.2898   0.428  0.66860    
## skinthickness   0.1329     0.1724   0.771  0.44074    
## insulin         0.1651     0.1468   1.124  0.26083    
## bmi             0.4364     0.1803   2.420  0.01553 *  
## dpf             0.3717     0.1360   2.733  0.00628 ** 
## age             0.4199     0.1769   2.374  0.01762 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 498.10  on 391  degrees of freedom
## Residual deviance: 372.41  on 383  degrees of freedom
## AIC: 390.41
## 
## Number of Fisher Scoring iterations: 4

####1. Reporta los OR completamente ajustados para el outcome y la nueva variable binaria hgl, usando una interpretación causal de este estimado (es decir, en términos de una intervención hipotética sobre los niveles de glucosa (hgl)).####

exp(1.5064)

## [1] 4.510464

Este estimado sugiere que si pudiéramos intervenir y hacer que una persona pasara de glucosa baja a glucosa alta, manteniendo iguales las demás variables del modelo, las odds del desenlace serían aproximadamente 4.5 veces mayores.

####2. Estima los efectos marginales promedio: RR (riesgo relativo), RD (diferencia de riesgos) y OR (odds ratio) usando el modelo frecuentista del Ejercicio 3.1. #### Sugerencia: puedes usar el paquete StdGLM Compara los estimados obtenidos mediante estandarización marginal con aquellos obtenidos directamente del modelo de regresión logística del Ejercicio 3.1.

library(dplyr)
library(broom)

# Modelo del punto 3.1
table(diabetes$outcome, useNA = "ifany")

## 
##   0   1 
## 262 130

modelo_p3 <- glm(
  outcome ~ pregnancies + hgl + hbp + skinthickness + insulin + bmi + dpf + age,
  data = diabetes,
  family = binomial(link = "logit")
)

summary(modelo_p3)

## 
## Call:
## glm(formula = outcome ~ pregnancies + hgl + hbp + skinthickness + 
##     insulin + bmi + dpf + age, family = binomial(link = "logit"), 
##     data = diabetes)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.3420     0.1809  -7.418 1.19e-13 ***
## pregnancies     0.1998     0.1703   1.173  0.24082    
## hgl             1.5064     0.3094   4.869 1.12e-06 ***
## hbp             0.1240     0.2898   0.428  0.66860    
## skinthickness   0.1329     0.1724   0.771  0.44074    
## insulin         0.1651     0.1468   1.124  0.26083    
## bmi             0.4364     0.1803   2.420  0.01553 *  
## dpf             0.3717     0.1360   2.733  0.00628 ** 
## age             0.4199     0.1769   2.374  0.01762 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 498.10  on 391  degrees of freedom
## Residual deviance: 372.41  on 383  degrees of freedom
## AIC: 390.41
## 
## Number of Fisher Scoring iterations: 4

# Riesgo predicho si todos tuvieran hgl = 1
p_hgl1 <- predict(
  modelo_p3,
  newdata = mutate(diabetes, hgl = 1),
  type = "response"
)

# Riesgo predicho si todos tuvieran hgl = 0
p_hgl0 <- predict(
  modelo_p3,
  newdata = mutate(diabetes, hgl = 0),
  type = "response"
)

# Riesgos promedio estandarizados
risk_hgl1 <- mean(p_hgl1)
risk_hgl0 <- mean(p_hgl0)

# Efectos marginales promedio
RD_marginal <- risk_hgl1 - risk_hgl0
RR_marginal <- risk_hgl1 / risk_hgl0
OR_marginal <- (risk_hgl1 / (1 - risk_hgl1)) / (risk_hgl0 / (1 - risk_hgl0))

# Resultados
resultados_marginales <- tibble(
  riesgo_hgl1 = risk_hgl1,
  riesgo_hgl0 = risk_hgl0,
  RD = RD_marginal,
  RR = RR_marginal,
  OR = OR_marginal
)

resultados_marginales

## # A tibble: 1 × 5
##   riesgo_hgl1 riesgo_hgl0    RD    RR    OR
##         <dbl>       <dbl> <dbl> <dbl> <dbl>
## 1       0.536       0.252 0.284  2.13  3.43

# OR ajustado directo del modelo para hgl
OR_modelo <- exp(coef(modelo_p3)["hgl"])
OR_modelo

##      hgl 
## 4.510472

Mediante estandarización marginal, el riesgo promedio estimado de desenlace fue 0.536 bajo glucosa alta y 0.252 bajoglucosa baja. A partir de estos riesgos se obtuvo una diferencia de riesgos de 0.284, una razón de riesgos de 2.13 y un OR marginal de 3.43. En comparación, el OR ajustado directo del modelo logístico para hgl fue 4.51. Esto muestra que el OR marginal y el OR del modelo no son iguales: el OR del modelo es un efecto condicional, mientras que el obtenido por estandarización es un efecto marginal promedio.