##############################################################################################
##############################################################################################
##############################################################################################
##############################################################################################
###  main analysis
### R --no-save < report.R > report.lst
##############################################################################################
##############################################################################################
##############################################################################################
##############################################################################################
 
source("/home/edensk/myStuff/mySources/statMisc.R")
source("/home/edensk/myStuff/mySources/skeMisc.R")
source("/home/edensk/myStuff/mySources/logranktimevariate/timevaryinglogrank20090930.R")
source("/home/edensk/myStuff/mySources/latextablezoo.s")
source("usefulFun.R")

library(chron)
library(Design)
library(Hmisc)
library(lattice)
library(survival)
outcomeLevels = c("lactAcidosis", "clinPancreatitis", "Both")

###----------------------------------LOAD-----------------------------------------
load("../data/data.rda")
load("../data/spssdata.rda")
load("../data/lacdata.rda")
load("../data/coxdata.rda")
load("../data/podata.rda")
load("../data/perid.rda")

coxdata$bmi3 = NA
coxdata$bmi3 = ifelse(coxdata$bmi<18.5, "<18.5", coxdata$bmi3)
coxdata$bmi3 = ifelse(coxdata$bmi>=18.5 & coxdata$bmi<=25, "18.5-25", coxdata$bmi3)
coxdata$bmi3 = ifelse(coxdata$bmi>25, ">25", coxdata$bmi3)
coxdata$bmi2 = NA
coxdata$bmi2 = ifelse(coxdata$bmi<=25, "<=25", coxdata$bmi2)
coxdata$bmi2 = ifelse(coxdata$bmi>25, ">25", coxdata$bmi2)
coxdata$cd4_3 = NA
coxdata$cd4_3 = ifelse(coxdata$cd4<100, "<100", coxdata$cd4_3)
coxdata$cd4_3 = ifelse(coxdata$cd4>=100 & coxdata$cd4<=200, "100-200", coxdata$cd4_3)
coxdata$cd4_3 = ifelse(coxdata$cd4>200, ">200", coxdata$cd4_3)
origCoxData = coxdata

###----------------------------------SINGLE IMPUTATION----------------------------
coxdata$hemglbn[is.na(coxdata$hemglbn)] = median(coxdata$hemglbn, na.rm=TRUE)
coxdata$bmi[is.na(coxdata$bmi)] = median(coxdata$bmi, na.rm=TRUE)
coxdata$logbmi = log(coxdata$bmi)
coxdata$viralload[is.na(coxdata$viralload)] = median(coxdata$viralload, na.rm=TRUE)
coxdata$logviralload10 = log(coxdata$viralload, 10)
coxdata$logage = log(coxdata$age)
coxdata$sqrtcd4 = sqrt(coxdata$cd4)
coxdata$loghemglbn = log(coxdata$hemglbn)
coxdata$logviralload10 = log(coxdata$viralload, 10)
coxdata$txnucl = factor(coxdata$txnucl)
coxdata$txnonnucl = factor(coxdata$txnonnucl)
coxdata$txadher = factor(coxdata$txadher)

###-----------------------------------Histograms-----------------------------------
pdftex(pathToPdf="pdf", pathToTex="gentex", "dataHist", "Histograms. Compare Transformations.", longCaption="Histograms. Compare Transformations.", label=NULL, width=6, height=7, fontsize="normalsize")
par(mfrow=c(5, 2), mar=c(4,2,2,1))
for (n in c("hemglbn", "loghemglbn", "bmi", "logbmi", "viralload", "logviralload10", "age", "logage", "cd4", "sqrtcd4")){
  hist(coxdata[[n]], main=n, breaks=50, xlab="")
}
dev.off()

coxdata$hemglbn[is.na(coxdata$hemglbn)] = median(coxdata$hemglbn, na.rm=TRUE)
coxdata$bmi[is.na(coxdata$bmi)] = median(coxdata$bmi, na.rm=TRUE)
coxdata$viralload[is.na(coxdata$viralload)] = median(coxdata$viralload, na.rm=TRUE)
coxdata$logviralload10 = log(coxdata$viralload, 10)
coxdata$logage = log(coxdata$age)
coxdata$sqrtcd4 = sqrt(coxdata$cd4)
coxdata$loghemglbn = log(coxdata$hemglbn)
coxdata$logviralload10 = log(coxdata$viralload, 10)
coxdata$txnucl = factor(coxdata$txnucl)
coxdata$txnonnucl = factor(coxdata$txnonnucl)
coxdata$txadher = factor(coxdata$txadher)

perid$start = 0
perid$cd4_3 = NA
perid$cd4_3 = ifelse(perid$cd4<200, "<200", perid$cd4_3)
perid$cd4_3 = ifelse(perid$cd4>=200 & perid$cd4<=250, "200-250", perid$cd4_3)
perid$cd4_3 = ifelse(perid$cd4>250, ">250", perid$cd4_3)
perid$sexcd4_3 = paste(perid$sex, perid$cd4_3)
perid$timetoskintoxicity = perid$timetoskint
perid$timetoskintoxicity[perid$skintoxi==0] = NA
perid$sktoxlabels = ifelse(perid$skintoxi==1, "Toxicity", "No toxicity")

###----------------------------------summarize basics-----------------------------------------
summ = summary( ~ age + sex + cd4 + cd4_3 + viralload + hemglbn + karnof + bmi + bmi2 + bmi3, data = coxdata, method="reverse", continuous = 5, nmin=1)
trash=latex(summ, prtest='P', file=paste('gentex/', "summary1", '.tex', sep=''), exclude1=FALSE, prmsd=TRUE, caption=paste("Basic Summary"), dotchart=FALSE, digits=2, where="hbp!")

table(data$drugrsn)
table(data$lactacid)
table(data$drugrsn[data$outcome==1])

###----------------------------------missing variables-----------------------------------------
sinktex(paste("gentex/", "numOfMissingVar", '.tex', sep=''))
  round(missingVariables(origCoxData), 1)
sinktex()

#################################################################################
#################################################################################
#################################################################################
################ SURVIVAL ANALYSIS FOR LACTIC ACIDOSIS
#################################################################################
#################################################################################
#################################################################################
dd = datadist(coxdata)
options(datadist="dd")
survObj = Surv(time=coxdata$timetoevent, event=coxdata$coxevent)
coximpdata = multipleImpForCPH(coxdata, timeVar="timetoevent", eventVar="coxevent", 10)

#################################################################################
################ unadjusted HRs for Lactic Acidosis For Not Imputed
#################################################################################
dd = datadist(origCoxData)
options(datadist="dd")

unitVector = c(txnucl="", txnonnucl="", txadher="", age="year", sex="", cd4="cell", cd4_3="", logviralload10="log unit", hemglbn="unit", karnof="unit", bmi="unit", bmi2="")

varListOfLabel= c(txnucl="Nucleosides", txnonnucl="Non-Nucleosides", txadher="Adherence", age="Age", sex="Gender", bmi="BMI, continuous", bmi2="BMI, 2 categories", cd4="CD4, continuous", cd4_3="CD4, 3 categories", logviralload10="Log viral load", hemglbn="Hemoglobin", karnof="KARNOFSKY")

tmp = subset(origCoxData, sex=="female")

hrMatrix = unadjustedHRsIMPNew(list(tmp), timeVar="timetoevent", eventVar="coxevent", varListOfInput=list(txnucl="ZDV/3TC", txnonnucl="EFV", txadher="Com-DOT", age=c(60, 70), bmi=c(19,20), bmi2="<=25", cd4=c(100, 200), cd4_3=">200", logviralload10=c(5, 6), hemglbn=c(8, 9), karnof=c(10, 20)), varListOfLabel = varListOfLabel, unitVector=unitVector)
hrTable = prepareHRorORtableForZebraTable(hrMatrix)

# zebraTable("gentex/unadjustedHRsLactNotImp.tex",
#            caption = paste("Unadjusted HR for Lactic Acidosis Using Original (not Imputed) Data.", nrow(tmp), "subjects"), orderGroups=FALSE, dataframe=hrTable, zebraPattern="none", fontSize="normalsize", by ="id", vars=setdiff(names(hrTable), c("id")), colNames = latexTextMode(c("", "HR", "95% CIs", "P-value")), landscape=FALSE, toLatexChar=FALSE)

