# Analysis code for MicronutTB by Schwalb et al. 2023, 
# Distributed under CC BY 4.0 license: https://creativecommons.org/licenses/by/4.0/

# Packages ---
library(rio)
library(here)
library(haven)
library(gmodels)
library(tidyverse)
library(ggplot2)
library(survival)
library(survminer)
library(lmtest)
library(Greg)
library(RColorBrewer)

# Data ---
micronutDB <- read_dta("data/cr_3.final_dataset_two.dta") # streamlined DB

# Fixes
micronutDB_extra <- read_dta("data/cr_3.final_dataset_one.dta") # extra variables
micronutDB$fu <- as.double(micronutDB$date_exit - micronutDB$date_recr)
micronutDB$dob <- micronutDB_extra$dob # add DOB
micronutDB$tbdeath <- micronutDB_extra$newtb_death # add TB death
micronutDB$time_arv <- as.double(micronutDB$date_exit - micronutDB$date_arv) # create time to ART
micronutDB$time_arv <- replace(micronutDB$time_arv, which(micronutDB$time_arv < 0), NA) # clean time to ART
micronutDB$wk8 <- as.numeric(micronutDB$fu >= 56)

# Descriptive analysis ---
# Follow-up and deaths
sum(micronutDB$fu) # total FU
sum(micronutDB$fu)/365.25 # total FU (person-years)
sum(micronutDB$fu[micronutDB$arm == 0]) # total FU (arm = LNS)
sum(micronutDB$fu[micronutDB$arm == 0])/365.25 # total FU (person-years) (arm = LNS)
sum(micronutDB$fu[micronutDB$arm == 1]) # total FU (arm = LNS-VM)
sum(micronutDB$fu[micronutDB$arm == 1])/365.25 # total FU (person-years) (arm = LNS-VM)


summary(micronutDB$fu) # median FU
table(micronutDB$xreason3) # lost to follow-up

table(!is.na(micronutDB$date_arv)) # counts all who started ART
table(micronutDB$admitted) # hospital admission
table(micronutDB$died) # total deaths

sum(micronutDB$died)/(sum(micronutDB$fu)/365)*10 # mortality rate per 10 person-year
poisson.test(sum(micronutDB$died), (sum(micronutDB$fu)/365))$conf.int[1]*10
poisson.test(sum(micronutDB$died), (sum(micronutDB$fu)/365))$conf.int[2]*10

# Demographic characteristics
table(micronutDB$site) # site frequencies

table(micronutDB$arm) # arm frequencies

micronutDB$age <- (as.numeric(micronutDB$date_recr - micronutDB$dob)/365) # create age variable
hist(micronutDB$age) # check histogram
summary(micronutDB$age) # median age
table(micronutDB$agegp2) # age group frequencies
CrossTable(micronutDB$agegp2, micronutDB$arm, chisq = TRUE) # crosstab age group x arm

table(micronutDB$sex) # sex frequencies
CrossTable(micronutDB$sex, micronutDB$arm, chisq = TRUE) # crosstab sex x arm

table(micronutDB$ses3) # SES tertile group frequencies
CrossTable(micronutDB$ses3, micronutDB$arm, chisq = TRUE) # crosstab SES tertile x arm

table(micronutDB$bmigp) # BMI group frequencies
CrossTable(micronutDB$bmigp, micronutDB$arm, chisq = TRUE) # crosstab BMI group x arm

table(micronutDB$cd4100) # CD4 binary group frequencies
CrossTable(micronutDB$cd4100, micronutDB$arm, chisq = TRUE) # crosstab CD4 group x arm

table(micronutDB$hbgp2) # Hb binary group frequencies
CrossTable(micronutDB$hbgp2, micronutDB$arm, chisq = TRUE) # crosstab Hb group x arm

table(micronutDB$crpgp) # CRP group frequencies
CrossTable(micronutDB$crpgp, micronutDB$arm, chisq = TRUE) # crosstab CD4 group x arm

# TB incidence and mortality ---
sum(micronutDB$tbdeath, na.rm = TRUE) # TB mortality rate

sum(micronutDB$newtb) # TB incident cases
sum(micronutDB$newtb)/(sum(micronutDB$fu)/365) # TB incidence rate

aggregate(micronutDB$fu, by = list(micronutDB$arm), FUN = sum) # FU per arm: 0(LNS) and 1 (LNS-VM)
aggregate(micronutDB$newtb, by = list(micronutDB$arm), FUN = sum) # TB incidence per arm

# LNS = 49713/365 = 136.2 // 137 TB 
LNSev <- 137
LNSt <- 136.2
poisson.test(LNSev, LNSt)$conf.int[1]*10
poisson.test(LNSev, LNSt)$conf.int[2]*10

