---
title: "ANALYSIS OF THE 28CN TRIAL"
author: "Nhat Le"
date: "Aug 7, 2018"
output: html_document
---
```{r, include=FALSE}
rm(list=ls()) # clear all
library(knitr)
opts_chunk$set(comment="",cache=FALSE,message=FALSE,warning=FALSE,echo=FALSE)
dir <- "V:/BIOSTATISTICS/RCT/CN_CNS_HIV/28CN_Tamoxifen_RCT"

dummy <- F # set to true (T) if the randomisation list is just a dummy list, otherwise false (F)

randolist_loc <- "Data/RandomizationList/28CN_Randlist.csv" # location of randomization list 
#randolist_loc <- "Data/RandomizationList/dummyRandomizationList.csv" # location of randomization list 
source(paste(dir,"Analysis/Ranalysis/Rfunctions/miscRfunctions28cn.R",sep="/"))

##----------------------  Load all libraries  ----------------------## 
library(Hmisc)
library(readr)
library(xtable)
library(survival)
library(lme4)
library(ggplot2)
library(pander)
library(rmarkdown)
library(plyr)
library(mitools)
library(compareGroups)
library(lmec)
require(JMbayes)
require(rstan)
library(kableExtra)
library(magrittr)
options(width=250)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
set.seed(123456)

##---------------------- Import data, merge randomization code ---------------------- 

#+ results='hide'
## Data
tmp <- sasxport.get(file=paste(dir,"Data/sasvadcsv/",sep="/"),method="csv",as.is=F)
#tmp <- lapply(tmp,function(x){x$usubjid <- as.character(x$usubjid);x}) # usubjid must be a character
for (i in 1:length(names(tmp))){
  assign(names(tmp[i]),tmp[[names(tmp[i])]])
}
rm(i,tmp)

## Randomization list
randolist <- read.csv(file=paste(dir,randolist_loc,sep="/"),as.is=T)

### We coorect information for two patients
### Patients 003-041 patient die at day 8 not censored
### Patient 003-037 patient removed from pp group because he didn't use drug since he was too severe.

# endpoint$evdeath.week10[endpoint$usubjid=="003-041"] <-1
# endpoint$evdeath[endpoint$usubjid=="003-041"] <-1
# endpoint$evtypetxt.irisdeath[endpoint$usubjid=="003-041"] <-1
# endpoint$evtypetxt.irisdeath[endpoint$usubjid=="003-041"] <-1



if (dummy) randolist$arm <- factor(randolist$arm,levels=c("Dummy Group A","Dummy Group B"))
if (!dummy){
  randolist<-reshape::rename(randolist,c("pat.id"="rando.number","r.arm"="arm"))
  randolist$arm <- factor(randolist$arm,levels=c("Control (Fluconazole)","Intervention (Tamoxifen + Fluconazole)"),labels = c("Control","Tamox"))
} 
lev.arm<-levels(randolist$arm)

#remove the arm column in population dataset
population<-population[,-9]

## Create minimal dataset with patient id, randomization arm, bl.date, and population info
pop <- merge(randolist,subset(bl.clinical,select=c(usubjid,rando.number,bl.date)),by="rando.number")
pop <- pop[,c("usubjid","arm","bl.date")]
pop <- merge(pop,subset(population,select=c(usubjid,itt,pp)),by="usubjid")
pop.itt <- subset(pop,itt==1)
bl.clinical$bl.date <- NULL # drop bl.date from bl.clinical (as it's now in pop)
bl.ecg$qtc.bl<-as.numeric(bl.ecg$qtc.bl)
ecg$qtc<-as.numeric(ecg$qtc)


## Order selected factors, convert selected factors to character etc.
ae$ae.name <- tolower(as.character(ae$ae.name))
ae$ae.soc <- toupper(as.character(ae$ae.soc))
labae$ae <- as.character(labae$ae)
labae$ae.code <- as.numeric(factor(labae$ae,levels=c(
  'Haemoglobin - low','White cell count - low','Neutrophils - low','Platelets - low',
  'ALT - high','AST - high','BlGlucose - high','BlGlucose - low','Creatinine - high',
  'Potassium - high','Potassium - low','Sodium - high','Sodium - low')))


levels.vis <- c("Normal for this patient","Blurred","Finger counting","Movement perception","Light perception",
                "No light perception","Unable to assess")
bl.clinical$visual.acuity.worst <- factor(bl.clinical$visual.acuity.worst,levels=levels.vis)
bl.clinical$visual.acuity.left <- factor(bl.clinical$visual.acuity.left,levels=levels.vis)
bl.clinical$visual.acuity.right <- factor(bl.clinical$visual.acuity.right,levels=levels.vis)
endpoint$visual.acuity.worst.week10 <- factor(endpoint$visual.acuity.worst.week10,levels=levels.vis)
endpoint$visual.acuity.left.week10 <- factor(endpoint$visual.acuity.left.week10,levels=levels.vis)
endpoint$visual.acuity.right.week10 <- factor(endpoint$visual.acuity.right.week10,levels=levels.vis)

levels.fundos <- c("Normal","Abnormal","Unable to visualize fundos","Not done")
bl.clinical$fundoscopy <- factor(bl.clinical$fundoscopy,levels=levels.fundos)
endpoint$fundoscopy.week10 <- factor(endpoint$fundoscopy.week10,levels=levels.fundos)

levels.disability <- c("Good","Intermediate","Severe disability","Death")
endpoint$disability.week10 <- factor(endpoint$disability.week10,levels=levels.disability)
endpoint$disability.rankin.week10 <- factor(endpoint$disability.rankin.week10,levels=levels.disability)
endpoint$disability.simpleq.week10 <- factor(endpoint$disability.simpleq.week10,levels=levels.disability)


```

## EXCLUSIONS FROM THE ITT AND PER-PROTOCOL POPULATION, STUDY DRUG AND STUDY COMPLETIONS
### Exclusions from the ITT

```{r echo=FALSE, message=FALSE, warning=FALSE, paged.print=TRUE}
tmp <- merge(subset(pop,select=c(usubjid,arm)),population,by="usubjid")
tmp <- merge(tmp,endpoint,by="usubjid")
tab <-mySummary.allvar(blvars=tmp[,3,drop=F],group=tmp$arm,contSummary="med.IQR",pval.comparison=F)
colnames(tab)<-c(tab[2,])
rownames(tab)<-NULL
kable(as.data.frame(tab[-2,])) %>% kable_styling("striped", full_width = F)
```



### Exclusions from the per-protocol population, completion of study drug and study (ITT population)
```{r echo=FALSE, paged.print=TRUE}
rm(pop) # as only pop.itt is needed subsequently
tab <-mySummary.allvar(blvars=tmp[tmp$itt==1,c(4:7,10,12,22,24)],group=tmp$arm[tmp$itt==1],contSummary="med.IQR",pval.comparison=F)
colnames(tab)<-c(tab[2,])
rownames(tab)<-NULL
kable(as.data.frame(tab[-2,])) %>% kable_styling("striped", full_width = F)
```



## BASELINE PATIENT CHARACTERISTICS AT STUDY ENTRY 


```{r, echo=FALSE}
##---------------------- PREPARE DATASET d USED FOR MOST EFFICACY ANALYSES  ---------------------- 
# Create dataset with all covariates used for subgroup analyses
cov1 <- subset(bl.clinical, select=c(usubjid,hiv.status,sex,age,resp.exam,resp.abnorm,gcs.category,visual.acuity.worst,fundoscopy))
cov2 <- subset(bl.imaging, select=c(usubjid,scan.cryptococcoma))
cov3 <- subset(bl.lp,select=c(usubjid,log10.csf.fungalcount.bl,openpres.bl,csf.wcc.bl))
cov4 <- subset(bl.lab,select=c(usubjid,cd4.bl))
cov5 <- subset(arv,select=c(usubjid,arv.status.bl,duration.arv.bl))
cov6 <- subset(bl.ecg,select=c(usubjid,qtc.bl,qtc.cat.bl))

cov.all <- merge(pop.itt,cov1,by="usubjid")
cov.all <- merge(cov.all,cov2,by="usubjid"); cov.all <- merge(cov.all,cov3,by="usubjid") 
cov.all <- merge(cov.all,cov4,by="usubjid"); cov.all <- merge(cov.all,cov5,by="usubjid");cov.all <- merge(cov.all,cov6,by="usubjid")
# tab<-compareGroups(arm~age+sex+gcs.category+log10.csf.fungalcount.bl+qtc.bl+qtc.cat.bl,data=d,method=2,show.p.overall = TRUE)
# createTable(tab)

tab <-mySummary.allvar(blvars=cov.all[c(8,7,11,15,21,22)],group=cov.all$arm,contSummary="med.IQR",pval.comparison=F)
colnames(tab)<-c(tab[2,])
rownames(tab)<-NULL
kable(as.data.frame(tab[-2,]))%>%
  kable_styling("striped", full_width = TRUE)
```




