-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathobesity_figure3-1.R
115 lines (98 loc) · 8.4 KB
/
obesity_figure3-1.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
107
108
109
110
111
112
113
114
115
library(DatabaseConnector)
library(ggplot2)
library(dplyr)
library(interactions)
library(ggpubr)
library(lme4)
library(tibble)
library(sjPlot)
library(nlme)
library(emmeans)
#################################
connectionDetails <- DatabaseConnector::createConnectionDetails(dbms = "sql server",
server = '',
user = '',
password = '',
pathToDriver = 'C:/jdbc_driver')
conn <- DatabaseConnector::connect(connectionDetails)
# #################################
person <- dbGetQuery(conn,"SELECT person_id, birth_datetime, gender_source_value FROM [CDMPv535_ABMI].[dbo].[person]")
StratPop_l1_s1_p1_t3316_c3317_s1_o2580 <- readRDS("~/output/cmOutput/StratPop_l1_s1_p1_t3316_c3317_s1_o2580.rds")
table(StratPop_l1_s1_p1_t3316_c3317_s1_o2580$treatment)
metaDF <- merge(StratPop_l1_s1_p1_t3316_c3317_s1_o2580, connected_db) %>% select(person_id, treatment, cohortStartDate, daysToCohortEnd)
colnames(metaDF) <- c('person_id', 'treatment', 'cohortStartDate', 'cohortEndDate')
metaDF <- metaDF %>% mutate(cohortEndDate = cohortStartDate + cohortEndDate)
target_cohort <- metaDF %>% filter(treatment == 1)
comparator_cohort <- metaDF %>% filter(treatment == 0)
#################################
# Total Eosinophils, Total Neutrophils 경향성을 보기 (+ 외래 측정치만 활용) (+ 나이, 성별, 반복측정치 보정)
### Eosinophil count (cells/mcL) : 3006504 (0, 1000) / Neutrophil count (cells/mcL) : 3018010 (0, 10000)
extractcode <- paste("select m.person_id, m.measurement_date, m.value_as_number from [CDMPv535_ABMI].[dbo].[measurement] m
join (select * from [CDMPv535_ABMI].[dbo].[visit_occurrence] where visit_concept_id in (9202)) opv
on m.visit_occurrence_id = opv.visit_occurrence_id
where m.measurement_concept_id = 3010813 and m.person_id in (", paste(target_cohort$person_id, collapse = ', '), ')')
target_leukocyte <- dbGetQuery(conn, extractcode)
target_leukocyte <- merge(target_leukocyte, person, by = 'person_id')
colnames(target_leukocyte) <- c('person_id', 'measurement_date', 'leukocyte', 'birth', 'gender')
target_leukocyte <- merge(target_cohort, target_leukocyte, by = 'person_id')
target_leukocyte <- target_leukocyte %>% filter(measurement_date >= cohortStartDate) %>% filter(cohortEndDate >= measurement_date) %>% mutate(followup_years = as.numeric(measurement_date - cohortStartDate)/365) %>% mutate(age = as.numeric(measurement_date - as.Date(substr(as.character(birth),1,10)))/365) %>% mutate(obesity = 'yes')
extractcode <- paste("select m.person_id, m.measurement_date, m.value_as_number from [CDMPv535_ABMI].[dbo].[measurement] m
join (select * from [CDMPv535_ABMI].[dbo].[visit_occurrence] where visit_concept_id in (9202)) opv
on m.visit_occurrence_id = opv.visit_occurrence_id
where m.measurement_concept_id = 3006504 and m.person_id in (", paste(target_cohort$person_id, collapse = ', '), ')')
target_eosinophil_percent <- dbGetQuery(conn, extractcode)
target_eosinophil_percent <- merge(target_eosinophil_percent, person, by = 'person_id')
colnames(target_eosinophil_percent) <- c('person_id', 'measurement_date', 'eosinophil_percent', 'birth', 'gender')
target_eosinophil_percent <- merge(target_cohort, target_eosinophil_percent, by = 'person_id')
target_eosinophil_percent <- target_eosinophil_percent %>% filter(measurement_date >= cohortStartDate) %>% filter(cohortEndDate >= measurement_date) %>% mutate(followup_years = as.numeric(measurement_date - cohortStartDate)/365) %>% mutate(age = as.numeric(measurement_date - as.Date(substr(as.character(birth),1,10)))/365) %>% mutate(obesity = 'yes')
target_eosinophil <- merge(target_leukocyte, target_eosinophil_percent %>% select(person_id, measurement_date, eosinophil_percent), by=c('person_id', 'measurement_date'))
target_eosinophil <- target_eosinophil %>% mutate(total_eosinophil_count = leukocyte * eosinophil_percent * 10)
extractcode <- paste("select m.person_id, m.measurement_date, m.value_as_number from [CDMPv535_ABMI].[dbo].[measurement] m
join (select * from [CDMPv535_ABMI].[dbo].[visit_occurrence] where visit_concept_id in (9202)) opv
on m.visit_occurrence_id = opv.visit_occurrence_id
where m.measurement_concept_id = 3010813 and m.person_id in (", paste(comparator_cohort$person_id, collapse = ', '), ')')
comparator_leukocyte <- dbGetQuery(conn, extractcode)
comparator_leukocyte <- merge(comparator_leukocyte, person, by = 'person_id')
colnames(comparator_leukocyte) <- c('person_id', 'measurement_date', 'leukocyte', 'birth', 'gender')
comparator_leukocyte <- merge(comparator_cohort, comparator_leukocyte, by = 'person_id')
comparator_leukocyte <- comparator_leukocyte %>% filter(measurement_date >= cohortStartDate) %>% filter(cohortEndDate >= measurement_date) %>% mutate(followup_years = as.numeric(measurement_date - cohortStartDate)/365) %>% mutate(age = as.numeric(measurement_date - as.Date(substr(as.character(birth),1,10)))/365) %>% mutate(obesity = 'no')
extractcode <- paste("select m.person_id, m.measurement_date, m.value_as_number from [CDMPv535_ABMI].[dbo].[measurement] m
join (select * from [CDMPv535_ABMI].[dbo].[visit_occurrence] where visit_concept_id in (9202)) opv
on m.visit_occurrence_id = opv.visit_occurrence_id
where m.measurement_concept_id = 3006504 and m.person_id in (", paste(comparator_cohort$person_id, collapse = ', '), ')')
comparator_eosinophil_percent <- dbGetQuery(conn, extractcode)
comparator_eosinophil_percent <- merge(comparator_eosinophil_percent, person, by = 'person_id')
colnames(comparator_eosinophil_percent) <- c('person_id', 'measurement_date', 'eosinophil_percent', 'birth', 'gender')
comparator_eosinophil_percent <- merge(comparator_cohort, comparator_eosinophil_percent, by = 'person_id')
comparator_eosinophil_percent <- comparator_eosinophil_percent %>% filter(measurement_date >= cohortStartDate) %>% filter(cohortEndDate >= measurement_date) %>% mutate(followup_years = as.numeric(measurement_date - cohortStartDate)/365) %>% mutate(age = as.numeric(measurement_date - as.Date(substr(as.character(birth),1,10)))/365) %>% mutate(obesity = 'no')
comparator_eosinophil <- merge(comparator_leukocyte, comparator_eosinophil_percent %>% select(person_id, measurement_date, eosinophil_percent), by=c('person_id', 'measurement_date'))
comparator_eosinophil <- comparator_eosinophil %>% mutate(total_eosinophil_count = leukocyte * eosinophil_percent * 10)
target_eosinophil <- target_eosinophil %>% filter(followup_years<=10)
comparator_eosinophil <- comparator_eosinophil %>% filter(followup_years<=10)
target_eosinophil <- target_eosinophil %>% filter(total_eosinophil_count != 0)
comparator_eosinophil <- comparator_eosinophil %>% filter(total_eosinophil_count != 0)
total <- rbind(target_eosinophil, comparator_eosinophil)
lmm_model <- lmer(
total_eosinophil_count ~ followup_years*obesity + age + gender
+ (1 | person_id),
data = total
)
#################################
predicted <- emmeans(lmm_model, specs = ~ obesity, at = list(followup_years = 10), pbkrtest.limit = 5000, lmerTest.limit = 5000)
conf_int <- confint(predicted)
print(conf_int)
print(pairs(predicted))
#################################
sjPlot::plot_model(lmm_model, type = "pred", terms = c("followup_years", 'obesity'), axis.title = c('Years of follow-up', 'Eosinophil count (cells/mcL)'), title = '', alpha = 0) +
theme(panel.background = element_blank(), axis.line = element_line(colour = "grey50")) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .7, position = position_dodge(width = 0.2), linetype = 1) +
scale_linetype_manual(values = c("dashed", "solid")) +
scale_fill_manual(values=c("white", "black")) +
scale_color_manual(values=c("black", "black")) +
#scale_y_continuous(limits = c(20, 32)) +
aes(linetype=group, color=group) +
geom_point(size = 4, stroke = .5, shape = 21, position = position_dodge(width = 0.05)) +
scale_x_discrete(limits=c(0, 2, 4, 6, 8, 10)) # 0, 1, 2, 3, 4, 5 ## 0, 2, 4, 6, 8, 10
sjPlot::tab_model(lmm_model)
# nrow(unique(total %>% filter(treatment == 1) %>% select(person_id)))
# nrow(unique(total %>% filter(treatment == 0) %>% select(person_id)))