%#############################################################
%#
%#   Program: analysis-scripts-017.nw
%#	  
%#   Project: Correlates of suboptimal entry into early infant 
%#            diagnosis in rural north central Nigeria
%# 
%#   Co-I:    Muktar H. Aliyu*; Meridith Blevins; 
%#            Karen M. Megazzini; Julie Dunlap; 
%#            Carolyn M. Audet; Ibrahim S. Sodangi;  
%#            Usman I. Gebi; Bryan E. Shepherd; C. 
%#            William Wester; Sten H. Vermund; 
%#       * corresponding author, muktar.aliyu at vanderbilt.edu
%#	  
%#   Biostatistician/Programmer: Meridith Blevins, MS**
%#   ** contact programmer: meridith.blevins at vanderbilt.edu
%#	  
%#   Purpose: Reproducible research.
%#	      
%#   Notes: This '.nw' file is a Sweave file which allows
%#          for results from R statistical software to be
%#          nested within a LaTeX document.
%#          
%#   Created: 12 October 2012
%#
%#############################################################


% Hidden R code chunk --- reading in the data
<<echo = false, results = hide>>=
rm(list=ls()) # SHOULD BE FIRST LINE IN ALL R PROGRAMS - CLEARS NAMESPACE
library(Hmisc)

## READ IN DATA FOR ANALYSIS
setwd("~/Projects/Nigeria")
care <- read.csv(file="Data/Analysis_20131108.csv",na.strings = c("","NA"),stringsAsFactors=FALSE)
care <- cleanup.import(care, lowernames=TRUE,datevars=c("HAART_Start_Date","DOB","lastclin","lastpharm","lastlab","AC_PC_Enrl_Dt"), dateformat="%m/%d/%Y")
care$deathdt <- as.Date(care$deathdt,origin="1960-01-01")
care <- care[!duplicated(care$unique_id),]
care <- care[care$clinic!="Kuta",]
care <- droplevels(care)

eid <- read.csv(file="Data/EID_Data_All-mb.csv",na.strings = c("","NA"),stringsAsFactors=FALSE)
eid <- cleanup.import(eid, lowernames=TRUE,datevars=c("Date","DOB"), dateformat="%m/%d/%Y")
eid <- eid[!is.na(eid$site) & eid$site != "IBB",]
warning(paste(c("DROPPING 94 WOMEN WITH MULTIPLE BIRTHS FROM DATASET:",eid$unique.id[duplicated(eid$unique.id)]),collapse=" "))
## 94 women had 2 births; no women had >2 births
## duplicated(paste(eid$unique.id,eid$dob))
## KEEP MOST RECENT BIRTH FOR OUTCOMES ANALYSIS
eid <- eid[order(eid$date),]
eid <- eid[!duplicated(eid$unique.id,fromLast = TRUE),]
eid$EID <- "yes"
names(eid)[names(eid) %in% "dob"] <- "child_dob"
## REPLACE BAD IDs IDENTIFIED BY SODANGI (11/25)
eid_correct_id <- read.csv("Data/eid_correct_id.csv",na.strings = c("","NA"),stringsAsFactors=FALSE)
eid <- merge(eid,eid_correct_id,all.x=TRUE)
eid$unique.id <- ifelse(!is.na(eid$correct_id),eid$correct_id,eid$unique.id)
eid <- eid[!duplicated(eid$unique.id),]  ## INTRODUCED 3 MORE DUPLICATES

queries <- read.csv("Data/pregnancy_queries_20131108_IS.csv",na.strings = c("","NA"),stringsAsFactors=FALSE)
queries$Revise.sex..yes.no[is.na(queries$Revise.sex..yes.no.) & queries$gender=="Female" & queries$Pregnant..yes.no=="male"] <- "yes"
queries$Revise.sex..yes.no[is.na(queries$Revise.sex..yes.no.) & queries$gender=="Male" & queries$Pregnant..yes.no=="male"] <- "no"
queries$Revise.sex..yes.no[is.na(queries$Revise.sex..yes.no.) & queries$gender=="Male" & queries$Pregnant..yes.no=="no"] <- "yes"
queries$Revise.sex..yes.no[is.na(queries$Revise.sex..yes.no.) & queries$gender=="Male" & queries$Pregnant..yes.no=="yes"] <- "yes"
query1 <- queries[queries$Revise.sex..yes.no=="yes",]
query2 <- queries[queries$Revise.sex..yes.no!="yes" & queries$pregsource=="Never pregnant" & queries$Pregnant..yes.no.=="yes",]
query3 <- queries[queries$Revise.sex..yes.no!="yes" & queries$pregsource!="Never pregnant" & queries$Pregnant..yes.no.=="no",]
dim(query3) ## THIS NEVER OCCURRED

## TRUST EID ABOVE ALL ELSE (EVEN GENDER) -- FROM SODANGI
## ANALYSIS PREGNANT VARIABLES
## PREGNANT: BASED ON ART REGIMEN
## ENPREG: If AC_Crnt_Prgn="Yes" within 90 days before or 30 days 
## EVERPREG: If AC_Crnt_Prgn="Yes"
eidnocare <- eid[!(eid$unique.id %in% care$unique_id),]

analysis <- merge(care,eid,by.x="unique_id",by.y="unique.id",all.x=TRUE)
analysis$EID[is.na(analysis$EID)] <- "no"
#analysis <- analysis[analysis$gender %in% "Female",]
analysis$preg <- rep("no",nrow(analysis))
analysis$pregsource <- ifelse(!is.na(analysis$pregnant) & analysis$pregnant=="yes","ARV regimen","")
analysis$preg[!is.na(analysis$pregnant) & analysis$pregnant=="yes"] <- "yes"
analysis$pregsource[analysis$pregsource=="" & !is.na(analysis$everpreg) & analysis$everpreg=="yes"] <- "CAREware"
analysis$pregsource[analysis$pregsource=="ARV regimen" & !is.na(analysis$everpreg) & analysis$everpreg=="yes"] <- "ARV regimen and CAREware"
analysis$preg[!is.na(analysis$everpreg) & analysis$everpreg=="yes"] <- "yes"
analysis$preg[analysis$EID=="yes"] <- "yes"
analysis$pregsource[analysis$pregsource=="" & analysis$EID=="yes"] <- "EID"
analysis$pregsource[analysis$pregsource=="ARV regimen" & analysis$EID=="yes"] <- "EID and ARV regimen"
analysis$pregsource[analysis$pregsource=="CAREware" & analysis$EID=="yes"] <- "EID and CAREware"
analysis$pregsource[analysis$pregsource=="ARV regimen and CAREware" & analysis$EID=="yes"] <- "EID and ARV regimen and CAREware"
analysis$pregsource[analysis$preg=="no"] <- "Never pregnant"

