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)
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()
[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()
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()
[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()
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 |