```{r, include=FALSE}

##---------------------- PREPARE DATASET d USED FOR MOST EFFICACY ANALYSES  ---------------------- 
# Create dataset with all covariates used for subgroup analyses
cov1 <- subset(bl.clinical, select=c(usubjid,hiv.status,sex,age,resp.exam,resp.abnorm,gcs,visual.acuity.worst,fundoscopy))
cov2 <- subset(bl.imaging, select=c(usubjid,scan.cryptococcoma))
cov3 <- subset(bl.lp,select=c(usubjid,log10.csf.fungalcount.bl,openpres.bl,csf.wcc.bl))
cov4 <- subset(bl.lab,select=c(usubjid,cd4.bl))
cov5 <- subset(arv,select=c(usubjid,arv.status.bl,duration.arv.bl))
cov6 <- subset(bl.ecg,select=c(usubjid,qtc.bl,qtc.cat.bl))

cov.all <- merge(pop.itt,cov1,by="usubjid")
cov.all <- merge(cov.all,cov2,by="usubjid"); cov.all <- merge(cov.all,cov3,by="usubjid") 
cov.all <- merge(cov.all,cov4,by="usubjid"); cov.all <- merge(cov.all,cov5,by="usubjid");cov.all <- merge(cov.all,cov6,by="usubjid")
cov.all$gcs.category <- factor(ifelse(cov.all$gcs<15,"<15","15"),levels=c("15","<15"))
cov.all$age.category <- factor(ifelse(cov.all$age<=35,"<=35",">35"),levels=c("<=35",">35"))
cov.all$fungalCount.category <- factor(ifelse(cov.all$log10.csf.fungalcount.bl<=5,"<=10^5",">10^5"),levels=c("<=10^5",">10^5"))
cov.all$cd4.category <- factor(ifelse(cov.all$cd4<=25,"<=25",">25"),levels=c("<=25",">25"))
cov.all$randomized.up.to.21Jul2014 <- factor(ifelse(cov.all$bl.date<=as.Date("21jul2014", "%d%b%Y"),"Yes","No"))
cov.all$openpres.category <- factor(ifelse(cov.all$openpres.bl>18,">18","<=18"),levels=c("<=18",">18"))
cov.all$csf.wcc.category <- factor(ifelse(cov.all$csf.wcc.bl<5,"<5",">=5"),levels=c("<5",">=5"))

d <- merge(cov.all,endpoint,by="usubjid")
d <- merge(d,subset(population,select=c(usubjid,ttdeath.pp,evdeath.pp,ttdeath.week10.pp,evdeath.week10.pp)),by="usubjid")
```




## RATE OF CSF STERILISATION DURING THE FIRST 2 WEEKS (PRIMARY ENDPOINT) 
##### Graphical display (including days 0-17, as for the modeling)

```{r, echo=FALSE, fig.height=4, fig.width=8}
#load fitted model
load(file="V:/BIOSTATISTICS/RCT/CN_CNS_HIV/28CN_Tamoxifen_RCT/Analysis/Ranalysis/Results/unbias.Rdata")
fc <- merge(pop.itt,subset(lp,(!is.na(csf.fungalcount))&(day>=0)&(day<=17),select=c(usubjid,day,csf.fungalcount,log10.csf.fungalcount)),by="usubjid")

fc <- merge(fc,subset(bl.lp,select=c(usubjid,log10.csf.fungalcount.bl)),by="usubjid",all.x=T)
fc <-merge(fc,subset(bl.clinical,select=c(usubjid,hiv.status)),by="usubjid",all.x=T)

dat<-data.frame(day=rep(0:17,each=2),arm=rep(levels(fc$arm),by=14))
mm<-model.matrix(~day+day:arm,dat)
pvar1 <- diag(mm %*% tcrossprod(vcov_fit(declines.all$obj),mm))
dat$pred<-mm%*%c(summary(declines.all$obj)$summary[1:3,"mean"]) #predict(m,newdat,re.form=NA) would give the same results
dat$plo<-dat$pred-qnorm(0.975,0,1)*sqrt(pvar1)
dat$phi <- dat$pred+qnorm(0.975,0,1)*sqrt(pvar1)
dat$arm <- factor(dat$arm, levels=c("Control","Tamox"))

ggplot(data=dat,aes(x=day,y=pred))+facet_grid(.~arm)+
  geom_line(data=fc, aes(day,pmax(log10.csf.fungalcount,log10(4.5)),group =usubjid),color=grey(0.8),size=0.6)+geom_point(data=fc, aes(day,pmax(log10.csf.fungalcount,log10(4.5))),color="darkgray",size=1)+
  geom_line(colour="#3182bd",size=1)+geom_ribbon(aes(x=day,ymin=plo,ymax=phi),fill="#3182bd",alpha=0.3)+
  scale_y_continuous("CSF Quantitative fungal count  [CFU/ml]",
                     breaks=c(log10(4.5),1:8),
                     labels=c("4.5","10",expression(10^2),expression(10^3),expression(10^4),expression(10^5),expression(10^6),expression(10^7),expression(10^8)))+
  scale_x_continuous("Study day",limits = c(0,17),breaks=c(1,8,15),minor_breaks=0:17)+
  theme_bw()+
  theme(panel.background=element_rect(fill = 'white', colour = 'grey'),
        axis.title.x = element_text(face="bold", colour="black", size=12, margin = margin(t = 20,r = 20,b =20, l = 20)),
        axis.title.y = element_text(face="bold", colour="black", size=12, margin = margin(t = 20,r = 20,b = 20, l = 20)),
        axis.text.y = element_text(face="plain", colour="black", size=10),
        axis.text.x = element_text(face="plain", colour="black", size=10),
        axis.ticks.x=element_blank(),
        legend.text = element_text(face="plain",size = 12),
        legend.title = element_blank(),
        legend.key = element_rect(colour = NA),
        legend.position = "right",
        strip.text=element_text(face = "bold",size=12))+coord_cartesian(ylim=c(log10(4.5),7))
ggsave(file = "V:/BIOSTATISTICS/RCT/CN_CNS_HIV/28CN_Tamoxifen_RCT/Analysis/Ranalysis/Figures/Figure1.png",   
    width = 7, 
    height = 5, dpi=400) 
```

```{r, echo=FALSE, fig.height=4, fig.width=8}

p <- ggplot(data=subset(fc,pp==1), aes(day,pmax(csf.fungalcount,1),group =usubjid))

p+geom_line(color="gray")+geom_point(color="darkgray",size=1)+facet_grid(.~arm)+
  scale_y_log10("Quantitative fungal count of PP [CFU/ml]",
                breaks=10^c(0:8),
                labels=c("0","10",expression(10^2),expression(10^3),expression(10^4),expression(10^5),expression(10^6),expression(10^7),expression(10^8)),
                minor_breaks=5*(10^c(0:7)))+
  scale_x_continuous("Study day",limits = c(0,17),breaks=c(1,8,15),minor_breaks=0:17)+
  geom_smooth(aes(group=arm),method="loess",se=F,lwd=1.2)+ theme_bw()+
  theme(panel.background=element_rect(fill = 'white', colour = 'grey'),
    axis.title.x = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b =20, l = 20)),
    axis.title.y = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b = 20, l = 20)),
    axis.text.y = element_text(face="plain", colour="black", size=10),
    axis.text.x = element_text(face="plain", colour="black", size=10),
    axis.ticks.x=element_blank(),
    legend.text = element_text(face="plain",size = 15),
    legend.title = element_blank(),
    legend.key = element_rect(colour = NA),
    legend.position = "right",
    strip.text=element_text(face = "bold",size=15))

dat<-data.frame(day=rep(0:17,each=2),arm=rep(levels(fc$arm),by=14))
mm<-model.matrix(~day+day:arm,dat)
pvar1 <- diag(mm %*% tcrossprod(vcov_fit(declines.all$obj),mm))
dat$pred<-mm%*%c(summary(declines.all$obj)$summary[1:3,"mean"]) #predict(m,newdat,re.form=NA) would give the same results
dat$plo<-dat$pred-qnorm(0.975,0,1)*sqrt(pvar1)
dat$phi <- dat$pred+qnorm(0.975,0,1)*sqrt(pvar1)
dat$arm <- factor(dat$arm, levels=c("Control","Tamox"))

ggplot(data=dat,aes(x=day,y=pred))+facet_grid(.~arm)+
  geom_line(data=fc, aes(day,pmax(log10.csf.fungalcount,log10(4.5)),group =usubjid),color=grey(0.8),size=0.6)+geom_point(data=fc, aes(day,pmax(log10.csf.fungalcount,log10(4.5))),color="darkgray",size=1)+
  geom_line(colour="#3182bd",size=1)+geom_ribbon(aes(x=day,ymin=plo,ymax=phi),fill="#3182bd",alpha=0.3)+
  scale_y_continuous("CSF Quantitative fungal count  [CFU/ml]",
                     breaks=c(log10(4.5),1:8),
                     labels=c("4.5","10",expression(10^2),expression(10^3),expression(10^4),expression(10^5),expression(10^6),expression(10^7),expression(10^8)))+
  scale_x_continuous("Study day",limits = c(0,17),breaks=c(1,8,15),minor_breaks=0:17)+
  theme_bw()+
  theme(panel.background=element_rect(fill = 'white', colour = 'grey'),
        axis.title.x = element_text(face="bold", colour="black", size=12, margin = margin(t = 20,r = 20,b =20, l = 20)),
        axis.title.y = element_text(face="bold", colour="black", size=12, margin = margin(t = 20,r = 20,b = 20, l = 20)),
        axis.text.y = element_text(face="plain", colour="black", size=10),
        axis.text.x = element_text(face="plain", colour="black", size=10),
        axis.ticks.x=element_blank(),
        legend.text = element_text(face="plain",size = 12),
        legend.title = element_blank(),
        legend.key = element_rect(colour = NA),
        legend.position = "right",
        strip.text=element_text(face = "bold",size=12))+coord_cartesian(ylim=c(log10(4.5),7))
ggsave(file = "V:/BIOSTATISTICS/RCT/CN_CNS_HIV/28CN_Tamoxifen_RCT/Analysis/Ranalysis/Figures/Figure1.png",   
    width = 7, 
    height = 5, dpi=400) 









```

```{r, echo=FALSE, fig.height=4, fig.width=8}

p <- ggplot(data=subset(fc,!is.na(hiv.status)), aes(day,pmax(csf.fungalcount,1),group =usubjid))

p+geom_line(color="gray")+geom_point(color="darkgray",size=1)+facet_grid(.~hiv.status+arm)+
  scale_y_log10("Quantitative fungal count of ITT with HIV classification [CFU/ml]",
                breaks=10^c(0:8),
                labels=c("0","10",expression(10^2),expression(10^3),expression(10^4),expression(10^5),expression(10^6),expression(10^7),expression(10^8)),
                minor_breaks=5*(10^c(0:7)))+
  scale_x_continuous("Study day",limits = c(0,17),breaks=c(1,8,15),minor_breaks=0:17)+
  geom_smooth(aes(group=arm),method="loess",se=F,lwd=1.2)+ theme_bw()+
  theme(panel.background=element_rect(fill = 'white', colour = 'grey'),
    axis.title.x = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b =20, l = 20)),
    axis.title.y = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b = 20, l = 20)),
    axis.text.y = element_text(face="plain", colour="black", size=10),
    axis.text.x = element_text(face="plain", colour="black", size=10),
    axis.ticks.x=element_blank(),
    legend.text = element_text(face="plain",size = 15),
    legend.title = element_blank(),
    legend.key = element_rect(colour = NA),
    legend.position = "right",
    strip.text=element_text(face = "bold",size=15))
```