#################################################################################
################ unadjusted HRs for Lactic Acidosis for Imputed
#################################################################################
coximpdatawomen = coximpdata
for (i in 1:length(coximpdata)){
  coximpdatawomen[[i]] = subset(coximpdata[[i]], sex=="female")
  tmp = subset(coximpdata[[i]], sex=="female")
}

dd = datadist(coximpdata[[1]])
options(datadist="dd")

unitVector = c(txnucl="", txnonnucl="", txadher="", age="year", sex="", imbmi="unit", bmi2="unit", cd4="cell", cd4_3="", log10viralload="log unit", imhemglbn="unit", karnof="unit")

varListOfLabel= c(txnucl="Nucleosides", txnonnucl="Non-Nucleosides", txadher="Adherence", age="Age", sex="Gender", imbmi="BMI, continuous", bmi2="BMI, 2 categories", cd4="CD4, continuous", cd4_3="CD4, 3 categories", log10viralload="Log viral load", imhemglbn="Hemoglobin", karnof="KARNOFSKY")

imhrMatrix = unadjustedHRsIMPNew(coximpdatawomen, timeVar="timetoevent", eventVar="coxevent", varListOfInput=list(txnucl="ZDV/3TC", txnonnucl="EFV", txadher="Com-DOT", age=c(60, 70), imbmi=c(19, 20), bmi2="<=25", cd4=c(100, 200), cd4_3=">200", log10viralload=c(5, 6), imhemglbn=c(9, 8), karnof=c(10, 20)), varListOfLabel = varListOfLabel, unitVector=unitVector)
hrTable = prepareHRorORtableForZebraTable(imhrMatrix)

# zebraTable("gentex/unadjustedHRsLact.tex",
#            caption = paste("Unadjusted HR for Lactic Acidosis Using Imputed Data.", nrow(coximpdatawomen[[1]]), "subjects."), orderGroups=FALSE, dataframe=hrTable, zebraPattern="none", fontSize="normalsize", by ="id", vars=setdiff(names(hrTable), c("id")), colNames = latexTextMode(c("", "HR", "95% CIs", "P-value")), landscape=FALSE, toLatexChar=FALSE)


#################################################################################
################ checking non-linearity for lactic acidosis
#################################################################################
coxdatawomen = subset(coxdata, sex=="female")
nonlinvars = c("age", "bmi", "cd4", "hemglbn")
nonlinpvalies = pvaluesForNonLinear(coxdatawomen, c("age", "bmi", "cd4", "hemglbn"), timeVar="timetoevent", eventVar="coxevent")
sinktex(paste("gentex/", "checknonlin", '.tex', sep=''))
  cat("Non linearity was checked in women's subgroup for the following variables:\n", paste(nonlinvars, collapse=", "),"\n")
  if (min(nonlinpvalies) >= .05){
    cat("\nNone of the variables displayed non-linearity\n")
    cat(paste("Minimum P-value for non-linear term was ", round(min(nonlinpvalies), 2), "\n"))
  }else{
    varStr = paste(names(nonlinpvalies)[nonlinpvalies < .05], " (P-value=", pvalueStr(nonlinpvalies[nonlinpvalies < .05]), ")", sep="")
    cat("\nThe following variables displayed non-linearity:\n")
    varStr
  }
sinktex()
#################################################################################
################ cox regression on nucleosides
#################################################################################

cphmodels4 = runCPHModels(coximpdatawomen, formula=as.formula("Surv(timetoevent, coxevent) ~ rcs(age, 3) + imbmi + txnucl"))

unitVector = c(txnucl="", age="year", imbmi="unit")
varListOfLabel= c(age="Age", imbmi="BMI", sex="Gender", txnucl="Nucleosides")
varListOfInput=list(age=c(50, 60), imbmi=c(19, 20), txnucl="ZDV/3TC")

imhrMatrix = getNiceOutputORorHRfromSummary(models = cphmodels4, data=coximpdata[[1]], varListOfInput=varListOfInput, varListOfLabel=varListOfLabel, unitVector=unitVector)
hrTable = prepareHRorORtableForZebraTable(imhrMatrix)

# zebraTable("gentex/coxmodelforlactivasidosis.tex",
#            caption = "Survival Model Output for Lactic Acidosis.", orderGroups=FALSE, dataframe=hrTable, zebraPattern="none", fontSize="small", by ="id", vars=setdiff(names(hrTable), c("id")), colNames = latexTextMode(c("", "HR", "95% CIs", "P-value")), landscape=FALSE, toLatexChar=FALSE)

#################################################################################
#################################################################################
#################################################################################
################ proportional odds regression for lactic acidosis
#################################################################################
#################################################################################
#################################################################################
podata$logage = log(podata$age)
podata$logbmi = log(podata$bmi)
podata$sqrtcd4 = sqrt(podata$cd4)
podata$sex = factor(podata$sex)
podata$log10viralload = log(podata$viralload, base=10)
podata$loghemglbn = log(podata$hemglbn)
dd = datadist(podata)
options(datadist="dd")

pdftex(pathToPdf="pdf", pathToTex="gentex", "followuptimehist", "Follow-up Time Histogram", longCaption="Follow-up Time Histogram", label=NULL, width=8, height=4, fontsize="normalsize")
  par(mfrow=c(1, 1))
  hist(podata$timetolastalive, freq=TRUE, ylab="Count", xlab="Length of follow-up time, days", breaks=50, main="Histogram of follow-up time")
  #plot(podata$timetolastalive[order(podata$timetolastalive)], cex=.5, xlab="Subject Index", ylab="Follow-up Time")
dev.off()
sinktex(paste("gentex/", "POLogisticRegressionModel", '.tex', sep=''))
  cat(paste("\n\nEffective sample size for Proportional Odds regression is ", round(nrow(podata) - (1/(nrow(podata)^2))*sum(table(podata$poevent)^3)), ", which means that\nwe can explore the effect of about ", round((nrow(podata) - (1/(nrow(podata)^2))*sum(table(podata$poevent)^3))/15),  " covariates without concern of overfitting", sep=""), "\n")
sinktex()

#################################################################################
################ proportional odds regression with imputed data for lactic acidosis
#################################################################################
set.seed(1)
impdata = multipleImpForLRM(podata, 10)
pomodels = runPOLRModels(impdata, as.formula("poevent ~ age + imbmi + cd4 + log10viralload + imhemglbn + sex + txnucl + txnonnucl + txadher"))
dd = datadist(impdata[[1]])
options(datadist="dd")

#################################################################################
################ proportional odds regression with imputed data for lactic acidosis
#################################################################################
impdatawomen = impdata
for (i in 1:length(impdata)){
  impdatawomen[[i]] = subset(impdata[[i]], sex=="female")
}
dd = datadist(impdatawomen[[1]])
options(datadist="dd")
pomodelswomen = runPOLRModels(impdatawomen, as.formula("poevent ~ age + imbmi + cd4 + log10viralload + imhemglbn + txnucl + txnonnucl + txadher"))

unitVector = c(txnucl="", txnonnucl="", txadher="", age="year", imbmi="unit", cd4="cell", log10viralload="log unit", imhemglbn="unit")
varListOfLabel= c(txnucl="Nucleosides", txnonnucl="Non-Nucleosides", txadher="Adherence", age="Age", imbmi="BMI, continuous", cd4="CD4, continuous", log10viralload="Log viral load", imhemglbn="Hemoglobin")
#varListOfInput=list(txnucl="ZDV/3TC", txnonnucl="EFV", txadher="Com-DOT", age=c(60, 70), imbmi=c(19, 20), cd4=c(100, 200), log10viralload=c(5, 6), imhemglbn=c(8, 9))
varListOfInput=list(txnucl="ZDV/3TC", txnonnucl="EFV", txadher="SOC", age=c(60, 70), imbmi=c(19, 20), cd4=c(100, 200), log10viralload=c(5, 6), imhemglbn=c(9, 8))

