%#############################################################
%#
%#   Program: analysis-scripts-001.nw
%#	  
%#   Project: ART in rural Mozambique: Patient characteristics,
%#            trends, and considerations as donor philosophy
%#            transitions away from an emergency response  
%# 
%#   PI: Janeen Burlison, MPH, MSN, WHNP-BC
%#       Mohsin Sidat, MD, PhD
%#       Alfredo Vergara, PhD
%# 	 Sten Vermund, MD, PhD
%#	 Troy Moon, MD, MPH*
%#       * corresponding author, troy.moon at fgh.org.mz
%#	  
%#   Biostatistician/Programmer: Meridith Blevins, MS**
%#   ** contact programmer: meridith.blevins at vanderbilt.edu
%#	  
%#   Statistician: Bryan E Shepherd, PhD
%#	  
%#   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: 17 May 2010
%#	  
%#   Revisions: 
%#   
%#
%#############################################################


% RUN SWEAVE (WITH TIME STAMP):
%   system.time(Sweave("analysis-scripts-001.nw"))

\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{graphicx} 
\usepackage{graphicx, subfigure}
\usepackage{threeparttable}
\usepackage{rotating}
\newcommand{\HRule}{\rule{\linewidth}{0.5mm}}


\begin{document}
\SweaveOpts{prefix.string=graphics/plot} % Created a "graphics" subdirectory to save graph files in
\begin{titlepage}
  \begin{center}
     \textsc{\Large Statistical Report:}\\[0.5cm]
    % Title
    \HRule \\[0.4cm]
    { \Large \bfseries ART in rural Mozambique:\\ Patient characteristics, trends, and considerations}\\[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}

% Hidden R code chunk --- reading in the data
<<echo = false, results = hide>>=
##########################################
# READ IN FULL DATASET "ADULTMOZAMBIQUE" #
##########################################
source("read_data.R")
attach(adultmozambique)

