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:

library(data.table)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between()     masks data.table::between()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ dplyr::first()       masks data.table::first()
## ✖ lubridate::hour()    masks data.table::hour()
## ✖ lubridate::isoweek() masks data.table::isoweek()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ dplyr::last()        masks data.table::last()
## ✖ lubridate::mday()    masks data.table::mday()
## ✖ lubridate::minute()  masks data.table::minute()
## ✖ lubridate::month()   masks data.table::month()
## ✖ lubridate::quarter() masks data.table::quarter()
## ✖ lubridate::second()  masks data.table::second()
## ✖ purrr::transpose()   masks data.table::transpose()
## ✖ lubridate::wday()    masks data.table::wday()
## ✖ lubridate::week()    masks data.table::week()
## ✖ lubridate::yday()    masks data.table::yday()
## ✖ lubridate::year()    masks data.table::year()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(broom)
library(vcd)
## Warning: package 'vcd' was built under R version 4.4.3
## Cargando paquete requerido: grid
library(sjPlot)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
library(epiR)
## Cargando paquete requerido: survival
## Package epiR 2.0.77 is loaded
## Type help(epi.about) for summary information
## Type browseVignettes(package = 'epiR') to learn how to use epiR for applied epidemiological analyses
library(gtsummary)
## Warning: package 'gtsummary' was built under R version 4.4.3
library(tableone)
library(ggpubr)
library(vcd)
library(vcdExtra)
## Warning: package 'vcdExtra' was built under R version 4.4.3
## Cargando paquete requerido: gnm
## Warning: package 'gnm' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'vcdExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     summarise
library(CalibrationCurves)
## Warning: package 'CalibrationCurves' was built under R version 4.4.3
## Cargando paquete requerido: rms
## Warning: package 'rms' was built under R version 4.4.3
## Cargando paquete requerido: Hmisc
## Warning: package 'Hmisc' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'Hmisc'
## 
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(pROC)
## Warning: package 'pROC' was built under R version 4.4.3
## Type 'citation("pROC")' for a citation.
## 
## Adjuntando el paquete: 'pROC'
## 
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
library(glmnet)
## Warning: package 'glmnet' was built under R version 4.4.3
## Cargando paquete requerido: Matrix
## 
## Adjuntando el paquete: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Loaded glmnet 4.1-9
library(caret)
## Warning: package 'caret' was built under R version 4.4.3
## Cargando paquete requerido: lattice
## 
## Adjuntando el paquete: 'lattice'
## 
## The following object is masked from 'package:gnm':
## 
##     barley
## 
## 
## Adjuntando el paquete: 'caret'
## 
## The following object is masked from 'package:survival':
## 
##     cluster
## 
## The following object is masked from 'package:purrr':
## 
##     lift
library(broom.helpers)
## Warning: package 'broom.helpers' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'broom.helpers'
## 
## The following objects are masked from 'package:gtsummary':
## 
##     all_categorical, all_continuous, all_contrasts, all_dichotomous,
##     all_interaction, all_intercepts
options(scipen = 999, digits = 3, encoding = 'UTF-8')

Se crea datos y se eliminan columna 7 a 10

#Se crea datos y se eliminan columna 7 a 10
library(readr)
fram <- read_csv("C:/Users/Administrador/Downloads/fram.csv")
## Rows: 1363 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (10): age, sbp, dbp, chol, sex, newchd, male, cig, packs, muerte
## 
## ℹ 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(fram)
datos <- fread("fram.csv")
datos$smoke <- ifelse(datos$cig >0,1,0)
datos <- datos[,-c(7:10)]
datos$smoke <- as.factor (datos$smoke)

Trasformo en factor aunque se va a usar datos para el análisis dado que es más eficiente este formato del dataframe

fram <- fram %>% mutate_at(vars(sex, newchd, male, muerte),~as.factor(.))

Predictivo: Modelo Parsimonioso Categorizo las variables de datos y hago la tabla 1 para evaluar las diferencias entre las variables clasificadas según newchd

datos [, c("sex", "smoke") := lapply(.SD, as.factor), .SDcols = c("sex", "smoke")]

catvars = c("sex", "smoke")
vars = c("age","sbp","dbp","chol","sex","smoke")

tabla_1 <- CreateTableOne(vars = vars, strata = "newchd", factorVars = catvars, data = datos)
kableone(tabla_1)
0 1 p test
n 1095 268
age (mean (SD)) 52.08 (4.75) 53.58 (4.76) <0.001
sbp (mean (SD)) 145.11 (26.28) 158.57 (31.14) <0.001
dbp (mean (SD)) 88.91 (13.53) 94.58 (15.44) <0.001
chol (mean (SD)) 232.99 (45.71) 241.30 (48.42) 0.008
sex = 1 (%) 479 (43.7) 164 (61.2) <0.001
smoke = 1 (%) 480 (43.9) 134 (50.0) 0.082

Análisis bivariado. Tabla

datos %>%
  tbl_uvregression(
    method = glm,
    y = newchd,
    method.args = list(family = binomial),
    exponentiate = TRUE,
    pvalue_fun = ~style_pvalue(.x, digits = 2)) %>% bold_p() %>% bold_labels()
Characteristic N OR 95% CI p-value
age 1,363 1.07 1.04, 1.10 <0.001
sbp 1,363 1.02 1.01, 1.02 <0.001
dbp 1,363 1.03 1.02, 1.04 <0.001
chol 1,363 1.00 1.00, 1.01 0.009
sex 1,363


    0

    1
2.03 1.55, 2.67 <0.001
smoke 1,362


    0

    1
1.28 0.98, 1.67 0.071
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

COLESTEROL

#colesterol. Lo paso a quintiles y hago el gráfico
datos<- datos %>% mutate(chol_q = ntile(chol, 5))