imhrMatrix = getNiceOutputORorHRfromSummary(models = pomodelswomen, data=impdatawomen[[1]], varListOfInput=varListOfInput, varListOfLabel=varListOfLabel, unitVector=unitVector)
hrTable = prepareHRorORtableForZebraTable(imhrMatrix)

# zebraTable("gentex/adjustedORsAsidosisWomen.tex",
#            caption = "Proportional Odds Model Output for Lactic Acidosis in Female Population.", orderGroups=FALSE, dataframe=hrTable, zebraPattern="none", fontSize="normalsize", by ="id", vars=setdiff(names(hrTable), c("id")), colNames = latexTextMode(c("", "OR", "95% CIs", "P-value")), landscape=FALSE, toLatexChar=FALSE)

tmpdata = impdatawomen[[1]]
dd = datadist(tmpdata)
options(datadist="dd")
tmppomod = lrm(poevent ~ age + imbmi + cd4 + log10viralload + imhemglbn + txnucl + txnonnucl + txadher, data=tmpdata)
summary(tmppomod, txnucl="ZDV/3TC", txnonnucl="EFV", txadher="Com-DOT", age=c(60, 70), imbmi=c(19, 20), cd4=c(100, 200), log10viralload=c(5, 6), imhemglbn=c(9, 8))

#################################################################################
#################################################################################
################ additional text for lactic acidosis paper ######################
#################################################################################
#################################################################################
sinktex(paste("gentex/", "additionalTextLacticAcidosis", '.tex', sep=''), verbatim = FALSE)
  cat("The follow-up time is ", round(sum(origCoxData$timetoevent/365.25)), " person-years", "\n", sep="")
  cat("With the median follow-up time of ", round(quantile((origCoxData$timetoevent/7), probs=c(0.5))), " weeks [IQR ", paste(round(quantile(origCoxData$timetoevent/7, probs=c(.25, .75))), collapse="-"),"].\n", sep="")
  tmpids=podata$id[podata$poevent != 0]
  laeventsdata = lacdata[lacdata$id %in% tmpids,]
  cat(kama(laeventsdata$id), " study participants had a total of ", nrow(laeventsdata), " serum lactate levels drawn.\n", sep="")
  tmpids=podata$id[podata$poevent == 2]
  laeventsdata = lacdata[lacdata$id %in% tmpids,]
  cat(kama(laeventsdata$id), " having a defined study event had a total of ", nrow(laeventsdata), " serum lactates drawn.\n", sep="")
  tmpids=podata$id[podata$poevent == 2]
  laeventsdata = data[data$id %in% tmpids,]
  cat("All ", kama(laeventsdata$id), " were female and their median age was ", round(quantile((laeventsdata$age), probs=c(0.5))), " years [IQR ", paste(round(quantile(laeventsdata$age, probs=c(.25, .75))), collapse="-"),"],\n", sep="")
  tmpids=podata$id[podata$poevent == 2]
  laeventsdata = origCoxData[origCoxData$id %in% tmpids,]
  cat("and were on cART for a median of ", round(quantile(laeventsdata$timetoevent/30.4375, probs=c(0.5)), 1), " months [IQR ", paste(round(quantile(laeventsdata$timetoevent/30.4375, probs=c(.25, .75)), 1), collapse="-"),"].\n", sep="")
  tmpids=podata$id[podata$poevent == 2]
  laeventsdata = lacdata[lacdata$id %in% tmpids,]
  cat("At the time of their moderate-severe symptomatic hyperlactemia or lactic acidosis event these ", kama(laeventsdata$id), " patients had a median serum lactate of ", round(quantile((laeventsdata$lactate), probs=c(0.5))), " [IQR ", paste(round(quantile(laeventsdata$lactate, probs=c(.25, .75))), collapse="-"),"].\n", sep="")
  tmpids=podata$id[podata$poevent == 2]
  laeventsdata = data[data$id %in% tmpids,]
  cat("and a median BMI of ", round(quantile((laeventsdata$basebmi), probs=c(0.5))), " [IQR ", paste(round(quantile(laeventsdata$basebmi, probs=c(.25, .75))), collapse="-"),"].\n", sep="")
sinktex(verbatim = FALSE)

#################################################################################
################ checking nonlinearity for proportional odds regression
#################################################################################
nonlinpvalies = pvaluesForNonLinearPOModel(podata, c("age", "bmi", "cd4", "hemglbn", "log10viralload"), outcomeVar="poevent")
sinktex(paste("gentex/", "checknonlinPO", '.tex', sep=''))
  cat("Non linearity was checked for the following variables:\n", paste(nonlinvars, collapse=", "),"\n")
  if (min(nonlinpvalies) >= .05){
    cat("\nNone of the variables displayed non-linearity\n")
    cat(paste("Minimum P-value for non-linear term was ", round(min(nonlinpvalies), 2), "\n"))
  }else{
    varStr = paste(names(nonlinpvalies)[nonlinpvalies < .05], " (P-value=", pvalueStr(nonlinpvalies[nonlinpvalies < .05]), ")", sep="")
    cat("\nThe following variables are candidates for non-linear association with the outcome:\n")
    varStr
  }
sinktex()


unitVector = c(txnucl="", txnonnucl="", txadher="", age="year", imbmi="unit", cd4="cell", log10viralload="log unit", imhemglbn="unit", sex="")
varListOfLabel= c(txnucl="Nucleosides", txnonnucl="Non-Nucleosides", txadher="Adherence", age="Age", imbmi="BMI, continuous", cd4="CD4, continuous", log10viralload="Log viral load", imhemglbn="Hemoglobin", sex="Gender")
varListOfInput=list(txnucl="ZDV/3TC", txnonnucl="EFV", txadher="SOC", age=c(60, 70), imbmi=c(19, 20), cd4=c(100, 200), log10viralload=c(5, 6), imhemglbn=c(9, 8), sex="male")
imhrMatrix = getNiceOutputORorHRfromSummary(models = pomodels, data=impdata[[1]], varListOfInput=varListOfInput, varListOfLabel=varListOfLabel, unitVector=unitVector)
hrTable = prepareHRorORtableForZebraTable(imhrMatrix)

# zebraTable("gentex/adjustedORsAsidosis.tex",
#            caption = "Proportional Odds Model Output for Lactic Acidosis.", orderGroups=FALSE, dataframe=hrTable, zebraPattern="none", fontSize="normalsize", by ="id", vars=setdiff(names(hrTable), c("id")), colNames = latexTextMode(c("", "OR", "95% CIs", "P-value")), landscape=FALSE, toLatexChar=FALSE)

allStuffFromSummary = retrStuffFromSummaryLrmOrCphForAllVars(list(models=pomodels, data=impdata[[1]]), age=c(50, 60), cd4=c(300, 400), log10viralload=c(5, 6), imhemglbn=c(9,8))

 
##############################################################################################
##############################################################################################
##############################################################################################
##############################################################################################
###  data preparation
### R --no-save < import.R > import.lst
##############################################################################################
##############################################################################################
##############################################################################################
##############################################################################################

source("/home/edensk/myStuff/mySources/skeMisc.R")
source("/home/edensk/myStuff/mySources/statMisc.R")
source("/home/edensk/myStuff/mySources/latextablezoo.s")
library(Hmisc)
library(Design)
library(chron)

pancdata = read.csv("../data/csv/timetopancreatitis.csv", as.is=TRUE)
pancdata$bid = gsub(" +", "", pancdata$bid)
pancdata$date = gsub("[()]", "", pancdata$date)
pancdata$date = gsub(" ", "/", pancdata$date)
pancdata$date = chron(pancdata$date, format="d/m/year")

skintdata = read.csv("../data/csv/skintoxicity.csv", as.is=TRUE)
skintdata$bid = gsub(" +", "", skintdata$bid)
skintdata$date = gsub("[()]", "", skintdata$date)
skintdata$date = gsub(" ", "/", skintdata$date)
skintdata$date = chron(skintdata$date, format="d/m/year")

data = read.csv("../data/csv/tshepo_fulldata_final.csv", as.is=TRUE)
load("../data/tshepo_fulldata_final.rdata")
data$id = data$bid
data$id = gsub("^B(0+)|-", "", data$id)