## INCORPORATE QUERY RESULTS 
analysis$gender[analysis$gender=="Male" & analysis$unique_id %in% query1$unique_id] <- "Female"
analysis$gender[analysis$gender=="Female" & analysis$unique_id %in% query1$unique_id] <- "Male"
analysis$preg[analysis$unique_id %in% query2$unique_id] <- "yes"
analysis$pregsource[analysis$unique_id %in% query2$unique_id] <- "HIVCare"
analysis$preg[analysis$unique_id %in% query3$unique_id] <- "no"
analysis$pregsource[analysis$unique_id %in% query3$unique_id] <- "Never pregnant"
## KEEP ONLY ENROLLMENT DATES OCCURING ON OR BEFORE 12/31/2012
analysis <- analysis[analysis$ac_pc_enrl_dt < as.Date("01-01-2013",format="%m-%d-%Y"),]
## KEEP ONLY PREGNANT FEMALES
analysis <- analysis[analysis$preg=="yes" & analysis$gender=="Female",]
## IMPUTE GBRH FOR 58 WOMEN PER SODANGI 11/20 EMAIL -- FORMATTING ISSUE IN IDs HAD AUTOMATICALLY ASSIGNED INCORRECT SITE
analysis$clinic[!is.na(analysis$clinic) & analysis$clinic=="UMYMH" & analysis$site=="GBRH" & !is.na(analysis$site)] <- "GBRH"


# FUNCTIONS tabulate counts/proportions and median(q1,q3) overall and by groups
## PRINT ONE LINE FOR VARIABLES WITH TWO CATEGORIES OR ONLY ONE OF INTEREST 
## EXAMPLE (PRINT FEMALE (%) -- NO NEED TO PRINT MALE(%) AS WELL)
# var must be a variable that takes on 1 value and the rest are NA
# similarly, rowlabel must be length 1
getnp1line <- function(var,by,rowlabel="",indent=TRUE,inc.pval=FALSE,by.order=FALSE,formatper="myper0",rowper=FALSE,noper=FALSE){
  rowlabel <- ifelse(rowlabel=="",label(var),rowlabel)
  if(rowlabel==""){stop("rowlabel was not specified, either as an attribute of var, or in the getnp1line call.")}
  myper1 <- get(formatper)
  n <- table(var,by)
  d <- table(by)
  n1 <- table(var)
  if(!rowper){p <- sapply(100*n[1,]/d,myper1)}
  if(rowper){p <- sapply(100*n[1,]/n1,myper1)}
  pt <- sapply(100*n1/length(by),myper1)
  if(indent==TRUE){line <- paste(paste("\\\\hspace{.2in}",rowlabel),paste(paste(n[1,]," (",p,"\\\\%)",sep=""),collapse=" & "),paste(sum(n)," (",pt,"\\\\%)",sep=""),sep=" & ")}
  else{line <- paste(paste(rowlabel),paste(paste(n[1,]," (",p,"\\\\%)",sep=""),collapse=" & "),paste(sum(n)," (",pt,"\\\\%)",sep=""),sep=" & ")}
  if(indent & noper){line <- paste(paste("\\\\hspace{.2in}",rowlabel),paste(n[1,],collapse=" & "),sum(n),sep=" & ")}
  if(!indent & noper){line <- paste(paste(rowlabel),paste(n[1,],collapse=" & "),sum(n),sep=" & ")}
  if(inc.pval==TRUE){
    var1 <- ifelse(is.na(var),0,var)
    if(by.order==TRUE){pval <- myformat4(wilcox.test(as.numeric(by)~ var1)$p.value)}
    if(by.order==FALSE){pval <- myformat4(chisq.test(table(var1,by))$p.value)}
    line <- paste(line,pval,sep=" & ")
    }
  return(list(line=line))
}

## PRINT MED (IQR) FOR VAR WITHIN EACH BY COLUMN AND OVERALL
# options include 'by.order'' which will affect the hypothesis test
# note will include \tnote{note} by the label
# OUTPUT INCLUDES line and line.incmiss (if there is missing data)
getmedtex <- function(var,by,label="",note=note,inc.pval=TRUE,lower.label=TRUE,miss.label=TRUE,by.order=FALSE,digits=8,formatper="myper0",iqrsep=" -- "){
  label <- ifelse(label=="",label(var),label)
  if(label==""){stop("label was not specified, either as an attribute of var, or in the getmedtex call.")}
  myper1 <- get(formatper)
  if(!is.numeric(var)){stop("The variable is non-numeric")}
  s <- round(summary(var),digits=digits)
  text <- paste(s[3]," (",s[2],iqrsep,s[5],")",sep="")
  level <- sort(unique(by))
  sby <- matrix(nrow=length(level), ncol=length(s)); textby <- rep(NA,length(level))
  for(i in 1:length(level)){
      sby[i,1:5] <- round(summary(var[by==level[i]])[1:5],digits=digits)
      textby[i] <- ifelse(!is.na(sby[i,3]),paste(sby[i,3]," (",sby[i,2],iqrsep,sby[i,5],")",sep=""), "--")
  }
  if(inc.pval==TRUE){
      if(length(level)==1) {pval <- " "}
      else if(length(level)>1 & by.order==FALSE) {pval <- myformat4(kruskal.test(var ~ by)$p.value)}
      else if(length(level)>1 & by.order==TRUE)  {pval <- myformat4(cor.test(var, as.numeric(by),alternative="two.sided",method="spearman")$p.value)}
      
      if(missing(note)){line1 <- paste(label,paste(textby,collapse=' & '),text,pval,sep=' & ')
      } else{line1 <- paste(paste(label,"\\\\tnote{",note,"}",sep=""),paste(textby,collapse=' & '),text,pval,sep=' & ')
      }
  } else{
      if(missing(note)){line1 <- paste(label,paste(textby,collapse=' & '),text,sep=' & ')
      } else{line1 <- paste(paste(label,"\\\\tnote{",note,"}",sep=""),paste(textby,collapse=' & '),text,sep=' & ')
      }
  }
  if(sum(is.na(var))>0){
    no <- table(is.na(var),by)[2,]
    nop <- sapply(100*no/table(by),myper1)
    nopt <- sapply(100*sum(no)/sum(table(by)),myper1)
    if(lower.label==TRUE){line2 <- paste(paste("\\\\hspace{.1in}Missing",ifelse(miss.label,paste(" ",tolower(label),sep=""),""),", n(\\\\%)",sep=""),paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),paste(sum(no)," (",nopt,"\\\\%)",sep=""),sep=" & ")}
      else{line2 <- paste(paste("\\\\hspace{.1in}Missing",ifelse(miss.label,paste(" ",label,sep=""),""),", n(\\\\%)",sep=""),paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),paste(sum(no)," (",nopt,"\\\\%)",sep=""),sep=" & ")}
    line3 <- paste(line1,line2,sep=" \\\\\\\\ ")
    return(list(line=line1,line.incmiss=line3))
  }
  else{return(list(line=line1))}
}
   
