-
Notifications
You must be signed in to change notification settings - Fork 2
/
02-epi_analitica.R
106 lines (77 loc) · 2.91 KB
/
02-epi_analitica.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
# paquetes ----------------------------------------------------------------
library(mosaicData)
library(tidyverse)
library(skimr)
library(rlang)
library(naniar)
library(compareGroups)
library(janitor)
library(epiR)
library(epiDisplay)
library(broom)
library(writexl)
library(avallecam)
# importar ----------------------------------------------------------------
data("Whickham")
smoke <- Whickham %>% as_tibble()
smoke %>% miss_var_summary()
smoke %>% skim()
# limpieza ----------------------------------------------------------------
smoke_clean <- read_rds("data/smokeclean_20190906.rds")
# tabla 1 y 2 -------------------------------------------------------------
tab1 <- smoke_clean %>%
compareGroups(outcome~smoker+agegrp+age,data = .,byrow = T) %>%
createTable(show.all = T,sd.type = 2)
tab1
#exportar en XLSX
tab1 %>%
export2xls("table/tabla02.xls")
# medidas de asociación ---------------------------------------------------
#epiR
smoke_tab1 <- with(smoke_clean,table(smoker_2,outcome_2)) %>% print()
epi.2by2(smoke_tab1,method = "cohort.count")
#epiDisplay
smoke_tab2 <- with(smoke_clean,table(outcome,smoker)) %>% print()
cs(cctable = smoke_tab2)
# evaluar confusión -------------------------------------------------------
#Mantel-Haenszel
smoke_tab3 <- with(smoke_clean,table(smoker_2,outcome_2,agegrp)) %>% print()
epi.2by2(smoke_tab3,method = "cohort.count")
mhor(mhtable=smoke_tab3,graph = F,design = "cohort")
# modelo ------------------------------------------------------------------
smoke_clean %>% glimpse()
#simple
wm1 <- glm(outcome_1 ~ smoker,
data = smoke_clean,
family = poisson(link = "log"))
epi_tidymodel_rr(wm1)
#multiple: controlar por confusión
wm1 <- glm(outcome_1 ~ smoker + age,
data = smoke_clean,
family = poisson(link = "log"))
epi_tidymodel_rr(wm1)
#exportar
epi_tidymodel_rr(wm1) %>%
write_xlsx("table/tabla03.xlsx")
# calcular error estandar robusto
# cov.m1 <- sandwich::vcovHC(wm1, type="HC0")
# std.err <- sqrt(diag(cov.m1))
# r.est <- cbind(Estimate= exp(coef(wm1)), "Robust SE" = exp(std.err),
# "Pr(>|z|)" = 2 * pnorm(abs(coef(wm1)/std.err), lower.tail=FALSE),
# LL = exp(coef(wm1)) - 1.96 * exp(std.err),
# UL = exp(coef(wm1)) + 1.96 * exp(std.err))
# r.est %>% round(digits = 2)
# time to event -----------------------------------------------------------
library(survival)
library(survminer)
url <- "http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/valung.csv"
valung <- read_csv(url) %>%
mutate(dead=as.factor(dead),
dead_1=as.numeric(dead)-1)
valung %>% skim()
valung %>% miss_var_summary()
fit <- survfit(Surv(t, dead_1) ~ therapy, data = valung)
fit %>% broom::tidy() %>% print(n=Inf)
ggsurvplot(fit, data = valung, risk.table = TRUE,
conf.int = T,ggtheme = theme_bw(),#ggtheme = theme_minimal(),
pval = T,censor=FALSE,tables.theme = theme_cleantable())