You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
library(tidyverse)
library(readxl)
library(writexl)
library(tableone)
library(haven) # Read sas files
library(here)
library(kableExtra)
library(jtools) # for summ() and plot_summs
library(sjPlot) # for tab_model
library(ggplot2) # survival/TTE analyses and other graphs
library(ggsurvfit) # survival/TTE analyses
library(survival) # survival/TTE analyses
library(gtsummary) # survival/TTE analyses
library(ggfortify) # autoplot
library(tidycmprsk) # competing risk analysis
library(ordinal) # clinstatus ordinal regression
library(logistf) # Firth regression in case of rare events
library(finalfit) # missing data exploration
library(mice) # multiple imputation
library(jomo) # multiple imputation
library(mitools) # multiple imputation
Load Data
Reshape dataset and define ITT population
#reshapedf<-df %>% ## corresponds to N in publication
pivot_wider(id_cols= c("USUBJID","ARM","ACTARM","BCSOSN","d29dthe0","d29dthe1","ttrecov0","ttrecov1","or15scor","agec","SEX","BDURSYMP","COMORB2","REGION1","REGION2","BMI","RACE","ETHNIC","STRATAR","STRATAV","HYPTFL","CADFL","CHFFL","CRDFL","CORQFL","CLDFL","CKDFL","DIAB1FL","DIAB2FL","OBESIFL","CANCFL","IMMDFL","ASTHMAFL","CVDFL","COAGFL","NICFL","RFDVTFL","SURGFL","IMMOFL","CPDVTFL","COMORB1"),
names_from="adyc",
values_from="ordscor")
#reorderdf<-df %>%
select("USUBJID","ARM","ACTARM","BCSOSN","Baseline","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","36","37","45", everything())
#define ITT -> keep all as randomized# addmargins(table(df$ARM, df$ACTARM))
Baseline Characteristics
df$trial<- c("ACTT2")
df$JAKi<- c("Baricitinib")
df<-df %>% ## no missing data // no randdate // ARM == 1 includes remdesivir!
rename(id_pat=USUBJID,
sex=SEX,
age=agec,
region=REGION2
)
df<-df %>%
mutate(trt= case_when(ARM=="Baricitinib + Remdesivir"~1,
TRUE~0))
# fix ages ">89" using the original datasetdf_original<-df_original %>%
rename(id_pat=USUBJID)
elderly<- unique(df[which(df$age==">89"), "id_pat"])
for (iin1:dim(elderly)[[1]]){
df[which(df$id_pat== paste(elderly[i,])), "age"] <- paste(df_original[which(df_original$id_pat== paste(elderly[i,])), "AGE"])
}
df$age<- as.numeric(df$age)
df %>%
drop_na(age) %>%
ggplot(aes(x=age)) +
geom_density(fill="blue", color="black") +
labs(title="Density Plot of Age",
x="Age",
y="Density")
# Ethnicitydf<-df %>%
mutate(ethn= case_when(RACE=="UNKNOWN"ÐNIC=="HISPANIC OR LATINO"~"HISPANIC OR LATINO",
TRUE~ c(RACE)))
# Country# For Country only Region (Asia, Europe, North America) and Non-US Site (n=148) and US Site (n=885) available. According to publication: "There were 67 trial sites in 8 countries: the United States (55 sites), Singapore (4), South Korea (2), Mexico (2), Japan (1), Spain (1), the United Kingdom (1), and Denmark (1)." => add region to country and reclassify laterdf$country<-df$region# Days with symptoms prior to randomizationdf$sympdur<-df$BDURSYMP## 16 missing datadf %>%
drop_na(sympdur) %>%
ggplot(aes(x=sympdur)) +
geom_density(fill="blue", color="black") +
labs(title="Density Plot of Symptom Duration",
x="Symptom Duration",
y="Density")
# Severity of COVID-19 with respect to respiratory support at randomisation / from publication: "Scores on the ordinal scale are as follows: 1, not hospitalized, no limitations of activities; 2, not hospitalized, limitation of activities, home oxygen requirement, or both; 3, hospitalized, not requiring supplemental oxygen and no longer requiring ongoing medical care (used if hospitalization was extended for infection-control reasons); 4, hospitalized, not requiring supplemental oxygen but requiring ongoing medical care (Covid-19–related or other medical conditions); 5, hospitalized, requiring any supplemental oxygen; 6, hospitalized, receiving noninvasive ventilation or use of high-flow oxygen devices; 7, hospitalized, receiving invasive mechanical ventilation or ECMO; and 8, death. Five deaths (three in patients receiving baricitinib plus RDV and two in patients receiving placebo plus RDV) occurred within the day 15 visit window but after 14 days — these deaths are included in the outcome of the ordinal score at day 15 but not in the outcome of mortality over the first 14 days."
# transform all clinical scoresscore_transform<-function(df, clinstatus_var, score_var) {
df<-df %>%
mutate({{ clinstatus_var }} :=
case_when({{ score_var }} %in% c(1, 2, 3) ~1,
{{ score_var }} ==4~2,
{{ score_var }} ==5~3,
{{ score_var }} ==6~4,
{{ score_var }} ==7~5,
{{ score_var }} ==8~6)) %>%
mutate({{ clinstatus_var }} :=factor({{ clinstatus_var }}, levels=1:6))
}
df<- score_transform(df, clinstatus_baseline, Baseline)
addmargins(table(df$clinstatus_baseline, df$trt, useNA="always")) # no missing / # corresponds to publication
# Vaccinationdf$vacc<-0# trial happened before roll-out of any vaccine# Viremia // Viral load value <LOQ and/or undectectabledf_vl_baseline<-df_vl %>%
filter(`Study Visit`=="Day 1") %>%
rename(id_pat=`Subject ID`)
df_vl_baseline<-df_vl_baseline %>%
mutate(vl_baseline= case_when(`SARS-CoV-2 PCR Qual. Result`=="Negative"~1,
`SARS-CoV-2 PCR Qual. Result`=="Positive"~0))
df<- left_join(df, df_vl_baseline[, c("vl_baseline", "id_pat")], by= join_by(id_pat==id_pat)) ## merge# table(df$vl_baseline, useNA = "always")# Variant // not available# Serology // not available
Discussion points BASELINE data:
Endpoints
# time to event datadf$death_d<- as.numeric(df$d29dthe0) # everyone with a value here, dieddf$discharge_d<- as.numeric(df$ttrecov0) # everyone with a value here, was discharged/reached discharge criteria // CAVE: 2 were discharged at day 9 but died later at home (within 28 days)# (i) Primary outcome: Mortality at day 28df<-df %>% # 47 have no outcome data, reason for censoring unclear (see below); a few have info in daily score
mutate(mort_28= case_when(death_d<29~1,
discharge_d<29~0, # discharge to hospice did not occurclinstatus_28%in% c(2,3,4,5) ~0, # still at hospital but alive, clinstatus variables after 28d are all empty (see View below)discharge_d>28~0)) # discharged later, proof of still alive (but empty)# table(df$mort_28, df$trt, useNA = "always") # deaths correspond to publication (but they have -1 based on different counting of day 0) // However, the reasons for the 21 missing in int and the 26 in cont are hard to identify from the publication# table(df$ACTARM, df$mort_28, useNA = "always") # the 16 not treated were censored at day 0 or 1 => NA for death_d/discharge_d/clinstatus_28 // The remaining 14 in intervention (see Fig 1 in publication): 8 Withdrew, 1 was withdrawn by investigator, 1 became ineligible after enrollment, 2 Had severe adverse event or adverse event other than death (?!), 2 Had other reason (?!) // The remaining 17 in control, I guess (see Fig 1): 16 Withdrew, 2 Were withdrawn by investigator, 1 Became ineligible after enrollment, 1 Had severe adverse event or adverse event other than death, 1 Was transferred to another hospital, 3 Had other reason.# df %>%# select(mort_28, ACTARM, death_d, discharge_d, d29dthe1, ttrecov1, clinstatus_baseline, clinstatus_1, clinstatus_2, clinstatus_3, clinstatus_4, clinstatus_5, clinstatus_6, clinstatus_7, clinstatus_8, clinstatus_9, clinstatus_10, clinstatus_11, clinstatus_12, clinstatus_13, clinstatus_14, clinstatus_15,# clinstatus_15_imp,# clinstatus_16, clinstatus_17, clinstatus_18, clinstatus_19, clinstatus_20, clinstatus_21, clinstatus_22, clinstatus_23, clinstatus_24, clinstatus_25, clinstatus_26, clinstatus_27, clinstatus_28, clinstatus_29, clinstatus_30, clinstatus_31, clinstatus_32, clinstatus_33, clinstatus_34, clinstatus_36, clinstatus_37, clinstatus_45) %>%# filter(is.na(mort_28)) %>%# View()# First, keep mort_28 as complete case# Second, use multiple imputation (see below)# Third, apply a deterministic imputation (see notes): we use the same rules as ACTT2 => we can use clinstatus_15_impdf<- score_transform(df, clinstatus_15_imp, df$or15scor)
df<-df %>%
mutate(mort_28_dimp= case_when(is.na(mort_28) &clinstatus_15_imp%in% c("1","2","3","4","5") ~0,
TRUE~ c(mort_28)))
# (ii) Mortality at day 60df$mort_60<-df$mort_28# max fup time in ACTT-2 was 28 days; thus mort_60 imputed from mort_28 // stick to complete case# (iii) Time to death within max. follow-up timedf<-df %>% # don't bother about the NA in mort_28 since they are not dead, i.e. just censor (==0) everyone else
mutate(death_reached= case_when(death_d<29~1,
TRUE~0))
df<-df %>% # no missing and those that were discharged and afterwards died have correct time to event data
mutate(death_time= case_when(death_d>=0~ c(death_d), # time to death, if no time to death, then...d29dthe1>=0~ c(d29dthe1))) # time to death censoring data# df %>% # those died have correct time to event data. Those discharged before d28 (and not withdrawn/ltfu/hospiz/SAE) should receive max fup time? But some do not, they are censored but reason for censoring unclear.# select(death_reached, ACTARM, death_time, death_d, discharge_d, d29dthe1, ttrecov1, clinstatus_baseline, clinstatus_1, clinstatus_2, clinstatus_3, clinstatus_4, clinstatus_5, clinstatus_6, clinstatus_7, clinstatus_8, clinstatus_9, clinstatus_10, clinstatus_11, clinstatus_12, clinstatus_13, clinstatus_14, clinstatus_15, or15scor, clinstatus_16, clinstatus_17, clinstatus_18, clinstatus_19, clinstatus_20, clinstatus_21, clinstatus_22, clinstatus_23, clinstatus_24, clinstatus_25, clinstatus_26, clinstatus_27, clinstatus_28, clinstatus_29, clinstatus_30, clinstatus_31, clinstatus_32, clinstatus_33, clinstatus_34, clinstatus_36, clinstatus_37, clinstatus_45) %>%# filter(death_reached != 1) %>%# filter(discharge_d < 29 & death_time < 28) %>%# View()# (iv) New mechanical ventilation among survivors within 28 days.df<-df %>% # see View below
mutate(new_mv_28= case_when((clinstatus_baseline==2|clinstatus_baseline==3|clinstatus_baseline==4)
& (mort_28==0| is.na(mort_28))
& (clinstatus_1==5|clinstatus_2==5|clinstatus_3==5|clinstatus_4==5|clinstatus_5==5|clinstatus_6==5|clinstatus_7==5|clinstatus_8==5|clinstatus_9==5|clinstatus_10==5|clinstatus_11==5|clinstatus_12==5|clinstatus_13==5|clinstatus_14==5|clinstatus_15==5|clinstatus_15_imp==5|clinstatus_16==5|clinstatus_17==5|clinstatus_18==5|clinstatus_19==5|clinstatus_20==5|clinstatus_21==5|clinstatus_22==5|clinstatus_23==5|clinstatus_24==5|clinstatus_25==5|clinstatus_26==5|clinstatus_27==5|clinstatus_28==5)
~1,
(clinstatus_baseline==2|clinstatus_baseline==3|clinstatus_baseline==4)
&mort_28==0~0))
# df %>%# select(new_mv_28, mort_28, death_reached, death_time, death_d, discharge_d, d29dthe1, ttrecov1, clinstatus_baseline, clinstatus_1, clinstatus_2, clinstatus_3, clinstatus_4, clinstatus_5, clinstatus_6, clinstatus_7, clinstatus_8, clinstatus_9, clinstatus_10, clinstatus_11, clinstatus_12, clinstatus_13, clinstatus_14, clinstatus_15, clinstatus_15_imp, clinstatus_16, clinstatus_17, clinstatus_18, clinstatus_19, clinstatus_20, clinstatus_21, clinstatus_22, clinstatus_23, clinstatus_24, clinstatus_25, clinstatus_26, clinstatus_27, clinstatus_28, clinstatus_29, clinstatus_30, clinstatus_31, clinstatus_32, clinstatus_33, clinstatus_34, clinstatus_36, clinstatus_37, clinstatus_45) %>%# filter(is.na(new_mv_28)) %>%# # filter(clinstatus_baseline != 5) %>%# # filter(death_reached == 0) %>%# View()# (iv) Alternative definition/analysis: New mechanical ventilation OR death within 28 days => include all in denominator. df<-df %>%
mutate(new_mvd_28= case_when(new_mv_28==1|mort_28==1~1,
new_mv_28==0|mort_28==0~0))
# (v) Clinical status at day 28df<-df %>% # Adapt clinstatus_28, since currently NA for those discharged or died. Then, check if nothing in 28d window (there was nothing), then LVCF including use or15scor
mutate(clinstatus_28= case_when(clinstatus_28==5~5,
clinstatus_28==4~4,
clinstatus_28==3~3,
clinstatus_28==2~2,
mort_28==1~6, # died within 28dmort_28==0~1)) # discharged alive / reached discharge criteria within 28ddf$clinstatus_28<-factor(df$clinstatus_28, levels=1:6)
# table(df$clinstatus_28, useNA = "always") # again, the 47 missing## Imputation according to protocol: If there was daily data for the ordinal score available but with missing data for single days, then we carried last observed value forward unless for day 28, whereby we first considered data from the window (+/-3 days) but there was nothing in that window => LVCFdfcs<-df %>%
select(id_pat, clinstatus_baseline, clinstatus_1, clinstatus_2, clinstatus_3, clinstatus_4, clinstatus_5, clinstatus_6, clinstatus_7, clinstatus_8, clinstatus_9, clinstatus_10, clinstatus_11, clinstatus_12, clinstatus_13, clinstatus_14, clinstatus_15_imp, clinstatus_16, clinstatus_17, clinstatus_18, clinstatus_19, clinstatus_20, clinstatus_21, clinstatus_22, clinstatus_23, clinstatus_24, clinstatus_25, clinstatus_26, clinstatus_27, clinstatus_28)
impute_last_forw=function(df){
first= which(names(df)%in%c("clinstatus_baseline"))
last= which(names(df)%in%c("clinstatus_28"))
for (iin1:dim(df)[[1]]){
for (jinfirst[1]:last[1]){
p=df[i, j]
df[i,j] <-
ifelse(!is.na(df[i, j]), p, df[i, j-1])
}
}
df
}
dfcs<- impute_last_forw(dfcs)
dfcs<-dfcs %>% # To control, don't overwrite
rename(clinstatus_28_imp=clinstatus_28)
df<- left_join(df, dfcs[, c("clinstatus_28_imp", "id_pat")], by= join_by(id_pat==id_pat)) # Merge imputed variable back# table(df$clinstatus_28, useNA = "always")# table(df$clinstatus_28_imp, useNA = "always") # All imputed (i.e. carried forward from imputed day 15)# (vi) Time to discharge or reaching discharge criteria up to day 28df<-df %>%
mutate(discharge_reached= case_when(discharge_d<29~1,
TRUE~0))
df<-df %>% # 2 are left without any time to event data => impute max. follow-up time
mutate(discharge_time= case_when(discharge_d>=0~ c(discharge_d), # time to discharge, if no time to discharge, then...ttrecov1>=0~ c(ttrecov1)+1)) # time to discharge censoring time (see ACTT-2 documentation regarding time zero -> now matching d29dthe1)df<-df %>% # restrict to max fup time 28d
mutate(discharge_time= case_when(discharge_time>28~28,
TRUE~discharge_time))
df<-df %>% # add 28d for those that died // Patients who died prior to day 28 are assumed not having reached discharge, i.e. counted as 28 days (as someone who has been censored on day 28).
mutate(discharge_time_sens= case_when(mort_28==1~28,
TRUE~discharge_time))
# table(df$discharge_time, df$discharge_reached, useNA = "always")# (vi) Sens-analysis: Alternative definition/analysis of outcome: time to sustained discharge within 28 days. There are no re-admissions documented in ACTT-2df$discharge_reached_sus<-df$discharge_reacheddf$discharge_time_sus<-df$discharge_time# (vii) Viral clearance up to day 5, day 10, and day 15 (Viral load value <LOQ and/or undectectable)df_vl<-df_vl %>%
rename(visit=`Actual Study Day`,
id_pat=`Subject ID`,
result=`SARS-CoV-2 PCR Qual. Result`)
df_vl<-df_vl %>%
mutate(visit= case_when(visit==-1~1,
TRUE~visit))
df_vl<-df_vl %>%
arrange(id_pat, visit)
df_vl<-df_vl %>%
mutate(vir_clear_5= case_when(visit==5&result=="Negative"~1,
visit==5&result=="Positive"~0,
visit==4&result=="Negative"~1,
visit==4&result=="Positive"~0,
visit==3&result=="Negative"~1,
visit==3&result=="Positive"~0,
visit==2&result=="Negative"~1,
visit==2&result=="Positive"~0,
visit==1&result=="Negative"~1,
visit==1&result=="Positive"~0))
df_vl<-df_vl %>%
mutate(vir_clear_10= case_when(visit==10&result=="Negative"~1,
visit==10&result=="Positive"~0,
visit==9&result=="Negative"~1,
visit==9&result=="Positive"~0,
visit==8&result=="Negative"~1,
visit==8&result=="Positive"~0,
visit==7&result=="Negative"~1,
visit==7&result=="Positive"~0,
visit==6&result=="Negative"~1,
visit==6&result=="Positive"~0,
vir_clear_5==1~1,
vir_clear_5==0~0))
df_vl<-df_vl %>%
mutate(vir_clear_15= case_when(visit==15&result=="Negative"~1,
visit==15&result=="Positive"~0,
visit==14&result=="Negative"~1,
visit==14&result=="Positive"~0,
visit==13&result=="Negative"~1,
visit==13&result=="Positive"~0,
visit==12&result=="Negative"~1,
visit==12&result=="Positive"~0,
visit==11&result=="Negative"~1,
visit==11&result=="Positive"~0,
vir_clear_10==1~1,
vir_clear_10==0~0))
df_vl5<-df_vl %>%
filter(!is.na(vir_clear_5)) %>%
group_by(id_pat) %>%
arrange(abs(visit-5)) %>%
filter(abs(visit-5) == min(abs(visit-5))) %>%
ungroup()
df<- left_join(df, df_vl5[, c("vir_clear_5", "id_pat")], by= join_by(id_pat==id_pat)) ## mergedf_vl10<-df_vl %>%
filter(!is.na(vir_clear_10)) %>%
group_by(id_pat) %>%
arrange(abs(visit-10)) %>%
filter(abs(visit-10) == min(abs(visit-10))) %>%
ungroup()
df<- left_join(df, df_vl10[, c("vir_clear_10", "id_pat")], by= join_by(id_pat==id_pat)) ## mergedf_vl15<-df_vl %>%
filter(!is.na(vir_clear_15)) %>%
group_by(id_pat) %>%
arrange(abs(visit-15)) %>%
filter(abs(visit-15) == min(abs(visit-15))) %>%
ungroup()
df<- left_join(df, df_vl15[, c("vir_clear_15", "id_pat")], by= join_by(id_pat==id_pat)) ## merge# (viii) Quality of life at day 28: Not available in ACTT-2# (ix) Participants with an adverse event grade 3 or 4, or a serious adverse event, excluding death, by day 28df_ae_tot<-df_ae_tot %>%
rename(id_pat=`Subject ID`)
df_ae_tot$`Study Day`<- as.numeric(df_ae_tot$`Study Day`)
## Warning: NAs introduced by coercion
df_ae<-df_ae_tot %>%
filter(Severity%in% c("Grade 3", "Grade 4", "Grade 5") | (Severity=="Grade 2"&`Serious Adverse Event`=="Yes")) %>%
filter(`Study Day`<29| is.na(`Study Day`))
# Keep just 1 id_pat (-> ANY adverse event grade 3 (severe), 4 (serious)) df_ae34_unique<-df_ae %>% distinct(id_pat, .keep_all=TRUE)
# Assign the outcomedf_ae34_unique$ae_28<-1# table(df_ae34_unique$ae_28)# mergedf<- left_join(df, df_ae34_unique[, c("ae_28", "id_pat")], by= join_by(id_pat==id_pat)) ## merge variable to main df# the remaining missing have no AE grade 34 -> recode as 0 and exclude deathsdf<-df %>%
mutate(ae_28= case_when(is.na(ae_28) ~0, # the LTFU were dischargedmort_28==1~NA, # exclude the deathsTRUE~ae_28))
# table(df$ae_28, df$mort_28, useNA = "always")# addmargins(table(df$ae_28, df$trt, useNA = "always"))# (ix) Sens-analysis: Alternative definition/analysis of outcome: incidence rate ratio (Poisson regression) -> AE per person by d28ae_npp<-df_ae %>%
group_by(id_pat)%>%
summarise(ae_28_sev= n())
df<- left_join(df, ae_npp[, c("ae_28_sev", "id_pat")], by= join_by(id_pat==id_pat)) # merge variable to main df# the remaining missing have no AE grade 34 -> recode as 0 and exclude deathsdf<-df %>%
mutate(ae_28_sev= case_when(is.na(ae_28_sev) ~0, # the LTFU were dischargedmort_28==1~NA, # exclude the deathsTRUE~ae_28_sev))
# addmargins(table(df$ae_28_sev, df$trt, useNA = "always"))# (ix) Sens-analysis: Alternative definition/analysis of outcome: time to first (of these) adverse event, within 28 days, considering death as a competing risk (=> censor and set to 28 days): Not available# (x) Adverse events of special interest within 28 days: a) thromboembolic events (venous thromboembolism, pulmonary embolism, arterial thrombosis), b) secondary infections (bacterial pneumonia including ventilator-associated pneumonia, meningitis and encephalitis, endocarditis and bacteremia, invasive fungal infection including pulmonary aspergillosis), c) Reactivation of chronic infection including tuberculosis, herpes simplex, cytomegalovirus, herpes zoster and hepatitis B, d) serious cardiovascular and cardiac events (including stroke and myocardial infarction), e) events related to signs of bone marrow suppression (anemia, lymphocytopenia, thrombocytopenia, pancytopenia), f) malignancy, g) gastrointestinal perforation (incl. gastrointestinal bleeding/diverticulitis), h) liver dysfunction/hepatotoxicity (grade 3 and 4)df_ae_tot<-df_ae_tot %>%
rename(ae=`MedDRA Preferred Term`,
ae_class=`MedDRA System Organ Class`,
ae_desc=`Adverse Event`) %>%
filter(`Study Day`<29| is.na(`Study Day`))
df_ae_tot<- left_join(df_ae_tot, df[, c("trt", "id_pat")], by= join_by(id_pat==id_pat))
df_thrombo<-df_ae_tot %>% # a) thromboembolic events (venous thromboembolism, pulmonary embolism, arterial thrombosis)
filter(ae%in% c("Axillary vein thrombosis", "Brachiocephalic vein thrombosis", "Deep vein thrombosis", "Embolism venous", "Peripheral artery occlusion", "Peripheral ischaemia", "Thrombosis", "Pulmonary embolism")) %>%
mutate(aesi="thrombo")
df_sec_inf<-df_ae_tot %>% # b) secondary infections (bacterial pneumonia including ventilator-associated pneumonia, meningitis and encephalitis, endocarditis and bacteremia, invasive fungal infection including pulmonary aspergillosis), but not COVID-19 pneumonia!
filter(ae_class%in% c("Infections and infestations") &!ae=="Septic shock") %>%
mutate(aesi="sec_inf")
# df_reactivate <- NA# df_ae_tot %>% # filter(grepl("hepatitis b|zoster|herpes|cytome|tuber|tb", ae, ignore.case = TRUE)) %>%# View()df_cardiac<-df_ae_tot %>% # d) serious cardiovascular and cardiac events (including stroke and myocardial infarction) (excl. hypertension)
filter(ae_class%in% c("Cardiac disorders") |ae=="Cerebrovascular accident") %>%
mutate(aesi="cardiac")
df_penia<-df_ae_tot %>% # e) events related to signs of bone marrow suppression (anemia, lymphocytopenia, thrombocytopenia, pancytopenia)
filter(grepl("penia|anemia|anaemia", ae, ignore.case=TRUE)) %>%
mutate(aesi="penia")
# df_malig <- NA# df_malig <- df_ae_tot %>% # f) malignancy# filter(ae_class %in% c("Neoplasms benign, malignant and unspecified (incl cysts and polyps)")) %>% # mutate(aesi = "malig")df_git_bl<-df_ae_tot %>% # g) gastrointestinal perforation (incl. gastrointestinal bleeding/diverticulitis)
filter(ae_class%in% c("Hepatobiliary disorders","Gastrointestinal disorders") & grepl("hemor|haemor|bleed", ae, ignore.case=TRUE)) %>%
mutate(aesi="git_bl")
df_hepatox<-df_ae_tot %>% # h) liver dysfunction/hepatotoxicity (grade 3 and 4)
filter(ae%in% c("Hyperbilirubinaemia", "Hypertransaminasaemia", "Liver injury")) %>%
mutate(aesi="hepatox")
df_mods<-df_ae_tot %>% # i) Multiple organ dysfunction syndrome and septic shock
filter(ae%in% c("Multiple organ dysfunction syndrome", "Septic shock", "Shock")) %>%
mutate(aesi="mods")
df_aesi<- rbind(df_mods, df_hepatox, df_git_bl, df_penia, df_cardiac, df_sec_inf, df_thrombo)
df_aesi<-df_aesi %>%
select(id_pat, trt, aesi, ae, ae_desc, ae_class)
# table(df_aesi$trt, df_aesi$aesi)# double-check if there are any duplicate AEs within the same person and if it is the same event or distinct onesdf_aesi<-df_aesi %>%
group_by(id_pat) %>%
mutate(duplicate_id= duplicated(ae) &!is.na(ae)) %>%
ungroup()
df_aesi<-df_aesi %>%
filter(duplicate_id==F)
# Save
saveRDS(df_aesi, file="df_aesi_actt2.RData")
# (xi) Adverse events, any grade and serious adverse event, within 28 days, grouped by organ classesdf_ae<-df_ae_tot %>%
select(id_pat, trt, ae, ae_desc, ae_class)
# double-check if there are any duplicate AEs within the same person and if it is the same event or distinct onesdf_ae<-df_ae %>%
group_by(id_pat) %>%
mutate(duplicate_id= duplicated(ae) &!is.na(ae)) %>%
ungroup()
df_ae<-df_ae %>%
filter(duplicate_id==F)
# Save
saveRDS(df_ae, file="df_ae_actt2.RData")
Discussion points
Define final datasets
# the overall setdf_all<-df# reduce the df set to our standardized set across all trialsdf<-df %>%
select(id_pat, trt, sex, age, trial, JAKi,
ethn,
country,
# icu, sympdur,
vacc,
clinstatus_baseline, vbaseline,
comed_dexa,
comed_rdv,
comed_toci, comed_ab, comed_acoa, comed_interferon, comed_other,
comed_cat,
comorb_lung, comorb_liver, comorb_cvd, comorb_aht, comorb_dm, comorb_obese, comorb_smoker, immunosupp,
any_comorb, comorb_cat, comorb_any, comorb_count,
comorb_autoimm, comorb_cancer, comorb_kidney,
crp,
# sero, vl_baseline,
# variant,mort_28, mort_28_dimp,
mort_60, death_reached, death_time,
new_mv_28, new_mvd_28,
clinstatus_28_imp,
discharge_reached, discharge_time, discharge_time_sens, discharge_reached_sus, discharge_time_sus,
ae_28, ae_28_sev,
vir_clear_5, vir_clear_10, vir_clear_15
)
# export for one-stage model, i.e., add missing variables df_os<-dfdf_os$icu<-NAdf_os$sero<-NAdf_os$variant<-NA# Save
saveRDS(df_os, file="df_os_actt2.RData")
Missing data plot: One-stage dataset
# Bar plot, missing data, each data point, standardized one-stage datasetoriginal_order<- colnames(df_os)
missing_plot<-df_os %>%
summarise_all(~ mean(is.na(.))) %>%
gather() %>%
mutate(key=factor(key, levels=original_order)) %>%
ggplot(aes(x=key, y=value)) +
geom_bar(stat="identity") +
labs(x="Columns", y="Proportion of Missing Values", title="Missing Data - standardized one-stage dataset") +
theme(axis.text.x= element_text(angle=45, hjust=1)) +
ylim(0, 1)
print(missing_plot)
### Baseline table, by individuals with no missing data vs any missing data (or only in mort_28)# df_core <- df_core %>% mutate(complete = ifelse(rowSums(is.na(.)) > 0, 0, 1));table(df_core$complete) # ANY missing df_core$resp<-ifelse(is.na(df_core$mort_28), 0, 1);table(df_core$resp) # only mort_28 missing
# Third, let's explore if the variables from my substantive model plus auxiliary variables are associated with mort_28mort28.aux<- glm(mort_28~trt+age+clinstatus_baseline+sex+ethn+country+sympdur+comorb_cat+comed_dexa+comed_ab+comed_acoa+comed_other+crp# + vl_baseline, # + ae_28_sev
,family="binomial"
,data=df_imp)
summary(mort28.aux)
##
## Call:
## glm(formula = mort_28 ~ trt + age + clinstatus_baseline + sex +
## ethn + country + sympdur + comorb_cat + comed_dexa + comed_ab +
## comed_acoa + comed_other + crp, family = "binomial", data = df_imp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.52980 -0.34836 -0.18870 -0.07287 3.02379
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -4.126e+01 3.249e+03 -0.013
## trt -4.999e-01 3.110e-01 -1.607
## age 5.774e-02 1.197e-02 4.824
## clinstatus_baseline3 1.586e+01 8.347e+02 0.019
## clinstatus_baseline4 1.713e+01 8.347e+02 0.021
## clinstatus_baseline5 1.812e+01 8.347e+02 0.022
## sexM 2.113e-01 3.208e-01 0.659
## ethnASIAN -9.239e-01 1.548e+00 -0.597
## ethnBLACK OR AFRICAN AMERICAN -1.185e+00 1.426e+00 -0.831
## ethnHISPANIC OR LATINO -1.043e+00 1.406e+00 -0.742
## ethnNATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER -1.017e-01 1.773e+00 -0.057
## ethnUNKNOWN 4.045e-01 1.678e+00 0.241
## ethnWHITE -8.258e-01 1.385e+00 -0.596
## countryEurope -1.559e+01 2.706e+03 -0.006
## countryNorth America 1.166e+00 1.312e+00 0.889
## sympdur -2.387e-02 3.315e-02 -0.720
## comorb_cat2 1.061e+00 8.001e-01 1.327
## comorb_cat3 1.343e+00 7.660e-01 1.754
## comorb_cat4 1.587e+00 1.143e+00 1.389
## comed_dexa1 -3.298e-01 6.712e-01 -0.491
## comed_ab1 4.476e-02 3.442e-01 0.130
## comed_acoa1 1.332e+00 1.091e+00 1.221
## comed_other1 1.611e+01 3.139e+03 0.005
## crp -5.835e-04 8.983e-04 -0.650
## Pr(>|z|)
## (Intercept) 0.9899
## trt 0.1080
## age 1.41e-06 ***
## clinstatus_baseline3 0.9848
## clinstatus_baseline4 0.9836
## clinstatus_baseline5 0.9827
## sexM 0.5101
## ethnASIAN 0.5507
## ethnBLACK OR AFRICAN AMERICAN 0.4060
## ethnHISPANIC OR LATINO 0.4579
## ethnNATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER 0.9543
## ethnUNKNOWN 0.8095
## ethnWHITE 0.5509
## countryEurope 0.9954
## countryNorth America 0.3739
## sympdur 0.4714
## comorb_cat2 0.1847
## comorb_cat3 0.0795 .
## comorb_cat4 0.1649
## comed_dexa1 0.6232
## comed_ab1 0.8965
## comed_acoa1 0.2222
## comed_other1 0.9959
## crp 0.5160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 436.60 on 948 degrees of freedom
## Residual deviance: 329.43 on 925 degrees of freedom
## (84 observations deleted due to missingness)
## AIC: 377.43
##
## Number of Fisher Scoring iterations: 18
# Fourth, let's explore if they are associated with missingness of mort_28:# df_imp$resp<-ifelse(is.na(df_imp$mort_28), 0, 1);table(df_imp$resp) # only mort_28 missing # mort28.aux.resp <- glm(resp ~ trt# + age # + clinstatus_baseline# + sex# + ethn# + region# + sympdur# + comorb_cat# ,family="binomial"# ,data=df_imp)# summary(mort28.aux.resp)df_imp %>%
missing_compare(dependent, explanatory) %>%
knitr::kable(row.names=FALSE, align= c("l", "l", "r", "r", "r"))
Missing data analysis: mort_28
Not missing
Missing
p
age
Mean (SD)
55.2 (15.6)
59.2 (16.3)
0.087
clinstatus_baseline
2
138 (97.2)
4 (2.8)
0.221
3
542 (96.1)
22 (3.9)
4
203 (94.0)
13 (6.0)
5
103 (92.8)
8 (7.2)
sex
F
363 (95.3)
18 (4.7)
0.959
M
623 (95.6)
29 (4.4)
ethn
AMERICAN INDIAN OR ALASKA NATIVE
10 (100.0)
0 (0.0)
0.792
ASIAN
98 (97.0)
3 (3.0)
BLACK OR AFRICAN AMERICAN
151 (96.8)
5 (3.2)
HISPANIC OR LATINO
234 (95.1)
12 (4.9)
NATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER
11 (100.0)
0 (0.0)
UNKNOWN
12 (92.3)
1 (7.7)
WHITE
470 (94.8)
26 (5.2)
country
Asia
66 (98.5)
1 (1.5)
0.328
Europe
13 (100.0)
0 (0.0)
North America
907 (95.2)
46 (4.8)
Duration of Symptoms (days)
Mean (SD)
8.5 (4.5)
8.3 (4.7)
0.821
comorb_cat
1
153 (90.5)
16 (9.5)
0.008
2
280 (97.2)
8 (2.8)
3
524 (96.0)
22 (4.0)
4
29 (96.7)
1 (3.3)
comed_dexa
0
945 (95.4)
46 (4.6)
0.756
1
41 (97.6)
1 (2.4)
comed_ab
0
384 (92.5)
31 (7.5)
<0.001
1
602 (97.4)
16 (2.6)
comed_acoa
0
62 (81.6)
14 (18.4)
<0.001
1
924 (96.6)
33 (3.4)
comed_other
0
8 (40.0)
12 (60.0)
<0.001
1
978 (96.5)
35 (3.5)
crp
Mean (SD)
143.2 (158.5)
166.5 (83.5)
0.382
vl_baseline
0
504 (95.8)
22 (4.2)
0.911
1
161 (96.4)
6 (3.6)
ae_28_sev
Mean (SD)
1.1 (2.3)
1.6 (2.4)
0.110
# Fifth, check age
summary(df_imp$age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18.00 43.00 56.00 55.39 67.00 101.00
hist(df_imp$age, breaks=50) # looks fine
# Sixth, check sympdur
summary(df_imp$sympdur)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 5.000 8.000 8.459 10.000 35.000 16
df_imp$sqcrptrunc=sqrt(df_imp$crptrunc)
hist(df_imp$sqcrptrunc) # looks fine
### Reshape to long format# names(df_imp)# str(df_imp)df_imp<-df_imp %>% # rename to differentiate between baseline clinical status and follow-up clinical statuses
rename(clinicalstatus_baseline=clinstatus_baseline)
# reshapedf_imp_long<-df_imp %>%
pivot_longer(cols= starts_with("clinstatus"), names_to="time", values_to="clinstatus")
# names(df_imp_long)# Convert time to numericdf_imp_long$time<- as.numeric(gsub("clinstatus_", "", df_imp_long$time))
# class(df_imp_long$time)# summary(df_imp_long$time)# have clinstatus as numericdf_imp_long$clinstatus_n<- as.numeric(df_imp_long$clinstatus)
### We will impute separately by treatment arm, since we have to expect an effect modification between outcome x trt over timedf_imp_long_int<-df_imp_long %>%
filter(trt==1)
df_imp_long_cont<-df_imp_long %>%
filter(trt==0)
## Explore distribution of clinical status over time, by armplot_clinstat_int<- ggplot(df_imp_long_int, aes(x=time, y=clinstatus_n)) +
stat_summary(fun="mean", geom="point", color="blue", size=3, na.rm=TRUE) +
geom_smooth(method="loess", se=FALSE, color="red", na.rm=TRUE) +
labs(title="Distribution of Clinical Status Over Time / Intervention",
x="Time",
y="Mean Clinical Status") +
theme_minimal()
print(plot_clinstat_int)
plot_clinstat_cont<- ggplot(df_imp_long_cont, aes(x=time, y=clinstatus_n)) +
stat_summary(fun="mean", geom="point", color="blue", size=3, na.rm=TRUE) +
geom_smooth(method="loess", se=FALSE, color="red", na.rm=TRUE) +
labs(title="Distribution of Clinical Status Over Time / Control",
x="Time",
y="Mean Clinical Status") +
theme_minimal()
print(plot_clinstat_cont)
Multiple imputation
#### INTERVENTION group## jomo only accepts numeric or factors, check and adapt# str(df_imp_long_int)df_imp_long_int$timesq<- sqrt(df_imp_long_int$time) # see X below
attach(df_imp_long_int)
Y2<-data.frame(mort_28# level 2 variables (baseline patient characteristics)
, age
, sex
, ethn
, country
, sqsympdur
, comorb_cat
, sqcrptrunc
, ae_28_sev
, vl_baseline
, comed_dexa
, comed_ab
, comed_other
)
Y<-data.frame(clinstatus_n) # level 1 variable within clustering variableX<- cbind(1, data.frame(clinicalstatus_baseline, time, timesq)) # matrix modelling linearity of clinstatus throughout day 28clus<-data.frame(id_pat) # clustering variable (patient)Z<-data.frame(rep(1,dim(df_imp_long_int)[1]),df_imp_long_int[,c("time")]) # random intercept and random slope
colnames(Z) <- c("const", "time")
nimp<-30# set number of iterations## run jomo# dry runimputed_int_mcmc<-jomo.MCMCchain(Y=Y, Y2=Y2, X=X, Z=Z, clus=clus, nburn=2)
# plot(c(1:2),imputed_int_mcmc$collectbeta[1,1,1:2],type="l")# plot(c(1:2),imputed_int_mcmc$collectcovu[5,5,1:2],type="l")
set.seed(1569)
imputed_int<- jomo(Y=Y, Y2=Y2, X=X, Z=Z, clus=clus, nburn=1000, nbetween=1000, nimp=nimp)
# nburn<-1000# imputed_int_mcmc<-jomo.MCMCchain(Y=Y, Y2=Y2, X=X, Z=Z, clus=clus, nburn=nburn)# plot(c(1:nburn),imputed_int_mcmc$collectbeta[1,1,1:nburn],type="l")# plot(c(1:nburn),imputed_int_mcmc$collectcovu[5,5,1:nburn],type="l")# convert to jomo object, split imputations, and exclude original data (imputation "0")imp.list_int<- imputationList(split(imputed_int, imputed_int$Imputation)[-1])
# checks
round(prop.table(table(imp.list_int[[1]]$`1`$mort_28, useNA="always"))*100,1) # first imputed dataset
round(prop.table(table(imp.list_int[[1]]$`2`$mort_28, useNA="always"))*100,1) # second imputed dataset
round(prop.table(table(df_imp_long_int$mort_28, useNA="always"))*100,1) # original data
summary(imp.list_int[[1]]$`1`$comorb_cat)
summary(imp.list_int[[1]]$`2`$sqsympdur)
summary(imp.list_int[[1]]$`2`$vl_baseline)
#### CONTROL group## jomo only accepts numeric or factors, check and adapt# str(df_imp_long_cont)df_imp_long_cont$timesq<- sqrt(df_imp_long_cont$time) # see X below
attach(df_imp_long_cont)
Y2<-data.frame(mort_28# level 2 variables (baseline patient characteristics)
, age
, sex
, ethn
, country
, sqsympdur
, comorb_cat
, sqcrptrunc
, ae_28_sev
, vl_baseline
, comed_dexa
, comed_ab
, comed_other
)
Y<-data.frame(clinstatus_n) # level 1 variable within clustering variableX<- cbind(1, data.frame(clinicalstatus_baseline, time, timesq)) # matrix modelling linearity of clinstatus throughout day 28clus<-data.frame(id_pat) # clustering variable (patient)Z<-data.frame(rep(1,dim(df_imp_long_cont)[1]),df_imp_long_cont[,c("time")]) # random intercept and random slope
colnames(Z)<-c("const", "time")
nimp<-30# set number of iterations# run jomo
set.seed(1569)
imputed_cont<- jomo(Y=Y, Y2=Y2, X=X, Z=Z, clus=clus, nburn=1000, nbetween=1000, nimp=nimp)
# nburn<-1000# imputed_cont_mcmc<-jomo.MCMCchain(Y=Y, Y2=Y2, X=X, Z=Z, clus=clus, nburn=nburn)# plot(c(1:nburn),imputed_cont_mcmc$collectbeta[1,1,1:nburn],type="l")# plot(c(1:nburn),imputed_cont_mcmc$collectcovu[5,5,1:nburn],type="l")# convert to jomo object, split imputations, and exclude original data (imputation "0")imp.list_cont<- imputationList(split(imputed_cont, imputed_cont$Imputation)[-1])
# checks
round(prop.table(table(imp.list_cont[[1]]$`1`$mort_28, useNA="always"))*100,1) # first imputed dataset
round(prop.table(table(imp.list_cont[[1]]$`2`$mort_28, useNA="always"))*100,1) # second imputed dataset
round(prop.table(table(df_imp_long_cont$mort_28, useNA="always"))*100,1) # original data
summary(imp.list_cont[[1]]$`1`$comorb_cat)
summary(imp.list_cont[[1]]$`2`$sqsympdur)
#### Add trt back, change from long to wide format, and finally combine the two data framesimputed_int$trt<-1imputed_int_s<-imputed_int %>% # remove imputation variables, not needed anymore
select(trt, age, sqsympdur, mort_28, sex, ethn, country, comed_dexa, comed_ab, comed_other, comorb_cat, sqcrptrunc, vl_baseline, clinicalstatus_baseline, ae_28_sev, clus, Imputation)
imputed_int_wide<-imputed_int_s %>% # change from long to wide format, i.e. remove duplicates within Imputation sets
group_by(Imputation) %>%
distinct(clus, .keep_all=TRUE) # 515 * 51 = 26265, correctimputed_cont$trt<-0# treatment variableimputed_cont_s<-imputed_cont %>% # remove imputation variables, not needed anymore
select(trt, age, sqsympdur, mort_28, sex, ethn, country, comed_dexa, comed_ab, comed_other, comorb_cat, sqcrptrunc, vl_baseline, clinicalstatus_baseline, ae_28_sev, clus, Imputation)
imputed_cont_wide<-imputed_cont_s %>% # change from long to wide format, i.e. remove duplicates within Imputation sets
group_by(Imputation) %>%
distinct(clus, .keep_all=TRUE) # 518 * 51 = 26418, correctimputed_combined<- rbind(imputed_cont_wide, imputed_int_wide)
#### Convert combined df to jomo object, split imputations, and exclude original data (imputation "0")imp.list<- imputationList(split(imputed_combined, imputed_combined$Imputation)[-1])
### Checks
round(prop.table(table(imp.list[[1]]$`1`$mort_28, imp.list[[1]]$`1`$trt, useNA="always"),2)*100,1) # first imputed dataset
round(prop.table(table(imp.list[[1]]$`2`$mort_28, imp.list[[1]]$`2`$trt, useNA="always"),2)*100,1) # second imputed dataset
round(prop.table(table(imp.list[[1]]$`3`$mort_28, imp.list[[1]]$`3`$trt, useNA="always"),2)*100,1) # second imputed dataset
round(prop.table(table(imp.list[[1]]$`4`$mort_28, imp.list[[1]]$`4`$trt, useNA="always"),2)*100,1) # second imputed dataset
round(prop.table(table(imp.list[[1]]$`5`$mort_28, imp.list[[1]]$`5`$trt, useNA="always"),2)*100,1) # second imputed dataset
round(prop.table(table(df_imp$mort_28, df_imp$trt, useNA="always"),2)*100,1) # original data
summary(imp.list[[1]]$`1`$comorb_cat)
summary(imp.list[[1]]$`2`$sqsympdur)
# Multiple imputation analysis under MAR; use mitools package to fit imputed and combined data list and apply Rubin's rules# mort.28.mi <- imp.list %>%# with(glm(mort_28 ~ trt# + age# + clinicalstatus_baseline# , family = binomial)) %>%# pool() %>%# summary(conf.int = T, exponentiate = T)# mort.28.mi
(i.i) Covariate adjustment for primary endpoint: Mortality at day 28
# unadjusted estimator for the (absolute) risk differencemort.28.prop.test<- prop.test(x= with(df, table(trt, mort_28)))
# print(mort.28.prop.test)# Estimate-diff(mort.28.prop.test$estimate)
# Covariate-Adjusted Analysis# Fit the `glm` object# Same as Complete case analysis, substantive model // but don't use piping, otherwise problem in margins::marginsdf_mort28_comp<-df %>% filter(!is.na(mort_28))
mort.28.cov.adj<-
glm(formula=mort_28~trt+age+clinstatus_baseline,
data=df_mort28_comp,
family= binomial(link="logit")
)
# Print a summary of the `glm` object
summary(mort.28.cov.adj)
##
## Call:
## glm(formula = mort_28 ~ trt + age + clinstatus_baseline, family = binomial(link = "logit"),
## data = df_mort28_comp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3296 -0.3558 -0.2139 -0.1177 2.8297
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -22.63979 862.31663 -0.026 0.979
## trt -0.35087 0.28903 -1.214 0.225
## age 0.05467 0.01068 5.120 3.05e-07 ***
## clinstatus_baseline3 16.03047 862.31632 0.019 0.985
## clinstatus_baseline4 17.22120 862.31631 0.020 0.984
## clinstatus_baseline5 18.34404 862.31632 0.021 0.983
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 457.65 on 985 degrees of freedom
## Residual deviance: 361.06 on 980 degrees of freedom
## AIC: 373.06
##
## Number of Fisher Scoring iterations: 18
# Or, we can obtain the standard error of the estimate two ways. The first way is using the margins::margins() command, using the robust standard errors from sandwich::vcovHC // The second way to obtain these would be the bias corrected and accelerated (BCa) non-parametric bootstrap# You’ll see that we now have a standard error, p-value under the hypothesis that the marginal effect is 0, and a 95% Confidence Interval for the estimate.
library(sandwich)
library(margins)
mort.28.cov.adj.ame<-margins::margins(
model=mort.28.cov.adj,
# Specify treatment variablevariables="trt",
# Convert to outcome scale, not link scaletype="response",
# Obtain robust standard errorsvcov=sandwich::vcovHC(x=mort.28.cov.adj, type="HC3")
)
summary(object=mort.28.cov.adj.ame, level=0.95)
## factor AME SE z p lower upper
## trt -0.0177 0.0147 -1.2022 0.2293 -0.0465 0.0111
# (iv) Alternative definition/analysis: New mechanical ventilation OR death within 28 days => include all in denominator.
table(df$new_mvd_28, df$trt, useNA="always")
# Sub-distribution hazards, i.e., represents the rate per unit of time of the event as well as the influence of competing events.df<-df %>% # cuminc needs a factor variable with censored patients coded as 0, the event as 1 and the competing event as 2.
mutate(discharge_reached_comp= case_when (discharge_reached==0& (mort_28==0| is.na(mort_28)) ~0,
discharge_reached==1& (mort_28==0| is.na(mort_28)) ~1,
mort_28==1~2))
df$discharge_reached_comp<- as.factor(df$discharge_reached_comp)
# Cumulative incidence for competing risks
cuminc(Surv(discharge_time, discharge_reached_comp) ~1, data=df)
# 1: patients without Dexamethasone nor Tocilizumab => JAKi effect alone# 2: patients with Dexamethasone but no Tocilizumab => JAKi effect with Dexa only# 3: patients with Dexamethasone and Tocilizumab => JAKi effect with Dexa + Toci# 4: patients with Tocilizumab but no Dexamethasone (if exist) => JAKi effect with Toci only mort.28.comed<-df %>%
glm(mort_28~trt*comed_cat+age+clinstatus_baseline
, family="binomial", data=.)
summ(mort.28.comed, exp=T, confint=T, model.info=T, model.fit=F, digits=2)
Discussion points