## PRINT N (%) FOR VAR WITHIN EACH BY COLUMN AND OVERALL
# options include 'by.order' and 'var.order' which will affect the hypothesis test
# note will include \tnote{note} by the label
# OUTPUT INCLUDES line and line.incmiss (if there is missing data)
getnptex <- function(var,by,label="",rowlabel=rowlabel,note=note,inc.pval=TRUE,var.order=FALSE,by.order=FALSE,formatper="myper0",rowper=FALSE){
  label <- ifelse(label=="",label(var),label)
  if(label==""){stop("label was not specified, either as an attribute of var, or in the getnptex call.")}
  myper1 <- get(formatper)
  if(missing(rowlabel)){
    if(is.factor(var)){rowlabel <- levels(var)}
    if(is.character(var)){rowlabel <- names(table(var))}
  }
  if(length(rowlabel) != sum(!(unique(var)%in% NA))){warning("Function GETNPTEX: number of unique values does not match length of rowlabel")}
  n <- table(var,by)
  n1 <- rowSums(n)
  d <- colSums(table(var,by)[1:nrow(n),])
  p <- matrix(nrow=nrow(n), ncol=ncol(n))
  pt <- rep(NA,nrow(n))
  lines <- rep(NA,nrow(n))
  for(i in 1:nrow(n)){
    p[i,which(d==0)] <- 0
    if(!rowper){p[i,which(d!=0)] <- sapply(100*n[i,which(d!=0)]/d[which(d!=0)],myper1)}
    if(rowper){p[i,which(d!=0)] <- sapply(100*n[i,which(d!=0)]/n1[i],myper1)}
    pt[i] <- sapply(100*table(var)[i]/sum(table(var)[1:nrow(n)]),myper1)
    lines[i] <- paste(paste("\\\\hspace{.2in}",rowlabel[i]),paste(paste(n[i,]," (",p[i,],"\\\\%)",sep=""),collapse=" & "),paste(sum(n[i,])," (",pt[i],"\\\\%)",sep=""),sep=" & ")
  }
# ADD PVALUE TO COLUMN HEADING IF DESIRED
  if(inc.pval==TRUE){
    if(by.order==TRUE){pval <- myformat4(kruskal.test(as.numeric(by), var)$p.value)}
    if(var.order==TRUE){pval <- myformat4(kruskal.test(var, by)$p.value)}
    if(by.order==TRUE & var.order==TRUE){pval <- myformat4(cor.test(var, by, alternative="two.sided", method="spearman")$p.value)} # ADDED 12/30/2010
    if(by.order==FALSE & var.order==FALSE){pval <- myformat4(chisq.test(n[n1!=0,d!=0])$p.value)}
    spaces <- rep(" ",1+length(table(by)))
    line1 <- paste(paste(label,", n(\\\\%)",sep=""),paste(spaces,collapse=" & "),pval,sep=" & ")
  } else{
    line1 <- paste(label,", n(\\\\%)",sep="")
  }
  all.lines <- paste(line1,paste(lines[which(n1!=0)],collapse=" \\\\\\\\ "),sep=" \\\\\\\\ ")
  if(sum(is.na(var))>0){
    no <- table(is.na(var),by)[2,]
    nop <- sapply(100*no/table(by),myper1)
    nopt <- sapply(100*sum(no)/sum(table(by)),myper1)
    if(missing(note)){ line2 <- paste("\\\\hspace{.1in}Missing",paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),paste(sum(no)," (",nopt,"\\\\%)",sep=""),sep=" & ")
    } else{line2 <- paste(paste("\\\\hspace{.1in}Missing\\\\tnote{",note,"}",sep=""),paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),paste(sum(no)," (",nopt,"\\\\%)",sep=""),sep=" & ")
    }
    all.lines.miss <- paste(line1,line2,paste(lines,collapse=" \\\\\\\\ "),sep=" \\\\\\\\ ")
    return(list(line=all.lines,line.incmiss=all.lines.miss))
  }
  else{return(list(line=all.lines))}
}
           
getrangetex <- function(var,by,label="",by.order=FALSE,digits=8,formatper="myper0",rangesep=" -- "){
  label <- ifelse(label=="",label(var),label)
  if(label==""){stop("label was not specified, either as an attribute of var, or in the getrangetex call.")}
  myper1 <- get(formatper)
  if(!is.numeric(var)){stop("The variable is non-numeric")}
  s <- round(summary(var),digits=digits)
  text <- paste(s[1],rangesep,s[6],sep="")
  level <- sort(unique(by))
  sby <- matrix(nrow=length(level), ncol=length(s)); textby <- rep(NA,length(level))
  for(i in 1:length(level)){
      sby[i,1:6] <- round(summary(var[by==level[i]])[1:6],digits=digits)
      textby[i] <- ifelse(!is.na(sby[i,3]),paste(sby[i,1],rangesep,sby[i,6],sep=""), "--")
  }
  line1 <- paste(paste0(label,", Range"),paste(textby,collapse=' & '),text,sep=' & ')
  return(list(line=line1))
}
## PRINT LATEX HEADER FOR GETNPTEX, GETMEDTEX, and GETNP1LINE functions
printheader <- function(crosstab,total=TRUE,total.label="Combined",pval=FALSE,pval.label="P-value"){
if(total){
  x <- paste(c("",names(crosstab),total.label),collapse=" & ")
  y <- paste(c("",paste("(n=",c(crosstab,sum(crosstab)),")",sep="")),collapse=" & ")
  }
if(!total){
  x <- paste(c("",names(crosstab)),collapse=" & ")
  y <- paste(c("",paste("(n=",c(crosstab),")",sep="")),collapse=" & ")
  }
if(pval){x <- paste(c(x,pval.label),collapse=" & ")}
  line <- paste(c(x,y),collapse=" \\\\\\\\ ")
  return(line)
}
myformat2 <- function(x){format(round(x,2),nsmall=2)}
myformat4 <- function(x){
    y <- NULL
    if(x < 0.001){y <- "<0.001"}
    else{y <- format(round(x,3),nsmall=3)}
    return(y)
}
myper0 <- function(x){
    y <- NULL
    if(x <= 0.5 & x!=0){y <- "<1"}
    else{y <- format(round(x),nsmall=0)}
    return(y)
}