```{r, echo=FALSE, fig.height=4, fig.width=8}

p <- ggplot(data=subset(fc,log10.csf.fungalcount.bl<5), aes(day,pmax(csf.fungalcount,1),group =usubjid))

p+geom_line(color="gray")+geom_point(color="darkgray",size=1)+facet_grid(.~arm)+
  scale_y_log10("Quantitative fungal count of low fc burden [CFU/ml]",
                breaks=10^c(0:8),
                labels=c("0","10",expression(10^2),expression(10^3),expression(10^4),expression(10^5),expression(10^6),expression(10^7),expression(10^8)),
                minor_breaks=5*(10^c(0:7)))+
  scale_x_continuous("Study day",limits = c(0,17),breaks=c(1,8,15),minor_breaks=0:17)+
  geom_smooth(aes(group=arm),method="loess",se=F,lwd=1.2)+ theme_bw()+
  theme(panel.background=element_rect(fill = 'white', colour = 'grey'),
    axis.title.x = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b =20, l = 20)),
    axis.title.y = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b = 20, l = 20)),
    axis.text.y = element_text(face="plain", colour="black", size=10),
    axis.text.x = element_text(face="plain", colour="black", size=10),
    axis.ticks.x=element_blank(),
    legend.text = element_text(face="plain",size = 15),
    legend.title = element_blank(),
    legend.key = element_rect(colour = NA),
    legend.position = "right",
    strip.text=element_text(face = "bold",size=15))
```

```{r, echo=FALSE, fig.height=4, fig.width=8}

p <- ggplot(data=subset(fc,(log10.csf.fungalcount.bl>=5)&(log10.csf.fungalcount.bl<=6)), aes(day,pmax(csf.fungalcount,1),group =usubjid))

p+geom_line(color="gray")+geom_point(color="darkgray",size=1)+facet_grid(.~arm)+
  scale_y_log10("Quantitative fungal count of middle fc burden [CFU/ml]",
                breaks=10^c(0:8),
                labels=c("0","10",expression(10^2),expression(10^3),expression(10^4),expression(10^5),expression(10^6),expression(10^7),expression(10^8)),
                minor_breaks=5*(10^c(0:7)))+
  scale_x_continuous("Study day",limits = c(0,17),breaks=c(1,8,15),minor_breaks=0:17)+
  geom_smooth(aes(group=arm),method="loess",se=F,lwd=1.2)+ theme_bw()+
  theme(panel.background=element_rect(fill = 'white', colour = 'grey'),
    axis.title.x = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b =20, l = 20)),
    axis.title.y = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b = 20, l = 20)),
    axis.text.y = element_text(face="plain", colour="black", size=10),
    axis.text.x = element_text(face="plain", colour="black", size=10),
    axis.ticks.x=element_blank(),
    legend.text = element_text(face="plain",size = 15),
    legend.title = element_blank(),
    legend.key = element_rect(colour = NA),
    legend.position = "right",
    strip.text=element_text(face = "bold",size=15))
```

```{r, echo=FALSE, fig.height=4, fig.width=8}

p <- ggplot(data=subset(fc,log10.csf.fungalcount.bl>6), aes(day,pmax(csf.fungalcount,1),group =usubjid))

p+geom_line(color="gray")+geom_point(color="darkgray",size=1)+facet_grid(.~arm)+
  scale_y_log10("Quantitative fungal count of high fc burden [CFU/ml]",
                breaks=10^c(0:8),
                labels=c("0","10",expression(10^2),expression(10^3),expression(10^4),expression(10^5),expression(10^6),expression(10^7),expression(10^8)),
                minor_breaks=5*(10^c(0:7)))+
  scale_x_continuous("Study day",limits = c(0,17),breaks=c(1,8,15),minor_breaks=0:17)+
  geom_smooth(aes(group=arm),method="loess",se=F,lwd=1.2)+ theme_bw()+
  theme(panel.background=element_rect(fill = 'white', colour = 'grey'),
    axis.title.x = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b =20, l = 20)),
    axis.title.y = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b = 20, l = 20)),
    axis.text.y = element_text(face="plain", colour="black", size=10),
    axis.text.x = element_text(face="plain", colour="black", size=10),
    axis.ticks.x=element_blank(),
    legend.text = element_text(face="plain",size = 15),
    legend.title = element_blank(),
    legend.key = element_rect(colour = NA),
    legend.position = "right",
    strip.text=element_text(face = "bold",size=15))
```



* There was 1 patient in control arm was excluded from ITT because she has 4 days of amphotericin B antifungal therapy after randomization for reasons other than death and 1 patient in Tamox arm was excluded from ITT because he had transaministis and didn't receive Tamox according to PI's decision. There was 1 patient in Tamox arm has no CSF measurement after beginning tamoxifen and died at the day 4th after randomization, but this patient didn't violate protocol.





### Unbias analysis based on joint modeling with linear mixed model for fungal count including detection limit and cox model for time to death outcome
```{r, include=FALSE}
fc<-merge(fc,subset(bl.lp,select=c("usubjid")),by="usubjid",all.y=T)
id.missing<-fc$usubjid[is.na(fc$csf.fungalcount)]
fc[fc$usubjid==id.missing,c("usubjid","day","csf.fungalcount","log10.csf.fungalcount")]<-lp[lp$usubjid==id.missing,c("usubjid","day","csf.fungalcount","log10.csf.fungalcount")]
fc[fc$usubjid==id.missing,c("arm","bl.date","itt","pp")]<-pop.itt[pop.itt$usubjid==id.missing,c("arm","bl.date","itt","pp")]
fc[fc$usubjid==id.missing,"log10.csf.fungalcount.bl"]<-bl.lp[bl.lp$usubjid==id.missing,c("log10.csf.fungalcount.bl")]
fc[fc$usubjid==id.missing,"hiv.status"]<-bl.clinical[bl.clinical$usubjid==id.missing,c("hiv.status")]


fc$studyday <- pmax(fc$day-1,0)
fc$arm1.studyday <- fc$studyday*(fc$arm==lev.arm[1])
fc$arm2.studyday <- fc$studyday*(fc$arm==lev.arm[2])

write.csv(fc,file="V:/BIOSTATISTICS/RCT/CN_CNS_HIV/28CN_Tamoxifen_RCT/Data/fc.csv")
clinical<-merge(endpoint,subset(pop.itt,select=c(usubjid,arm,pp)),by="usubjid")
clinical <- merge(clinical,subset(bl.lp,select=c(usubjid,log10.csf.fungalcount.bl)),by="usubjid",all.x=T)
clinical <-merge(clinical,subset(bl.clinical,select=c(usubjid,hiv.status)),by="usubjid",all.x=T)


declines.all <- declineEstimates.unbias(fc,clinical)

declines.pp <- declineEstimates.unbias(subset(fc,pp==1),subset(clinical,pp==1))

declines.hiv.pos <- declineEstimates.unbias(subset(fc,hiv.status=="Positive"),subset(clinical,hiv.status=="Positive"))
declines.hiv.neg <- declineEstimates.unbias(subset(fc,hiv.status=="Negative"),subset(clinical,hiv.status=="Negative"),timepoint=c(0,17),n.adapt=3000,n.iter=6000)

declines.low <- declineEstimates.unbias(subset(fc,log10.csf.fungalcount.bl<5),subset(clinical,log10.csf.fungalcount.bl<5),timepoint=c(0,7,10,13,17),adapt_delta = 0.98,n.adapt=4000,n.iter=6000)
declines.mid <- declineEstimates.unbias(subset(fc,(log10.csf.fungalcount.bl>=5)&(log10.csf.fungalcount.bl<=6)),subset(clinical,(log10.csf.fungalcount.bl>=5)&(log10.csf.fungalcount.bl<=6)),timepoint=c(0,17),adapt_delta = 0.9851,4000,n.iter=6000)
declines.high <- declineEstimates.unbias(subset(fc,log10.csf.fungalcount.bl>6),subset(clinical,log10.csf.fungalcount.bl>6),timepoint=c(0,17),adapt_delta = 0.98,4000,n.iter=6000)
save(declines.all,declines.pp,declines.hiv.pos,declines.hiv.neg,declines.low,declines.mid,declines.high,file="V:/BIOSTATISTICS/RCT/CN_CNS_HIV/28CN_Tamoxifen_RCT/Analysis/Ranalysis/Results/unbias.Rdata")

load(file="V:/BIOSTATISTICS/RCT/CN_CNS_HIV/28CN_Tamoxifen_RCT/Analysis/Ranalysis/Results/unbias.Rdata")
fungal.tab<-matrix("",nrow=8,ncol=7)
fungal.tab[1,]<-c("Subgroup","n",paste(lev.arm[1],"Change/day (95% CI)",sep="-"),"n",paste(lev.arm[2],"Change/day (95% CI)",sep="-"),"Difference in change (95% CI); p-value","p-adjusted")
fungal.tab[2,1:6]<-c("ITT population",unlist(declines.all))
fungal.tab[3,1:6]<-c("Per Protocol population",unlist(declines.pp))
fungal.tab[4,1:6]<-c("HIV positive patients only",unlist(declines.hiv.pos))
fungal.tab[5,1:6]<-c("HIV negative patients only",unlist(declines.hiv.neg))
fungal.tab[6,1:6]<-c("Patients with baseline count <5 log10 CFU/ml",unlist(declines.low))
fungal.tab[7,1:6]<-c("Patients with baseline count 5-6 log10 CFU/ml",unlist(declines.mid))
fungal.tab[8,1:6]<-c("Patients with baseline count >6 log10 CFU/ml",unlist(declines.high))
p.value<-c(0.9532,0.9577,0.4078,0.3712,0.8185,0.0141,0.0007)
```


