#Author: Marcel Wolbers
#Last modification: 14Mar2015

## Summary for baseline characteristics

mySummary.allvar <- function(blvars,group,pooledGroup=F,contSummary="med.IQR",caption=NULL,filename=NA,pval.comparison=F,cont=NA){
  # contSummary can be median (90% range) "med.90" or median (IQR) "med.IQR" or median (range) "med.range" or "MIC"
  require(xtable); require(Hmisc); library(gdata)
  if (is.factor(group)) {
    levelOrder <- levels(group)
    was.factor <- T
  } else {was.factor <- F}
  group <- as.character(group)
  if (was.factor) {group <- factor(group,levels=levelOrder[!is.na(match(levelOrder,group))])
  }else group <- factor(group)
  gr.lev <- levels(group)
  if (pooledGroup){ # Add pooled summaries for all patients
    mylabels <- unlist(lapply(blvars,Hmisc::label)) # save them as rbind destroys some of them
    blvars <- rbind(blvars,blvars)
    for (i in 1:ncol(blvars)) label(blvars[,i]) <- mylabels[i] # add labels again
    group <- c(as.character(group),rep("All patients",length(group)))
    group <- factor(group,levels=c("All patients",gr.lev))
    gr.lev <- levels(group)
  }
  header1 <- c("",c(rbind(rep("",length(gr.lev)),paste(gr.lev," (N=",table(group),")",sep=""))))
  header2 <- c("Characteristic",rep(c("n","Summary statistic"),length(gr.lev)))
  result <-  rbind(header1,header2)
  if (pval.comparison) result <- cbind(result,c("Comparison","(p-value)"))
  if (is.na(cont)[1]) cont <- rep(cont,ncol(blvars))
  for (i in 1:ncol(blvars)){
    result.i <- mySummary.onevar(varname=ifelse(Hmisc::label(blvars[,i])!="",label(blvars[,i]),names(blvars)[i]),
                                 blvars[,i],group,contSummary=contSummary,pval.comparison=pval.comparison,cont=cont[i])
    result <- rbind(result,result.i)
  }
  rownames(result) <- rep("",nrow(result))
  if (!is.na(filename)){  # generate html table
    x <- print(xtable(result,caption=caption),type="html",file=filename,
               caption.placement="top",include.rownames=F,include.colnames=F,
               html.table.attributes="border=1")
    cont.summary.options <- c("med.IQR","median.90","median.range","MIC")      
    cont.summary.text <- c("median (IQR)","median (90% range)","median (range)","median, 90% quantile and range")      
    x <- paste(x,"Summary statistic is absolute count (%) for categorical variables and ",
               cont.summary.text[match(contSummary,cont.summary.options)],"for continuous data.",
               ifelse(pval.comparison,"Comparisons based on Fisher's exact test for categorical data and Wilcoxon test for continuous data.",""))
    write(x,file=filename)
  }
  result
}

mySummary.onevar <- function(varname,variable,group,cont=NA,contSummary="med.90",pval.comparison=F){
  require(gdata)
  if (is.na(cont)) cont <- ifelse(is.factor(variable)|length(unique(na.omit(variable)))<=5,F,T)
  ngroup <- length(levels(group))
  mycont.summary <- function(variable,group) {
    if (is.na(match(contSummary,c("med.90","med.IQR","med.range","MIC"))))
      stop("contSummary=",contSummary," not yet implemented")
    n <- c(by(variable,group,function(x) length(na.omit(x))))
    if (contSummary=="med.90") summarystat <- by(variable,group,function(x) quantile(x,c(0.05,0.5,0.95),na.rm=T))
    if (contSummary=="med.IQR") summarystat <- by(variable,group,function(x) quantile(x,c(0.25,0.5,0.75),na.rm=T))
    if (contSummary=="med.range") summarystat <- by(variable,group,function(x) quantile(x,c(0,0.5,1),na.rm=T))
    if (!is.na(match(contSummary,c("med.90","med.IQR","med.range")))){
      summarystat.nice <- lapply(summarystat,function(x){ x <- formatC(round(x,2),2,format="f");
                                                          paste(x[2],"(",x[1],",",x[3],")",sep="")})
      result <- matrix("",ncol=ngroup*2+1,nrow=1)
      result[1,seq(2,ncol(result),by=2)] <- n
      result[1,seq(3,ncol(result),by=2)] <- unlist(summarystat.nice)
    }
    if (contSummary=="MIC") {
      result <- matrix("",ncol=ngroup*2+1,nrow=4)
      result[,1] <- c("","- MIC 50","- MIC 90","- range")
      result[1,seq(2,ncol(result),by=2)] <- n
      medians <- by(variable,group,function(x) quantile(x,c(0.5),na.rm=T))
      q.90s <- by(variable,group,function(x) quantile(x,c(0.9),na.rm=T))
      ranges <- by(variable,group,function(x) quantile(x,c(0,1),na.rm=T))
      result[2,seq(3,ncol(result),by=2)] <- unlist(lapply(medians,function(x){x <- formatC(round(x,2),2,format="f")}))      
      result[3,seq(3,ncol(result),by=2)] <- unlist(lapply(q.90s,function(x){x <-   formatC(round(x,2),2,format="f")}))
      result[4,seq(3,ncol(result),by=2)] <- unlist(lapply(ranges,function(x) {x <- formatC(round(x,2),2,format="f"); 
                                                                              paste("(",x[1],",",x[2],")",sep="")}))
      # # Replace "257" by ">256" (and "33" by ">32" for cotrimoxazole) 
      # result <- gsub("257",">256",result)
      # if (varname=="cotrimoxazole") result <- gsub("33",">32",result)
    }     
    if (pval.comparison) {
      pval <- kruskal.test(variable[group!="All patients"]~group[group!="All patients"])$p.value
      if (pval<0.001) pval <- "<0.001" else pval <- formatC(round(pval,3))
      result <- cbind(result,"")
      result[1,ncol(result)] <- pval
    }
    result
  }
  mycat.summary <- function(variable,group) {
    ta <- table(group,variable)
    ta.n <- apply(ta,1,sum)
    ta.prop <- ta/ta.n
    ta.nice <- matrix(paste(ta,"/",ta.n," (",round(100*ta.prop),"%",")",sep=""),nrow=nrow(ta),ncol=ncol(ta))
    result <- matrix("",ncol=ngroup*2+1,nrow=ncol(ta)+1)
    result[2:nrow(result),1] <- paste("- ",colnames(ta))
    result[2:nrow(result),seq(3,ncol(result),by=2)] <- t(ta.nice)
    result[1,seq(2,ncol(result),by=2)] <- apply(ta,1,sum) # n's
    if (pval.comparison) {
      if (ncol(table(group[group!="All patients"],variable[group!="All patients"]))==1) {pval <- "NA"
      }else {
        # Use simulated p-value if normal fisher-test doesn't work
        options(show.error.messages=F) 
        tab<- table(group[group!="All patients"],variable[group!="All patients"])
        if(prod(as.numeric(chisq.test(tab)$expected>1))==0){
          ft <- fisher.test(table(group[group!="All patients"],variable[group!="All patients"]))
          options(show.error.messages=T)
          if (class(ft)!="try-error"){
            pval <- ft$p.value
          } else {
            warning("Simulated p-values for Fisher test used with B=10000")
            pval <- fisher.test(table(group[group!="All patients"],variable[group!="All patients"]),
                                simulate.p.value=T,B=10000)$p.value
          }
        }else{
          ft <- chisq.test(table(group[group!="All patients"],variable[group!="All patients"]))
          pval <- ft$p.value
        }

        if (pval<0.001) pval <- "<0.001" else pval <- formatC(round(pval,3))
      }
      result <- cbind(result,c(pval,rep("",nrow(result)-1)))
    }
    result
  }
  if (cont) {r <- mycont.summary(variable,group)
  }else {r <- mycat.summary(variable,group)}
  r[1,1] <- varname
  r
}