convertdate1 <- function(mydate){as.Date(mydate,"%Y-%m-%d")}
convertdate <- function(mydate){as.Date(mydate,"%m/%d/%Y")}

analysis$lastdate <- apply(with(analysis,data.frame(haart_start_date,dob,lastclin,lastpharm,lastlab,deathdt,ac_pc_enrl_dt)),1,max,na.rm=TRUE)
analysis$lastdate <- convertdate1(analysis$lastdate)
analysis$en2last <- as.numeric(analysis$lastdate - analysis$ac_pc_enrl_dt )
analysis$tx2last <- as.numeric(analysis$lastdate - analysis$haart_start_date )
analysis$en2tx <- as.numeric(analysis$haart_start_date - analysis$ac_pc_enrl_dt )
analysis$txin90 <- ifelse(is.na(analysis$en2tx) | analysis$en2tx < 0 | analysis$en2tx>90,0,1)
analysis$txbmi <- round(analysis$txwt / (analysis$txht/100)^2,1)
analysis$txbmi[analysis$txbmi<10 & !is.na(analysis$txbmi)] <- NA
analysis$txbmi[analysis$txbmi>50 & !is.na(analysis$txbmi)] <- NA
analysis$enbmi <- round(analysis$enwt / (analysis$enht/100)^2,1)
analysis$enbmi[analysis$enbmi<10 & !is.na(analysis$enbmi)] <- NA
analysis$enbmi[analysis$enbmi>50 & !is.na(analysis$enbmi)] <- NA
analysis$ph_entr_pnt[analysis$ph_entr_pnt %in% c("","NGR42900430","Negative","Positive")] <- NA
analysis$ph_entr_pnt <- as.factor(analysis$ph_entr_pnt)
levels(analysis$ph_entr_pnt)[2] <- "Outside clinic/provider"
analysis$ph_mrtl_stt[analysis$ph_mrtl_stt %in% c("","HIV-related","N/A","Negative","Not HIV-related","Positive","Unknown")] <- NA
analysis$ph_mrtl_stt[analysis$ph_mrtl_stt %in% c("Divorced","Separated")] <- "Divorced/Separated"
analysis$ph_pt_edu[analysis$ph_pt_edu %in% c("","N/A","Other")] <- NA
analysis$ph_pt_occp[analysis$ph_pt_occp %in% c("","N/A")] <- NA
analysis$cd4cat <- cut(analysis$encd4,breaks=c(0,50,200,350,1801))
analysis$hemcat <- cut(analysis$enhem,breaks=c(0,8,10,20))
## create an ART eligible variable
analysis$arteligible <- rep(NA,nrow(analysis))
analysis$arteligible[analysis$ac_pc_enrl_dt < convertdate("06/01/2010") & analysis$encd4 < 200 & !is.na(analysis$encd4)] <- 1
analysis$arteligible[analysis$ac_pc_enrl_dt < convertdate("06/01/2010") & analysis$encd4 < 350 & !is.na(analysis$encd4) & analysis$enwho==3 & !is.na(analysis$enwho)] <- 1
analysis$arteligible[analysis$ac_pc_enrl_dt < convertdate("06/01/2010") & analysis$enwho==4 & !is.na(analysis$enwho)] <- 1
analysis$arteligible[analysis$ac_pc_enrl_dt >= convertdate("06/01/2010") & analysis$encd4 <= 350 & !is.na(analysis$encd4)] <- 1
analysis$arteligible[analysis$ac_pc_enrl_dt >= convertdate("06/01/2010") & analysis$enwho==3 & !is.na(analysis$enwho)] <- 1
analysis$arteligible[analysis$ac_pc_enrl_dt >= convertdate("06/01/2010") & analysis$enwho==4 & !is.na(analysis$enwho)] <- 1
analysis$noassess <- rep(NA,nrow(analysis))
analysis$noassess <- ifelse(!is.na(analysis$encd4) | !is.na(analysis$arteligible),0,1)
analysis$ph_pt_edu[!is.na(analysis$ph_pt_edu) & analysis$ph_pt_edu=="Qur'anic"] <- "None"
analysis$ph_pt_edu[!is.na(analysis$ph_pt_edu) & analysis$ph_pt_edu %in% c("Completed primary","Started primary")] <- "Primary"
edu_codes <- data.frame(unique(analysis$ph_pt_edu)[c(3,4,2,5)],1:4)
analysis$education <- factor(edu_codes[match(analysis$ph_pt_edu,edu_codes[,1]),2])
levels(analysis$education) <- edu_codes[,1]
## MEDICATION
analysis$mednum <- unlist(lapply(strsplit(analysis$med,"+",fixed=TRUE),FUN=function(x) length(x[!is.na(x)])))
analysis$mednum[analysis$mednum==0 & analysis$art_status!="No HAART"] <- NA
analysis$mednum[analysis$mednum==1] <- 2 # collapse 1 and 2 to prophylaxis
analysis$mednum[analysis$mednum==4] <- 3 # collapse 3 and 4 to tripe
analysis$mednum <- factor(analysis$mednum,labels=c("No ART","Prophylaxis","HAART"))
## TRACKING DATA
analysis$endeath12mo<-ifelse(analysis$death=="Death" & analysis$en2last<365.25,1,0)
analysis$enfup12mo <- ifelse(analysis$en2last<365.25,analysis$en2last,365.25)
analysis$enltfu12mo <- ifelse(analysis$endeath12mo==0 & analysis$en2last<365.25 & convertdate("05/31/2011") - analysis$lastdate > 180,1,0)
analysis <- droplevels(analysis)

analysis$EID <- as.factor(analysis$EID)
PREG.N <- table(analysis$EID)
names(PREG.N) <- c("No EID record","EID entry")
PREG.HEAD <- printheader(PREG.N,pval=TRUE,pval.label="P-value\\\\tnote{b}")