data$rxnrti = substituteSKE(data$rxnrti, what=c(1,2,3), towhat=c("ZDV/3TC", "ZDV/ddI", "d4T/3TC"))
data$rxnnrti = substituteSKE(data$rxnnrti, what=c(1,2), towhat=c("NVP", "EFV"))
data$rxadh = substituteSKE(data$rxadh, what=c(1,2), towhat=c("SOC", "Com-DOT"))
data$alltreat = paste(data$rxnrti, data$rxnnrti, data$rxadh, sep="/")
data$treat = paste(data$rxnrti, data$rxnnrti, sep="/")

label(data$id) = "ID"
label(data$adhstrat) = "Observed Therapy"
label(data$age) = "Age"
label(data$basebmi) = "Baseline BMI"
label(data$cd4) = "CD4 Count"
label(data$dthcscd) = "Cause of Death"
label(data$height) = "Height"
label(data$hemglbn) = "Hemoglobin"
label(data$karnof) = "Karnofsky Score"
label(data$lastdate) = "Last Date Known Alive"
label(data$lastvsdt) = "Date of Last Visit"
label(data$lstdosdt) = "Date of Last Dose"
label(data$offstrs) = "End of Study Status"
label(data$randdt) = "Randomization Date"
label(data$sex) = "Sex"
label(data$stagecln) = "WHO clinical stage (1-4)"
label(data$stageprf) = "Karnofsky score (0-100)"
label(data$visitdt) = "First Visit Date"
label(data$vl_0) = "Baseline Viral Load"
label(data$weight) = "Weight"
label(data$mod_ane) = "Moderate Anemia (Yes No)"
label(data$wasting) = "HIV wasting diagnosis YES or NO"
label(data$tb_any) = "TB Diagnosis"
label(data$zoster_new) = "Incident Herpes Zoster (Shingles) Cases?"
label(data$ks_new) = "Incident Kaposi Sarcoma Cases?"
label(data$sev_ane) = "Severe Anemia Yes or No"
label(data$endptdt) = "Date of End Point"
label(data$deathdate) = "Death Date"
label(data$deathwk) = "Death Date (in Terms of Number of Weeks on Study)"
label(data$vfrescens) = "Virologic Failure with Resistance (Genotypic Drug Resistance)"
label(data$tmtdt_n) = "Date of Treatment Modifying Toxicity"
label(data$deathcens) = "Died"

label(data$cd4nnrti) = "median CD4 by NNRTI (those randomized to efavirenz versus nevirapine)??"
label(data$cd4strat) = "median CD4 by CD4 strata (upper (CD4 strata 201-350) versus lower CD4 strata (< 201)"
label(data$hemstrat) = "hemoglobin strata (not sure why this is here as we did not stratify by hemoglobin, but in analyses as low hemoglobin (anemia) is a predictor of death, the data may be stratified by hemoglobin cutoffs (i.e. < 8 grams/dL = severfe anemia, 8-10 moderate???) "
label(data$latedth) = "not sure???  (late death???) we did have one patient who died post coming off study?  Hwanhee will need to help us here??"
label(data$log_vl_0) = "baseline log viral load"
label(data$offstcd) = "CD4 cell count value off study??? not sure why we have this as once patients completed their 39th visit (7390) they completed 3 years of follow-up; hwanhee please assist??"
label(data$offstdt1) = "off study date (the date they completed 3 years of study follow-up)"
label(data$pcpdx) = "pneumocystis jiroveci (formerly carinii) pneumonia diagnosis?  yes or no??"
label(data$rxadh) = "treatment regimen by adherence strategy??"
label(data$rxnnrti) = "NNRTI assignment (i.e. were they randomized to nevirapine (NVP) or efavirenz (EFV)??"
label(data$rxnrti) = "NRTI assignment (zidovudine (ZDV) plus lamuvidine (3TC); zidovudine (ZDV) plus didanosine (ddI), OR stavudine (d4T) plus lamuvudine (3TC))??"
label(data$rxonly) = "???"
label(data$rxstrdt1) = "???"
label(data$rxstrdt2) = "???"
label(data$rxstrdt3) = "???"
label(data$stopdt1) = "???"
label(data$stopdt2) = "???"
label(data$stopdt3) = "???"
label(data$neuro) = "incident peripheral neuropathy (does this include all grades and does this exclude those few patients who had grade 1 peripheral neuropathy at baseline??)"
label(data$drugrsn) = "adverse drug reaction to one or more of the antiretroviral medications (or to any medication; including antibiotics (bactrim, penicillin, etc.)?  YES or NO"
label(data$drugrade) = "grade of adverse drug reaction?"
label(data$medcensdt) = "medication censor date?? is this for as-treated not ITT anlayses--please clarify??"
label(data$medcens) = "???"
label(data$medcenswk) = "???"
#label(data$deathcens) = "date death censored??"
label(data$vfresdt) = "???"
label(data$vfreswk) = "???"
label(data$deathdt_n) = "???"
label(data$deathwk_n) = "???"
label(data$tmtwk_n) = "???"
label(data$vfresdt_n) = "???"
label(data$vfreswk_n) = "???"

#tmp = latex(file=paste('gentex/', "describedata", ".tex", sep=""), describe(data))

data = subset(data, select=c("age", "basebmi", "cd4", "dthcscd", "height", "hemglbn", "karnof", "lastdate", "lastvsdt", "latedth", "lstdosdt", "offstdt1", "offstrs", "pcpdx", "randdt", "rxadh", "rxnnrti", "rxnrti", "rxonly", "rxstrdt1", "rxstrdt2", "rxstrdt3", "sex", "stagecln", "stopdt1", "stopdt2", "stopdt3", "visitdt", "vl_0", "weight", "mod_ane", "wasting", "tb_any", "zoster_new", "ks_new", "neuro", "sev_ane", "drugrsn", "drugrade", "medcens", "medcenswk", "endptdt", "deathdate", "deathwk", "vfresdt", "vfreswk", "deathdt_n", "deathwk_n", "tmtdt_n", "tmtwk_n", "vfresdt_n", "vfreswk_n", "id", "bid", "deathcens"))

data$sex = factor(substituteSKE(data$sex, what=c(1, 2), towhat=c("male", "female")))
label(data$sex) = "Sex"
data$followuptime = as.numeric(chron(dates. = as.character(data$lastdate), format = c(dates = "m/d/y")) - chron(dates. = as.character(data$randdt), format = c(dates = "m/d/y")))
data$followuptimeweeks = round(data$followuptime/7, 1)
label(data$followuptimeweeks) = "Follow-up time (weeks)"

data$timetolastvisitweeks = as.numeric(chron(dates. = as.character(data$lastvsdt), format = c(dates = "m/d/y")) - chron(dates. = as.character(data$randdt), format = c(dates = "m/d/y")))
data$timetolastvisitweeks = round(data$timetolastvisitweeks/7, 1)
label(data$timetolastvisitweeks) = "Time to last visit (weeks)"

data$timetolastdoseweeks = as.numeric(chron(dates. = as.character(data$lastdate), format = c(dates = "m/d/y")) - chron(dates. = as.character(data$randdt), format = c(dates = "m/d/y")))
data$timetolastdoseweeks = round(data$timetolastdoseweeks/7, 1)
label(data$timetolastdoseweeks) = "Time to last dose (weeks)"

data$timetodeathorcens = as.numeric(chron(dates. = as.character(data$deathdate), format = c(dates = "m/d/y")) - chron(dates. = as.character(data$randdt), format = c(dates = "m/d/y")))
data$timetodeathorcens[data$deathcens==0] = data$followuptime[data$deathcens==0]

lacticAcidosisSubj = c("B002233-7", "B002305-8", "B000896-3", "B001529-7", "B000872-9", "B001423-7", "B002120-0", "B002219-3", "B002306-4", "B001617-4", "B001103-5", "B002202-6", "B001542-5", "B001897-2", "B001054-2", "B001105-7", "B001378-8", "B001624-1", "B001847-7")
lacticAcidosisSubjIDs = gsub("^B(0+)|-", "", lacticAcidosisSubj)