## Summary for AE
mySummary.crude.ae <- function(ae,pt.arm,by.SOC=T){
  # create a crude AE table either by system organ class only (by.SOC=T) or by ae preferred term only (by.SOC=F)
  # function assumes the following structure
  # ae: columns usubjid, ae.soc and ae.soc.code (body system, if by.SOC=T), ae.pt and ae.pt.code (ae term, if by.SOC=F)
  # pt.arm: columns usubjid (unique patient identifier) and arm (randomized treatment group)

  if (nrow(ae)==0){stop("There are no adverse event to summarize")}

  # chose correct AE name and code
  if (by.SOC==T){
    ae$aename<-ae$ae.soc
    ae$aecode<-ae$ae.soc.code
  }
  if (by.SOC==F){
    ae$aename<-ae$ae.pt
    ae$aecode<-paste(ae$ae.soc.code,ae$ae.pt.code,sep="")
  }
  
  # add dummy rows to generate numbers of AEs of any type
  ae.any <- ae; ae.any$aecode <- "AAE"; ae.any$aename <- "All AEs combined"
  ae <- rbind(ae,ae.any)
  # add randomized arm to AE datasets and keep only subjects for which an arm was provided (useful for subgroup analyses)
  ae.arm <- merge(pt.arm,ae,by="usubjid")
  if(by.SOC==T){
    ae.arm$aecode<-factor(ae.arm$aecode,levels=sort(unique(as.character(ae.arm$aecode)))) 
  }

  if (nrow(ae.arm)==0){stop("There are no adverse event to summarize")}

  gr.lev <- levels(pt.arm$arm)
  n.1 <- sum(pt.arm$arm==gr.lev[1])
  n.2 <- sum(pt.arm$arm==gr.lev[2])

  # number of AE of each type
  all.aecodes <- sort(unique(ae.arm$aecode))#factor(sort(unique(as.character(ae.arm$aecode))))
  all.aenames.indexes <- match(all.aecodes,ae.arm$aecode)
  all.aenames <- ae.arm$aename[all.aenames.indexes]
  if (any(duplicated(all.aenames))) warning("Identical AE names with differnt AE codes: ",all.aenames[duplicated(all.aenames)])

  ae.count <- data.frame("ae.code"=all.aecodes,"ae.name"=all.aenames,"n.pt.1"=NA,"n.ae.1"=NA,"n.pt.2"=NA,"n.ae.2"=NA,"p.val"=NA,stringsAsFactors=F)
  for (i in 1:nrow(ae.count)){
    pt.ae.1 <- ae.arm$usubjid[((ae.arm$aecode==ae.count$ae.code[i])&(ae.arm$arm==gr.lev[1]))&(!is.na(ae.arm$aecode))]
    ae.count$n.ae.1[i] <- length(pt.ae.1)
    ae.count$n.pt.1[i] <- length(unique(pt.ae.1))
    pt.ae.2 <- ae.arm$usubjid[((ae.arm$aecode==ae.count$ae.code[i])&(ae.arm$arm==gr.lev[2]))&(!is.na(ae.arm$aecode))]
    ae.count$n.ae.2[i] <- length(pt.ae.2)
    ae.count$n.pt.2[i] <- length(unique(pt.ae.2))
    tab<-cbind(c(ae.count$n.pt.1[i],n.1-ae.count$n.pt.1[i]),c(ae.count$n.pt.2[i],n.2-ae.count$n.pt.2[i]))
    if(prod(as.numeric(chisq.test(tab)$expected>1))==0){
      ae.count$p.val[i] <- fisher.test(tab)$p.value
      
    }else{
      ae.count$p.val[i] <- chisq.test(tab)$p.value
    }
    
  }
  ae.count
}