```{r, echo=FALSE}
colnames(fungal.tab) <- fungal.tab[1,]; fungal.tab <- fungal.tab[-1,]
fungal.tab<-as.matrix(fungal.tab)
fungal.tab[,7]<-c(p.adjust(p.value,"hochberg"))
rownames(tab)<-NULL
kable(fungal.tab)%>%
  kable_styling("striped", full_width = TRUE)
```

*Of note, only subgroup "Patients with baseline count >6 log10 CFU/ml" has the significant beneficial effect after adjusted for multiple testing. But the sample in the Tamox group is only one.*

###   SURVIVAL UNTIL 10 WEEKS (SECONDARY ENDPOINT)
####  KAPLAN-MEIER CURVES OF OVERALL SURVIVAL AND ADDITIONAL EXPLORATORY FIGURES
##### Kaplan-Meier curves

```{r, echo=FALSE, fig.height=8.5, fig.width=10}
mykmplot(Surv(ttdeath.week10,evdeath.week10)~1,d,ylab="Time since randomization [days]",xlab="Survival probability",main="All patients")
```


```{r, echo=FALSE, fig.height=8.5, fig.width=10}
mykmplot(Surv(ttdeath.week10,evdeath.week10)~1,subset(d,hiv.status=="Positive"),ylab="Time since randomization [days]",xlab="Survival probability",main="Infected HIV patients")
```


```{r, echo=FALSE, fig.height=8.5, fig.width=10}
mykmplot(Surv(ttdeath.week10,evdeath.week10)~1,subset(d,hiv.status=="Negative"),ylab="Time since randomization [days]",xlab="Survival probability",main="Un-infected HIV patients")
```


#### 10-week mortality overall.

```{r, echo=FALSE}
tmp <- summary(survfit(Surv(ttdeath.week10,evdeath.week10)~arm,data=d),times=71,extend=T)
print(data.frame(group=rownames(tmp$table),day=tmp$time,mortality=1-tmp$surv,ci.lower=1-tmp$upper,ci.upper=1-tmp$lower),row.names=F)
```

#### 10-week mortality by HIV.

```{r, echo=FALSE}
tmp <- summary(survfit(Surv(ttdeath,evdeath)~hiv.status+arm,data=d),times=71,extend=T)
print(data.frame(group=rownames(tmp$table),day=tmp$time,mortality=1-tmp$surv,ci.lower=1-tmp$upper,ci.upper=1-tmp$lower),row.names=F)

```

### Exploratory figure - differences in mortality

The figures below show differences in mortality between the two groups (black lines), estimates+/-standard error (dark gray areas),and point-wise confidence intervals (light gray areas).

```{r, include=FALSE}
my.exploratory.survdiffplot <- function(data,main){
  tmp.1 <- summary(survfit(Surv(ttdeath,evdeath)~arm,data=subset(data,arm==levels(data$arm)[1])),times=1:72,extend=T)
  tmp.2 <- summary(survfit(Surv(ttdeath,evdeath)~arm,data=subset(data,arm==levels(data$arm)[2])),times=1:72,extend=T)
  diff <- data.frame(t=1:72,diff.surv=tmp.1$surv-tmp.2$surv,se.diff.surv=sqrt(tmp.1$std.err^2+tmp.2$std.err^2))
  plot(diff$t,100*diff$diff.surv,type="l",
       ylim=c(-25,25),ylab="Difference in mortality risks [%]",xlab="Days since randomization",axes=F,main=main)
  axis(1,at=seq(0,270,by=14),labels=NA)
  tp <- c(0,14,42,70)
  axis(1,at=tp)
  axis(2,at=seq(-30,30,by=10))
  tt <- c(1:72,72:1); 
  y1 <- 100*c(diff$diff.surv-diff$se.diff.surv,rev(diff$diff.surv+diff$se.diff.surv))
  y2 <- 100*c(diff$diff.surv-qnorm(0.975)*diff$se.diff.surv,rev(diff$diff.surv+qnorm(0.975)*diff$se.diff.surv))
  polygon(tt,y2,col="lightgray",border=NA)
  polygon(tt,y1,col="gray",border=NA)
  lines(diff$t,100*diff$diff.surv,type="l",ylim=c(-0.2,0.2),lwd=2)
  abline(h=0,lty=2)
}


```




```{r, echo=FALSE, fig.height=6, fig.width=18}
my.exploratory.survdiffplot(d,"All patients")

```

```{r, fig.height=6, fig.width=18}
my.exploratory.survdiffplot(subset(d,hiv.status=="Positive"),"HIV-Positive")
```

```{r, fig.height=6, fig.width=18}
my.exploratory.survdiffplot(subset(d,hiv.status=="Negative"),"HIV-Negative")
```

### Exploratory figure - hazards
The figures below shows estimated piecewise constant hazards in the two groups over time.

```{r, include=FALSE}
my.plothaz <- function(data,main){
  require(muhaz)
  d1 <- subset(data,arm==levels(data$arm)[1])
  d2 <- subset(data,arm==levels(data$arm)[2])
  sink("NUL")
  fit1 <- pehaz(d1$ttdeath,d1$evdeath,width=14) 
  fit2 <- pehaz(d2$ttdeath,d2$evdeath,width=14)
  sink() # sink is just to avoid annoying outputs in the html 
  plot(fit1,main=main,axes=F,ylab="Estimated hazard",xlab="Days since randomization")
  axis(1,at=seq(0,270,by=14),labels=NA)
  tp <- c(0,14,42,70,98,126,154,182)
  axis(1,at=tp)
  axis(2)
  lines(fit2,col=gray(0.7),lwd=2,lty="41")
  lines(fit1,lwd=2,lty="solid")
}
```

```{r, echo=FALSE, fig.height=6, fig.width=18}
my.plothaz(d,"All patients")
```

```{r, echo=FALSE, fig.height=6, fig.width=18}
my.plothaz(subset(d,hiv.status=="Positive"),"HIV-positive")
```

```{r, echo=FALSE, fig.height=6, fig.width=18}
my.plothaz(subset(d,hiv.status=="Negative"),"HIV-negative")
```



```{r, echo=FALSE}
#' ## SURVIVAL UNTIL 10 WEEKS (PRIMARY ENDPOINT) 
#' ### Overall analysis and pre-defined subgroups
#' Note: Displayed risks are based on the Kaplan-Meier method.

# base model for all analyses except for the per protocol and the analysis by HIV status
base.model <- Surv(ttdeath.week10,evdeath.week10)~arm+strata(hiv.status) 
# ITT, per protocol analysis
r1 <- surv.comparison(model=base.model,data=d)
r2 <- surv.comparison(model=Surv(ttdeath.week10.pp,evdeath.week10.pp)~arm+strata(hiv.status),data=subset(d,pp==T))
r <- cbind(c("Subgroup","","ITT","Per protocol"),
           rbind(r1,r2[-c(1,2),]),
           c("Test for heterogeneity","p-value","",""))
# Add subgroups
# by HIV status (not stratified by HIV status)
s <- surv.comparison.subgroup(base.model=Surv(ttdeath.week10,evdeath.week10)~arm,subgroup.model=~hiv.status,data=d)
r <- rbind(r,s[-c(1:3),])
# # all others (stratified by HIV.status) 
# s <- surv.comparison.subgroup(base.model=base.model,subgroup.model=~gcs.category+arv.status.bl+sex+age.category+fungalCount.category+cd4.category+openpres.category+csf.wcc.category,data=d)
# r <- rbind(r,s[-c(1:3),])
# # output all
# colnames(r) <- r[1,]; r <- r[-1,]
# kable(r)

#' ### Formal comparison of risks at 10 weeks 
r <- cbind(c("","ITT","Per protocol","HIV-positive","HIV-negative"),
           rbind(surv.comparison.abs.risk(model=Surv(ttdeath.week10,evdeath.week10)~.,data=d,time=71),
                 surv.comparison.abs.risk(model=Surv(ttdeath.week10.pp,evdeath.week10.pp)~.,data=subset(d,pp==T),time=71)[2,],
                 surv.comparison.abs.risk(model=Surv(ttdeath.week10,evdeath.week10)~.,data=subset(d,hiv.status=="Positive"),time=71)[2,],
                 surv.comparison.abs.risk(model=Surv(ttdeath.week10,evdeath.week10)~.,data=subset(d,hiv.status=="Negative"),time=71)[2,]))
rownames(r)<-NULL
kable(r)%>%
  kable_styling("striped", full_width = TRUE)
```


### Cox regression 

<!-- ```{r, eval=FALSE, include=FALSE} -->
<!-- fit.imp <- by(all.data.imp.long,all.data.imp.long$.imp, function(x)  -->
<!--   coxph(Surv(ttdeath.week10,evdeath.week10)~ -->
<!--           arm+strata(hiv.status)+ -->
<!--           log10.csf.fungalcount.bl+I(gcs<15)+arv.status.bl,data=x)) -->
<!-- #+ results='hide' -->
<!-- sfit <- summary(MIcombine(fit.imp),logeffect=T) -->
<!-- tmp <- data.frame(HR=sfit$expresults,lower.CI=sfit[,"(lower"],upper.CI=sfit[,"upper)"], -->
<!--                   pval=format.pval(sfit[,"pVal"])) -->
<!-- rownames(tmp) <- rownames(sfit) -->
<!-- #+ output=TRUE -->
<!-- tmp -->
<!-- #+ output=FALSE -->
<!-- ``` -->


