From 2e4892cf578cfd41dd35d847f65c44173bea1c33 Mon Sep 17 00:00:00 2001
From: twillia2 <thomas.christie.williams@ed.ac.uk>
Date: Thu, 30 Jan 2025 22:54:02 +0000
Subject: [PATCH] Upload New File
---
.../BronchSTOP_analysis_for_submission.qmd | 769 ++++++++++++++++++
1 file changed, 769 insertions(+)
create mode 100644 maternal_VE_study/BronchSTOP_analysis_for_submission.qmd
diff --git a/maternal_VE_study/BronchSTOP_analysis_for_submission.qmd b/maternal_VE_study/BronchSTOP_analysis_for_submission.qmd
new file mode 100644
index 0000000..9d8b5d4
--- /dev/null
+++ b/maternal_VE_study/BronchSTOP_analysis_for_submission.qmd
@@ -0,0 +1,769 @@
+---
+title: "BronchStop data quality report and analysis"
+format: pdf
+---
+
+## Analysis on `r format(Sys.Date(),'%d/%m/%Y')`
+
+```{r setup}
+#| echo: false
+#| message: false
+
+library(tidyverse)
+library(janitor)
+library(tidyr)
+library(glue)
+library(gt)
+library(survival)
+library(broom)
+library(gtsummary)
+
+load(file = "bronchstop_livedata.rdata")
+
+print(glue("Number of consented subjects in database: {data%>%filter(consented==TRUE)%>%nrow()}"))
+
+data<-data%>%
+ filter(consented==TRUE)%>%
+ mutate(vax_gest=case_when(vax_gest %in% c(1:18) ~ vax_gest+23,
+ vax_gest == 99 ~ 99))%>%
+ mutate(vax_gest_record=case_when(vax_gest_record %in% c(1:18) ~ vax_gest_record+23,
+ vax_gest_record == 99 ~ 99))%>%
+ mutate(age_gest=case_when(age_gest %in% c(1:15) ~ age_gest+21,
+ age_gest == 99 ~ 99))%>%
+ mutate(age_gest=replace_na(age_gest,40))%>%
+ mutate(age_gest_term_or_not=if_else(age_gest<36,FALSE,TRUE))%>%
+ mutate(vax_gest_born=age_gest-vax_gest)%>%
+ mutate(potential_excludes=if_else(vax_gest_born %in% c(-1,0,1,2),TRUE,FALSE))%>%
+ mutate(LOA = time_length(interval(ymd(attend_dat),ymd(dis_dat_final)), "day"))%>%
+ mutate(follow_up_form_complete=case_when(follow_up_form_complete==0 ~"incomplete",
+ follow_up_form_complete==1 ~"unverified",
+ follow_up_form_complete==2 ~"complete"))%>%
+ mutate(vax_yn_record=case_when(vax_yn_record==0 ~"confirmed_unvax",
+ vax_yn_record==1 ~"confirmed_vax",
+ vax_yn_record==99 ~"checked_not_known",
+ is.na(vax_yn_record) ~"not_done_yet"))%>%
+ mutate(mum_thinks_vaccinated=case_when(vax_yn==0 ~"thinks_unvax",
+ vax_yn==1 ~"thinks_vax",
+ vax_yn==99 ~"mum_unsure"))%>%
+ mutate(vax_status_combined=case_when(vax_yn_record=="confirmed_vax" ~ TRUE,
+ vax_yn_record=="confirmed_unvax" ~ FALSE,
+ mum_thinks_vaccinated=="thinks_vax" ~ TRUE,
+ mum_thinks_vaccinated=="thinks_unvax" ~ FALSE,
+ .default = NA))%>%
+ mutate(vax_status=case_when(vax_yn==1 ~ TRUE,
+ vax_yn==0 ~ FALSE))%>%
+ mutate(pertussis_status=case_when(pert_record==1 ~ TRUE,
+ pert_record==0 ~ FALSE))%>%
+ #mutate(sex=case_when(sex==2 ~"female",
+ # sex==1 ~"male"))%>%
+ mutate(male_gender=case_when(sex==2 ~FALSE,sex==1 ~TRUE))%>%
+ mutate(age_days = time_length(interval(ymd(dob_infant),ymd(attend_dat)), "day"))%>%
+ mutate(age_mnths_calc = time_length(interval(ymd(dob_infant),ymd(attend_dat)), "month"))%>%
+ mutate(month=month(attend_dat))%>%
+ mutate(age_under_3_months = if_else(age_mnths_calc<4,TRUE,FALSE))%>%
+ mutate(case=case_when(virus %in% c("rsv") ~ TRUE,
+ virus %in% c("none","other") ~ FALSE))%>%
+ rename("comorb-other"=comorb_othr)%>%
+ mutate(across(starts_with("comorb_"), ~ as.logical(.x)))%>%
+ mutate(across(starts_with("treat_hosp_"), ~ as.logical(.x)))%>%
+ mutate(imd_eng=case_when(imd_eng %in% 1:2 ~ 1,
+ imd_eng %in% 3:4 ~ 2,
+ imd_eng %in% 5:6 ~ 3,
+ imd_eng %in% 7:8 ~ 4,
+ imd_eng %in% 9:10 ~ 5,
+ .default = NA))%>%
+ mutate(imd_scot_2=if_else(imd_scot_2==99,NA,imd_scot_2))%>%
+ mutate(imd=coalesce(imd_eng,imd_scot_2))%>%
+ mutate(breastfed=if_else(breast==99,NA,as.logical(breast)))%>%
+ mutate(maternal_ethnicity=case_when(mat_eth==1 ~ "White",
+ mat_eth==2 ~ "Asian",
+ mat_eth==3 ~ "Black",
+ mat_eth==4 ~ "Mixed",
+ mat_eth==99 ~ "Other",
+ mat_eth==999 ~ NA,
+ .default = NA))%>%
+ mutate(maternal_white=case_when(mat_eth==1 ~ TRUE,
+ mat_eth %in% c(2:99) ~ FALSE,
+ mat_eth==999 ~ NA,
+ .default = NA))%>%
+ mutate(prem=comorb_1)%>%
+ filter(attend_dat<dmy("21-01-2025"))
+
+```
+## Data quality check
+### incomplete forms / data errors
+```{r incomplete}
+#| echo: false
+#| message: false
+
+ineligible_exclusions=c("gla-1203","gla-1188","add-580")
+
+#exclude GLA-1203 and GLA-1188
+#exclude add-580 (not admitted)
+cat(glue("Excluding {length(ineligible_exclusions)} who weren't eligible to recruit. ({paste(ineligible_exclusions, collapse=', ')})"))
+data<-data%>%filter(!record_id %in% ineligible_exclusions)
+
+print(glue("Number of remaining subjects: {data%>%filter(consented==TRUE)%>%nrow()}"))
+
+```
+
+\newpage
+### assessing maternal vaccine recall accuracy
+```{r mat_vax_recall_accuracy}
+#| echo: false
+#| message: false
+
+data%>%
+ tabyl(vax_yn_record,mum_thinks_vaccinated)%>%
+ adorn_totals(where=c("row", "col"))%>%gt()%>%
+ tab_header(title = md("validation of maternal recall"))
+
+library(epiR)
+vax_stats<-data%>%
+ mutate(pred = case_when(mum_thinks_vaccinated=="thinks_vax" ~ "vax",
+ mum_thinks_vaccinated=="thinks_unvax" ~ "unvax"),
+ truth = case_when(vax_yn_record=="confirmed_vax" ~ "vax",
+ vax_yn_record=="confirmed_unvax" ~ "unvax"))%>%
+ tabyl(pred,truth,show_na = F)%>%as.data.frame()%>%arrange(desc(pred))%>%select(vax,unvax)%>%
+ as.matrix()
+
+colnames(vax_stats)<-c("recall_vax","recall_unvax")
+rownames(vax_stats)<-c("conf_vax","conf_unvax")
+
+accuracy<-binom.test((vax_stats[1,1]+vax_stats[2,2]),sum(vax_stats))%>%
+ tidy()%>%
+ select(estimate,conf.low, conf.high)%>%
+ mutate(across(everything(),~sprintf("%0.1f%%", .x * 100)))
+
+cat(glue("Accuracy of maternal recall {accuracy$estimate} ({accuracy$conf.low}-{accuracy$conf.high})"))
+
+vax_stats%>%
+ epi.tests(., conf.level = 0.95)
+
+```
+
+\newpage
+### Exclusions data
+```{r exclusions}
+#| echo: false
+#| message: false
+
+exclusionTable<-data.frame(reason=NULL,count=NULL)
+# Had monoclonals:
+hadMabs<-data%>%filter(paliz==1 |
+ mabs_othr_1==1 |
+ mabs_othr_1==2)%>%
+ select(record_id)
+cat(glue("Excluding {nrow(hadMabs)} who had had monoclonal antibodies in last month"))
+if(nrow(hadMabs)>0){
+ data<-data%>%filter(!record_id %in% hadMabs$record_id)
+ exclusionTable<-rbind(exclusionTable, data.frame(reason="had mabs",count=nrow(hadMabs)))
+}
+
+# records not marked as complete:
+notComplete<-data%>%filter(follow_up_form_complete=="incomplete")%>%
+ select(record_id)
+
+cat(glue("Excluding {nrow(notComplete)} whos records were not marked as complete"))
+if(nrow(notComplete)>0){
+ usefuls<-data%>%filter(record_id %in% notComplete$record_id)%>%select(record_id,case,vax_yn)%>%drop_na(case)%>%filter(!vax_yn==99)
+ cat(glue("\t {nrow(usefuls)} of these have enough to be useful though. ({paste(usefuls$record_id, collapse=', ')})"))
+ data<-data%>%filter(!record_id %in% notComplete$record_id)
+ exclusionTable<-rbind(exclusionTable, data.frame(reason="not completed",count=nrow(notComplete)))
+}
+
+# incomplete testing data:
+incompleteTests<-data%>%filter(!virus %in% c("rsv","none","other"))%>%select(record_id)
+cat(glue("Excluding {nrow(incompleteTests)} who had had no PCR result"))
+if(nrow(incompleteTests)>0){
+ data<-data%>%filter(!record_id %in% incompleteTests$record_id)
+ exclusionTable<-rbind(exclusionTable, data.frame(reason="not tested",count=nrow(incompleteTests)))
+}
+
+# # mum did not know if vaccinated or not:
+ unknownVax<-data%>%filter(is.na(vax_status))%>%select(record_id)
+ cat(glue("Excluding {nrow(unknownVax)} with unknown vaccine status"))
+ if(nrow(unknownVax)>0){
+ data<-data%>%filter(!record_id %in% unknownVax$record_id)
+ exclusionTable<-rbind(exclusionTable, data.frame(reason="unknown vax status",count=nrow(unknownVax)))
+ }
+
+print(glue("Number of remaining subjects: {data%>%filter(consented==TRUE)%>%nrow()}"))
+
+exclusionTable%>%
+ gt()%>%tab_header(title = md("Exclusions"))
+
+data%>%select(record_id,age_mnths_calc,dob_infant,attend_dat)%>%arrange(age_mnths_calc)%>%top_n(5,age_mnths_calc)%>%
+ gt()%>%tab_header("oldest kids and their DoBs")
+
+cat("min & max recruitment dates....")
+summary(data$attend_dat)
+
+```
+
+```{r looking_at_gestations}
+#| echo: false
+#| message: false
+#| include: true
+
+#if use recall data as priority
+under2weeks_recall<-
+data%>%
+ #calculate differences with gestations
+ mutate(vax_gest=if_else(vax_gest==99,NA,vax_gest))%>%
+ mutate(vax_gest_record=if_else(vax_gest_record==99,NA,vax_gest_record))%>%
+ mutate(age_gest=if_else(age_gest==99,NA,age_gest))%>%
+ mutate(vax_gest_combined=coalesce(vax_gest,vax_gest_record))%>%
+ mutate(gestVax_to_gestBorn=age_gest-vax_gest_combined)%>%
+
+ #calculate differences with dates
+ mutate(vax_dat_combined=coalesce(vax_dat,vax_dat_record))%>%
+ mutate(dob_to_vaxdate=time_length(interval(ymd(vax_dat_combined),ymd(dob_infant)), "weeks"))%>%
+
+ filter(vax_yn==1)%>%
+ mutate(dob_to_vaxdate=if_else(dob_to_vaxdate<0,NA,dob_to_vaxdate))%>%
+ mutate(vax_diff=coalesce(dob_to_vaxdate,gestVax_to_gestBorn))%>%
+ mutate(under2weeks=if_else(vax_diff<2,TRUE,FALSE))%>%
+ select(record_id,under2weeks,vax_diff,dob_to_vaxdate,gestVax_to_gestBorn,vax_gest,vax_gest_record,age_gest,vax_dat,vax_dat_record,dob_infant)
+
+
+#if use vax data as priority
+under2weeks_records<-
+ data%>%
+ #calculate differences with gestations
+ mutate(vax_gest=if_else(vax_gest==99,NA,vax_gest))%>%
+ mutate(vax_gest_record=if_else(vax_gest_record==99,NA,vax_gest_record))%>%
+ mutate(age_gest=if_else(age_gest==99,NA,age_gest))%>%
+ mutate(vax_gest_combined=coalesce(vax_gest_record,vax_gest))%>%
+ mutate(gestVax_to_gestBorn=age_gest-vax_gest_combined)%>%
+
+ #calculate differences with dates
+ mutate(vax_dat_combined=coalesce(vax_dat_record,vax_dat))%>%
+ mutate(dob_to_vaxdate=time_length(interval(ymd(vax_dat_combined),ymd(dob_infant)), "weeks"))%>%
+
+ filter(vax_yn==1)%>%
+ mutate(dob_to_vaxdate=if_else(dob_to_vaxdate<0,NA,dob_to_vaxdate))%>%
+ mutate(vax_diff=coalesce(dob_to_vaxdate,gestVax_to_gestBorn))%>%
+ mutate(under2weeks=if_else(vax_diff<2,TRUE,FALSE))%>%
+ select(record_id,under2weeks,vax_diff,dob_to_vaxdate,gestVax_to_gestBorn,vax_gest,vax_gest_record,age_gest,vax_dat,vax_dat_record,dob_infant)
+
+
+under2weeks_recall%>%tabyl(under2weeks)%>%gt()%>%tab_header("prioritising recall to identify <2 weeks since vaccine")
+
+under2weeks_records%>%tabyl(under2weeks)%>%gt()%>%tab_header("prioritising records to identify <2 weeks since vaccine")
+
+data<-data%>%
+ left_join(under2weeks_records%>%select(record_id,under2weeks))
+
+```
+
+```{r plotting function}
+#| echo: false
+#| message: false
+#| include: false
+
+calculate_VE<-function(vaxData){
+
+ t1<-vaxData%>%
+ mutate(case=if_else(case,"case","control"))%>%
+ mutate(vax_status=if_else(vax_status,"vax","unvax"))%>%
+ tabyl(vax_status,case)%>%
+ adorn_totals(where=c("row", "col"))%>%
+ rename(" "=vax_status)%>%
+ gt()%>%tab_header(title = md("2x2 data feeding into regression"))
+
+tableOutput<-data.frame(formula=NULL,VE=NULL,additional=NULL)
+
+VE_calc<-clogit(formula=case ~ vax_status, data=vaxData)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(estimate=1-estimate,
+ conf_high=1-conf.low,
+ conf_low=1-conf.high)%>%
+ select(estimate,conf_low,conf_high)
+
+tableOutput<-rbind(tableOutput,data.frame(
+ formula="Unadjusted",
+ VE=glue("{100 * round(VE_calc$estimate,3)}% ({100 * round(VE_calc$conf_low,3)}% - {100 * round(VE_calc$conf_high,3)}%)"),
+ additional="")
+)
+
+VE_calc<-clogit(formula=case ~ vax_status + strata(month), data=vaxData)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(estimate=1-estimate,
+ conf_high=1-conf.low,
+ conf_low=1-conf.high)%>%
+ select(estimate,conf_low,conf_high)
+
+tableOutput<-rbind(tableOutput,data.frame(
+ formula="month stratified",
+ VE=glue("{100 * round(VE_calc$estimate,3)}% ({100 * round(VE_calc$conf_low,3)}% - {100 * round(VE_calc$conf_high,3)}%)"),
+ additional="")
+)
+
+VE_calc<-clogit(formula=case ~ vax_status + strata(siteID) + strata(month), data=vaxData)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(estimate=1-estimate,
+ conf_high=1-conf.low,
+ conf_low=1-conf.high)%>%
+ select(estimate,conf_low,conf_high)
+
+tableOutput<-rbind(tableOutput,data.frame(
+ formula="site/month stratified",
+ VE=glue("{100 * round(VE_calc$estimate,3)}% ({100 * round(VE_calc$conf_low,3)}% - {100 * round(VE_calc$conf_high,3)}%)"),
+ additional="")
+)
+
+VE_calc<-clogit(formula=case ~ vax_status + strata(siteID) + strata(month) + breast, data=vaxData)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(estimate=1-estimate,
+ conf_high=1-conf.low,
+ conf_low=1-conf.high)%>%
+ select(estimate,conf_low,conf_high)
+
+tableOutput<-rbind(tableOutput,data.frame(
+ formula="site/month stratified + BF",
+ VE=glue("{100 * round(VE_calc[1,]$estimate,3)}% ({100 * round(VE_calc[1,]$conf_low,3)}% - {100 * round(VE_calc[1,]$conf_high,3)}%)"),
+ additional=glue("Breastfeeding OR {1-round(VE_calc[2,]$estimate,2)} ({1-round(VE_calc[2,]$conf_high,2)} - {1-round(VE_calc[2,]$conf_low,2)})"))
+)
+
+VE_calc<-clogit(formula=case ~ vax_status + strata(siteID) + strata(month) + age_under_3_months, data=vaxData)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(estimate=1-estimate,
+ conf_high=1-conf.low,
+ conf_low=1-conf.high)%>%
+ select(estimate,conf_low,conf_high)
+
+tableOutput<-rbind(tableOutput,data.frame(
+ formula="site/month stratified + age",
+ VE=glue("{100 * round(VE_calc[1,]$estimate,3)}% ({100 * round(VE_calc[1,]$conf_low,3)}% - {100 * round(VE_calc[1,]$conf_high,3)}%)"),
+ additional=glue("<3m OR {1-round(VE_calc[2,]$estimate,2)} ({1-round(VE_calc[2,]$conf_high,2)} - {1-round(VE_calc[2,]$conf_low,2)})"))
+)
+
+VE_calc<-clogit(formula=case ~ vax_status + strata(siteID) + strata(month) + age_under_3_months, data=vaxData)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(estimate=1-estimate,
+ conf_high=1-conf.low,
+ conf_low=1-conf.high)%>%
+ select(estimate,conf_low,conf_high)
+
+tableOutput<-rbind(tableOutput,data.frame(
+ formula="site/month stratified + age",
+ VE=glue("{100 * round(VE_calc[1,]$estimate,3)}% ({100 * round(VE_calc[1,]$conf_low,3)}% - {100 * round(VE_calc[1,]$conf_high,3)}%)"),
+ additional=glue("<3m OR {1-round(VE_calc[2,]$estimate,2)} ({1-round(VE_calc[2,]$conf_high,2)} - {1-round(VE_calc[2,]$conf_low,2)})"))
+)
+
+
+
+VE_calc<-clogit(formula=case ~ pertussis_status, data=vaxData)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(estimate=1-estimate,
+ conf_high=1-conf.low,
+ conf_low=1-conf.high)%>%
+ select(estimate,conf_low,conf_high)
+
+tableOutput<-rbind(tableOutput,data.frame(
+ formula="pertussis VE vs RSV",
+ VE=glue("{100 * round(VE_calc$estimate,3)}% ({100 * round(VE_calc$conf_low,3)}% - {100 * round(VE_calc$conf_high,3)}%)"),
+ additional="")
+)
+
+VE_calc<-clogit(formula=case ~ pertussis_status + strata(siteID), data=vaxData)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(estimate=1-estimate,
+ conf_high=1-conf.low,
+ conf_low=1-conf.high)%>%
+ select(estimate,conf_low,conf_high)
+
+tableOutput<-rbind(tableOutput,data.frame(
+ formula="pertussis VE vs RSV (site stratified)",
+ VE=glue("{100 * round(VE_calc$estimate,3)}% ({100 * round(VE_calc$conf_low,3)}% - {100 * round(VE_calc$conf_high,3)}%)"),
+ additional="")
+)
+
+
+t2<-tableOutput%>%gt()
+
+ return(list(table1=t1,table2=t2))
+}
+
+calculate_adjusted_VE<-function(vaxData){
+
+ VE_calc<-clogit(formula=case ~ vax_status + strata(siteID) + strata(month) + age_under_3_months + male_gender + prem, data=vaxData)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(VEestimate=1-estimate, VEconf_high=1-conf.low, VEconf_low=1-conf.high)%>%
+ select(estimate,conf.low,conf.high,VEestimate,VEconf_high,VEconf_low)
+
+tableOutput<-data.frame(variable=NULL,estimate=NULL)
+tableOutput<-rbind(tableOutput,data.frame(
+ variable="vaccine effectiveness",
+ estimate=glue("{100 * round(VE_calc[1,]$VEestimate,3)}% ({100 * round(VE_calc[1,]$VEconf_low,3)}% - {100 * round(VE_calc[1,]$VEconf_high,3)}%)")))
+
+tableOutput<-rbind(tableOutput,data.frame(
+ variable="OR age <3m",
+ estimate=glue("{round(VE_calc[2,]$estimate,2)} ({round(VE_calc[2,]$conf.low,2)} - {round(VE_calc[2,]$conf.high,2)})")))
+
+tableOutput<-rbind(tableOutput,data.frame(
+ variable="OR male",
+ estimate=glue("{round(VE_calc[3,]$estimate,2)} ({round(VE_calc[3,]$conf.low,2)} - {round(VE_calc[3,]$conf.high,2)})")))
+
+tableOutput<-rbind(tableOutput,data.frame(
+ variable="OR <37w gestation",
+ estimate=glue("{round(VE_calc[4,]$estimate,2)} ({round(VE_calc[4,]$conf.low,2)} - {round(VE_calc[4,]$conf.high,2)})")))
+
+tableOutput%>%gt()%>%tab_header(title="fully adjusted model",subtitle = "case ~ vax_status + strata(siteID) + strata(month) + age_under_3_months + male_gender + prem")%>%return()
+
+}
+
+make_summary_tables<-function(vaxData){
+t1<-vaxData%>%
+ drop_na(case)%>%
+ mutate(case=if_else(case,"case","control"))%>%
+ mutate(male_gender=if_else(male_gender,"Male","Female"))%>%
+ mutate(prem=if_else(prem,"Born preterm (<37w)","Born at term"))%>%
+ mutate(maternal_white=if_else(maternal_white,"White","non-white"))%>%
+ mutate(breastfed=if_else(breastfed,"Yes","No"))%>%
+ mutate(any_oxygen=treat_hosp_3|treat_hosp_4|treat_hosp_5|treat_hosp_6)%>%
+ tbl_summary(include = c(age_mnths_calc, male_gender, prem,comorb_2:comorb_6 ,maternal_white,breastfed,imd,
+ any_oxygen,treat_hosp_4,treat_hosp_5,treat_hosp_6,admit_locn_3,admit_locn_4),
+ type=list(age_mnths_calc~"continuous",
+ male_gender~"categorical",
+ prem~"categorical",
+ breastfed~"categorical",
+ imd~"continuous",
+ maternal_white="categorical",
+ any_oxygen="dichotomous",
+ treat_hosp_4="dichotomous",
+ treat_hosp_5="dichotomous",
+ treat_hosp_6="dichotomous",
+ admit_locn_3="dichotomous",
+ admit_locn_4="dichotomous"
+ ),
+ label=(list(age_mnths_calc~"Age in months",
+ male_gender~"Male gender",
+ prem~"Gestation at birth",
+ breastfed ="Breastfed",
+ imd="Socioeconomic quintile",
+ maternal_white="Maternal ethnicity",
+ comorb_2="Chronic lung disease",
+ comorb_3="Congenital cardiac disease",
+ comorb_4="Neuromuscular disease",
+ comorb_5="Other",
+ comorb_6="None",
+ any_oxygen="Supplemental oxygen use",
+ treat_hosp_4="High flow nasal cannulae",
+ treat_hosp_5="Continuous positive airway pressure",
+ treat_hosp_6="Invasive ventilation",
+ admit_locn_3="HDU admission",
+ admit_locn_4="PICU admission"
+ )),
+ by=case)%>%
+ add_p()%>%
+ modify_caption("**Table 1. Characteristics of the included case and control patients **")
+
+t2<-vaxData%>%
+ filter(case==T)%>%
+ mutate(vax_status=if_else(vax_status,"vaccinated","un-vaccinated"))%>%
+ tbl_summary(include = c(LOA, treat_hosp_1:treat_hosp_6,admit_locn_3,admit_locn_4),
+ label=(list(
+ LOA="Length of admission",
+ treat_hosp_1="NG feeds",
+ treat_hosp_2="IV fluids",
+ treat_hosp_3="Low flow O2",
+ treat_hosp_4="HFNCO2",
+ treat_hosp_5="CPAP",
+ treat_hosp_6="Ventilation",
+ admit_locn_3="admitted to HDU",
+ admit_locn_4="admitted to PICU"
+ )),
+ by=vax_status)%>%
+ add_p()%>%
+ modify_caption("**Table 2. Comparison of clinical outcomes for RSV-positive infants with mothers vaccinated more than 14 days prior to delivery compared to infants with unvaccinated mothers **")
+
+
+
+
+ return(list(table1=t1,table2=t2))
+}
+
+```
+\newpage
+# Outputs using maternal recall (whole cohort)
+```{r VE_estimates_maternal}
+#| echo: false
+#| message: false
+
+output<-data%>%
+ select(case,vax_status,siteID,month,breast,age_under_3_months,age_mnths_calc,pertussis_status)%>%
+ drop_na(case, vax_status)%>%
+ calculate_VE()
+
+output$table1
+
+output$table2%>%tab_header(title = md("VE maternal recall"))
+
+data%>%
+ select(case,vax_status,siteID,month,breast,age_under_3_months,age_mnths_calc,pertussis_status, male_gender,prem)%>%
+ drop_na(case, vax_status)%>%
+ calculate_adjusted_VE()
+
+tables<-data%>%make_summary_tables()
+tables$table1
+tables$table2
+
+```
+\newpage
+# VE estimates excluding those who we think had it <2 weeks before birth
+```{r VE_estimates_under2weeks1}
+#| echo: false
+#| message: false
+
+remove_under_2weeks<-data%>%filter(under2weeks %in% c(NA,FALSE))
+
+remove_under_2weeks%>% mutate(vax_status=if_else(vax_status,"vax","unvax"))%>%
+ tabyl(vax_status,under2weeks)%>%rename(" "=vax_status)%>%
+ gt()%>%tab_header(title = md("subjects known to be vaccinated <2 weeks before birth"))
+
+remove_under_2weeks%>%
+ mutate(case=if_else(case,"case","control"))%>%
+ mutate(vax_status=if_else(vax_status,"vax","unvax"))%>%
+ tabyl(vax_status,case)%>%
+ adorn_totals(where=c("row", "col"))%>%
+ rename(" "=vax_status)%>%
+ gt()%>%tab_header(title = md("2x2 data feeding into regression"))
+
+remove_under_2weeks%>%
+ calculate_adjusted_VE()%>%tab_header("Fully adjusted model removing <2 weeks")
+
+tables<-remove_under_2weeks%>%make_summary_tables()
+tables$table1
+tables$table2
+
+```
+
+\newpage
+# VE estimates including only those who we know had it >2 weeks before birth
+```{r VE_estimates_under2weeks2}
+#| echo: false
+#| message: false
+
+only_over_2weeks<-data%>%filter((vax_status & !under2weeks) | !vax_status)
+
+only_over_2weeks%>% mutate(vax_status=if_else(vax_status,"vax","unvax"))%>%
+ tabyl(vax_status,under2weeks)%>%rename(" "=vax_status)%>%
+ gt()%>%tab_header(title = md("subjects known to be vaccinated <2 weeks before birth"))
+
+only_over_2weeks%>%
+ mutate(case=if_else(case,"case","control"))%>%
+ mutate(vax_status=if_else(vax_status,"vax","unvax"))%>%
+ tabyl(vax_status,case)%>%
+ adorn_totals(where=c("row", "col"))%>%
+ rename(" "=vax_status)%>%
+ gt()%>%tab_header(title = md("2x2 data feeding into regression"))
+
+output<-only_over_2weeks%>%
+ calculate_VE()
+output$table1
+output$table2%>%tab_header(title = md("VE only under 2 weeks"))
+
+only_over_2weeks%>%
+ calculate_adjusted_VE()%>%tab_header("Fully adjusted model for known >=2 weeks")
+
+tables<-only_over_2weeks%>%make_summary_tables()
+tables$table1
+tables$table2
+
+```
+
+\newpage
+# Fully Adjusted pertussis
+```{r VE_adj_pertussis}
+#| echo: false
+#| message: false
+
+#have to fudge month & swap vax for pertussus...
+pertussis<-data%>%
+ mutate(vax_status=pertussis_status)%>%
+ mutate(siteID=1)
+
+pertussis%>%
+ mutate(case=if_else(case,"case","control"))%>%
+ mutate(vax_status=if_else(vax_status,"vax","unvax"))%>%
+ tabyl(vax_status,case)%>%
+ adorn_totals(where=c("row", "col"))%>%
+ rename(" "=vax_status)%>%
+ gt()%>%tab_header(title = md("2x2 data feeding into regression"))
+
+pertussis%>%
+ calculate_adjusted_VE()%>%tab_header("Fully adjusted model for pertussis")
+
+```
+
+\newpage
+# VE estimates using the records data we have
+```{r VE_estimates_records}
+#| echo: false
+#| message: false
+
+records_only<-data%>%mutate(vax_status = case_when(vax_yn_record=="confirmed_vax" ~ TRUE,
+ vax_yn_record=="confirmed_unvax" ~ FALSE))
+
+records_only%>%
+ mutate(case=if_else(case,"case","control"))%>%
+ mutate(vax_status=if_else(vax_status,"vax","unvax"))%>%
+ tabyl(vax_status,case)%>%
+ adorn_totals(where=c("row", "col"))%>%
+ rename(" "=vax_status)%>%
+ gt()%>%tab_header(title = md("2x2 data feeding into regression"))
+
+records_only%>%
+ calculate_adjusted_VE()%>%tab_header("Just using the maternal records we have")
+
+tables<-records_only%>%make_summary_tables()
+tables$table1
+tables$table2
+
+records_only%>%filter((vax_status & !under2weeks) | !vax_status)%>%
+ mutate(case=if_else(case,"case","control"))%>%
+ mutate(vax_status=if_else(vax_status,"vax","unvax"))%>%
+ tabyl(vax_status,case)%>%
+ adorn_totals(where=c("row", "col"))%>%
+ rename(" "=vax_status)%>%
+ gt()%>%tab_header(title = md("2x2 data feeding into regression if we only included <2weeks"))
+
+records_only%>%filter((vax_status & !under2weeks) | !vax_status)%>%
+ calculate_adjusted_VE()%>%tab_header("Just using the maternal records we have confirmed over 2 weeks!")
+
+```
+
+```{r VE_estimates_combined}
+#| echo: false
+#| message: false
+#| include: false
+
+# VE estimates using combined maternal + any records data we have
+
+output<-data%>%
+ mutate(vax_status = vax_status_combined)%>%
+ select(case,vax_status,siteID,month,breast,age_under_3_months,age_mnths_calc,pertussis_status)%>%
+ drop_na(case, vax_status)%>%
+ calculate_VE()
+
+output$table1
+
+output$table2%>%tab_header(title = md("VE combined records data"))
+
+```
+
+\newpage
+## Screening method comparision of subgroups vs all controls (unadjusted)
+(nb I don't rank the other treatments at the moment.... so no cpap can include ventilation.)
+```{r screening_method}
+#| echo: false
+#| message: false
+
+subgroupVE<-function(subgroup,outcome){
+VE_calc<-clogit(formula=case ~ vax_status, data=subgroup)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(estimate=1-estimate,
+ conf_high=1-conf.low,
+ conf_low=1-conf.high)%>%
+ mutate(measure=glue("VE against {outcome}"))%>%
+ select(measure,estimate,conf_low,conf_high)
+cat(glue("{VE_calc$measure} {100 * round(VE_calc$estimate,3)}% ({100 * round(VE_calc$conf_low,3)}% - {100 * round(VE_calc$conf_high,3)}%)"))
+
+subgroup%>%
+ mutate(case=if_else(case,"case","control"))%>%
+ mutate(vax_status=if_else(vax_status,"vax","unvax"))%>%
+ tabyl(case,vax_status)%>%adorn_totals(where=c("row", "col"))%>%rename(" "=case)%>%
+ mutate(percVac=vax/Total)%>%mutate(percUnvac=unvax/Total)%>%
+ gt()%>%fmt_percent(starts_with("perc"),decimals = 1)%>%
+ tab_header(title = glue("{outcome} vs all controls"))%>%
+ return()
+}
+
+outcomes<-data%>%
+ select(case,vax_status,siteID,month,breast,age_under_3_months,age_mnths_calc,pertussis_status,treat_hosp_4,treat_hosp_5,treat_hosp_6,
+ admit_locn_3,admit_locn_4,under2weeks)%>%
+ rename(treat_hfno=treat_hosp_4,treat_cpap=treat_hosp_5,treat_vent=treat_hosp_6,
+ admit_hdu=admit_locn_3,admit_picu=admit_locn_4)%>%
+ drop_na(case, vax_status)
+
+outcomes%>%filter(case & treat_hfno | !case)%>%subgroupVE("HFNCO")
+
+outcomes%>%filter(case & treat_cpap | !case)%>%subgroupVE("CPAP")
+
+outcomes%>%filter(case & treat_vent | !case)%>%subgroupVE("vent")
+
+outcomes%>%filter(case & admit_hdu | !case)%>%subgroupVE("hdu")
+
+outcomes%>%filter(case & admit_picu | !case)%>%subgroupVE("picu")
+
+```
+\newpage
+## Screening method comparision of subgroups vs all controls and only under 2 weeks (unadjusted)
+```{r screening_method_over_2_weeks}
+#| echo: false
+#| message: false
+
+outcomes<-outcomes%>%
+ filter((vax_status & !under2weeks) | !vax_status)
+
+outcomes%>%filter(case & treat_hfno | !case)%>%subgroupVE("HFNCO")
+
+outcomes%>%filter(case & treat_cpap | !case)%>%subgroupVE("CPAP")
+
+outcomes%>%filter(case & treat_vent | !case)%>%subgroupVE("vent")
+
+outcomes%>%filter(case & admit_hdu | !case)%>%subgroupVE("hdu")
+
+outcomes%>%filter(case & admit_picu | !case)%>%subgroupVE("picu")
+
+```
+
+\newpage
+## Screening method comparision of subgroups vs all controls. (stratified by month and site)
+can't do more adjustments or it breaks
+```{r screening_method_stratified}
+#| echo: false
+#| message: false
+
+subgroupVE<-function(subgroup,outcome){
+VE_calc<-clogit(formula=case ~ vax_status + strata(siteID) + strata(month), data=subgroup)%>%
+ tidy(conf.int=TRUE, exponentiate=TRUE)%>%
+ mutate(estimate=1-estimate,
+ conf_high=1-conf.low,
+ conf_low=1-conf.high)%>%
+ mutate(measure=glue("VE against {outcome}"))%>%
+ select(measure,estimate,conf_low,conf_high)
+cat(glue("{VE_calc$measure} {100 * round(VE_calc$estimate,3)}% ({100 * round(VE_calc$conf_low,3)}% - {100 * round(VE_calc$conf_high,3)}%)"))
+
+subgroup%>%
+ mutate(case=if_else(case,"case","control"))%>%
+ mutate(vax_status=if_else(vax_status,"vax","unvax"))%>%
+ tabyl(case,vax_status)%>%adorn_totals(where=c("row", "col"))%>%
+ mutate(percVac=vax/Total)%>%mutate(percUnvac=unvax/Total)%>%
+ rename(" "=case)%>%gt()%>%fmt_percent(starts_with("perc"),decimals = 1)%>%
+ tab_header(title = glue("{outcome} vs all controls"))%>%
+ return()
+}
+
+
+
+outcomes<-data%>%
+ select(case,vax_status,siteID,month,breast,age_under_3_months,age_mnths_calc,pertussis_status,treat_hosp_4,treat_hosp_5,treat_hosp_6,
+ admit_locn_3,admit_locn_4,under2weeks,male_gender,prem)%>%
+ rename(treat_hfno=treat_hosp_4,treat_cpap=treat_hosp_5,treat_vent=treat_hosp_6,
+ admit_hdu=admit_locn_3,admit_picu=admit_locn_4)%>%
+ drop_na(case, vax_status)
+
+outcomes%>%filter(case & treat_hfno | !case)%>%subgroupVE("HFNCO")
+
+outcomes%>%filter(case & treat_cpap | !case)%>%subgroupVE("CPAP")
+
+outcomes%>%filter(case & treat_vent | !case)%>%subgroupVE("vent")
+
+outcomes%>%filter(case & admit_hdu | !case)%>%subgroupVE("hdu")
+
+outcomes%>%filter(case & admit_picu | !case)%>%subgroupVE("picu")
+
+```
--
GitLab