### not cases according to Bill's email from 02/01/2010
notLacticAcidosisIDs =  gsub("^B(0+)|-", "", c("B001532-0", "B001469-3"))

pancrSubj = pancdata$bid
pancrSubjIDs = gsub("^B(0+)|-", "", pancrSubj)

skintSubj = skintdata$bid
skintSubjIDs = gsub("^B(0+)|-", "", skintSubj)


#outcomes = c("lactAcidosis", "clinPancreatitis")

data$lactacid = 0
data$clinpanc = 0
data$skintoxi = 0
data$lactacid[data$bid %in% lacticAcidosisSubj] = 1
data$clinpanc[data$bid %in% pancrSubj] = 1
data$skintoxi[data$bid %in% skintSubj] = 1
data = merge(data, data.frame(bid=pancdata$bid, pancdt = pancdata$date), by="bid", all.x=TRUE)
data = merge(data, data.frame(bid=skintdata$bid, skintdt = skintdata$date), by="bid", all.x=TRUE)

save(data, file="../data/data.rda", compress=TRUE)

######################################################################################3
### download data for lactic acidosis, the one bill used in hist paper
#######################################################################################
dateOutFormat = "day-month-year"
spssdata = read.csv("../data/lactic_acidosis_spss_jaids_ms_final_dataset_05_june_2007_using_30_nov_2005_as_data_cutoff_cww.csv", as.is=TRUE)
names(spssdata) = tolower(names(spssdata))
names(spssdata) = gsub("^bid$", "id", names(spssdata))
spssdata$datelact = chron(dates.=as.character(spssdata$datelact), out.format = dateOutFormat)
save(spssdata, file="../data/spssdata.rda", compress=TRUE)

lacticdata1 = read.csv("../data/all_lact FINAL FILE FROM ANN 14 NOV 2006from30NOV2005.csv", as.is=TRUE)
names(lacticdata1) = tolower(names(lacticdata1))
lacticdata1$id = lacticdata1$bid
lacticdata1$id = gsub("^B(0+)|-", "", lacticdata1$id)
save(lacticdata1, file="../data/lacticdata1.rda", compress=TRUE)

lacticdata2 = read.csv("../data/all_lact FINAL FILE FROM ANN 14 NOV 2006ALLLACTATES.csv", as.is=TRUE)
names(lacticdata2) = tolower(names(lacticdata2))
lacticdata2$id = lacticdata2$bid
lacticdata2$id = gsub("^B(0+)|-", "", lacticdata2$id)
save(lacticdata2, file="../data/lacticdata2.rda", compress=TRUE)

lacticdata3 = read.csv("../data/all_lact30NOV(1) AMT.csv", as.is=TRUE)
names(lacticdata3) = tolower(names(lacticdata3))
lacticdata3$id = lacticdata3$bid
lacticdata3$id = gsub("^B(0+)|-", "", lacticdata3$id)
save(lacticdata3, file="../data/lacticdata3.rda", compress=TRUE)

tmp = merge(lacticdata1, lacticdata3, by=c("id", "bid", "specdt", "lactate"), all=TRUE)
lacdata = merge(tmp, lacticdata2, by=c("id", "bid", "specdt", "lactate", "grade"), all=TRUE)

spssnotinexcel = spssdata[spssdata$id %in% setdiff(spssdata$id, lacdata$id),]
lacdata$specdt = chron(dates.=as.character(lacdata$specdt), format = "d-month-y", out.format = dateOutFormat)
lacdata = rbind(lacdata, data.frame(id =  spssnotinexcel$id, bid=NA, specdt=spssnotinexcel$datelact, lactate = spssnotinexcel$lacresult, grade=NA))

###---------------------------------find time to 4.4 or more of lactate
lacdata = lacdata[order(lacdata$id, lacdata$specdt), ]
lacdata$num = 1:nrow(lacdata)
tmpld = lacdata[lacdata$lactate>=4.4,]
tmp = tapply(tmpld$num, tmpld$id, function(x) as.character(min(x)))
tmpd = subset(lacdata, num %in% tmp, select=c("id", "specdt"))
names(tmpd) = c("id", "eventdt")
lacdata = merge(lacdata, tmpd, by="id", all.x=TRUE)
lacdata$eventdt = chron(lacdata$eventdt, format="day-month-year")
lacdata$num=NULL
lacdata$screened = TRUE

save(lacdata, file="../data/lacdata.rda", compress=TRUE)


#################################################################################
################ data per id
#################################################################################
perid = data.frame(id = data$id,
                   #tx = data$rxnrti,
                   tx = data$rxonly,
                   txnucl = data$rxnrti,
                   txnonnucl = data$rxnnrti,
                   txadher = data$rxadh,
                   age = data$age,
                   bmi = data$basebmi,
                   sex = data$sex,
                   cd4 = data$cd4,
                   viralload = data$vl_0,
                   hemglbn = data$hemglbn,
                   karnof = data$karnof,
                   clinpanc = data$clinpanc,
                   pancdt = data$pancdt,
                   skintoxi = data$skintoxi,
                   skintdt = data$skintdt,
                   firsttvistitdt = chron(dates. = as.character(data$visitdt), format = c(dates = "m/d/y"), out.format=dateOutFormat),
                   deathdt = chron(dates. = as.character(data$deathdt), format = c(dates = "m/d/y"), out.format=dateOutFormat),
                   timetodeathorcens = data$timetodeathorcens,
                   deathcens = data$deathcens,
                   lastalivedt = chron(dates. = as.character(data$lastdate), format = c(dates = "m/d/y"), out.format=dateOutFormat))
perid = merge(perid, unique(lacdata[,c("id", "eventdt", "screened")]), by=c("id"), all.x=TRUE)
perid$timetodeath = as.numeric(perid$deathdt - perid$firsttvistitdt)
perid$timetolastalive = as.numeric(perid$lastalivedt - perid$firsttvistitdt)
perid$timetoevent = as.numeric(perid$eventdt - perid$firsttvistitdt)
perid$timetoevent[is.na(perid$eventdt)] = perid$timetolastalive[is.na(perid$eventdt)]
perid$timetopanc = as.numeric(perid$pancdt - perid$firsttvistitdt)
perid$timetopanc[is.na(perid$pancdt)] = perid$timetolastalive[is.na(perid$pancdt)]
perid$timetoskint = as.numeric(perid$skintdt - perid$firsttvistitdt)
perid$timetoskint[is.na(perid$skintdt)] = perid$timetolastalive[is.na(perid$skintdt)]
perid$coxevent = 0
perid$coxevent[!is.na(perid$eventdt)] = 1
perid$coxevent[perid$id %in% notLacticAcidosisIDs] = 0
perid$poevent = 0
perid$poevent[!is.na(perid$screened)] = 1
perid$poevent[!is.na(perid$eventdt)] = 2
perid$poevent[perid$id %in% notLacticAcidosisIDs] = 0
perid$screened[is.na(perid$screened)] = FALSE
perid = unique(perid)
perid$tx1 = ifelse(perid$tx==1, 1, 0)
perid$tx2 = ifelse(perid$tx==2, 1, 0)
perid$tx3 = ifelse(perid$tx==3, 1, 0)
perid$tx4 = ifelse(perid$tx==4, 1, 0)
perid$tx5 = ifelse(perid$tx==5, 1, 0)
perid$sexf = ifelse(perid$sex=="female", 1, 0)
perid$logviralload10 = log(perid$viralload)

#################################################################################
################ data for cox regression
#################################################################################
coxdata = perid
save(coxdata, file="../data/coxdata.rda", compress=TRUE)

#################################################################################
################ data for proportional odds regression
#################################################################################
podata = perid
save(podata, file="../data/podata.rda", compress=TRUE)

#################################################################################
################ data for proportional odds regression
#################################################################################
pancrdata = perid
save(pancrdata, file="../data/pancrdata.rda", compress=TRUE)
#################################################################################
################ data
#################################################################################
save(perid, file="../data/perid.rda", compress=TRUE)

##############################################################################################
##############################################################################################
##############################################################################################
##############################################################################################
###  imputation source code
##############################################################################################
##############################################################################################
##############################################################################################
##############################################################################################