mySummary.ae <- function(ae,pt.arm,SOC.only=T){
  # create AE summary either by body system only (SOC.only=T) or by body system and preferred term (SOC.only=F)
  # see mySummary.crude.ae for the assumed structure of arguments
  gr.lev <- levels(pt.arm$arm)
  n.1 <- sum(pt.arm$arm==gr.lev[1])
  n.2 <- sum(pt.arm$arm==gr.lev[2])
  if (SOC.only==T){
    ae.count <- mySummary.crude.ae(ae,pt.arm,by.SOC=T)
  }
  if (SOC.only==F){
    ae.count.SOC <- mySummary.crude.ae(ae,pt.arm,by.SOC=T)
    ae.count.PT <- mySummary.crude.ae(ae,pt.arm,by.SOC=F)
    ae.count.PT <- ae.count.PT[-2,] # drop overall count
    ae.count.PT$ae.name <- paste(" -",ae.count.PT$ae.name,sep=" ")
    ae.count <- rbind(ae.count.SOC,ae.count.PT)
    ae.count$ae.code<-as.character(ae.count$ae.code)
    ae.count<-arrange(ae.count,ae.code)
    ae.count$ae.code<-as.factor(ae.count$ae.code)
  }
  ae.count.final <- cbind("Adverse event name"=ae.count$ae.name,
                          "n.pt.1"=paste(ae.count$n.pt.1," (",round(100*ae.count$n.pt.1/n.1,2),"%)",sep=""),
                          "n.ae.1"=as.character(ae.count$n.ae.1),
                          "n.pt.2"=paste(ae.count$n.pt.2," (",round(100*ae.count$n.pt.2/n.2,2),"%)",sep=""),
                          "n.ae.2"=as.character(ae.count$n.ae.2),
                          "(p value)"=as.character(round(ae.count$p.val,3)))
  ae.count.final <- rbind(c("",paste(gr.lev[1]," (n=",n.1,")",sep=""),"",paste(gr.lev[2]," (n=",n.2,")",sep=""),"","Comparison"),
                          colnames(ae.count.final),ae.count.final)
  rownames(ae.count.final)<-NULL
  ae.count.final
}


## KM curves
## KM curves
mykmplot <- function(survmodel,data,main="",xlab="",ylab="Time since randomization [days]",x.max=71){
  tp <- seq(0,x.max,by=14)
  lev <- levels(data$arm)
  lev.labels <- lev
  attr(survmodel,".Environment") <- environment()
  survfit.group1 <- survfit(survmodel,data=data,subset=(arm==lev[1]))
  at.risk.group1 <- summary(survfit.group1,times=tp,extend=T)$n.risk
  survfit.group2 <- survfit(survmodel,data=data,subset=(arm==lev[2]))
  at.risk.group2 <- summary(survfit.group2,times=tp,extend=T)$n.risk
  plot(survfit.group1,main=main,mark.time=F,conf.int=F,
       ylim=c(-0.5,1),xlim=c(-15,x.max),
       xlab="",ylab="",col="black",lty=1,axes=F,lwd=2)
  lines(survfit.group2,mark.time=F,conf.int=F,col=gray(0.7),lty="41",lwd=2)
  axis(1,at=seq(0,x.max,by=7),labels=NA,pos=0,adj=1)
  axis(1,at=tp,pos=0)
  axis(2,at=c(0,0.25,0.5,0.75,1),pos=0)
  mtext(text=xlab,side=2,at=0.5,line=2,cex=1.2,padj=4.5)
  mtext(text=ylab,side=1,at=round(x.max/2), line=2,cex=1.2,padj = -15)
  mtext(text="No. at risk",side=1,at=-round(0.25*x.max), line=3,cex=1,padj=-14.8)
  mtext(text=lev.labels[1],side=1,at=-round(0.25*x.max),line=4,cex=1,adj=0.25,padj=-14.3)
  mtext(text=lev.labels[2],side=1,at=-round(0.25*x.max),line=5,cex=1,adj=0.25,padj=-13.8)
  mtext(text=at.risk.group1,sid=1,at=tp, line=4,cex=1,adj=0.5,padj=-14.3)
  mtext(text=at.risk.group2,sid=1,at=tp, line=5,cex=1,adj=0.5,padj=-13.8)
  legend(x=round(0.6*x.max),y=0.25,legend=c(lev.labels[1],lev.labels[2]),lty=c("solid","41"),lwd=2,col=c("black",gray(0.7)),cex=1.2)
  invisible(NULL)
}



## Treatment comparisons overall and in subgroups - Cox regression
format.pval <- function(p){
  ifelse(p<0.0001,"<0.0001",formatC(p,4,format="f"))
}

surv.comparison.death <- function(model,data,add.risk=T,add.prop.haz.test=T){
  ##---------------------------------------------------------------------------------------------------------
  ## Purpose: Summarize results for a Cox survival model with the treatment arm (variable "arm") as the main covariate
  ## !! model formula can include other covariates than arm BUT arm must be the first covariate in the model !!
  ## add.risk: if T, the event probability ("absolute risk") at time "infinity" is also displayed
  ## add.prop.haz.test: if T, a test for proportional hazards is also added
  ## Author: Marcel Wolbers, 9March2015
  ##---------------------------------------------------------------------------------------------------------
  arm.names <- levels(data$arm)
  # Table header
  header1 <- c(paste(arm.names,"(n=",table(data$arm),")",sep=""),"Comparison")
  header2 <- c(rep(ifelse(add.risk,"Deaths/n (risk [%])","Deaths/n"),2),"HR (95%CI);p-value")
  header <- rbind(header1,header2)
  result <- rbind(header,"")
  # add number of events and risks
  fit.surv <- summary(survfit(update(model,.~arm),data=data),times=1e10,extend=T)
  events.n <- paste(fit.surv$table[,"events"],fit.surv$table[,"n.max"],sep="/")
  if (add.risk) events.n <- paste(events.n," (",formatC(100*(1-fit.surv$surv),2,format="f"),")",sep="")
  result[3,1:2] <- events.n
  # add HR, CI, p-value
  fit.coxph <- summary(coxph(model,data))
  hr <- formatC(fit.coxph$coef[1,"exp(coef)"],2,format="f")
  pval <- format.pval(fit.coxph$coef[1,"Pr(>|z|)"])
  ci <- paste(formatC(fit.coxph$conf.int[1,c("lower .95")],2,format="f"),
              formatC(fit.coxph$conf.int[1,c("upper .95")],2,format="f"),sep="-")
  hr.ci.p <- paste(hr,"(",ci,"); p=",pval,sep="")
  result[3,3] <- hr.ci.p
  # add test for proportional hazards
  if (add.prop.haz.test){
    attr(model,".Environment") <- environment() # needed for cox.zph to work 
    p.prop.haz <- cox.zph(coxph(model,data))$table[1,"p"]
    result <- cbind(result,c("Test for proportional hazards","p-value",format.pval(p.prop.haz)))  
  }
  rownames(result) <- NULL
  result
}


