%#############################################################
%#
%#   Program: analysis-scripts-016.nw
%#	  
%#   Project: Comprehensive site assessment for the 
%#            International Epidemiologic Databases to 
%#            Evaluate AIDS (IeDEA) Collaboration
%# 
%#   Co-I:    SN Duda*, AM Farr, ML Lindegren, M Blevins, 
%#            CW Wester, R McKaig, K. Wools-Kaloustian, 
%#            D. Ekouevi, M. Egger, J. Hemingway-Foday, 
%#            D. Cooper, R. Moore, and D. Nash for the 
%#       * corresponding author, stephany.duda 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: 10 October 2012
%#
%#############################################################


\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[utf8]{inputenc}
\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}
\usepackage{float}
\usepackage[breaklinks,colorlinks=true]{hyperref}  %% ,citecolor=blue
\usepackage[all]{hypcap}
\newcommand{\squishlist}{
 \begin{list}{$\ast$}
  { \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}  }
\newcommand{\HRule}{\rule{\linewidth}{0.5mm}}
\title{Attitudes toward migration}

\begin{document}
\SweaveOpts{concordance=TRUE}
\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]{graphics/biostat-logo-mb.png}\\[2cm]

    \textsc{\Large Analysis Supplement:}\\[0.5cm]
 
    % Title
    \HRule \\[0.4cm]
    { \Large \bfseries IeDEAS Site Assessment:\\[0.8cm]Comprehensive HIV Care and Treatment}\\[0.4cm]
    \HRule \\[1.5cm]
   
    % Biostatistician and Faculty Member
\begin{minipage}{0.4\textwidth}
\begin{center} \large
\emph{Biostatistician:}\\
Meridith Blevins, M.S.
\end{center}
\end{minipage}

\vfill
 
% Bottom of the page
{\large \today}
 
\end{center}
 

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

\clearpage

% Hidden R code chunk 
<<echo = false, results = hide>>=
rm(list=ls()) # SHOULD BE FIRST LINE IN ALL R PROGRAMS - CLEARS NAMESPACE
setwd("/home/blevinml/Projects/IeDEAS/SiteAssessment"); 
set.seed(1)

#Read Data
data <- read.csv("Data/IeDEA_SiteAssessment_v1_SectionsA-C.csv",na.strings = "")

## HDI CLASSIFICATION
data$hdi <- rep(NA,nrow(data))
data$hdi[data$country %in% c("USA","United States","Canada","Argentina","Brazil","Chile","Mexico","Peru","Australia","Japan","Republic of Korea","Korea","Malaysia","Singapore","Taiwan")] <- 3
data$hdi[data$country %in% c("South Africa","Botswana","Honduras","Cambodia","China","India","Indonesia","Philippines","Thailand")] <- 2
data$hdi[data$country %in% c("Mozambique","Malawi","Zimbabwe","Rwanda",
			     "Zambia","Haiti","Benin","Burkina Faso","C\xf4te d'Ivoire",
			     "Cote d’Ivoire","Gambia","Mali","Nigeria","Senegal","Burundi","Cameroon","DRC","Congo","Kenya","Tanzania","Uganda")] <- 1
data$hdi=factor(data$hdi,levels=c(1:3),ordered=TRUE)
levels(data$hdi)=c("UN HDI Low-income","UN HDI Middle-income","UN HDI High-income")

## PEPFAR IN 2008 CLASSIFICATION
data$pepfar08 <- rep(NA,nrow(data))
data$pepfar08[data$country %in% c("Argentina","Australia","Benin","Brazil","Burkina Faso","Burundi","Cameroon","Canada",
                                  "Chile","Gambia","Honduras","Japan","Malaysia","Mali","Mexico","Peru","Philippines",
                                  "Republic of Korea","Senegal","Singapore","Taiwan","USA",
                                  "Cambodia","China","DRC","India","Indonesia","Malawi","Thailand","Zimbabwe")] <- 2
data$pepfar08[data$country %in% c("Botswana","C\xf4te d'Ivoire","Haiti",
                                  "Kenya","Mozambique","Nigeria","Rwanda","South Africa","Tanzania",
                                  "Uganda","Zambia")] <- 1
data$pepfar08=factor(data$pepfar08,levels=c(1:2))
levels(data$pepfar08)=c("PEPFAR","No PEPFAR")

## CLINIC DATA 
data$type <- as.factor(ifelse(!is.na(data$a06a1),"Primary",ifelse(!is.na(data$a06a4),"Secondary",ifelse(!is.na(data$a06a7),"Tertiary",NA))))
data$public <- ifelse(data$type=="Primary",as.character(data$a06a1),ifelse(data$type=="Secondary",as.character(data$a06a4),ifelse(data$type=="Tertiary",as.character(data$a06a7),NA)))
data$adult <- ifelse(!is.na(data$a05a1) & is.na(data$a05a3) & is.na(data$a05a4),"Adult only",ifelse(!is.na(data$a05a3) | !is.na(data$a05a4),"Adults and Children",NA))
data$yearcat <- cut(data$a13,breaks=c(0,1990,1995,2000,2005,2012),right=FALSE)
data$urban <- rep(NA,nrow(data))
data$urban[(!is.na(data$a04a1) | !is.na(data$a04a2)) & is.na(data$a04a3) & is.na(data$a04a4) & is.na(data$a04a5)] <- 1
data$urban[is.na(data$a04a1) & is.na(data$a04a2) & (!is.na(data$a04a3) | !is.na(data$a04a4)) & is.na(data$a04a5)] <- 2
data$urban[is.na(data$urban)] <- 3
data$urban=factor(data$urban,levels=c(1:3))
levels(data$urban)=c("Urban","Rural","Mixed")