# LNS-VM = 50549/365 = 138.5 // 126 TB
LNSVMev <- 126
LNSVMt <- 138.5
poisson.test(LNSVMev, LNSVMt)$conf.int[1]*10
poisson.test(LNSVMev, LNSVMt)$conf.int[2]*10

aggregate(micronutDB$fu, by = list(micronutDB$bmigp), FUN = sum) # FU per BMI: 0 (17-18.5), 1 (16-16.9), 2 (<16)
aggregate(micronutDB$newtb, by = list(micronutDB$bmigp), FUN = sum) # TB incidence per BMI group
# 17-18.5 = 47246/365 = 129.4 // 89
# 16-16.9 = 25300/365 = 69.3 // 64
# <16 = 27716/365 = 75.9 // 110

# Survival analysis ---
surv <- Surv(time = micronutDB$fu, event = micronutDB$newtb)

coxraw <- coxph(surv ~ factor(micronutDB$arm)) # Cox regression per arm (raw)
summary(coxraw)

coxadj <- coxph(surv ~ factor(micronutDB$arm) + factor(micronutDB$agegp2) + factor(micronutDB$sex) + 
                  factor(micronutDB$bmigp) + factor(micronutDB$cd4100) + 
                  factor(micronutDB$hbgp2) + factor(micronutDB$crpgp)) # Cox regression per arm (adjusted)
summary(coxadj)

# BMI association with TB incidence ---
coxbmiraw <- coxph(surv ~ factor(micronutDB$bmigp)) # Cox regression per BMI (raw)
summary(coxbmiraw)

coxbmiadj <- coxph(surv ~ factor(micronutDB$bmigp) + factor(micronutDB$agegp2) + factor(micronutDB$sex) + 
                     factor(micronutDB$cd4100) + factor(micronutDB$hbgp2) + factor(micronutDB$crpgp)) # Cox regression per BMI (adjusted)
summary(coxbmiadj)

# LRT - BMI association
bmiestA <- coxph(surv ~ factor(micronutDB$bmigp) + factor(micronutDB$agegp2) + factor(micronutDB$sex) + 
                   factor(micronutDB$cd4100) + factor(micronutDB$hbgp2) + factor(micronutDB$crpgp)) 
bmiestB <- coxph(surv ~ factor(micronutDB$agegp2) + factor(micronutDB$sex) + 
                   factor(micronutDB$cd4100) + factor(micronutDB$hbgp2) + factor(micronutDB$crpgp)) 
lrtest(bmiestA, bmiestB)

# BMI - Interaction with sex
bmiestC <- coxph(surv ~ factor(micronutDB$bmigp)*factor(micronutDB$sex) + factor(micronutDB$agegp2) +
                   factor(micronutDB$cd4100) + factor(micronutDB$hbgp2) + factor(micronutDB$crpgp))
lrtest(bmiestA, bmiestC)

# BMI - Interaction with CD4 count
bmiestD <- coxph(surv ~ factor(micronutDB$bmigp)*factor(micronutDB$cd4100) + factor(micronutDB$agegp2) +
                   factor(micronutDB$sex) + factor(micronutDB$hbgp2) + factor(micronutDB$crpgp))
lrtest(bmiestA, bmiestD)

# BMI (linearity) association with TB incidence ---
coxbmiadjlin <- coxph(surv ~ micronutDB$bmigp + factor(micronutDB$agegp2) + factor(micronutDB$sex) + 
                     factor(micronutDB$cd4100) + factor(micronutDB$hbgp2) + factor(micronutDB$crpgp)) # Cox regression per BMI (adjusted)
summary(coxbmiadjlin)

# Survival curves ---
display.brewer.all(colorblindFriendly = T)
RColorBrewer::brewer.pal(8, 'Dark2')
RColorBrewer::brewer.pal(8, 'Set2')

# Fig 1A - Survival plot per arm (650w x 450h)
fig1a <- survfit(Surv(micronutDB$fu, micronutDB$newtb) ~ micronutDB$arm, data = micronutDB)

tiff(here("Fig2.tiff"),  width = 8, height = 6, units = 'in', res = 200) 
ggsurvplot(fig1a, ggtheme = theme_classic(), fun = "event", censor = FALSE, pval = FALSE,
           linetype = "strata", palette = c("#66C2A5", "#FC8D62"), conf.int = TRUE,
           xlab = "Time in study (months)", ylab = "Cumulative incident TB",
           legend = "bottom", legend.title = "Arm",  legend.labs = c("LNS", "LNS-VM"),
           risk.table = "abs_pct", fontsize = 2.25, ylim = c(0, 0.4), surv.scale = c("percent"),
           xscale = 30, break.x.by = 30, xlim = c(0,120), axes.offset = TRUE)
dev.off()