surv.comparison <- function(model,data,add.risk=T,add.prop.haz.test=T){
  ##---------------------------------------------------------------------------------------------------------
  ## Purpose: Summarize results for a Cox survival model with the treatment arm (variable "arm") as the main covariate
  ## !! model formula can include other covariates than arm BUT arm must be the first covariate in the model !!
  ## add.risk: if T, the event probability ("absolute risk") at time "infinity" is also displayed
  ## add.prop.haz.test: if T, a test for proportional hazards is also added
  ## Author: Marcel Wolbers, 9March2015
  ##---------------------------------------------------------------------------------------------------------
  arm.names <- levels(data$arm)
  # Table header
  header1 <- c(paste(arm.names,"(n=",table(data$arm),")",sep=""),"Comparison")
  header2 <- c(rep(ifelse(add.risk,"events/n (risk [%])","events/n"),2),"HR (95%CI);p-value")
  header <- rbind(header1,header2)
  result <- rbind(header,"")
  # add number of events and risks
  fit.surv <- summary(survfit(update(model,.~arm),data=data),time=1e10,extend=T)
  events.n <- paste(fit.surv$table[,"events"],fit.surv$table[,"n.max"],sep="/")
  if (add.risk) events.n <- paste(events.n," (",formatC(100*(1-fit.surv$surv),2,format="f"),")",sep="")
  result[3,1:2] <- events.n
  # add HR, CI, p-value
  fit.coxph <- summary(coxph(model,data))
  hr <- formatC(fit.coxph$coef[1,"exp(coef)"],2,format="f")
  pval <- format.pval(fit.coxph$coef[1,"Pr(>|z|)"])
  ci <- paste(formatC(fit.coxph$conf.int[1,c("lower .95")],2,format="f"),
              formatC(fit.coxph$conf.int[1,c("upper .95")],2,format="f"),sep="-")
  hr.ci.p <- paste(hr,"(",ci,"); p=",pval,sep="")
  result[3,3] <- hr.ci.p
  # add test for proportional hazards
  if (add.prop.haz.test){
    attr(model,".Environment") <- environment() # needed for cox.zph to work 
    p.prop.haz <- cox.zph(coxph(model,data))$table[1,"p"]
    result <- cbind(result,c("Test for proportional hazards","p-value",format.pval(p.prop.haz)))  
  }
  rownames(result) <- NULL
  result
}

surv.comparison.subgroup <- function(base.model,subgroup.model,data,...){
  ##---------------------------------------------------------------------------------------------------------
  ## Purpose: Summarize results for a Cox survival model by treatment arm (variable "arm") and subgroup
  ## base.model: Model from which sub-group specific estimates are extracted (!! arm must be the first covariate in the model)
  ## subgroup.model: model of the form "~subgrouping.variable1+subgrouping.variable2" 
  ##                 (!! subgrouping.variable must be factors and there should be nothing on the left-hand side of the formula)
  ## ...: arguments that are passed to surv.comparison
  ## Author: Marcel Wolbers, 20March2015
  ##---------------------------------------------------------------------------------------------------------
  # result in entire population
  result <- surv.comparison(base.model,data,...)
  result <- cbind(c("Subgroup","","All patients"),result,c("Test for heterogeneity","p-value",""))
  # Preparation of models and data
  subgroup.char <- attr(terms(subgroup.model),"term.labels")
  for (k in 1:length(subgroup.char)){
    main.model <- update(base.model,as.formula(paste(".~.+",subgroup.char[k],sep=""))) 
    ia.model <- update(base.model,as.formula(paste(".~.+arm*",subgroup.char[k],sep=""))) 
    data$.subgroup.var <- data[,subgroup.char[k]]
    factor.levels <- levels(data[,subgroup.char[k]])
    # Add interaction test for heterogeneity
    result <- rbind(result,"")
    result[nrow(result),1] <- subgroup.char[k]
    ia.pval <- anova(coxph(ia.model,data=data),coxph(main.model,data=data),test="Chisq")[2,"P(>|Chi|)"]
    result[nrow(result),ncol(result)] <- format.pval(ia.pval)
    # Add results for each subgroup level
    for (j in 1:length(factor.levels)){
      result <- rbind(result,"")
      result[nrow(result),1] <- paste("-",factor.levels[j])
      d.subgroup <- subset(data,.subgroup.var==factor.levels[j])
      result[nrow(result),2:(ncol(result)-1)] <- surv.comparison(model=base.model,data=d.subgroup,...)[3,]
    }
  }
  result
}




surv.comparison.abs.risk <- function(model,data,time){
  ##---------------------------------------------------------------------------------------------------------
  ## Purpose: Formally compare risks at a specific time point based on Greenwood estimates
  ## !! model formula should only include the survival model
  ## add.risk: if T, the event probability ("absolute risk") at time "infinity" is also displayed
  ## add.prop.haz.test: if T, a test for proportional hazards is also added
  ## Author: Marcel Wolbers, 9March2015
  ##---------------------------------------------------------------------------------------------------------
  arm.names <- levels(data$arm)
  # Table header
  header1 <- c(paste(arm.names,"(n=",table(data$arm),")",sep=""),"Comparison")
  header2 <- c(rep("events/n (risk [%])",2),"Absolute risk difference [%] (95%CI);p-value")
  header <- rbind(header1,header2)
  result <- rbind(header,"")
  # add number of events and risks
  fit.surv <- summary(survfit(update(model,.~arm),data=data),time=time,extend=T)
  events.n <- paste(fit.surv$table[,"events"],fit.surv$table[,"n.max"],sep="/")
  events.n <- paste(events.n," (",formatC(100*(1-fit.surv$surv),2,format="f"),")",sep="")
  result[3,1:2] <- events.n
  # add comparison
  diff <- (1-fit.surv$surv[2])-(1-fit.surv$surv[1])
  se.diff <- sqrt(sum(fit.surv$std.err^2))
  ci.diff <- c(diff-qnorm(0.975)*se.diff,diff+qnorm(0.975)*se.diff)
  pval.diff <- 2*pnorm(-abs(diff)/se.diff)
  result[3,3] <- paste(formatC(100*diff,2,format="f"),
                       "(",formatC(100*ci.diff[1],2,format="f")," to ",formatC(100*ci.diff[2],2,format="f"),
                       "); p=",format.pval(pval.diff),sep="")
  colnames(result) <- result[1,]; rownames(result) <- NULL
  result[-1,]
}