## COUNTS MATCH STEPHANY
## FORCE NUMBER will take a factor/character to a number without issuing the warnings
forcenumber <- function(var){
    options(warn = -1)
    neg <- ifelse(substr(var,1,1)=="-","-","")
    x <- as.numeric(paste0(neg,gsub("[^0-9.]","",as.character(var))))
    options(warn = 1)
    return(x)
}
data$PR1ONSITE <- apply(data.frame(data$a11a1>0,data$a11a2>0,data$a11a3>0,data$a11a4>0,data$a11a5>0,data$a11a6>0,forcenumber(data$a11a7)>0),1,sum,na.rm=TRUE) > 0
data$PR2ONSITE <- apply(data.frame(data$a11b1>0,data$a11b2>0,data$a11b3>0,data$a11b4>0,data$a11b5>0,data$a11b6>0,forcenumber(data$a11b7)>0),1,sum,na.rm=TRUE) > 0
data$PR3ONSITE <- apply(data.frame(data$a11c1>0,data$a11c2>0,data$a11c3>0,data$a11c4>0,data$a11c5>0,data$a11c6>0,forcenumber(data$a11c7)>0),1,sum,na.rm=TRUE) > 0
data$PR4ONSITE <- apply(data.frame(data$a11d1>0,data$a11d2>0,data$a11d3>0,data$a11d4>0,data$a11d5>0,data$a11d6>0,forcenumber(data$a11d7)>0),1,sum,na.rm=TRUE) > 0
data$PR5ONSITE <- apply(data.frame(data$a11e1>0,data$a11e2>0,data$a11e3>0,data$a11e4>0,data$a11e5>0,data$a11e6>0,forcenumber(data$a11e7)>0),1,sum,na.rm=TRUE) > 0
data$PR6ONSITE <- apply(data.frame(data$a11f1>0,data$a11f2>0,data$a11f3>0,data$a11f4>0,data$a11f5>0,data$a11f6>0,forcenumber(data$a11f7)>0),1,sum,na.rm=TRUE) > 0
data$PR7ONSITE <- apply(data.frame(data$a11g1>0,data$a11g2>0,data$a11g3>0,data$a11g4>0,data$a11g5>0,data$a11g6>0,forcenumber(data$a11g7)>0),1,sum,na.rm=TRUE) > 0
data$PR8ONSITE <- apply(data.frame(data$a11h1>0,data$a11h2>0,data$a11h3>0,data$a11h4>0,data$a11h5>0,data$a11h6>0,forcenumber(data$a11h7)>0),1,sum,na.rm=TRUE) > 0
data$PR9ONSITE <- apply(data.frame(data$a11i1>0,data$a11i2>0,data$a11i3>0,data$a11i4>0,data$a11i5>0,data$a11i6>0,forcenumber(data$a11i7)>0),1,sum,na.rm=TRUE) > 0
data$PR10ONSITE <- apply(data.frame(data$a11j1>0,data$a11j2>0,data$a11j3>0,data$a11j4>0,data$a11j5>0,data$a11j6>0,forcenumber(data$a11j7)>0),1,sum,na.rm=TRUE) > 0


library(Hmisc)

#(1) ART adherence
#B06 -> 1 (1-on-1 counseling), 6 (calendars, checklist), 10 (routine review of pickup).
#AND B07 -> 2-6 (provide services more than once)
B07_2to6 <- c("At ART initiation and infrequently thereafter","At ART initiation and more than once a month thereafter","At ART initiation and once every 1-3 months thereafter","At ART initiation and once every 4-6 months thereafter","Less frequently")
B07_SPEC <- c("AT ART INITIATION AND AT EVERY APPOINTMENT","At ARV initiation and every visit afterwards (scheduled every 3-4 months).","at every clinic visit","at every visit","every visit","Every visit","Every Visit")
data$comp1 <- (!is.na(data$b06a1) | !is.na(data$b06a6) | !is.na(data$b06a10)) & (data$b07a1 %in% B07_2to6 | data$b07a2 %in% B07_SPEC)
B06MISS <- apply(data.frame(data$b06a1,data$b06a2,data$b06a3,data$b06a4,data$b06a5,data$b06a6,data$b06a7,data$b06a8,data$b06a9,data$b06a10,data$b06a11,data$b06a12,data$b06a13,data$b06a14),1,FUN=function(x) sum(!is.na(x)))==0
B07MISS <- apply(data.frame(data$b07a1,data$b07a2),1,FUN=function(x) sum(!is.na(x)))==0
data$comp1[!data$comp1 & (B06MISS | B07MISS)] <- NA
label(data$comp1) <- "ART adherence"

#(2) Nutritional support
#B09 -> (1 (nutritional counseling) AND one of 2-8) OR Nutritionist on-site
NUTONSITE <- apply(data.frame(data$a11i1>0,data$a11i2>0,data$a11i3>0,data$a11i4>0,data$a11i5>0,data$a11i6>0,forcenumber(data$a11i7)>0),1,sum,na.rm=TRUE) > 0
data$comp2 <- (!is.na(data$b09a1) & (!is.na(data$b09a2) | !is.na(data$b09a3) | !is.na(data$b09a4) | !is.na(data$b09a5) | !is.na(data$b09a6) | !is.na(data$b09a7) | !is.na(data$b09a8))) | (NUTONSITE)
B09MISS <- apply(data.frame(data$b09a1,   data$b09a3,   data$b09a5,   data$b09a7,   data$b09a9,   data$b09a11,  data$b09a13,  data$b09a15,  
                 data$b09a2,   data$b09a4,   data$b09a6,   data$b09a8,   data$b09a10,  data$b09a12,  data$b09a14,  data$b09a16),1,FUN=function(x) sum(!is.na(x)))==0
A11MISS <- apply(data.frame(data$a11i1,data$a11i2,data$a11i3,data$a11i4,data$a11i5,data$a11i6,data$a11i7),1,FUN=function(x) sum(!is.na(x)))==0
data$comp2[!data$comp2 & (B09MISS | A11MISS)] <- NA
label(data$comp2) <- "Nutritional support"

#(3) PMTCT
#B59 -> 2,3,4, Backup B03 if B59 missing. if B59 =missing, then true if B03 = 1-4, maybe 5 (check answer in "other" field).
B59MISS <- is.na(data$b59a2) & is.na(data$b59a3) & is.na(data$b59a4) 
B03MISS <- apply(data.frame(data$b03a1,data$b03a2,data$b03a3,data$b03a4,data$b03a5,data$b03a6,data$b03a7),1,FUN=function(x) sum(!is.na(x)))==0
data$comp3 <- ifelse(!B59MISS,!is.na(data$b59a2) | !is.na(data$b59a3) | !is.na(data$b59a4),!is.na(data$b03a1) | !is.na(data$b03a2) | !is.na(data$b03a3) | !is.na(data$b03a4) | data$b03a6 %in% c("PMTCT/No ANC","PMTCT within ART  clinic","ARV for PMTCT","PMTCT with HIV"))
data$comp3[B03MISS & !data$comp3] <- NA
label(data$comp3) <- "PMTCT"

#(4) CD4 available on or off-site from B34
data$comp4 <- data$b34a12 %in% c("Off Site","On Site")
data$comp4[is.na(data$b34a12) & !data$comp4] <- NA
label(data$comp4) <- "CD4 testing"

#(5) TB. Screening 
#C04 -> YES, everyone screened 
data$comp5 <- data$c04=="YES, for everyone screened for TB"
data$comp5[is.na(data$c04) & !data$comp5] <- NA  ## DON'T REALLY NEED THIS LINE
label(data$comp5) <- "TB screening"