# Fig 1B - Survival plot per arm (650w x 450h)
fig1b <- survfit(Surv(micronutDB$fu, micronutDB$newtb) ~ micronutDB$bmigp, data = micronutDB)

tiff(here("Fig3.tiff"),  width = 8, height = 6, units = 'in', res = 200) 
ggsurvplot(fig1b, ggtheme = theme_classic(), fun = "event", censor = FALSE, pval = FALSE,
           linetype = "strata", palette = c("#66C2A5", "#FC8D62", "#8DA0CB"), conf.int = TRUE,
           xlab = "Time in study (months)", ylab = "Cumulative incident TB",
           legend = "bottom", legend.title = "BMI",  legend.labs = c("17-18.5", "16-16.9", "<16"),
           risk.table = "abs_pct", fontsize = 2.25, ylim = c(0, 0.4), surv.scale = c("percent"),
           xscale = 30, break.x.by = 30, xlim = c(0,120), axes.offset = TRUE)
dev.off()

# ART as time-dependent variable
micronutDB_timedep <- tmerge(data1 = micronutDB %>% dplyr::select(id, fu, newtb),
                             data2 = micronutDB %>% dplyr::select(id, fu, newtb, time_arv, arv),
                             id = id, incidTB = event(fu, newtb), art = tdc(arv))

cox_timedep <- coxph(Surv(time = tstart, time2 = tstop, event = micronutDB_timedep$newtb) ~ art, data = micronutDB_timedep)
summary(cox_timedep)
ggforest(cox_timedep, data = micronutDB_timedep)

# Sensitivity analysis ---
micronutDBsa <- micronutDB %>%
  filter(wk8 == 1)

sum(micronutDBsa$fu[micronutDBsa$arm == 0]) # total FU (arm = LNS)
sum(micronutDBsa$fu[micronutDBsa$arm == 0])/365.25 # total FU (person-years) (arm = LNS)
sum(micronutDBsa$fu[micronutDBsa$arm == 1]) # total FU (arm = LNS-VM)
sum(micronutDBsa$fu[micronutDBsa$arm == 1])/365.25 # total FU (person-years) (arm = LNS-VM)

micronutDBsa2 <- micronutDB %>%
  filter(wk8 == 0)

# Demographic characteristics
table(micronutDBsa$arm) # arm frequencies

micronutDBsa$age <- (as.numeric(micronutDBsa$date_recr - micronutDBsa$dob)/365) # create age variable
hist(micronutDBsa$age) # check histogram
summary(micronutDBsa$age) # median age
table(micronutDBsa$agegp2) # age group frequencies
CrossTable(micronutDBsa$agegp2, micronutDBsa$arm, chisq = TRUE) # crosstab age group x arm

table(micronutDBsa$sex) # sex frequencies
CrossTable(micronutDBsa$sex, micronutDBsa$arm, chisq = TRUE) # crosstab sex x arm

table(micronutDBsa$ses3) # SES tertile group frequencies
CrossTable(micronutDBsa$ses3, micronutDBsa$arm, chisq = TRUE) # crosstab SES tertile x arm

table(micronutDBsa$bmigp) # BMI group frequencies
CrossTable(micronutDBsa$bmigp, micronutDBsa$arm, chisq = TRUE) # crosstab BMI group x arm

table(micronutDBsa$cd4100) # CD4 binary group frequencies
CrossTable(micronutDBsa$cd4100, micronutDBsa$arm, chisq = TRUE) # crosstab CD4 group x arm

table(micronutDBsa$hbgp2) # Hb binary group frequencies
CrossTable(micronutDBsa$hbgp2, micronutDBsa$arm, chisq = TRUE) # crosstab Hb group x arm

table(micronutDBsa$crpgp) # CRP group frequencies
CrossTable(micronutDBsa$crpgp, micronutDBsa$arm, chisq = TRUE) # crosstab CD4 group x arm

aggregate(micronutDBsa$fu, by = list(micronutDBsa$arm), FUN = sum) # FU per arm: 0(LNS) and 1 (LNS-VM)
aggregate(micronutDBsa$newtb, by = list(micronutDBsa$arm), FUN = sum) # TB incidence per arm

# LNS = 42435/365 = 116.3 // 28 TB 
LNSev <- 28
LNSt <- 116.3
poisson.test(LNSev, LNSt)$conf.int[1]*10
poisson.test(LNSev, LNSt)$conf.int[2]*10

# LNS-VM = 44517/365 = 121.9 // 21 TB
LNSVMev <- 21
LNSVMt <- 121.9
poisson.test(LNSVMev, LNSVMt)$conf.int[1]*10
poisson.test(LNSVMev, LNSVMt)$conf.int[2]*10