## Treatment comparisons overall and in subgroups - Logistic regression
logist.comparison <- function(model,data){
  ##---------------------------------------------------------------------------------------------------------
  ## Purpose: Summarize results for a logistic regression model ith the treatment arm (variable "arm") as the main covariate
  ## !! model formula can include other covariates than arm BUT arm must be the first covariate in the model !!
  ## Author: Marcel Wolbers, 16March2015
  ##---------------------------------------------------------------------------------------------------------
  arm.names <- levels(data$arm)
  # Table header
  header1 <- c(paste(arm.names,"(n=",table(data$arm),")",sep=""),"Comparison")
  header2 <- c(rep("events/n (%)",2),"OR (95%CI);p-value")
  header <- rbind(header1,header2)
  result <- rbind(header,"")
  # add number of events and risks
  tmp <- model.frame(model,data)
  tab <- table(model.response(tmp),tmp$arm)
  events.n <- paste(tab[2,],tab[1,]+tab[2,],sep="/")
  events.n <- paste(events.n," (",formatC(100*tab[2,]/(tab[1,]+tab[2,]),2,format="f"),")",sep="")
  result[3,1:2] <- events.n
  # add OR, CI, p-value
  fit.glm <- glm(model,data,family=binomial)
  or <- formatC(exp(fit.glm$coef[2]),2,format="f") # first coef is intercept, second is treatment arm
  pval <- format.pval(anova(fit.glm,update(fit.glm,.~.-arm),test="Chisq")[2,"Pr(>Chi)"]) # LR test
  profile.ci <- confint(fit.glm)
  ci <- paste(formatC(exp(profile.ci[2,1]),2,format="f"),
              formatC(exp(profile.ci[2,2]),2,format="f"),sep="-")
  or.ci.p <- paste(or,"(",ci,"); p=",pval,sep="")
  result[3,3] <- or.ci.p
  rownames(result) <- NULL
  result
}

logist.comparison.subgroup <- function(base.model,subgroup.model,data){
  ##---------------------------------------------------------------------------------------------------------
  ## Purpose: Summarize results for a logistic model by treatment arm (variable "arm") and subgroup
  ## base.model: Model from which sub-group specific estimates are extracted (!! arm must be the first covariate in the model)
  ## subgroup.model: model of the form "~subgrouping.variable1+subgrouping.variable2" 
  ##                 (!! subgrouping.variable must be factors and there should be nothing on the left-hand side of the formula)
  ## Author: Marcel Wolbers, 20March2015
  ##---------------------------------------------------------------------------------------------------------
  # result in entire population
  result <- logist.comparison(base.model,data)
  result <- cbind(c("Subgroup","","All patients"),result,c("Test for heterogeneity","p-value",""))
  # Preparation of models and data
  subgroup.char <- attr(terms(subgroup.model),"term.labels")
  for (k in 1:length(subgroup.char)){
    main.model <- update(base.model,as.formula(paste(".~.+",subgroup.char[k],sep=""))) 
    ia.model <- update(base.model,as.formula(paste(".~.+arm*",subgroup.char[k],sep=""))) 
    data$.subgroup.var <- data[,subgroup.char[k]]
    factor.levels <- levels(data[,subgroup.char[k]])
    # Add interaction test for heterogeneity
    result <- rbind(result,"")
    result[nrow(result),1] <- subgroup.char[k]
    ia.pval <- anova(glm(ia.model,data=data,family=binomial),
                     glm(main.model,data=data,family=binomial),test="Chisq")[2,"Pr(>Chi)"]
    result[nrow(result),ncol(result)] <- format.pval(ia.pval)
    # Add results for each subgroup level
    for (j in 1:length(factor.levels)){
      result <- rbind(result,"")
      result[nrow(result),1] <- paste("-",factor.levels[j])
      d.subgroup <- subset(data,.subgroup.var==factor.levels[j])
      result[nrow(result),2:(ncol(result)-1)] <- logist.comparison(model=base.model,data=d.subgroup)[3,]
    }
  }
  result
}

summary.MIresult <- function (object, ..., alpha = 0.05, logeffect = FALSE){
  #Overwrite the summary.MIresult from library(mitools):
  # add p-value and slightly modify reporting in case of logeffect=T
  #Last change: 8Aug2007 (corrected p-value calculation; use t instead of normal distribution)
  cat("Multiple imputation results:\n")
  lapply(object$call, function(a) {
    cat("      ")
    print(a)
  })
  out <- data.frame(results = coef(object), se = sqrt(diag(vcov(object))))
  crit <- qt(alpha/2, object$df, lower = FALSE)
  out$"pVal" <- round(2*pt(-abs(out$results)/out$se,df=object$df),5)
  if (logeffect) { out$expresults <- exp(out$results) }
  out$"(lower" <- out$results - crit * out$se
  out$"upper)" <- out$results + crit * out$se
  if (logeffect) {
    out$"(lower" <- exp(out$"(lower")
    out$"upper)" <- exp(out$"upper)")
  }
  out$missInfo <- paste(round(100 * object$missinfo), "%")
  print(out, ...)
}


# Logistic regressions  
myextractlogit <- function(model,data){
  # OR, CI, p-value
  fit <- glm(model,data,family=binomial)
  OR <- round(exp(coef(fit)[-1]),2)
  pval <- round(summary(fit)$coef[-1,"Pr(>|z|)"],5)
  ci <- round(exp(confint.default(fit)[-1,]),2)            
  or.ci.p <- paste(OR,"(",ci[1],",",ci[2],"); p=",pval,sep="")
  or.ci.p
}