sr = function(){
  source("usefulFun.R")
}

firstTimeReachVal = function(time, allValue, critValue){
  min(time[value >= critValue])
}

firstTimeReachVal = function(data, timeVarName, allValueVarName, critValue){
  min(data[,timeVarName][data[,allValueVarName] >= critValue])
}

findFirstDateOfBadLactate = function(lacdata){
  lacdata = lacdata[order(lacdata$id, lacdata$specdt), ]
  lacdata$num = 1:nrow(lacdata)
  tmp1 = tapply(lacdata$num, lacdata$id, min)
  tmp2 = tapply(lacdata$num, lacdata$id, max)
  indexd = merge(data.frame(id=names(tmp1), start=tmp1), data.frame(id=names(tmp2), stop=tmp2), by="id")
  indexd$firstdate = NA
  for (i in 1:nrow(indexd)){
    indexd$firstdate[i] = as.character(firstTimeReachVal(lacdata[indexd$start[i]:indexd$stop[i],], "specdt", "lactate", 4.4))
  }
  indexd$firstdate[indexd$firstdate==Inf] = NA
  lacdata = merge(lacdata, indexd[, c("id", "firstdate")], by="id")
  #lacdata$firstdate = chron(lacdata$firstdate, format="day-month-year")
  lacdata$firstdate
}

findFirstDateOfBadLactate = function(lacdata){
  lacdata = lacdata[order(lacdata$id, lacdata$specdt), ]
  lacdata$num = 1:nrow(lacdata)
  lacdata$minnum = tapply(lacdata$num, lacdata$id, min)[as.character(lacdata$id)]
  lacdata$firstdate = lacdata$specdt
  lacdata$firstdate[lacdata$minnum != lacdata$num] = NA
  lacdata$firstdate
}


### ------------------------------------ multiple imputation by for PO LR model
multipleImpForLRM = function(data, n=5){
  set.seed(1)
  ### n - number of times to run the loop of imputation
  data$imbmi = data$bmi
  data$log10viralload = log(data$viralload, base=10)
  data$imhemglbn = data$hemglbn
  impdata = list()
  ### --------------------------------------- impute viral load
  data$log10viralload[is.na(data$viralload)] = median(log(data$viralload, base=10), na.rm=TRUE)
  ### --------------------------------------- impute bmi
  data$imbmi[is.na(data$bmi)] = median(data$bmi, na.rm=TRUE)
  for (i in 1:n) {
    ### --------------------------------------- predict hemoglobin
    hemglFit = lm(hemglbn ~ age + sex + cd4 + log10viralload + imbmi + poevent + txnucl + txnonnucl + txadher, data=data)
    predictloghemglbn = predict(hemglFit, newdata = data)
    hemglResid = sample(residuals(hemglFit), nrow(data), replace=TRUE)
    ### --------------------------------------- impute hemoglobin using predicted hemoglobin
    data$imhemglbn = ifelse(is.na(data$hemglbn), (predictloghemglbn + hemglResid)*(predictloghemglbn + hemglResid > 0), data$hemglbn)

    ### --------------------------------------- predict bmi
    bmiFit = lm(bmi ~ age + sex + cd4 + log10viralload + imhemglbn + poevent + txnucl + txnonnucl + txadher, data=data)
    predictbmi = predict(bmiFit, newdata = data)
    bmiResid = sample(residuals(bmiFit), nrow(data), replace=TRUE)
    ### --------------------------------------- impute bmi using predicted bmi
    data$imbmi = ifelse(is.na(data$bmi), (predictbmi + bmiResid)*(predictbmi + bmiResid > 0), data$bmi)

    ### --------------------------------------- predict viral load
    vloadFit = lm(log10viralload ~ age + sex + cd4 + imbmi + imhemglbn + poevent + txnucl + txnonnucl + txadher, data=data)
    predictlog10vload = predict(vloadFit, newdata = data)
    vloadResid = sample(residuals(vloadFit), nrow(data), replace=TRUE)
    ### --------------------------------------- impute viral load using predicted viral load
    data$log10viralload = ifelse(is.na(data$viralload), (predictlog10vload + vloadResid) * (predictlog10vload + vloadResid > 0), data$log10viralload)

    impdata[[i]] = data
  }
  impdata
}

multipleImpForCPH = function(data, timeVar, eventVar, n=5){
  set.seed(1)
  ### n - number of times to run the loop of imputation
  data$imbmi = data$bmi
  data$log10viralload = log(data$viralload, base=10)
  data$imhemglbn = data$hemglbn
  impdata = list()
  ### --------------------------------------- impute viral load
  data$log10viralload[is.na(data$viralload)] = median(log(data$viralload, base=10), na.rm=TRUE)
  ### --------------------------------------- impute bmi
  data$imbmi[is.na(data$bmi)] = median(data$bmi, na.rm=TRUE)
  for (i in 1:n) {
    ### --------------------------------------- predict hemoglobin
    hemglFit = lm(as.formula(paste("hemglbn ~ age + sex + cd4 + log10viralload + imbmi + txnucl + txnonnucl + txadher", eventVar, timeVar, sep=" + ")), data=data)
    predictloghemglbn = predict(hemglFit, newdata = data)
    hemglResid = sample(residuals(hemglFit), nrow(data), replace=TRUE)
    ### --------------------------------------- impute hemoglobin using predicted hemoglobin
    data$imhemglbn = ifelse(is.na(data$hemglbn), (predictloghemglbn + hemglResid)*(predictloghemglbn + hemglResid > 0), data$hemglbn)

    ### --------------------------------------- predict bmi
    bmiFit = lm(as.formula(paste("bmi ~ age + sex + cd4 + log10viralload + imhemglbn + txnucl + txnonnucl + txadher", eventVar, timeVar, sep=" + ")), data=data)
    predictbmi = predict(bmiFit, newdata = data)
    bmiResid = sample(residuals(bmiFit), nrow(data), replace=TRUE)
    ### --------------------------------------- impute bmi using predicted bmi
    data$imbmi = ifelse(is.na(data$bmi), (predictbmi + bmiResid)*(predictbmi + bmiResid > 0), data$bmi)

    ### --------------------------------------- predict viral load
    vloadFit = lm(as.formula(paste("log10viralload ~ age + sex + cd4 + imbmi + imhemglbn + txnucl + txnonnucl + txadher", eventVar, timeVar, sep=" + ")), data=data)
    predictlog10vload = predict(vloadFit, newdata = data)
    vloadResid = sample(residuals(vloadFit), nrow(data), replace=TRUE)
    ### --------------------------------------- impute viral load using predicted viral load
    data$log10viralload = ifelse(is.na(data$viralload), (predictlog10vload + vloadResid) * (predictlog10vload + vloadResid > 0), data$log10viralload)
    data$bmi3 = NA
    data$bmi3 = ifelse(data$imbmi<18.5, "<18.5", data$bmi3)
    data$bmi3 = ifelse(data$imbmi>=18.5 & data$imbmi<=25, "18.5-25", data$bmi3)
    data$bmi3 = ifelse(data$imbmi>25, ">25", data$bmi3)
    data$bmi2 = NA
    data$bmi2 = ifelse(data$imbmi<=25, "<=25", data$bmi2)
    data$bmi2 = ifelse(data$imbmi>25, ">25", data$bmi2)

    impdata[[i]] = data
  }
  impdata
}