aggregate(micronutDBsa$fu, by = list(micronutDBsa$bmigp), FUN = sum) # FU per BMI: 0 (17-18.5), 1 (16-16.9), 2 (<16)
aggregate(micronutDBsa$newtb, by = list(micronutDBsa$bmigp), FUN = sum) # TB incidence per BMI group
# 17-18.5 = 42581/365 = 116.6 // 20
BMI1ev <- 20
BMI1t <- 116.6
poisson.test(BMI1ev, BMI1t)$conf.int[1]*10
poisson.test(BMI1ev, BMI1t)$conf.int[2]*10

# 16-16.9 = 21958/365 = 60.2 // 15
BMI2ev <- 15
BMI2t <- 60.2
poisson.test(BMI2ev, BMI2t)$conf.int[1]*10
poisson.test(BMI2ev, BMI2t)$conf.int[2]*10

# <16 = 22413/365 = 61.4 // 14
BMI3ev <- 14
BMI3t <- 61.4
poisson.test(BMI3ev, BMI3t)$conf.int[1]*10
poisson.test(BMI3ev, BMI3t)$conf.int[2]*10

# Survival analysis
survsa <- Surv(time = micronutDBsa$fu, event = micronutDBsa$newtb)

coxrawsa <- coxph(survsa ~ factor(micronutDBsa$arm)) # Cox regression per arm (raw)
summary(coxrawsa)

coxadjsa <- coxph(survsa ~ factor(micronutDBsa$arm) + factor(micronutDBsa$agegp2) + factor(micronutDBsa$sex) + 
                  factor(micronutDBsa$bmigp) + factor(micronutDBsa$cd4100) + 
                  factor(micronutDBsa$hbgp2) + factor(micronutDBsa$crpgp)) # Cox regression per arm (adjusted)
summary(coxadjsa)

# Fig S1A - Survival plot per arm (650w x 450h)
figs1a <- survfit(Surv(micronutDBsa$fu, micronutDBsa$newtb) ~ micronutDBsa$arm, data = micronutDBsa)

ggsurvplot(figs1a, ggtheme = theme_classic(), fun = "event", censor = FALSE, pval = FALSE,
           linetype = "strata", palette = c("#66C2A5", "#FC8D62"), conf.int = TRUE,
           xlab = "Time in study (months)", ylab = "Cumulative incident TB",
           legend = "bottom", legend.title = "Arm",  legend.labs = c("LNS", "LNS-VM"),
           risk.table = "abs_pct", fontsize = 2.25, ylim = c(0, 0.4), surv.scale = c("percent"),
           xscale = 30, break.x.by = 30, xlim = c(0,120), axes.offset = TRUE)

# BMI association with TB incidence ---
coxbmirawsa <- coxph(survsa ~ factor(micronutDBsa$bmigp)) # Cox regression per BMI (raw)
summary(coxbmirawsa)

coxbmiadjsa <- coxph(survsa ~ factor(micronutDBsa$bmigp) + factor(micronutDBsa$agegp2) + factor(micronutDBsa$sex) + 
                     factor(micronutDBsa$cd4100) + factor(micronutDBsa$hbgp2) + factor(micronutDBsa$crpgp)) # Cox regression per BMI (adjusted)
summary(coxbmiadjsa)

bmiestAsa <- coxph(survsa ~ factor(micronutDBsa$bmigp) + factor(micronutDBsa$agegp2) + factor(micronutDBsa$sex) + 
                   factor(micronutDBsa$cd4100) + factor(micronutDBsa$hbgp2) + factor(micronutDBsa$crpgp)) 
bmiestBsa <- coxph(survsa ~ factor(micronutDBsa$agegp2) + factor(micronutDBsa$sex) + 
                   factor(micronutDBsa$cd4100) + factor(micronutDBsa$hbgp2) + factor(micronutDBsa$crpgp)) 
lrtest(bmiestAsa, bmiestBsa)

# Fig S1B - Survival plot per arm (650w x 450h)
figs1b <- survfit(Surv(micronutDBsa$fu, micronutDBsa$newtb) ~ micronutDBsa$bmigp, data = micronutDBsa)

ggsurvplot(figs1b, ggtheme = theme_classic(), fun = "event", censor = FALSE, pval = FALSE,
           linetype = "strata", palette = c("#66C2A5", "#FC8D62", "#8DA0CB"), conf.int = TRUE,
           xlab = "Time in study (months)", ylab = "Cumulative incident TB",
           legend = "bottom", legend.title = "BMI",  legend.labs = c("17-18.5", "16-16.9", "<16"),
           risk.table = "abs_pct", fontsize = 2.25, ylim = c(0, 0.4), surv.scale = c("percent"),
           xscale = 30, break.x.by = 30, xlim = c(0,120), axes.offset = TRUE)

table(micronutDBsa2$xreason1, micronutDBsa2$arm)