generated from opensafely/research-template
/
02b_models_coxme.R
67 lines (50 loc) · 2.44 KB
/
02b_models_coxme.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
######################################
# This script:
# - imports processed data
# - fits a mixed effects cox model using the coxme package to a subset of the data
# - saves model summaries (tables and figures)
######################################
# Preliminaries ----
## Import libraries
library('tidyverse')
library('lubridate')
library('survival')
library('coxme')
library('gtsummary')
library('gt')
#library('ehahelper')
## Create output directory
dir.create(here::here("output", "models", "testing"), showWarnings = FALSE, recursive=TRUE)
## Import processed data
data_tte <- read_rds(here::here("output", "data", "data_modelling.rds"))
## Converts logical to integer so that model coefficients print nicely in gtsummary methods
data_cox <- data_tte %>%
mutate(
across(
where(is.logical),
~.x*1L
)
)
# Subset data
## Registered patients counts
practice_counts <- data_cox %>%
group_by(practice_id_latest_active_registration) %>%
summarise(`Number_of_registered_patients` = n())
## Exclude
# data_cox_reduced <- data_cox %>%
# filter(#practice_id_latest_active_registration %in%
# # subset(practice_counts, Number_of_registered_patients >= 100)$practice_id_latest_active_registration,
# practice_id_latest_active_registration %in%
# subset(practice_counts, Number_of_registered_patients >= 100)$practice_id_latest_active_registration[1:1000]) %>%
# droplevels()
# MODELS ----
# Mixed effects Cox model - adjusted; baseline demographics, comorbs, geographical, flu, shielding & pracice as random effect
mod.coxme.adj <- coxme(Surv(follow_up_time, covid_vax) ~
ageband + sex + ethnicity + morbid_obesity +
chronic_heart_disease + diabetes + chronic_kidney_disease_diagnostic + chronic_kidney_disease_all_stages +
chronic_kidney_disease_all_stages_3_5 + sev_mental_ill + learning_disability + chronic_neuro_dis_inc_sig_learn_dis +
asplenia + chronic_liver_disease + chronis_respiratory_disease + immunosuppression_diagnosis +
immunosuppression_medication + imd + rural_urban + prior_covid + flu_vaccine + shielded + shielded_since_feb_15 +
rural_urban + region + (1 | practice_id_latest_active_registration),
data = data_cox)
write_rds(mod.coxme.adj, here::here("output", "models", "testing", "mod_test_coxme_adj.rds"), compress="gz")