multipleImpForCPH_short = function(data, timeVar, eventVar, n=5){
  set.seed(1)
  ### n - number of times to run the loop of imputation
  data$imbmi = data$bmi
  data$imhemglbn = data$hemglbn
  impdata = list()
  ### --------------------------------------- impute bmi
  data$imbmi[is.na(data$bmi)] = median(data$bmi, na.rm=TRUE)
  for (i in 1:n) {
    ### --------------------------------------- predict hemoglobin
    hemglFit = lm(as.formula(paste("hemglbn ~ age + sex + cd4 + imbmi", eventVar, timeVar, sep=" + ")), data=data)
    predictloghemglbn = predict(hemglFit, newdata = data)
    hemglResid = sample(residuals(hemglFit), nrow(data), replace=TRUE)
    ### --------------------------------------- impute hemoglobin using predicted hemoglobin
    data$imhemglbn = ifelse(is.na(data$hemglbn), (predictloghemglbn + hemglResid)*(predictloghemglbn + hemglResid > 0), data$hemglbn)
    ### --------------------------------------- predict bmi
    bmiFit = lm(as.formula(paste("bmi ~ age + sex + cd4 + imhemglbn", eventVar, timeVar, sep=" + ")), data=data)
    predictbmi = predict(bmiFit, newdata = data)
    bmiResid = sample(residuals(bmiFit), nrow(data), replace=TRUE)
    ### --------------------------------------- impute bmi using predicted bmi
    data$imbmi = ifelse(is.na(data$bmi), (predictbmi + bmiResid)*(predictbmi + bmiResid > 0), data$bmi)
    data$imbmi3 = NA
    data$imbmi3 = ifelse(data$imbmi<18.5, "<18.5", data$imbmi3)
    data$imbmi3 = ifelse(data$imbmi>=18.5 & data$imbmi<=25, "18.5-25", data$imbmi3)
    data$imbmi3 = ifelse(data$imbmi>25, ">25", data$imbmi3)
    data$imbmi2 = NA
    data$imbmi2 = ifelse(data$imbmi<=25, "<=25", data$imbmi2)
    data$imbmi2 = ifelse(data$imbmi>25, ">25", data$imbmi2)

    impdata[[i]] = data
    #models[[i]] = lrm(poevent ~ age + bmi + cd4 + log10viralload + hemglbn + sexf + txnucl + txnonnucl + txadher, data = data)
  }
  impdata
}

runPOLRModels = function(impdata, formula){
### PROPORTIONAL ODDS MODELS
  models = list()
  for (i in 1:length(impdata)){
    models[[i]] = lrm(formula, data = impdata[[i]])
  }
  models
}

runCPHModels = function(impdata, formula){
### SURVIVAL MODELS
  models = list()
  for (i in 1:length(impdata)){
    models[[i]] = cph(formula, data = impdata[[i]])
  }
  models
}

plotOddsRatios = function(oddsMatrix, xRange=NULL, yLabelLine=0, main=""){
  ####################################################
  ### matrix has to have three columns for OR and CIs
  ####################################################
  xscale = c(.5, 1, 2, 3, 4)
  roundn = 2
  height = .16
  mainY = nrow(oddsMatrix):1
  logOddsMatrix = log(oddsMatrix)
  par(mar=c(3, 16, 3, 2))
  yRange = c(0+1, nrow(logOddsMatrix)+1)
  if (is.null(xRange)){
    xRange = range(logOddsMatrix)+c(-.6, 1.2)
  }
  plot(x=0, y=0, type="n", ylim = yRange, xlim=xRange, axes=FALSE, ann=FALSE, main=main)
  lines(x=rep(log(1), 2), y=yRange+c(-.5, -.5), lty=3, col=gray(.5))
  rect(xleft=logOddsMatrix[,2], ybottom=mainY, xright=logOddsMatrix[,3], ytop=mainY + height, density = NULL, col="gray", border="gray")
  points(x=logOddsMatrix[,1], y=mainY+height/2, pch=18, cex=1.3)
  mtext(text=rownames(logOddsMatrix), at=mainY, side=2, line=yLabelLine, las=2, cex=.9)
  mtext(text=xscale, at=log(xscale), side=3, line=-1, las=1, cex=.7)
  mtext(text=main, at=log(1), side=3, line= 0, las=1, cex=1)
  text(x=max(logOddsMatrix[,3]), y=mainY, labels=paste(round(oddsMatrix[,1], roundn), " (",round(oddsMatrix[,2], roundn),", ", round(oddsMatrix[,3], roundn), ")", sep=""), pos=4, cex=.7)
  par(mar=c(5, 4, 4, 2))
}

retrStuffFromSummaryLrmOrCphForAllVars = function(modelsAndData, ...){
  ### LIMITATION! WILL NOT CALCULATE THE RIGHT PVALUES FOR NON-LINEAR !
  ####################################################
  ### get model results from lrm or cph summary for ALL variables
  ####################################################
  models = modelsAndData[["models"]]
  data = modelsAndData[["data"]]
  dd <<- datadist(data)
  options(datadist="dd")
  n = length(models)
  nvars = dim(summary(models[[1]]))[1]/2
  vars = rownames(summary(models[[1]], ...))[seq(1, 2*nvars, 2)]
  variances = matrix(NA, nrow=nvars, ncol=n)
  effects = matrix(NA, nrow=nvars, ncol=n)
  stuff =  matrix(NA, nrow=nvars, ncol=4)
  for (i in 1:n){
    summOut = summary(models[[i]], ...)
    #summOut = summary(models[[i]])
    effects[, i] = summOut[seq(1, 2*nvars, 2), c("Effect")]
    variances[, i] = summOut[seq(1, 2*nvars, 2), c("S.E.")]^2
  }
  meanEffect = apply(effects, 1, mean)
  ### bootstrap variance estimate
  if (n==1){
    totalVariance = variances
  }else{
    totalVariance = apply(variances, 1, mean) + (n+1) * apply(effects, 1, var) / n
  }
  stuff[,1:3] = cbind(exp(meanEffect), exp(meanEffect - 1.96*sqrt(totalVariance)), exp(meanEffect + 1.96*sqrt(totalVariance)))
  ### Z, Wald
  z = meanEffect/sqrt(totalVariance)
  stuff[, 4] = 1-pchisq(z^2, 1)
  ### add difference
  stuff = cbind(matrix(summOut[seq(1, 2*nvars, 2), c("Low", "High")], ncol=2), stuff)
  rownames(stuff) = vars
  colnames(stuff) = c("Low", "High", "Effect", "Lower 0.95", "Upper 0.95", "P")
  stuff
}

retrieveStuffFromPOLRWithoutSummary = function(modelsAndData){
  ####################################################
  ### get model results WITHOUT summary for ALL variables
  ####################################################
  models = modelsAndData[["models"]]
  data = modelsAndData[["data"]]
  n = length(models)
  mainvar = rownames(anova(models[[1]]))
  mainvar = mainvar[1:(length(mainvar)-1)]
  variances = matrix(NA, nrow=n, ncol=length(mainvar))
  colnames(variances) = mainvar
  effects<-matrix(NA, nrow=n, ncol=length(mainvar))
  colnames(effects) = mainvar
  stuff = list()
  for (i in 1:n){
    m = models[[i]]
    effects[i, mainvar] = m$coef[mainvar]
    variances[i, mainvar] = diag(m$var[mainvar, mainvar])
  }
  meanEffect = apply(effects, 2, mean)
  ### bootstrap variance estimate
  totalVariance = apply(variances, 2, mean) + (n+1)*apply(effects, 2, var)/n
  stuff[["estimates"]] = cbind(exp(meanEffect),
                               exp(meanEffect - 1.96*sqrt(totalVariance)),
                               exp(meanEffect + 1.96*sqrt(totalVariance)))
  colnames(stuff[["estimates"]]) = c("Effect", "Lower 0.95", "Upper 0.95")
  ### Z, Wald
  z = meanEffect/sqrt(totalVariance)
  stuff[["pvalues"]] = 1-pchisq(z^2, 1)
  stuff
}

missingVariables = function(data){
  tmp = sapply(data, function(x){sum(is.na(x))})
  tmp = tmp[tmp>0]
  res = matrix(tmp, nrow=length(tmp))
  rownames(res) = names(tmp)
  res = cbind(res, 100*tmp/nrow(data))
  colnames(res) = c(paste("# of NA out of", nrow(data)), "% of NA")
  res
}