age.m1 <- getmedtex(var=analysis$age,by=analysis$EID,label="Age",note="a",inc.pval=TRUE,digits=0)  
ager.m1 <- getrangetex(var=analysis$age,by=analysis$EID,label="Age",digits=0)  
psc.n1 <- getnptex(var=analysis$pregsource,by=analysis$EID,label="Source of Pregnancy Information",inc.pval=FALSE)
edu.n1 <- getnptex(var=analysis$education,by=analysis$EID,label="Education",inc.pval=TRUE)
mar.n1 <- getnptex(var=analysis$ph_mrtl_stt,by=analysis$EID,label="Marital status",inc.pval=TRUE)
occ.n1 <- getnptex(var=analysis$ph_pt_occp,by=analysis$EID,label="Occupation",inc.pval=TRUE)
cli.n1 <- getnptex(var=analysis$clinic,by=analysis$EID,label="Clinic",inc.pval=TRUE)
analysis$med[analysis$med %in% c("3TC + AZT + TDF","3TC + ABC + AZT","3TC + AZT + FTC + TDF","3TC + AZT + EFV + NVP","3TC + AZT + FTC + NVP + TDF","3TC + AZT + EFV + FTC + TDF","3TC + ABC + EFV + FTC","3TC + ABC + AZT + NVP","")] <- NA
reg.n1 <- getnptex(var=analysis$med,by=analysis$EID,label="analysis ART regimen",inc.pval=FALSE)
ref.n1 <- getnptex(var=analysis$ph_entr_pnt,by=analysis$EID,label="Referral type",inc.pval=TRUE)
hei.m1 <- getmedtex(var=analysis$enht,by=analysis$EID,label="Height (cm)",note="ac",inc.pval=TRUE,digits=0)
wei.m1 <- getmedtex(var=analysis$enwt,by=analysis$EID,label="Weight (kg)",note="ac",inc.pval=TRUE,digits=0)
bmi.m1 <- getmedtex(var=analysis$enbmi   ,by=analysis$EID,label="BMI (kg/m$^2$)",note="ac",inc.pval=TRUE,lower.label=FALSE,digits=1)
encd4.m1 <- getmedtex(var=analysis$encd4,by=analysis$EID,label="CD4 count",note="a",inc.pval=TRUE,lower.label=FALSE,digits=0)
encd4.n1 <- getnptex(var=analysis$cd4cat,by=analysis$EID,label="CD4 count category",inc.pval=FALSE,rowlabel=c("<50","51-200","201-350",">350"))
enhem.m1 <- getmedtex(var=analysis$enhem,by=analysis$EID,label="Hemoglobin",note="a",inc.pval=TRUE,digits=1)
enhem.n1 <- getnptex(var=analysis$hemcat,by=analysis$EID,label="Hemoglobin category",inc.pval=FALSE,rowlabel=c("<8","8-10",">10"))
encre.m1 <- getmedtex(var=analysis$encre,by=analysis$EID,label="Creatinine",note="a",inc.pval=TRUE,digits=2)
enwho.n1 <- getnptex(var=analysis$enwho,by=analysis$EID,label="WHO stage",inc.pval=TRUE,rowlabel=c("I","II","III","IV"))
enfs.n1 <- getnptex(var=analysis$enfs,by=analysis$EID,label="Functional status",inc.pval=TRUE)
art.n1 <- getnptex(var=analysis$mednum,by=analysis$EID,label="First ART receipt",inc.pval=TRUE)
table1  <- paste(age.m1$line,ager.m1$line,psc.n1$line,mar.n1$line.incmiss,occ.n1$line.incmiss,cli.n1$line,
                 ref.n1$line.incmiss,sep=" \\\\\\\\ ")
table2  <- paste(hei.m1$line.incmiss,wei.m1$line.incmiss,bmi.m1$line.incmiss,enfs.n1$line.incmiss,
		             encd4.m1$line.incmiss,encd4.n1$line.incmiss,enhem.m1$line.incmiss,encre.m1$line.incmiss,enwho.n1$line.incmiss, 
                 art.n1$line,sep=" \\\\\\\\ ")
setwd("~/Projects/Nigeria/Aim3")
firsten <- sapply(as.Date(unlist(lapply(split(analysis$ac_pc_enrl_dt,analysis$clinic),min)), origin = "1970-01-01"),format,"%e %B %Y")
firsttx <- sapply(as.Date(unlist(lapply(split(analysis$haart_start_date,analysis$clinic),min,na.rm=TRUE)), origin = "1970-01-01"),format,"%e %B %Y")
firstpt <- paste(paste(names(firsten),firsten,firsttx,sep= " & "),collapse=" \\\\\\\\ ")

## OUTCOMES
table(analysis$pcr.2.result,useNA='always',analysis$EID)
table(analysis$pcr.2.result,useNA='always',analysis$EID)
table(analysis$rapid.test.result,useNA='always',analysis$EID)
analysis$ctx.initiation.rev <- analysis$ctx.initiation
analysis$ctx.initiation.rev[!(analysis$ctx.initiation.rev %in% c(1,2,3)) & !is.na(analysis$ctx.initiation.rev)] <- NA
table(analysis$ctx.initiation.rev,useNA='always',analysis$EID)
table(analysis$feeding.option.in.1st.3.months,useNA='always',analysis$EID)

## VITAL STATUS DATA IS NOT RELIABLE
dead <- paste(analysis$dead.at.6mo,analysis$dead.at.12mo,analysis$dead.at.12mo.1)
well <- paste(analysis$well.at.6mo,analysis$well.at.12mo,analysis$well.at.12mo.1)
sick <- paste(analysis$sick.at.6mo,analysis$sick.at.12mo,analysis$sick.at.12mo.1)
write.csv(analysis[analysis$EID=="yes",],"~/Projects/Nigeria/Output/eid_20140107.csv",row.names=FALSE)
@


\documentclass[11pt]{article}
\usepackage[margin=0.75in]{geometry}
\usepackage{setspace,relsize}   % needed for latex(describe()), \code
\usepackage{moreverb}           % for verbatimtabinput
\usepackage{url}
\usepackage{amssymb}
\usepackage{color}
\usepackage[pdftex]{lscape}   % allows tables to be landscape when
\usepackage{longtable}  % for subject listings 
\usepackage{listings}
\lstset{breaklines=true}
\lstset{basicstyle=\small\ttfamily}
\usepackage{graphicx} 
\usepackage{graphicx, subfigure}
\usepackage{threeparttable}
\usepackage[breaklinks,colorlinks=true]{hyperref}  %% ,citecolor=blue
\usepackage{rotating}
\usepackage{tikz}
\usepackage[utf8]{inputenc}
\usetikzlibrary{shapes,shadows,arrows,trees}
\usepackage{caption}
\newcommand{\HRule}{\rule{\linewidth}{0.5mm}}
\newcommand{\squishlist}{
 \begin{list}{$\bullet$}
  { \setlength{\itemsep}{0pt}
     \setlength{\parsep}{3pt}
     \setlength{\topsep}{3pt}
     \setlength{\partopsep}{0pt}
     \setlength{\leftmargin}{1.5em}
     \setlength{\labelwidth}{1em}
     \setlength{\labelsep}{0.5em} } }