#(6) Prevention
#B04 -> 2 or more of 1-10 (=b04a1 - b04a11 in REDCap)
data$comp6 <- !is.na(data$b04a1) & apply(data.frame(data$b04a2,data$b04a3,data$b04a4,data$b04a5,data$b04a6,data$b04a7,data$b04a8,data$b04a9,data$b04a10,data$b04a11),1,FUN=function(x) sum(!is.na(x)))>0
B04MISS <- apply(data.frame(data$b04a1,data$b04a2,data$b04a3,data$b04a4,data$b04a5,data$b04a6,data$b04a7,data$b04a8,data$b04a9,data$b04a10,data$b04a11,data$b04a12,data$b04a13,data$b04a14,data$b04a15,data$b04a16),1,FUN=function(x) sum(!is.na(x)))==0
data$comp6[!is.na(data$pg09_b04_b05_complete) & data$pg09_b04_b05_complete=="Incomplete" & !data$comp6] <- NA
label(data$comp6) <- "Prevention"

#(7) Outreach
#B05=3 (only outreach post-ART)
data$comp7 <- data$b05a3=="Outreach program for patients who miss appointments" & !is.na(data$b05a3)
B05MISS <- apply(data.frame(data$b05a1,data$b05a2,data$b05a3,data$b05a4,data$b05a5,data$b05a6,data$b05a7),1,FUN=function(x) sum(!is.na(x)))==0
data$comp7[B05MISS & !data$comp7] <- NA
label(data$comp7) <- "Outreach"

data$compcomp <- apply(with(data,data.frame(comp1,comp2,comp3,comp4,comp5,comp6,comp7)),1,FUN=function(x) sum(!is.na(x)))
data$compmiss <- apply(with(data,data.frame(comp1,comp2,comp3,comp4,comp5,comp6,comp7)),1,FUN=function(x) sum(is.na(x)))
## CASEWISE DELETION
data$compsum1 <- apply(with(data,data.frame(comp1,comp2,comp3,comp4,comp5,comp6,comp7)),1,sum)
## ASSUME NO SERVICE
data$compsum3 <- apply(with(data,data.frame(comp1,comp2,comp3,comp4,comp5,comp6,comp7)),1,sum,na.rm=TRUE)
## PERCENT OF NON-MISSING
data$compsum2 <- round(100*data$compsum3/data$compcomp)  
## ASSUME SERVICE
data$compsum4 <- apply(with(data,data.frame(comp1,comp2,comp3,comp4,comp5,comp6,comp7)),1,FUN=function(x) sum(ifelse(is.na(x),1,x)))
## LOOK AT SUMMARY DERIVATIONS
with(data,data.frame(compsum1,compsum2,compsum3,compsum4,compmiss))


comprehensive <- with(data,data.frame(comp1,comp2,comp3,comp4,comp5,comp6,comp7))
library(stringr)

data$region <- str_trim(str_replace_all(data$region, "[^[:alpha:]]", " "))
library(lattice)
allcomps <- c(data$compsum1,data$compsum2,data$compsum3,data$compsum4)
allgrps <- c(rep("Complete case",nrow(data)),rep("Percent of nonmissing",nrow(data)),rep("Assume no service",nrow(data)),rep("Assume service",nrow(data)))