```{r, echo=FALSE}
fit <- coxph(Surv(ttdeath.week10,evdeath.week10)~arm+strata(hiv.status)+log10.csf.fungalcount.bl+gcs.category+arv.status.bl,data=d) 
#Only one value was missing so we dont perfomed the imputation data.
sfit <- summary(fit)
#+ output=TRUE
tab<-data.frame(HR=sfit$coefficients[,"exp(coef)"],lower.CI=sfit$conf.int[,"lower .95"],upper.CI=sfit$conf.int[,"upper .95"],
           pval=format.pval(sfit$coefficients[,"Pr(>|z|)"]))
#+ output=FALSE

#' ### Additional exploratory analysis: Hazard ratios during days 1-22, 23-43, and 44-71
#' This splits the time axis into first half of study treatment period, second half, and remaining time.
# base.model <- Surv(ttdeath.week10,evdeath.week10)~arm+strata(hiv.status) 
# # ITT, per protocol analysis
# r1 <- surv.comparison(model=Surv(pmin(ttdeath.week10,22),ifelse(ttdeath.week10<=22,evdeath.week10,0))~arm+strata(hiv.status),data=d)
# r2 <- surv.comparison(model=Surv(pmin(ttdeath.week10,43),ifelse(ttdeath.week10<=43,evdeath.week10,0))~arm+strata(hiv.status),data=subset(d,ttdeath.week10>22))
# r3 <- surv.comparison(model=Surv(ttdeath.week10,evdeath.week10)~arm+strata(hiv.status),data=subset(d,ttdeath.week10>43))
# 
# r <- cbind(c("Time period","","Days 1-22","Days 23-43","Days 44-71"),
#            rbind(r1,r2[-c(1,2),],r3[-c(1,2),]))
# colnames(r) <- r[1,]; r <- r[-1,]
kable(tab)%>%
  kable_styling("striped", full_width = TRUE)
```


### DISABILITY AT 10 WEEKS
##### Descriptive table

```{r, echo=FALSE}
d.tmp <- d
d.tmp$disability.week10.hiv.pos <- d.tmp$disability.week10.hiv.neg <- d.tmp$disability.week10
d.tmp$disability.week10.hiv.pos[d.tmp$hiv.status!="Positive"] <- NA
d.tmp$disability.week10.hiv.neg[d.tmp$hiv.status!="Negative"] <- NA

d.tmp$disability.month6.hiv.pos <- d.tmp$disability.month6.hiv.neg <- d.tmp$disability.month6
d.tmp$disability.month6.hiv.pos[d.tmp$hiv.status!="Positive"] <- NA
d.tmp$disability.month6.hiv.neg[d.tmp$hiv.status!="Negative"] <- NA

tab <-mySummary.allvar(blvars=d.tmp[,c("disability.week10","disability.week10.hiv.pos","disability.week10.hiv.neg",
                                       "disability.rankin.week10","disability.simpleq.week10")],
                       group=d.tmp$arm,contSummary="med.IQR",pval.comparison=T)
colnames(tab)<-c(tab[2,])
rownames(tab)<-NULL
kable(as.data.frame(tab[-2,])) %>%
  kable_styling("striped", full_width = TRUE)

```



##### Overall and subgroup analyses for the probability of a "good" outcome excluding missing disability assessments

Of note, we cannot perfom this table when we dont have the ongoing status of all these patients
```{r echo=FALSE}
#base model for all analysis except for the analysis by hiv.status
base.model <- I(disability.week10=="Good")~arm+hiv.status
# ITT and per protocol analysis
r1 <- logist.comparison(model=base.model,data=d)
r2 <- logist.comparison(model=base.model,data=subset(d,pp==T))
r <- cbind(c("Subgroup","","ITT","Per protocol"),
           rbind(r1,r2[-c(1,2),]),
           c("Test for heterogeneity","p-value","",""))
# Add subgroups
# by hiv.status 
s <- logist.comparison.subgroup(base.model=I(disability.week10=="Good")~arm,subgroup.model=~hiv.status,data=d)
r <- rbind(r,s[-c(1:3),])
# output all
colnames(r) <- r[1,]; r <- r[-1,]
rownames(r)<-NULL
kable(r)%>%
  kable_styling("striped", full_width = TRUE)
```


## CLINICAL ADVERSE EVENTS AND NEW LABORATORY ABNORMALITIES
In all adverse event summaries, n.pt.1 and n.pt.2 refer to the number of patients with at least one adverse event in each study arm. n.ae.1 and n.ae.2 refer to the number of adverse events in each study arm. P-values are based on Chi-squared test but descriptive only (due to multiple testing).

```{r, echo=FALSE}
pt.arm <- subset(pop.itt,select=c("usubjid","arm")) 
pt.arm <- plyr::join_all(list(pt.arm,subset(bl.clinical,select=c("usubjid","hiv.status"))))
pt.arm.Hpos <- subset(pt.arm,hiv.status=="Positive",select=c("usubjid","arm")) 
pt.arm.Hneg <- subset(pt.arm,hiv.status=="Negative",select=c("usubjid","arm")) 
```



### Number of adverse events by patient
```{r, echo=FALSE}
num.ae <- ddply(ae,"usubjid",summarise,num.ae=length(ae.name))
tmp <- merge(pt.arm,num.ae,by="usubjid",all.x=T) 
tmp$num.ae[is.na(tmp$num.ae)] <- 0
tmp$num.ae.cat <- cut(tmp$num.ae,
                      breaks=c(-.5,.5,1.5,2.5,5.5,9.5,Inf),
                      labels=c("0","1","2","3-5","6-9","10 or more"))
num.ae.grade34 <- ddply(subset(ae,(ae.grade==""|ae.grade=="Grade 3 or 4")),"usubjid",summarise,num.grade34.ae=length(ae.name))
tmp <- merge(tmp,num.ae.grade34,by="usubjid",all.x=T) 
tmp$num.grade34.ae[is.na(tmp$num.grade34.ae)] <- 0
tmp$num.grade34.ae.cat <- cut(tmp$num.grade34.ae,
                      breaks=c(-.5,.5,1.5,2.5,5.5,9.5,Inf),
                      labels=c("0","1","2","3-5","6-9","10 or more"))
tab <-mySummary.allvar(blvars=tmp[,-c(1:2)],group=tmp$arm,contSummary="med.IQR",pval.comparison=T)
colnames(tab)<-c(tab[2,])
rownames(tab)<-NULL
kable(as.data.frame(tab[-2,]))%>%
  kable_styling("striped", full_width = TRUE)
```

### Adverse events by type
```{r, echo=FALSE}
kable(as.data.frame(mySummary.ae(ae,pt.arm=pt.arm,SOC.only=T))[-2,])%>%
  kable_styling("striped", full_width = TRUE)
```

### Adverse events by type and subtype
```{r, echo=FALSE}
kable(as.data.frame(mySummary.ae(ae,pt.arm=pt.arm,SOC.only=F)))%>%
  kable_styling("striped", full_width = TRUE)
```

### Adverse events by type and subtype - HIV positive only
```{r, echo=FALSE}
kable(as.data.frame(mySummary.ae(ae,pt.arm=pt.arm.Hpos,SOC.only=F))[-2,])%>%
  kable_styling("striped", full_width = TRUE)
```


### Adverse events by type and subtype - HIV negative only
```{r, echo=FALSE}
kable(as.data.frame(mySummary.ae(ae,pt.arm=pt.arm.Hneg,SOC.only=F)))%>%
  kable_styling("striped", full_width = TRUE)
```

### Grade 3 & 4 adverse events by type
```{r, echo=FALSE}
kable(as.data.frame(mySummary.ae(subset(ae,(ae.grade==""|ae.grade=="Grade 3 or 4")),pt.arm=pt.arm,SOC.only=T))[-2,])%>%
  kable_styling("striped", full_width = TRUE)
```

### Grade 3 & 4 adverse events by type and subtype
```{r, echo=FALSE}
kable(as.data.frame(mySummary.ae(subset(ae,(ae.grade==""|ae.grade=="Grade 3 or 4")),pt.arm=pt.arm,SOC.only=F))[-2,])%>%
  kable_styling("striped", full_width = TRUE)
```

### Grade 3 & 4 adverse events by type and subtype - HIV positive only
```{r, echo=FALSE}
kable(as.data.frame(mySummary.ae(subset(ae,(ae.grade==""|ae.grade=="Grade 3 or 4")),pt.arm=pt.arm.Hpos,SOC.only=F))[-2,])%>%
  kable_styling("striped", full_width = TRUE)
```

### Grade 3 & 4 adverse events by type and subtype - HIV negative only
```{r, echo=FALSE}
kable(as.data.frame(mySummary.ae(subset(ae,(ae.grade==""|ae.grade=="Grade 3 or 4")),pt.arm=pt.arm.Hneg,SOC.only=F)))%>%
  kable_styling("striped", full_width = TRUE)
```


### New grade 3 & 4 laboratory abnormalities
```{r, echo=FALSE}
labae$ae.soc <- labae$ae; labae$ae.soc.code <- labae$ae.code
tmp <- mySummary.ae(ae=labae,pt.arm)
colnames(tmp)[1] <- tmp[2,1] <- "Lab abormality name"
tmp[3,1] <- "Any selected lab abormality"
kable(as.data.frame(tmp)[-2,])%>%
  kable_styling("striped", full_width = TRUE)
```


#### Add number of new lab abnormalitites per patient and comparison between arms

```{r, echo=FALSE}
num.ae <- ddply(labae,"usubjid",summarise,num.ae=length(ae))
tmp <- merge(pt.arm,num.ae,by="usubjid",all.x=T) 
tmp$num.ae[is.na(tmp$num.ae)] <- 0
tmp$num.ae.cat <- cut(tmp$num.ae,
                      breaks=c(-.5,.5,1.5,2.5,5.5,9.5,Inf),
                      labels=c("0","1","2","3-5","6-9","10 or more"))
tab <-mySummary.allvar(blvars=tmp[,-c(1:2)],group=tmp$arm,contSummary="med.IQR",pval.comparison=T)
colnames(tab)<-c(tab[2,])
rownames(tab)<-NULL
kable(as.data.frame(tab[-2,]))%>%
  kable_styling("striped", full_width = TRUE)
```



## QT PROLONGATION-SECONDARY ENDPOINT
##### Graphical display (including days 0-28, as for the modeling)


##### Non-linear pattern
```{r, fig.height=15, fig.width=17}
qt<-merge(pop.itt,subset(ecg,(!is.na(qtc))&(day>=0)&(tamofluco.measure.time%in%c("2 hours after","before")),select=c(usubjid,day,qtc,tamofluco.measure.time)),by="usubjid")
```


