clear

import delimited "/Users/rachaelburke/Documents/2020 CASTLE/CASTLE database/CASTLE_data/castle_trial_sharing.csv"

** Make strings into factors in stata
encode a310_tbsus, generate(tb_sus)
encode z999_cd4under100, generate(cd4_u100)


*** Primary outcome; glm w. robust SEs ***
tab e999_treatment b004_arm, col
glm e999_treatment b004_arm, family(binomial) link(log) vce(cluster cluster_n) eform
glm e999_treatment b004_arm, family(binomial) link(identity) vce(cluster cluster_n)

*** Primary outcome by subgroup; glm w. robust SEs ***
table e999_treatment b004_arm, by (tb_sus) col 

glm e999_treatment b004_arm##i.tb_sus, link(log) vce(cluster cluster_n) eform
glm e999_treatment b004_arm##i.tb_sus, link(identity) vce(cluster cluster_n)
* ^ Coded so that "tb_sus"==no is baseline

glm e999_treatment b004_arm##ib2.tb_sus, link(log) vce(cluster cluster_n) eform
glm e999_treatment b004_arm##ib2.tb_sus, link(identity) vce(cluster cluster_n)
* ^ Coded so that "tb_sus"==yes is baseline

table e999_treatment b004_arm, by (cd4_u100) col 
glm e999_treatment b004_arm##i.cd4_u100, link(log) vce(cluster cluster_n) eform
glm e999_treatment b004_arm##i.cd4_u100, link(identity) vce(cluster cluster_n)

glm e999_treatment b004_arm##ib2.cd4_u100, link(log) vce(cluster cluster_n) eform
glm e999_treatment b004_arm##ib2.cd4_u100, link(identity) vce(cluster cluster_n)

glm e999_treatment b004_arm##ib3.cd4_u100, link(log) vce(cluster cluster_n) eform
glm e999_treatment b004_arm##ib3.cd4_u100, link(identity) vce(cluster cluster_n)

testparm 1.b004_arm##2.cd4 1.b004_arm##3.cd4

*** Secondary outcomes

*** 1. TB treatment within 24 hours (same analysis approach as primary outcome)
tab z999_tbrx24h b004_arm, col
glm z999_tbrx24h_num b004_arm, family(binomial) link(log) vce(cluster cluster_n) eform
glm  z999_tbrx24h_num b004_arm, family(binomial) link(identity) vce(cluster cluster_n)

*** 2. Death up to 56 days (hazard ratio)
tab z999_outcome b004_arm, col
tab z999_failure b004_arm, col
stset z999_time2out_num, failure(z999_failure==1)
stcox b004_arm, vce(cluster cluster_n)
sts graph, by(b004_arm) fail risktable(,failevents) xlabel(0(14)56)
* ^ graph looks same as one drawn in R (for publication)glm e999_treatment b004_arm, family(binomial) link(identity) vce(cluster cluster_n)


*** 3. Undiagnosed TB at discharge (no events in one arm)
tab z999_undiag b004_arm
glm z999_undiag b004_arm, family(binomial) link(log) vce(cluster cluster_n) eform
glm z999_undiag b004_arm, family(binomial) link(identity) vce(cluster cluster_n)
* ^ code for glm for undiagnosed TB, not meaningful as none in one arm

csi 2 0 206 207, exact
* ^ use csi instead to get p value for difference

**** ICC / amount of clustering (as per CONSORT cluster RCT extension)
melogit e999_treatment b004_arm || cluster_n:, eform
melogit e999_treatment || cluster_n:


* ^ NB. logit link not log link, so can't compare to glm with log-binomial function
estat icc

melogit e999_treatment if b004_arm==1 || cluster_n:
estat icc