setwd("/home/blevinml/Projects/IeDEAS/SiteAssessment/Code")
myper1 <- function(x){
    y <- NULL
    if(x < 0.05 & x!=0){y <- "<0.1"}
    else{y <- format(round(x,1),nsmall=1)}
    return(y)
}
printline <- function(var,subset=TRUE){
if(sum(var[subset],na.rm=TRUE)>0){
  line <- paste(c("",label(var),table(!var[subset],useNA='always'),paste(myper1(100*sum(var[subset],na.rm=TRUE)/sum(!is.na(var[subset]))),"\\\\%",sep="")),collapse=" & ")
  }
if(sum(var[subset],na.rm=TRUE)==length(var[subset])){
  line <- paste(c("",label(var),length(var[subset]),"0",sum(is.na(var[subset])),paste(myper1(100*sum(var[subset],na.rm=TRUE)/sum(!is.na(var[subset]))),"\\\\%",sep="")),collapse=" & ")
  }
if(sum(var[subset],na.rm=TRUE)>0 & !any(!var[subset],na.rm=TRUE) & sum(var[subset],na.rm=TRUE)!=length(var[subset])){
  line <- paste(c("",label(var),length(var[subset]),"0",sum(is.na(var[subset])),paste(myper1(100*sum(var[subset],na.rm=TRUE)/sum(!is.na(var[subset]))),"\\\\%",sep="")),collapse=" & ")
  }
if(sum(!var[subset],na.rm=TRUE)==length(var[subset])){
  line <- paste(c("",label(var),"0",length(var[subset]),sum(is.na(var[subset])),paste(myper1(100*sum(var[subset],na.rm=TRUE)/sum(!is.na(var[subset]))),"\\\\%",sep="")),collapse=" & ")
  }
if(sum(var[subset],na.rm=TRUE)==0){
  line <- paste(c("",label(var),"-","-",sum(is.na(var[subset])),"-"),collapse=" & ")
  }
return(line)
}
printlinecc <- function(var){
line <- paste(c(label(var),table(!var),paste(myper1(100*sum(var,na.rm=TRUE)/sum(!is.na(var))),"\\\\%",sep="")),collapse=" & ")
if(all(var)){line <- paste(c(label(var),table(!var),0,paste(myper1(100*sum(var,na.rm=TRUE)/sum(!is.na(var))),"\\\\%",sep="")),collapse=" & ")}
return(line)
}
sumtable <- paste(printline(data$comp1),printline(data$comp2),printline(data$comp3),printline(data$comp4),printline(data$comp5),printline(data$comp6),printline(data$comp7),sep=" \\\\\\\\ ")
sumtable1 <- paste(printline(data$comp1,subset=(data$region=="NA ACCORD")),printline(data$comp2,subset=(data$region=="NA ACCORD")),printline(data$comp3,subset=(data$region=="NA ACCORD")),printline(data$comp4,subset=(data$region=="NA ACCORD")),printline(data$comp5,subset=(data$region=="NA ACCORD")),printline(data$comp6,subset=(data$region=="NA ACCORD")),printline(data$comp7,subset=(data$region=="NA ACCORD")),sep=" \\\\\\\\ ")
sumtable2 <- paste(printline(data$comp1,subset=(data$region=="CCASAnet")),printline(data$comp2,subset=(data$region=="CCASAnet")),printline(data$comp3,subset=(data$region=="CCASAnet")),printline(data$comp4,subset=(data$region=="CCASAnet")),printline(data$comp5,subset=(data$region=="CCASAnet")),printline(data$comp6,subset=(data$region=="CCASAnet")),printline(data$comp7,subset=(data$region=="CCASAnet")),sep=" \\\\\\\\ ")
sumtable3 <- paste(printline(data$comp1,subset=(data$region=="Asia Pacific")),printline(data$comp2,subset=(data$region=="Asia Pacific")),printline(data$comp3,subset=(data$region=="Asia Pacific")),printline(data$comp4,subset=(data$region=="Asia Pacific")),printline(data$comp5,subset=(data$region=="Asia Pacific")),printline(data$comp6,subset=(data$region=="Asia Pacific")),printline(data$comp7,subset=(data$region=="Asia Pacific")),sep=" \\\\\\\\ ")
sumtable4 <- paste(printline(data$comp1,subset=(data$region=="Central Africa")),printline(data$comp2,subset=(data$region=="Central Africa")),printline(data$comp3,subset=(data$region=="Central Africa")),printline(data$comp4,subset=(data$region=="Central Africa")),printline(data$comp5,subset=(data$region=="Central Africa")),printline(data$comp6,subset=(data$region=="Central Africa")),printline(data$comp7,subset=(data$region=="Central Africa")),sep=" \\\\\\\\ ")
sumtable5 <- paste(printline(data$comp1,subset=(data$region=="East Africa")),printline(data$comp2,subset=(data$region=="East Africa")),printline(data$comp3,subset=(data$region=="East Africa")),printline(data$comp4,subset=(data$region=="East Africa")),printline(data$comp5,subset=(data$region=="East Africa")),printline(data$comp6,subset=(data$region=="East Africa")),printline(data$comp7,subset=(data$region=="East Africa")),sep=" \\\\\\\\ ")
sumtable6 <- paste(printline(data$comp1,subset=(data$region=="Southern Africa")),printline(data$comp2,subset=(data$region=="Southern Africa")),printline(data$comp3,subset=(data$region=="Southern Africa")),printline(data$comp4,subset=(data$region=="Southern Africa")),printline(data$comp5,subset=(data$region=="Southern Africa")),printline(data$comp6,subset=(data$region=="Southern Africa")),printline(data$comp7,subset=(data$region=="Southern Africa")),sep=" \\\\\\\\ ")
sumtable7 <- paste(printline(data$comp1,subset=(data$region=="West Africa")),printline(data$comp2,subset=(data$region=="West Africa")),printline(data$comp3,subset=(data$region=="West Africa")),printline(data$comp4,subset=(data$region=="West Africa")),printline(data$comp5,subset=(data$region=="West Africa")),printline(data$comp6,subset=(data$region=="West Africa")),printline(data$comp7,subset=(data$region=="West Africa")),sep=" \\\\\\\\ ")
sumtablecc <- paste(printlinecc(data$comp1[!is.na(data$compsum1)]),printlinecc(data$comp2[!is.na(data$compsum1)]),printlinecc(data$comp3[!is.na(data$compsum1)]),printlinecc(data$comp4[!is.na(data$compsum1)]),printlinecc(data$comp5[!is.na(data$compsum1)]),printlinecc(data$comp6[!is.na(data$compsum1)]),printlinecc(data$comp7[!is.na(data$compsum1)]),sep=" \\\\\\\\ ")
#[1] public/private
#[2] primary/secondary/tertiary
#[3] adult-only vs. adult/peds
#[4] diversity of staff types


## WRITE COUNTRY LEVEL DATA FOR MAP
j <- 1
countryinfo <- matrix(NA,nrow=80,ncol=4)
for(i in unique(data$country)){
countryinfo[j,] <- c(i,length(data$compsum1[data$country==i & !is.na(data$compsum1)]),mean(data$compsum1[data$country==i & !is.na(data$compsum1)]),median(data$compsum1[data$country==i & !is.na(data$compsum1)]))
j <- j + 1
}  
countryinfo <- as.data.frame(countryinfo[!is.na(countryinfo[,1]),])
names(countryinfo) <- c("Country","Site Count","Mean of Comprehensive Count","Median of Comprehensive Count")
write.csv(countryinfo,"countryinfo.csv",row.names=FALSE)

j <- 1
regioninfo <- matrix(NA,nrow=80,ncol=4)
for(i in unique(data$region)){
regioninfo[j,] <- c(i,length(data$compsum1[data$region==i & !is.na(data$compsum1)]),mean(data$compsum1[data$region==i & !is.na(data$compsum1)]),median(data$compsum1[data$region==i & !is.na(data$compsum1)]))
j <- j + 1
}  
regioninfo <- as.data.frame(regioninfo[!is.na(regioninfo[,1]),])
names(regioninfo) <- c("Region","Site Count","Mean of Comprehensive Count","Median of Comprehensive Count")
write.csv(regioninfo,"regioninfo.csv",row.names=FALSE)

siteinfo <- with(data,data.frame(site_id,country,region,compsum1))
names(siteinfo) <- c("Site ID","Country","Region","Comprehensive Count")
write.csv(siteinfo,"siteinfo.csv",row.names=FALSE)
@

\section{Comprehensive Services}

\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Comprehensive Service Categories}
\label{tab:attlab}
\begin{tabular}{llp{8cm}}
\hline
  Number & Label & Survey Inputs \\
\hline 
 1 & ART adherence & B06: 1 (1-on-1 counseling), 6 (calendars, checklist), 10 (routine review of pickup). AND B07: 2-6 (provide services more than once)\\
 2 & Nutritional support & B09: (1 (nutritional counseling) AND one of 2-8) OR Nutritionist on-site\\
 3 & PMTCT & B59: 2,3,4. If B59 missing, then true if B03: 1-4, maybe 5 (check answer in "other" field).\\
 4 & CD4 testing & CD4 available on or off-site from B34\\
 5 & TB screening & C04: YES, everyone screened \\
 6 & Prevention & B04: 1 (HIV counseling and testing) AND B04: 1 or more of 2-10\\
 7 & Outreach & B05: 3 (only outreach post-ART) \\
%\\
\hline\hline
\end{tabular}
}
\end{threeparttable}
\end{center} 


\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Services offered at all sites}
\label{tab:attlab}
\begin{tabular}{lllllll}
\hline
  Region & Service & Offered & Not Offered & Missing & Percent offered  \\
  & &  &   &  & (of nonmissing) \\