\newcommand{\squishlisttwo}{
 \begin{list}{$\bullet$}
  { \setlength{\itemsep}{0pt}
     \setlength{\parsep}{0pt}
    \setlength{\topsep}{0pt}
    \setlength{\partopsep}{0pt}
    \setlength{\leftmargin}{2em}
    \setlength{\labelwidth}{1.5em}
    \setlength{\labelsep}{0.5em} } }

\newcommand{\squishend}{
  \end{list}  }


\usepackage{Sweave}
\begin{document}
 \SweaveOpts{prefix.string=graphics/plot} % Created a "graphics" subdirectory to save graph files in
\begin{titlepage}
  \begin{center}
    % Upper part of the page
    \includegraphics[scale=0.15]{/home/blevinml/Projects/MetaCytology/Data/biostat-logo-mb.png}\\[2cm]

    \textsc{\Large Statistical Analysis Report:}\\[0.5cm]
 
    % Title
    \HRule \\[0.4cm]
    { \Large \bfseries Pregnant patient characteristics by Early Infant Diagnosis entry  \\[0.15cm] in rural Kwara and Niger states of Nigeria}\\[0.4cm]
    \HRule \\[1.5cm]

    % Biostatistician and Faculty Member
\begin{minipage}{0.4\textwidth}
\begin{flushleft} \large
\emph{Biostatistician:}\\
Meridith Blevins, M.S.
\end{flushleft}
\end{minipage}
\begin{minipage}{0.4\textwidth}
\begin{flushright} \large
\emph{Faculty Statistician:} \\
Bryan Shepherd, Ph.D.
\end{flushright}
\end{minipage}
 
\vfill
 
% Bottom of the page
{\large \today}
 
\end{center}
 
\end{titlepage}

\section{Introduction}
Niger and Kwara states are located in Nigeria's North Central region, the region with the highest HIV prevalence
 in the country (2010 adult HIV prevalence of 7.5\% vs. national HIV prevalence of 3.6\%).
 Both Ilorin, the capital of Kwara State and Minna, capital of Niger State are located on major
 highways that link the north and south of Nigeria, with substantial populations of HIV at-risk groups.
 The prevalence of HIV among adults in Niger and Kwara states is estimated to be 4.0\% and 2.2\%, respectively.\\

VIGH/FGH-supported activities in Nigeria utilize an integrated and holistic health system strengthening approach,
 organized around two core principles: 1) Direct technical assistance to government health facilities to implement
 integrated HIV clinical services at the state level; and 2) Health workforce capacity development.
 VIGH/FGH-supported activities include the following program areas: adult and pediatric HIV care and treatment
 (HCT), prevention of mother-to-child HIV transmission (PMTCT), HIV testing and counseling (HTC), TB-HIV,
 orphans and vulnerable children (OVC), pharmaceutical logistics including procurement and provision
 of antiretroviral drugs (ARVs), strengthening of laboratory infrastructure, and improvement of
 strategic information (SI) services. At the time of this study FGH was supporting HIV treatment services
 in five clinics, namely: Sobi Specialist Hospital (Ilorin) and Lafiagi General Hospital in Kwara state;
 and Gawu Babangida Rural Hospital, Kuta Rural Hospital, and Umaru Yar Adua Hospital Sabon Wuse in Niger State.
 These five clinics are served by nine HTC/PMTCT feeder sites (‘satellite’ sites) in a hub and spoke model.
 
\subsection{Research Aims}
\vspace{-0.07in}
\subsubsection{To describe pregnant patients enrolled in HIV care and treatment.}
\vspace{-0.07in}
\subsubsection{To describe women who bring their babies to Early Infant Diagnosis (EID).}
\vspace{-0.07in}
\subsubsection{To determine characteristics associated with EID initiation.}

\section{Methods}
\subsection{Participants}
This analysis uses data collected through September 30, 2013 and include HIV-infected pregnant patients 
 entering HIV care and treatment programs at aged 12 and older. 
 In order to ensure that the mothers had enough time to carry the baby to term and go for EID, we 
 include only women enrolled on or before December 31, 2012.   Only four of five clinics will be included
 for analysis: Sobi, Lafiagi, Gaw, and Umaru Yar Adua; Kuta is excluded because there EID records are not kept 
 for that population.  Additionally, EID records are kept for IBB, but CAREWare data is not collected for that site;
 IBB EID data will not be included for analysis.  

\subsection{Outcomes}
The primary outcome will be entry into the EID program.  The unit of analysis is a woman documented as pregnant at any time
 following enrollment into HIV care and treatment.  The HIV care and treatment database does not have due dates,
 so the time of pregnancy cannot be assessed for women who did not being their baby for EID, nor can we
 account for any one mom who bring in babies from separate pregnancies (i.e. she will be credited once for EID
 entry as long as one or more babes entered the EID program).    \\

\subsection{Multiple centers}
There are five clinics in two states.  Sobi and Lafiagi Hospitals are located in Kwara state.
 Gawu, Kuta, and Umaru YarAdua Hospitals are located in Niger State.  Clinic is a covariate of interest; thus
 we will not account for clustering within clinics.  % Significant heterogeneity
 % is not anticipated among these five hospitals; however, Cox models will be stratified by clinic.

\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Dates of Care and Treatment Initiation by Clinic}
\label{tab:base}
\begin{tabular}{llllllll}
\hline
Clinic & Care Initiation & Treatment Initiation\\
\hline 
\Sexpr{firstpt}\\
%\\
\hline
\end{tabular}
}
\begin{tablenotes}
\item[a] Note that the database close is 30 September 2013, but the inclusion criteria
 was women enrolled by 31 December 2012.
\end{tablenotes}
\end{threeparttable}
\end{center}

\subsection{Data Sources and Measurements}

The data extraction includes all patients enrolled on or before 30 September 2013. 
 Clinical characteristics collected closest to the date of enrollment up to a 90 day window before or
 after enrollment are deemed as enrollment status indicators (eg, weight, WHO stage, CD4)
 with the exception of height which was allowed a 365 day window in either direction. \\

VU/FGH uses CAREWare to record client level data at all HIV care and treatment sites. 
  VU/FGH performs routine audits of medical records to ensure forms are completed accurately and laboratory data is entered correctly. 
 VU/FGH is in the process of improving data management system for CAREWare to better ensure quality data. 
 A data management plan has been developed and is in the final stage of being fully implemented in the coming months.  \\

Demographic information includes: sex, age, marital status, educational status, occupational status, service entry into the program.
 Routine clinical data are as follows: weight, height, functional status, WHO stage, TB status.  Labs that may be collected include:
 CD4, hepatitis B, VDRL, pregnancy, hepatitis C, hematocrit,
 hemoglobin, ALT, creatinine, total cholesterol, HDL, LDL, and triglycerides.