```{r, include=FALSE}
require(splines)
require(pbkrtest)
qt$day<-qt$day-1# we set the day1 as day 0 in the modeling.
qt<-subset(qt,day>=0)
qt$tamofluco.measure.time<-factor(qt$tamofluco.measure.time,labels=c("before","2 hours after"))
fit_treatment<-lmer(qtc~(ns(day,df=3)+ns(day,df=3):arm)*tamofluco.measure.time+(1+day|usubjid),data=qt)
fit_withouttreatment<-lmer(qtc~(ns(day,df=3))*tamofluco.measure.time+(1+day|usubjid),data=qt)
summary(fit_treatment)$coefficients
a<-KRmodcomp(fit_treatment,fit_withouttreatment)
```

```{r echo=FALSE,fig.height=10, fig.width=12}
fit_treatment<-lmer(qtc~(ns(day,df=3)+ns(day,df=3):arm)*tamofluco.measure.time+(1+day|usubjid),data=qt)
dat<-data.frame(day=rep(0:13,each=4),arm=rep(rep(levels(qt$arm),each=2),by=14),tamofluco.measure.time=rep(rep(levels(qt$tamofluco.measure.time),by=2),by=14))

mm<-model.matrix(~(ns(day,df=3)+ns(day,df=3):arm)*tamofluco.measure.time,dat)
pvar1 <- diag(mm %*% tcrossprod(vcov(fit_treatment),mm))
dat$pred<-mm%*%fixef(fit_treatment) #predict(m,newdat,re.form=NA) would give the same results
dat$plo<-dat$pred-qnorm(0.975,0,1)*sqrt(pvar1)
dat$phi <- dat$pred+qnorm(0.975,0,1)*sqrt(pvar1)
dat$arm <- factor(dat$arm, levels=c("Control","Tamox"))
dat$tamofluco.measure.time <- factor(dat$tamofluco.measure.time, levels=c("before","2 hours after"))  

qt$id<-paste(qt$usubjid,qt$tamofluco.measure.time,sep="-")
ggplot(data=dat,aes(x=day,y=pred))+geom_line(aes(col=tamofluco.measure.time),size=1)+geom_ribbon(aes(x=day,ymin=plo,ymax=phi,fill=tamofluco.measure.time),alpha=0.2)+facet_grid(.~arm)+
  geom_line(data=qt,aes(day,qtc,group =id,color=tamofluco.measure.time),alpha=0.2,size=0.5,position = "dodge")+geom_point(data=qt,aes(day,qtc,color=tamofluco.measure.time),size=1,position="dodge")+scale_y_continuous("QT correction [ms]",breaks=seq(0,600,by=25))+
  scale_x_continuous("Study day",limits = c(0,13),minor_breaks=0:13)+
  theme_bw()+
  theme(panel.background=element_rect(fill = 'white', colour = 'grey'),
        axis.title.x = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b =20, l = 20)),
        axis.title.y = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b = 20, l = 20)),
        axis.text.y = element_text(face="plain", colour="black", size=10),
        axis.text.x = element_text(face="plain", colour="black", size=10),
        axis.ticks.x=element_blank(),
        legend.text = element_text(face="plain",size = 15),
        legend.title = element_blank(),
        legend.key = element_rect(colour = NA),
        legend.position = "right",
        strip.text=element_text(face = "bold",size=15))+scale_color_manual(values=c("#3182bd","#f03b20"))
```



```{r, echo=FALSE}
KRmodcomp(fit_treatment,fit_withouttreatment)
```
 
* Of note, p.value <0.05 says that treatment has effect on QTc.*

### AUC OF QTC OF THE FIRST TWO WEEKS (days 1-14)
```{r, include=FALSE}


AUC_compute<-function(data.auc,main=" "){
  patid<-unique(data.auc$usubjid)
  data<-NULL
  for(i in 1:length(patid)){
    d<-data.auc[data.auc$usubjid==patid[i],]
    auc<-ifelse(nrow(d)>=2,auc.cal(x=d$day,y=d$qtc,min=0,max=16,detectionLimit=log10(4.5)),NA)
    data.c<-data.frame(usubjid=patid[i],auc=auc,arm=d$arm[1])
    data<-rbind(data,data.c)
  }
  
  #data <- merge(bl.lp,data,by="usubjid")
  # Comparion
  tab2 <-mySummary.onevar(varname=main,variable=data$auc,group=data$arm,contSummary="med.IQR",pval.comparison=T)
  # fit<-lm(auc~arm +log10.csf.fungalcount.bl,data=data)
  # tab1<-c("Tamox",round(summary(fit)$coefficients[2,],4))
  # conf<-round(confint(fit),2)
  # # replace crude p-value by adjusted p-value (and estimate, CI)
  # tab2[1,6]<-paste(paste(round(fit$coefficients[2],2),"(",conf[2,1],",",conf[2,2],")",";",sep=""),format.pval(summary(fit)$coefficients[2,"Pr(>|t|)"]))
  tab2<-rbind(c(" ","n","Med(IQR)","n","Med(IQR)","p.value"),tab2)
  colnames(tab2)<-c(" "," ",levels(data$arm)[1]," ",levels(data$arm)[2],"Comparision")
  return(tab2)
}

AUC.all_before <- AUC_compute(subset(qt,tamofluco.measure.time=="before"),main="AUC of QTC before of ITT ")
AUC.all_2_hours <- AUC_compute(subset(qt,tamofluco.measure.time=="2 hours after"),main="AUC of QTC 2 hours after of ITT ")[-1,]
AUC.pp_before <- AUC_compute(subset(qt,pp==1 & tamofluco.measure.time=="before"),main="AUC of QTC before of PP")[-1,]
AUC.pp_2_hours <- AUC_compute(subset(qt,pp==1 & tamofluco.measure.time=="2 hours after"),main="AUC of QTC 2 hours after of PP")[-1,]
QTC.tab<-rbind(AUC.all_before,unname(AUC.all_2_hours),unname(AUC.pp_before),unname(AUC.pp_2_hours))
```

```{r, echo=FALSE}
kable(QTC.tab) %>%
  kable_styling("striped", full_width = TRUE)
```




## RATE OF IRIS
### Rate of IRIS until 10 weeks (main analysis)

```{r, echo=FALSE}
r <- cbind(c("Subgroup","","ITT"),
           surv.comparison(model=Surv(tt.irisdeath.week10,evtype.irisdeath.week10==1)~arm+strata(hiv.status),
                           add.risk=F,add.prop.haz.test=F,data=d),
           c("Test for heterogeneity","p-value",""))
#s <- surv.comparison.subgroup(base.model=Surv(tt.irisdeath.week10,evtype.irisdeath.week10==1)~arm,                              subgroup.model=~hiv.status,data=d,add.risk=F,add.prop.haz.test=F)
#r <- rbind(r,s[-c(1:3),])
colnames(r) <- r[1,]; r <- r[-1,]
kable(r)%>%
  kable_styling("striped", full_width = TRUE)
```
Note: Total number of events
```{r, echo=FALSE}
table(d$evtypetxt.irisdeath.week10,d$arm)
```


## TIME TO NEW AIDS-DEFINING ILLNESS OR DEATH
### Time to new AIDS-defining illness or death until 10 weeks (main analysis)
```{r, echo=FALSE}
r <- cbind(c("Subgroup","","ITT"),
           surv.comparison(model=Surv(tt.aidsdeath.week10,evtype.aidsdeath.week10!=0)~arm+strata(hiv.status),
                           add.risk=T,add.prop.haz.test=F,data=d),
           c("Test for heterogeneity","p-value",""))
#s <- surv.comparison.subgroup(base.model=Surv(tt.aidsdeath.week10,evtype.aidsdeath.week10!=0)~arm, subgroup.model=~hiv.status,data=d,add.risk=T,add.prop.haz.test=F)
#r <- rbind(r,s[-c(1:3),])
colnames(r) <- r[1,]; r <- r[-1,]
kable(r)%>%
  kable_styling("striped", full_width = TRUE)
```


Note: Number of events
```{r, echo=FALSE}
table(d$evtypetxt.aidsdeath.week10,d$arm)
```

## VISUAL DEFICIT AT 10 WEEKS 
#### Desriptive table 
```{r, echo=FALSE}
tab <-mySummary.allvar(blvars=d[,c("visual.acuity.worst.week10","visual.acuity.left.week10","visual.acuity.right.week10","fundoscopy.week10")],group=d$arm,contSummary="med.IQR",pval.comparison=F)
colnames(tab)<-c(tab[2,])
rownames(tab)<-NULL
kable(as.data.frame(tab[-2,])) %>%
  kable_styling("striped", full_width = TRUE)
```



#### Analyses for the probability of a "normal" acuity
Note: Subjects without an acuity assessments or those recorded as "Unable to assess" were excluded from the analysis.

All the data haven't update yet! No result for this endpoint. 

```{r echo=FALSE}
# base model for all analysis except for the analysis by hiv status
d$outcome <- d$visual.acuity.worst.week10
d$outcome[d$outcome=="Unable to assess"] <- NA
base.model <- I(outcome=="Normal for this patient") ~ arm + hiv.status

# ITT 
r <- cbind(c("Subgroup","","ITT","Patients with normal acuity at baseline"),
           rbind(logist.comparison(model=base.model,data=d),
                 logist.comparison(model=base.model,data=subset(d,visual.acuity.worst=="Normal for this patient"))[-c(1:2),]),c("Test for heterogeneity","p-value","",""))

# Add subgroups
# by 
s <- logist.comparison.subgroup(base.model=I(outcome=="Normal for this patient")~arm,subgroup.model=~hiv.status,data=d)
r <- rbind(r,s[-c(1:3),])

# output all
colnames(r) <- r[1,]; r <- r[-1,]
kable(r)%>%
  kable_styling("striped", full_width = TRUE)
```


#### Analyses for the probability of a "normal" fundoscopy result 
Note: only 2 patients have available data.
```{r echo=FALSE}
# base model for all analysis except for the analysis by hiv status
d$outcome <- factor(d$fundoscopy.week10,levels=c("Normal","Abnormal"))
base.model <- I(outcome=="Normal")~arm # do not adjust for continent as no events in Africa
# ITT and asian subjects
r <- cbind(c("","","ITT","HIV"),
           logist.comparison(model=base.model,data=d))
# output all
colnames(r) <- r[1,]; r <- r[-1,]
kable(r)%>%
  kable_styling("striped", full_width = TRUE)
```