auc.cal <- function(x,y,min,max,detectionLimit=-Inf){
  # Caculate AUC for an x-range from min to max based on linear interpolation/extrapolation of given x and y values
  # !! x, min and max must be integers for this AUC calculation to be correct
  val <- approxExtrap(x=x, y=y,xout=min:max)
  y <- pmax(val$y,detectionLimit)
  sol <- sum(y)-y[1]/2-tail(y,1)/2 
  c("AUC"=sol)
}
fit.S.sd<-function(data.long,data.surv,n.adapt=1000,n.iter=2000,n.chain=3,max_treedepth=13,adapt_delta = 0.9,seed=10){
  set.seed(seed)
  find_index<-function(ttdeath,time_spec){
    index<-rep(1,length(ttdeath))
    for(i in 1:length(ttdeath)){
      for(j in 1:(length(time_spec)-1)){
        if((ttdeath[i]<=time_spec[j+1])&(ttdeath[i]>time_spec[j])){
          index[i]<-j
          
        }
      }
    }
    return(index)
  }
  # longitudinal dataset
  N_long<-nrow(data.long)# number of patients
  data.long$y2.censInd<-(data.long$log10.csf.fungalcount<=log10(4.5))+0
  data.long$N_order<- c(1:N_long)
  index_cens<-data.long$N_order[data.long$y2.censInd==1]
  ncens<-length(index_cens)
  index_obs<-data.long$N_order[data.long$y2.censInd==0]
  nobs<-length(index_obs)
  data.longitudinal_obs<-data.long[data.long$N_order%in%index_obs,]
  data.longitudinal_cens<-data.long[data.long$N_order%in%index_cens,]
  data.long<-rbind(data.longitudinal_obs,data.longitudinal_cens)
  I<-nrow(data.surv)
  C<-log10(4.5)
  log10fc_obs<-data.long$log10.csf.fungalcount[1:nobs]
  time<-data.long$new.day
  patid<-as.integer(as.factor(data.long$usubjid)) 
  X_inter_long<-model.matrix(~1,data=data.long)
  p_X_intercept<-ncol(X_inter_long)
  X_slope_long<-model.matrix(~arm,data=data.long)
  p_X_slope<-ncol(X_slope_long)
  # survival dataset
  #X<-model.matrix(formula.longitudinal,data=data.long)
  X_intercept<-model.matrix(~1,data=data.surv)
  X_slope<-model.matrix(~I(arm=="Control"),data=data.surv)
  Z <- model.matrix(~1,data=data.surv)
  p_Z<-ncol(Z)
  death<-data.surv$evdeath.week10
  ttdeath<-data.surv$ttdeath.week10
  time_spec<-c(0,7,14,21,28,42,71) #quantile(ttdeath,c(0.0025,0.005,0.025,0.5,0.75,0.95,0.975))
  n_time_interval<-length(time_spec)-1
  index_interval<-find_index(ttdeath,time_spec)
  zeros<-matrix(rep(0,n_time_interval*I),nrow=I,ncol=n_time_interval)
  time_mean<-mean(data.long$studyday)
  time_sd<-sd(data.long$studyday)
  CM.data<-c("time_mean","time_sd","N_long","ncens","nobs","I","patid","time","log10fc_obs","C","p_X_intercept","p_X_slope",
             "X_intercept","X_slope","p_Z","Z","ttdeath","n_time_interval","time_spec","death","zeros","index_interval")
  fit <- stan(file ='V:/BIOSTATISTICS/RCT/CN_CNS_HIV/28CN_Tamoxifen_RCT/Analysis/Ranalysis/fit_S_scale.stan', data = CM.data,warmup=n.adapt,thin=10,
              iter = n.iter, chains = n.chain,pars=c("a_fix","b_fix","b_fix_unscale","b_2"),control=list(max_treedepth=max_treedepth,adapt_delta = adapt_delta))
  return(fit)
}