\subsubsection{Data Cleaning}

Data queries were generated for out of range and missing data.   Each site addressed data queries and the clean data was
 extracted for final analyses.  In order to ensure that our denominator (pregnant women in HIV care and treatment) was 
 complete, we performed 100\% audit of pregnancy status of all women aged 12-55 against an external database (HIVCare) maintained
 by FGH personnel (Sodangi) for internal reporting.  \\

Along with data query resolution, values of hemoglobin less than 1 and greater than 18 (were imputed
as missing. Similarly, values of CD4 count greater than 1500 were imputed as missing.
Creatinine larger than 100 were imputed as missing. If height was recorded less than 100 
or larger than 220, it is imputed as missing. If weight was recorded less than 20 or larger than 120, it is
imputed as missing. BMI is calculated as weight in kilograms divided by the square of height in meters.
Extreme BMI records below 10 kg/m$^2$ and above kg/m$^2$ were imputed as missing.

\subsection{Statistical methods}

{\it 1. To describe pregnant patients enrolled in HIV care and treatment.}\\
Summary characteristics will be tabulated for all pregnant women.\\

\noindent{\it 2. To describe women who bring their babies to Early Infant Diagnosis (EID).}\\
Summary characteristics will be tabulated for all pregnant women by entry into the Early Infant Diagnosis program.\\

\noindent{\it 3. To determine characteristics associated with EID initiation.}\\
We modeled the probability of EID initiation using multivariable logistic regression.  Covariates identified by the study primary investigator {\it apriori}
 include: age, marital status, clinic, BMI, functional status, CD4 count, WHO stage, ART status, date of enrollment. 
% Results from the model
 % showed little incremental association between covariates and EID entry after adjusting for program level covariates (clinic and date of
 % enrollment).  A second model was established  {\it pos hoc} which dropped clinic and date of enrollment.  Both models are presented.
 Linearity assumptions were met for all continuous covariates. 
  Multiple imputation was used to account for missing values of covariates and to prevent case-wise deletion of missing data;
  180 (25\%) patients had complete data for all covariates.
  We used predictive mean matching to take random draws from imputation models; 25 imputation data sets were used in the analysis.\footnote[1]{Harrell Jr FE, Dupont MC. Package `Hmisc’. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. 2012.}
  

\noindent R-software 3.0.2 (www.r-project.org) will be used for data analyses.

    \listoftables\addcontentsline{toc}{section}{Index of Tables}
    \listoffigures\addcontentsline{toc}{section}{Index of Figures}


\section{Results}
Data include all patients enrolled into care on or before 30 September 2013.


%%%%%%%%%%%%%%%%%%%%
%     TABLE 1      %
%%%%%%%%%%%%%%%%%%%%
\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Summary of Patient Demographics by EID entry}
\label{tab:base}
\begin{tabular}{llllllll}
\hline
\Sexpr{PREG.HEAD }\\
\hline 
\Sexpr{table1}\\
%\\
\hline
\end{tabular}
}
\begin{tablenotes}
\item[a] Continuous variables are reported as medians (interquartile range).
\item[b] To compare the distribution of study characteristics for participants by sex,
 we employ chi-square tests.   Similarly, we use a Wilcoxon rank sum test for continuous variables by sex.
\end{tablenotes}
\end{threeparttable}
\end{center}


%%%%%%%%%%%%%%%%%%%%
%     TABLE 1      %
%%%%%%%%%%%%%%%%%%%%
\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Summary of Patient Clinical Characteristics by EID entry}
\label{tab:base}
\begin{tabular}{llllllll}
\hline
\Sexpr{PREG.HEAD}\\
\hline 
\Sexpr{table2}\\
%\\
\hline
\end{tabular}
}
\begin{tablenotes}
\item[a] Continuous variables are reported as medians (interquartile range).
\item[b] To compare the distribution of study characteristics for participants by sex,
 we employ chi-square tests.   Similarly, we use a Wilcoxon rank sum test for continuous variables by sex.
\item[c] Weight, height, functional status, CD4, hemoglobin, creatinine, WHO stage are collected at enrollment. Enrollment data is collected in a window of $\pm$ 90 days from date of enrollment.
\end{tablenotes}
\end{threeparttable}
\end{center}


% Hidden R code chunk --- reading in the data
<<echo = false, results = hide>>=
library(rms)
set.seed(2013111)
lrm1.d1 <- analysis
lrm1.d1$enroll_dt <- as.numeric(lrm1.d1$ac_pc_enrl_dt) - 13000
lrm1.d1$enwho[lrm1.d1$enwho %in% 3:4] <- "3/4"
lrm1.d1$enwho <- as.factor(lrm1.d1$enwho)
lrm1.d1$enfs[lrm1.d1$enfs %in% c("Ambulatory","Bedridden")] <- "Ambulatory/Bedridden"
table(lrm1.d1$art_status,lrm1.d1$med)
lrm1.d1$EID <- ifelse(lrm1.d1$EID=="yes",1,0)
lrm1.d1$clinic <- as.factor(lrm1.d1$clinic)
lrm1.d1$enfs <- as.factor(lrm1.d1$enfs)
lrm1.d1$ph_mrtl_stt <- as.factor(lrm1.d1$ph_mrtl_stt)
lrm1.d1$education <- as.factor(lrm1.d1$education)
lrm1.dd1 <- datadist(with(lrm1.d1,data.frame(encd4, enhem, enwho, enfs, clinic, ph_mrtl_stt, education, enroll_dt, age, enbmi,mednum,art_status,EID)))

lrm1.mi <- aregImpute(~ encd4 + enhem + enwho + enfs + ph_mrtl_stt + education + enroll_dt + age + enbmi + mednum + EID,n.impute=25,data=lrm1.d1) 
options(datadist='lrm1.dd1')

########## MODEL 1: LOGISITIC REGRESSIOn (ALL COVARIATES SELECTED APRIORI)
lrm1.nl <- fit.mult.impute(EID ~ clinic + rcs(age,4) + education + ph_mrtl_stt + rcs(enbmi,4) + enfs + rcs(encd4,4) + rcs(enhem,4) + enwho + rcs(enroll_dt,5) + mednum, lrm, lrm1.mi, lrm1.d1)
lrm1 <- fit.mult.impute(EID ~ clinic + age + education + ph_mrtl_stt + enbmi + enfs + encd4 + enhem + enwho + rcs(enroll_dt,4) + mednum, lrm, lrm1.mi, lrm1.d1)
lrm1.nomi <- lrm(EID ~ clinic + age + education + ph_mrtl_stt + enbmi + enfs + encd4 + enhem + enwho + enroll_dt + mednum, lrm1.d1)
getlrm1 <- as.data.frame(anova(lrm1))$P


cdates <- as.numeric(as.Date(c("01/01/2010","01/01/2011","01/01/2012"), "%m/%d/%Y"))  - 13000
indexsum1 <-c(2,18,20,22,24,26,28,4,30,6,8,32,34,36,38,12,14,16,10)
lrmsum1 <- summary(lrm1,clinic="UMYMH",age=c(20,25),education="None",ph_mrtl_stt="Married",enbmi=c(20,21),enfs="Working",
                        encd4=c(50,100),enhem=c(10,11),enwho=1,enroll_dt=c(cdates[1],cdates[2]),mednum="No ART")[indexsum1,c(4,6,7)]
lrmsum2 <- summary(lrm1,clinic="UMYMH",age=c(20,25),education="None",ph_mrtl_stt="Married",enbmi=c(20,21),enfs="Working",
                        encd4=c(50,100),enhem=c(10,11),enwho=1,enroll_dt=c(cdates[1],cdates[3]),mednum="No ART")[10,c(4,6,7)]
lrmsum <- rbind(lrmsum1,lrmsum2)
lrm1.labels <- c("Age (per 5 yrs)","Education","\\\\hspace{0.25in} None (ref)","\\\\hspace{0.25in} Primary",
	                                       "\\\\hspace{0.25in} Secondary","\\\\hspace{0.25in} Post secondary",
	         "Marital status","\\\\hspace{0.25in} Married (ref)","\\\\hspace{0.25in} Divorced/Separated","\\\\hspace{0.25in} Single","\\\\hspace{0.25in} Widowed",
                 "BMI (per 1 kg/m$^2$)",
	         "Functional status","\\\\hspace{0.25in} Working (ref)","\\\\hspace{0.25in} Ambulatory/Bedridden",
                 "CD4 count (per 50 cells/$\\\\mu$L)","Hemoglobin (per 1 g/dL)",
		 "WHO stage","\\\\hspace{0.25in} I (ref)","\\\\hspace{0.25in} II","\\\\hspace{0.25in} III or IV",
                 "First ART receipt","\\\\hspace{0.25in} No ART (ref)","\\\\hspace{0.25in} Prophylaxis","\\\\hspace{0.25in} HAART",
                 "Clinic","\\\\hspace{0.25in} UMYMH (ref)","\\\\hspace{0.25in} GBRH","\\\\hspace{0.25in} LGH","\\\\hspace{0.25in} SBSH",
                 "Date of enrollment","\\\\hspace{0.25in} January 2010 (ref)","\\\\hspace{0.25in} January 2011","\\\\hspace{0.25in} January 2012"
		  )
lrmsumlines1 <- paste(paste(sapply(lrmsum[,1],myformat2)," (",sapply(lrmsum[,2],myformat2),", ",sapply(lrmsum[,3],myformat2),")",sep=""),sep=" & ")
lrmsumlines2 <- c(lrmsumlines1[1],"","1",lrmsumlines1[2:4],"","1",lrmsumlines1[5:8],"","1",lrmsumlines1[9:11],"","1",lrmsumlines1[12:13],
                  "","1",lrmsumlines1[14:15],"","1",lrmsumlines1[16:18],"","1",lrmsumlines1[19:20])
lrm1pval <- c(sapply(getlrm1[2:3],myformat4),rep("",4),sapply(getlrm1[c(4)],myformat4),rep("",4),sapply(getlrm1[5:6],myformat4),rep("",2),
                     sapply(getlrm1[7:9],myformat4),rep("",3),myformat4(getlrm1[12]),rep("",3),
                     myformat4(getlrm1[1]),rep("",4),sapply(getlrm1[c(10)],myformat4),rep("",3))
lrm1table <- paste(lrm1.labels,lrmsumlines2,lrm1pval,sep=" & ")
lrm1tabletex <- paste(lrm1table,collapse=" \\\\\\\\ ")

p1 <- Predict(lrm1,enroll_dt)

@

\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Logistic Regression Model: Factors associated with EID uptake}
\label{tab:lrm1}
\begin{tabular}{llll}
\hline
            & OR (95\% CI)      & P-value     \\
\hline
\Sexpr{lrm1tabletex}\\
\hline
\end{tabular}
}
%\begin{tablenotes}
%\item[a]  
%\end{tablenotes}
\end{threeparttable}
\end{center}



\setkeys{Gin}{width=.9\textwidth} %default is 0.7
\begin{figure}[h!]
   \caption{Predicted log-odds of EID entry by enrollment in HIV care and treatment}
   \begin{center}
      \label{fig:lrm1}
<<fig=TRUE, echo=FALSE,height=3.5,width=7>>=
  myadj1 <- 1
  myadj2 <- -0.5

#png("graphics/vitDbdate.png",res=325,width=1750,height=1050, bg="transparent") 
days1 <- c("01/15/","2/15/","03/15/","04/15/","05/15/","06/15/","07/15/","08/15/","09/15/","10/15/","11/15/","12/15/")
months1 <- c("Jan-","Feb-","Mar-","Apr-","May-","Jun-","Jul-","Aug-","Sep-","Oct-","Nov-","Dec-")
datelabs1 <- as.numeric(as.Date(c(paste0(days1,"2009"),paste0(days1,"2010"),paste0(days1,"2011"),paste0(days1,"2012"),paste0(days1,"2013")), "%m/%d/%Y") ) - 13000
datelabs2 <- c(paste0(months1,"09"),paste0(months1,"10"),paste0(months1,"11"),paste0(months1,"12"),paste0(months1,"13"))
#datelabs2 <- c("Jan-10  ","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"," Jan-11")
print(plot(p1, ylab='Log-odds of EID entry',xlab="Date of enrollment",adj.subtitle=FALSE,data=lrm1.d1,lwd=2,scales=list(x=list(at=datelabs1[seq(1,60,by=6)], labels=datelabs2[seq(1,60,by=4)])),scat1d.opts=list(col="darkred",lwd=0.4,frac=0.05),
      par.settings = list(
        axis.components = list(
           left   = list(pad1 = myadj1, pad2 = myadj2)
          ,bottom = list(pad1 = myadj1, pad2 = myadj2)
          ,top    = list(pad1 = myadj1, pad2 = myadj2)
          ,right  = list(pad1 = myadj1, pad2 = myadj2)
        ))))
#dev.off()
@
  \end{center}
\begin{lstlisting}
Adjusted to: \Sexpr{attr(p1,"info")$adjust}
\end{lstlisting}
\end{figure} 


\end{document}