\hline 
North America    \Sexpr{sumtable1}\\
\hline
Latin America    \Sexpr{sumtable2}\\
\hline
Asia-Pacific     \Sexpr{sumtable3}\\
\hline
Central Africa   \Sexpr{sumtable4}\\
\hline
East Africa      \Sexpr{sumtable5}\\
\hline
Southern Africa  \Sexpr{sumtable6}\\
\hline
West Africa      \Sexpr{sumtable7}\\
\hline
All Regions      \Sexpr{sumtable}\\
%\\
\hline\hline
\end{tabular}
}
\end{threeparttable}
\end{center} 


 \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Cumulative distribution of ART provision by region}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
plotcumsum1 <- data.frame(table(data$a13,data$region))
for(i in 1:7){
plotcumsum1$cumsum[plotcumsum1$Var2==unique(data$region)[i]] <- cumsum(plotcumsum1$Freq[plotcumsum1$Var2==unique(data$region)[i]])
}
year <- seq(1985,2010,by=5)
plotcumsum1$Var1 <- as.numeric(as.character(plotcumsum1$Var1))

#png("Code/graphics/cumsum_cohort_year.png",res=250,width=1200,height=1200, bg="transparent") 
p=xyplot(cumsum~Var1, plotcumsum1, lwd=1.5,type='l', group=Var2, auto.key = list(points = FALSE, lines = TRUE, columns = 2, cex=.8),scales=list(x=list(at=year, labels=year)),xlab="Year ART services began",ylab="Number of sites",aspect=1)
print(p)
#dev.off()
@
    \end{center}
  \end{figure} 


 \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Cumulative proportion of ART provision by region}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
for(i in 1:7){
plotcumsum1$total[plotcumsum1$Var2==unique(data$region)[i]] <- sum(plotcumsum1$Freq[plotcumsum1$Var2==unique(data$region)[i]])
}
plotcumsum1$cumprop <- 100*plotcumsum1$cumsum/plotcumsum1$total

#png("Code/graphics/cumsum_cohort_year.png",res=250,width=1200,height=1200, bg="transparent") 
p=xyplot(cumprop~Var1, plotcumsum1, lwd=1.5,type='l', group=Var2, auto.key = list(points = FALSE, lines = TRUE, columns = 2, cex=.8),scales=list(x=list(at=year, labels=year)),xlab="Year ART services began",ylab="Percentage of sites",aspect=1)
print(p)
#dev.off()
@
    \end{center}
  \end{figure} 


  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Comprehensive Services Histograms}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
 print(histogram(~allcomps | allgrps,scales=list(x=list(relation="free")),breaks=NULL,xlab="Comprehensive Services"))
@ 
\\
{\footnotesize There are 128 sites in the database and 93 have complete information (upper left) on comprehensive service offerings. If we take the percent of nonmissing (upper right), the distribution is similar
 to assuming the site offers the comprehensive service (lower right); whereas, assuming no service may overpenalize those sites that did not report the information (lower left).  
The remaining figures in this report will utilize only those 93 sites with complete information.}
    \end{center}
  \end{figure} 



\vspace{-.35in}
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Comprehensive Services by Year ART services began with Region}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
allyears <- rep(data$a13,4)
allregion <- rep(data$region,4)

print(xyplot(allcomps~allyears | allgrps, jitter.x=TRUE, jitter.y=TRUE,groups=allregion, scales=list(y=list(relation="free")),breaks=NULL,ylab="Comprehensive Services",xlab="Year ART services began",
 auto.key = list(points =TRUE, columns = 2)))
@
    \end{center}
  \end{figure}
 
% Hidden R code chunk 
<<echo = false, results = hide>>=
# 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 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))}
}
## 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)
}
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)
}


