rm(list=ls())
# ref is the reference level for HAART time period.  Values can be 'noHAART' or 'gt180'.
# If 'noHAART' is used, then the marginal structural model uses total_HAART_indicator;
# if 'gt180' is used, then the marginal structural model uses indicator_haart_month and
# indHAART180I_month
ref <- "noHAART" 
#-------------------------------#
#       Loading libraries       #
#-------------------------------#
library(Design)
library(sandwich) # necessary to calculate appropriate SEs

#-------------------------------#
#        Loading data           #
#-------------------------------#
# Data formatted as month-to-month
load("timedfImonth.rda")


#-------------------------------#
#     Calculating weights       #
#-------------------------------#
# Definition of variables in the model:
# lag_indicator_haart_month = an indicator for whether the subject had taken HAART lagged by a month
# combined_term_month = (CD4 + lag_indicator_haart_month*CD4)*no_na_ind_cd4_month
#    where no_na_ind_cd4_month = 0, if CD4 at time t is missing; 1, if otherwise.
# na_ind_cd4_month = 0, if CD4 is non-missing at time t; 1, if otherwise.
wt.mod.denom.full <- glm(indicator_haart_month ~ rcs(month,3) + lag_indicator_haart_month + 
                                      combined_term_month + na_ind_cd4_month + age_at_first_visit + 
                                      sex + ivdu + birth_country.factor + 
                                      race.factor, data=time.df.month, family="binomial")
denom.full <- exp(predict(wt.mod.denom.full))/(1+exp(predict(wt.mod.denom.full)))

wt.mod.num.full <- glm(indicator_haart_month ~ rcs(month,3) + age_at_first_visit + sex + ivdu + birth_country.factor + 
                                               race.factor, data=time.df.month, family="binomial")
num.full <- exp(predict(wt.mod.num.full))/(1+exp(predict(wt.mod.num.full)))

time.df.month$wt.month.full <- ifelse(time.df.month$indicator_haart_month == 1, num.full/denom.full, (1-num.full)/(1-denom.full))

# Truncated weights (< 1% and > 99%)
quant.wt.full <- quantile(time.df.month$wt.month.full, probs=c(0.01,0.99))
time.df.month$wt3.month.full <- ifelse(time.df.month$wt.month.full < quant.wt.full[1], quant.wt.full[1],
                                ifelse(time.df.month$wt.month.full > quant.wt.full[2], quant.wt.full[2],
                                time.df.month$wt.month.full))

#--------------------------------#
# Logistic regression model with #
#  inverse probability weights   #
#--------------------------------#
if(ref=="noHAART"){
   # total_HAART_indicator is a 3-level categorical variable with levels -- no HAART, <= 180 days of HAART, and > 180 days of HAART
   msm.mod4.month.trunc.full <- glm(tb_month ~ rcs(month,3) + total_HAART_indicator, data=time.df.month, 
                                    weights=wt3.month.full,family="binomial")
} else {
   # indicator_haart_month = 1, if took HAART in month; 0, otherwise.
   # indHAART180I_month = 1, if in the first 180 days of a HAART regimen; 0, if more than 180 days in HAART regimen.
   msm.mod4.month.trunc.full <- glm(tb_month ~ rcs(month,3) + indicator_haart_month + indHAART180I_month, data=time.df.month, 
                                    weights=wt3.month.full,family="binomial")
}

# Computing appropriate SE given repeated measures on individuals
sand.var.mod4.month.trunc.full <- sandwich(msm.mod4.month.trunc.full)


#--------------------------------#
#   Formatting data for output   #
#--------------------------------#
# If embedding this code in a Sweave (*.nw) file, can use the latex function to create a pdf of the data frame created
# below.
month.I.trunc2.coeff.full <- msm.mod4.month.trunc.full$coeff[-c(1:3)]
month.I.trunc2.OR.full <- exp(msm.mod4.month.trunc.full$coeff[-c(1:3)])
month.I.trunc2.SE.full <- sqrt(diag(sand.var.mod4.month.trunc.full)[-c(1:3)])
month.I.trunc2.lci.full <- exp(month.I.trunc2.coeff.full - qnorm(0.975)*month.I.trunc2.SE.full)
month.I.trunc2.uci.full <- exp(month.I.trunc2.coeff.full + qnorm(0.975)*month.I.trunc2.SE.full)
month.I.trunc2.ci.full <- paste("(",round(month.I.trunc2.lci.full,2),", ",round(month.I.trunc2.uci.full,2),")",sep="")
wald.stat.full <- month.I.trunc2.coeff.full/month.I.trunc2.SE.full
month.I.trunc2.p.full <- 1 - pchisq(wald.stat.full^2,1)
month.I.trunc2.p.full <- ifelse(month.I.trunc2.p.full < 0.0001,"< 0.0001",
                          ifelse(month.I.trunc2.p.full > 0.0001 & month.I.trunc2.p.full < 0.001,round(month.I.trunc2.p.full,3),
                           round(month.I.trunc2.p.full,2)))
df.trunc2.full <- data.frame("OR"=round(month.I.trunc2.OR.full,2),
                             "CI"=month.I.trunc2.ci.full,
                             "P"=month.I.trunc2.p.full)