unadjustedHRsIMPNew = function(coximpdata, timeVar="timetoevent", eventVar="coxevent", varListOfInput, varListOfLabel=NULL, unitVector, ...){
### runs cox models, for each variable and creates a table of unadjusted HRs
  hrMatrix = matrix(NA, nrow=1, ncol=5)
  survObjectFormula = paste("Surv(time=",timeVar,", event=",eventVar,") ~")
  counter = 1
  for (n in names(varListOfInput)){
    vlabel = if (is.null(varListOfLabel)) n else varListOfLabel[[n]]
    cphmod = runCPHModels(coximpdata, formula=as.formula(paste(survObjectFormula, n)))
    tmpMatrix = getStuffForOneVarFromFranksModelSummaryImp(cphmod, coximpdata[[1]], varListOfInput[[n]], n, vlabel, unit=unitVector[n])
    cat(n, "\n")
    tmpMatrix = cbind(tmpMatrix, matrix(rep(counter, nrow(tmpMatrix)), ncol=1))
    counter = counter + 1
    hrMatrix = rbind(hrMatrix, tmpMatrix)
  }
  hrMatrix[2:nrow(hrMatrix),]
}

pvaluesForNonLinear = function(coxdata, vars, timeVar="timetoevent", eventVar="coxevent", ...){
  dd = datadist(coxdata)
  options(datadist="dd")
  pvalues = rep(NA, length(vars))
  names(pvalues) = vars
  for (v in vars){
    survObj = Surv(time=coxdata[[timeVar]], event=coxdata[[eventVar]])
    cphmod = cph(as.formula(paste("survObj ~ rcs(", v,",3)")), data = coxdata, surv=TRUE)
    pvalues[v] = anova(cphmod)[" Nonlinear", "P"]
  }
  pvalues
}

pvaluesForNonLinearPOModel = function(data, vars, outcomeVar="poevent", ...){
  dd <<- datadist(data)
  options(datadist="dd")
  pvalues = rep(NA, length(vars))
  names(pvalues) = vars
  for (v in vars){
    pomod = lrm(as.formula(paste(outcomeVar, " ~ rcs(", v,",3)")), data = data)
    pvalues[v] = anova(pomod)[" Nonlinear", "P"]
  }
  pvalues
}

somethingAndCIsStr = function(somethingCilowCihigh, roundn=3){
  paste(round(somethingCilowCihigh[1], roundn), " (",round(somethingCilowCihigh[2], roundn), ",", round(somethingCilowCihigh[3], roundn), ")", sep="")
}

getNiceOutputORorHRfromSummary = function(models, data, varListOfInput, varListOfLabel=NULL, unitVector, varsToShowRange=c()){
  hrMatrix = matrix(NA, nrow=1, ncol=5)
  counter = 1
  for (n in names(varListOfInput)){
    vlabel = if (is.null(varListOfLabel)) n else varListOfLabel[[n]]
    tmpMatrix = getStuffForOneVarFromFranksModelSummaryImp(models, data, varListOfInput[[n]], n, vlabel, unit=unitVector[n], showRange = (n %in% varsToShowRange))
    cat(n, "\n")
    tmpMatrix = cbind(tmpMatrix, matrix(rep(counter, nrow(tmpMatrix)), ncol=1))
    counter = counter + 1
    hrMatrix = rbind(hrMatrix, tmpMatrix)
  }
  hrMatrix[2:nrow(hrMatrix),]
}

prepareHRorORtableForZebraTable = function(hrMatrix, roundn=3){
### hrMatrix is an output of either getNiceOutputORorHRfromSummary or unadjustedHRsIMPNew
  hrTable = data.frame(info=rownames(hrMatrix), hr=round(hrMatrix[,1], roundn), hrci=paste("(", round(hrMatrix[,2], roundn), ", ", round(hrMatrix[,3], roundn), ")", sep=""), pvalue=pvalueStr(hrMatrix[,4]), id=hrMatrix[,5], stringsAsFactors = FALSE)
  hrTable$hr[is.na(hrTable$hr)]=""
  hrTable$hrci[hrTable$hrci == "(NA, NA)"]=""
  hrTable$info = latexTextMode(hrTable$info)
  hrTable$info = gsub("\\.\\.\\.\\.", "~~~~", hrTable$info)
  hrTable$pvalue = latexTextMode(hrTable$pvalue)
  hrTable
}

getStuffForOneVarFromFranksModelSummaryImp = function(models, data, varInput, varName, varLabel=varName, unit="unit", showRange=FALSE){
  ### model - model (cph, lrm), already exponentiated
  ### varName - the name of the variable
  ### varInput - if varInput is one number(or string), then it is assumed
  ###    that the variable is categorical, and varInput is a reference value
  ###    if varInput is a matrix of Nx2 then it is assumed that the variable is numeric(or ordinal)
  ### units - units of numeric(or ordinal) variable
  ### showRange - if TRUE show for which values of the variable the HR or OR was calculated the default is FALSE
  # ###    but if variable has been transformed then showRange will be TRUE - impossible to implement (formula string doesn't carry over)
  ### RETURNS - a character matrix of ORs or HRs with CIs and Pvalue
  ### for example:
  ###
  ### "BMI"                    ""        ""        ""        ""
  ### "....18.5 <= BMI <= 25"     "1"        ""        ""        ""
  ### "....BMI < 18.5"          "1.1"     "0.4"     "3.1"     "0.8"
  ### "....BMI > 25"            "0.8"     "0.2"     "2.5"     "0.7"
  ###
  ### another example:
  ### "....CD4, continuous"        ""        ""        ""        ""
  ### "....(per 100 cells)"     "1.4"     "0.8"     "2.3"     "0.2"
  plural = function(val, str){
    ifelse(val==1, str, paste(str, "s", sep=""))
  }
  blankspstr = "...."
  if (length(varInput)==2) varInput = matrix(varInput, ncol=2)
  if (length(varInput)>1){
    if(class(varInput)!="matrix" | ncol(varInput)!=2) stop("varInput should be a matrix with 2 columns (see function comments)")
    ## if variable appears in a model formula as a non-linear term, set showRange to TRUE
    resMatr = matrix(NA, nrow = nrow(varInput) + 1, ncol=4)
    colnames(resMatr) = c("Effect", "Lower 0.95", "Upper 0.95", "P")
    varDiff = varInput[,2]-varInput[,1]
    if (showRange){
      rownames(resMatr) = c(varLabel, paste(blankspstr, "change from ", varInput[,1], " to ", varInput[,2], " ", plural(varDiff, unit), sep=""))
    }else{
      rownames(resMatr) = c(varLabel, paste(blankspstr, "per", " ", varDiff, " ", plural(varDiff, unit), sep=""))
    }
    for (i in 1:nrow(varInput)){
      args = list(list(models=models, data=data), varInput[i, ])
      names(args) = c("modelsAndData", varName)
      sumOutput = do.call("retrStuffFromSummaryLrmOrCphForAllVars", args)
      resMatr[1+i, c("Effect", "Lower 0.95", "Upper 0.95", "P")] = sumOutput[varName, c("Effect", "Lower 0.95", "Upper 0.95", "P")]
    }
  }else{
    args = list(list(models=models, data=data), varInput)
    names(args) = c("modelsAndData", varName)
    sumOutput = do.call("retrStuffFromSummaryLrmOrCphForAllVars", args)
    varNum = length(grep(varName, rownames(sumOutput)))
    resMatr = matrix(NA, nrow = varNum + 2, ncol=4)
    smallSumOut = sumOutput[grep(varName, rownames(sumOutput), value=TRUE),]
    if (is.null(dim(smallSumOut))){
      varRowNames = grep(varName, rownames(sumOutput), value=TRUE)
    }else{
      varRowNames = rownames(smallSumOut)
    }
    colnames(resMatr) = c("Effect", "Lower 0.95", "Upper 0.95", "P")
    rownames(resMatr) = c(varLabel, paste(blankspstr, varInput, sep=""), paste(blankspstr, varRowNames, sep=""))
    resMatr[2,1] = 1
    resMatr[paste(blankspstr, varRowNames, sep=""), c("Effect", "Lower 0.95", "Upper 0.95", "P")] = sumOutput[varRowNames, c("Effect", "Lower 0.95", "Upper 0.95", "P")]
    removeStr = paste(paste(varName, " - ", sep=""), "|", paste(":", varInput, sep=""), sep="")
    rownames(resMatr) = gsub(removeStr, "", rownames(resMatr))
  }
  resMatr
}