fit.S.cohort<-function(formula.longitudinal,data.long,formula.survival,data.surv,n.adapt=3000,n.iter=4000,n.chain=3,max_treedepth=14,adapt_delta = 0.95,seed=10,timepoint=c(0,2,4,7,10,13,17)){
  data.long<-arrange(data.long,usubjid,day)
  data.surv<-arrange(data.surv,usubjid)
  set.seed(seed)
  find_index<-function(ttdeath,time_spec){
    index<-rep(1,length(ttdeath))
    for(i in 1:length(ttdeath)){
      for(j in 1:(length(time_spec)-1)){
        if((ttdeath[i]<=time_spec[j+1])&(ttdeath[i]>time_spec[j])){
          index[i]<-j
          
        }
      }
    }
    return(index)
  }
  scale_matrix<-function(a){
    mean.a<-apply(t(a),1,mean)
    sd.a<-apply(t(a),1,sd)
    b<-a
    for(i in 2:ncol(a)){
      b[,i]<-(a[,i]-mean.a[i])/sd.a[i]
    }
    return(b)
  }
  # longitudinal dataset
  N_long<-nrow(data.long)# number of patients
  y2_censInd<-(data.long$log10.csf.fungalcount<=log10(4.5))+0
  I<-nrow(data.surv)
  C<-log10(4.5)
  log10fc_obs<-data.long$log10.csf.fungalcount
  patid<-as.integer(as.factor(data.long$usubjid)) 
  X_long_scale<-scale_matrix(model.matrix(formula.longitudinal,data=data.long))
  p_X<-ncol(X_long_scale)
  time_scale<-X_long_scale[,"studyday"]
  X_intercept_unscale<-model.matrix(~1,data=data.surv)
  X_slope_unscale<-model.matrix(~arm,data=data.surv)
  if(nlevels(data.surv$arm)==1){
    X_slope_unscale<-model.matrix(~1,data=data.surv)
  }
  p_X_intercept<-ncol(X_intercept_unscale)
  p_X_slope<-ncol(X_slope_unscale)
  Z <- model.matrix(formula.survival,data=data.surv)
  p_Z<-ncol(Z)
  data.surv$evdeath.17d<-ifelse(data.surv$ttdeath.week10<17 & data.surv$evdeath.week10==1,1,0)
  data.surv$ttdeath.17d<-pmin(data.surv$ttdeath.week10,17)
  death<-data.surv$evdeath.17d
  ttdeath<-data.surv$ttdeath.17d
  time_spec<-timepoint
  n_time_interval<-length(time_spec)-1
  index_interval<-find_index(ttdeath,time_spec)
  zeros<-matrix(rep(0,n_time_interval*I),nrow=I,ncol=n_time_interval)
  mean_X<-apply(t(model.matrix(formula.longitudinal,data=data.long)),1,mean)
  sd_X<-apply(t(model.matrix(formula.longitudinal,data=data.long)),1,sd)
  
  CM.data<-c("mean_X","sd_X","N_long","y2_censInd","I","patid","time_scale","log10fc_obs","C","p_X","p_X_intercept","p_X_slope",
             "X_long_scale","X_intercept_unscale","X_slope_unscale","p_Z","Z","ttdeath","n_time_interval","time_spec","death","zeros","index_interval")
  fit <- stan(file = "V:/BIOSTATISTICS/RCT/CN_CNS_HIV/28CN_Tamoxifen_RCT/Analysis/Ranalysis/fit_S_ran_eff_12_Jan19.stan", data = CM.data,warmup=n.adapt,thin=1,
              iter = n.iter, chains = n.chain,pars=c("a_fix_unscale","b_fix_unscale","b_fix_unscale_true","beta","lambda","eta","sigma","lp__","a_unscale","b_unscale"),control=list(max_treedepth=max_treedepth,adapt_delta = adapt_delta))

  # fit <- stan(file = "V:/BIOSTATISTICS/RCT/CN_CNS_HIV/28CN_Tamoxifen_RCT/Analysis/Ranalysis/fit_S_ran_eff_12_Jan19.stan", data = CM.data,warmup=n.adapt,thin=10,
  #             iter = n.iter, chains = n.chain,pars=c("a_fix_unscale","b_fix_unscale","b_fix_unscale_true","beta"),control=list(max_treedepth=max_treedepth,adapt_delta = adapt_delta))
  
    return(fit)
}
declineEstimates.unbias <- function(data.long,data.surv,n.adapt=2500,n.iter=5000,n.chain=3,max_treedepth=15,adapt_delta = 0.97,timepoint=c(0,2,4,7,10,13,17) ){
  data.long$usubjid<-droplevels(data.long$usubjid)
  data.surv$usubjid<-droplevels(data.surv$usubjid)
  data.long$new.day<-(data.long$studyday-mean(data.long$studyday))/sd(data.long$studyday)
  # fit.1<-fit.S.sd(data.long=data.long,data.surv=data.surv,n.adapt=n.adapt,n.iter=n.iter,n.chain=n.chain,max_treedepth=max_treedepth,adapt_delta = adapt_delta)
  fit.1<-fit.S.cohort(formula.longitudinal=~studyday+studyday:arm, data.long=data.long,
                      formula.survival=~1, data.surv=data.surv, n.adapt=n.adapt, n.iter=n.iter, n.chain=n.chain,
                      max_treedepth=max_treedepth, adapt_delta = adapt_delta,timepoint = timepoint) 
  res.1<-as.data.frame(summary(fit.1)$summary)
  n.arm1 <- length(unique(data.long$usubjid[data.long$arm==lev.arm[1]])) 
  n.arm2 <- length(unique(data.long$usubjid[data.long$arm==lev.arm[2]])) 
  change.arm1 <- paste(formatC(res.1[4,"mean"],3,format="f")," (",
                       formatC(res.1[4,"2.5%"],3,format="f"),",",
                       formatC(res.1[4,"97.5%"],3,format="f"),")",sep="")
  change.arm2 <- paste(formatC(res.1[5,"mean"],3,format="f")," (",
                       formatC(res.1[5,"2.5%"],3,format="f"),",",
                       formatC(res.1[5,"97.5%"],3,format="f"),")",sep="")  
  diff.1vs2 <-   paste(formatC(res.1[3,"mean"],3,format="f")," (",
                       formatC(res.1[3,"2.5%"],3,format="f"),",",
                       formatC(res.1[3,"97.5%"],3,format="f"),")",
                       format.pval(2*pnorm(-abs(res.1[3,"mean"]/res.1[3,"sd"]))),sep="")
  list(n.arm1=n.arm1,change.arm1=change.arm1,n.arm2=n.arm2,change.arm2=change.arm2,diff.1vs2=diff.1vs2,obj=fit.1)
}
declineEstimates.DL <- function(outcome.model,data){
  data$usubjid<-droplevels(data$usubjid)
  y <- data$log10.csf.fungalcount
  cens<-(y<=log10(4.5))+0
  X<-model.matrix(update(outcome.model,.~studyday+arm1.studyday),data=data)
  Z <- model.matrix(log10.csf.fungalcount~1+studyday,data=data)
  cluster = as.numeric(as.factor(data$usubjid))
  fit.1<-lmec(yL=y,cens=cens, X=X, Z=Z, cluster=cluster, method='ML',abspmv=1e-6, maxstep=1000000)
  n.arm1 <- length(unique(data$usubjid[data$arm==lev.arm[1]])) 
  change.arm2 <- paste(formatC(fit.1$beta[2],3,format="f")," (",
                       formatC(fit.1$beta[2]-qnorm(0.975,0,1)* sqrt(fit.1$varFix[2,2]),3,format="f"),",",
                       formatC(fit.1$beta[2]+qnorm(0.975,0,1)* sqrt(fit.1$varFix[2,2]),3,format="f"),")",sep="")
  
  X<-model.matrix(update(outcome.model,.~studyday+arm2.studyday),data=data)
  fit.2<-lmec(yL=y,cens=cens, X=X, Z=Z, cluster=cluster, method='ML',abspmv=1e-6, maxstep=1000)
  n.arm2 <- length(unique(data$usubjid[data$arm==lev.arm[2]])) 
  change.arm1 <- paste(formatC(fit.2$beta[2],3,format="f")," (",
                       formatC(fit.2$beta[2]-qnorm(0.975,0,1)* sqrt(fit.2$varFix[2,2]),3,format="f"),",",
                       formatC(fit.2$beta[2]+qnorm(0.975,0,1)* sqrt(fit.2$varFix[2,2]),3,format="f"),")",sep="")
  diff.1vs2 <-   paste(formatC(fit.2$beta[3],3,format="f")," (",
                       formatC(fit.2$beta[3]-qnorm(0.975,0,1)* sqrt(fit.2$varFix[3,3]),3,format="f"),",",
                       formatC(fit.2$beta[3]+qnorm(0.975,0,1)* sqrt(fit.2$varFix[3,3]),3,format="f"),"); p value ",
                       format.pval(2*pnorm(-abs(fit.2$beta[3]/sqrt(fit.2$varFix[3,3])))),sep="")
  list(n.arm1=n.arm1,change.arm1=change.arm1,n.arm2=n.arm2,change.arm2=change.arm2,diff.1vs2=diff.1vs2)
}
declineEstimates <- function(outcome.model,data){
  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(coef(fit.lme2)["studyday","Estimate"],3,format="f")," (",
                       formatC(coef(fit.lme2)["studyday","Estimate"]-qnorm(0.975)*coef(fit.lme2)["studyday","Std. Error"],3,format="f"),",",
                       formatC(coef(fit.lme2)["studyday","Estimate"]+qnorm(0.975)*coef(fit.lme2)["studyday","Std. Error"],3,format="f"),")",sep="")
  change.arm2 <- paste(formatC(coef(fit.lme1)["studyday","Estimate"],3,format="f")," (",
                       formatC(coef(fit.lme1)["studyday","Estimate"]-qnorm(0.975)*coef(fit.lme1)["studyday","Std. Error"],3,format="f"),",",
                       formatC(coef(fit.lme1)["studyday","Estimate"]+qnorm(0.975)*coef(fit.lme1)["studyday","Std. Error"],3,format="f"),")",sep="")
  diff.1vs2 <-   paste(formatC(coef(fit.lme2)["arm2.studyday","Estimate"],3,format="f")," (",
                       formatC(coef(fit.lme2)["arm2.studyday","Estimate"]-qnorm(0.975)*coef(fit.lme2)["arm2.studyday","Std. Error"],3,format="f"),",",
                       formatC(coef(fit.lme2)["arm2.studyday","Estimate"]+qnorm(0.975)*coef(fit.lme2)["arm2.studyday","Std. Error"],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)
}