###################################
# CUSTOMIZE DATA FORMAT FUNCTIONS #
###################################
myformat2 <- function(x){format(round(x,2),nsmall=2)}
myformat3 <- function(x){
    y <- NULL
    if(x < 0.01){y <- "<0.01"}
    else{y <- format(round(x,2),nsmall=2)}
    return(y)
}
myformat4 <- function(x){
    y <- NULL
    if(x < 0.001){y <- "<0.001"}
    else{y <- format(round(x,3),nsmall=3)}
    return(y)
}
myper1 <- function(x){
    y <- NULL
    if(x < 0.05 & x!=0){y <- "<0.1"}
    else{y <- format(round(x,1),nsmall=1)}
    return(y)
}
#############################################
# CUSTOMIZE FUNCTIONS TO WRITE LATEX TABLES #
#############################################
# WRITE 3 FUNCTIONS WHICH tabulate counts/proportions and median(q1,q3) overall and by covariate
getnp1line <- function(var,by,rowlabel,indent=TRUE){
  n <- table(var,by)
  d <- table(by)
  p <- sapply(100*n[1,]/d,myper1)
  pt <- sapply(100*table(var)/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=" & ")}
  return(list(line=line))
}
getmedtex <- function(var,by,label,note=note,inc.pval=TRUE,lower.label=TRUE){
  s <- round(summary(var),1)
  text <- paste(s[3]," (",s[2],", ",s[5],")",sep="")
  level <- sort(unique(by))
  sby <- matrix(nrow=length(level), ncol=length(s)); textby <- rep(NA,length(level))
  if(length(level)==1) {pval <- " "}
  else if(length(level)>1) {pval <- myformat3(kruskal.test(var ~ by)$p.value)}
  for(i in 1:length(level)){
      sby[i,1:5] <- round(summary(var[by==level[i]])[1:5],1)
      textby[i] <- paste(sby[i,3]," (",sby[i,2],", ",sby[i,5],")",sep="")    
  }
  if(inc.pval==TRUE){
      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 ",tolower(label),", n(\\\\%)",sep=""),paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),paste(sum(no)," (",nopt,"\\\\%)",sep=""),sep=" & ")}
      else{line2 <- paste(paste("\\\\hspace{.1in}Missing ",label,", 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))}
}
getnptex <- function(var,by,label,rowlabel,note=note,inc.pval=TRUE){
  n <- table(var,by)
  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,] <- sapply(100*n[i,]/d,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 PLANNED
  if(inc.pval==TRUE){
    pval <- myformat3(chisq.test(n)$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,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))}
}

##################################
# TABULATE/PROPORTIONS/QUARTILES #
##################################

# CREATE DATASETS OF ALL SUBGROUPS
  allmoz <- adultmozambique
# NOW WE SUBSET TO THOSE PARTICIPANTS WHO WERE ENROLLED BY 12/31/2008
# THEN THOSE WHO HAD TREATMENT WITHIN 90 DAYS OF ENROLLMENT ARE IN TX GROUP
  enrolledby <- as.Date(adultmozambique$dataabertura, "%m/%d/%Y") < as.Date("01/01/2009", "%m/%d/%Y")
  mymoz <- cbind(adultmozambique[enrolledby,],en2tx[enrolledby],education[enrolledby])
  mymoz$txin90 <- ifelse(mymoz$en2tx < 91 & !is.na(mymoz$en2tx),1,0)
# THEN SUBSET TO THOSE WHO HAD AT LEAST ONE VISIT
  fupmoz <- mymoz[!(as.character(mymoz$dataabertura)==as.character(mymoz$lastdate)),]
# FINALLY IDENTIFY ALL ELIGIBLE ART CANDIDATES BASED ON CD4 and/or WHO STAGE
  elig1 <- ifelse(mymoz$encd4<200 & !is.na(mymoz$encd4),1,0)
  elig2 <- ifelse(mymoz$enwho==4 & !is.na(mymoz$enwho),1,0)
  elig3 <- ifelse(mymoz$encd4>=200 & mymoz$encd4<350 & !is.na(mymoz$encd4) & mymoz$enwho==3 & !is.na(mymoz$enwho),1,0)
  eligmoz <- mymoz[rowSums(data.frame(elig1,elig2,elig3))>0,]

bysex.N <- table(allmoz$sexo)
bysex.age <- getmedtex(allmoz$idade,allmoz$sexo,label="Age",note="a",inc.pval=FALSE,lower.label=FALSE)
bysex.coddistrito <- getnptex(allmoz$coddistrito,allmoz$sexo,label="District",rowlabel=levels(allmoz$coddistrito),inc.pval=FALSE)
bysex.ed <- getnptex(allmoz$education,allmoz$sexo,label="Education Level",rowlabel=c("Illiterate","Primary","Secondary","Technical","University Level"),inc.pval=FALSE,note="b")
bysex.enwho <- getnptex(allmoz$enwho,allmoz$sexo,label="WHO Staging at Enrollment\\\\tnote{c}",rowlabel=c("I","II","III","IV"),inc.pval=FALSE)
bysex.enweight <- getmedtex(allmoz$enweight,allmoz$sexo,label="Weight at Enrollment",inc.pval=FALSE,lower.label=FALSE)
bysex.enheight <- getmedtex(allmoz$enheight,allmoz$sexo,label="Height at Enrollment",inc.pval=FALSE,lower.label=FALSE)
bysex.enBMI <- getmedtex(allmoz$enBMI,allmoz$sexo,label="BMI at Enrollment",inc.pval=FALSE,lower.label=FALSE)
bysex.encd4 <- getmedtex(allmoz$encd4,allmoz$sexo,label="CD4 at Enrollment",inc.pval=FALSE,lower.label=FALSE)
bysex.enhgb <- getmedtex(allmoz$enhgb,allmoz$sexo,label="HGB at Enrollment",inc.pval=FALSE,lower.label=FALSE)
bysex.emtarv <- getnptex(allmoz$emtarv,allmoz$sexo,label="ART Ever Initiated",rowlabel=c("No","Yes"),inc.pval=FALSE)

mymoz.N <- table(mymoz$txin90)
mymoz.sex <- getnp1line(mymoz$sexo,mymoz$txin90,rowlabel="Female",indent=FALSE)
mymoz.age <- getmedtex(mymoz$idade,mymoz$txin90,label="Age",note="a",inc.pval=FALSE,lower.label=FALSE)
mymoz.coddistrito <- getnptex(mymoz$coddistrito,mymoz$txin90,label="District",rowlabel=levels(allmoz$coddistrito),inc.pval=FALSE)
mymoz.ed <- getnptex(mymoz$education,mymoz$txin90,label="Education Level",rowlabel=c("Illiterate","Primary","Secondary","Technical","University Level"),inc.pval=FALSE,note="b")
mymoz.enwho <- getnptex(mymoz$enwho,mymoz$txin90,label="WHO Staging at Enrollment\\\\tnote{c}",rowlabel=c("I","II","III","IV"),inc.pval=FALSE)
mymoz.enweight <- getmedtex(mymoz$enweight,mymoz$txin90,label="Weight at Enrollment",inc.pval=FALSE,lower.label=FALSE)
mymoz.enheight <- getmedtex(mymoz$enheight,mymoz$txin90,label="Height at Enrollment",inc.pval=FALSE,lower.label=FALSE)
mymoz.enBMI <- getmedtex(mymoz$enBMI,mymoz$txin90,label="BMI at Enrollment",inc.pval=FALSE,lower.label=FALSE)
mymoz.encd4 <- getmedtex(mymoz$encd4,mymoz$txin90,label="CD4 at Enrollment",inc.pval=FALSE,lower.label=FALSE)
mymoz.enhgb <- getmedtex(mymoz$enhgb,mymoz$txin90,label="HGB at Enrollment",inc.pval=FALSE,lower.label=FALSE)

fupmoz.N <- table(fupmoz$txin90)
fupmoz.sex <- getnp1line(fupmoz$sexo,fupmoz$txin90,rowlabel="Female",indent=FALSE)
fupmoz.age <- getmedtex(fupmoz$idade,fupmoz$txin90,label="Age",note="a",inc.pval=FALSE,lower.label=FALSE)
fupmoz.coddistrito <- getnptex(fupmoz$coddistrito,fupmoz$txin90,label="District",rowlabel=levels(allmoz$coddistrito),inc.pval=FALSE)
fupmoz.ed <- getnptex(fupmoz$education,fupmoz$txin90,label="Education Level",rowlabel=c("Illiterate","Primary","Secondary","Technical","University Level"),inc.pval=FALSE,note="b")
fupmoz.enwho <- getnptex(fupmoz$enwho,fupmoz$txin90,label="WHO Staging at Enrollment\\\\tnote{c}",rowlabel=c("I","II","III","IV"),inc.pval=FALSE)
fupmoz.enweight <- getmedtex(fupmoz$enweight,fupmoz$txin90,label="Weight at Enrollment",inc.pval=FALSE,lower.label=FALSE)
fupmoz.enheight <- getmedtex(fupmoz$enheight,fupmoz$txin90,label="Height at Enrollment",inc.pval=FALSE,lower.label=FALSE)
fupmoz.enBMI <- getmedtex(fupmoz$enBMI,fupmoz$txin90,label="BMI at Enrollment",inc.pval=FALSE,lower.label=FALSE)
fupmoz.encd4 <- getmedtex(fupmoz$encd4,fupmoz$txin90,label="CD4 at Enrollment",inc.pval=FALSE,lower.label=FALSE)
fupmoz.enhgb <- getmedtex(fupmoz$enhgb,fupmoz$txin90,label="HGB at Enrollment",inc.pval=FALSE,lower.label=FALSE)

eligmoz.N <- table(eligmoz$txin90)
eligmoz.sex <- getnp1line(eligmoz$sexo,eligmoz$txin90,rowlabel="Female",indent=FALSE)
eligmoz.age <- getmedtex(eligmoz$idade,eligmoz$txin90,label="Age",note="a",inc.pval=FALSE,lower.label=FALSE)
eligmoz.coddistrito <- getnptex(eligmoz$coddistrito,eligmoz$txin90,label="District",rowlabel=levels(allmoz$coddistrito),inc.pval=FALSE)
eligmoz.ed <- getnptex(eligmoz$education,eligmoz$txin90,label="Education Level",rowlabel=c("Illiterate","Primary","Secondary","Technical","University Level"),inc.pval=FALSE,note="b")
eligmoz.enwho <- getnptex(eligmoz$enwho,eligmoz$txin90,label="WHO Staging at Enrollment\\\\tnote{c}",rowlabel=c("I","II","III","IV"),inc.pval=FALSE)
eligmoz.enweight <- getmedtex(eligmoz$enweight,eligmoz$txin90,label="Weight at Enrollment",inc.pval=FALSE,lower.label=FALSE)
eligmoz.enheight <- getmedtex(eligmoz$enheight,eligmoz$txin90,label="Height at Enrollment",inc.pval=FALSE,lower.label=FALSE)
eligmoz.enBMI <- getmedtex(eligmoz$enBMI,eligmoz$txin90,label="BMI at Enrollment",inc.pval=FALSE,lower.label=FALSE)
eligmoz.encd4 <- getmedtex(eligmoz$encd4,eligmoz$txin90,label="CD4 at Enrollment",inc.pval=FALSE,lower.label=FALSE)
eligmoz.enhgb <- getmedtex(eligmoz$enhgb,eligmoz$txin90,label="HGB at Enrollment",inc.pval=FALSE,lower.label=FALSE)
@

\section{Introduction}
Mozambique's Zamb\'ezia province is the country's second most populace province with a population of 
 approximately 4 million persons and an overall estimated HIV/AIDS prevalence of nearly 20\%.  VIGH and its 
 affiliate non-governmental organization Friends in Global Health, LLC (FGH) have been working 
 exclusively in Zamb\'ezia province since February 2007 as a PEPFAR funded partner to the Mozambican 
 Ministry of Health.  FGH is assisting local health authorities with the development and 
 implementation of a comprehensive model for rural antiretroviral therapy (ART) expansion that 
 integrates health and social services, addresses human resources and infrastructure constraints, 
 and develops a sustainable health care system.  

\subsection{Research Aims}
\subsubsection{To summarize the characteristics of adult patients entering HIV care and treatment programs in six districts of rural Zamb\'ezia, Mozambique.}
\subsubsection{To identify enrollment trends in sex and WHO staging.}
\subsubsection{To identify missing data trends for CD4 and Hemoglobin.}

\section{Methods}
\subsection{Participants}
The planned analysis will use data collected through March 31, 2009 and include HIV-infected patients 
 entering HIV care and treatment programs at age 15 years or older.  

\subsection{Outcomes}
The primary outcome will be ART initiation in the first 90 days.  In all 
 analyses, we plan to use an `intent-to-continue-treatment' approach and thus ignore subsequent 
 changes to treatment, including treatment interruptions and terminations.  \\

Of secondary interest are the associations between quarter of enrollment and baseline characteristics or data collection.  
 It is anticipated that as the program matured and enrollment numbers increased, then new participants 
 enrolled earlier in disease course and data collection improved.


\subsection{Statistical Methods}
{\it 1. To summarize the characteristics of adult patients entering HIV care and treatment
 programs in six districts of rural Zamb\'ezia, Mozambique.}\\
Summary characteristics will be tabulated over total enrollment and by those who initiated
 ART and those who did not.  A logistic regression may identify whether baseline demographics,
 lab results, and clinical assessment are predictive of starting therapy or not.  Specifically,
 year of enrollment, age, sex, education level, CD4 count, weight, hemoglobin, and WHO stage at entry
 have been identified as predictors of interest. If needed, missing data methods will be
 considered. Options include complete case analysis, dummy variable adjustment, or multiple
 imputation.\\
\\
{\it 2. To identify enrollment trends in sex and WHO staging.}\\
The relationship between sex and enrollment may be modeled using a logistic regression
 with sex as the outcome and date of enrollment as a predictor.  Similarly, the
 relationship between staging and enrollment may be modeled using proportional odds
 (ordinal logistic regression) with stage as the outcome and date of enrollment as a
 predictor.  Figures depicting proportion of female and proportions of WHO staged I, II,
 I, or IV participants by quarter of enrollment will also summarize any trend.  \\
\\
{\it 3. To identify data collection trends in CD4 and Hemoglobin.}\\
It is expected that as the program matured and more laboratory/data collection resources
 were provided for that collection of CD4 and Hemoglobin improved.
The relationship between presence of lab value and enrollment may be modeled
 using a logistic regression with missingness as the outcome and date of enrollment
 as a predictor. Figures depicting proportion of CD4 and proportions of hemoglobin
 collected by quarter of enrollment will also summarize any trend.  \\
\\
R-software 2.10.0 (www.r-project.org) is used for data analyses. 
 Analysis plan dated 22 July 2009; {\it post hoc} incorporation of data collection trends.


\section{Results}
Between June 2006 and December 2006, each of the six districts initiated enrollment into the
 HIV Care and Treatment program. By the end of March 2009 after nearly 3 years of enrollment,
 there were a total of \Sexpr{nrow(adultmozambique)} adult patients (aged 15 years or older) in the care and treatment program. 

\subsection{Summary Characteristics}
Table \ref{tab:allsex} summarizes all participants enrolled at all six districts by 31 March 2009 by sex. 
 Table \ref{tab:allsites90} summarizes all participants enrolled at all six districts before 1 January 2009 by ART initiation within 90 days of enrollment. 
 Table \ref{tab:allsites90fup} summarizes all participants enrolled at all six districts before 1 January 2009 who had at least one day of follow-up 
 by ART initiation within 90 days of enrollment.  
 Table \ref{tab:allsites90elig} summarizes all participants enrolled at all six districts before 1 January 2009 who are recommended for therapy based on CD4/WHO Stage
 by ART initiation within 90 days of enrollment (if a patient had a CD4 count <200/$\mu$L independent of WHO Stage, was in WHO Stage IV independent of CD4 count,
 or had a CD4 count between 200-349/$\mu$L and was documented in WHO Stage III, they were recommended for treatment
 initiation). 
 There are \Sexpr{nrow(eligmoz)} patients who met the criteria for ART initiation; \Sexpr{sum(eligmoz$txin90)} individuals
 initiated ART within the first 90 days of enrollment.


%%%%%%%%%%%%%%%%%%%%
% TABLE 1 - BY SEX %
%%%%%%%%%%%%%%%%%%%%

\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Summary of Patient Characteristics by Sex}
\label{tab:allsex}
\begin{tabular}{p{6cm}p{2.5cm}p{2.5cm}p{2.5cm}}
\hline
   & Female             & Male                  & Combined   \\
   & (n=\Sexpr{bysex.N[1]})  & (n=\Sexpr{bysex.N[2]})  & (n=\Sexpr{sum(bysex.N)}) \\
\hline 
\Sexpr{bysex.age$line}\\
\Sexpr{bysex.coddistrito$line}\\
\Sexpr{bysex.ed$line.incmiss}\\
\Sexpr{bysex.enwho$line.incmiss}\\
\Sexpr{bysex.enweight$line.incmiss}\\
\Sexpr{bysex.enheight$line.incmiss}\\
\Sexpr{bysex.enBMI$line.incmiss}\\
\Sexpr{bysex.encd4$line.incmiss}\\
\Sexpr{bysex.enhgb$line.incmiss}\\
\Sexpr{bysex.emtarv$line}\\
\hline
\end{tabular}
}
\begin{tablenotes}
\item[a] Continuous variables are reported as medians (interquartile range).
\item[b] Percentages are computed using the number of patients with a non-missing value.
\item[c] Enrollment data is collected in a window of $\pm$ 90 days from date of enrollment.
\end{tablenotes}
\end{threeparttable}
\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOW A TABLE FOR PPTS REGISTERED %
% BY 12/31/2008 AND STARTING ART  %
% WITHIN 90 DAYS OF REGISTRATION  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Summary of Patient Characteristics by ART Initiation}
\label{tab:allsites90}
\begin{tabular*}{0.80\textwidth}{lp{2.5cm}p{2.5cm}l}
\hline
  & Not starting ART within 90 days of Registration      & Starting ART within 90 days of Registration       & Combined \\
  & (n=\Sexpr{mymoz.N[1]})  & (n=\Sexpr{mymoz.N[2]})    & (n=\Sexpr{sum(mymoz.N)}) \\
\hline \\
\Sexpr{mymoz.sex$line}\\
\Sexpr{mymoz.age$line}\\
\Sexpr{mymoz.coddistrito$line}\\
\Sexpr{mymoz.ed$line.incmiss}\\
\Sexpr{mymoz.enwho$line.incmiss}\\
\Sexpr{mymoz.enweight$line.incmiss}\\
\Sexpr{mymoz.enheight$line.incmiss}\\
\Sexpr{mymoz.enBMI$line.incmiss}\\
\Sexpr{mymoz.encd4$line.incmiss}\\
\Sexpr{mymoz.enhgb$line.incmiss}\\
\hline
\end{tabular*}
}
\begin{tablenotes}
\item[a] Continuous variables are reported as medians (interquartile range).
\item[b] Percentages are computed using the number of patients with a non-missing value.
\item[c] Enrollment data is collected in a window of $\pm$ 90 days from date of enrollment.
\end{tablenotes}
\end{threeparttable}
\end{center}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOW A TABLE FOR PPTS REGISTERED %
% BY 12/31/2008 AND STARTING ART  %
% WITHIN 90 DAYS OF REGISTRATION  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Summary of Patient Characteristics for Participants with at least One Day of Follow-up by ART Initiation}
\label{tab:allsites90fup}
\begin{tabular*}{0.80\textwidth}{lp{2.5cm}p{2.5cm}l}
\hline
  & Not starting ART within 90 days of Registration      & Starting ART within 90 days of Registration       & Combined \\
  & (n=\Sexpr{fupmoz.N[1]})  & (n=\Sexpr{fupmoz.N[2]})    & (n=\Sexpr{sum(fupmoz.N)}) \\
\hline \\
\Sexpr{fupmoz.sex$line}\\
\Sexpr{fupmoz.age$line}\\
\Sexpr{fupmoz.coddistrito$line}\\
\Sexpr{fupmoz.ed$line.incmiss}\\
\Sexpr{fupmoz.enwho$line.incmiss}\\
\Sexpr{fupmoz.enweight$line.incmiss}\\
\Sexpr{fupmoz.enheight$line.incmiss}\\
\Sexpr{fupmoz.enBMI$line.incmiss}\\
\Sexpr{fupmoz.encd4$line.incmiss}\\
\Sexpr{fupmoz.enhgb$line.incmiss}\\
\hline
\end{tabular*}
}
\begin{tablenotes}
\item[a] Continuous variables are reported as medians (interquartile range).
\item[b] Percentages are computed using the number of patients with a non-missing value.
\item[c] Enrollment data is collected in a window of $\pm$ 90 days from date of enrollment.
\end{tablenotes}
\end{threeparttable}
\end{center}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  NOW A TABLE FOR PPTS ELIGIBLE     %
% FOR ART (CD4/WHO) AND STARTING ART %
%  WITHIN 90 DAYS OF REGISTRATION    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Summary of Patient Characteristics for Participants Recommended for Treatment Initiation}
\label{tab:allsites90elig}
\begin{tabular*}{0.80\textwidth}{lp{2.5cm}p{2.5cm}l}
\hline
  & Not starting ART within 90 days of Registration      & Starting ART within 90 days of Registration       & Combined \\
  & (n=\Sexpr{eligmoz.N[1]})  & (n=\Sexpr{eligmoz.N[2]})    & (n=\Sexpr{sum(eligmoz.N)}) \\
\hline \\
\Sexpr{eligmoz.sex$line}\\
\Sexpr{eligmoz.age$line}\\
\Sexpr{eligmoz.coddistrito$line}\\
\Sexpr{eligmoz.ed$line.incmiss}\\
\Sexpr{eligmoz.enwho$line.incmiss}\\
\Sexpr{eligmoz.enweight$line.incmiss}\\
\Sexpr{eligmoz.enheight$line.incmiss}\\
\Sexpr{eligmoz.enBMI$line.incmiss}\\
\Sexpr{eligmoz.encd4$line.incmiss}\\
\Sexpr{eligmoz.enhgb$line.incmiss}\\
\hline
\end{tabular*}
}
\begin{tablenotes}
\item[a] Continuous variables are reported as medians (interquartile range).
\item[b] Percentages are computed using the number of patients with a non-missing value.
\item[c] Enrollment data is collected in a window of $\pm$ 90 days from date of enrollment.
\end{tablenotes}
\end{threeparttable}
\end{center}


<<echo = false, results = hide>>=
# RUN LOGISTIC REGRESSION MODEL WITH ART ENTRY AS OUTCOME - USE MULTIPLE IMPUTATION FOR MISSING DATA
  library(Design)
  vcov.lrm <- vcov.default # FIX BUG IN DESIGN PACKAGE
# PREDICTORS OF ART INITIATION IN FIRST 90 DAYS FOLLOWING ENROLLMENT
  set.seed(2)
  myenwho <- as.factor(mymoz$enwho)
  myenroll  <- as.numeric(as.Date(mymoz$dataabertura, "%m/%d/%Y")) - 13000
  myeducation <- as.factor(mymoz$education)
# SET VARIABLES IN DATA FRAME
  lrm1.d1 <- data.frame(mymoz$txin90,mymoz$coddistrito,mymoz$sexo,
                               mymoz$idade,mymoz$encd4,mymoz$enweight,
                               myenroll,mymoz$enhgb,myenwho,myeducation)
  colnames(lrm1.d1) <- c("txin90","site","sexo","idade","encd4","enweight","myenroll","enhgb","myenwho","education")
  lrm1.dd1 <- datadist(lrm1.d1)
# IMPUTE MISSING CD4, HGB, WEIGHT and WHO STAGE
  lrm1.mi <- aregImpute(~txin90+site+sexo+idade+education+encd4+enweight+enhgb+myenwho+myenroll,n.impute=10,data=lrm1.d1)
  options(datadist='lrm1.dd1')
# FIT LOGISTIC REGRESSION WITH DEFAULT N=10 IMPUTATIONS
  lrm1 <- fit.mult.impute(txin90 ~ site + sexo + rcs(idade,4) + education +  rcs(sqrt(encd4),4) +  
                        rcs(enweight,4) + enhgb + myenwho + rcs(myenroll,4), lrm, lrm1.mi, lrm1.d1)
  (getanova <- as.data.frame(anova(lrm1),row.names=c("site","sexo","idade","nonlin1","education","encd4",
             "nonlin2","enweight","nonlin3","enhgb","myenwho","myenroll","nonlin4","totnon","total"))[c(1:3,5:6,8,10:12),])
  (getnonlrm1 <- as.data.frame(anova(lrm1),row.names=c("site","sexo","idade","age nonlinear","education","encd4",
             "CD4 nonlinear","enweight","Weight nonlinear","enhgb","myenwho","myenroll","Date of Initiation nonlinear","totnon","total"))[c(4,7,9,13),])
  lrm1.stats <- lrm1$stats
  cdates <- as.numeric(as.Date(c("01/01/2006","01/01/2007","01/01/2008","01/01/2009"), "%m/%d/%Y")) - 13000
  mymoda <- summary(lrm1,myenroll=c(cdates[1],cdates[2]),encd4=c(50,100),idade=c(20,30),sex='F',enweight=c(25,50),enhgb=c(8,9),education=1)
  mymodb <- summary(lrm1,myenroll=c(cdates[1],cdates[3]),encd4=c(50,200),idade=c(20,40),enweight=c(25,75))
  mymodc <- summary(lrm1,myenroll=c(cdates[1],cdates[4]),encd4=c(50,350),idade=c(20,50))
# CREATE LaTeX TABLE OF MODEL EFFECTS AND 95% CI
  lrm1table <- rep(NA,8)
  lrm1table[1] <- paste("Female",      myformat3(mymoda[22,4]),myformat3(mymoda[22,6]),myformat3(mymoda[22,7]),myformat4(getanova$P[2]),sep=" & ")
  lrm1table[2] <- paste("20  (ref) ", 1 , " " , " ", myformat4(getanova$P[3]),sep=" & ")
  lrm1table[3] <- paste("30 ",myformat3(mymoda[2,4]), myformat3(mymoda[2,6]), myformat3(mymoda[2,7]),sep=" & ")
  lrm1table[4] <- paste("40 ",myformat3(mymodb[2,4]), myformat3(mymodb[2,6]), myformat3(mymodb[2,7]),sep=" & ")
  lrm1table[5] <- "DELETE"
  lrm1table[6] <- paste("Namacurra (ref) ",1 ," " ," ",myformat4(getanova$P[1]),sep=" & ")
  lrm1table[7] <- paste("Alto Molocue",myformat3(mymoda[12,4]), myformat3(mymoda[12,6]), myformat3(mymoda[12,7]),sep=" & ")
  lrm1table[8] <- paste("Gile",        myformat3(mymoda[14,4]), myformat3(mymoda[14,6]), myformat3(mymoda[14,7]),sep=" & ")
  lrm1table[9] <- paste("Ile",         myformat3(mymoda[16,4]), myformat3(mymoda[16,6]), myformat3(mymoda[16,7]),sep=" & ")
  lrm1table[10] <- paste("Inhassunge",  myformat3(mymoda[18,4]), myformat3(mymoda[18,6]), myformat3(mymoda[18,7]),sep=" & ")
  lrm1table[11] <- paste("Lugela",      myformat3(mymoda[20,4]), myformat3(mymoda[20,6]), myformat3(mymoda[20,7]),sep=" & ")
  lrm1table[12] <- paste("Illiterate (ref)", 1 , sep=" & ")
  lrm1table[13] <- paste("Primary School",   myformat3(mymoda[24,4]), myformat3(mymoda[24,6]), myformat3(mymoda[24,7]),myformat4(getanova$P[4]),sep=" & ")
  lrm1table[14] <- paste("Secondary School", myformat3(mymoda[26,4]), myformat3(mymoda[26,6]), myformat3(mymoda[26,7]),sep=" & ")
  lrm1table[15] <- paste("Technical School", myformat3(mymoda[28,4]), myformat3(mymoda[28,6]), myformat3(mymoda[28,7]),sep=" & ")
  lrm1table[16] <- paste("University Level", myformat3(mymoda[30,4]), myformat3(mymoda[30,6]), myformat3(mymoda[30,7]),sep=" & ")
  lrm1table[17] <- paste("50 vs. 25 ",        myformat3(mymoda[6,4]), myformat3(mymoda[6,6]), myformat3(mymoda[6,7]), myformat4(getanova$P[6]),sep=" & ")
  lrm1table[19] <- paste("75 vs. 25",        myformat3(mymodb[6,4]), myformat3(mymodb[6,6]), myformat3(mymodb[6,7]),sep=" & ")
  lrm1table[21] <- paste("100 vs. 50",myformat3(mymoda[4,4]), myformat3(mymoda[4,6]), myformat3(mymoda[4,7]), myformat4(getanova$P[5]),sep=" & ")
  lrm1table[22] <- paste("200 vs. 50",myformat3(mymodb[4,4]), myformat3(mymodb[4,6]), myformat3(mymodb[4,7]),sep=" & ")
  lrm1table[23] <- paste("350 vs. 50",myformat3(mymodc[4,4]), myformat3(mymodc[4,6]), myformat3(mymodc[4,7]),sep=" & ")
  lrm1table[24] <- paste("Hemoglobin (per g/dL)", myformat3(mymoda[8,4]), myformat3(mymoda[8,6]), myformat3(mymoda[8,7]),myformat4(getanova$P[7]),sep=" & ")
  lrm1table[25] <- paste(" I (ref)", 1 , " ", " ", myformat4(getanova$P[8]),sep=" & ")
  lrm1table[26] <- paste(" II",  myformat3(mymoda[32,4]), myformat3(mymoda[32,6]), myformat3(mymoda[32,7]),sep=" & ")
  lrm1table[27] <- paste(" III", myformat3(mymoda[34,4]), myformat3(mymoda[34,6]), myformat3(mymoda[34,7]),sep=" & ")
  lrm1table[28] <- paste(" IV",  myformat3(mymoda[36,4]), myformat3(mymoda[36,6]), myformat3(mymoda[36,7]),sep=" & ")
  lrm1table[29] <- paste("2006 (ref)", 1 , " ", " ", myformat4(getanova$P[9]),sep=" & ")
  lrm1table[30] <- paste("2007", myformat3(mymoda[10,4]), myformat3(mymoda[10,6]), myformat3(mymoda[10,7]),sep=" & ")
  lrm1table[31] <- paste("2008", myformat3(mymodb[10,4]), myformat3(mymodb[10,6]), myformat3(mymodb[10,7]),sep=" & ")
  lrm1table[32] <- paste("2009", myformat3(mymodc[10,4]), myformat3(mymodc[10,6]), myformat3(mymodc[10,7]),sep=" & ")
# get nonlinear p-value for hemoglobin
  lrm.hgb1 <- fit.mult.impute(txin90 ~  rcs(enhgb,4), lrm, lrm1.mi, lrm1.d1)
  lrm.hgb1.pval <- myformat3(as.data.frame(anova(lrm.hgb1))$P[2])
@

\newpage
\subsection{Baseline Characteristics as Predictors of ART Initiation}
We use a logistic regression (with multiple imputation) to identify whether baseline demographics,
 lab results, and clinical assessment are predictive of starting therapy.  
 The outcome is ART initiation within the first 90 days of enrollment.  Please reference table \ref{tab:allsites90}
 for summary statistics of this subpopulation.
 Please note that both CD4 and hemoglobin have 
 >50\% missing data, but they are included as known predictors of ART initiation.

Some continuous predictors are nonlinear with log-odds of ART initiation; these are, age, CD4 count,
 weight, and date of enrollment.
 To relax our linearity assumption, we model each one using a 4 knot restricted cubic spline function.  
 Predictive ability may be quantified using the $c$ index. A value of $c$ of 0.5 indicates random
 predictions, while a value of 1 indicates perfect prediction$^1$.  This model has good predictive
 ability ($c=\Sexpr{format(round(lrm1.stats[6],2),nsmall=2)}$).
\\

\footnotetext[1]{Harrell, FE. {\it Regression Modeling Strategies}. New York, NY: Springer; 2001.}

\vspace{-.15in}

\begin{center}
\begin{threeparttable}
\footnotesize
\caption{Model Effects: ART entry in first 90 days}
\centering{
\begin{tabular}{llllll}
\hline
  & Odds Ratio  & Lower 95\% & Upper 95\% & P-value \\
\hline 
\\
\Sexpr{lrm1table[1]} \\
Age (yrs)\\
\hspace{.2in}\Sexpr{lrm1table[2]} \\
\hspace{.2in}\Sexpr{lrm1table[3]} \\
\hspace{.2in}\Sexpr{lrm1table[4]} \\
District \\
\hspace{.2in}\Sexpr{lrm1table[6]} \\
\hspace{.2in}\Sexpr{lrm1table[7]} \\
\hspace{.2in}\Sexpr{lrm1table[8]} \\
\hspace{.2in}\Sexpr{lrm1table[9]} \\
\hspace{.2in}\Sexpr{lrm1table[10]} \\
\hspace{.2in}\Sexpr{lrm1table[11]} \\
Education \\
\hspace{.2in}\Sexpr{lrm1table[12]} \\
\hspace{.2in}\Sexpr{lrm1table[13]} \\
\hspace{.2in}\Sexpr{lrm1table[14]} \\
\hspace{.2in}\Sexpr{lrm1table[15]} \\
\hspace{.2in}\Sexpr{lrm1table[16]} \\
Weight (kg)\\
\hspace{.2in}\Sexpr{lrm1table[17]} \\
\hspace{.2in}\Sexpr{lrm1table[19]} \\
CD4 count (cells/$\mu$L) \\
\hspace{.2in}\Sexpr{lrm1table[21]} \\
\hspace{.2in}\Sexpr{lrm1table[22]} \\
\hspace{.2in}\Sexpr{lrm1table[23]} \\
% Hemoglobin
\Sexpr{lrm1table[24]} \\
WHO Stage \\
\hspace{.2in}\Sexpr{lrm1table[25]} \\
\hspace{.2in}\Sexpr{lrm1table[26]} \\
\hspace{.2in}\Sexpr{lrm1table[27]} \\
\hspace{.2in}\Sexpr{lrm1table[28]} \\
Year of Enrollment \\
\hspace{.2in}\Sexpr{lrm1table[29]} \\
\hspace{.2in}\Sexpr{lrm1table[30]} \\
\hspace{.2in}\Sexpr{lrm1table[31]} \\
\hspace{.2in}\Sexpr{lrm1table[32]} \\
\\
\hline
\end{tabular}
}
\begin{tablenotes}
\item There is strong evidence that the association between age and the log-odds of ART initiation is non-linear (p\Sexpr{myformat3(getnonlrm1[1,3])}).
\item There is strong evidence that the association between CD4 count and the log-odds of ART initiation is non-linear (p\Sexpr{myformat3(getnonlrm1[2,3])}).
\item There is some evidence that the association between weight and the log-odds of ART initiation is non-linear (p=\Sexpr{myformat3(getnonlrm1[3,3])}).
\item There is strong evidence that the association between date of ART initiation and the log-odds of ART initiation is non-linear (p\Sexpr{myformat3(getnonlrm1[4,3])}).
\item There is no evidence that the association between hemoglobin and the log-odds of ART initiation is non-linear (p=\Sexpr{lrm.hgb1.pval}).
\end{tablenotes}
\end{threeparttable}
\label{tab:lrm1}
\end{center}
\newpage



<<echo = false, results = hide>>=
# RUN LOGISTIC REGRESSION MODEL WITH ART ENTRY AS OUTCOME FOR PARTICIPANTS ELIGIBLE FOR ART - USE MULTIPLE IMPUTATION FOR MISSING DATA
# PREDICTORS OF ART INITIATOIN AMONG THOSE PATIENTS ELIGIBLE TO BEGIN TREATMENT
# If a patient had a CD4 count <200/μL independent of WHO Stage, was in WHO Stage IV independent of CD4 count,
# or had a CD4 count between 200-349/μL and was documented in WHO Stage III, they were recommended for treatment
# initiation with the first line regimen of lamivudine (3TC), stavudine (d4T), and nevirapine (NVP). 
  elig1 <- ifelse(mymoz$encd4<200 & !is.na(mymoz$encd4),1,0)
  elig2 <- ifelse(mymoz$enwho==4 & !is.na(mymoz$enwho),1,0)
  elig3 <- ifelse(mymoz$encd4>=200 & mymoz$encd4<350 & !is.na(mymoz$encd4) & mymoz$enwho==3 & !is.na(mymoz$enwho),1,0)
  eligmoz <- mymoz[rowSums(data.frame(elig1,elig2,elig3))>0,]
  eligenwho <- as.factor(eligmoz$enwho)
  eligstart  <- as.numeric(as.Date(eligmoz$dataabertura, "%m/%d/%Y")) - 13000
  eligmoz$education[eligmoz$education==6]<-4
  eligeducation <- as.factor(eligmoz$education)
# SET VARIABLES IN DATA FRAME
  lrm1.eligd1 <- data.frame(eligmoz$txin90,eligmoz$coddistrito,eligmoz$sexo,
                               eligmoz$idade,eligmoz$encd4,eligmoz$enweight,
                               eligstart,eligmoz$enhgb,eligenwho,eligeducation)
  colnames(lrm1.eligd1) <- c("txin90","site","sexo","idade","encd4","enweight","eligstart","enhgb","eligenwho","education")
  lrm1.eligdd1 <- datadist(lrm1.eligd1)
# IMPUTE MISSING CD4, HGB, WEIGHT and WHO STAGE
  set.seed(1000)
  lrm1.eligmi <- aregImpute(~txin90+site+sexo+idade+education+encd4+enweight+enhgb+eligenwho+eligstart,n.impute=10,data=lrm1.eligd1)
  options(datadist='lrm1.eligdd1')
# FIT LOGISTIC REGRESSION WITH DEFAULT N=10 IMPUTATIONS
  lrmelig1 <- fit.mult.impute(txin90 ~ site + sexo + rcs(idade,4) + education +  rcs(sqrt(encd4),4) +  
                        rcs(enweight,4) + enhgb + eligenwho + rcs(eligstart,4), lrm, lrm1.eligmi, lrm1.eligd1)
  (getanovaelig <- as.data.frame(anova(lrmelig1),row.names=c("site","sexo","idade","nonlin1","education","encd4",
             "nonlin2","enweight","nonlin3","enhgb","eligenwho","eligstart","nonlin4","totnon","total"))[c(1:3,5:6,8,10:12),])
  (getnonlrmelig1 <- as.data.frame(anova(lrmelig1),row.names=c("site","sexo","idade","age nonlinear","education","encd4",
             "CD4 nonlinear","enweight","Weight nonlinear","enhgb","eligenwho","eligstart","Date of Initiation nonlinear","totnon","total"))[c(4,7,9,13),])
  lrmelig1.stats <- lrmelig1$stats
  cdates <- as.numeric(as.Date(c("01/01/2006","01/01/2007","01/01/2008","01/01/2009"), "%m/%d/%Y")) - 13000
  eligmoda <- summary(lrmelig1,eligstart=c(cdates[1],cdates[2]),encd4=c(50,100),idade=c(20,30),sex='F',enweight=c(25,50),enhgb=c(8,9),education=1)
  eligmodb <- summary(lrmelig1,eligstart=c(cdates[1],cdates[3]),encd4=c(50,200),idade=c(20,40),enweight=c(25,75))
  eligmodc <- summary(lrmelig1,eligstart=c(cdates[1],cdates[4]),encd4=c(50,350),idade=c(20,50))
# CREATE LaTeX TABLE OF MODEL EFFECTS AND 95% CI
  lrmelig1table <- rep(NA,8)
  lrmelig1table[1] <- paste("Female",      myformat3(eligmoda[22,4]),myformat3(eligmoda[22,6]),myformat3(eligmoda[22,7]),myformat4(getanovaelig$P[2]),sep=" & ")
  lrmelig1table[2] <- paste("20  (ref) ", 1 , " " , " ", myformat4(getanovaelig$P[3]),sep=" & ")
  lrmelig1table[3] <- paste("30 ",myformat3(eligmoda[2,4]), myformat3(eligmoda[2,6]), myformat3(eligmoda[2,7]),sep=" & ")
  lrmelig1table[4] <- paste("40 ",myformat3(eligmodb[2,4]), myformat3(eligmodb[2,6]), myformat3(eligmodb[2,7]),sep=" & ")
  lrmelig1table[6] <- paste("Namacurra (ref) ",1 ," " ," ",myformat4(getanovaelig$P[1]),sep=" & ")
  lrmelig1table[7] <- paste("Alto Molocue",myformat3(eligmoda[12,4]), myformat3(eligmoda[12,6]), myformat3(eligmoda[12,7]),sep=" & ")
  lrmelig1table[8] <- paste("Gile",        myformat3(eligmoda[14,4]), myformat3(eligmoda[14,6]), myformat3(eligmoda[14,7]),sep=" & ")
  lrmelig1table[9] <- paste("Ile",         myformat3(eligmoda[16,4]), myformat3(eligmoda[16,6]), myformat3(eligmoda[16,7]),sep=" & ")
  lrmelig1table[10] <- paste("Inhassunge",  myformat3(eligmoda[18,4]), myformat3(eligmoda[18,6]), myformat3(eligmoda[18,7]),sep=" & ")
  lrmelig1table[11] <- paste("Lugela",      myformat3(eligmoda[20,4]), myformat3(eligmoda[20,6]), myformat3(eligmoda[20,7]),sep=" & ")
  lrmelig1table[12] <- paste("Illiterate (ref)", 1 , sep=" & ")
  lrmelig1table[13] <- paste("Primary School",   myformat3(eligmoda[24,4]), myformat3(eligmoda[24,6]), myformat3(eligmoda[24,7]),myformat4(getanovaelig$P[4]),sep=" & ")
  lrmelig1table[14] <- paste("Secondary School", myformat3(eligmoda[26,4]), myformat3(eligmoda[26,6]), myformat3(eligmoda[26,7]),sep=" & ")
  lrmelig1table[15] <- paste("Technical/University", myformat3(eligmoda[28,4]), myformat3(eligmoda[28,6]), myformat3(eligmoda[28,7]),sep=" & ")
  lrmelig1table[17] <- paste("50 vs. 25 ",        myformat3(eligmoda[6,4]), myformat3(eligmoda[6,6]), myformat3(eligmoda[6,7]), myformat4(getanovaelig$P[6]),sep=" & ")
  lrmelig1table[19] <- paste("75 vs. 25",        myformat3(eligmodb[6,4]), myformat3(eligmodb[6,6]), myformat3(eligmodb[6,7]),sep=" & ")
  lrmelig1table[21] <- paste("100 vs. 50",myformat3(eligmoda[4,4]), myformat3(eligmoda[4,6]), myformat3(eligmoda[4,7]), myformat4(getanovaelig$P[5]),sep=" & ")
  lrmelig1table[22] <- paste("200 vs. 50",myformat3(eligmodb[4,4]), myformat3(eligmodb[4,6]), myformat3(eligmodb[4,7]),sep=" & ")
  lrmelig1table[23] <- paste("350 vs. 50",myformat3(eligmodc[4,4]), myformat3(eligmodc[4,6]), myformat3(eligmodc[4,7]),sep=" & ")
  lrmelig1table[24] <- paste("Hemoglobin (per g/dL)", myformat3(eligmoda[8,4]), myformat3(eligmoda[8,6]), myformat3(eligmoda[8,7]),myformat4(getanovaelig$P[7]),sep=" & ")
  lrmelig1table[25] <- paste(" I (ref)", 1 , " ", " ", myformat4(getanovaelig$P[8]),sep=" & ")
  lrmelig1table[26] <- paste(" II",  myformat3(eligmoda[30,4]), myformat3(eligmoda[30,6]), myformat3(eligmoda[30,7]),sep=" & ")
  lrmelig1table[27] <- paste(" III", myformat3(eligmoda[32,4]), myformat3(eligmoda[32,6]), myformat3(eligmoda[32,7]),sep=" & ")
  lrmelig1table[28] <- paste(" IV",  myformat3(eligmoda[34,4]), myformat3(eligmoda[34,6]), myformat3(eligmoda[34,7]),sep=" & ")
  lrmelig1table[29] <- paste("2006 (ref)", 1 , " ", " ", myformat4(getanovaelig$P[9]),sep=" & ")              
  lrmelig1table[30] <- paste("2007", myformat3(eligmoda[10,4]), myformat3(eligmoda[10,6]), myformat3(eligmoda[10,7]),sep=" & ")
  lrmelig1table[31] <- paste("2008", myformat3(eligmodb[10,4]), myformat3(eligmodb[10,6]), myformat3(eligmodb[10,7]),sep=" & ")
  lrmelig1table[32] <- paste("2009", myformat3(eligmodc[10,4]), myformat3(eligmodc[10,6]), myformat3(eligmodc[10,7]),sep=" & ")
# get nonlinear p-value for hemoglobin
  lrm.hgbelig1 <- fit.mult.impute(txin90 ~  rcs(enhgb,4), lrm, lrm1.eligmi, lrm1.eligd1)
  lrm.hgbelig1.pval <- myformat3(as.data.frame(anova(lrm.hgbelig1))$P[2])
@

\newpage
\subsection{Baseline Characteristics as Predictors of ART Initiation for Patients Recommended for Treatment Initiation}
 If a patient had a CD4 count <200/$\mu$L independent of WHO Stage, was in WHO Stage IV independent of CD4 count,
 or had a CD4 count between 200-349/$\mu$L and was documented in WHO Stage III, they were recommended for treatment
 initiation with the first line regimen of lamivudine (3TC), stavudine (d4T), and nevirapine (NVP). 
 There are \Sexpr{nrow(eligmoz)} patients who met the criteria for ART initiation; \Sexpr{sum(eligmoz$txin90)} individuals
 initiated ART within the first 90 days of enrollment.
 We use a logistic regression (with multiple imputation) to identify whether baseline demographics,
 lab results, and clinical assessment are predictive of starting therapy among patients eligible for ART initiation.  
 The outcome is ART initiation within the first 90 days of enrollment.
 Please reference table \ref{tab:allsites90elig} for summary statistics of this subpopulation.
 Please note that hemoglobin has
 >50\% missing data, but it is included as known predictor of ART initiation.

Some continuous predictors are nonlinear with log-odds of ART initiation; these are, age, CD4 count,
 weight, and date of enrollment.
 To relax our linearity assumption, we model each one using a 4 knot restricted cubic spline function.  
 Predictive ability may be quantified using the $c$ index. A value of $c$ of 0.5 indicates random
 predictions, while a value of 1 indicates perfect prediction$^1$.  This model has good predictive
 ability ($c=\Sexpr{format(round(lrmelig1.stats[6],2),nsmall=2)}$).\\

\footnotetext[1]{Harrell, FE. {\it Regression Modeling Strategies}. New York, NY: Springer; 2001.}

\vspace{-.15in}

\begin{center}
\begin{threeparttable}
\footnotesize
\caption{Model Effects: ART entry in first 90 days among eligible patients}
\centering{
\begin{tabular}{llllll}
\hline
  & Odds Ratio  & Lower 95\% & Upper 95\% & P-value \\
\hline 
\\
\Sexpr{lrmelig1table[1]} \\
Age (yrs)\\
\hspace{.2in}\Sexpr{lrmelig1table[2]} \\
\hspace{.2in}\Sexpr{lrmelig1table[3]} \\
\hspace{.2in}\Sexpr{lrmelig1table[4]} \\
District \\
\hspace{.2in}\Sexpr{lrmelig1table[6]} \\
\hspace{.2in}\Sexpr{lrmelig1table[7]} \\
\hspace{.2in}\Sexpr{lrmelig1table[8]} \\
\hspace{.2in}\Sexpr{lrmelig1table[9]} \\
\hspace{.2in}\Sexpr{lrmelig1table[10]} \\
\hspace{.2in}\Sexpr{lrmelig1table[11]} \\
Education \\
\hspace{.2in}\Sexpr{lrmelig1table[12]} \\
\hspace{.2in}\Sexpr{lrmelig1table[13]} \\
\hspace{.2in}\Sexpr{lrmelig1table[14]} \\
\hspace{.2in}\Sexpr{lrmelig1table[15]} \\
Weight (kg)\\
\hspace{.2in}\Sexpr{lrmelig1table[17]} \\
\hspace{.2in}\Sexpr{lrmelig1table[19]} \\
CD4 count (cells/$\mu$L) \\
\hspace{.2in}\Sexpr{lrmelig1table[21]} \\
\hspace{.2in}\Sexpr{lrmelig1table[22]} \\
\hspace{.2in}\Sexpr{lrmelig1table[23]} \\
% Hemoglobin
\Sexpr{lrmelig1table[24]} \\
WHO Stage \\
\hspace{.2in}\Sexpr{lrmelig1table[25]} \\
\hspace{.2in}\Sexpr{lrmelig1table[26]} \\
\hspace{.2in}\Sexpr{lrmelig1table[27]} \\
\hspace{.2in}\Sexpr{lrmelig1table[28]} \\
Year of Enrollment \\
\hspace{.2in}\Sexpr{lrmelig1table[29]} \\
\hspace{.2in}\Sexpr{lrmelig1table[30]} \\
\hspace{.2in}\Sexpr{lrmelig1table[31]} \\
\hspace{.2in}\Sexpr{lrmelig1table[32]} \\
\\
\hline
\end{tabular}
}
\begin{tablenotes}
\item There is strong evidence that the association between age and the log-odds of ART initiation is non-linear (p\Sexpr{myformat3(getnonlrmelig1[1,3])}).
\item There is strong evidence that the association between CD4 count and the log-odds of ART initiation is non-linear (p\Sexpr{myformat3(getnonlrmelig1[2,3])}).
\item There is strong evidence that the association between date of ART initiation and the log-odds of ART initiation is non-linear (p\Sexpr{myformat3(getnonlrmelig1[4,3])}).
\item There is no evidence that the association between weight and the log-odds of ART initiation is non-linear (p=\Sexpr{myformat3(getnonlrmelig1[3,3])}).
\item There is no evidence that the association between hemoglobin and the log-odds of ART initiation is non-linear (p=\Sexpr{lrm.hgbelig1.pval}).
\end{tablenotes}
\end{threeparttable}
\label{tab:lrmelig1}
\end{center}
\newpage

<<echo = false, results = hide>>=
set.seed(10)

# CREATE MYENROLL
  myenroll  <- as.numeric(as.Date(dataabertura, "%m/%d/%Y")) - 13000
  cdates <- as.numeric(as.Date(c("07/01/2006","07/01/2007","07/01/2008","07/01/2009"), "%m/%d/%Y")) - 13000
#########################
#  LOGISTIC REGRESSION  #
#   SEX BY ENROLLMENT   #
#########################
# SET VARIABLES IN DATA FRAME
  newsexo <- ifelse(sexo=='F',1,0)
  lrm2.d2 <- data.frame(newsexo,myenroll)
  lrm2.dd2 <- datadist(lrm2.d2)
  options(datadist='lrm2.dd2')
# There are zero missing covariates
  lrm2 <- lrm(newsexo ~ rcs(myenroll,5), x=TRUE, y=TRUE, data=lrm2.d2) 
  (getanova2 <- as.data.frame(anova(lrm2)))
# MODEL ESTIMATES FOR FEMALE ENROLLMENT
  mymod2a <- summary(lrm2,myenroll=c(cdates[1],cdates[2]))
  mymod2b <- summary(lrm2,myenroll=c(cdates[1],cdates[3]))
  mymod2c <- summary(lrm2,myenroll=c(cdates[1],cdates[4]))
# LaTeX TABLE FOR FEMALE ENROLLMENT
  lrm2table <- rep(NA,4)
  lrm2table[1] <- paste("2006 (ref)", 1 , " ", " ", myformat4(getanova2$P[1]),sep=" & ")
  lrm2table[2] <- paste("2007", myformat3(mymod2a[2,4]), myformat3(mymod2a[2,6]), myformat3(mymod2a[2,7]),sep=" & ")
  lrm2table[3] <- paste("2008", myformat3(mymod2b[2,4]), myformat3(mymod2b[2,6]), myformat3(mymod2b[2,7]),sep=" & ")
  lrm2table[4] <- paste("2009", myformat3(mymod2c[2,4]), myformat3(mymod2c[2,6]), myformat3(mymod2c[2,7]),sep=" & ")

#############################
#  PROPORTIONAL ODDS MODEL  #
#  WHO STAGE BY ENROLLMENT  #
#############################
# SET VARIABLES IN DATA FRAME
  myenwho3 <- as.factor(enwho)
  lrm3.d3 <- data.frame(myenwho3,myenroll)
  lrm3.dd3 <- datadist(lrm3.d3)
  options(datadist='lrm3.dd3')
# There are missing covariates
  lrm3 <- lrm(myenwho3 ~ rcs(myenroll,5), x=TRUE, y=TRUE, data=lrm3.d3) 
  (getanova3 <- as.data.frame(anova(lrm3)))
# MODEL ESTIMATES FOR FEMALE ENROLLMENT
  mymod3a <- summary(lrm3,myenroll=c(cdates[1],cdates[2]))
  mymod3b <- summary(lrm3,myenroll=c(cdates[1],cdates[3]))
  mymod3c <- summary(lrm3,myenroll=c(cdates[1],cdates[4]))
# LaTeX TABLE FOR WHO STAGE ENROLLMENT
  lrm3table <- rep(NA,4)
  lrm3table[1] <- paste("2006 (ref)", 1 , " ", " ", myformat4(getanova3$P[1]),sep=" & ")
  lrm3table[2] <- paste("2007", myformat3(mymod3a[2,4]), myformat3(mymod3a[2,6]), myformat3(mymod3a[2,7]),sep=" & ")
  lrm3table[3] <- paste("2008", myformat3(mymod3b[2,4]), myformat3(mymod3b[2,6]), myformat3(mymod3b[2,7]),sep=" & ")
  lrm3table[4] <- paste("2009", myformat3(mymod3c[2,4]), myformat3(mymod3c[2,6]), myformat3(mymod3c[2,7]),sep=" & ")
@
\clearpage 
\subsection{Enrollment Trends}
This section summarizes trends in enrollment (regardless of ART initiation) by sex
 and WHO staging to the HIV Care and Treatment program on or before 31 March 2009. 
 The relationship between sex and enrollment is modeled using a logistic regression
 with sex as the outcome and date of enrollment as a predictor.  
 Similarly, the relationship between staging and enrollment is modeled using proportional odds
 (ordinal logistic regression) with stage as the outcome and date of enrollment as a
 predictor. 
 For both models, date of enrollment is a continuous predictor
 that is nonlinear with the log-odds. To relax our linearity assumption, we model date of enrollment with
 a 5 knot restricted cubic spline function. 

\begin{center}
\begin{threeparttable}
\footnotesize
\caption{Model Effects: Sex enrollment trend for HIV Care and Treatment Programs}
\centering{
\begin{tabular}{llllll}
\hline
  & Odds Ratio  & Lower 95\% & Upper 95\% & P-value \\
\hline 
\\
Year of Enrollment (F vs. M)\\
\hspace{.2in}\Sexpr{lrm2table[1]} \\
\hspace{.2in}\Sexpr{lrm2table[2]} \\
\hspace{.2in}\Sexpr{lrm2table[3]} \\
\hspace{.2in}\Sexpr{lrm2table[4]} \\
\\
\hline
\end{tabular}
}
\end{threeparttable}
\label{tab:lrm2}
\end{center}

\vspace{-.25in}

<<echo = false, results = hide>>=
whonew <- quarter <- rep(NA,length(dataabertura))
for(i in 1:length(dataabertura)){
  if(as.Date("01/01/2006", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("10/01/2006", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 1}  # 01/1/2006-3/31/2006 - NOTE: only 3 ppts before 7/1/06
  if(as.Date("10/01/2006", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("01/01/2007", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 2}  # 10/1/2006-3/31/2007
  if(as.Date("01/01/2007", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("04/01/2007", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 3}  # 01/1/2007-3/31/2007
  if(as.Date("04/01/2007", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("07/01/2007", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 4}  # 03/1/2007-6/30/2007
  if(as.Date("07/01/2007", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("10/01/2007", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 5}  # 07/1/2007-9/30/2007
  if(as.Date("10/01/2007", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("01/01/2008", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 6}  # 10/1/2007-1/31/2008
  if(as.Date("01/01/2008", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("04/01/2008", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 7}  # 01/1/2008-3/31/2008
  if(as.Date("04/01/2008", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("07/01/2008", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 8}  # 04/1/2008-6/30/2008
  if(as.Date("07/01/2008", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("10/01/2008", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 9}  # 07/1/2008-9/30/2008
  if(as.Date("10/01/2008", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("01/01/2009", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 10} # 10/1/2009-1/31/2009
  if(as.Date("01/01/2009", "%m/%d/%Y") <= as.Date(dataabertura[i], "%m/%d/%Y") & as.Date("04/01/2009", "%m/%d/%Y") > as.Date(dataabertura[i], "%m/%d/%Y")){quarter[i] <- 11} # 01/1/2009-3/31/2009
}

qdates <- format(as.Date(c("06/01/2006","10/01/2006",
		    "01/01/2007","04/01/2007","07/01/2007","10/01/2007",
                    "01/01/2008","04/01/2008","07/01/2008","10/01/2008","01/01/2009"), "%m/%d/%Y"), "%m/%d/%Y")
qmmyy <- format(as.Date(c("06/01/2006","10/01/2006",
		    "01/01/2007","04/01/2007","07/01/2007","10/01/2007",
                    "01/01/2008","04/01/2008","07/01/2008","10/01/2008","01/01/2009"), "%m/%d/%Y"), "%m/%y")
qlabel <- c("Q1","Q2","Q3","Q4","Q5","Q6","Q7","Q8","Q9","Q10","Q11")
whonew[enwho==1 | enwho==2] <- 1
whonew[enwho==3 | enwho==4] <- 2
propwho <- table(whonew,quarter)[2,]/(table(whonew,quarter)[1,] + table(whonew,quarter)[2,]) # STAGES III/IV  /  STAGES I/II
propwho1 <- table(enwho,quarter)[1,]/(table(enwho,quarter)[1,] + table(enwho,quarter)[2,] + table(enwho,quarter)[3,] + table(enwho,quarter)[4,]) # STAGES I  /  STAGES I/II/III/IV
propwho2 <- table(enwho,quarter)[2,]/(table(enwho,quarter)[1,] + table(enwho,quarter)[2,] + table(enwho,quarter)[3,] + table(enwho,quarter)[4,]) # STAGES I  /  STAGES I/II/III/IV
propwho3 <- table(enwho,quarter)[3,]/(table(enwho,quarter)[1,] + table(enwho,quarter)[2,] + table(enwho,quarter)[3,] + table(enwho,quarter)[4,]) # STAGES I  /  STAGES I/II/III/IV
propwho4 <- table(enwho,quarter)[4,]/(table(enwho,quarter)[1,] + table(enwho,quarter)[2,] + table(enwho,quarter)[3,] + table(enwho,quarter)[4,]) # STAGES I  /  STAGES I/II/III/IV
@

\setkeys{Gin}{width=.55\textwidth} %default is 0.7
\begin{figure}[h!]
   \caption{Proportion of Female Enrollees by Quarter}
   \begin{center}
      \label{fig:rawfem}
<<fig=TRUE, echo=FALSE>>=
propsexo <- table(sexo,quarter)[1,]/(table(sexo,quarter)[1,] + table(sexo,quarter)[2,]) # Proportion FEMALE 
plot(1:11,propsexo,axes=F,col=0,ylim=c(0.5,0.8),ylab="Proportion of Enrollees",xlab="Quarter of Enrollment")
box()
axis(side=2)
axis(side=1,at=c(1:11),labels=qlabel)
axis(side=3,at=c(1:11),labels=qmmyy,tick=F)
abline(h=0.5,lty=2)
lines(1:11,propsexo,lwd=1.5)
points(1:11,propsexo,cex=0.7)
legend(1,.8,"Female",col=1,lty=1,lwd=1.5,bty="n",cex=.8)
@
  \end{center}
\end{figure} 


\clearpage
\begin{center}
%\setlength{\tabcolsep}{4pt}
\begin{threeparttable}
\footnotesize
\caption{Model Effects: WHO staging enrollment trend for HIV Care and Treatment Programs}
\centering{
\begin{tabular}{llllll}
\hline
  & Odds Ratio  & Lower 95\% & Upper 95\% & P-value \\
\hline 
\\
Year of Enrollment (Higher vs. Lower Stage)\\
\hspace{.2in}\Sexpr{lrm3table[1]} \\
\hspace{.2in}\Sexpr{lrm3table[2]} \\
\hspace{.2in}\Sexpr{lrm3table[3]} \\
\hspace{.2in}\Sexpr{lrm3table[4]} \\
\\
\hline
\end{tabular}
}
\end{threeparttable}
\label{tab:lrm3}
\end{center}

\vspace{-.25in}



\setkeys{Gin}{width=.55\textwidth} %default is 0.7
\begin{figure}[h!]
   \caption{Proportion of WHO Stages by Quarter}
   \begin{center}
      \label{fig:rawstage}
<<fig=TRUE, echo=FALSE>>=
plot(1:11,propwho,axes=F,col=0,ylim=c(0,0.5),ylab="Proportion of Enrollees",xlab="Quarter of Enrollment")
box()
axis(side=2)
axis(side=1,at=c(1:11),labels=qlabel)
axis(side=3,at=c(1:11),labels=qmmyy,tick=F)
lines(1:11,propwho1,lty=1,col=gray(.4), lwd=1.5)
points(1:11,propwho1,cex=0.7,col=gray(.4))
lines(1:11,propwho2,lty=2,col=gray(.4), lwd=1.5)
points(1:11,propwho2,cex=0.7,col=gray(.4))
lines(1:11,propwho3,lty=3)
points(1:11,propwho3,cex=0.7)
lines(1:11,propwho4,lty=4)
points(1:11,propwho4,cex=0.7)
abline(h=0.5,lty=2)
abline(h=0.25,lty=6)
legend(1,0.08,c("Stage I","Stage II","Stage III","Stage IV"),
       col=c(gray(.4),gray(.4),1,1),lty=c(1,2,3,4),
       lwd=c(1.5,1.5,1,1),bty="n",cex=.8)
@
  \end{center}
\end{figure} 




<<echo = false, results = hide>>=
# RUN LOGISTIC REGRESSION MODEL WITH MISSING DATA AS OUTCOME
# SET VARIABLES IN DATA FRAME
  miss.encd4 <- ifelse(is.na(encd4),0,1)
  miss.enhgb <- ifelse(is.na(enhgb),0,1)
  lrm.miss.d <- data.frame(miss.encd4,miss.enhgb,myenroll)
  lrm.miss.dd <- datadist(lrm.miss.d)
  options(datadist='lrm.miss.dd')
#######################
# LOGISTIC REGRESSION #
#  CD4 BY ENROLLMENT  #
#######################
  lrm.miss.cd4 <- lrm(miss.encd4 ~ rcs(myenroll,5), x=TRUE, y=TRUE, data=lrm.miss.d) 
  (miss.cd4 <- as.data.frame(anova(lrm.miss.cd4)))
  lrm.miss.cd4.stats <- lrm.miss.cd4$stats
# MODEL ESTIMATES FOR MISSING CD4 DATA
  mymod2a <- summary(lrm.miss.cd4,myenroll=c(cdates[1],cdates[2]))
  mymod2b <- summary(lrm.miss.cd4,myenroll=c(cdates[1],cdates[3]))
  mymod2c <- summary(lrm.miss.cd4,myenroll=c(cdates[1],cdates[4]))
# LaTeX TABLE FOR MISSING CD4 DATA BY ENROLLMENT
  misscd4.lrmtable <- rep(NA,4)
  misscd4.lrmtable[1] <- paste("2006 (ref)", 1 , " ", " ", myformat4(miss.cd4$P[1]),sep=" & ")
  misscd4.lrmtable[2] <- paste("2007", myformat3(mymod2a[2,4]), myformat3(mymod2a[2,6]), myformat3(mymod2a[2,7]),sep=" & ")
  misscd4.lrmtable[3] <- paste("2008", myformat3(mymod2b[2,4]), myformat3(mymod2b[2,6]), myformat3(mymod2b[2,7]),sep=" & ")
  misscd4.lrmtable[4] <- paste("2009", myformat3(mymod2c[2,4]), myformat3(mymod2c[2,6]), myformat3(mymod2c[2,7]),sep=" & ")
#######################
# LOGISTIC REGRESSION #
#  HGB BY ENROLLMENT  #
#######################
  lrm.miss.hgb <- lrm(miss.enhgb ~ rcs(myenroll,5), x=TRUE, y=TRUE, data=lrm.miss.d) 
  (miss.hgb <- as.data.frame(anova(lrm.miss.hgb)))
# MODEL ESTIMATES FOR MISSING HGB DATA
  mymod2a <- summary(lrm.miss.hgb,myenroll=c(cdates[1],cdates[2]))
  mymod2b <- summary(lrm.miss.hgb,myenroll=c(cdates[1],cdates[3]))
  mymod2c <- summary(lrm.miss.hgb,myenroll=c(cdates[1],cdates[4]))
# LaTeX TABLE FOR MISSING HGB DATA BY ENROLLMENT
  misshgb.lrmtable <- rep(NA,4)
  misshgb.lrmtable[1] <- paste("2006 (ref)", 1 , " ", " ", myformat4(miss.hgb$P[1]),sep=" & ")
  misshgb.lrmtable[2] <- paste("2007", myformat3(mymod2a[2,4]), myformat3(mymod2a[2,6]), myformat3(mymod2a[2,7]),sep=" & ")
  misshgb.lrmtable[3] <- paste("2008", myformat3(mymod2b[2,4]), myformat3(mymod2b[2,6]), myformat3(mymod2b[2,7]),sep=" & ")
  misshgb.lrmtable[4] <- paste("2009", myformat3(mymod2c[2,4]), myformat3(mymod2c[2,6]), myformat3(mymod2c[2,7]),sep=" & ")
#######################
#  GRAPHIC - CD4 AND  #
#    HGB BY QUARTER   #
#######################
propmiss.encd4 <- table(miss.encd4,quarter)[1,]/(table(miss.encd4,quarter)[1,] + table(miss.encd4,quarter)[2,]) # Proportionn NOT MISSING 
propmiss.enhgb <- table(miss.enhgb,quarter)[1,]/(table(miss.enhgb,quarter)[1,] + table(miss.enhgb,quarter)[2,]) # Proportion NOT MISSING 

@
\clearpage
\subsection{Data Collection Trends}
This section summarizes trends in data collection for all enrollees (regardless of ART initiation)
 to the HIV Care and Treatment program on or before 31 March 2009. 
 The relationship between presence of CD4 or hemoglobin records and enrollment is
 modeled using a logistic regression
 with missingness as the outcome and date of enrollment as a predictor.  
 For both models, date of enrollment is a continuous predictor
 that is nonlinear with the log-odds. To relax our linearity assumption, we model date of enrollment with
 a 5 knot restricted cubic spline function. \\


\begin{center}
\begin{threeparttable}
\footnotesize
\caption{Model Effects: CD4 data collection trends for HIV Care and Treatment Programs}
\centering{
\begin{tabular}{llllll}
\hline
  & Odds Ratio  & Lower 95\% & Upper 95\% & P-value \\
\hline 
\\
Year of Enrollment (Present vs. Missing)\\
\hspace{.2in}\Sexpr{misscd4.lrmtable[1]} \\
\hspace{.2in}\Sexpr{misscd4.lrmtable[2]} \\
\hspace{.2in}\Sexpr{misscd4.lrmtable[3]} \\
\hspace{.2in}\Sexpr{misscd4.lrmtable[4]} \\
\\
\hline
\end{tabular}
}
\end{threeparttable}
\label{tab:lrm2}
\end{center}

\vspace{-.1in}


\begin{center}
\begin{threeparttable}
\footnotesize
\caption{Model Effects: Hemoglobin data collection trends for HIV Care and Treatment Programs}
\centering{
\begin{tabular}{llllll}
\hline
  & Odds Ratio  & Lower 95\% & Upper 95\% & P-value \\
\hline 
\\
Year of Enrollment (Present vs. Missing)\\
\hspace{.2in}\Sexpr{misshgb.lrmtable[1]} \\
\hspace{.2in}\Sexpr{misshgb.lrmtable[2]} \\
\hspace{.2in}\Sexpr{misshgb.lrmtable[3]} \\
\hspace{.2in}\Sexpr{misshgb.lrmtable[4]} \\
\\
\hline
\end{tabular}
}
\end{threeparttable}
\label{tab:lrm2}
\end{center}


\vspace{-.25in}


\setkeys{Gin}{width=.45\textwidth} %default is 0.7
\begin{figure}[h!]
   \caption{Proportion of Participants with Lab Record by Quarter}
\vspace{-.15in}
   \begin{center}
      \label{fig:rawcd4}
<<fig=TRUE, echo=FALSE>>=
plot(1:11,propmiss.encd4,axes=F,col=0,ylim=c(0,0.61),ylab="Proportion of Enrollees",xlab="Quarter of Enrollment")
box()
axis(side=2)
axis(side=1,at=c(1:11),labels=qlabel)
axis(side=3,at=c(1:11),labels=qmmyy,tick=F)
lines(1:11,1-propmiss.encd4,lwd=1.5)
points(1:11,1-propmiss.encd4,cex=0.7)
lines(1:11,1-propmiss.enhgb,lwd=1.5,col=gray(0.4),lty=2)
points(1:11,1-propmiss.enhgb,cex=0.7)
legend(8,.05,c("CD4","Hemoglobin"),col=c(1,gray(0.4)),lty=c(1,2),lwd=1.5,bty="n",cex=.8)
@
  \end{center}
\end{figure} 



\end{document}