%############################################################# %# %# Program: analysis-scripts-004.nw %# Project: Patient Characteristics, Enrollment Trends and %# Delayed Initiation of Combination Antiretroviral %# Therapy in Rural Northcentral Nigeria %# %# PI: Muktar Aliyu, MD PHD MPH* %# Amy Fenoglio %# * corresponding author, muktar.aliyu at vanderbilt.edu %# %# Biostatistician/Programmer: Meridith Blevins, MS %# ** contact programmer: meridith.blevins at vanderbilt.edu %# %# Purpose: Reproducible resarch. %# %# Notes: This '.nw' file is a Sweave file which allows %# for results from R statistical software to be %# nested within a LaTeX document. %# %# Published: 5 October 2012 %# %# Revisions: %# %# %# %############################################################# % RUN SWEAVE: % setwd("/home/blevinml/Projects/Nigeria/Code");Sweave("analysis-scripts-005.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[breaklinks,colorlinks=true]{hyperref} %% ,citecolor=blue \usepackage{rotating} \usepackage{tikz} \usetikzlibrary{shapes,shadows,arrows,trees} \usepackage{caption} \newcommand{\HRule}{\rule{\linewidth}{0.5mm}} \newcommand{\squishlist}{ \begin{list}{$\bullet$} { \setlength{\itemsep}{0pt} \setlength{\parsep}{3pt} \setlength{\topsep}{3pt} \setlength{\partopsep}{0pt} \setlength{\leftmargin}{1.5em} \setlength{\labelwidth}{1em} \setlength{\labelsep}{0.5em} } } \newcommand{\squishlisttwo}{ \begin{list}{$\bullet$} { \setlength{\itemsep}{0pt} \setlength{\parsep}{0pt} \setlength{\topsep}{0pt} \setlength{\partopsep}{0pt} \setlength{\leftmargin}{2em} \setlength{\labelwidth}{1.5em} \setlength{\labelsep}{0.5em} } } \newcommand{\squishend}{ \end{list} } \usepackage{Sweave} \begin{document} \SweaveOpts{prefix.string=graphics/plot} % Created a "graphics" subdirectory to save graph files in \begin{titlepage} \begin{center} % Upper part of the page \textsc{\Large Statistical Analysis Report:}\\[0.5cm] % Title \HRule \\[0.4cm] { \Large \bfseries ART initiation, enrollment and programmatic trends, mortality and retention among HIV infected patients\\[0.15cm] in rural Kwara and Niger states of Nigeria}\\[0.4cm] \HRule \\[1.5cm] % Biostatistician and Faculty Member \begin{minipage}{0.4\textwidth} \begin{flushleft} \large \emph{Biostatistician:}\\ Meridith Blevins, M.S. \end{flushleft} \end{minipage} \begin{minipage}{0.4\textwidth} \begin{flushright} \large \emph{Faculty Statistician:} \\ Bryan Shepherd, Ph.D. \end{flushright} \end{minipage} \vfill % Bottom of the page {\large \today} \end{center} \end{titlepage} \section{Introduction} The Vanderbilt University Institute for Global Health (VIGH) and its affiliate Friends in Global Health (FGH) currently supports HIV services in 13 sites in Kwara and Niger States. All 13 sites support HIV testing and counseling (HTC) services and 12 sites support prevention of mother-to-child transmission (PMTCT) of HIV services. Five of the sites provide HIV care and treatment services in addition to HTC and PMTCT services. \\ In a highly populated, diverse, and decentralized country like Nigeria, the Government of Nigeria, Presidents Emergency Plan for AIDS Relief (PEPFAR), and other programs have made enormous strides in addressing the HIV/AIDS epidemic, but much work remains. While the prevalence of HIV in Nigeria (estimated at 3.6\%) is relatively low for sub-Saharan Africa, due to its vast population, Nigeria has the second-largest number of persons living with HIV in the world, estimated at 3,300,000 people (UNAIDS Report on the Global AIDS Epidemic, 2010). Women are disproportionately affected by the epidemic, with a prevalence rate of 2.3\% (for ages 15-24 or higher) compared to young men at 0.8\%. As of 2008, PEPFAR programs have provided care and support to 1,043,100 Nigerians, only about 40\% of those in need (The President's Emergency Plan for AIDS Relief: Sixth Annual Report to Congress).\\ While PEPFAR has made a solid start in addressing HIV/AIDS in many cities, services are not accessible to all in need, particularly in rural areas. The ability to access care and treatment services beyond larger cities is very limited, underscoring a clear need for significant expansion of services. As of 2009, ARV coverage rate in Nigeria was only 21\% with 302,973 individuals receiving treatment, and 1,400,000 in need (WHO/UNAIDS/UNICEF, Towards Universal Access: Scaling Up Priority HIV/AIDS Interventions in the Health Sector, September 2010). With the planned analysis, we will use routinely collected patient monitoring data to describe the program and patient outcomes. \\ \subsection{Research Aims} \vspace{-0.07in} \subsubsection{To describe patients enrolled in HIV care and treatment.} \vspace{-0.07in} \subsubsection{To describe patients who initiate ART among those eligible for ART.} \vspace{-0.07in} \subsubsection{To identify enrollment trends in sex and WHO staging.} \vspace{-0.07in} \subsubsection{To identify missing data trends for CD4 and hemoglobin.} \vspace{-0.07in} \subsubsection{To examine trends from enrollment or ART inintiation of CD4 and hemoglobin.} \vspace{-0.07in} \subsubsection{To estimate all-cause mortality and LTFU during the first year of care and treatment.} \vspace{-0.07in} \subsubsection{To correct mortality for patients on ART.} \section{Methods} \subsection{Participants} This analysis uses data collected through May 31, 2011 and include HIV-infected patients entering HIV care and treatment programs at aged 15 and older. Only those participants who initiated treatment will be included in the all-cause mortality analysis. \\ All HIV-infected clients (pre-ART and ART) enrolled into care at VU/FGH-supported clinics receive baseline labs (including CD4, hematology, chemistry), basic care kit (BKC) containing an insecticide treated net, water vessel, and information, educational and communication IEC materials, along with monthly refills of water guard, condoms, and soap. Distribution of the BCK monthly refills encouraged retention into the program. Pre-ART and ART clients are also screened for TB and other opportunistic infections (OIs) and receive treatment as indicated. VU/FGH supports provision of cotrimoxazole prophylaxis per national guidelines as well as diagnosis and treatment of OIs, including fungal, bacterial and protozoal conditions. Pain is assessed during clinic visits and medication is provided to manage pain as indicated. \subsection{Outcomes} The primary outcome will be ART initiation within 90 days of enrollment among those eligible for ART (as definied below). 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. \\ Loss to follow-up (LTFU) is defined as those patients that are not deceased and have not had contact in 180 days prior to database close; this adheres to the Nigerian standard of more than 90 days late to scheduled visit. Because date of missed visit is not captured in the database, we select the missed visit date based on a maximum 3 month window between visits, so as not to overcount LTFU. \\ \noindent\textbf{\large Treatment eligibility in Nigeria}\\ \noindent Period 1: ART eligibility during June 2009 to May 2010 \squishlist \item WHO stage IV irrespective of CD4 cell count \item WHO stage III if CD4 cell count <350 \item WHO Stage I or II if CD4 cell count <200 \squishend Period 2: ART eligibility during June 2010 to present: \squishlist \item CD4 cell count $\le$ 350 regardless of WHO stage \item CD4 count >350 if WHO stage 3 or 4 \squishend \subsection{Multiple centers} There are five clinics in two states. Sobi and Lafiagi Hospitals are located in Kwara state. Gawu, Kuta, and Umaru YarAdua Hospitals are located in Niger State. \subsection{Data Sources and Measurements} The data extraction includes all patients enrolled on or before 31 May 2011. For analysis, we will include only those patients enrolled by 28 February 2011. Clinical characteristics collected closest to the date of enrollment up to a 90 day window before or after enrollment are deemed as enrollment status indicators (eg, weight, WHO stage, CD4) with the exception of height which was allowed a 365 day window in either direction. \\ VU/FGH uses CAREWare to record client level data at all HIV care and treatment sites. VU/FGH performs routine audits of medical records to ensure forms are completed accurately and laboratory data is entered correctly. VU/FGH is in the process of improving data management system for CAREWare to better ensure quality data. A data management plan has been developed and is in the final stage of being fully implemented in the coming months. \\ Demographic information includes: sex, age, marital status, educational status, occupational status, service entry into the program. Routine clinical data are as follows: weight, height, functional status, WHO stage, TB status. Labs that may be collected include: CD4, hepatitis B, VDRL, pregnancy, hepatitis C, hematocrit, hemoglobin, ALT, creatinine, total cholesterol, HDL, LDL, and triglycerides. \subsubsection{Data Cleaning} Data queries were generated for out of range and missing data. Each site addressed data queries and the clean data was extracted for final analyses. \\ Along with data query resolution, values of hemoglobin less than 1 (N=46) and greater than 18 (N=60) were imputed as missing. Similarly, values of CD4 count less than 0 (N=0) or greater than 1500 (N=7) were imputed as missing. Creatinine larger than 100 (N=10) and less than 0.1 (N=0) were imputed as missing. If height was recorded less than 100 (N=21) or larger than 220 (N=4), it is imputed as missing. If weight was recorded less than 20 (N=23) or larger than 120 (N=19), it is imputed as missing. BMI is calculated as weight in kilograms divided by the square of height in meters. Extreme BMI records below 10 kg/m$^2$ (N=11) and above kg/m$^2$ (N=38) were imputed as missing. \subsection{Statistical methods} {\it 1. To describe patients enrolled in HIV care and treatment.}\\ Summary characteristics will be tabulated over total enrollment and by sex.\\ \\ \noindent{\it 2. To describe patients who initiate ART among those eligible for ART.}\\ We will subset to those patients who were eligible for ART at enrollment into HIV care and treatment according to Nigerian national guidelines at the date of enrollment. Summary characteristics will be tabulated for patients eligible for ART by ART initiation within 90 days of enrollment 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, BMI, hemoglobin, and WHO stage at entry have been identified as predictors of interest. If needed, missing data methods will be considered. Multiple imputation is used to account for missing values of baseline predictors in order to prevent casewise deletion of missing data. \\ \\ \noindent{\it 3. 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, III, or IV participants by quarter of enrollment will summarize any trend. \\ \\ \noindent{\it 4. 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. \\ \\ \noindent{\it 5. To examine trends from enrollment or ART initiation of CD4 and hemoglobin.}\\ The trajectory of CD4 count and hemoglobin will be depicted using spaghetti plots overlaid with a lowess curve for two groups of patients: entire cohort (time since enrollment) and ART initiators (time since initiation). No formal comparisons are planned as this is a big picture summary and will not address issues of timely initiation and patient dropout.\\ \\ \noindent{\it 6. To estimate all-cause mortality and LTFU during the first year of ART.}\\ Kaplan-Meier estimates will be used to compute mortality and mortality/LTFU during the first 12 months on ART. We will estimate the cumulative incidence of LTFU treating death as a competing risk during the first 12 months. These estimates will include 95\% confidence intervals and they will be calculated only for those patients who initiated ART. \\ \\ \noindent{\it 7. To correct mortality for patients on ART.}\\ A PLoS Medicine article recently published$^1$ attempts to correct for mortality not observed in those lost to follow-up. The mortality observed among patients remaining in care can be multiplied by a correction factor $C$ to obtain an estimate of program-level mortality that takes death among patients lost to follow-up into account. This correction factor is obtained from a nomogram. This method applies only to one year estimates and for patients who have initiated ART.\\ \footnotetext[1]{Egger M, Spycher BD, Sidle J, Weigel R, Geng EH, et al. (2011) Correcting Mortality for Loss to Follow-Up: A Nomogram Applied to Antiretroviral Treatment Programmes in Sub-Saharan Africa. PLoS Med 8(1): e1000390.} \\ \noindent R-software 2.13.1 (www.r-project.org) will be used for data analyses. \listoftables\addcontentsline{toc}{section}{Index of Tables} \listoffigures\addcontentsline{toc}{section}{Index of Figures} \clearpage % Hidden R code chunk --- reading in the data <>= rm(list=ls()) # SHOULD BE FIRST LINE IN ALL R PROGRAMS - CLEARS NAMESPACE # Programs to tabulate counts/proportions and median(q1,q3) overall and by groups convertdate1 <- function(mydate){as.Date(mydate,"%Y-%m-%d")} convertdate <- function(mydate){as.Date(mydate,"%m/%d/%Y")} myformat1 <- function(x){format(round(x,1),nsmall=1)} 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) } ## FORCE NUMBER will take a factor/character to a number without issuing the warnings forcenumber <- function(var){ options(warn = -1) x <- as.numeric(gsub("[^0-9.]","",as.character(var))) options(warn = 1) return(x) } getnptex <- function(var,by,label,rowlabel,note=note,inc.pval=TRUE,var.order=FALSE,by.order=FALSE){ 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 p[i,which(d!=0)] <- sapply(100*n[i,which(d!=0)]/d[which(d!=0)],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))} } getmedtex <- function(var,by,label,note=note,inc.pval=TRUE,lower.label=TRUE,by.order=FALSE,digits=8){ if(!is.numeric(var)){stop("The variable is non-numeric")} s <- round(summary(var,digits=digits),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)) for(i in 1:length(level)){ sby[i,1:5] <- round(summary(var[by==level[i]],digits=digits)[1:5],1) textby[i] <- ifelse(!is.na(sby[i,3]),paste(sby[i,3]," (",sby[i,2],", ",sby[i,5],")",sep=""), "--") } if(inc.pval==TRUE){ if(length(level)==1) {pval <- " "} else if(length(level)>1 & by.order==FALSE) {pval <- myformat4(kruskal.test(var ~ by)$p.value)} else if(length(level)>1 & by.order==TRUE) {pval <- myformat4(cor.test(var, as.numeric(by),alternative="two.sided",method="spearman")$p.value)} if(missing(note)){line1 <- paste(label,paste(textby,collapse=' & '),text,pval,sep=' & ') } else{line1 <- paste(paste(label,"\\\\tnote{",note,"}",sep=""),paste(textby,collapse=' & '),text,pval,sep=' & ') } } else{ if(missing(note)){line1 <- paste(label,paste(textby,collapse=' & '),text,sep=' & ') } else{line1 <- paste(paste(label,"\\\\tnote{",note,"}",sep=""),paste(textby,collapse=' & '),text,sep=' & ') } } if(sum(is.na(var))>0){ no <- table(is.na(var),by)[2,] nop <- sapply(100*no/table(by),myper1) nopt <- sapply(100*sum(no)/sum(table(by)),myper1) if(lower.label==TRUE){line2 <- paste(paste("\\\\hspace{.1in}Missing ",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))} } getnp1line <- function(var,by,rowlabel,indent=TRUE,inc.pval=FALSE,by.order=FALSE){ 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=" & ")} if(inc.pval==TRUE){ warning("Function NP1LINE: assuming if var!=1 then var=0 for hypothesis test") 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)) } library(Hmisc) ## READ IN DATA FOR ANALYSIS source("read_data.R") ## TABLE 1 STUFF SEX.N <- table(analysis$gender) age.m1 <- getmedtex(var=analysis$age,by=analysis$gender,label="Age",note="a",inc.pval=TRUE) edu_codes <- levels(analysis$ph_pt_edu)[c(2,6,1,5,3,4)] edu.n1 <- getnptex(var=match(analysis$ph_pt_edu,edu_codes),by=analysis$gender,label="Education",inc.pval=TRUE,rowlabel=edu_codes) mar.n1 <- getnptex(var=analysis$ph_mrtl_stt,by=analysis$gender,label="Marital status",inc.pval=TRUE,rowlabel=levels(analysis$ph_mrtl_stt)) occ.n1 <- getnptex(var=analysis$ph_pt_occp,by=analysis$gender,label="Occupation",inc.pval=TRUE,rowlabel=levels(analysis$ph_pt_occp)) cli.n1 <- getnptex(var=analysis$clinic,by=analysis$gender,label="Clinic",inc.pval=TRUE,rowlabel=levels(analysis$clinic)) reg.n1 <- getnptex(var=analysis$med,by=analysis$gender,label="Baseline ART regimen",inc.pval=FALSE,rowlabel=levels(analysis$med)) ref.n1 <- getnptex(var=analysis$ph_entr_pnt,by=analysis$gender,label="Referral type",inc.pval=TRUE,rowlabel=levels(analysis$ph_entr_pnt)) preg.n1 <- getnp1line(var=ifelse(analysis$enpreg=="yes",1,NA),by=analysis$gender,rowlabel="Pregnant (enrollment)",indent=FALSE,inc.pval=FALSE) hei.m1 <- getmedtex(var=analysis$enht,by=analysis$gender,label="Height (cm)",note="ac",inc.pval=TRUE) wei.m1 <- getmedtex(var=analysis$enwt,by=analysis$gender,label="Weight (kg)",note="ac",inc.pval=TRUE) bmi.m1 <- getmedtex(var=analysis$enbmi ,by=analysis$gender,label="BMI (kg/m$^2$)",note="ac",inc.pval=TRUE,lower.label=FALSE) encd4.m1 <- getmedtex(var=analysis$encd4,by=analysis$gender,label="CD4 count",note="a",inc.pval=TRUE,lower.label=FALSE) encd4.n1 <- getnptex(var=analysis$cd4cat,by=analysis$gender,label="CD4 count category",inc.pval=FALSE,rowlabel=c("<50","51-200","201-350",">350")) enhem.m1 <- getmedtex(var=analysis$enhem,by=analysis$gender,label="Hemoglobin",note="a",inc.pval=TRUE) enhem.n1 <- getnptex(var=analysis$hemcat,by=analysis$gender,label="Hemoglobin category",inc.pval=FALSE,rowlabel=c("<8","8-10",">10")) encre.m1 <- getmedtex(var=analysis$encre,by=analysis$gender,label="Creatinine",note="a",inc.pval=TRUE) enwho.n1 <- getnptex(var=analysis$enwho,by=analysis$gender,label="WHO stage",inc.pval=TRUE,rowlabel=c("I","II","III","IV")) fs_codes <- levels(analysis$enfs)[c(2,1,3)] enfs.n1 <- getnptex(var=match(analysis$enfs,fs_codes),by=analysis$gender,label="Functional status",inc.pval=TRUE,rowlabel=fs_codes) art.n1 <- getnptex(var=analysis$art_status,by=analysis$gender,label="HAART status",inc.pval=TRUE,rowlabel=levels(analysis$art_status)) death.n1 <- getnp1line(var=ifelse(analysis$endeath12mo==1,1,NA),by=analysis$gender,rowlabel="Death in 12 months (enrollment)",indent=FALSE,inc.pval=TRUE) lost.n1 <- getnp1line(var=ifelse(analysis$enltfu12mo==1,1,NA),by=analysis$gender,rowlabel="Lost in 12 months (enrollment)",indent=FALSE,inc.pval=TRUE) table1 <- paste(cli.n1$line,age.m1$line,edu.n1$line.incmiss,mar.n1$line.incmiss,occ.n1$line.incmiss, ref.n1$line.incmiss,preg.n1,sep=" \\\\\\\\ ") table2 <- paste(hei.m1$line.incmiss,wei.m1$line.incmiss,bmi.m1$line.incmiss,enfs.n1$line.incmiss, encd4.m1$line.incmiss,encd4.n1$line.incmiss,enhem.m1$line.incmiss,enhem.n1$line.incmiss,encre.m1$line.incmiss,enwho.n1$line.incmiss, art.n1$line,death.n1,lost.n1,sep=" \\\\\\\\ ") txin90text <- paste("Includes ",sum(is.na(analysis$en2tx))," patients not initiating treatment prior to database cut date, and ", sum(analysis$en2tx>90,na.rm=T)," initiating >90 days from enrollment. Additionally, ",sum(analysis$art_status=="Not ART-naive"), " not ART-naive patients were dropped from this table and subsequent analysis.",sep="") eligtxin90text <- paste("Includes ",sum(is.na(arteligible$en2tx))," patients not initiating treatment prior to database cut date, and ", sum(arteligible$en2tx>90,na.rm=T)," initiating >90 days from enrollment.",sep="") ELIGTXIN90.N <- table(arteligible$txin90) age.m3 <- getmedtex(var=arteligible$age,by=arteligible$txin90,label="Age",note="a",inc.pval=TRUE) sex.n3 <- getnp1line(var=ifelse(arteligible$gender=="Female",1,NA),by=arteligible$txin90,rowlabel="Female",indent=FALSE,inc.pval=TRUE) edu.n3 <- getnptex(var=match(arteligible$ph_pt_edu,edu_codes),by=arteligible$txin90,label="Education",inc.pval=TRUE,rowlabel=edu_codes) mar.n3 <- getnptex(var=arteligible$ph_mrtl_stt,by=arteligible$txin90,label="Marital status",inc.pval=TRUE,rowlabel=levels(arteligible$ph_mrtl_stt)) occ.n3 <- getnptex(var=arteligible$ph_pt_occp,by=arteligible$txin90,label="Occupation",inc.pval=TRUE,rowlabel=levels(arteligible$ph_pt_occp)) cli.n3 <- getnptex(var=arteligible$clinic,by=arteligible$txin90,label="Clinic",inc.pval=TRUE,rowlabel=levels(arteligible$clinic)) reg.n3 <- getnptex(var=arteligible$med,by=arteligible$txin90,label="Baseline ART regimen",inc.pval=FALSE,rowlabel=levels(arteligible$med)) ref.n3 <- getnptex(var=arteligible$ph_entr_pnt,by=arteligible$txin90,label="Referral type",inc.pval=TRUE,rowlabel=levels(arteligible$ph_entr_pnt)) preg.n3 <- getnp1line(var=ifelse(arteligible$enpreg=="yes",1,NA),by=arteligible$txin90,rowlabel="Pregnant (enrollment)",indent=FALSE,inc.pval=TRUE) hei.m3 <- getmedtex(var=arteligible$enht,by=arteligible$txin90,label="Height (cm)",note="ac",inc.pval=TRUE) wei.m3 <- getmedtex(var=arteligible$enwt,by=arteligible$txin90,label="Weight (kg)",note="ac",inc.pval=TRUE) bmi.m3 <- getmedtex(var=arteligible$enbmi ,by=arteligible$txin90,label="BMI (kg/m$^2$)",note="ac",inc.pval=TRUE,lower.label=FALSE) encd4.m3 <- getmedtex(var=arteligible$encd4,by=arteligible$txin90,label="CD4 count",note="a",inc.pval=TRUE,lower.label=FALSE) enhem.m3 <- getmedtex(var=arteligible$enhem,by=arteligible$txin90,label="Hemoglobin",note="a",inc.pval=TRUE) encd4.n3 <- getnptex(var=arteligible$cd4cat,by=arteligible$txin90,label="CD4 count category",inc.pval=FALSE,rowlabel=c("<50","51-200","201-350",">350")) enhem.n3 <- getnptex(var=arteligible$hemcat,by=arteligible$txin90,label="Hemoglobin category",inc.pval=FALSE,rowlabel=c("<8","8-10",">10")) encre.m3 <- getmedtex(var=arteligible$encre,by=arteligible$txin90,label="Creatinine",note="a",inc.pval=TRUE) enwho.n3 <- getnptex(var=arteligible$enwho,by=arteligible$txin90,label="WHO stage",inc.pval=TRUE,rowlabel=c("I","II","III","IV")) enfs.n3 <- getnptex(var=match(arteligible$enfs,fs_codes),by=arteligible$txin90,label="Functional status",inc.pval=TRUE,rowlabel=fs_codes) death.n3 <- getnp1line(var=ifelse(arteligible$endeath12mo==1,1,NA),by=arteligible$txin90,rowlabel="Death in 12 months (enrollment)",indent=FALSE,inc.pval=TRUE) lost.n3 <- getnp1line(var=ifelse(arteligible$enltfu12mo==1,1,NA),by=arteligible$txin90,rowlabel="Lost in 12 months (enrollment)",indent=FALSE,inc.pval=TRUE) table9 <- paste(cli.n3$line,age.m3$line,sex.n3,edu.n3$line.incmiss,mar.n3$line.incmiss,occ.n3$line.incmiss, ref.n3$line.incmiss,preg.n3, sep=" \\\\\\\\ ") table10 <- paste(hei.m3$line.incmiss,wei.m3$line.incmiss,bmi.m3$line.incmiss,enfs.n3$line.incmiss, encd4.m3$line.incmiss,encd4.n3$line.incmiss,enhem.m3$line.incmiss,enhem.n3$line.incmiss,encre.m3$line.incmiss,enwho.n3$line.incmiss, # reg.n3$line.incmiss, death.n3,lost.n3,sep=" \\\\\\\\ ") ## FLOWCHART ## ADULTS >= 15 flow1 <- nrow(analysis) ## Not ART-naive flow3 <- sum(analysis$art_status=="Not ART-naive") ## Not able to assess eligiblity flow4 <- sum(analysis$noassess==1 & analysis$art_status!="Not ART-naive") ## Not eligible flow5 <- sum(is.na(analysis$arteligible) & analysis$noassess!=1 & analysis$art_status!="Not ART-naive") ## Eligible flow6 <- nrow(arteligible) ## Initiate ART flow7 <- sum(arteligible$txin90==1) ## FAIL to initiate ART flow8 <- sum(arteligible$txin90==0) @ \clearpage \section{Results} Data include all patients enrolled into care on or before 28 February 2011. \subsection{Cohort profile} \begin{center} \tikzstyle{every node}=[draw=black,thick,anchor=west] \tikzstyle{block} = [draw,rectangle,fill=blue!30,text width=11em,text centered, minimum height=10mm] \tikzstyle{block2} = [draw,rectangle,fill=green!30,text width=15em,text centered, minimum height=10mm] \tikzstyle{block3} = [draw,rectangle,fill=red!50,text width=10em,text centered, minimum height=10mm, node distance=5em] \tikzstyle{line} = [draw,-stealth,thick] \captionof{figure}{Flowchart of enrollment and ART eligibility} \begin{tikzpicture}[% grow via three points={one child at (0.5,-1.9) and two children at (0.5,-1.9) and (0.5,-3.55)}, edge from parent path={(\tikzparentnode.south) |- (\tikzchildnode.west)}] \node [block](p1) {\Sexpr{formatC(flow1,big.mark=",")} HIV-seropositive adults enrolled into HIV care by 28 February 2011} child { node [block2] {\Sexpr{formatC(flow3,big.mark=",")} (\Sexpr{myformat1(100*flow3/flow1)}\%) persons were not ART-na\"ive}} child { node [block2] {\Sexpr{formatC(flow4,big.mark=",")} (\Sexpr{myformat1(100*flow4/flow1)}\%) persons were not able to be evalutated for ART eligibility}}; \node [block,yshift=-15.0em](p2) {\Sexpr{formatC(flow5+flow6,big.mark=",")} HIV-seropositive adults assessed for ART eligibility within 90 days of enrollment} child { node [block2] {\Sexpr{formatC(flow5,big.mark=",")} (\Sexpr{myformat1(100*flow5/(flow5+flow6))}\%) Not eligible for ART}}; \node [block,yshift=-25.1em](p3) {\Sexpr{formatC(flow6,big.mark=",")} HIV-seropositive adults eligible for ART within 90 days of enrollment}; \node[block3,below of=p3,xshift=-7.5em](art1){\Sexpr{formatC(flow7,big.mark=",")} (\Sexpr{myformat1(100*flow7/flow6)}\%) Initiate ART}; \node[block3,below of=p3,xshift=7.5em](art2){\Sexpr{formatC(flow8,big.mark=",")} (\Sexpr{myformat1(100*flow8/flow6)}\%)\\ Fail to initiate ART}; \path[line](p1) -- (p2); \path[line](p2) -- (p3); \path[line](p3) -- (art1); \path[line](p3) -- (art2); \end{tikzpicture} \\[.5in] Because patients with WHO stage IV (and III after June 2010) are eligible for treatment regardless of CD4 count, there may be some bias in concluding that \Sexpr{myformat1(100*flow6/(flow5+flow6))}\% of patients were ART eligible. This could be an overestimate; it is owing to the fact that we cannot determine eligibility for patients missing CD4 count with WHO stage I and II (and III before June 2010). \end{center} %%%%%%%%%%%%%%%%%%%% % TABLE 1 % %%%%%%%%%%%%%%%%%%%% \begin{center} \setlength{\tabcolsep}{3pt} \begin{threeparttable} \footnotesize \centering{ \caption{Summary of Patient Demographics by Sex} \label{tab:base} \begin{tabular}{llllllll} \hline & Female & Male & Combined & P-value\tnote{b} \\ & (n=\Sexpr{SEX.N[1]}) & (n=\Sexpr{SEX.N[2]}) & (n=\Sexpr{sum(SEX.N)}) \\ \hline \Sexpr{table1}\\ %\\ \hline \end{tabular} } \begin{tablenotes} \item[a] Continuous variables are reported as medians (interquartile range). \item[b] To compare the distribution of study characteristics for participants by sex, we employ chi-square tests. Similarly, we use a Wilcoxon rank sum test for continuous variables by sex. \end{tablenotes} \end{threeparttable} \end{center} %%%%%%%%%%%%%%%%%%%% % TABLE 1 % %%%%%%%%%%%%%%%%%%%% \begin{center} \setlength{\tabcolsep}{3pt} \begin{threeparttable} \footnotesize \centering{ \caption{Summary of Patient Clinical Characteristics by Sex} \label{tab:base} \begin{tabular}{llllllll} \hline & Female & Male & Combined & P-value\tnote{b} \\ & (n=\Sexpr{SEX.N[1]}) & (n=\Sexpr{SEX.N[2]}) & (n=\Sexpr{sum(SEX.N)}) \\ \hline \Sexpr{table2}\\ %\\ \hline \end{tabular} } \begin{tablenotes} \item[a] Continuous variables are reported as medians (interquartile range). \item[b] To compare the distribution of study characteristics for participants by sex, we employ chi-square tests. Similarly, we use a Wilcoxon rank sum test for continuous variables by sex. \item[c] Weight, height, functional status, CD4, hemoglobin, creatinine, WHO stage are collected at enrollment. Enrollment data is collected in a window of $\pm$ 90 days from date of enrollment. \end{tablenotes} \end{threeparttable} \end{center} <>= ## LOGISTIC REGRESSION MODEL library(rms) set.seed(20120611) lrm1.d1 <- baseline[!is.na(baseline$arteligible) & baseline$arteligible==1,] lrm1.d1$notxin90 <- abs(lrm1.d1$txin90 - 1) lrm1.d1$ph_pt_occp[lrm1.d1$ph_pt_occp=="Retired"] <- NA lrm1.d1$ph_pt_occp <- droplevels(lrm1.d1$ph_pt_occp) lrm1.d1$enroll_dt <- as.numeric(lrm1.d1$ac_pc_enrl_dt) - 13000 lrm1.d1$enwho <- as.factor(lrm1.d1$enwho) lrm1.d1$sexandpreg <- as.character(lrm1.d1$gender) lrm1.d1$sexandpreg[lrm1.d1$enpreg=="yes"] <- "Pregnant Female" lrm1.d1$sexandpreg <- as.factor(lrm1.d1$sexandpreg) lrm1.dd1 <- datadist(with(lrm1.d1,data.frame(encd4, enhem, enwho, enfs, clinic, sexandpreg, ph_mrtl_stt, ph_pt_edu, ph_pt_occp, enroll_dt, age, notxin90, enbmi))) lrm1.mi <- aregImpute(~ encd4 + enhem + enwho + enfs + clinic + sexandpreg + ph_mrtl_stt + ph_pt_edu + ph_pt_occp + enroll_dt + age + notxin90 + enbmi,n.impute=25,data=lrm1.d1) options(datadist='lrm1.dd1') ########## MODEL 1: LOGISITIC REGRESSION lrm1.nl <- fit.mult.impute(notxin90 ~ clinic + rcs(age,4) + sexandpreg + ph_pt_edu + ph_mrtl_stt + ph_pt_occp + rcs(enbmi,4) + enfs + rcs(encd4,4) + rcs(enhem,4) + enwho + rcs(enroll_dt,5), lrm, lrm1.mi, lrm1.d1) lrm1 <- fit.mult.impute(notxin90 ~ clinic + age + sexandpreg + ph_pt_edu + ph_mrtl_stt + ph_pt_occp + enbmi + enfs + rcs(encd4,5) + enhem + enwho + rcs(enroll_dt,3), lrm, lrm1.mi, lrm1.d1) getlrm1 <- as.data.frame(anova(lrm1))$P cdates <- as.numeric(as.Date(c("06/30/2009","12/31/2009","06/30/2010","12/31/2010","05/31/2011"), "%m/%d/%Y")) - 13000 lrmsum1 <- summary(lrm1,clinic="GBRH",age=c(20,25),sexandpreg="Male",ph_pt_edu="None",ph_mrtl_stt="Married",ph_pt_occp="Employed",enbmi=c(20,21),enfs="Working", encd4=c(350,50),enhem=c(10,11),enwho=1,enroll_dt=c(cdates[1],cdates[2]))[c(12,14,16,18,2,20,22,32,24,30,26,28,34,36,38,40,42,44,46,4,48,50,6),c(4,6,7)] lrmsum2 <- summary(lrm1,clinic="GBRH",age=c(20,25),sexandpreg="Male",ph_pt_edu="None",ph_mrtl_stt="Married",ph_pt_occp="Employed",enbmi=c(20,21),enfs="Working", encd4=c(350,200),enhem=c(10,11),enwho=1,enroll_dt=c(cdates[1],cdates[2]))[c(6),c(4,6,7)] lrmsum3 <- summary(lrm1,clinic="GBRH",age=c(20,25),sexandpreg="Male",ph_pt_edu="None",ph_mrtl_stt="Married",ph_pt_occp="Employed",enbmi=c(20,21),enfs="Working", encd4=c(350,500),enhem=c(10,11),enwho=1,enroll_dt=c(cdates[1],cdates[2]))[c(6,8,52,54,56,10),c(4,6,7)] lrmsum4 <- summary(lrm1,clinic="GBRH",age=c(20,25),sexandpreg="Male",ph_pt_edu="None",ph_mrtl_stt="Married",ph_pt_occp="Employed",enbmi=c(20,21),enfs="Working", encd4=c(350,500),enhem=c(10,11),enwho=1,enroll_dt=c(cdates[1],cdates[3]))[c(10),c(4,6,7)] lrmsum5 <- summary(lrm1,clinic="GBRH",age=c(20,25),sexandpreg="Male",ph_pt_edu="None",ph_mrtl_stt="Married",ph_pt_occp="Employed",enbmi=c(20,21),enfs="Working", encd4=c(350,500),enhem=c(10,11),enwho=1,enroll_dt=c(cdates[1],cdates[4]))[c(10),c(4,6,7)] lrmsum6 <- summary(lrm1,clinic="GBRH",age=c(20,25),sexandpreg="Male",ph_pt_edu="None",ph_mrtl_stt="Married",ph_pt_occp="Employed",enbmi=c(20,21),enfs="Working", encd4=c(350,500),enhem=c(10,11),enwho=1,enroll_dt=c(cdates[1],cdates[5]))[c(10),c(4,6,7)] lrmsum <- rbind(lrmsum1,lrmsum2,lrmsum3,lrmsum4,lrmsum5) lrm1.labels <- c("Clinic","\\\\hspace{0.25in} GBRH (ref)","\\\\hspace{0.25in} Kuta","\\\\hspace{0.25in} LGH","\\\\hspace{0.25in} SBSH","\\\\hspace{0.25in} UMYMH", "Age (per 5 yrs)","Sex and pregnancy","\\\\hspace{0.25in} Male (ref)","\\\\hspace{0.25in} Female","\\\\hspace{0.25in} Pregnant Female", "Education","\\\\hspace{0.25in} None (ref)","\\\\hspace{0.25in} Started primary","\\\\hspace{0.25in} Completed primary","\\\\hspace{0.25in} Secondary","\\\\hspace{0.25in} Post secondary","\\\\hspace{0.25in} Qur'anic", "Marital status","\\\\hspace{0.25in} Married (ref)","\\\\hspace{0.25in} Divorced","\\\\hspace{0.25in} Separated","\\\\hspace{0.25in} Single","\\\\hspace{0.25in} Widowed", "Occupation","\\\\hspace{0.25in} Employed (ref)","\\\\hspace{0.25in} Other","\\\\hspace{0.25in} Student","\\\\hspace{0.25in} Unemployed","BMI (per 1 kg/m$^2$)", "Functional status","\\\\hspace{0.25in} Working (ref)","\\\\hspace{0.25in} Ambulatory","\\\\hspace{0.25in} Bedridden", "CD4 count","\\\\hspace{0.25in} 50","\\\\hspace{0.25in} 200","\\\\hspace{0.25in} 350 (ref)","\\\\hspace{0.25in} 500","Hemoglobin (per 1 g/dL)", "WHO stage","\\\\hspace{0.25in} I (ref)","\\\\hspace{0.25in} II","\\\\hspace{0.25in} III","\\\\hspace{0.25in} IV", "Month of enrollment","\\\\hspace{0.25in} June 2009 (ref)","\\\\hspace{0.25in} December 2009","\\\\hspace{0.25in} June 2010","\\\\hspace{0.25in} December 2010" ) myformat <- function(x){format(round(x,1),nsmall=1)} lrmsumlines1 <- paste(paste(sapply(lrmsum[,1],myformat2)," (",sapply(lrmsum[,2],myformat2),", ",sapply(lrmsum[,3],myformat2),")",sep=""),sep=" & ") lrmsumlines2 <- c("","1",lrmsumlines1[1:5],"","1",lrmsumlines1[6:7],"","1",lrmsumlines1[8:12],"","1",lrmsumlines1[13:16],"","1",lrmsumlines1[17:20],"","1",lrmsumlines1[21:22],"",lrmsumlines1[23:24],"1",lrmsumlines1[25:26],"","1",lrmsumlines1[27:29],"","1",lrmsumlines1[30:32]) lrm1pval <- c(myformat4(getlrm1[1]),rep("",5),sapply(getlrm1[c(2,3)],myformat4),rep("",3),myformat4(getlrm1[4]),rep("",6),myformat4(getlrm1[5]),rep("",5),myformat4(getlrm1[6]),rep("",4),sapply(getlrm1[c(7,8)],myformat4),rep("",3),myformat4(getlrm1[9]),rep("",4),sapply(getlrm1[c(11,12)],myformat4),rep("",4),myformat4(getlrm1[13]),rep("",4)) lrm1table <- paste(lrm1.labels,lrmsumlines2,lrm1pval,sep=" & ") lrm1tabletex <- paste(lrm1table,collapse=" \\\\\\\\ ") @ \clearpage \subsection{ART initiation among ART eligible} %%%%%%%%%%%%%%%%%%%% % TABLE 1 % %%%%%%%%%%%%%%%%%%%% \begin{center} \setlength{\tabcolsep}{3pt} \begin{threeparttable} \footnotesize \centering{ \caption{Summary of Patient Demographics by HAART Initiation among those eligible to start HAART} \label{tab:base} \begin{tabular}{llllllll} \hline & No HAART\tnote{d} & HAART in 90 days & Combined & P-value\tnote{b} \\ & (n=\Sexpr{ELIGTXIN90.N[1]}) & (n=\Sexpr{ELIGTXIN90.N[2]}) & (n=\Sexpr{sum(ELIGTXIN90.N)}) \\ \hline \Sexpr{table9}\\ \hline \end{tabular} } \begin{tablenotes} \item[a] Continuous variables are reported as medians (interquartile range). \item[b] To compare the distribution of study characteristics for participants by ART initiation in 90 days, we employ chi-square tests. Similarly, we use a Wilcoxon rank sum test for continuous variables by ART initiation in 90 days. \item[d] \Sexpr{eligtxin90text} \end{tablenotes} \end{threeparttable} \end{center} %%%%%%%%%%%%%%%%%%%% % TABLE 1 % %%%%%%%%%%%%%%%%%%%% \begin{center} \setlength{\tabcolsep}{3pt} \begin{threeparttable} \footnotesize \centering{ \caption{Summary of Patient Clinical Characteristics by HAART Initiation among those eligible to start HAART} \label{tab:base} \begin{tabular}{llllllll} \hline & No HAART\tnote{d} & HAART in 90 days & Combined & P-value\tnote{b} \\ & (n=\Sexpr{ELIGTXIN90.N[1]}) & (n=\Sexpr{ELIGTXIN90.N[2]}) & (n=\Sexpr{sum(ELIGTXIN90.N)}) \\ \hline \Sexpr{table10}\\ \hline \end{tabular} } \begin{tablenotes} \item[a] Continuous variables are reported as medians (interquartile range). \item[b] To compare the distribution of study characteristics for participants by ART initiation in 90 days, we employ chi-square tests. Similarly, we use a Wilcoxon rank sum test for continuous variables by ART initiation in 90 days. \item[c] Weight, height, functional status, CD4, hemoglobin, creatinine, WHO stage are collected at enrollment. Enrollment data is collected in a window of $\pm$ 90 days from date of enrollment. \item[d] \Sexpr{eligtxin90text} \end{tablenotes} \end{threeparttable} \end{center} \begin{center} \setlength{\tabcolsep}{3pt} \begin{threeparttable} \footnotesize \centering{ \caption{Logistic Regression: Failure to initiate HAART in 90 days among those eligible to start HAART} \begin{tabular}{llllllll} \hline & Odds Ratio (95\% CI) & P-value \\ \hline \Sexpr{lrm1tabletex}\\ \hline \end{tabular} } \begin{tablenotes} \item[a] There is evidence that the association between the log-odds of delayed HAART initiation is non-linear with CD4 count (p=\Sexpr{myformat4(getlrm1[10])}) and date of enrollment (p=\Sexpr{myformat4(getlrm1[14])}). \item[b] There are \Sexpr{sum(lrm1$freq)} patients included in this model. Missing values of baseline predictors are accounted for using multiple imputation. \end{tablenotes} \end{threeparttable} \end{center} \setkeys{Gin}{width=.45\textwidth} %default is 0.7 \begin{figure}[h!] \caption{Predicted Log Odds: Failure to initiate HAART in 90 days by CD4 count among those eligible to start HAART} \begin{center} \label{fig:ors} <>= p <- Predict(lrm1,encd4) print(plot(p,xlab="CD4 count at enrollment",adj.subtitle=FALSE)) @ \tiny\begin{verbatim}Adjusted to: \Sexpr{gsub("1845","2010-08-24",attr(p,"info")$adjust)}\end{verbatim} \end{center} \end{figure} \setkeys{Gin}{width=.45\textwidth} %default is 0.7 \begin{figure}[h!] \caption{Predicted Log Odds: Failure to initiate HAART in 90 days by date of enrollment among those eligible to start HAART} \begin{center} \label{fig:ors} <>= # as.Date(1842+13000, origin="1970-01-01") clabels <- c("Jun 2009","Dec 2009","Jun 2010","Dec 2010","May 2011") cdates <- as.numeric(as.Date(c("06/30/2009","12/31/2009","06/30/2010","12/31/2010","05/31/2011"), "%m/%d/%Y")) - 13000 p <- Predict(lrm1,enroll_dt) print(plot(p,scales=list(x=list(at=cdates, labels=clabels)),xlab="Date of enrollment",adj.subtitle=FALSE)) @ \tiny\begin{verbatim}Adjusted to: \Sexpr{gsub("1848","2010-08-27",attr(p,"info")$adjust)}\end{verbatim} \end{center} \end{figure} \clearpage \setkeys{Gin}{width=.55\textwidth} %default is 0.7 \begin{figure}[h!] \caption{Predicted Probability: Failure to initiate HAART in 90 days by date of enrollment among those eligible to start HAART} \begin{center} \label{fig:ors} <>= # as.Date(1848+13000, origin="1970-01-01") clabels <- c("Jun 2009","Dec 2009","Jun 2010","Dec 2010","May 2011") cdates <- as.numeric(as.Date(c("06/30/2009","12/31/2009","06/30/2010","12/31/2010","05/31/2011"), "%m/%d/%Y")) - 13000 p <- Predict(lrm1,enroll_dt,fun=plogis) print(plot(p,scales=list(x=list(at=cdates, labels=clabels)),xlab="Date of enrollment",ylab="Probability of HAART in 90 days",adj.subtitle=FALSE)) @ \tiny\begin{verbatim}Adjusted to: \Sexpr{gsub("1804","2010-07-14",attr(p,"info")$adjust)}\end{verbatim} \end{center} \end{figure} <>= clabels <- c("Jun 2009","August 2009","October 2009","Dec 2009","Feb 2010","Apr 2010","Jun 2010","August 2010","October 2010","Dec 2010","Feb 2011") cdates <- as.numeric(as.Date(c("06/30/2009","08/31/2009","10/31/2009","12/31/2009","02/28/2010","04/30/2010","06/30/2010","08/31/2010","10/31/2010","12/31/2010","02/28/2011"), "%m/%d/%Y")) - 13000 getp <- Predict(lrm1,enroll_dt=cdates,fun=plogis) getplines1 <- paste(clabels,paste(sapply(getp$yhat,myformat2)," (",sapply(getp$lower,myformat2),", ",sapply(getp$upper,myformat2),")",sep=""),sep=" & ") getptabletex <- paste(getplines1,collapse=" \\\\\\\\ ") @ \begin{center} \setlength{\tabcolsep}{3pt} \begin{threeparttable} \footnotesize \centering{ \caption{Predicted Probability: Failure to initiate HAART in 90 days among those eligible to start HAART} \begin{tabular}{llllllll} \hline & Probability (95\% CI) \\ \hline \Sexpr{getptabletex}\\ %\\ \hline \end{tabular} } %\begin{tablenotes} %\item[a] %\end{tablenotes} \end{threeparttable} \tiny\begin{verbatim}Adjusted to: \Sexpr{gsub("1848","2010-08-27",attr(getp,"info")$adjust)}\end{verbatim} \end{center} \begin{figure}[!h] \caption{ART initiation during the first 90 days after enrollment by sex} \begin{center} <>= arteligible$sexandpreg <- as.character(arteligible$gender) arteligible$sexandpreg[arteligible$enpreg=="yes"] <- "Pregnant Female" arteligible$sexandpreg <- as.factor(arteligible$sexandpreg) arteligible$en2tx90 <- arteligible$en2tx arteligible$en2tx90[is.na(arteligible$en2tx90)] <- 91 S<-Surv(arteligible$en2tx90,arteligible$txin90) Sx.1<-Surv(arteligible$en2tx90[arteligible$sexandpreg=="Female"],arteligible$txin90[arteligible$sexandpreg=="Female"]) Sx.2<-Surv(arteligible$en2tx90[arteligible$sexandpreg=="Pregnant Female"],arteligible$txin90[arteligible$sexandpreg=="Pregnant Female"]) Sx.3<-Surv(arteligible$en2tx90[arteligible$sexandpreg=="Male"],arteligible$txin90[arteligible$sexandpreg=="Male"]) Sx.pval <- as.data.frame(anova(cph(S ~ arteligible$sexandpreg, x=TRUE, y=TRUE)))$P[1] #png("graphics/InitiationBySexPreg.png",res=250,width=1000,height=1125, bg="transparent") par(mgp=c(2,1,0)) plot(survfit(Sx.2~1),xlab="Time (days)",ylab="Proportion on ART",axes=FALSE,conf.int=FALSE,mark.time=FALSE, xmax=90,fun="event",lty=3,lwd=1.5,col="#CC00CC",ylim=c(0,0.82)) box() axis(side=2) axis(side=1,at=c(0,15,30,45,60,75,90),labels=c(0,15,30,45,60,75,90)) lines(survfit(Sx.1~1),lty=1,lwd=1,fun="event",mark.time=FALSE,xmax=90,col="#FF0033") lines(survfit(Sx.3~1),lty=4,lwd=1.5,fun="event",mark.time=FALSE,xmax=90,col="#3333CC") legend("bottomright",c("Non-Pregnant Females", "Pregnant Females", "Males"), col=c("#FF0033","#CC00CC","#3333CC"),lty=c(1,3,4), lwd=c(1,1.5,1.5),bty="n",cex=.65) legend("bottomleft",c("Log-rank test:",paste("p",myformat4(Sx.pval),sep="")), ,bty="n",cex=.65) #dev.off() @ \end{center} \end{figure} \begin{figure}[!h] \caption{ART initiation during the first 90 days after enrollment by clinic} \begin{center} <>= library(RColorBrewer) colvec <- brewer.pal(7, "Set1") Sx.C1<-Surv(arteligible$en2tx90[arteligible$clinic=="GBRH"],arteligible$txin90[arteligible$clinic=="GBRH"]) Sx.C2<-Surv(arteligible$en2tx90[arteligible$clinic=="Kuta"],arteligible$txin90[arteligible$clinic=="Kuta"]) Sx.C3<-Surv(arteligible$en2tx90[arteligible$clinic=="LGH"],arteligible$txin90[arteligible$clinic=="LGH"]) Sx.C4<-Surv(arteligible$en2tx90[arteligible$clinic=="SBSH"],arteligible$txin90[arteligible$clinic=="SBSH"]) Sx.C5<-Surv(arteligible$en2tx90[arteligible$clinic=="UMYMH"],arteligible$txin90[arteligible$clinic=="UMYMH"]) Sx.Cpval <- as.data.frame(anova(cph(S ~ arteligible$clinic, x=TRUE, y=TRUE)))$P[1] #png("graphics/InitiationBySexPreg.png",res=250,width=1000,height=1125, bg="transparent") par(mgp=c(2,1,0)) plot(survfit(S~1),xlab="Time (days)",ylab="Proportion on ART",axes=FALSE,conf.int=FALSE,mark.time=FALSE, xmax=90,fun="event",lty=3,lwd=1.5,col=0,ylim=c(0,0.82)) box() axis(side=2) axis(side=1,at=c(0,15,30,45,60,75,90),labels=c(0,15,30,45,60,75,90)) lines(survfit(Sx.C1~1),lty=1,lwd=1.5,fun="event",mark.time=FALSE,xmax=90,col=colvec[1]) lines(survfit(Sx.C2~1),lty=2,lwd=1.5,fun="event",mark.time=FALSE,xmax=90,col=colvec[2]) lines(survfit(Sx.C3~1),lty=6,lwd=1.5,fun="event",mark.time=FALSE,xmax=90,col=colvec[3]) lines(survfit(Sx.C4~1),lty=4,lwd=1.5,fun="event",mark.time=FALSE,xmax=90,col=colvec[4]) lines(survfit(Sx.C5~1),lty=5,lwd=1.5,fun="event",mark.time=FALSE,xmax=90,col=colvec[5]) legend("bottomright",c("GBRH","Kuta","LGH","SBSH","UMYMH"), col=c(colvec[1:5]),lty=c(1,2,6,4,5), lwd=c(1.5),bty="n",cex=.65) legend("bottomleft",c("Log-rank test:",paste("p",myformat4(Sx.Cpval),sep="")), ,bty="n",cex=.65) #dev.off() @ \end{center} \end{figure} <>= ## TREND ANALYSIS library(latticeExtra) quarter <- rep(NA,length(analysis$ac_pc_enrl_dt)) qdates <- c("06/01/2009","10/01/2009","02/01/2010","06/01/2010","10/01/2010","02/01/2011","06/01/2011") qdatesnum <- as.Date(qdates,"%m/%d/%Y") qmmyy <- format(qdatesnum, "%m/%y") qlabel <- c("Q1","Q2","Q3","Q4","Q5","Q6") for(i in 1:length(analysis$ac_pc_enrl_dt)){ if(qdatesnum[1] <= analysis$ac_pc_enrl_dt[i] & qdatesnum[2] > analysis$ac_pc_enrl_dt[i]){quarter[i] <- 1} if(qdatesnum[2] <= analysis$ac_pc_enrl_dt[i] & qdatesnum[3] > analysis$ac_pc_enrl_dt[i]){quarter[i] <- 2} if(qdatesnum[3] <= analysis$ac_pc_enrl_dt[i] & qdatesnum[4] > analysis$ac_pc_enrl_dt[i]){quarter[i] <- 3} if(qdatesnum[4] <= analysis$ac_pc_enrl_dt[i] & qdatesnum[5] > analysis$ac_pc_enrl_dt[i]){quarter[i] <- 4} if(qdatesnum[5] <= analysis$ac_pc_enrl_dt[i] & qdatesnum[6] > analysis$ac_pc_enrl_dt[i]){quarter[i] <- 5} if(qdatesnum[6] <= analysis$ac_pc_enrl_dt[i] & qdatesnum[7] > analysis$ac_pc_enrl_dt[i]){quarter[i] <- 6} } propsex <- table(analysis$gender,quarter)[1,]/(table(analysis$gender,quarter)[1,] + table(analysis$gender,quarter)[2,]) # Proportion FEMALE propwho1 <- table(analysis$enwho,quarter)[1,]/table(quarter[!is.na(analysis$enwho)]) # Proportion WHO stage 1 propwho2 <- table(analysis$enwho,quarter)[2,]/table(quarter[!is.na(analysis$enwho)]) # Proportion WHO stage 2 propwho3 <- table(analysis$enwho,quarter)[3,]/table(quarter[!is.na(analysis$enwho)]) # Proportion WHO stage 3 propwho4 <- table(analysis$enwho,quarter)[4,]/table(quarter[!is.na(analysis$enwho)]) # Proportion WHO stage 4 trend.d1 <- analysis trend.d1$enroll_dt <- as.numeric(trend.d1$ac_pc_enrl_dt) - 13000 trend.d1$enwho <- as.factor(trend.d1$enwho) trend.d1$enwho1 <- ifelse(trend.d1$enwho==1,1,0) trend.d1$enwho2 <- ifelse(trend.d1$enwho==2,1,0) trend.d1$enwho3 <- ifelse(trend.d1$enwho==3,1,0) trend.d1$enwho4 <- ifelse(trend.d1$enwho==4,1,0) trend.d1$hascd4 <- !is.na(trend.d1$encd4) trend.d1$hashem <- !is.na(trend.d1$enhem) trend.d1$female <- ifelse(trend.d1$gender=="Female",1,0) trend.dd1 <- datadist(with(trend.d1,data.frame(enroll_dt,encd4, enhem,hascd4,hashem, enwho,enwho1,enwho2,enwho3,enwho4, enfs, clinic, gender,ph_mrtl_stt, ph_pt_edu, ph_pt_occp, age, txin90, enbmi))) options(datadist='trend.dd1') tr1 <- lrm(female~rcs(enroll_dt,5),data=trend.d1) tr2 <- lrm(enwho~rcs(enroll_dt,5),data=trend.d1) tr2a <- lrm(enwho1~rcs(enroll_dt,5),data=trend.d1) tr2b <- lrm(enwho2~rcs(enroll_dt,5),data=trend.d1) tr2c <- lrm(enwho3~rcs(enroll_dt,5),data=trend.d1) tr2d <- lrm(enwho4~rcs(enroll_dt,5),data=trend.d1) tr3 <- lrm(hascd4~rcs(enroll_dt,5),data=trend.d1) tr4 <- lrm(hashem~rcs(enroll_dt,5),data=trend.d1) clabels <- c("Jun 2009","Dec 2009","Jun 2010","Dec 2010","May 2011") cdates <- as.numeric(as.Date(c("06/30/2009","12/31/2009","06/30/2010","12/31/2010","05/31/2011"), "%m/%d/%Y")) - 13000 myadj1 <- 1 myadj2 <- -0.5 colvec <- brewer.pal(7, "Set1") civec <- brewer.pal(7, "Pastel1") prepanel.ci <- function(x, y, ly, uy, subscripts, ...) { x <- as.numeric(x) ly <- as.numeric(ly[subscripts]) uy <- as.numeric(uy[subscripts]) list(ylim = range(y, uy, ly, finite = TRUE)) } panel.ci <- function(x, y, ly, uy, subscripts, pch = 16, density=4, angle=45, col = 'black', ...) { x <- as.numeric(x) y <- as.numeric(y) ly <- as.numeric(ly[subscripts]) uy <- as.numeric(uy[subscripts]) panel.polygon(c(x,rev(x)), c(uy, rev(ly)), col=col, density=density, angle=angle, border=FALSE, ...) } p1 <- Predict(tr2a,enroll_dt,fun=plogis) p2 <- Predict(tr2b,enroll_dt,fun=plogis) p3 <- Predict(tr2c,enroll_dt,fun=plogis) p4 <- Predict(tr2d,enroll_dt,fun=plogis) p1$group <- 1 p2$group <- 2 p3$group <- 3 p4$group <- 4 p <- rbind(p1,p2,p3,p4) trendwho <- xyplot(yhat ~ enroll_dt, groups=group,lwd=2, ly=p$lower,uy=p$upper,data=p,lty=c(1,2,3,4),type="l",col.line=colvec,col=civec, aspect = 1, xlab="Date of Enrollment",ylab="Probability",scales=list(x=list(at=cdates, labels=clabels)),prepanel=prepanel.ci, panel=function(x,y,...){ panel.superpose(x,y,panel.groups=panel.ci, ...) panel.xyplot(x, y, ...)}, par.settings = list( axis.components = list( left = list(pad1 = myadj1, pad2 = myadj2) ,bottom = list(pad1 = myadj1, pad2 = myadj2) ,top = list(pad1 = myadj1, pad2 = myadj2) ,right = list(pad1 = myadj1, pad2 = myadj2) ))) trendwho <- trendwho + layer(panel.key(c("Stage I","Stage II","Stage III","Stage IV"), corner = c(.02,0.98), lines = TRUE, points = FALSE, packets = 1,cex=0.9)) # FUNCTION BELOW CHANGES THE COLORS IN THE KEY TO MATCH THE COLORS IN MY PLOT update(trendwho, par.settings = custom.theme(symbol = colvec, fill = civec, lwd=2, lty=c(1:4))) clabels <- c("Jun 2009","Dec 2009","Jun 2010","Dec 2010","May 2011") cdates <- as.numeric(as.Date(c("06/30/2009","12/31/2009","06/30/2010","12/31/2010","05/31/2011"), "%m/%d/%Y")) - 13000 p <- Predict(tr1,enroll_dt,fun=plogis) p$group <- 1 trendsex <- xyplot(yhat ~ enroll_dt, groups=group,lwd=2, ly=p$lower,uy=p$upper,data=p,lty=c(1),type="l",col.line=colvec,col=civec, aspect = 1, xlab="Date of Enrollment",ylab="Probability",scales=list(x=list(at=cdates, labels=clabels)),prepanel=prepanel.ci, panel=function(x,y,...){ panel.superpose(x,y,panel.groups=panel.ci, ...) panel.xyplot(x, y, ...)}, par.settings = list( axis.components = list( left = list(pad1 = myadj1, pad2 = myadj2) ,bottom = list(pad1 = myadj1, pad2 = myadj2) ,top = list(pad1 = myadj1, pad2 = myadj2) ,right = list(pad1 = myadj1, pad2 = myadj2) ))) trendsex <- trendsex + layer(panel.key(c("Females"), corner = c(.98,0.02), lines = TRUE, points = FALSE, packets = 1,cex=0.9)) # FUNCTION BELOW CHANGES THE COLORS IN THE KEY TO MATCH THE COLORS IN MY PLOT update(trendsex, par.settings = custom.theme(symbol = colvec, fill = civec, lwd=2, lty=c(1))) propcd4 <- table(trend.d1$hascd4,quarter)[2,]/(table(trend.d1$hascd4,quarter)[1,] + table(trend.d1$hascd4,quarter)[2,]) # Proportion FEMALE prophem <- table(trend.d1$hashem,quarter)[2,]/(table(trend.d1$hashem,quarter)[1,] + table(trend.d1$hashem,quarter)[2,]) # Proportion FEMALE p1 <- Predict(tr3,enroll_dt,fun=plogis) p2 <- Predict(tr4,enroll_dt,fun=plogis) p1$group <- 1 p2$group <- 2 p <- rbind(p1,p2) trenddatcol <- xyplot(yhat ~ enroll_dt, groups=group,lwd=2, ly=p$lower,uy=p$upper,data=p,lty=c(1,2),type="l",col.line=colvec,col=civec, aspect = 1, xlab="Date of Enrollment",ylab="Probability",scales=list(x=list(at=cdates, labels=clabels)),prepanel=prepanel.ci, panel=function(x,y,...){ panel.superpose(x,y,panel.groups=panel.ci, ...) panel.xyplot(x, y, ...)}, par.settings = list( axis.components = list( left = list(pad1 = myadj1, pad2 = myadj2) ,bottom = list(pad1 = myadj1, pad2 = myadj2) ,top = list(pad1 = myadj1, pad2 = myadj2) ,right = list(pad1 = myadj1, pad2 = myadj2) ))) trenddatcol <- trenddatcol + layer(panel.key(c("CD4 collection","Hemoglobin colleciton"), corner = c(.02,0.02), lines = TRUE, points = FALSE, packets = 1,cex=0.9)) # FUNCTION BELOW CHANGES THE COLORS IN THE KEY TO MATCH THE COLORS IN MY PLOT update(trenddatcol, par.settings = custom.theme(symbol = colvec, fill = civec, lwd=2, lty=c(1:2))) ## CREATE SUMMARY TABLE QUART.N <- table(quarter) cumen.t2 <- paste(c("Cumulative enrollment",cumsum(table(quarter))),collapse=" & ") sex.t2 <- getnp1line(var=ifelse(trend.d1$gender=="Female",1,NA),by=quarter,rowlabel="Female",indent=FALSE,inc.pval=FALSE) enwho.t2 <- getnptex(var=trend.d1$enwho,by=quarter,label="WHO stage",inc.pval=FALSE,rowlabel=c("I","II","III","IV")) encd4.t2 <- getmedtex(var=trend.d1$encd4,by=quarter,label="CD4 count",note="a",inc.pval=FALSE) hascd4.t2 <- getnp1line(var=ifelse(trend.d1$hascd4,1,NA),by=quarter,rowlabel="CD4 data collection",indent=FALSE,inc.pval=FALSE) hashem.t2 <- getnp1line(var=ifelse(trend.d1$hashem,1,NA),by=quarter,rowlabel="Hemoglobin data collection",indent=FALSE,inc.pval=FALSE) trendtable <- paste(cumen.t2,sex.t2,enwho.t2$line,encd4.t2$line.incmiss,hascd4.t2,hashem.t2,sep=" \\\\\\\\ ") @ \begin{center} \setlength{\tabcolsep}{3pt} \begin{threeparttable} \footnotesize \centering{ \caption{Quarterly Key Indicators} \label{tab:base} \begin{tabular}{p{3.5cm}p{1.75cm}p{1.75cm}p{1.75cm}p{1.75cm}p{1.75cm}p{1.75cm}p{1.75cm}} \hline & \Sexpr{paste(qlabel,collapse=" & ")} & Combined \\ & \Sexpr{paste(paste("(n=",QUART.N,")",sep=""),collapse=" & ")} & (n=\Sexpr{sum(QUART.N)}) \\ \hline \Sexpr{trendtable}\\ %\\ \hline \end{tabular} } \begin{tablenotes} \item[a] Continuous variables are reported as medians (interquartile range). \end{tablenotes} \end{threeparttable} \end{center} \clearpage \subsection{Enrollment trends} \setkeys{Gin}{width=.55\textwidth} %default is 0.7 \begin{figure}[h!] \caption{Proportion of Female Enrollees by Quarter} \vspace{-.25in} \begin{center} \label{fig:rawfem} <>= plot(1:6,seq(0,1,.2),axes=F,col=0,xlim=c(0.5,6.5),ylim=c(0.5,0.8),ylab="Proportion of Enrollees",xlab="Quarter of Enrollment") box() axis(side=2) axis(side=1,at=c(1:6),labels=qlabel) axis(side=3,at=c(1:7)-0.5,labels=qmmyy,tick=F) abline(h=0.5,lty=2) lines(1:6,propsex,lwd=1.5) points(1:6,propsex,cex=0.7) legend(1,.8,"Female",col=1,lty=1,lwd=1.5,bty="n",cex=.8) @ \end{center} \end{figure} \setkeys{Gin}{width=.55\textwidth} %default is 0.7 \begin{figure}[h!] \caption{Predicted Probability: Female enrollment by date of enrollment} \vspace{-.25in} \begin{center} \label{fig:ors} <>= #png("graphics/lrm_sex.png",res=250,width=1000,height=1125, bg="transparent") update(trendsex, par.settings = custom.theme(symbol = colvec, fill = civec, lwd=2, lty=c(1))) #dev.off() @ \\ \tiny Test of association between date of enrollment and odds of female enrollee: p=\Sexpr{myformat4(as.data.frame(anova(tr1))$P[1])}. \end{center} \end{figure} \clearpage \setkeys{Gin}{width=.55\textwidth} %default is 0.7 \begin{figure}[h!] \caption{Proportion of WHO Stages by Quarter} \vspace{-.25in} \begin{center} \label{fig:rawstage} <>= plot(1:6,seq(0,1,.2),axes=F,col=0,xlim=c(0.5,6.5),ylim=c(0,0.6),ylab="Proportion of Enrollees",xlab="Quarter of Enrollment") box() axis(side=2) axis(side=1,at=c(1:6),labels=qlabel) axis(side=3,at=c(1:7)-0.5,labels=qmmyy,tick=F) lines(1:6,propwho1,lty=1,col=gray(.4), lwd=1.5) points(1:6,propwho1,cex=0.7,col=gray(.4)) lines(1:6,propwho2,lty=2,col=gray(.4), lwd=1.5) points(1:6,propwho2,cex=0.7,col=gray(.4)) lines(1:6,propwho3,lty=3) points(1:6,propwho3,cex=0.7) lines(1:6,propwho4,lty=4) points(1:6,propwho4,cex=0.7) legend("topleft",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} \setkeys{Gin}{width=.55\textwidth} %default is 0.7 \begin{figure}[h!] \caption{Predicted Probability: WHO staging by date of enrollment} \vspace{-.25in} \begin{center} \label{fig:ors} <>= #png("graphics/lrm_who.png",res=250,width=1000,height=1125, bg="transparent") update(trendwho, par.settings = custom.theme(symbol = colvec, fill = civec, lwd=2, lty=c(1:4))) #dev.off() @ \\ \tiny Test of association between date of enrollment and proprtional odds of higher staged enrollee: p\Sexpr{myformat4(as.data.frame(anova(tr2))$P[1])}. \end{center} \end{figure} \clearpage \subsection{Lab collection and trajectory} \setkeys{Gin}{width=.55\textwidth} %default is 0.7 \begin{figure}[h!] \caption{Proportion of patients with enrollment CD4 or hemoglobin by Quarter} \vspace{-.25in} \begin{center} \label{fig:rawstage} <>= plot(1:6,seq(0,1,.2),axes=F,col=0,xlim=c(0.5,6.5),ylim=c(0.5,1),ylab="Proportion of Enrollees",xlab="Quarter of Enrollment") box() axis(side=2) axis(side=1,at=c(1:6),labels=qlabel) axis(side=3,at=c(1:7)-0.5,labels=qmmyy,tick=F) lines(1:6,propcd4,lty=1,col=gray(.4), lwd=1.5) points(1:6,propcd4,cex=0.7,col=gray(.4)) lines(1:6,prophem,lty=2) points(1:6,prophem,cex=0.7) legend("topleft",c("CD4 collection","Hemoglobin collection"), col=c(gray(.4),1),lty=c(1,2), lwd=c(1.5,1),bty="n",cex=.8) @ \end{center} \end{figure} \setkeys{Gin}{width=.55\textwidth} %default is 0.7 \begin{figure}[h!] \caption{Predicted Probability: CD4 and hemoglobin collection by date of enrollment} \vspace{-.25in} \begin{center} \label{fig:ors} <>= #png("graphics/lrm_collection.png",res=250,width=1000,height=1125, bg="transparent") update(trenddatcol, par.settings = custom.theme(symbol = colvec, fill = civec, lwd=2, lty=c(1:2))) #dev.off() @ \\ \tiny Test of association between date of enrollment and probability of CD4 and hemoglobin collection, respectively: p\Sexpr{myformat4(as.data.frame(anova(tr3))$P[1])} and p\Sexpr{myformat4(as.data.frame(anova(tr4))$P[1])}. \end{center} \end{figure} \begin{figure}[htp] \caption{CD4 and hemoglobin trajectories} \vspace{-.5in} \begin{center} \hspace{-0.05\textwidth} \subfigure[CD4, entire cohort]{\includegraphics[width=3.7in]{graphics/cd4long.png}} \hspace{-0.05\textwidth} \subfigure[CD4, ART initiators]{\includegraphics[width=3.7in]{graphics/cd4longtx.png}} \hspace{-0.1\textwidth}\\ \hspace{-0.05\textwidth} \subfigure[Hemoglobin, entire cohort]{\includegraphics[width=3.7in]{graphics/hemlong.png}} \hspace{-0.05\textwidth} \subfigure[Hemoglobin, ART initiators]{\includegraphics[width=3.7in]{graphics/hemlongtx.png}} \hspace{-0.1\textwidth} \end{center} \end{figure} <>= onart <- analysis[!is.na(analysis$haart_start_date),] library(cmprsk) getltfucmprsk <- function(censorindicator,subset=TRUE){ time2death <- as.numeric(onart$deathdt[subset] - onart$haart_start_date[subset]) time2event <- as.numeric(onart$lastdate[subset] - onart$haart_start_date[subset]) time <- as.numeric(onart$lastdate[subset] - onart$haart_start_date[subset]) time2event[time2event < 0] <- NA time2death[time2death < 0] <- NA time[time < 0] <- NA incevent <- 1*(censorindicator[subset] & is.na(onart$deathdt[subset]) & time2event <= time) + 2*(!is.na(onart$deathdt[subset]) & time2death<=time) I <- cuminc(time,incevent) Ionly <- list(x=I$"1 1"$time,y=I$"1 1"$est) return(list(incplot=Ionly,inc=I,event=incevent,time=time)) } testltfucmprsk <- function(censorindicator,group,subset=TRUE){ time2death <- as.numeric(onart$deathdt[subset] - onart$haart_start_date[subset]) time2event <- as.numeric(onart$lastdate[subset] - onart$haart_start_date[subset]) time <- as.numeric(onart$lastdate[subset] - onart$haart_start_date[subset]) time2event[time2event < 0] <- NA time2death[time2death < 0] <- NA time[time < 0] <- NA incevent <- 1*(censorindicator[subset] & is.na(onart$deathdt[subset]) & time2event <= time) + 2*(!is.na(onart$deathdt[subset]) & time2death<=time) I <- cuminc(time,incevent,group=group[subset]) Ionly <- list(x=I$"1 1"$time,y=I$"1 1"$est) return(list(incplot=Ionly,inc=I,event=incevent,time=time)) } getcumincest <- function(itime,inc,label,alpha=0.05){ ## Aalen's variance estimator -- Delta method ## CI is based on 1-F(t) = S ## (S)^exp(+-A) ## A=zalpha*sqrt(Var(S) / (S*logS) attime=timepoints(inc$inc,times=itime) est=1-attime$est[1,] se=sqrt(attime$var[1,]) zalpha=qnorm(1-alpha/2) A=zalpha*se/(est*log(est)) L1=est^(exp(A)) L2=est^(exp(-A)) inc.sum1 <- paste(label,paste(itime,"days"),paste(sapply(1-est,myper2)," & (",sapply(1-L1,myper2),", ",sapply(1-L2,myper2),")",sep=""),sep=" & ") return(inc.sum1) } onart$death12mo<-ifelse(onart$death=="Death" & onart$tx2last<365.25,1,0) onart$fup12mo <- ifelse(onart$tx2last<365.25,onart$tx2last,365.25) onart$ltfu12mo <- ifelse(onart$death12mo==0 & onart$tx2last<365.25 & convertdate("05/31/2011") - onart$lastdate > 180,1,0) S<-Surv(onart$fup12mo,onart$death12mo) L<-getltfucmprsk(onart$ltfu12mo==1) onart$deathlost12mo <- with(onart,apply(data.frame(death12mo,ltfu12mo),1,max)) SL<-Surv(onart$fup12mo,onart$deathlost12mo) ## SUMMARIZE CRUDE RATE OF ADHERENCE BY DURATION myper2 <- function(x){ y <- paste(format(round(100*x,1),nsmall=1),"\\\\%",sep="") return(y) } S.sum1 <- summary(survfit(S~1),times=c(90,180,365)) S.sum2 <- paste(c("Mortality","",""),paste(S.sum1$time,"days"),paste(sapply(1-S.sum1$surv,myper2)," & (",sapply(1-S.sum1$lower,myper2),", ",sapply(1-S.sum1$upper,myper2),")",sep=""),sep=" & ") S.sum3 <- paste(S.sum2,collapse=" \\\\\\\\ ") L.sum3 <- paste(paste(getcumincest(itime=c(90,180,365),inc=L,label=c("LTFU\\\\tnote{a}","",""))),collapse=" \\\\\\\\ ") SL.sum1 <- summary(survfit(SL~1),times=c(90,180,365)) SL.sum2 <- paste(c("Mortality/LTFU","",""),paste(SL.sum1$time,"days"),paste(sapply(1-SL.sum1$surv,myper2)," & (",sapply(1-SL.sum1$lower,myper2),", ",sapply(1-SL.sum1$upper,myper2),")",sep=""),sep=" & ") SL.sum3 <- paste(SL.sum2,collapse=" \\\\\\\\ ") @ \subsection{Mortality and loss to follow-up on patient on ART} \begin{center} \begin{threeparttable} \footnotesize \caption{Kaplan-Meier estimates of Mortality/LTFU} \centering{ \begin{tabular}{lllllll} \hline Event type & Treatment duration & Rate & (95\% CI) \\ \hline\hline %\\ \Sexpr{S.sum3}\\ \hline \Sexpr{SL.sum3}\\ \hline \Sexpr{L.sum3}\\ \hline\hline \end{tabular} } \begin{tablenotes} \item[a] Deaths are censored at date of death; thus, deceased patients are no longer included in the risk set for LTFU. \end{tablenotes} \end{threeparttable} \end{center} \subsection{Correcting for mortality in LTFU} A PLoS Medicine article recently published$^1$ attempts to correct for mortality not observed in those lost to follow-up. The mortality observed among patients remaining in care can be multiplied by a correction factor $C$ to obtain an estimate of program-level mortality that takes death among patients lost to follow-up into account. This correction factor is obtained from a nomogram. This method applies only to one year estimates and for patients who have initiated ART.\\ \footnotetext[1]{Egger M, Spycher BD, Sidle J, Weigel R, Geng EH, et al. (2011) Correcting Mortality for Loss to Follow-Up: A Nomogram Applied to Antiretroviral Treatment Programmes in Sub-Saharan Africa. PLoS Med 8(1): e1000390.} <>= ## GET INPUTS FROM THOSE WHO WERE NOT LOST S.EL <-Surv(onart$fup12mo[onart$ltfu12mo!=1],onart$death12mo[onart$ltfu12mo!=1]) mysum.S <- summary(survfit(S.EL~1),time=365) mysum.KM <- 1 - with(mysum.S,data.frame(surv,lower,upper)) M_nltfu <- round(mysum.KM[1]*100,1) # estimated mortality from a Kaplan Meier M_nltfu_lb <- round(mysum.KM[3]*100,1) # 95% CI lower bound M_nltfu_ub <- round(mysum.KM[2]*100,1) # 95% CI upper bound N_ltfu <- nrow(onart) # N potential LTFU n_ltfu <- sum(onart$ltfu12mo) # N true LTFU ## NOMOGRAM FUNCTION #'GNU S' R Code compiled by R2WASP v. 1.0.44 () #To cite this work: Martin WG Brinkhof, Ben D Spycher, 2010, Correcting Mortality Estimates for Loss to Follow-up in sub-Saharan African ART Programmes (v1.0) in Free Statistics Software (v1.1.23), Office for Research Development and Education, URL http://www.wessa.net/rwasp_mortality_mc.wasp/ #Source of accompanying publication: Egger M, Spycher BD, et al. under review (source updated ASAP). #Technical description: S_nltfu <- (100-M_nltfu)/100 # estimated survival at 1 year from a Kaplan Meier S_nltfu_lb <- (100-M_nltfu_ub)/100 # 95%CI lower bound S_nltfu_ub <- (100-M_nltfu_lb)/100 # 95%CI upper bound n_rep<-10^5 set.seed(12345) beta<-as.vector(c(-.04044087,.57287318)) V<-rbind(c(.00021312,-.0044198),c(-.0044198,.12357747)) tau2<-.3682446984157812 seu_S_nltfu<-as.numeric((log(-log(S_nltfu))-log(-log(S_nltfu_ub)))/1.96) sel_S_nltfu<-as.numeric((log(-log(S_nltfu_lb))-log(-log(S_nltfu)))/1.96) se_S_nltfu<-mean(c(seu_S_nltfu,sel_S_nltfu)) m_S_nltfu<-as.numeric(log(-log(S_nltfu))) dev_S_nltfu<-rnorm(n_rep, mean = m_S_nltfu, sd = se_S_nltfu) dev_S_nltfu<-exp(-exp(dev_S_nltfu)) dev_R_ltfu<-rbinom(n_rep, N_ltfu, n_ltfu/N_ltfu)/N_ltfu x_val<-cbind(100*dev_R_ltfu,1) # percent LTFU sd_M_ltfu<-sqrt(apply((x_val%*%V)*x_val,1,sum)+tau2) m_M_ltfu<-x_val%*%beta logitM_ltfu<-rnorm(rep(1,n_rep),mean=m_M_ltfu, sd=sd_M_ltfu) dev_M_ltfu<-1/(1+exp(-logitM_ltfu)) dev_M_corr<-(1-dev_R_ltfu)*(1-dev_S_nltfu)+dev_R_ltfu*dev_M_ltfu M_corr<-mean(dev_M_corr)*100 M_corr_ci<-quantile(dev_M_corr,c(0.025,0.975))*100 M_corr M_corr_ci @ Among all patients who initiated ART, the numbers at risk and lost to follow-up are \Sexpr{nrow(onart)} and \Sexpr{sum(onart$ltfu12mo)}, respectively. The Kaplan-Meier (95\% CI) estimate for mortality at 1 year of patients NOT lost to follow-up is \Sexpr{myformat1(mysum.KM[1]*100)}\% (\Sexpr{myformat1(mysum.KM[3]*100)}\% to \Sexpr{myformat1(mysum.KM[2]*100)}\%). The overall estimate (95\% CI) of program-level mortality at 1 year is \Sexpr{myformat1(M_corr)}\% (\Sexpr{myformat1(M_corr_ci[1])}\% to \Sexpr{myformat1(M_corr_ci[2])}\%).\\ \end{document}