## TIME TO NEW NEUROLOGICAL EVENT OR DEATH
### Time to new neurological event or death until 10 weeks (main analysis)
```{r, echo=TRUE}
r <- cbind(c("Subgroup","","ITT"),
           surv.comparison(model=Surv(tt.neurodeath.week10,evtype.neurodeath.week10!=0)~arm+strata(hiv.status),
                           add.risk=T,add.prop.haz.test=F,data=d),
           c("Test for heterogeneity","p-value",""))
s <- surv.comparison.subgroup(base.model=Surv(tt.neurodeath.week10,evtype.neurodeath.week10!=0)~arm,
                              subgroup.model=~hiv.status,data=d,add.risk=T,add.prop.haz.test=F)
r <- rbind(r,s[-c(1:3),])
colnames(r) <- r[1,]; r <- r[-1,]
kable(r)%>%
  kable_styling("striped", full_width = TRUE)
```


Note: Number of events
```{r, echo=TRUE}
table(d$evtypetxt.neurodeath.week10,d$arm)
```


### Rate of new neurological events until 10 weeks (supplementary)
```{r, echo=TRUE}
r <- cbind(c("Subgroup","","ITT"),
           surv.comparison(model=Surv(tt.neurodeath.week10,evtype.neurodeath.week10==1)~arm+strata(hiv.status),
                           add.risk=F,add.prop.haz.test=F,data=d),
           c("Test for heterogeneity","p-value",""))
s <- surv.comparison.subgroup(base.model=Surv(tt.neurodeath.week10,evtype.neurodeath.week10==1)~arm,
                              subgroup.model=~hiv.status,data=d,add.risk=F,add.prop.haz.test=F)
r <- rbind(r,s[-c(1:3),])
colnames(r) <- r[1,]; r <- r[-1,]
kable(r)%>%
  kable_styling("striped", full_width = TRUE)
```




Note: Number of events
```{r, echo=TRUE}
table(d$evtypetxt.neurodeath,d$arm)
```

### Cumulative incidence of new neurological events at 10 weeks (supplementary)

```{r, echo=TRUE}
fit <- survfit(Surv(tt.neurodeath,evtype.neurodeath,type="mstate")~arm,data=d)
tmp <- summary(fit,times=50)
print(data.frame(group=tmp$strata,time=tmp$time,risk.neuroevent=tmp$pstate[,1],risk.death.without.neuroevent=tmp$pstate[,2]),row.names=F)
```


## LONGITUDINAL INTRACRANIAL PRESSURE DURING THE FIRST 2 WEEKS
### Graphical display
```{r, echo=FALSE, fig.height=4, fig.width=8}
op <- merge(pop.itt,subset(lp,(!is.na(openpres))&(day>=0)&(day<=17),select=c(usubjid,day,openpres)),by="usubjid")
op<-merge(subset(bl.clinical,select=c(usubjid,hiv.status)),op,by="usubjid")
op<-merge(subset(bl.lp,select=c(usubjid,log10.csf.fungalcount.bl)),op,by="usubjid")
p <- ggplot(data=op, aes(day,pmax(openpres,1),group =usubjid))

p+geom_line(color="gray")+geom_point(color="darkgray",size=1)+facet_grid(.~arm)+
  scale_y_log10("CSF opening pressure (cm)",breaks=c(2.5,5,10,20,40,80,160))+
  scale_x_continuous("Study day",limits = c(0,17),breaks=c(1,8,15),minor_breaks=0:17)+
  geom_smooth(aes(group=arm),method="loess",se=F,lwd=1.2)+ theme_bw()+
  theme(panel.background=element_rect(fill = 'white', colour = 'grey'),
    axis.title.x = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b =20, l = 20)),
    axis.title.y = element_text(face="bold", colour="black", size=15, margin = margin(t = 20,r = 20,b = 20, l = 20)),
    axis.text.y = element_text(face="plain", colour="black", size=10),
    axis.text.x = element_text(face="plain", colour="black", size=10),
    axis.ticks.x=element_blank(),
    legend.text = element_text(face="plain",size = 15),
    legend.title = element_blank(),
    legend.key = element_rect(colour = NA),
    legend.position = "right",
    strip.text=element_text(face = "bold",size=15))
```



### Naive analysis based on mixed model (ignoring truncation by death) - analysis on log scale
```{r, echo=FALSE}
declineEstimates.percChange.14days <- function(outcome.model,data){
  # calculate percentage change over 14 days assuming in original unit (assuming outcome is modeled on the log10-scale)
  fit.lme1 <- summary(lmer(update(outcome.model,.~studyday+arm1.studyday+(1+studyday|usubjid)),data=data))
  fit.lme2 <- summary(lmer(update(outcome.model,.~studyday+arm2.studyday+(1+studyday|usubjid)),data=data)) 
  n.arm1 <- length(unique(data$usubjid[data$arm==lev.arm[1]])) 
  n.arm2 <- length(unique(data$usubjid[data$arm==lev.arm[2]])) 
  change.arm1 <- paste(formatC(100*10^(14*coef(fit.lme2)["studyday","Estimate"])-100,3,format="f")," (",
                       formatC(100*10^(14*(coef(fit.lme2)["studyday","Estimate"]-qnorm(0.975)*coef(fit.lme2)["studyday","Std. Error"]))-100,3,format="f"),",",
                       formatC(100*10^(14*(coef(fit.lme2)["studyday","Estimate"]+qnorm(0.975)*coef(fit.lme2)["studyday","Std. Error"]))-100,3,format="f"),")",sep="")
  change.arm2 <- paste(formatC(100*10^(14*coef(fit.lme1)["studyday","Estimate"])-100,3,format="f")," (",
                       formatC(100*10^(14*(coef(fit.lme1)["studyday","Estimate"]-qnorm(0.975)*coef(fit.lme1)["studyday","Std. Error"]))-100,3,format="f"),",",
                       formatC(100*10^(14*(coef(fit.lme1)["studyday","Estimate"]+qnorm(0.975)*coef(fit.lme1)["studyday","Std. Error"]))-100,3,format="f"),")",sep="")
  diff.1vs2 <-   paste(formatC(100*10^(14*coef(fit.lme2)["arm2.studyday","Estimate"])-100,3,format="f")," (",
                       formatC(100*10^(14*(coef(fit.lme2)["arm2.studyday","Estimate"]-qnorm(0.975)*coef(fit.lme2)["arm2.studyday","Std. Error"]))-100,3,format="f"),",",
                       formatC(100*10^(14*(coef(fit.lme2)["arm2.studyday","Estimate"]+qnorm(0.975)*coef(fit.lme2)["arm2.studyday","Std. Error"]))-100,3,format="f"),"); p value ",
                       format.pval(2*pnorm(-abs(coef(fit.lme1)["arm1.studyday","t value"]))),sep="")
  list(n.arm1=n.arm1,change.arm1=change.arm1,n.arm2=n.arm2,change.arm2=change.arm2,diff.1vs2=diff.1vs2)
}

```



```{r, echo=FALSE}
# covariates for treatment arm - day interactions (which is only in effect after day 1)
op$studyday <- pmax(op$day-1,0)
op$arm1.studyday <- op$studyday*(op$arm==lev.arm[1])
op$arm2.studyday <- op$studyday*(op$arm==lev.arm[2])
```




```{r, echo=FALSE}
# Estimate declines
declines.all <- declineEstimates.percChange.14days(log10(pmax(openpres,1))~.,op)
declines.pp <- declineEstimates.percChange.14days(log10(pmax(openpres,1))~.,subset(op,pp==1))
declines.hiv.pos <- declineEstimates.percChange.14days(log10(pmax(openpres,1))~.,subset(op,hiv.status=="Positive"))
declines.hiv.neg <- declineEstimates.percChange.14days(log10(pmax(openpres,1))~.,subset(op,hiv.status=="Negative"))
declines.low <- declineEstimates.percChange.14days(log10(pmax(openpres,1))~.,subset(op,log10.csf.fungalcount.bl<5))
declines.mid <- declineEstimates.percChange.14days(log10(pmax(openpres,1))~.,subset(op,(log10.csf.fungalcount.bl>=5)&(log10.csf.fungalcount.bl<=6)))
declines.high <- declineEstimates.percChange.14days(log10(pmax(openpres,1))~.,subset(op,log10.csf.fungalcount.bl>6))


op.tab<-matrix("",nrow=9,ncol=6)
op.tab[1,]<-c("","",lev.arm[1],"",lev.arm[2],"Comparison")
op.tab[2,]<-c("Continent","n","% change over 14 days (95% CI)","n","% change over 14 days (95% CI)","Rel. difference in % change over 14 days")
op.tab[3,]<-c("All",unlist(declines.all))
op.tab[4,]<-c("PP",unlist(declines.pp))
op.tab[5,]<-c("HIV positive patients only",unlist(declines.hiv.pos))
op.tab[6,]<-c("HIV negative patients only",unlist(declines.hiv.neg))
op.tab[7,]<-c("Low fungal burden ",unlist(declines.low))
op.tab[8,]<-c("Middle fungal burden ",unlist(declines.mid))
op.tab[9,]<-c("High fungal burden ",unlist(declines.high))
colnames(op.tab) <- op.tab[1,]; op.tab <- op.tab[-1,]
kable(op.tab)%>%
  kable_styling("striped", full_width = TRUE)

```


<!-- ##   RELAPSE IN THE 10 WEEKS AFTER RANDOMIZATION -->
<!-- There were no relapse record so the analysis for relapse will be obmited. -->
<!-- ```{r eval=FALSE, include=FALSE} -->
<!-- r <- cbind(c("Subgroup","","ITT"), -->
<!--            surv.comparison(model=Surv(tt.relapse,evtype.relapse==1)~arm+strata(hiv.status), -->
<!--                            add.risk=F,add.prop.haz.test=F,data=d), -->
<!--            c("Test for heterogeneity","p-value","")) -->
<!-- #s <- surv.comparison.subgroup(base.model=Surv(tt.relapse,evtype.relapse==1)~arm,                              subgroup.model=~hiv.status,data=d,add.risk=F,add.prop.haz.test=F) -->
<!-- #r <- rbind(r,s[-c(1:3),]) -->
<!-- colnames(r) <- r[1,]; r <- r[-1,] -->
<!-- kable(r)%>% -->
<!--   kable_styling("striped", full_width = TRUE) -->
<!-- ``` -->