data$pephdi <- paste(data$pepfar08,": ",data$hdi,sep="")
data <- data[!is.na(data$compsum1),]
data$compsum1cat <- cut(data$compsum1,breaks=c(0,5,6,7))
levels(data$compsum1cat) <- c("Low (3-5)","Medium (6)","High (7)")
COMP.N <- printheader(table(data$compsum1cat),total.label="Overall")#,pval=TRUE)
COMP1.N <- printheader(table(data$compsum1),total.label="Overall")
comp1.n1 <- getnp1line(var=ifelse(data$comp1,1,NA),by=data$compsum1cat,rowlabel=label(data$comp1),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
comp2.n1 <- getnp1line(var=ifelse(data$comp2,1,NA),by=data$compsum1cat,rowlabel=label(data$comp2),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
comp3.n1 <- getnp1line(var=ifelse(data$comp3,1,NA),by=data$compsum1cat,rowlabel=label(data$comp3),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
comp4.n1 <- getnp1line(var=ifelse(data$comp4,1,NA),by=data$compsum1cat,rowlabel=label(data$comp4),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
comp5.n1 <- getnp1line(var=ifelse(data$comp5,1,NA),by=data$compsum1cat,rowlabel=label(data$comp5),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
comp6.n1 <- getnp1line(var=ifelse(data$comp6,1,NA),by=data$compsum1cat,rowlabel=label(data$comp6),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
comp7.n1 <- getnp1line(var=ifelse(data$comp7,1,NA),by=data$compsum1cat,rowlabel=label(data$comp7),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)

att1.n1 <- getnptex(var=data$region,by=data$compsum1cat,label="Region",inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
att2.n1 <- getnptex(var=data$urban,by=data$compsum1cat,label="Patient population",inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
att3.n1 <- getnptex(var=data$yearcat,by=data$compsum1cat,label="Year of ART Provision",inc.pval=FALSE,by.order=TRUE,rowlabel=c("1985-1989 ","1990-1994 ","1995-1999 ","2000-2004 ","2005-2009 "),rowper=TRUE)
att4.n1 <- getnptex(var=data$public,by=data$compsum1cat,label="Type of Facility",inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
att5.n1 <- getnptex(var=data$type,by=data$compsum1cat,label="Level of Facility",inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
att6.n1 <- getnptex(var=data$adult,by=data$compsum1cat,label="Patients seen in clinic",inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
att7.n1 <- getnptex(var=data$pepfar08,by=data$compsum1cat,label="PEPFAR Country (2008)",inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
att8.n1 <- getnptex(var=data$hdi,by=data$compsum1cat,label="HDI Income Category (2010)",inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
att9.n1 <- getnptex(var=data$pephdi,by=data$compsum1cat,label="PEPFAR (2008) and HDI (2010)",inc.pval=FALSE,by.order=TRUE,rowper=TRUE)


ONSITE1.n1 <- getnp1line(var=ifelse(data$PR1ONSITE ,1,NA),by=data$compsum1cat,rowlabel="Physicians",indent=TRUE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
ONSITE2.n1 <- getnp1line(var=ifelse(data$PR2ONSITE ,1,NA),by=data$compsum1cat,rowlabel="Pediatricians",indent=TRUE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
ONSITE3.n1 <- getnp1line(var=ifelse(data$PR3ONSITE ,1,NA),by=data$compsum1cat,rowlabel="Mid-level providers",indent=TRUE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
ONSITE4.n1 <- getnp1line(var=ifelse(data$PR4ONSITE ,1,NA),by=data$compsum1cat,rowlabel="Nurses/Midwives",indent=TRUE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
ONSITE5.n1 <- getnp1line(var=ifelse(data$PR5ONSITE ,1,NA),by=data$compsum1cat,rowlabel="Nursing assistants",indent=TRUE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
ONSITE6.n1 <- getnp1line(var=ifelse(data$PR6ONSITE ,1,NA),by=data$compsum1cat,rowlabel="Lay health workers, adherence counselors or outreach workers",indent=TRUE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
ONSITE7.n1 <- getnp1line(var=ifelse(data$PR7ONSITE ,1,NA),by=data$compsum1cat,rowlabel="Pharmacists",indent=TRUE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
ONSITE8.n1 <- getnp1line(var=ifelse(data$PR8ONSITE ,1,NA),by=data$compsum1cat,rowlabel="Pharmacy assistants",indent=TRUE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
ONSITE9.n1 <- getnp1line(var=ifelse(data$PR9ONSITE ,1,NA),by=data$compsum1cat,rowlabel="Nutritionists",indent=TRUE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)
ONSITE10.n1 <- getnp1line(var=ifelse(data$PR10ONSITE,1,NA),by=data$compsum1cat,rowlabel="Data capturers",indent=TRUE,inc.pval=FALSE,by.order=TRUE,rowper=TRUE)


table1 <- paste(comp1.n1,comp2.n1,comp3.n1,comp4.n1,comp5.n1,comp6.n1,comp7.n1,sep=" \\\\\\\\ ")
table2 <- paste(att1.n1,att2.n1,att3.n1$line.incmiss,att4.n1$line.incmiss,att5.n1$line.incmiss,att6.n1,att7.n1,att8.n1,att9.n1,"Number of sites with provider category on staff",ONSITE1.n1,ONSITE2.n1 ,ONSITE3.n1 ,ONSITE4.n1 ,ONSITE5.n1 ,ONSITE6.n1 ,ONSITE7.n1 ,ONSITE8.n1 ,ONSITE9.n1 ,ONSITE10.n1,sep=" \\\\\\\\ ")


## GET COUNTS FOR ADHERENCE
table(data$b06a2[!(data$comp1)])
table(data$b06a3[!(data$comp1)])
table(data$b06a4[!(data$comp1)])
table(data$b06a5[!(data$comp1)])
table(data$b06a7[!(data$comp1)])
table(data$b06a8[!(data$comp1)])
table(data$b06a9[!(data$comp1)])
table(data$b06a11[!(data$comp1)])
table(data$b06a12[!(data$comp1)])
table(data$b06a13[!(data$comp1)])
table(data$b06a1[!(data$comp1)])
table(data$b06a6[!(data$comp1)])
table(data$b06a10[!(data$comp1)])

data$b07a1[(!data$comp1 & !is.na(data$b06a1))]
data$b07a2[(!data$comp1 & !is.na(data$b06a1))]

adherencetext <- "Among those sites without comprehensive ART adherence support, other offerings included: group counseling (8), appointment slips (8), education material (6),
                  pill boxes (5), educational videotape (1), alarm clocks (1), pharmacist on multidisciplinary team (7). There are sites that offered comprehensive services, but
                   not at regular intervals: 1-on-1 counseling (10), calendars (4), and routine review of pick-up (7)."



HDI.N <- printheader(table(data$hdi),total.label="Overall")#,pval=TRUE)
comp1.n2 <- getnp1line(var=ifelse(data$comp1,1,NA),by=data$hdi,rowlabel=label(data$comp1),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=FALSE)
comp2.n2 <- getnp1line(var=ifelse(data$comp2,1,NA),by=data$hdi,rowlabel=label(data$comp2),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=FALSE)
comp3.n2 <- getnp1line(var=ifelse(data$comp3,1,NA),by=data$hdi,rowlabel=label(data$comp3),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=FALSE)
comp4.n2 <- getnp1line(var=ifelse(data$comp4,1,NA),by=data$hdi,rowlabel=label(data$comp4),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=FALSE)
comp5.n2 <- getnp1line(var=ifelse(data$comp5,1,NA),by=data$hdi,rowlabel=label(data$comp5),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=FALSE)
comp6.n2 <- getnp1line(var=ifelse(data$comp6,1,NA),by=data$hdi,rowlabel=label(data$comp6),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=FALSE)
comp7.n2 <- getnp1line(var=ifelse(data$comp7,1,NA),by=data$hdi,rowlabel=label(data$comp7),indent=FALSE,inc.pval=FALSE,by.order=TRUE,rowper=FALSE)
table3 <- paste(comp1.n2,comp2.n2,comp3.n2,comp4.n2,comp5.n2,comp6.n2,comp7.n2,sep=" \\\\\\\\ ")
@
\section{Comprehensive Services: Complete Case only}

\begin{comment}
\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Service Summary (Complete Cases, n=93)}
\label{tab:attlab}
\begin{tabular}{lllll}
\hline
  Service & Offered & Not Offered & Percent offered  \\
\hline 
\Sexpr{sumtablecc}\\
\hline\hline
\end{tabular}
}
\end{threeparttable}
\end{center} 
\end{comment}

\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Distribution of comprehensive services}
\label{tab:attlab}
\begin{tabular}{llllllll}
\hline
 & \multicolumn{4}{c}{Number of Services}\\
 \Sexpr{COMP.N}  \\
\hline 
\Sexpr{table1}\\
\hline\hline
\end{tabular}
}
\end{threeparttable}
\end{center} 

\Sexpr{adherencetext}


\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Distribution of comprehensive services by UND HDI Income}
\label{tab:attlab}
\begin{tabular}{llllllll}
\hline
 & \multicolumn{4}{c}{UN HDI Income}\\
 \Sexpr{HDI.N}  \\
\hline 
\Sexpr{table3}\\
\hline\hline
\end{tabular}
}
	  \begin{tablenotes}
	  \item[a] Column percentages in order to understand which essential services are lacking among high-income countries.  Outreach, TB screening, and nutritional support appear to have the lowest percentages.
	  \end{tablenotes}
\end{threeparttable}
\end{center} 

\Sexpr{adherencetext}
  
  
\vspace{-.35in}
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Comprehensive Services by Year ART services began with Region}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
print(xyplot(data$compsum1~data$a13, jitter.x=TRUE, jitter.y=TRUE,groups=data$region, scales=list(x=list(relation="free")),breaks=NULL,ylab="Comprehensive Services",xlab="Year ART services began",
 auto.key = list(points =TRUE, columns = 2)))
@
    \end{center}
  \end{figure} 


  \setkeys{Gin}{width=.99\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Distribution of Comprehensive Services by Site (Complete Cases, n=93)}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, height=7, width=14, echo=FALSE>>=
compsubset <- with(data[order(data$compsum1),],c(comp1,comp2,comp3,comp4,comp5,comp6,comp7))
labels <- c(sapply(label(comprehensive),rep,nrow(data)))
sitenum <- as.factor(rep(1:nrow(data),7))
yield <- rep(data$compsum1[order(data$compsum1)],7)
sitenum1 <- sitenum[!is.na(yield) & compsubset & !is.na(compsubset)]
labels1 <- labels[!is.na(yield) & compsubset & !is.na(compsubset)]
yield1 <- yield[!is.na(yield) & compsubset & !is.na(compsubset)]/7
yield2 <- rep(1,length(sitenum1))
print(barchart(yield2~sitenum1, groups=labels1, ylab="Service offered",xlab="Sites",stack=TRUE,auto.key=list(columns=4),scales=list(x=list(rot=45,cex=0.7))))
@
This figure allows us to assess the combination of services offered by each site.
    \end{center}
  \end{figure} 


\begin{center}
\setlength{\tabcolsep}{3pt}
\begin{threeparttable}
\footnotesize
\centering{
\caption{Facility characteristics}
\label{tab:attlab}
\begin{tabular}{p{6cm}lllllll}
\hline
 & \multicolumn{4}{c}{Number of Services}\\
 \Sexpr{COMP.N}  \\
\hline 
\Sexpr{table2}\\
%\\
\hline\hline
\end{tabular}
}
\end{threeparttable}
\end{center} 


  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Percentage of Sites with Comprehensive Services by Year ART services began}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
getcompper <- function(var,by,subset=TRUE){
  if(all(by)){x <- rep(100,length(unique(by)))}
  if(!any(by,na.rm=TRUE)){x <- rep(0,length(unique(by)))}
  if(any(by,na.rm=TRUE) & !all(by)){x <- 100*table(var[subset],by[subset])[,2]/table(var[subset & !is.na(by)])}
  return(x)
}

allcomp3 <- getcompper(data$yearcat,data$compsum1<4)
allcomp4 <- getcompper(data$yearcat,data$compsum1==4)
allcomp5 <- getcompper(data$yearcat,data$compsum1==5)
allcomp6 <- getcompper(data$yearcat,data$compsum1==6)
allcomp7 <- getcompper(data$yearcat,data$compsum1>6)
allyears <- rep(1:5,5)
labels1 <- c(rep("3",5),rep("4",5),rep("5",5),rep("6",5),rep("7",5))
#tiff(filename = "graphics/panel1.tiff",res = 300,height=1200,width=1200) #, pointsize = 12)
print(barchart(c(allcomp3,allcomp4,allcomp5,allcomp6,allcomp7) ~ as.factor(allyears), groups=labels1, stack=TRUE,auto.key=list(columns=2), scales=list(x=list(at=1:5, labels=c("1985-\n1989 ","1990-\n1994 ","1995-\n1999 ","2000-\n2004 ","2005-\n2009 "))),type="l",breaks=NULL,xlab="Year ART services began",ylab="Percentage"))
#dev.off()
@
    \end{center}
  \end{figure} 


\clearpage
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Comprehensive Services by Year ART services began with Patient population}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
print(xyplot(data$compsum1~data$a13, jitter.x=TRUE, jitter.y=TRUE, groups=data$urban, ylab="Comprehensive Services",xlab="Year ART services began",type=c("g","p"),
 auto.key = list(points =TRUE, columns = 1)))
@
    \end{center}
  \end{figure} 
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Distribution of Patient population by Comprehensive Services}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
allpublic3 <- getcompper(data$compsum1,data$urban=="Urban")
allpublic4 <- getcompper(data$compsum1,data$urban=="Rural")
allpublic5 <- getcompper(data$compsum1,data$urban=="Mixed")
allyears <- rep(3:7,3)
labels1=factor(c(rep(1,5),rep(2,5),rep(3,5)),levels=c(1:3),ordered=TRUE)
levels(labels1)=c("Urban","Rural","Mixed")
#tiff(filename = "graphics/panel2.tiff",res = 300,height=1200,width=1200) #, pointsize = 12)
print(barchart(c(allpublic3,allpublic4,allpublic5) ~ as.factor(allyears), groups=labels1, stack=TRUE,auto.key=list(columns=1), scales=list(x=list(at=1:5, labels=c(3:7))),type="l",breaks=NULL,xlab="Number of comprehensive services",ylab="Percentage"))
#dev.off()
@
    \end{center}
  \end{figure} 


\clearpage
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Comprehensive Services by Year ART services began with Patients seen in clinic}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
print(xyplot(data$compsum1~data$a13, jitter.x=TRUE, jitter.y=TRUE, groups=data$adult, ylab="Comprehensive Services",xlab="Year ART services began",type=c("g","p"), #,"smooth"),
 auto.key = list(points =TRUE, columns = 1)))
@
    \end{center}
  \end{figure} 
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Distribution of Patients seen in clinic by Comprehensive Services}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
allpublic3 <- getcompper(data$compsum1,data$adult=="Adult only")
allpublic4 <- getcompper(data$compsum1,data$adult=="Adults and Children")
allyears <- rep(3:7,2)
labels1 <- c(rep("Adult only",5),rep("Adults and Children",5))
print(barchart(c(allpublic3,allpublic4) ~ as.factor(allyears), groups=labels1, stack=TRUE,auto.key=list(columns=1), scales=list(x=list(at=1:5, labels=c(3:7))),type="l",breaks=NULL,xlab="Number of comprehensive services",ylab="Percentage"))
@
    \end{center}
  \end{figure} 



\clearpage
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Comprehensive Services by Year ART services began with Type of Facility}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
print(xyplot(data$compsum1~data$a13, jitter.x=TRUE, jitter.y=TRUE, groups=data$public, ylab="Comprehensive Services",xlab="Year ART services began",type=c("g","p"), #,"smooth"),
 auto.key = list(points =TRUE, columns = 1)))
@
    \end{center}
  \end{figure} 
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Distribution of Type of Facility by Comprehensive Services}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
allpublic3 <- getcompper(data$compsum1,data$public=="Public or Government")
allpublic4 <- getcompper(data$compsum1,data$public=="Private Clinic")
allyears <- rep(3:7,2)
labels1 <- c(rep("Public or Government",5),rep("Private Clinic",5))
print(barchart(c(allpublic3,allpublic4) ~ as.factor(allyears), groups=labels1, stack=TRUE,auto.key=list(columns=1), scales=list(x=list(at=1:5, labels=c(3:7))),type="l",breaks=NULL,xlab="Number of comprehensive services",ylab="Percentage"))
@
    \end{center}
  \end{figure} 


\clearpage
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Comprehensive Services by Year ART services began with Level of Facility}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
print(xyplot(data$compsum1~data$a13, jitter.x=TRUE, jitter.y=TRUE, groups=data$type, ylab="Comprehensive Services",xlab="Year ART services began",type=c("g","p"), #), #,"smooth"),
 auto.key = list(points =TRUE, columns = 1)))
@
    \end{center}
  \end{figure} 
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Distribution of Level of Facility by Comprehensive Services}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
alltype3 <- getcompper(data$compsum1,data$type=="Primary")
alltype4 <- getcompper(data$compsum1,data$type=="Secondary")
alltype5 <- getcompper(data$compsum1,data$type=="Tertiary")
allyears <- rep(3:7,3)
labels1 <- c(rep("Primary",5),rep("Secondary",5),rep("Tertiary",5))
#tiff(filename = "graphics/panel3.tiff",res = 300,height=1200,width=1200) #, pointsize = 12)
print(barchart(c(alltype3,alltype4,alltype5) ~ as.factor(allyears), groups=labels1, stack=TRUE,auto.key=list(columns=1), scales=list(x=list(at=1:5, labels=c(3:7))),type="l",breaks=NULL,xlab="Number of comprehensive services",ylab="Percentage"))
#dev.off()
@
    \end{center}
  \end{figure} 


\clearpage
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Comprehensive Services by Year ART services began with HDI}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=

print(xyplot(data$compsum1~data$a13, jitter.x=TRUE, jitter.y=TRUE, groups=data$hdi, ylab="Comprehensive Services",xlab="Year ART services began",type=c("g","p"), #,"smooth"),
 auto.key = list(points =TRUE, columns = 1)))
#cor.test(data$compsum1,data$a13,method="spearman")

@
    \end{center}
  \end{figure} 

  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Distribution of HDI Income by Comprehensive Services}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
allhdi3 <- getcompper(data$compsum1,data$hdi=="UN HDI Low-income")
allhdi4 <- getcompper(data$compsum1,data$hdi=="UN HDI Middle-income")
allhdi5 <- getcompper(data$compsum1,data$hdi=="UN HDI High-income")
allyears <- rep(3:7,3)

labels1=factor(c(rep(1,5),rep(2,5),rep(3,5)),levels=c(1:3),ordered=TRUE)
levels(labels1)=c("UN HDI Low-income","UN HDI Middle-income","UN HDI High-income")
#tiff(filename = "graphics/panel3.tiff",res = 300,height=1200,width=1200) #, pointsize = 12)
print(barchart(c(allhdi3,allhdi4,allhdi5) ~ as.factor(allyears), groups=labels1, stack=TRUE,auto.key=list(columns=1), scales=list(x=list(at=1:5, labels=c(3:7))),type="l",breaks=NULL,xlab="Number of comprehensive services",ylab="Percentage"))
#dev.off()
@
    \end{center}
  \end{figure} 


\clearpage
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Comprehensive Services by Year ART services began with PEPFAR country status (2008)}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
print(xyplot(data$compsum1~data$a13, jitter.x=TRUE, jitter.y=TRUE, groups=data$pepfar08, ylab="Comprehensive Services",xlab="Year ART services began",type=c("g","p"), #,"smooth"),
 auto.key = list(points =TRUE, columns = 1)))
@
    \end{center}
  \end{figure} 
  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Distribution of PEPFAR country status (2008) by Comprehensive Services}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=
allpublic3 <- getcompper(data$compsum1,data$pepfar08=="PEPFAR")
allpublic4 <- getcompper(data$compsum1,data$pepfar08=="No PEPFAR")
allyears <- rep(3:7,2)
labels1=factor(c(rep(1,5),rep(2,5)),levels=c(1:2),ordered=TRUE)
levels(labels1)=c("PEPFAR","No PEPFAR")
#tiff(filename = "graphics/panel4.tiff",res = 300,height=1200,width=1200) #, pointsize = 12)
print(barchart(c(allpublic3,allpublic4) ~ as.factor(allyears), groups=labels1, stack=TRUE,auto.key=list(columns=1), scales=list(x=list(at=1:5, labels=c(3:7))),type="l",breaks=NULL,xlab="Number of comprehensive services",ylab="Percentage"))
#dev.off()
@
    \end{center}
  \end{figure} 

  \setkeys{Gin}{width=.55\textwidth} %default is 0.7
  \begin{figure}[H]
    \caption{Distribution of PEPFAR country status (2008) by Comprehensive Services}
    \vspace{-0.2in}
    \begin{center}
	\label{fig:attsvc}
<<fig=TRUE, echo=FALSE>>=

allpephdi1 <- getcompper(data$compsum1,data$pephdi=="No PEPFAR: UN HDI Low-income"   )
allpephdi2 <- getcompper(data$compsum1,data$pephdi=="No PEPFAR: UN HDI Middle-income")
allpephdi3 <- getcompper(data$compsum1,data$pephdi=="No PEPFAR: UN HDI High-income"  )
allpephdi4 <- getcompper(data$compsum1,data$pephdi=="PEPFAR: UN HDI Low-income"      )
allyears <- rep(3:7,4)
labels1=factor(c(rep(1,5),rep(2,5),rep(3,5),rep(4,5)),levels=c(1:4),ordered=TRUE)
levels(labels1)=c("PEPFAR: UN HDI Low-income","No PEPFAR: UN HDI Low-income" ,"No PEPFAR: UN HDI Middle-income","No PEPFAR: UN HDI High-income" )
#tiff(filename = "graphics/panel5.tiff",res = 300,height=1200,width=1200) #, pointsize = 12)
print(barchart(c(allpephdi4,allpephdi1,allpephdi2,allpephdi3) ~ as.factor(allyears), groups=labels1, stack=TRUE,auto.key=list(columns=1, cex=.75), scales=list(x=list(at=1:5, labels=c(3:7))),type="l",breaks=NULL,xlab="Number of comprehensive services",ylab="Percentage"))
#dev.off()
@
    \end{center}
  \end{figure} 

\end{document}