datos %>%
  dplyr::group_by(chol_q) %>%
  dplyr::summarise(tar = mean(newchd)) %>%
  ggplot(aes(x = chol_q, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", col = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## `geom_smooth()` using formula = 'y ~ x'

datos %>%
  dplyr::group_by(chol_q) %>%
  dplyr::summarise(tar = mean(newchd)) %>%
  ggplot(aes(x = chol_q, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "Cuartiles de colesterol", y = "Proporción de eventos CHD") +
  theme_minimal()

#MODELO COLESTEROL CUANTI
modelocol <- glm(newchd ~ chol, family = "binomial", data = fram)
t1<-tab_model(modelocol, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - col cuanti")
knitr::asis_output(t1$knitr)
Newchd - col cuanti
  newchd
Predictors Odds Ratios std. Error CI p
chol 1.00 0.00 1.00 – 1.01 0.009
Observations 1363
Deviance 1344.419
log-Likelihood -672.210
modelocol
## 
## Call:  glm(formula = newchd ~ chol, family = "binomial", data = fram)
## 
## Coefficients:
## (Intercept)         chol  
##    -2.30698      0.00379  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1361 Residual
## Null Deviance:       1350 
## Residual Deviance: 1340  AIC: 1350
tab_model(modelocol)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.10 0.05 – 0.20 <0.001
chol 1.00 1.00 – 1.01 0.009
Observations 1363
R2 Tjur 0.005
#MODELO COLESTEROL QUINTILOS
modelocolq2<- glm(newchd ~ chol_q, family = "binomial", data = datos)
t2<-tab_model(modelocolq2, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - col quantil")
knitr::asis_output(t2$knitr)
Newchd - col quantil
  newchd
Predictors Odds Ratios std. Error CI p
chol_q 1.09 0.05 0.99 – 1.20 0.070
Observations 1363
Deviance 1347.959
log-Likelihood -673.979
modelocolq2
## 
## Call:  glm(formula = newchd ~ chol_q, family = "binomial", data = datos)
## 
## Coefficients:
## (Intercept)       chol_q  
##     -1.6746       0.0875  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1361 Residual
## Null Deviance:       1350 
## Residual Deviance: 1350  AIC: 1350
tab_model(modelocolq2)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.19 0.13 – 0.26 <0.001
chol q 1.09 0.99 – 1.20 0.070
Observations 1363
R2 Tjur 0.002
#MODELO COLESTEROL Q1 COMO DUMMY. PREGUNTAR DIFERENCIAS
datos$chol_q1 <- datos$chol_q
datos$chol_q1 <- as.factor(datos$chol_q1)
modelocolq11<- glm(newchd ~ chol_q1, family = "binomial", data = datos)
tq1<-tab_model(modelocolq11, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - col quantil")
knitr::asis_output(tq1$knitr)
Newchd - col quantil
  newchd
Predictors Odds Ratios std. Error CI p
chol_q12 1.08 0.24 0.70 – 1.66 0.741
chol_q13 1.02 0.23 0.66 – 1.59 0.912
chol_q14 0.95 0.21 0.61 – 1.48 0.838
chol_q15 1.58 0.33 1.05 – 2.40 0.028
Observations 1363
Deviance 1343.405
log-Likelihood -671.702
modelocolq11
## 
## Call:  glm(formula = newchd ~ chol_q1, family = "binomial", data = datos)
## 
## Coefficients:
## (Intercept)     chol_q12     chol_q13     chol_q14     chol_q15  
##     -1.5198       0.0729       0.0247      -0.0461       0.4601  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1358 Residual
## Null Deviance:       1350 
## Residual Deviance: 1340  AIC: 1350
tab_model(modelocolq11)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.22 0.16 – 0.30 <0.001
chol q1 [2] 1.08 0.70 – 1.66 0.741
chol q1 [3] 1.02 0.66 – 1.59 0.912
chol q1 [4] 0.95 0.61 – 1.48 0.838
chol q1 [5] 1.58 1.05 – 2.40 0.028
Observations 1363
R2 Tjur 0.006
# COLESTEROL la deviance de los modelos cuanti y cuali son similares 1350 a 1340- VER POR QUE TIENEN EL MISMO GRADO DE LIBERTAD. 

EDAD

#edad. Paso edad a quintilos
datos<- datos %>% mutate(age_q = ntile(age, 5))


datos %>%
  dplyr::group_by(age_q) %>%
  dplyr::summarise(tar = mean(newchd)) %>%
  ggplot(aes(x = age_q, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", age = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## Warning in geom_smooth(method = "lm", age = "blue"): Ignoring unknown
## parameters: `age`
## `geom_smooth()` using formula = 'y ~ x'

datos %>%
  dplyr::group_by(age_q) %>%
  dplyr::summarise(tar = mean(newchd)) %>%
  ggplot(aes(x = age_q, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "Cuartiles de edad", y = "Proporción de eventos CHD") +
  theme_minimal()

#modelo edad continua
modeloedad <- glm(newchd ~ age, family = "binomial", data = fram)
t3<-tab_model(modeloedad, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - edad cuanti")
knitr::asis_output(t3$knitr)
Newchd - edad cuanti
  newchd
Predictors Odds Ratios std. Error CI p
age 1.07 0.02 1.04 – 1.10 <0.001
Observations 1363
Deviance 1330.165
log-Likelihood -665.082
modeloedad
## 
## Call:  glm(formula = newchd ~ age, family = "binomial", data = fram)
## 
## Coefficients:
## (Intercept)          age  
##     -4.8843       0.0658  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1361 Residual
## Null Deviance:       1350 
## Residual Deviance: 1330  AIC: 1330
tab_model(modeloedad)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.01 0.00 – 0.03 <0.001
age 1.07 1.04 – 1.10 <0.001
Observations 1363
R2 Tjur 0.016
#modelo edad quintiles

modeloedadq<- glm(newchd ~ age_q, family = "binomial", data = datos)
t4<-tab_model(modeloedadq, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - edad quantil")
knitr::asis_output(t4$knitr)
Newchd - edad quantil
  newchd
Predictors Odds Ratios std. Error CI p
age_q 1.26 0.06 1.15 – 1.39 <0.001
Observations 1363
Deviance 1328.464
log-Likelihood -664.232
modeloedadq
## 
## Call:  glm(formula = newchd ~ age_q, family = "binomial", data = datos)
## 
## Coefficients:
## (Intercept)        age_q  
##      -2.139        0.233  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1361 Residual
## Null Deviance:       1350 
## Residual Deviance: 1330  AIC: 1330
tab_model(modeloedadq)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.12 0.08 – 0.17 <0.001
age q 1.26 1.15 – 1.39 <0.001
Observations 1363
R2 Tjur 0.017
# edad : la deviance de los modelos cuanti y cuali son diferentes: 1330 para edad cuanti y 1328 para edad quintilos. 

PRESION SISTOLICA

#presión sistólica. Se crea variable sbp_q quintilos de la presión sistólica
datos<- datos %>% mutate(sbp_q = ntile(sbp, 5))
datos$sbp_q1<-datos$sbp_q
datos$sbp_q1<- as.factor(datos$sbp_q1)

datos %>%
  dplyr::group_by(sbp_q) %>%
  dplyr::summarise(tar = mean(newchd)) %>%
  ggplot(aes(x = sbp_q, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", sbp = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## Warning in geom_smooth(method = "lm", sbp = "blue"): Ignoring unknown
## parameters: `sbp`
## `geom_smooth()` using formula = 'y ~ x'

datos %>%
  dplyr::group_by(sbp_q) %>%
  dplyr::summarise(tar = mean(newchd)) %>%
  ggplot(aes(x = sbp_q, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "Cuartiles de sbp", y = "Proporción de eventos CHD") +
  theme_minimal()

#modelo con spb continua
modelosbp <- glm(newchd ~ sbp, family = "binomial", data = fram)
t3<-tab_model(modelosbp, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - edad cuanti")
knitr::asis_output(t3$knitr)
Newchd - edad cuanti
  newchd
Predictors Odds Ratios std. Error CI p
sbp 1.02 0.00 1.01 – 1.02 <0.001
Observations 1363
Deviance 1305.048
log-Likelihood -652.524
modelosbp
## 
## Call:  glm(formula = newchd ~ sbp, family = "binomial", data = fram)
## 
## Coefficients:
## (Intercept)          sbp  
##     -3.7691       0.0156  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1361 Residual
## Null Deviance:       1350 
## Residual Deviance: 1310  AIC: 1310
tab_model(modelosbp)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.02 0.01 – 0.05 <0.001
sbp 1.02 1.01 – 1.02 <0.001
Observations 1363
R2 Tjur 0.036
#modelo con sbp quintilos

modelosbpq<- glm(newchd ~ sbp_q, family = "binomial", data = datos)
t4<-tab_model(modelosbpq, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - sbp quantil")
knitr::asis_output(t4$knitr)
Newchd - sbp quantil
  newchd
Predictors Odds Ratios std. Error CI p
sbp_q 1.39 0.07 1.26 – 1.54 <0.001
Observations 1363
Deviance 1306.486
log-Likelihood -653.243
modelosbpq
## 
## Call:  glm(formula = newchd ~ sbp_q, family = "binomial", data = datos)
## 
## Coefficients:
## (Intercept)        sbp_q  
##      -2.466        0.331  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1361 Residual
## Null Deviance:       1350 
## Residual Deviance: 1310  AIC: 1310
tab_model(modelosbpq)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.08 0.06 – 0.12 <0.001
sbp q 1.39 1.26 – 1.54 <0.001
Observations 1363
R2 Tjur 0.033
# sbp: la deviance de los modelos cuanti y cuali  no son diferentes: 1305 para sbp cuanti y 1306 para sbp quintilos. 
#Variable dbp. Creo dbp_q quintiles
datos<- datos %>% mutate(dbp_q = ntile(dbp, 5))
datos$dbp_q1<-datos$dbp_q
datos$dbp_q1<- as.factor(datos$dbp_q1)

datos %>%
  dplyr::group_by(dbp_q) %>%
  dplyr::summarise(tar = mean(newchd)) %>%
  ggplot(aes(x = dbp_q, y = tar)) +
  geom_point() +
  geom_smooth(method = "lm", dbp = "blue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.5)
## Warning in geom_smooth(method = "lm", dbp = "blue"): Ignoring unknown
## parameters: `dbp`
## `geom_smooth()` using formula = 'y ~ x'

datos %>%
  dplyr::group_by(dbp_q) %>%
  dplyr::summarise(tar = mean(newchd)) %>%
  ggplot(aes(x = dbp_q, y = tar)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
  labs(x = "Cuartiles de dbp", y = "Proporción de eventos CHD") +
  theme_minimal()

#modelo dbp continua
modelodias <- glm(newchd ~ dbp, family = "binomial", data = fram)
t8<-tab_model(modelodias, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - dbp cuanti")
knitr::asis_output(t8$knitr)
Newchd - dbp cuanti
  newchd
Predictors Odds Ratios std. Error CI p
dbp 1.03 0.00 1.02 – 1.04 <0.001
Observations 1363
Deviance 1318.082
log-Likelihood -659.041
modelodias
## 
## Call:  glm(formula = newchd ~ dbp, family = "binomial", data = fram)
## 
## Coefficients:
## (Intercept)          dbp  
##     -3.8583       0.0268  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1361 Residual
## Null Deviance:       1350 
## Residual Deviance: 1320  AIC: 1320
tab_model(modelodias)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.02 0.01 – 0.05 <0.001
dbp 1.03 1.02 – 1.04 <0.001
Observations 1363
R2 Tjur 0.026
#modelo dbp quintiles
modelodiasq<- glm(newchd ~ dbp_q, family = "binomial", data = datos)
t9<-tab_model(modelodiasq, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - dbp quantil")
knitr::asis_output(t9$knitr)
Newchd - dbp quantil
  newchd
Predictors Odds Ratios std. Error CI p
dbp_q 1.32 0.07 1.20 – 1.46 <0.001
Observations 1363
Deviance 1318.712
log-Likelihood -659.356
modelodiasq
## 
## Call:  glm(formula = newchd ~ dbp_q, family = "binomial", data = datos)
## 
## Coefficients:
## (Intercept)        dbp_q  
##       -2.29         0.28  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1361 Residual
## Null Deviance:       1350 
## Residual Deviance: 1320  AIC: 1320
tab_model(modelodiasq)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.10 0.07 – 0.14 <0.001
dbp q 1.32 1.20 – 1.46 <0.001
Observations 1363
R2 Tjur 0.025
# diastolica: la deviance de los modelos cuanti y cuali son diferentes: 1318 para sbp cuanti y 1312 para dbp quintilos. 
#modelo sbp + dbp
modelopresq<- glm(newchd ~ dbp_q + sbp_q, family = "binomial", data = datos)
t10<-tab_model(modelopresq, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - dbp sbp quantil")
knitr::asis_output(t10$knitr)
Newchd - dbp sbp quantil
  newchd
Predictors Odds Ratios std. Error CI p
dbp_q 1.08 0.08 0.94 – 1.25 0.269
sbp_q 1.31 0.10 1.13 – 1.52 <0.001
Observations 1363
Deviance 1305.256
log-Likelihood -652.628
modelopresq
## 
## Call:  glm(formula = newchd ~ dbp_q + sbp_q, family = "binomial", data = datos)
## 
## Coefficients:
## (Intercept)        dbp_q        sbp_q  
##     -2.5301       0.0806       0.2714  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1360 Residual
## Null Deviance:       1350 
## Residual Deviance: 1310  AIC: 1310
tab_model(modelopresq)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.08 0.05 – 0.12 <0.001
dbp q 1.08 0.94 – 1.25 0.269
sbp q 1.31 1.13 – 1.52 <0.001
Observations 1363
R2 Tjur 0.034
vif(modelopresq)
## dbp_q sbp_q 
##  2.12  2.12
# la multicolinealidad es aceptable cuando se evaluan los quintiles: NO SE HACE MULTICOLINEALIDAD EN REGRESION LOGISTICA DIJO EN LA CLASE
modelopres<- glm(newchd ~ dbp + sbp, family = "binomial", data = datos)
t11<-tab_model(modelopres, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Newchd - dbp sbp cuanti")
knitr::asis_output(t11$knitr)
Newchd - dbp sbp cuanti
  newchd
Predictors Odds Ratios std. Error CI p
dbp 1.00 0.01 0.99 – 1.02 0.555
sbp 1.01 0.00 1.01 – 1.02 <0.001
Observations 1363
Deviance 1304.699
log-Likelihood -652.350
modelopres
## 
## Call:  glm(formula = newchd ~ dbp + sbp, family = "binomial", data = datos)
## 
## Coefficients:
## (Intercept)          dbp          sbp  
##    -3.91271      0.00449      0.01385  
## 
## Degrees of Freedom: 1362 Total (i.e. Null);  1360 Residual
## Null Deviance:       1350 
## Residual Deviance: 1300  AIC: 1310
tab_model(modelopres)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.02 0.01 – 0.05 <0.001
dbp 1.00 0.99 – 1.02 0.555
sbp 1.01 1.01 – 1.02 <0.001
Observations 1363
R2 Tjur 0.036
vif(modelopres)
##  dbp  sbp 
## 2.68 2.68
# el vif tambien es aceptable cuando la variable es cuantitativa aunque se modifica la p de dbp y deja de ser significativo. SE CONSIDERA QUE SON COLINEALES Y SE ELIMINA DBP DEL MODELO

CONSTRUCCIÓN DE LOS MODELOS

#MODELO EDAD (1) AGREGANDO SEXO (2). EL COEFICIENTE DE LA EDAD NO SE MMODIFICA CON EL AGREGADO DE SEXO.AMBOS DISMINUYEN LA DEVIANCE EN FORMA ESTADÍSTICAMENTE SIGNIFICATIVA
modelo1 <- glm(newchd~age, data = datos, family = binomial)
modelo2 <- glm(newchd~age+sex, data = datos, family = binomial)

t<-tab_model(modelo1,modelo2, transform=NULL,show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F) 

knitr::asis_output(t$knitr)
  newchd newchd
Predictors Log-Odds std. Error CI p Log-Odds std. Error CI p
age 0.07 0.01 0.04 – 0.09 <0.001 0.07 0.01 0.04 – 0.10 <0.001
sex1 0.72 0.14 0.44 – 0.99 <0.001
Observations 1363 1363
Deviance 1330.165 1303.531
log-Likelihood -665.082 -651.766
anova(modelo2, test="LRT")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: newchd
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev   Pr(>Chi)    
## NULL                  1362       1351               
## age   1     21.1      1361       1330 0.00000440 ***
## sex   1     26.6      1360       1304 0.00000025 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelo3 <- glm(newchd~age+sex+smoke, data = datos, family = binomial)
modelo4 <- glm(newchd~age+sex+smoke+chol, data = datos, family = binomial)
modelo5 <- glm(newchd~age+sex+smoke+chol+sbp, data = datos, family = binomial)
modelo6 <- glm(newchd~age+sex+smoke+chol+sbp+dbp, data = datos, family = binomial)
t15<-tab_model(modelo3,modelo4,modelo5,modelo6, transform=NULL,show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F) 

knitr::asis_output(t15$knitr)
  newchd newchd newchd newchd
Predictors Log-Odds std. Error CI p Log-Odds std. Error CI p Log-Odds std. Error CI p Log-Odds std. Error CI p
age 0.07 0.01 0.04 – 0.10 <0.001 0.07 0.01 0.04 – 0.10 <0.001 0.05 0.02 0.02 – 0.08 <0.001 0.06 0.02 0.03 – 0.09 <0.001
chol 0.01 0.00 0.00 – 0.01 <0.001 0.00 0.00 0.00 – 0.01 0.003 0.00 0.00 0.00 – 0.01 0.003
dbp 0.00 0.01 -0.01 – 0.02 0.535
sbp 0.02 0.00 0.01 – 0.02 <0.001 0.02 0.00 0.01 – 0.02 <0.001
sex1 0.67 0.15 0.38 – 0.96 <0.001 0.79 0.15 0.50 – 1.10 <0.001 0.96 0.16 0.65 – 1.27 <0.001 0.95 0.16 0.64 – 1.27 <0.001
smoke1 0.14 0.15 -0.15 – 0.43 0.348 0.13 0.15 -0.16 – 0.42 0.373 0.17 0.15 -0.12 – 0.47 0.251 0.18 0.15 -0.12 – 0.48 0.241
Observations 1362 1362 1362 1362
Deviance 1302.316 1289.881 1241.541 1241.157
log-Likelihood -651.158 -644.941 -620.771 -620.579
anova(modelo6, test="LRT")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: newchd
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev        Pr(>Chi)    
## NULL                   1361       1351                    
## age    1     21.1      1360       1330 0.0000043316375 ***
## sex    1     26.5      1359       1303 0.0000002637865 ***
## smoke  1      0.9      1358       1302         0.34792    
## chol   1     12.4      1357       1290         0.00042 ***
## sbp    1     48.3      1356       1242 0.0000000000036 ***
## dbp    1      0.4      1355       1241         0.53535    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#smoke no es significativo. dbp no es significativo en el modelo. ni smoke ni dbp mejoran la deviance del modelo
modelo7 <- glm(newchd~age+sex+chol+sbp+smoke, data = datos, family = binomial)
t16<-tab_model(modelo7,transform=NULL,show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F) 

knitr::asis_output(t16$knitr)
  newchd
Predictors Log-Odds std. Error CI p
age 0.05 0.02 0.02 – 0.08 <0.001
chol 0.00 0.00 0.00 – 0.01 0.003
sbp 0.02 0.00 0.01 – 0.02 <0.001
sex1 0.96 0.16 0.65 – 1.27 <0.001
smoke1 0.17 0.15 -0.12 – 0.47 0.251
Observations 1362
Deviance 1241.541
log-Likelihood -620.771
anova(modelo7, test="LRT")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: newchd
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev        Pr(>Chi)    
## NULL                   1361       1351                    
## age    1     21.1      1360       1330 0.0000043316375 ***
## sex    1     26.5      1359       1303 0.0000002637865 ***
## chol   1     12.5      1358       1291          0.0004 ***
## sbp    1     47.8      1357       1243 0.0000000000047 ***
## smoke  1      1.3      1356       1242          0.2511    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#mi elección es el modelo 7q SIN dbp con SMOKE tiene deviance 1249

#modelo 7q tiene deviance 1249 con edad q. ESTE MODELO ELIJO
modelo7q <- glm(newchd~age_q+sex+chol_q+sbp_q+smoke, data = datos, family = binomial)
t17<-tab_model(modelo7q,transform=NULL,show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F)

knitr::asis_output(t17$knitr)
  newchd
Predictors Log-Odds std. Error CI p
age_q 0.19 0.05 0.09 – 0.29 <0.001
chol_q 0.10 0.05 0.00 – 0.20 0.048
sbp_q 0.35 0.05 0.24 – 0.46 <0.001
sex1 0.89 0.16 0.58 – 1.19 <0.001
smoke1 0.17 0.15 -0.12 – 0.47 0.253
Observations 1362
Deviance 1249.781
log-Likelihood -624.890
anova(modelo7q, test="LRT")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: newchd
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev       Pr(>Chi)    
## NULL                    1361       1351                   
## age_q   1     22.9      1360       1328 0.000001735634 ***
## sex     1     26.9      1359       1301 0.000000217485 ***
## chol_q  1      6.3      1358       1295          0.012 *  
## sbp_q   1     43.6      1357       1251 0.000000000039 ***
## smoke   1      1.3      1356       1250          0.253    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(modelo7)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.00 0.00 – 0.00 <0.001
age 1.06 1.02 – 1.09 <0.001
sex [1] 2.61 1.91 – 3.58 <0.001
chol 1.00 1.00 – 1.01 0.003
sbp 1.02 1.01 – 1.02 <0.001
smoke [1] 1.19 0.88 – 1.60 0.251
Observations 1362
R2 Tjur 0.080
tab_model(modelo7q)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.02 0.01 – 0.04 <0.001
age q 1.21 1.09 – 1.34 <0.001
sex [1] 2.42 1.79 – 3.30 <0.001
chol q 1.11 1.00 – 1.23 0.048
sbp q 1.42 1.28 – 1.58 <0.001
smoke [1] 1.19 0.88 – 1.60 0.253
Observations 1362
R2 Tjur 0.073
#Modelo 7q1 deviance 1249 edad sin q
modelo7q1 <- glm(newchd~age+sex+chol_q+sbp_q+smoke, data = datos, family = binomial)
t77<-tab_model(modelo7q1,transform=NULL,show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F)

knitr::asis_output(t77$knitr)
  newchd
Predictors Log-Odds std. Error CI p
age 0.05 0.02 0.02 – 0.08 <0.001
chol_q 0.10 0.05 0.00 – 0.20 0.049
sbp_q 0.35 0.05 0.25 – 0.46 <0.001
sex1 0.88 0.16 0.58 – 1.19 <0.001
smoke1 0.18 0.15 -0.12 – 0.48 0.234
Observations 1362
Deviance 1249.909
log-Likelihood -624.955
anova(modelo7q1, test="LRT")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: newchd
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev       Pr(>Chi)    
## NULL                    1361       1351                   
## age     1     21.1      1360       1330 0.000004331638 ***
## sex     1     26.5      1359       1303 0.000000263786 ***
## chol_q  1      6.3      1358       1297          0.012 *  
## sbp_q   1     45.6      1357       1251 0.000000000015 ***
## smoke   1      1.4      1356       1250          0.234    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(modelo7q1)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.00 0.00 – 0.01 <0.001
age 1.06 1.03 – 1.09 <0.001
sex [1] 2.42 1.78 – 3.29 <0.001
chol q 1.11 1.00 – 1.23 0.049
sbp q 1.43 1.28 – 1.59 <0.001
smoke [1] 1.20 0.89 – 1.61 0.234
Observations 1362
R2 Tjur 0.072
#Quiero probar interaccioón entre edadq y sbpq
modelo7qint <- glm(newchd~age_q+sex+chol_q+sbp_q+smoke+age_q:sbp_q, data = datos, family = binomial)
t17i<-tab_model(modelo7qint,transform=NULL,show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F)

knitr::asis_output(t17i$knitr)
  newchd
Predictors Log-Odds std. Error CI p
age_q 0.10 0.14 -0.16 – 0.37 0.443
age_q:sbp_q 0.03 0.04 -0.05 – 0.10 0.505
chol_q 0.10 0.05 0.00 – 0.20 0.047
sbp_q 0.27 0.13 0.01 – 0.53 0.042
sex1 0.89 0.16 0.59 – 1.20 <0.001
smoke1 0.18 0.15 -0.12 – 0.47 0.246
Observations 1362
Deviance 1249.337
log-Likelihood -624.669
anova(modelo7qint, test="LRT")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: newchd
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev       Pr(>Chi)    
## NULL                         1361       1351                   
## age_q        1     22.9      1360       1328 0.000001735634 ***
## sex          1     26.9      1359       1301 0.000000217485 ***
## chol_q       1      6.3      1358       1295          0.012 *  
## sbp_q        1     43.6      1357       1251 0.000000000039 ***
## smoke        1      1.3      1356       1250          0.253    
## age_q:sbp_q  1      0.4      1355       1249          0.506    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(modelo7qint)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.02 0.01 – 0.06 <0.001
age q 1.11 0.85 – 1.45 0.443
sex [1] 2.44 1.80 – 3.32 <0.001
chol q 1.11 1.00 – 1.23 0.047
sbp q 1.31 1.01 – 1.70 0.042
smoke [1] 1.19 0.89 – 1.61 0.246
age q × sbp q 1.03 0.95 – 1.10 0.505
Observations 1362
R2 Tjur 0.073
#edad : sbp no significativo. No lo dejo. no agrega deviance
#HOSMER & LEMERSHOW. LA P ME DA 0.08. NO RECHAZA LA HIPÓTESIS NULA.ES COMO QUIERO QUE ME DE 
pred_modelo <- as.data.frame(predict(modelo7q, datos, type = "response")) %>% rename(pred = `predict(modelo7q, datos, type = "response")`)

pred_modelo$verdad <- datos$newchd

H_L <- HLtest(modelo7q,10)
HL_DF <- cbind(as.data.frame(H_L$table), indice = as.factor(c(1,2,3,4,5,6,7,8,9,10)))

H_L
## Hosmer and Lemeshow Goodness-of-Fit Test 
## 
## Call:
## glm(formula = newchd ~ age_q + sex + chol_q + sbp_q + smoke, 
##     family = binomial, data = datos)
##  ChiSquare df P_value
##       13.8  8  0.0881
HL_DF[,1:4]
##                cut total obs   exp
## 1  [0.0345,0.0708]   140 135 132.4
## 2  (0.0708,0.0995]   133 128 121.5
## 3   (0.0995,0.127]   137 117 121.5
## 4    (0.127,0.152]   135 115 116.2
## 5    (0.152,0.177]   139 109 116.2
## 6    (0.177,0.207]   133 114 107.7
## 7    (0.207,0.239]   136 101 105.8
## 8    (0.239,0.281]   137  97 101.3
## 9    (0.281,0.353]   135  94  92.9
## 10   (0.353,0.569]   137  84  78.6
HL_DF %>% knitr::kable()
cut total obs exp chi indice
[0.0345,0.0708] 140 135 132.4 0.228 1
(0.0708,0.0995] 133 128 121.5 0.587 2
(0.0995,0.127] 137 117 121.5 -0.404 3
(0.127,0.152] 135 115 116.2 -0.116 4
(0.152,0.177] 139 109 116.2 -0.669 5
(0.177,0.207] 133 114 107.7 0.610 6
(0.207,0.239] 136 101 105.8 -0.463 7
(0.239,0.281] 137 97 101.3 -0.428 8
(0.281,0.353] 135 94 92.9 0.116 9
(0.353,0.569] 137 84 78.6 0.613 10
ggplot(HL_DF, aes(x=indice))+geom_point(aes(y=exp, color = "blue"))+geom_point(aes(y=obs, color = "red"))+scale_color_manual(labels = c("Esperado", "Observado"), values = c("blue", "red"))

#slope e intercept
pred_modelo <- pred_modelo[-1258,]# esta parece ser una observacion duplicada

a <- valProbggplot(pred_modelo$pred, pred_modelo$verdad)
a$ggPlot

#a$Calibration
#curva roc. AUC 0.69
roc<-roc(pred_modelo$verdad, pred_modelo$pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc, main = "Curva ROC", col = "red", lwd = 3)

text(0.8, 0.2, paste("AUC ROC =", round(roc$auc, 3)), adj = c(0, 1))

#punto de corte 40%
pred_modelo <- pred_modelo %>% dplyr::mutate(prediccion = ifelse(pred > 0.4, 1,0))

tab <- tab_xtab(pred_modelo$prediccion, pred_modelo$verdad ,tdcol.row = "black",tdcol.n = "black", title = "Matriz de Confusion - cutoff =30") ## ojo si usan el RStudio con fondo blanco cambien el color de los numeros o no van a ver nada...

knitr::asis_output(tab$knitr)
Matriz de Confusion - cutoff =30
prediccion verdad Total
0 1
0 1043 233 1276
1 51 35 86
Total 1094 268 1362
χ2=24.264 · df=1 · &phi=0.137 · p=0.000
cat("Sensibilidad:", paste0(round(caret::sensitivity(factor(pred_modelo$verdad), factor(pred_modelo$prediccion)),2)))
## Sensibilidad: 0.82
cat("Especificidad:", paste0(round(caret::specificity(factor(pred_modelo$verdad), factor(pred_modelo$prediccion)),2)))
## Especificidad: 0.41
SENS=roc$sensitivities
ESPEC=roc$specificities
CUTOFF=roc$thresholds

d <- cbind(SENS,ESPEC,CUTOFF)
d%>% knitr::kable()
SENS ESPEC CUTOFF
1.000 0.000 -Inf
1.000 0.007 0.036
1.000 0.015 0.039
1.000 0.018 0.041
1.000 0.021 0.042
0.996 0.027 0.043
0.996 0.031 0.045
0.996 0.032 0.046
0.996 0.034 0.047
0.996 0.039 0.048
0.996 0.042 0.049
0.996 0.047 0.050
0.993 0.050 0.050
0.993 0.052 0.051
0.993 0.053 0.052
0.993 0.055 0.053
0.993 0.058 0.054
0.993 0.060 0.056
0.993 0.065 0.057
0.985 0.067 0.058
0.985 0.069 0.058
0.985 0.074 0.059
0.985 0.075 0.059
0.985 0.077 0.060
0.985 0.080 0.060
0.985 0.080 0.061
0.985 0.081 0.062
0.985 0.085 0.063
0.985 0.090 0.064
0.985 0.091 0.064
0.985 0.096 0.065
0.985 0.099 0.065
0.985 0.101 0.066
0.985 0.105 0.066
0.985 0.109 0.067
0.985 0.110 0.068
0.985 0.112 0.069
0.981 0.116 0.070
0.981 0.117 0.070
0.981 0.120 0.071
0.981 0.123 0.071
0.981 0.124 0.072
0.981 0.125 0.072
0.981 0.131 0.073
0.981 0.133 0.074
0.981 0.135 0.075
0.981 0.138 0.076
0.981 0.142 0.076
0.978 0.144 0.077
0.978 0.146 0.078
0.978 0.148 0.079
0.978 0.153 0.080
0.978 0.156 0.080
0.978 0.157 0.081
0.978 0.161 0.081
0.978 0.162 0.082
0.978 0.164 0.082
0.978 0.168 0.083
0.978 0.170 0.083
0.978 0.171 0.084
0.974 0.174 0.085
0.974 0.176 0.086
0.970 0.180 0.086
0.970 0.183 0.087
0.970 0.186 0.088
0.970 0.190 0.088
0.970 0.193 0.089
0.970 0.194 0.089
0.970 0.196 0.090
0.970 0.199 0.090
0.970 0.200 0.091
0.970 0.201 0.092
0.970 0.207 0.093
0.970 0.208 0.093
0.970 0.209 0.093
0.970 0.218 0.094
0.970 0.218 0.095
0.970 0.223 0.095
0.970 0.224 0.095
0.970 0.229 0.096
0.970 0.230 0.096
0.966 0.231 0.097
0.966 0.233 0.097
0.966 0.236 0.098
0.963 0.237 0.098
0.963 0.239 0.099
0.963 0.240 0.100
0.955 0.244 0.100
0.955 0.247 0.101
0.951 0.250 0.102
0.951 0.251 0.102
0.940 0.252 0.102
0.933 0.257 0.103
0.933 0.258 0.104
0.933 0.259 0.104
0.933 0.260 0.104
0.933 0.261 0.105
0.933 0.265 0.106
0.933 0.267 0.106
0.933 0.271 0.107
0.933 0.273 0.109
0.929 0.277 0.109
0.925 0.277 0.110
0.925 0.280 0.111
0.925 0.282 0.111
0.925 0.282 0.112
0.925 0.284 0.112
0.922 0.287 0.113
0.922 0.291 0.114
0.922 0.295 0.114
0.918 0.297 0.114
0.918 0.298 0.115
0.918 0.301 0.115
0.914 0.301 0.116
0.914 0.304 0.116
0.914 0.308 0.117
0.914 0.311 0.118
0.914 0.312 0.119
0.914 0.314 0.119
0.914 0.314 0.120
0.914 0.316 0.121
0.903 0.320 0.121
0.899 0.323 0.122
0.896 0.326 0.123
0.896 0.327 0.123
0.896 0.329 0.123
0.896 0.335 0.124
0.892 0.337 0.124
0.892 0.339 0.125
0.892 0.340 0.126
0.888 0.340 0.127
0.888 0.346 0.127
0.888 0.347 0.127
0.884 0.348 0.128
0.884 0.351 0.129
0.884 0.354 0.130
0.884 0.355 0.130
0.884 0.362 0.131
0.881 0.365 0.131
0.881 0.367 0.132
0.881 0.373 0.132
0.881 0.375 0.132
0.881 0.378 0.133
0.877 0.378 0.133
0.877 0.382 0.134
0.873 0.384 0.134
0.869 0.387 0.134
0.866 0.389 0.135
0.866 0.390 0.136
0.866 0.393 0.136
0.862 0.396 0.137
0.862 0.398 0.137
0.862 0.401 0.138
0.851 0.403 0.139
0.851 0.405 0.139
0.851 0.407 0.140
0.847 0.410 0.141
0.843 0.411 0.141
0.843 0.412 0.142
0.836 0.417 0.143
0.836 0.419 0.143
0.836 0.422 0.144
0.836 0.425 0.144
0.836 0.426 0.144
0.836 0.428 0.144
0.828 0.430 0.145
0.828 0.432 0.146
0.828 0.434 0.146
0.828 0.437 0.147
0.828 0.439 0.148
0.825 0.441 0.148
0.825 0.443 0.149
0.817 0.444 0.150
0.813 0.446 0.151
0.813 0.448 0.151
0.813 0.450 0.152
0.813 0.452 0.152
0.810 0.458 0.153
0.810 0.460 0.153
0.810 0.464 0.154
0.806 0.467 0.154
0.806 0.468 0.154
0.799 0.469 0.155
0.799 0.470 0.155
0.795 0.472 0.156
0.795 0.474 0.156
0.787 0.477 0.157
0.787 0.480 0.157
0.787 0.483 0.158
0.787 0.484 0.159
0.787 0.491 0.159
0.787 0.492 0.159
0.776 0.492 0.160
0.776 0.493 0.161
0.772 0.495 0.161
0.769 0.498 0.162
0.765 0.500 0.163
0.754 0.502 0.163
0.746 0.503 0.164
0.746 0.505 0.165
0.746 0.506 0.165
0.743 0.509 0.166
0.743 0.510 0.167
0.739 0.512 0.167
0.739 0.516 0.167
0.739 0.516 0.168
0.739 0.520 0.169
0.739 0.523 0.169
0.739 0.527 0.170
0.728 0.533 0.171
0.720 0.536 0.172
0.713 0.537 0.173
0.713 0.539 0.174
0.713 0.540 0.174
0.709 0.541 0.175
0.705 0.546 0.175
0.705 0.548 0.176
0.705 0.549 0.176
0.701 0.552 0.177
0.698 0.554 0.177
0.698 0.556 0.178
0.694 0.558 0.178
0.690 0.559 0.179
0.687 0.564 0.180
0.687 0.567 0.180
0.679 0.569 0.180
0.675 0.572 0.181
0.675 0.573 0.182
0.675 0.576 0.182
0.668 0.581 0.182
0.668 0.583 0.183
0.660 0.583 0.184
0.657 0.585 0.185
0.657 0.588 0.186
0.653 0.588 0.186
0.653 0.591 0.187
0.653 0.594 0.188
0.653 0.596 0.189
0.653 0.597 0.190
0.653 0.601 0.191
0.649 0.608 0.191
0.649 0.609 0.192
0.649 0.611 0.193
0.646 0.611 0.193
0.646 0.612 0.193
0.646 0.614 0.194
0.646 0.621 0.195
0.646 0.623 0.195
0.646 0.625 0.196
0.646 0.626 0.197
0.646 0.628 0.198
0.646 0.630 0.199
0.646 0.631 0.200
0.646 0.633 0.200
0.638 0.634 0.201
0.634 0.634 0.202
0.634 0.636 0.202
0.634 0.638 0.202
0.634 0.642 0.203
0.634 0.643 0.204
0.634 0.645 0.204
0.631 0.647 0.205
0.631 0.650 0.205
0.631 0.652 0.206
0.631 0.654 0.206
0.631 0.655 0.207
0.631 0.656 0.207
0.623 0.664 0.207
0.623 0.665 0.208
0.619 0.668 0.209
0.616 0.671 0.209
0.616 0.673 0.209
0.616 0.674 0.210
0.608 0.675 0.211
0.604 0.682 0.212
0.601 0.684 0.213
0.601 0.686 0.214
0.601 0.686 0.215
0.601 0.687 0.216
0.597 0.687 0.216
0.597 0.688 0.218
0.593 0.693 0.219
0.590 0.694 0.219
0.575 0.699 0.220
0.571 0.699 0.221
0.563 0.700 0.221
0.563 0.703 0.222
0.563 0.706 0.223
0.563 0.707 0.224
0.560 0.708 0.224
0.556 0.709 0.225
0.541 0.716 0.226
0.541 0.718 0.227
0.534 0.719 0.229
0.534 0.721 0.231
0.534 0.728 0.232
0.526 0.729 0.234
0.522 0.730 0.234
0.519 0.732 0.235
0.519 0.737 0.236
0.507 0.739 0.237
0.504 0.741 0.237
0.504 0.747 0.237
0.504 0.749 0.238
0.500 0.749 0.239
0.500 0.753 0.239
0.496 0.755 0.240
0.496 0.756 0.241
0.496 0.757 0.242
0.493 0.758 0.242
0.485 0.759 0.243
0.470 0.764 0.245
0.470 0.765 0.245
0.463 0.766 0.246
0.459 0.766 0.249
0.459 0.768 0.250
0.444 0.771 0.251
0.440 0.776 0.253
0.440 0.778 0.253
0.437 0.780 0.253
0.437 0.782 0.254
0.429 0.790 0.257
0.429 0.791 0.258
0.425 0.791 0.260
0.425 0.792 0.261
0.422 0.793 0.263
0.418 0.796 0.264
0.414 0.799 0.265
0.403 0.803 0.268
0.403 0.807 0.269
0.403 0.809 0.270
0.403 0.810 0.271
0.396 0.817 0.272
0.392 0.820 0.273
0.388 0.821 0.274
0.381 0.824 0.275
0.362 0.831 0.277
0.362 0.833 0.279
0.358 0.834 0.280
0.354 0.835 0.281
0.351 0.837 0.283
0.351 0.838 0.285
0.347 0.843 0.286
0.340 0.843 0.287
0.340 0.844 0.287
0.340 0.846 0.289
0.336 0.848 0.290
0.332 0.851 0.290
0.325 0.856 0.292
0.310 0.862 0.293
0.310 0.863 0.293
0.302 0.863 0.295
0.299 0.865 0.296
0.295 0.865 0.298
0.295 0.867 0.300
0.295 0.868 0.302
0.295 0.872 0.303
0.295 0.876 0.305
0.287 0.877 0.307
0.284 0.878 0.308
0.276 0.878 0.309
0.276 0.879 0.310
0.272 0.881 0.311
0.261 0.882 0.312
0.261 0.883 0.313
0.250 0.886 0.315
0.250 0.887 0.315
0.250 0.888 0.316
0.246 0.890 0.318
0.243 0.891 0.319
0.235 0.892 0.321
0.235 0.893 0.323
0.235 0.897 0.324
0.231 0.899 0.326
0.228 0.900 0.327
0.224 0.900 0.329
0.224 0.902 0.330
0.224 0.904 0.332
0.224 0.907 0.334
0.224 0.909 0.335
0.220 0.912 0.338
0.216 0.912 0.340
0.213 0.913 0.342
0.209 0.916 0.343
0.209 0.918 0.344
0.209 0.919 0.346
0.209 0.920 0.348
0.201 0.923 0.350
0.198 0.923 0.352
0.194 0.924 0.353
0.187 0.925 0.355
0.183 0.925 0.357
0.183 0.926 0.360
0.179 0.927 0.363
0.179 0.928 0.365
0.168 0.931 0.367
0.168 0.935 0.369
0.164 0.936 0.373
0.160 0.936 0.378
0.157 0.939 0.380
0.153 0.940 0.382
0.149 0.942 0.383
0.142 0.943 0.384
0.142 0.944 0.386
0.142 0.946 0.387
0.138 0.948 0.389
0.138 0.949 0.391
0.138 0.952 0.393
0.138 0.953 0.396
0.131 0.953 0.399
0.127 0.955 0.403
0.123 0.958 0.406
0.123 0.960 0.409
0.116 0.963 0.413
0.116 0.964 0.418
0.108 0.966 0.423
0.101 0.968 0.425
0.097 0.968 0.427
0.093 0.970 0.429
0.093 0.972 0.431
0.090 0.973 0.433
0.090 0.975 0.435
0.090 0.979 0.438
0.086 0.979 0.443
0.086 0.981 0.448
0.086 0.984 0.450
0.086 0.985 0.453
0.086 0.986 0.456
0.078 0.987 0.463
0.075 0.988 0.470
0.067 0.989 0.474
0.060 0.989 0.476
0.056 0.989 0.478
0.052 0.991 0.481
0.049 0.993 0.488
0.041 0.993 0.495
0.034 0.995 0.499
0.022 0.996 0.510
0.022 0.998 0.521
0.011 0.998 0.525
0.007 0.999 0.535
0.004 0.999 0.557
0.000 1.000 Inf

Otro modelo. Voy a poner colesterol dicotómica según el último quintilo. Modelo col cat con 2 categorias mayor y menor a 271 correspondiente al último quintilo MODELOCOLCAT7 DEVIANCE 1245

summary(datos$chol_q1)
##   1   2   3   4   5 
## 273 273 273 272 272
describe(datos$chol)
## datos$chol 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##     1363        0      179        1    234.6      233     51.8    166.0 
##      .10      .25      .50      .75      .90      .95 
##    178.0    200.0    230.0    264.5    296.0    317.0 
## 
## lowest :  96 121 126 133 134, highest: 386 400 405 409 430
quantile(datos$chol,0.80) #da 271
## 80% 
## 271
datos <- datos %>% mutate(chol_cat=as.factor(case_when(chol>=271 ~"1",
chol<271 ~"0")))
datos$chol_cat <- factor (datos$chol_cat, levels= c ("0","1"))
datos$chol_cat
##    [1] 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0
##   [38] 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 1 0 1 0 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0
##   [75] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
##  [112] 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 0 0 0
##  [149] 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
##  [186] 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
##  [223] 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 0 0 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0
##  [260] 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0
##  [297] 0 0 0 0 1 1 0 1 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 1 0 1 0 0 1
##  [334] 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
##  [371] 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 1 0 0 0 1 0 0 0 0
##  [408] 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
##  [445] 1 1 0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [482] 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 0 0
##  [519] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1
##  [556] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0
##  [593] 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0
##  [630] 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1
##  [667] 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0 1 0 0
##  [704] 1 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [741] 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 0
##  [778] 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0
##  [815] 0 0 0 0 0 0 0 0 0 1 0 0 1 1 1 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1
##  [852] 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 1 0 0 1 0 1 1 0 0
##  [889] 0 1 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0
##  [926] 1 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
##  [963] 0 0 1 0 1 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0
## [1000] 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## [1037] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 1 0 0 0
## [1074] 0 1 0 0 1 0 0 1 0 1 0 1 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1
## [1111] 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 1 0 0 1 1 1 0 1 0 0 0 1 0 1 0
## [1148] 1 0 0 0 0 0 1 0 1 0 0 1 0 0 1 0 1 0 1 1 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1
## [1185] 0 0 1 0 1 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 1 0 1 0
## [1222] 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## [1259] 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 1 1 0 0 1 0 0 1 0 0 0 1 0 0 0
## [1296] 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [1333] 0 1 0 0 0 1 0 0 1 0 1 1 0 0 0 0 0 0 1 1 0 1 0 0 0 1 1 0 0 0 0
## Levels: 0 1
describe(datos$chol_cat)
## datos$chol_cat 
##        n  missing distinct 
##     1363        0        2 
##                       
## Value          0     1
## Frequency   1088   275
## Proportion 0.798 0.202
modelocolcat<-glm (newchd~chol_cat, data=datos, family=binomial)#bivariado colesterol cat
tab_model(modelocolcat)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.22 0.19 – 0.26 <0.001
chol cat [1] 1.53 1.12 – 2.09 0.007
Observations 1363
R2 Tjur 0.005
summary(modelocolcat)
## 
## Call:
## glm(formula = newchd ~ chol_cat, family = binomial, data = datos)
## 
## Coefficients:
##             Estimate Std. Error z value            Pr(>|z|)    
## (Intercept)  -1.5030     0.0786  -19.13 <0.0000000000000002 ***
## chol_cat1     0.4284     0.1592    2.69              0.0071 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1351.2  on 1362  degrees of freedom
## Residual deviance: 1344.3  on 1361  degrees of freedom
## AIC: 1348
## 
## Number of Fisher Scoring iterations: 4
modelocolcat7<- glm(newchd~age+sex+chol_cat+sbp_q+smoke, data = datos, family = binomial)
tcolcat<-tab_model(modelocolcat7,transform=NULL,show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F)

knitr::asis_output(tcolcat$knitr)
  newchd
Predictors Log-Odds std. Error CI p
age 0.06 0.02 0.03 – 0.09 <0.001
chol_cat1 0.52 0.17 0.19 – 0.85 0.002
sbp_q 0.36 0.05 0.25 – 0.46 <0.001
sex1 0.90 0.16 0.60 – 1.21 <0.001
smoke1 0.18 0.15 -0.12 – 0.48 0.239
Observations 1362
Deviance 1244.653
log-Likelihood -622.326
anova(modelocolcat7, test="LRT")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: newchd
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev      Pr(>Chi)    
## NULL                      1361       1351                  
## age       1     21.1      1360       1330 0.00000433164 ***
## sex       1     26.5      1359       1303 0.00000026379 ***
## chol_cat  1     10.9      1358       1292       0.00096 ***
## sbp_q     1     46.3      1357       1246 0.00000000001 ***
## smoke     1      1.4      1356       1245       0.23933    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(modelocolcat7)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.00 0.00 – 0.01 <0.001
age 1.06 1.03 – 1.09 <0.001
sex [1] 2.46 1.81 – 3.34 <0.001
chol cat [1] 1.68 1.20 – 2.34 0.002
sbp q 1.43 1.29 – 1.59 <0.001
smoke [1] 1.20 0.89 – 1.61 0.239
Observations 1362
R2 Tjur 0.076

VOY A CATEGORIZAR DBP POR EL ÚLTIMO CUARTIL QUE AUMENTA A 0.33 Y VOY A INTENTAR METERLO EN EL MODELO. no da significaivo ni es confundidor de sbp

summary(datos$dbp_q1)
##   1   2   3   4   5 
## 273 273 273 272 272
describe(datos$dbp_q1)
## datos$dbp_q1 
##        n  missing distinct 
##     1363        0        5 
##                               
## Value        1   2   3   4   5
## Frequency  273 273 273 272 272
## Proportion 0.2 0.2 0.2 0.2 0.2
quantile(datos$dbp,0.80) #da 100
## 80% 
## 100
datos <- datos %>% mutate(dbp_cat=as.factor(case_when(dbp>=100 ~"1",
dbp<100 ~"0")))
datos$dbp_cat <- factor (datos$dbp_cat, levels= c ("0","1"))
describe(datos$dbp_cat)
## datos$dbp_cat 
##        n  missing distinct 
##     1363        0        2 
##                     
## Value         0    1
## Frequency  1036  327
## Proportion 0.76 0.24
modelodbpcat<-glm (newchd~dbp_cat, data=datos, family=binomial)#bivariado colesterol cat
tab_model(modelodbpcat)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.20 0.17 – 0.23 <0.001
dbp cat [1] 2.18 1.63 – 2.90 <0.001
Observations 1363
R2 Tjur 0.021
summary(modelodbpcat)
## 
## Call:
## glm(formula = newchd ~ dbp_cat, family = binomial, data = datos)
## 
## Coefficients:
##             Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)  -1.6281     0.0839   -19.4 < 0.0000000000000002 ***
## dbp_cat1      0.7793     0.1470     5.3           0.00000011 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1351.2  on 1362  degrees of freedom
## Residual deviance: 1324.3  on 1361  degrees of freedom
## AIC: 1328
## 
## Number of Fisher Scoring iterations: 4
modelodbpcat7<- glm(newchd~age+sex+chol_cat+sbp_q+smoke+dbp_cat, data = datos, family = binomial)
tdbpcat<-tab_model(modelodbpcat7,transform=NULL,show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F)

knitr::asis_output(tdbpcat$knitr)
  newchd
Predictors Log-Odds std. Error CI p
age 0.06 0.02 0.03 – 0.09 <0.001
chol_cat1 0.51 0.17 0.18 – 0.84 0.003
dbp_cat1 0.24 0.19 -0.14 – 0.62 0.218
sbp_q 0.31 0.07 0.18 – 0.44 <0.001
sex1 0.89 0.16 0.59 – 1.20 <0.001
smoke1 0.18 0.15 -0.12 – 0.48 0.242
Observations 1362
Deviance 1243.135
log-Likelihood -621.568
anova(modelodbpcat7, test="LRT")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: newchd
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev      Pr(>Chi)    
## NULL                      1361       1351                  
## age       1     21.1      1360       1330 0.00000433164 ***
## sex       1     26.5      1359       1303 0.00000026379 ***
## chol_cat  1     10.9      1358       1292       0.00096 ***
## sbp_q     1     46.3      1357       1246 0.00000000001 ***
## smoke     1      1.4      1356       1245       0.23933    
## dbp_cat   1      1.5      1355       1243       0.21799    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(modelodbpcat7)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.00 0.00 – 0.01 <0.001
age 1.06 1.03 – 1.09 <0.001
sex [1] 2.44 1.80 – 3.32 <0.001
chol cat [1] 1.67 1.19 – 2.32 0.003
sbp q 1.36 1.19 – 1.55 <0.001
smoke [1] 1.19 0.89 – 1.61 0.242
dbp cat [1] 1.27 0.87 – 1.85 0.218
Observations 1362
R2 Tjur 0.078

Ahora voy a probar edad segun el quantil 60 dado que ahi sube el % de enf cardiovascular. Da 54 años. deviance 1248. Pero me quedo con esta

summary(datos$age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    45.0    48.0    52.0    52.4    56.0    62.0
describe(datos$age)
## datos$age 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##     1363        0       18    0.996    52.38     52.5     5.51       45 
##      .10      .25      .50      .75      .90      .95 
##       46       48       52       56       59       60 
##                                                                             
## Value         45    46    47    48    49    50    51    52    53    54    55
## Frequency    109    87    84    79    81    92    87    99    67    85    79
## Proportion 0.080 0.064 0.062 0.058 0.059 0.067 0.064 0.073 0.049 0.062 0.058
##                                                     
## Value         56    57    58    59    60    61    62
## Frequency     87    77    74    66    65    34    11
## Proportion 0.064 0.056 0.054 0.048 0.048 0.025 0.008
## 
## For the frequency table, variable is rounded to the nearest 0
quantile(datos$age,0.60) #da 54
## 60% 
##  54
datos <- datos %>% mutate(age_cat=as.factor(case_when(age>=54 ~"1",
age<54 ~"0")))
datos$age_cat <- factor (datos$age_cat, levels= c ("0","1"))
describe(datos$age_cat)
## datos$age_cat 
##        n  missing distinct 
##     1363        0        2 
##                       
## Value          0     1
## Frequency    785   578
## Proportion 0.576 0.424
modeloagecat<-glm (newchd~age_cat, data=datos, family=binomial)#bivariado edad cat
tab_model(modeloagecat)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.19 0.15 – 0.23 <0.001
age cat [1] 1.77 1.35 – 2.32 <0.001
Observations 1363
R2 Tjur 0.013
summary(modeloagecat)
## 
## Call:
## glm(formula = newchd ~ age_cat, family = binomial, data = datos)
## 
## Coefficients:
##             Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)  -1.6735     0.0979  -17.10 < 0.0000000000000002 ***
## age_cat1      0.5702     0.1372    4.16             0.000032 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1351.2  on 1362  degrees of freedom
## Residual deviance: 1333.9  on 1361  degrees of freedom
## AIC: 1338
## 
## Number of Fisher Scoring iterations: 4
modeloagecat7<- glm(newchd~age_cat+sex+chol_cat+sbp_q+smoke, data = datos, family = binomial)
tagecat<-tab_model(modeloagecat7,transform=NULL,show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F)

knitr::asis_output(tagecat$knitr)
  newchd
Predictors Log-Odds std. Error CI p
age_cat1 0.46 0.14 0.18 – 0.74 0.001
chol_cat1 0.52 0.17 0.18 – 0.85 0.002
sbp_q 0.36 0.05 0.25 – 0.46 <0.001
sex1 0.90 0.16 0.60 – 1.21 <0.001
smoke1 0.15 0.15 -0.14 – 0.45 0.315
Observations 1362
Deviance 1247.689
log-Likelihood -623.844
anova(modeloagecat7, test="LRT")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: newchd
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev        Pr(>Chi)    
## NULL                      1361       1351                    
## age_cat   1     17.5      1360       1333 0.0000291615794 ***
## sex       1     26.7      1359       1307 0.0000002427154 ***
## chol_cat  1     10.8      1358       1296         0.00099 ***
## sbp_q     1     47.1      1357       1249 0.0000000000066 ***
## smoke     1      1.0      1356       1248         0.31563    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(modeloagecat7)
  newchd
Predictors Odds Ratios CI p
(Intercept) 0.03 0.02 – 0.05 <0.001
age cat [1] 1.59 1.20 – 2.11 0.001
sex [1] 2.46 1.82 – 3.35 <0.001
chol cat [1] 1.68 1.20 – 2.33 0.002
sbp q 1.43 1.29 – 1.59 <0.001
smoke [1] 1.16 0.87 – 1.56 0.315
Observations 1362
R2 Tjur 0.075

LE VOY A HACER LAS COMPROBACIONES A MODELOAGECAT7

#HOSMER & LEMERSHOW. . p 0.18 NO RECHAZA LA HIPÓTESIS NULA.ES COMO QUIERO QUE ME DE 
pred_modeloagecat7 <- as.data.frame(predict(modeloagecat7, datos, type = "response")) %>% rename(pred = `predict(modeloagecat7, datos, type = "response")`)

pred_modeloagecat7$verdad <- datos$newchd

H_L <- HLtest(modeloagecat7,10)
HL_DF <- cbind(as.data.frame(H_L$table), indice = as.factor(c(1,2,3,4,5,6,7,8,9,10)))

H_L
## Hosmer and Lemeshow Goodness-of-Fit Test 
## 
## Call:
## glm(formula = newchd ~ age_cat + sex + chol_cat + sbp_q + smoke, 
##     family = binomial, data = datos)
##  ChiSquare df P_value
##       11.3  8   0.185
HL_DF[,1:4]
##                cut total obs   exp
## 1  [0.0444,0.0718]   163 157 153.5
## 2   (0.0718,0.103]   137 128 124.2
## 3     (0.103,0.12]   135 119 119.2
## 4     (0.12,0.152]   110  94  94.9
## 5    (0.152,0.174]   147 113 123.1
## 6    (0.174,0.206]   130 108 105.4
## 7    (0.206,0.236]   153 112 118.5
## 8    (0.236,0.281]   146 113 107.0
## 9    (0.281,0.358]   121  80  81.2
## 10   (0.358,0.597]   120  70  67.0
HL_DF %>% knitr::kable()
cut total obs exp chi indice
[0.0444,0.0718] 163 157 153.5 0.279 1
(0.0718,0.103] 137 128 124.2 0.342 2
(0.103,0.12] 135 119 119.2 -0.018 3
(0.12,0.152] 110 94 94.9 -0.088 4
(0.152,0.174] 147 113 123.1 -0.913 5
(0.174,0.206] 130 108 105.4 0.256 6
(0.206,0.236] 153 112 118.5 -0.596 7
(0.236,0.281] 146 113 107.0 0.578 8
(0.281,0.358] 121 80 81.2 -0.131 9
(0.358,0.597] 120 70 67.0 0.363 10
ggplot(HL_DF, aes(x=indice))+geom_point(aes(y=exp, color = "blue"))+geom_point(aes(y=obs, color = "red"))+scale_color_manual(labels = c("Esperado", "Observado"), values = c("blue", "red"))

#slope e intercept
pred_modeloagecat7 <- pred_modeloagecat7[-1258,]# esta parece ser una observacion duplicada

a <- valProbggplot(pred_modeloagecat7$pred, pred_modeloagecat7$verdad)
a$ggPlot

#a$Calibration
#curva roc. AUC 0.69
roc<-roc(pred_modeloagecat7$verdad, pred_modeloagecat7$pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc, main = "Curva ROC", col = "red", lwd = 3)

text(0.8, 0.2, paste("AUC ROC =", round(roc$auc, 3)), adj = c(0, 1))

#punto de corte 40%
pred_modeloagecat7 <- pred_modeloagecat7 %>% dplyr::mutate(prediccion = ifelse(pred > 0.4, 1,0))

tab <- tab_xtab(pred_modeloagecat7$prediccion, pred_modeloagecat7$verdad ,tdcol.row = "black",tdcol.n = "black", title = "Matriz de Confusion - cutoff =30") ## ojo si usan el RStudio con fondo blanco cambien el color de los numeros o no van a ver nada...

knitr::asis_output(tab$knitr)
Matriz de Confusion - cutoff =30
prediccion verdad Total
0 1
0 1047 234 1281
1 47 34 81
Total 1094 268 1362
χ2=25.614 · df=1 · &phi=0.141 · p=0.000
cat("Sensibilidad:", paste0(round(caret::sensitivity(factor(pred_modeloagecat7$verdad), factor(pred_modeloagecat7$prediccion)),2)))
## Sensibilidad: 0.82
cat("Especificidad:", paste0(round(caret::specificity(factor(pred_modeloagecat7$verdad), factor(pred_modeloagecat7$prediccion)),2)))
## Especificidad: 0.42
SENS=roc$sensitivities
ESPEC=roc$specificities
CUTOFF=roc$thresholds

d <- cbind(SENS,ESPEC,CUTOFF)
d%>% knitr::kable()
SENS ESPEC CUTOFF
1.000 0.000 -Inf
0.996 0.038 0.048
0.993 0.063 0.057
0.978 0.102 0.065
0.978 0.116 0.070
0.978 0.144 0.072
0.978 0.149 0.076
0.978 0.153 0.081
0.978 0.157 0.085
0.974 0.192 0.091
0.966 0.213 0.098
0.955 0.226 0.100
0.948 0.234 0.102
0.944 0.261 0.106
0.944 0.270 0.110
0.940 0.277 0.112
0.940 0.284 0.116
0.903 0.331 0.119
0.884 0.369 0.123
0.884 0.374 0.128
0.862 0.400 0.134
0.851 0.415 0.137
0.851 0.423 0.139
0.840 0.441 0.145
0.840 0.449 0.150
0.825 0.455 0.152
0.810 0.472 0.155
0.802 0.479 0.158
0.765 0.505 0.161
0.750 0.509 0.162
0.716 0.537 0.167
0.716 0.540 0.173
0.698 0.559 0.176
0.672 0.585 0.180
0.668 0.590 0.184
0.664 0.599 0.185
0.653 0.607 0.188
0.642 0.627 0.195
0.627 0.634 0.201
0.627 0.644 0.204
0.616 0.657 0.208
0.608 0.663 0.212
0.560 0.698 0.215
0.556 0.702 0.221
0.552 0.704 0.230
0.534 0.724 0.233
0.522 0.724 0.235
0.463 0.760 0.239
0.451 0.764 0.244
0.440 0.773 0.249
0.429 0.782 0.256
0.429 0.785 0.263
0.414 0.796 0.265
0.396 0.809 0.268
0.384 0.824 0.273
0.381 0.834 0.278
0.340 0.863 0.281
0.340 0.865 0.289
0.336 0.866 0.299
0.299 0.885 0.303
0.295 0.887 0.309
0.291 0.890 0.319
0.265 0.896 0.330
0.261 0.899 0.339
0.231 0.910 0.344
0.213 0.924 0.353
0.187 0.936 0.359
0.183 0.937 0.368
0.168 0.939 0.379
0.138 0.952 0.383
0.134 0.952 0.390
0.127 0.957 0.408
0.112 0.958 0.426
0.075 0.973 0.439
0.071 0.978 0.458
0.034 0.987 0.470
0.030 0.988 0.477
0.026 0.991 0.496
0.022 0.996 0.535
0.015 0.999 0.579
0.000 1.000 Inf