generated from opensafely/research-template
/
Table_postopcovid_tv.R
176 lines (137 loc) · 11.7 KB
/
Table_postopcovid_tv.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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
library(foreach)
library(data.table)
ncores <- parallel::detectCores(logical = T)
data.table::setDTthreads(ncores)
source(here::here("analysis","Utils.R"))
###########################################################
dt.tv.splits <- data.table::setDT(feather::read_feather(here::here("output","cohort_postdisch_week_splits.feather")))
procedures <- c('Abdominal','Cardiac','Obstetrics','Orthopaedic','Thoracic', 'Vascular')
n.type.events <- 1:2 #sort(unique(dt.tv[(postop.covid.cohort) ,event]))[-1]
data.table::setkey(dt.tv.splits,patient_id,tstart,tstop)
dt.tv.splits[event == 3, event := 2]
data.table::setkey(dt.tv.splits, patient_id, end.fu, start)
post.op.covid.model.split <-
lapply(n.type.events, function(i) survival::coxph(survival::Surv(start,end,event==i) ~ Abdominal+ Cardiac + Obstetrics + Orthopaedic + Thoracic + Vascular + age.cat + sex + bmi.cat + imd5 + wave + vaccination.status.factor + region + Current.Cancer + Emergency*week.post.disch + Charl12 + recentCOVID + previousCOVID, id = patient_id,data = dt.tv.splits[(postop.covid.cohort) & start <=end], model = T))
data.table::fwrite(broom::tidy(post.op.covid.model.split[[1]], exponentiate= T, conf.int = T), file = here::here("output","postopcovidmodelsplit.csv"))
newdata.rows <- length(levels(dt.tv.splits$week.post.disch)) - 1
newdata.pred <- data.table::data.table('start' = c(-7,0,7,14,21),
'end' = c(0,7,14,21,28),
'event' = rep(F,newdata.rows),
'week.post.disch' = paste(0:(newdata.rows - 1)),
'patient_id' = 1:newdata.rows,
'Abdominal' = rep(T,newdata.rows),
'Cardiac'= rep(F,newdata.rows),
'Obstetrics'=rep(F,newdata.rows),
'Orthopaedic'=rep(F,newdata.rows),
'Thoracic'=rep(F,newdata.rows),
'Vascular'=rep(F,newdata.rows),
'age.cat' = rep('(50,70]',newdata.rows),
'sex' = rep('F',newdata.rows),
'bmi.cat' = rep(levels(dt.tv.splits$bmi.cat)[2],newdata.rows),
'imd5' = rep(levels(dt.tv.splits$imd5)[3], newdata.rows),
'wave' = rep(paste0('Wave_',4),times = newdata.rows),
'vaccination.status.factor' = rep('3',newdata.rows),
'region' = rep("East Midlands",newdata.rows),
'Current.Cancer' = rep(T,newdata.rows),
'Emergency' = rep(F,newdata.rows),
'Charl12' = rep('Single',newdata.rows),
'recentCOVID' = rep(F,newdata.rows),
'previousCOVID' = rep(F,newdata.rows)
)
#newdata.pred[,(procedures) := lapply(procedures, function(x) x == procedures[which.max(dt.tv[,lapply(.SD,sum,na.rm = T), .SDcols = c(procedures)])])]
# newdata.pred[,(covariates[-c(1:length(procedures))]) := lapply(((length(procedures)+1):length(covariates)), function(i.c) {
# if(is.factor(dt.tv[!is.na(get(covariates[i.c])),get(covariates[i.c])])) {
# as.character(rep(max.category(i.c),newdata.rows))
# } else if(is.logical(dt.tv[!is.na(get(covariates[i.c])),get(covariates[i.c])])) {
# as.logical(rep(max.category(i.c),newdata.rows))
# } else if(is.numeric(dt.tv[!is.na(get(covariates[i.c])),get(covariates[i.c])])) {
# is.numeric(rep(max.category(i.c),newdata.rows))
# } else {
# rep(max.category(i.c),newdata.rows)
# }
# })]
times <- lapply(n.type.events, function(i) { survival::basehaz(post.op.covid.model.split[[i]],centered = F)$time })
base.haz <- lapply(n.type.events, function(i) { data.table::data.table('time' = times[[i]],
'base.haz' = survival::basehaz(post.op.covid.model.split[[i]],centered = F)[,1] -
c(0,head(survival::basehaz(post.op.covid.model.split[[i]],centered = F)[,1],-1)))})
lp <- lapply(n.type.events, function(i) { data.table::data.table('time' = seq(-7,21,7),
'risk' = exp(predict(object = post.op.covid.model.split[[i]],
type = 'lp',
newdata = newdata.pred))) })
base.haz.merge <- Reduce(x =base.haz,f = function(x,y) merge(x,y,by = 'time', no.dups = T, suffixes = c(".x",".y"), all = T, sort = T))
weekly.post.op.risk <-
unlist(round(100*apply(exp(apply(safelog(1 - Reduce('+',lapply(n.type.events, function(i) {
lp[[i]][base.haz.merge[order(time),.SD,.SDcols = c(1,i+1)],,roll =Inf,on = 'time', rollends = c(T,T)][time >= -7][order(time),.(.SD[,1]*.SD[,2]),.SDcols = c(2:3)]
}))),2,cumsum))*
lp[[1]][base.haz.merge[order(time),.SD,.SDcols = c(1,2)],,roll =Inf,on = 'time', rollends = c(T,T)][time >= -7][order(time),.(.SD[,1]*.SD[,2]),.SDcols = c(2:3)] ,2,cumsum), digits = 3))
weekly.post.op.risk[!is.finite(weekly.post.op.risk)] <- 0
times.comb <- unique(sort(unlist(times)))[unique(sort(unlist(times))) >= -7]
weekly.post.op.risk <- c(weekly.post.op.risk[max(which(times.comb <= 0))],
weekly.post.op.risk[max(which(times.comb <= 7))],
weekly.post.op.risk[max(which(times.comb <= 14))],
weekly.post.op.risk[max(which(times.comb <= 21))],
weekly.post.op.risk[max(which(times.comb <= 28))])
weekly.post.op.risk <- data.table::data.table("Risk" = weekly.post.op.risk - c(0,weekly.post.op.risk[-length(weekly.post.op.risk)]),
"Risk period" = c("Week pre discharge","1st week","2nd week","3rd week","4th week"))
##############
post.op.VTE.model.split <-
lapply(n.type.events, function(i) survival::coxph(survival::Surv(start,end,event.VTE==i) ~ Abdominal + survival::strata(postcovid) +
Cardiac + Obstetrics + Orthopaedic + Thoracic + Vascular + age.cat +
sex + bmi.cat + imd5 + wave + vaccination.status.factor + region +
Current.Cancer + Emergency*week.post.disch + Charl12 + recentCOVID + previousCOVID,
id = patient_id,
data = dt.tv.splits[(postcovid.VTE.cohort) & start <=end], model = T))
data.table::fwrite(broom::tidy(post.op.VTE.model.split[[1]], exponentiate= T, conf.int = T), file = here::here("output","postopVTEmodelsplit.csv"))
newdata.rows <- length(levels(dt.tv.splits$week.post.disch)) - 1
newdata.pred <- data.table::data.table('start' = rep(c(-7,0,7,14,21), times = 2),
'end' = rep(c(0,7,14,21,28),times = 2),
'event.VTE' = rep(F,newdata.rows*2),
'week.post.disch' = rep(paste(0:(newdata.rows - 1)), times = 2),
'patient_id' = rep(1:2,each = newdata.rows),
'Abdominal' = rep(T,newdata.rows*2),
'Cardiac'= rep(F,newdata.rows*2),
'Obstetrics'=rep(F,newdata.rows*2),
'Orthopaedic'=rep(F,newdata.rows*2),
'Thoracic'=rep(F,newdata.rows*2),
'postcovid'=rep(c(0,1), each = newdata.rows),
'Vascular'=rep(F,newdata.rows*2),
'age.cat' = rep('(50,70]',newdata.rows*2),
'sex' = rep('F',newdata.rows*2),
'bmi.cat' = rep(levels(dt.tv.splits$bmi.cat)[2],newdata.rows*2),
'imd5' = rep(levels(dt.tv.splits$imd5)[3], newdata.rows*2),
'wave' = rep(paste0('Wave_',4),times = newdata.rows*2),
'vaccination.status.factor' = rep('3',newdata.rows*2),
'region' = rep("East Midlands",newdata.rows*2),
'Current.Cancer' = rep(T,newdata.rows*2),
'Emergency' = rep(F,newdata.rows*2),
'Charl12' = rep('Single',newdata.rows*2),
'recentCOVID' = rep(F,newdata.rows*2),
'previousCOVID' = rep(F,newdata.rows*2)
)
times <- lapply(n.type.events, function(i) { survival::basehaz(post.op.VTE.model.split[[i]],centered = F)$time })
base.haz <- lapply(n.type.events, function(i) { data.table::data.table('time' = times[[i]],
'base.haz' = survival::basehaz(post.op.VTE.model.split[[i]],centered = F)[,1] -
c(0,head(survival::basehaz(post.op.VTE.model.split[[i]],centered = F)[,1],-1)))})
lp <- lapply(n.type.events, function(i) { data.table::dcast(data.table::data.table('patient_id' = rep(1:2, each = 5),
'time' = rep(seq(-7,21,7),2),
'risk' = exp(predict(object = post.op.VTE.model.split[[i]],
type = 'lp',
newdata = newdata.pred))),time ~patient_id, value.var = 'risk')})
base.haz.merge <- Reduce(x =base.haz,f = function(x,y) merge(x,y,by = 'time', no.dups = T, suffixes = c(".x",".y"), all = T, sort = T))
weekly.post.op.VTE.risk <-
unlist(round(100*apply(exp(apply(safelog(1 - Reduce('+',lapply(n.type.events, function(i) {
lp[[i]][base.haz.merge[order(time),.SD,.SDcols = c(1,i+1)],,roll =Inf,on = 'time', rollends = c(T,T)][time >= -7][order(time),.(.SD[,1]*.SD[,3], .SD[,2]*.SD[,3]),.SDcols = c(2:4)]
}))),2,cumsum))*
lp[[1]][base.haz.merge[order(time),.SD,.SDcols = c(1,2)],,roll =Inf,on = 'time', rollends = c(T,T)][time >= -7][order(time),.(.SD[,1]*.SD[,3], .SD[,2]*.SD[,3]),.SDcols = c(2:4)], 2, cumsum ), digits = 3))
weekly.post.op.VTE.risk[!is.finite(weekly.post.op.VTE.risk)] <- 0
times.comb <- unique(sort(unlist(times)))[unique(sort(unlist(times))) >= -7]
weekly.post.op.VTE.risk <- rbind(weekly.post.op.VTE.risk[max(which(times.comb <= 0)),],
weekly.post.op.VTE.risk[max(which(times.comb <= 7)),],
weekly.post.op.VTE.risk[max(which(times.comb <= 14)),],
weekly.post.op.VTE.risk[max(which(times.comb <= 21)),],
weekly.post.op.VTE.risk[max(which(times.comb <= 28)),])
weekly.post.op.VTE.risk <- data.table::data.table("COVID"= rep(c(F,T), each = 5),
"Risk" = as.vector(weekly.post.op.VTE.risk - rbind(c(0,0),weekly.post.op.VTE.risk[-nrow(weekly.post.op.VTE.risk),])),
"Risk period" = c("Week pre discharge","1st week","2nd week","3rd week","4th week"))
##################################
save(weekly.post.op.risk,weekly.post.op.VTE.risk, file = here::here("output","postopcovid_tv.RData"))