vcov_fit<-function(obj){
  list_of_draws <- rstan::extract(obj)
  M<-cbind(a=list_of_draws$a_fix_unscale,b=list_of_draws$b_fix_unscale)
  k <- ncol(M) #number of variables
  n <- nrow(M) #number of subjects
  
  #create means for each column
  M_mean <- matrix(data=1, nrow=n) %*% colMeans(M) 
  
  #creates a difference matrix
  D <- M - M_mean
  
  #creates the covariance matrix
  C <- (n-1)^-1*t(D) %*% D
  return(C)
}

pred.lmer.diff_msm<-function(obj,newdata){
  library(msm)
  dat<-newdata
  
  for(i in 1:nrow(dat)){
    mm<-model.matrix(~(ns(day,kn=knots.time,Bo = c(knots.min,knots.max))+ns(day,kn=knots.time,Bo = c(knots.min,knots.max)):arm)*tamofluco.measure.time,data=dat[i,])
    a1<-mm[1]
    a2<-mm[2]
    a3<-mm[3]
    a4<-mm[4]
    a5<-mm[5]
    a6<-mm[6]
    a7<-mm[7]
    a8<-mm[8]
    a9<-mm[9]
    a10<-mm[10]
    a11<-mm[11]
    a12<-mm[12]
    a13<-mm[13]
    a14<-mm[14]

    pse_before<-deltamethod(~(a6*x6+a7*x7+a8*x8),
                      mean=fixef(fit_treatment), cov=vcov(fit_treatment))
    pse_after<-deltamethod(~(a12*x12+a13*x13+a14*x14),
                           mean=fixef(fit_treatment), cov=vcov(fit_treatment))
    
    pse_Control<-deltamethod(~(a5*x5+a9*x9+a10*x10+a11*x11),
                            mean=fixef(fit_treatment), cov=vcov(fit_treatment))
    pse_Tamox<-deltamethod(~(a5*x5+a12*x12+a13*x13+a14*x14),
                             mean=fixef(fit_treatment), cov=vcov(fit_treatment))

    dat$absolute_diff_before[i]<-mm[1,6:8]%*%fixef(fit_treatment)[6:8]
    dat$absolute_diff_before_lo[i]<-dat$absolute_diff_before[i]-qnorm(0.975,0,1)*pse_before
    dat$absolute_diff_before_hi[i]<-dat$absolute_diff_before[i]+qnorm(0.975,0,1)*pse_before
    dat$absolute_diff_after[i]<-mm[1,12:14]%*%fixef(fit_treatment)[12:14]
    dat$absolute_diff_after_lo[i]<-dat$absolute_diff_after[i]-qnorm(0.975,0,1)*pse_after
    dat$absolute_diff_after_hi[i]<-dat$absolute_diff_after[i]+qnorm(0.975,0,1)*pse_after
    
    dat$absolute_diff_Control[i]<-mm[1,c(5,9:11)]%*%fixef(fit_treatment)[c(5,9:11)]
    dat$absolute_diff_Control_lo[i]<-dat$absolute_diff_Control[i]-qnorm(0.975,0,1)*pse_Control
    dat$absolute_diff_Control_hi[i]<-dat$absolute_diff_Control[i]+qnorm(0.975,0,1)*pse_Control
    dat$absolute_diff_Tamox[i]<-mm[1,c(5,12:14)]%*%fixef(fit_treatment)[c(5,12:14)]
    dat$absolute_diff_Tamox_lo[i]<-dat$absolute_diff_Tamox[i]-qnorm(0.975,0,1)*pse_Tamox
    dat$absolute_diff_Tamox_hi[i]<-dat$absolute_diff_Tamox[i]+qnorm(0.975,0,1)*pse_Tamox
    
    
  }
  return(dat)
}