<!-- Note: Total number of events -->
<!-- ```{r eval=FALSE, include=FALSE} -->
<!-- table(d$tt.relapse,d$evtype.relapse) -->
<!-- ``` -->


<!-- ### Cumulative incidence of relapses at 10 weeks and 6 months (supplementary) -->
<!-- ```{r, echo=FALSE} -->
<!-- fit <- survfit(Surv(tt.relapse,evtype.relapse,type="mstate")~arm,data=d) -->
<!-- tmp <- summary(fit,times=70) -->
<!-- print(data.frame(group=tmp$strata,time=tmp$time,risk.relapse=tmp$pstate[,1],risk.death.without.relapse=tmp$pstate[,2]),row.names=F) -->
<!-- ``` -->

## ADDITIONAL PLANNED AUXILIARY ANALYSES
### ARV status, time to ARV after enrolment, duration of AmB treatment (post enrolment)
```{r, echo=FALSE}
tmp <- merge(pop.itt[,-4],arv,by="usubjid",all.x=T)
tmp <- merge(tmp,subset(population,select=c(usubjid,days.amphob)),by="usubjid")
tmp$days.amphob.cat <- cut(tmp$days.amphob,breaks=c(-0.5,0.5,5.5,9.5,13.5,14.5,Inf),
                           labels=c("0 days","1-5 days","6-9 days","10-13 days","14 days","More than 14 days"))
tab <-mySummary.allvar(blvars=tmp[,-c(1:6)],group=tmp$arm,contSummary="med.IQR",pval.comparison=F)
colnames(tab)<-c(tab[2,])
rownames(tab)<-NULL
kable(as.data.frame(tab[-2,]))%>%
  kable_styling("striped", full_width = TRUE)
```


### Number of X-rays, CT, MRI performed
```{r, echo=FALSE}
# Xrays
tmp <- merge(pop.itt,xray,by="usubjid")
num.xray <- ddply(tmp,.(arm),summarise,
  num.Xray=length(xray.result),
  num.Xray.abnormal=sum(xray.result=="Abnormal", na.rm = T),
  num.Xray.after.bl=sum(day>3, na.rm = T),
  num.Xray.abnormal.after.bl=sum((day>3)&(xray.result=="Abnormal"), na.rm = T))
# CT
tmp <- merge(pop.itt,subset(brainscan,scan.method=="CT"),by="usubjid")
num.ct <- ddply(tmp,.(arm),summarise,
                num.CT=length(scan.result),
                num.CT.abnormal=sum(scan.result=="Abnormal", na.rm = T),
                num.CT.after.bl=sum(day>3, na.rm = T),
                num.CT.abnormal.after.bl=sum((day>3)&(scan.result=="Abnormal"), na.rm = T))
# MRI
tmp <- merge(pop.itt,subset(brainscan,scan.method=="MRI"),by="usubjid")
num.mri <- ddply(tmp,.(arm),summarise,
                 num.MRI=length(scan.result),
                 num.MRI.abnormal=sum(scan.result=="Abnormal", na.rm = T),
                 num.MRI.after.bl=sum(day>3, na.rm = T),
                 num.MRI.abnormal.after.bl=sum((day>3)&(scan.result=="Abnormal"), na.rm = T))

#+ output=TRUE
print(num.xray,row.names=F)
print(num.ct,row.names=F)
print(num.mri,row.names=F)
```


### CD4 comparision
```{r}
library(plyr)
data.cd4<-cd4and8
data.cd4<-merge(data.cd4,subset(bl.lab,select=c(usubjid,cd4.bl)))
dat<-ddply(data.cd4,.(usubjid),summarize,day.first=min(day),day.last=max(day))
data.cd4<-merge(data.cd4,dat,by="usubjid")
data.cd4<-subset(data.cd4,day==day.last)
data.cd4$change.cd4<-(data.cd4$cd4)-(data.cd4$cd4.bl)
data.cd4<-merge(data.cd4,subset(pop.itt,select=c(usubjid,arm,pp)),by="usubjid")
data.cd4<-merge(data.cd4,bl.clinical,by="usubjid")
data.cd4<-merge(data.cd4,subset(bl.lp,select=c(usubjid,log10.csf.fungalcount.bl)),by="usubjid")

cd4.change<-function(data,main=" "){
# Comparion
  
tab2 <-mySummary.onevar(varname=main,variable=data$change.cd4,group=data$arm,contSummary="med.IQR",pval.comparison=T)
# fit<-lm(change~arm +log10.lacto.load.bl,data=data.day7)
# tab1<-c("armdummy.B",round(summary(fit)$coefficients[2,],4))
# conf<-round(confint(fit),2)
# # replace crude p-value by adjusted p-value (and estimate, CI)
#tab2[1,6]<-paste(paste(round(fit$coefficients[2],2),"(",conf[2,1],",",conf[2,2],")",";",sep=""),format.pval(summary(fit)$coefficients[2,"Pr(>|t|)"]))
tab2<-rbind(c(" ","n","Med(IQR)","n","Med(IQR)","p.value"),tab2)
colnames(tab2)<-c(" "," ",levels(data$arm)[1]," ",levels(data$arm)[2],"Comparision")
return(tab2)
}

CD4.change.all <- cd4.change(data.cd4,main="CD4 change after treatment ITT")
CD4.change.pp <- cd4.change(subset(data.cd4,pp==1),main="CD4 change after treatment PP")[-1,]
CD4.change.hiv.pos <- cd4.change(subset(data.cd4,hiv.status=="Positive"),main="CD4 change after treatment HIV +")[-1,]
CD4.change.hiv.neg <- cd4.change(subset(data.cd4,hiv.status=="Negative"),main="CD4 change after treatment HIV -")[-1,]

CD4.change.low <- cd4.change(subset(data.cd4,log10.csf.fungalcount.bl<5),main="CD4 change after treatment low burden")[-1,]
CD4.change.middle.n.high <- cd4.change(subset(data.cd4,(log10.csf.fungalcount.bl>=5)&(log10.csf.fungalcount.bl<=6)),main="CD4 change after treatment middle & high burden")[-1,]
CD4.tab<-rbind(CD4.change.all,unname(CD4.change.pp),unname(CD4.change.hiv.pos),unname(CD4.change.hiv.neg),unname(CD4.change.low),unname(CD4.change.middle.n.high))
```

```{r, echo=FALSE}
kable(CD4.tab)%>%
  kable_styling("striped", full_width = TRUE)
```


### AUC OF CSF FUNGAL COUNT ACTIVITY OF THE FIRST TWO WEEKS (days 1-17)
```{r, include=FALSE}


AUC_compute<-function(data.auc,main=" "){
  patid<-unique(data.auc$usubjid)
  data<-NULL
  for(i in 1:length(patid)){
  d<-data.auc[data.auc$usubjid==patid[i],]
  auc<-ifelse(nrow(d)>=2,auc.cal(x=d$day,y=d$log10.csf.fungalcount,min=0,max=16,detectionLimit=log10(4.5)),NA)
  data.c<-data.frame(usubjid=patid[i],auc=auc,arm=d$arm[1])
  data<-rbind(data,data.c)
}

#data <- merge(bl.lp,data,by="usubjid")
# Comparion
tab2 <-mySummary.onevar(varname=main,variable=data$auc,group=data$arm,contSummary="med.IQR",pval.comparison=T)
# fit<-lm(auc~arm +log10.csf.fungalcount.bl,data=data)
# tab1<-c("Tamox",round(summary(fit)$coefficients[2,],4))
# conf<-round(confint(fit),2)
# # replace crude p-value by adjusted p-value (and estimate, CI)
# tab2[1,6]<-paste(paste(round(fit$coefficients[2],2),"(",conf[2,1],",",conf[2,2],")",";",sep=""),format.pval(summary(fit)$coefficients[2,"Pr(>|t|)"]))
tab2<-rbind(c(" ","n","Med(IQR)","n","Med(IQR)","p.value"),tab2)
colnames(tab2)<-c(" "," ",levels(data$arm)[1]," ",levels(data$arm)[2],"Comparision")
return(tab2)
}

AUC.all <- AUC_compute(fc,main="AUC of CSF fungal of ITT")
AUC.pp <- AUC_compute(subset(fc,pp==1),main="AUC of CSF fungal of PP")[-1,]
AUC.hiv.pos <- AUC_compute(subset(fc,hiv.status=="Positive"),main="AUC of CSF fungal of HIV +")[-1,]
AUC.hiv.neg <- AUC_compute(subset(fc,hiv.status=="Negative"),main="AUC of CSF fungal of HIV -")[-1,]

AUC.low <- AUC_compute(subset(fc,log10.csf.fungalcount.bl<5),main="AUC of CSF fungal of low burden")[-1,]
AUC.mid <- AUC_compute(subset(fc,(log10.csf.fungalcount.bl>=5)&(log10.csf.fungalcount.bl<=6)),main="AUC of CSF fungal of middle burden")[-1,]
AUC.high <- AUC_compute(subset(fc,log10.csf.fungalcount.bl>6),main="AUC of CSF fungal of high burden")[-1,]
fungal.tab<-rbind(AUC.all,unname(AUC.pp),unname(AUC.hiv.pos),unname(AUC.hiv.neg),unname(AUC.low),unname(AUC.mid),unname(AUC.high))
```

```{r, echo=FALSE}
kable(fungal.tab) %>%
  kable_styling("striped", full_width = TRUE)
```


```{r eval=FALSE, include=FALSE, results='asis'}
cat('<style>.toc-content{z-index: 1000} td{font-size:10.2px; padding:5.2px!important} th{font-size:11px}</style>')
```