-
Notifications
You must be signed in to change notification settings - Fork 0
/
obesity_subgroup3-1.R
121 lines (101 loc) · 8.92 KB
/
obesity_subgroup3-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
116
117
118
119
120
121
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_s2_p1_t3316_c3317_s1_o3954 <- readRDS("~/output/cmOutput/StratPop_l1_s2_p1_t3316_c3317_s1_o3954.rds")
Obesity <- StratPop_l1_s2_p1_t3316_c3317_s1_o3954 %>% filter(treatment == 1)
Nonobesity <- StratPop_l1_s2_p1_t3316_c3317_s1_o3954 %>% filter(treatment == 0)
Obesity <- merge(Obesity, connected_db) %>% select(person_id, treatment, cohortStartDate, daysToCohortEnd, daysToEvent)
Nonobesity <- merge(Nonobesity, connected_db) %>% select(person_id, treatment, cohortStartDate, daysToCohortEnd, daysToEvent)
colnames(Obesity) <- c('person_id', 'treatment', 'cohortStartDate', 'cohortEndDate', 'daysToEvent')
colnames(Nonobesity) <- c('person_id', 'treatment', 'cohortStartDate', 'cohortEndDate', 'daysToEvent')
Obesity <- Obesity %>% mutate(cohortEndDate = cohortStartDate + cohortEndDate)
Nonobesity <- Nonobesity %>% mutate(cohortEndDate = cohortStartDate + cohortEndDate)
Obesity_WCT <- Obesity %>% filter(is.na(daysToEvent)) %>% select(person_id, cohortStartDate, cohortEndDate)
Obesity_UCT <- Obesity %>% filter(!is.na(daysToEvent)) %>% select(person_id, cohortStartDate, cohortEndDate)
Nonobesity_WCT <- Nonobesity %>% filter(is.na(daysToEvent)) %>% select(person_id, cohortStartDate, cohortEndDate)
Nonobesity_UCT <- Nonobesity %>% filter(!is.na(daysToEvent)) %>% select(person_id, cohortStartDate, cohortEndDate)
### 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(Obesity_UCT$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(Obesity_UCT, 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 = 3018010 and m.person_id in (", paste(Obesity_UCT$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(Obesity_UCT, 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(Nonobesity_UCT$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(Nonobesity_UCT, 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 = 3018010 and m.person_id in (", paste(Nonobesity_UCT$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(Nonobesity_UCT, 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)))