%############################################################# %# %# Program: analysis-scripts-002.nw %# %# Project: Nutritional Causes for Early ART Mortality in %# Zambia: Relationships between Dietary Intake, %# Appetite, and Early Mortality among Low BMI %# Adults Initiating HIV Treatment in Zambia %# %# PI: John R. Koethe, MD %# Meridith Blevins, MS %# Claire Bosire, MSPH, RD %# Christopher Nyirenda, MD %# Edmond K. Kabagambe, DVM, PhD %# Albert Mwango, MD %# Webster Kasongo, MPH %# Isaac Zulu, MD, MPH %# Bryan E. Shepherd, PhD %# Douglas C. Heimburger, MD, MS* %# * corresponding author, douglas.heimburger at vanderbilt.edu %# %# Biostatistician/Programmer: Meridith Blevins, MS** %# ** contact programmer: meridith.blevins at vanderbilt.edu %# %# Statistician: Bryan E Shepherd, PhD %# %# Purpose: Reproducible research. %# %# Notes: This '.nw' file is a Sweave file which allows %# for results from R statistical software to be %# nested within a LaTeX document. %# %# Created: 28 February 2011 %# %# Revisions: %# %# %############################################################# % setwd("/home/blevinml/Projects/NEMART/Aim2"); Sweave("analysis-scripts-002.nw") \documentclass[11pt]{article} \usepackage[margin=0.75in]{geometry} \usepackage{setspace,relsize} % needed for latex(describe()), \code \usepackage{moreverb} % for verbatimtabinput \usepackage{url} \usepackage{amssymb} \usepackage{color} \usepackage[pdftex]{lscape} % allows tables to be landscape when \usepackage{longtable} % for subject listings \usepackage{graphicx} \usepackage{graphicx, subfigure} \usepackage{threeparttable} \usepackage{rotating} \newcommand{\HRule}{\rule{\linewidth}{0.5mm}} \begin{document} \SweaveOpts{prefix.string=graphics/plot} % Created a "graphics" subdirectory to save graph files in \begin{titlepage} \begin{center} % Upper part of the page \textsc{\Large Statistical Report:}\\[0.5cm] % Title \HRule \\[0.4cm] { \Large \bfseries NEMART:\\ Nutritional Causes for Early ART Mortality in Zambia\\[12pt] Dietary Analysis}\\[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 recent introduction of highly active anti-retroviral therapy (HAART) on a large scale among persons with HIV/AIDS in sub-Saharan Africa has been attended by a surprisingly high early mortality rate. At the University of Alabama at Birmingham's (UAB) Centre for Infectious Disease Research in Zambia (CIDRZ), 13.3\% of patients in the highest risk group died in the first 90 days of therapy, and most of these deaths occurred in the first 30 days. It is probable that the high prevalence of malnutrition in the region, worsened by the presence of HIV disease, is a factor in early ART mortality. It is possible that the refeeding syndrome (RS), which has been observed in multiple settings in which severely cachectic or marasmic patients are treated and fed, is responsible for a portion of the early deaths, perhaps even a substantial portion. Patients are at particular risk for RS when feeding is introduced too rapidly, when calories come predominantly from carbohydrate, and when serum phosphate levels drop significantly. Within several days, HAART can improve appetite and intestinal mucosal integrity, potentially resulting in increased food intake and/or nutrient absorption that exceeds the body's ability to meet the high phosphorus demands of carbohydrate metabolism. The resulting hypophosphatemia could lead to adverse outcomes, including cardiopulmonary decompensation and death. \subsection{Research Aims} \subsubsection{To summarize the change in dietary intake and appetite during the course of ART.} \subsubsection{To estimate the correlation between appetite scale and dietary intake.} \subsubsection{To estimate the association of 90 day mortality with dietary intake and appetite during the course of ART.} \subsubsection{To estimate the association of 90 day mortality/LTFU with dietary intake and appetite during the course of ART.} \section{Methods} \subsection{Participants} This is an observational cohort study of 142 adults who initiated ART between November 6, 2006 and November 12, 2007 in a Zambian clinic. Consecutive persons were enrolled on the day they initiated ART if they had a BMI < 16 kg/m$^2$ or CD4+ count < 50 cells/$\mu$L. ART was initiated according to Zambian Ministry of Health national guidelines, which parallel those of the World Health Organization. Participants or family members of participants who missed clinic visits were contacted to determine their vital status and encourage them to return to the clinic. Participants were followed for 12 weeks or until all follow-up was discontinued on December 27, 2007, whichever occurred first; after this, they continued in the routine ART program. \subsection{Outcomes} The primary outcomes are longitudinal measures of appetite and dietary intake (specifically energy, carbohydrate, protein, and fat). Predictors of interest are time, 90 day outcome, and BMI. \subsection{Data Sources and Measurements} At 0, 1, 2, 4, 8, and 12 weeks of ART, study staff conducted a symptom review and physical examination focused on symptoms and signs of malnutrition and refeeding syndrome. Dietary intake was estimated with 24-hour recalls at weeks 0, 1, 8, and 12. \subsection{Statistical Methods} {\it 1. To summarize the change in dietary intake and appetite.}\\ To describe the relationship of dietary intake and appetite over time, we will produce spaghetti plots. Similarly we will produce descriptive statistics by time point and test for change in dietary intake during the course of ART using a linear mixed effects model with random intercepts for patient. \\ \noindent {\it 2. To estimate the correlation between appetite scale and dietary intake.} \\ Graphical representations of appetite and dietary intake by visit will be produced. We are planning a mixed effects model for the $i^{th}$ patient at the $j^{th}$ time point, \begin{eqnarray*} Y_{ij} & = & \beta_1 + \beta_2t_{ij} + \beta_3\mathrm{APPETITE}_{ij} + \beta_4t_{ij}\times\mathrm{APPETITE}_{ij} + b_{1i} + e_{ij}. \end{eqnarray*} The null hypothesis of no appetite differences in the changes in dietary intake can be expressed as $\mathrm{H_0\hspace{-3pt}:}\hspace{3pt}\beta_4=0$. All interaction terms that do not significantly impact the model ($p>0.1$) are dropped from each model.\\ \noindent {\it 3. To estimate the association of 90 day mortality with dietary intake and appetite.}\\ To illustrate the relationship of dietary intake and appetite over time, we will produce spaghetti plots by 90 day outcome. Cox regression with baseline and time-varying covariates will be used to assess the relationship between time to death, baseline status and dietary intake or appetite. Baseline variables include age, sex, hemoglobin, baseline serum phosphate, CD4 count and BMI. We will use principal components analysis with age, phosphate, hemoglobin CD4 count and BMI as continuous variables to extract 1 principal factor to represent baseline status. The resulting principal factor(s) will be used to adjust the models. We use single imputation to fill in missing data for <6\% of predictors prior to extracting the principal components. Energy, protein, fat, carbohydrate intake, and appetite have been identified as post-ART time-dependent predictors of interest. Participants with no assessment will be dropped from analysis. Only those intake or appetite corresponding to the time between visits will be included for analysis. \\ \noindent {\it 4. To estimate the association of 90 day mortality/LTFU with dietary intake and appetite.}\\ We repeat the above Cox regression, the outcome is redefined as dead or lost. \\ R-software 2.11.1 (www.r-project.org) will be used for data analyses. \\ % Hidden R code chunk --- reading in the data <>= rm(list=ls()) # SHOULD BE FIRST LINE IN ALL R PROGRAMS - CLEARS NAMESPACE ################################################# # READ IN FULL DATASET "BASELINE" and "LONGFUP" # ################################################# source("read_data.R") ################################### # CUSTOMIZE DATA FORMAT FUNCTIONS # ################################### myformat2 <- function(x){format(round(x,2),nsmall=2)} myformat3 <- function(x){ y <- NULL if(x < 0.01){y <- "<0.01"} else{y <- format(round(x,2),nsmall=2)} return(y) } myformat4 <- function(x){ y <- NULL if(x < 0.001){y <- "<0.001"} else{y <- format(round(x,3),nsmall=3)} return(y) } myper1 <- function(x){ y <- NULL if(x < 0.05 & x!=0){y <- "<0.1"} else{y <- format(round(x,1),nsmall=1)} return(y) } format2 <- function(x){ y <- format(round(x,2),nsmall=2,scientific=FALSE) return(y) } ############################################# # CUSTOMIZE FUNCTIONS TO WRITE LATEX TABLES # ############################################# # WRITE FUNCTIONS WHICH tabulate counts/proportions and median(q1,q3) overall and by covariate getmedtex <- function(var,by,label,note=note,inc.pval=TRUE,lower.label=TRUE){ s <- round(summary(var),1) text <- paste(s[3]," (",s[2],", ",s[5],")",sep="") level <- sort(unique(by)) sby <- matrix(nrow=length(level), ncol=length(s)); textby <- rep(NA,length(level)) if(length(level)==1) {pval <- " "} else if(length(level)>1) {pval <- myformat3(kruskal.test(var ~ by)$p.value)} for(i in 1:length(level)){ sby[i,1:5] <- round(summary(var[by==level[i]])[1:5],1) textby[i] <- paste(sby[i,3]," (",sby[i,2],", ",sby[i,5],")",sep="") } if(inc.pval==TRUE){ if(missing(note)){line1 <- paste(label,paste(textby,collapse=' & '),text,pval,sep=' & ') } else{line1 <- paste(paste(label,"\\\\tnote{",note,"}",sep=""),paste(textby,collapse=' & '),text,pval,sep=' & ') } } else{ if(missing(note)){line1 <- paste(label,paste(textby,collapse=' & '),text,sep=' & ') } else{line1 <- paste(paste(label,"\\\\tnote{",note,"}",sep=""),paste(textby,collapse=' & '),text,sep=' & ') } } if(sum(is.na(var))>0){ no <- table(is.na(var),by)[2,] nop <- sapply(100*no/table(by),myper1) nopt <- sapply(100*sum(no)/sum(table(by)),myper1) if(lower.label==TRUE){line2 <- paste(paste("\\\\hspace{.1in}Missing ",tolower(label),", n(\\\\%)",sep=""),paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),paste(sum(no)," (",nopt,"\\\\%)",sep=""),sep=" & ")} else{line2 <- paste(paste("\\\\hspace{.1in}Missing ",label,", n(\\\\%)",sep=""),paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),paste(sum(no)," (",nopt,"\\\\%)",sep=""),sep=" & ")} line3 <- paste(line1,line2,sep=" \\\\\\\\ ") return(list(line=line1,line.incmiss=line3)) } else{return(list(line=line1))} } getmeantex <- function(var,by,label,note=note,inc.pval=TRUE,lower.label=TRUE){ m <- round(mean(var,na.rm=TRUE)) s <- round(sd(var,na.rm=TRUE)) text <- paste(m," (",s,")",sep="") level <- sort(unique(by)) sby <- matrix(nrow=length(level), ncol=2); textby <- rep(NA,length(level)) for(i in 1:length(level)){ sby[i,1] <- round(mean(var[by==level[i]],na.rm=TRUE)) sby[i,2] <- round(sd(var[by==level[i]],na.rm=TRUE)) textby[i] <- paste(sby[i,1]," (",sby[i,2],")",sep="") } if(inc.pval==TRUE){ if(length(level)==1) {pval <- " "} else if(length(level)>1) {pval <- myformat3(kruskal.test(var ~ by)$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))} } getnptex <- function(var,by,label,rowlabel,note=note,inc.pval=TRUE){ n <- table(var,by) d <- colSums(table(var,by)[1:nrow(n),]) p <- matrix(nrow=nrow(n), ncol=ncol(n)) pt <- rep(NA,nrow(n)) lines <- rep(NA,nrow(n)) for(i in 1:nrow(n)){ p[i,] <- sapply(100*n[i,]/d,myper1) pt[i] <- sapply(100*table(var)[i]/sum(table(var)[1:nrow(n)]),myper1) lines[i] <- paste(paste("\\\\hspace{.2in}",rowlabel[i]),paste(paste(n[i,]," (",p[i,],"\\\\%)",sep=""),collapse=" & "),paste(sum(n[i,])," (",pt[i],"\\\\%)",sep=""),sep=" & ") } # ADD PVALUE TO COLUMN HEADING IF PLANNED if(inc.pval==TRUE){ pval <- myformat3(chisq.test(n)$p.value) spaces <- rep(" ",1+length(table(by))) line1 <- paste(paste(label,", n(\\\\%)",sep=""),paste(spaces,collapse=" & "),pval,sep=" & ") } else{ line1 <- paste(label,", n(\\\\%)",sep="") } all.lines <- paste(line1,paste(lines,collapse=" \\\\\\\\ "),sep=" \\\\\\\\ ") if(sum(is.na(var))>0){ no <- table(is.na(var),by)[2,] nop <- sapply(100*no/table(by),myper1) nopt <- sapply(100*sum(no)/sum(table(by)),myper1) if(missing(note)){ line2 <- paste("\\\\hspace{.1in}Missing",paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),paste(sum(no)," (",nopt,"\\\\%)",sep=""),sep=" & ") } else{line2 <- paste(paste("\\\\hspace{.1in}Missing\\\\tnote{",note,"}",sep=""),paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),paste(sum(no)," (",nopt,"\\\\%)",sep=""),sep=" & ") } all.lines.miss <- paste(line1,line2,paste(lines,collapse=" \\\\\\\\ "),sep=" \\\\\\\\ ") return(list(line=all.lines,line.incmiss=all.lines.miss)) } else{return(list(line=all.lines))} } longmedtex <- function(var,by,caseid,label,note=note,inc.pval=TRUE,lower.label=TRUE){ s <- round(summary(var),1) level <- sort(unique(by)) sby <- matrix(nrow=length(level), ncol=length(s)); textby <- rep(NA,length(level)) if(length(level)==1) {pval <- " "} for(i in 1:length(level)){ sby[i,1:5] <- round(summary(var[by==level[i]])[1:5],1) if(!is.na(sby[i,3])){textby[i] <- paste(sby[i,3]," (",sby[i,2],", ",sby[i,5],")",sep="")} else if(is.na(sby[i,3])){textby[i] <- "-"} } if(inc.pval==TRUE & length(level)>1){ newdata1 <- data.frame(var,by,caseid)[!is.na(var),] newdata <- groupedData(var ~ by | caseid, data=newdata1) mylme <- lme(var ~ by, data=newdata, random=~1|caseid, na.action=na.omit,method="ML") pval <- myformat4(as.data.frame(anova(mylme))$'p-value'[2]) if(missing(note)){line1 <- paste(label,paste(textby,collapse=' & '),pval,sep=' & ') } else{line1 <- paste(paste(label,"\\\\tnote{",note,"}",sep=""),paste(textby,collapse=' & '),pval,sep=' & ') } } else{ if(missing(note)){line1 <- paste(label,paste(textby,collapse=' & '),sep=' & ') } else{line1 <- paste(paste(label,"\\\\tnote{",note,"}",sep=""),paste(textby,collapse=' & '),sep=' & ') } } if(sum(is.na(var))>0){ no <- table(is.na(var),by)[2,] nop <- sapply(100*no/table(by),myper1) if(lower.label==TRUE){line2 <- paste(paste("\\\\hspace{.1in}Missing ",tolower(label),", n(\\\\%)",sep=""),paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),sep=" & ")} else{line2 <- paste(paste("\\\\hspace{.1in}Missing ",label,", n(\\\\%)",sep=""),paste(paste(no," (",nop,"\\\\%)",sep=""),collapse=" & "),sep=" & ")} line3 <- paste(line1,line2,sep=" \\\\\\\\ ") return(list(line=line1,line.incmiss=line3)) } else{return(list(line=line1))} } ################################## # TABULATE/PROPORTIONS/QUARTILES # ################################## baseline$CarbPer <- 100 * with(baseline,4 * carbohydrate/energy) baseline$ProPer <- 100 * with(baseline,4 * TotalProtein/energy) baseline$FatPer <- 100 * with(baseline,9 * TotalFat/energy) ## BASELINE ONLY outcome90.N <- table(baseline$outcome90) sex.n1 <- getnptex(var=baseline$sex,by=baseline$outcome90,label="Sex",inc.pval=TRUE,rowlabel=c("Female","Male")) age.m1 <- getmedtex(var=baseline$age,by=baseline$outcome90,label="Age",note="a",inc.pval=TRUE) wgt.m1 <- getmedtex(var=baseline$weight,by=baseline$outcome90,label="Weight",note="a",inc.pval=TRUE) bmi.m1 <- getmedtex(var=baseline$bmi,by=baseline$outcome90,label="BMI",inc.pval=TRUE) cd4.m1 <- getmedtex(var=baseline$cd4,by=baseline$outcome90,label="CD4",inc.pval=TRUE) hem.m1 <- getmedtex(var=baseline$h30,by=baseline$outcome90,label="Hemoglobin",inc.pval=TRUE) app.m1 <- getmedtex(var=baseline$appetite,by=baseline$outcome90,label="Appetite",inc.pval=TRUE) app.n1 <- getnptex(var=baseline$appetite,by=baseline$outcome90,label="Appetite",inc.pval=TRUE,rowlabel=c("None","Little","Normal","Hungry")) nut1.m1 <- getmedtex(var=baseline$energy,by=baseline$outcome90,label="Total Energy (kcal)",inc.pval=TRUE) nut2.m1 <- getmedtex(var=baseline$carbohydrate,by=baseline$outcome90,label="Total Carbohydrate (g)",inc.pval=TRUE) nut2.p1 <- getmeantex(var=baseline$CarbPer,by=baseline$outcome90,label="Carbohydrate (\\\\% kcal)",inc.pval=FALSE,lower.label=FALSE) nut3.m1 <- getmedtex(var=baseline$TotalProtein,by=baseline$outcome90,label="Total Protein (g)",inc.pval=TRUE) nut3.p1 <- getmeantex(var=baseline$ProPer,by=baseline$outcome90,label="Protein (\\\\% kcal)",inc.pval=FALSE,lower.label=FALSE) nut4.m1 <- getmedtex(var=baseline$TotalFat,by=baseline$outcome90,label="Total Fat (g)",inc.pval=TRUE) nut4.p1 <- getmeantex(var=baseline$FatPer,by=baseline$outcome90,label="Fat (\\\\% kcal)",inc.pval=FALSE,lower.label=FALSE) reg.n1 <- getnptex(var=baseline$ARVregimen,by=baseline$outcome90,label="Regimen",inc.pval=FALSE,rowlabel=levels(baseline$ARVregimen)) baseline.table1 <- paste(wgt.m1$line.incmiss,age.m1$line.incmiss,wgt.m1$line.incmiss,bmi.m1$line.incmiss,cd4.m1$line.incmiss,hem.m1$line.incmiss,app.m1$line,app.n1$line,nut1.m1$line.incmiss, nut2.m1$line.incmiss,nut2.p1$line,nut3.m1$line.incmiss,nut3.p1$line,nut4.m1$line.incmiss,nut4.p1$line,reg.n1$line,sep=" \\\\\\\\ ") ## SIX STUDY VISITS attach(longfup) library(nlme) visit.N <- table(visit) set.seed(12345) wgt.m1 <- longmedtex(var=weight,by=visit,caseid=pepfarid,label="Weight",note="a",inc.pval=TRUE) bmi.m1 <- longmedtex(var=bmi_visit,by=visit,caseid=pepfarid,label="BMI",inc.pval=TRUE) app.m1 <- longmedtex(var=appetite,by=visit,caseid=pepfarid,label="Appetite",inc.pval=TRUE) nut1.m1 <- longmedtex(var=energy,by=visit,caseid=pepfarid,label="Total Energy (kcal)",inc.pval=TRUE) nut2.m1 <- longmedtex(var=carbohydrate,by=visit,caseid=pepfarid,label="Total Carbohydrate (g)",inc.pval=TRUE) nut3.m1 <- longmedtex(var=TotalProtein,by=visit,caseid=pepfarid,label="Total Protein (g)",inc.pval=TRUE) nut4.m1 <- longmedtex(var=TotalFat,by=visit,caseid=pepfarid,label="Total Fat (g)",inc.pval=TRUE) table1 <- paste(wgt.m1$line.incmiss,bmi.m1$line.incmiss,app.m1$line.incmiss,nut1.m1$line.incmiss, nut2.m1$line.incmiss,nut3.m1$line.incmiss,nut4.m1$line.incmiss,sep=" \\\\\\\\ ") library(ICC) eneicc <- ICCest(pepfarid,energy,longfup)$ICC caricc <- ICCest(pepfarid,carbohydrate,longfup)$ICC faticc <- ICCest(pepfarid,TotalFat,longfup)$ICC proicc <- ICCest(pepfarid,TotalProtein,longfup)$ICC appicc <- ICCest(pepfarid,appetite,longfup)$ICC @ \begin{landscape} \section{Results} \subsection{Baseline Patient Characteristics} %%%%%%%%%%%%%%%%%%%% % TABLE 1 % %%%%%%%%%%%%%%%%%%%% \begin{center} \setlength{\tabcolsep}{3pt} \begin{threeparttable} \footnotesize \centering{ \caption{Summary of Patient Visit Data by 90-day Outcome} \label{tab:nut1} \begin{tabular}{lllllllll} \hline & Alive & Dead & Lost & Combined & P-value$^b$ \\ & (n=\Sexpr{outcome90.N[1]}) & (n=\Sexpr{outcome90.N[2]}) & (n=\Sexpr{outcome90.N[3]}) & (n=\Sexpr{sum(outcome90.N)}) \\ \hline \Sexpr{baseline.table1}\\ %\\ \hline \end{tabular} } \begin{tablenotes} \item[a] Continuous variables are reported as medians (interquartile range). \item[b] To determine linear trend in the distribution of study characteristics over time, we test whether the coefficient of time (days from ART initiation) in a linear mixed effects model is different from zero. This approach ignores missing data. \end{tablenotes} \end{threeparttable} \end{center} \subsection{Post-ART Patient Characteristics} Table \ref{tab:nut1} summarizes patient weight and dietary recall by visit for all patients. Spaghetti plots are color-coded to show results for the alive, deceased and lost cohorts. A patient must have a minimum of two records on any of the 6 scheduled visits in order to appear in these plots. Each line represents one study participant. %%%%%%%%%%%%%%%%%%%% % TABLE 2 % %%%%%%%%%%%%%%%%%%%% \begin{center} \setlength{\tabcolsep}{3pt} \begin{threeparttable} \footnotesize \centering{ \caption{Summary of Patient Visit Data} \label{tab:nut1} \begin{tabular}{lllllllll} \hline & Visit 1 & Visit 2 & Visit 3 & Visit 4 & Visit 5 & Visit 6 & P-value$^b$ \\ & (n=\Sexpr{visit.N[1]}) & (n=\Sexpr{visit.N[2]}) & (n=\Sexpr{visit.N[3]}) & (n=\Sexpr{visit.N[4]}) & (n=\Sexpr{visit.N[5]}) & (n=\Sexpr{visit.N[6]}) \\ \hline \Sexpr{table1}\\ %\\ \hline \end{tabular} } \begin{tablenotes} \item[a] Continuous variables are reported as medians (interquartile range). \item[b] To determine linear trend in the distribution of study characteristics over time, we test whether the coefficient of time (days from ART initiation) in a linear mixed effects model is different from zero. This approach ignores missing data. \item[c] The intraclass correlation coefficient describes how strongly nutritional intake measures from the same individual resemble each other. The ICC for energy is \Sexpr{myformat2(eneicc)}, carbohydrate is \Sexpr{myformat2(caricc)}, protein is \Sexpr{myformat2(proicc)}, fat is \Sexpr{myformat2(faticc)}, and appetite is \Sexpr{myformat2(appicc)}. \end{tablenotes} \end{threeparttable} \end{center} \end{landscape} % Hidden R code chunk --- reading in the data <>= ## EXAMPLE CODE FOR SPAGHETTI PLOTS USING LATTICE ## library(lattice) my.panel <- function(x,y,...){ panel.xyplot(x,y,type=c('g','l'),col="#FFCC00") } panel.smoother <-function(x,y,...){ panel.superpose(x,y,panel.groups=my.panel,...) panel.loess(x,y,col="#000099",lwd=6,span=0.8) } time4 <- 30*(0:3) zeroout <- function(var){ newvar <- rep(NA,nrow(longfup)) for(i in unique(longfup$pepfarid)){ if(sum(!is.na(var[longfup$pepfarid==i]))>1){ newvar[which(longfup$pepfarid==i & !is.na(var))] <- 1 } } return(newvar) } longfup$enout <- zeroout(longfup$energy) ## PLOT PREDICTED ENERGY INTAKE OVER TIME (OUTPUT AS LMM_APP_ENERGY.PDF FOR LATEX INCLUSION) #### Energy (kcal) PANEL energymax <- max(energy,na.rm=TRUE) energymin <- min(energy,na.rm=TRUE) xyplot(energy ~ time, data = longfup, groups = pepfarid, ylim=c(energymin,energymax),xlim=c(0,90), scales=list(x=list(at=time4, labels=time4)), aspect = 1, subset=(outcome90==0 & enout==1),panel=panel.smoother, ylab="",xlab="") ## REPEAT THIS FOR DIFFERENT OUTCOMES AND INTAKE/APPETITE @ % Hidden R code chunk --- reading in the data <>= ## HOW IS APPETITE ASSOCIATED WITH DIETARY INTAKE? ## LINEAR MIXED MODELS WITH DIETARY AS OUTCOME myapp_lmm <- function(var,label){ longfup$var <- var lmm_app <- lme(var ~ time*as.factor(appetite), data=longfup, random=~1|pepfarid, na.action=na.omit,method="ML") sum_app <- summary(lmm_app) ### TOOK OUT REFERENCE GROUPS FOR EASE OF INTERPRETATION # DERIVATION OF ESTIMATES WITHOUT REFERENCE GETS COMPLICATED WITH COVARIANCE sum1.est <- sapply(sum_app$tTable[1:8,1],format2) sum1.est <- sapply(c(sum_app$tTable[1,1], sum_app$tTable[3:5,1]+sum_app$tTable[1,1], sum_app$tTable[2,1], sum_app$tTable[6:8,1]+sum_app$tTable[2,1]),format2) sum1.ll <- sapply(c(sum_app$tTable[1,1] - 1.96*sum_app$tTable[1,2], sum_app$tTable[3:5,1]+sum_app$tTable[1,1] - 1.96*sqrt(sum_app$tTable[3:5,2]^2 + sum_app$tTable[1,2]^2 + 2*lmm_app$varFix[3:5,1]), sum_app$tTable[2,1] - 1.96*sum_app$tTable[2,2], sum_app$tTable[6:8,1]+sum_app$tTable[2,1] - 1.96*sqrt(sum_app$tTable[6:8,2]^2 + sum_app$tTable[2,2]^2 + 2*lmm_app$varFix[6:8,2])),format2) sum1.ul <- sapply(c(sum_app$tTable[1,1] + 1.96*sum_app$tTable[1,2], sum_app$tTable[3:5,1]+sum_app$tTable[1,1] + 1.96*sqrt(sum_app$tTable[3:5,2]^2 + sum_app$tTable[1,2]^2 + 2*lmm_app$varFix[3:5,1]), sum_app$tTable[2,1] + 1.96*sum_app$tTable[2,2], sum_app$tTable[6:8,1]+sum_app$tTable[2,1] + 1.96*sqrt(sum_app$tTable[6:8,2]^2 + sum_app$tTable[2,2]^2 + 2*lmm_app$varFix[6:8,2])),format2) an1 <- anova(lmm_app) sum1.ests1 <- paste(sum1.est," (",sum1.ll,", ",sum1.ul,")",sep="") pval1 <- c(myformat4(an1$"p-value"[3]),"","","",myformat4(an1$"p-value"[4]),"","","") labels1 <- c("\\\\hspace{0.2in}None","\\\\hspace{0.2in}Little","\\\\hspace{0.2in}Normal","\\\\hspace{0.2in}Hungry","\\\\hspace{0.2in}Time $\\\\times$ None", "\\\\hspace{0.2in}Time $\\\\times$ Little","\\\\hspace{0.2in}Time $\\\\times$ Normal","\\\\hspace{0.2in}Time $\\\\times$ Hungry") lmm.table1 <- paste("\\\\hline ",label," \\\\\\\\",paste(paste(labels1,sum1.ests1,pval1,sep=" & "),collapse=" \\\\\\\\ "),sep="") ## HERE WE CHECK THE APPROPRIATENESS OF INCLUDING THE INTERACTION TERM, RETAINING THOSE WITH P < 0.10 (somewhat arbitrary level of significance) if(an1$"p-value"[4] > 0.1){ lmm_app <- lme(var ~ time + as.factor(appetite), data=longfup, random=~1|pepfarid, na.action=na.omit,method="ML") sum_app <- summary(lmm_app) sum1.est <- sapply(c(sum_app$tTable[1,1], sum_app$tTable[3:5,1]+sum_app$tTable[1,1], sum_app$tTable[2,1]),format2) sum1.ll <- sapply(c(sum_app$tTable[1,1] - 1.96*sum_app$tTable[1,2], sum_app$tTable[3:5,1]+sum_app$tTable[1,1] - 1.96*sqrt(sum_app$tTable[3:5,2]^2 + sum_app$tTable[1,2]^2 + 2*lmm_app$varFix[3:5,1]), sum_app$tTable[2,1] - 1.96*sum_app$tTable[2,2]),format2) sum1.ul <- sapply(c(sum_app$tTable[1,1] + 1.96*sum_app$tTable[1,2], sum_app$tTable[3:5,1]+sum_app$tTable[1,1] + 1.96*sqrt(sum_app$tTable[3:5,2]^2 + sum_app$tTable[1,2]^2 + 2*lmm_app$varFix[3:5,1]), sum_app$tTable[2,1] + 1.96*sum_app$tTable[2,2]),format2) an1 <- anova(lmm_app) sum1.ests1 <- paste(sum1.est," (",sum1.ll,", ",sum1.ul,")",sep="") pval1 <- c(myformat4(an1$"p-value"[3]),"","","",myformat4(an1$"p-value"[2])) labels1 <- c("\\\\hspace{0.2in}None","\\\\hspace{0.2in}Little","\\\\hspace{0.2in}Normal","\\\\hspace{0.2in}Hungry","\\\\hspace{0.2in}Time") lmm.table1 <- paste("\\\\hline ",label," \\\\\\\\",paste(paste(labels1,sum1.ests1,pval1,sep=" & "),collapse=" \\\\\\\\ "),sep="") } return(lmm.table1) } ene.applmm <- myapp_lmm(energy,"Energy (kcal)") car.applmm <- myapp_lmm(carbohydrate,"Carbohydrate (g)") pro.applmm <- myapp_lmm(TotalProtein,"Protein (g)") fat.applmm <- myapp_lmm(TotalFat,"Fat (g)") applmm.table <- paste(ene.applmm,car.applmm,pro.applmm,fat.applmm,sep=" \\\\\\\\ ") ## extract p-value for test of interaction for Protein and Fat lmm_app_proint <- anova(lme(TotalProtein ~ time*as.factor(appetite), data=longfup, random=~1|pepfarid, na.action=na.omit,method="ML")) lmm_app_fatint <- anova(lme(TotalFat ~ time*as.factor(appetite), data=longfup, random=~1|pepfarid, na.action=na.omit,method="ML")) @ \clearpage \subsection{Mixed Effects Models for Dietary Intake and Appetite} A mixed effects model is basically a generalization of the generalized linear model to allow repeated measures per subject. \begin{center} \begin{threeparttable} \footnotesize \caption{Fixed and Interaction Effects from Linear Mixed Models} \centering{ \begin{tabular}{lll} \hline & Effect (95\% CI) & P-value$^a$ \\ \Sexpr{applmm.table} \\ \hline \end{tabular} } \begin{tablenotes} \item[a] Wald F-tests for time effect (time + time by appetite interaction), appetite effect (appetite + time by appetite interaction), and test of interaction term. Non-significant interaction terms were dropped from the protein (p=\Sexpr{myformat2(lmm_app_proint$'p-value'[4])}) and fat (p=\Sexpr{myformat2(lmm_app_fatint$'p-value'[4])}) models. \item[b] Only one patient on one occasion gave ``very hungry'' response to appetite scale, this response was reclassified as ``hungry'' to avoid singularity in modeling. \end{tablenotes} \end{threeparttable} \end{center} % Hidden R code chunk --- reading in the data <>= ## CREATE PLOTS OF PREDICTED DIETARY INTAKE BY APPETITE FOR MANUSCRIPT USING LATTICEEXTRA ## library(latticeExtra) 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, col.line = 'black', ...) { x <- as.numeric(x) y <- as.numeric(y) ly <- as.numeric(ly[subscripts]) uy <- as.numeric(uy[subscripts]) panel.lines(x, uy, col = col.line, lty=5, lwd=1) panel.lines(x, ly, col = col.line, lty=5, lwd=1) panel.xyplot(x, y, pch = pch,col.line=col.line, ...) } longfup$APPETITE <- as.character(longfup$appetite) getpred <- function(var,int=TRUE){ longfup$var <- var if(int){lmm1 <- lme(var ~ time*as.factor(APPETITE), data=longfup, random=~1|pepfarid, na.action=na.omit,method="ML")} if(!int){lmm1 <- lme(var ~ time + as.factor(APPETITE), data=longfup, random=~1|pepfarid, na.action=na.omit,method="ML")} sum1 <- summary(lmm1) dat1 <- expand.grid(time=seq(0,90,length=91),APPETITE=as.factor(1:4)) dat1$linear.predictors <- predict(lmm1, dat1, level=c(0)) Designmat <- model.matrix(eval(eval(lmm1$call$fixed)[-2]),dat1[-3]) dat1$se.fit <- sqrt(diag(Designmat %*% lmm1$varFix %*% t(Designmat))) dat1$lower <- dat1$linear.predictors-1.96*dat1$se.fit dat1$upper <- dat1$linear.predictors+1.96*dat1$se.fit return(dat1) } ## PLOT PREDICTED ENERGY INTAKE OVER TIME (OUTPUT AS LMM_APP_ENERGY.PDF FOR LATEX INCLUSION) dat1 <- getpred(longfup$energy) energy1 <- xyplot(linear.predictors ~ time, groups=APPETITE,lwd=2, ly=dat1$lower,uy=dat1$upper,data=dat1,prepanel=prepanel.ci,panel=panel.superpose,panel.groups=panel.ci, type="l",col.line=c("#E41A1C","#377EB8","#4DAF4A","#EEC900"), aspect = 1, xlab="Time (days)",ylab="Energy (kcal)",scales=list(x=list(at=c(0,30,60,90), labels=c(0,30,60,90)))) energy1 <- energy1 + layer(panel.key(c("None", "Little", "Normal","Hungry"), 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(energy1, par.settings = custom.theme(symbol = c("#E41A1C","#377EB8","#4DAF4A","#EEC900"), fill = c("#E41A1C","#377EB8","#4DAF4A","#EEC900"), lwd=2)) ## REPEAT from dat1 on for carbohydrate, protein and fat @ \hspace{-0.13\textwidth} \begin{minipage}[b]{0.45\textwidth} \centering \includegraphics{graphics/lmm_app_energy} \end{minipage} \hspace{-0.1\textwidth} \begin{minipage}[b]{0.45\textwidth} \centering \includegraphics{graphics/lmm_app_carb} \end{minipage} \hspace{-0.1\textwidth} \begin{minipage}[b]{0.45\textwidth} \centering \includegraphics{graphics/lmm_app_protein} \end{minipage} \hspace{-0.1\textwidth} \begin{minipage}[b]{0.45\textwidth} \centering \includegraphics{graphics/lmm_app_fat} \end{minipage} % Hidden R code chunk --- reading in the data <>= ## COX REGRESSIONS WITH TIME-DEPENDENT COVARIATES ## # data includes one line per time interval (between visits), # so there are multiple lines per patient # note that 'V2' is a first principal component fit as follows: # V2 <- with(baseline, princomp(~age+sqrt(cd4)+bmi+h30+phosphate)) library(survival) library(Design) ################### ## TIME TO DEATH ## ################### dd<-datadist(COXTDDATA2) ## COXTDDATA2 does not have last observation carried forward options(datadist='dd') mod.ene1.cph <- cph(Surv(start.time,stop.time,event90)~sex+V2+energy,x=T,y=T,surv=T,data=COXTDDATA2) mysum.ene1 <- summary(mod.ene1.cph,energy=c(1000,1100))[c(2,6,4),c(4,6:7)] myanova.ene1 <- anova(mod.ene1.cph) mylab.ene1 <- "Energy (per 100 kcal increase)" mod.car1.cph <- cph(Surv(start.time,stop.time,event90)~sex+V2+carbohydrate,x=T,y=T,surv=T,data=COXTDDATA2) mysum.car1 <- summary(mod.car1.cph,carbohydrate=c(200,225))[c(2,6,4),c(4,6:7)] myanova.car1 <- anova(mod.car1.cph) mylab.car1 <- "Carbohydrate (per 25 g increase)" mod.pro1.cph <- cph(Surv(start.time,stop.time,event90)~sex+V2+TotalProtein,x=T,y=T,surv=T,data=COXTDDATA2) mysum.pro1 <- summary(mod.pro1.cph,TotalProtein=c(50,55))[c(2,6,4),c(4,6:7)] myanova.pro1 <- anova(mod.pro1.cph) mylab.pro1 <- "Protein (per 5 g increase)" mod.fat1.cph <- cph(Surv(start.time,stop.time,event90)~sex+V2+TotalFat,x=T,y=T,surv=T,data=COXTDDATA2) mysum.fat1 <- summary(mod.fat1.cph,TotalFat=c(50,55))[c(2,6,4),c(4,6:7)] myanova.fat1 <- anova(mod.fat1.cph) mylab.fat1 <- "Fat (per 5 g increase)" mod.app1.cph <- cph(Surv(start.time,stop.time,event90)~sex+V2+as.factor(appetite),x=T,y=T,surv=T,data=COXTDDATA2) mod.app1.cph.num <- cph(Surv(start.time,stop.time,event90)~sex+V2+as.numeric(appetite),x=T,y=T,surv=T,data=COXTDDATA2) mod.app1.cph.numvis <- cph(Surv(start.time,stop.time,event90)~sex+V2+as.numeric(appetite),x=T,y=T,surv=T,data=COXTDDATA2,subset=(visit %in% c(1,2,4,6))) mod.app1.cph.vis <- cph(Surv(start.time,stop.time,event90)~sex+V2+as.factor(appetite),x=T,y=T,surv=T,data=COXTDDATA2,subset=(visit %in% c(1,2,4,6))) ## all are similar mysum.app1 <- summary(mod.app1.cph,appetite=3)[c(2,4,6,8,10),c(4,6:7)] myanova.app1 <- anova(mod.app1.cph) mylab.app1 <- "Appetite" getTDest <- function(mysum,myanova,mylab){ estimates <- data.frame(c("1$^{\\\\mathrm{st}}$ Principal Component","Male",mylab),paste(sapply((mysum[,1]),myformat3), " (",sapply((mysum[,2]),myformat3),", ",sapply((mysum[,3]),myformat3),")",sep=""), sapply(myanova[c(2,1,3),3],myformat4)) print.latex <- paste(paste(estimates[,1],estimates[,2],estimates[,3],sep=" & "),collapse=" \\\\\\\\ ") return(print.latex) } estimates.ene1 <- getTDest(mysum.ene1,myanova.ene1,mylab.ene1) estimates.car1 <- getTDest(mysum.car1,myanova.car1,mylab.car1) estimates.pro1 <- getTDest(mysum.pro1,myanova.pro1,mylab.pro1) estimates.fat1 <- getTDest(mysum.fat1,myanova.fat1,mylab.fat1) labels.app1 <- c("1$^{\\\\mathrm{st}}$ Principal Component","Male","Appetite","\\\\hspace{0.2in}None vs. Normal","\\\\hspace{0.2in}Little vs. Normal","\\\\hspace{0.2in}Hungry vs. Normal") estimates <- paste(sapply((mysum.app1[,1]),myformat3)," (",sapply((mysum.app1[,2]),myformat3),", ",sapply((mysum.app1[,3]),myformat3),")",sep="") pvals.app1 <- sapply(myanova.app1[c(2,1,3),3],myformat4) estimates.app1a <- data.frame(labels.app1,c(estimates[1:2],"",estimates[3:5]),c(pvals.app1[1:3],"","","")) estimates.app1 <- paste(paste(estimates.app1a[,1],estimates.app1a[,2],estimates.app1a[,3],sep=" & "),collapse=" \\\\\\\\ ") ######################## ## TIME TO DEATH/LTFU ## ######################## mod.ene2.cph <- cph(Surv(start.time,stop.time,deadlost90)~sex+V2+energy,x=T,y=T,surv=T,data=COXTDDATA2) mysum.ene2 <- summary(mod.ene2.cph,energy=c(1000,1100))[c(2,6,4),c(4,6:7)] myanova.ene2 <- anova(mod.ene2.cph) mylab.ene2 <- "Energy (per 100 kcal increase)" mod.car2.cph <- cph(Surv(start.time,stop.time,deadlost90)~sex+V2+carbohydrate,x=T,y=T,surv=T,data=COXTDDATA2) mysum.car2 <- summary(mod.car2.cph,carbohydrate=c(200,225))[c(2,6,4),c(4,6:7)] myanova.car2 <- anova(mod.car2.cph) mylab.car2 <- "Carbohydrate (per 25 g increase)" mod.pro2.cph <- cph(Surv(start.time,stop.time,deadlost90)~sex+V2+TotalProtein,x=T,y=T,surv=T,data=COXTDDATA2) mysum.pro2 <- summary(mod.pro2.cph,TotalProtein=c(50,55))[c(2,6,4),c(4,6:7)] myanova.pro2 <- anova(mod.pro2.cph) mylab.pro2 <- "Protein (per 5 g increase)" mod.fat2.cph <- cph(Surv(start.time,stop.time,deadlost90)~sex+V2+TotalFat,x=T,y=T,surv=T,data=COXTDDATA2) mysum.fat2 <- summary(mod.fat2.cph,TotalFat=c(50,55))[c(2,6,4),c(4,6:7)] myanova.fat2 <- anova(mod.fat2.cph) mylab.fat2 <- "Fat (per 5 g increase)" mod.app2.cph <- cph(Surv(start.time,stop.time,deadlost90)~sex+V2+as.factor(appetite),x=T,y=T,surv=T,data=COXTDDATA2) mysum.app2 <- summary(mod.app2.cph,appetite=3)[c(2,4,6,8,10),c(4,6:7)] myanova.app2 <- anova(mod.app2.cph) mylab.app2 <- "Appetite" estimates.ene2 <- getTDest(mysum.ene2,myanova.ene2,mylab.ene2) estimates.car2 <- getTDest(mysum.car2,myanova.car2,mylab.car2) estimates.pro2 <- getTDest(mysum.pro2,myanova.pro2,mylab.pro2) estimates.fat2 <- getTDest(mysum.fat2,myanova.fat2,mylab.fat2) labels.app2 <- c("1$^{\\\\mathrm{st}}$ Principal Component","Male","Appetite","\\\\hspace{0.2in}None vs. Normal","\\\\hspace{0.2in}Little vs. Normal","\\\\hspace{0.2in}Hungry vs. Normal") estimates <- paste(sapply((mysum.app2[,1]),myformat3)," (",sapply((mysum.app2[,2]),myformat3),", ",sapply((mysum.app2[,3]),myformat3),")",sep="") pvals.app2 <- sapply(myanova.app2[c(2,1,3),3],myformat4) estimates.app2a <- data.frame(labels.app2,c(estimates[1:2],"",estimates[3:5]),c(pvals.app2[1:3],"","","")) estimates.app2 <- paste(paste(estimates.app2a[,1],estimates.app2a[,2],estimates.app2a[,3],sep=" & "),collapse=" \\\\\\\\ ") @ \subsection{Time-Dependent Cox Models for Dietary Intake, Appetite and 90 Day Mortality} Cox regression models with baseline and time-dependent covariates are used to assess the relationship between time to death, baseline variables and dietary intake or appetite. We use the first principal component of age, hemoglobin, CD4 count, BMI, and baseline phosphate marker to adjust for baseline status. Energy, carbohydrate, protein, and fat intake and appetite have been identified as post-ART time-dependent predictors of interest. Time intervals with no corresponding intake or appetite are excluded from analysis; rather than LOCF, the patient is dropped from the risk set for that interval.$^1$ \footnotetext[1]{Therneau, TM, Grambsch PM. {\it Modeling Survival Data: Extending the Cox Model}. New York, NY: Springer; 2000.} \begin{center} \begin{threeparttable} \footnotesize \caption{Adjusted Hazard Ratios for Mortality in First 90 Days - Exclude Missing} \centering{ \begin{tabular}{llll} \hline & HR (95\% CI) & P-value \\ \hline \\ \Sexpr{estimates.ene1} \\ \hline \Sexpr{estimates.car1} \\ \hline \Sexpr{estimates.pro1} \\ \hline \Sexpr{estimates.fat1} \\ \hline \Sexpr{estimates.app1} \\ \\ \hline \\ \end{tabular} } \end{threeparttable} \end{center} \subsection{Time-Dependent Cox Models for Dietary Intake, Appetite and 90 Day Mortality/LTFU} Cox regression models with baseline and time-dependent covariates are used to assess the relationship between time to death or LTFU, baseline variables and dietary intake or appetite. \begin{center} \begin{threeparttable} \footnotesize \caption{Adjusted Hazard Ratios for Mortality/LTFU in First 90 Days - Exclude Missing} \centering{ \begin{tabular}{llll} \hline & HR (95\% CI) & P-value \\ \hline \\ \Sexpr{estimates.ene2} \\ \hline \Sexpr{estimates.car2} \\ \hline \Sexpr{estimates.pro2} \\ \hline \Sexpr{estimates.fat2} \\ \hline \Sexpr{estimates.app2} \\ \\ \hline \\ \end{tabular} } \end{threeparttable} \end{center} <>= ############################################################# # BASELINE ONLY # ############################################################# S<-Surv(tx2lastknown,dead90) SL<-Surv(tx2lastknown,deadlost90) ################### ## TIME TO DEATH ## ################### library(survival) library(Design) dd<-datadist(with(baseline,data.frame(sex,V2,energy,carbohydrate,TotalProtein,TotalFat,appetite))) options(datadist='dd') mod.ene3.cph <- cph(S~sex+V2+energy,x=T,y=T,surv=T,data=baseline) mysum.ene3 <- summary(mod.ene3.cph,energy=c(1000,1100))[c(2,6,4),c(4,6:7)] myanova.ene3 <- anova(mod.ene3.cph) mylab.ene3 <- "Baseline Energy (per 100 kcal increase)" mod.car3.cph <- cph(S~sex+V2+carbohydrate,x=T,y=T,surv=T,data=baseline) mysum.car3 <- summary(mod.car3.cph,carbohydrate=c(200,225))[c(2,6,4),c(4,6:7)] myanova.car3 <- anova(mod.car3.cph) mylab.car3 <- "Baseline Carbohydrate (per 25 g increase)" mod.pro3.cph <- cph(S~sex+V2+TotalProtein,x=T,y=T,surv=T,data=baseline) mysum.pro3 <- summary(mod.pro3.cph,TotalProtein=c(50,55))[c(2,6,4),c(4,6:7)] myanova.pro3 <- anova(mod.pro3.cph) mylab.pro3 <- "Baseline Protein (per 5 g increase)" mod.fat3.cph <- cph(S~sex+V2+TotalFat,x=T,y=T,surv=T,data=baseline) mysum.fat3 <- summary(mod.fat3.cph,TotalFat=c(50,55))[c(2,6,4),c(4,6:7)] myanova.fat3 <- anova(mod.fat3.cph) mylab.fat3 <- "Baseline Fat (per 5 g increase)" mod.app3.cph <- cph(S~sex+V2+as.factor(appetite),x=T,y=T,surv=T,data=baseline) mod.app3.cph.num <- cph(S~sex+V2+as.numeric(appetite),x=T,y=T,surv=T,data=baseline) mod.app3.cph.numvis <- cph(S~sex+V2+as.numeric(appetite),x=T,y=T,surv=T,data=baseline,subset=(visit %in% c(1,2,4,6))) mod.app3.cph.vis <- cph(S~sex+V2+as.factor(appetite),x=T,y=T,surv=T,data=baseline,subset=(visit %in% c(1,2,4,6))) ## all are similar mysum.app3 <- summary(mod.app3.cph,appetite=3)[c(2,4,6,8,10),c(4,6:7)] myanova.app3 <- anova(mod.app3.cph) mylab.app3 <- "Baseline Appetite" getTDest <- function(mysum,myanova,mylab){ estimates <- data.frame(c("1$^{\\\\mathrm{st}}$ Principal Component","Male",mylab),paste(sapply((mysum[,1]),myformat3), " (",sapply((mysum[,2]),myformat3),", ",sapply((mysum[,3]),myformat3),")",sep=""), sapply(myanova[c(2,1,3),3],myformat4)) print.latex <- paste(paste(estimates[,1],estimates[,2],estimates[,3],sep=" & "),collapse=" \\\\\\\\ ") return(print.latex) } estimates.ene3 <- getTDest(mysum.ene3,myanova.ene3,mylab.ene3) estimates.car3 <- getTDest(mysum.car3,myanova.car3,mylab.car3) estimates.pro3 <- getTDest(mysum.pro3,myanova.pro3,mylab.pro3) estimates.fat3 <- getTDest(mysum.fat3,myanova.fat3,mylab.fat3) labels.app3 <- c("1$^{\\\\mathrm{st}}$ Principal Component","Male","Appetite","\\\\hspace{0.2in}None vs. Normal","\\\\hspace{0.2in}Little vs. Normal","\\\\hspace{0.2in}Hungry vs. Normal") estimates <- paste(sapply((mysum.app3[,1]),myformat3)," (",sapply((mysum.app3[,2]),myformat3),", ",sapply((mysum.app3[,3]),myformat3),")",sep="") pvals.app3 <- sapply(myanova.app3[c(2,1,3),3],myformat4) estimates.app3a <- data.frame(labels.app3,c(estimates[1:2],"",estimates[3:5]),c(pvals.app3[1:3],"","","")) estimates.app3 <- paste(paste(estimates.app3a[,1],estimates.app3a[,2],estimates.app3a[,3],sep=" & "),collapse=" \\\\\\\\ ") ######################## ## TIME TO DEATH/LTFU ## ######################## mod.ene4.cph <- cph(SL~sex+V2+energy,x=T,y=T,surv=T,data=baseline) mysum.ene4 <- summary(mod.ene4.cph,energy=c(1000,1100))[c(2,6,4),c(4,6:7)] myanova.ene4 <- anova(mod.ene4.cph) mylab.ene4 <- "Baseline Energy (per 100 kcal increase)" mod.car4.cph <- cph(SL~sex+V2+carbohydrate,x=T,y=T,surv=T,data=baseline) mysum.car4 <- summary(mod.car4.cph,carbohydrate=c(200,225))[c(2,6,4),c(4,6:7)] myanova.car4 <- anova(mod.car4.cph) mylab.car4 <- "Baseline Carbohydrate (per 25 g increase)" mod.pro4.cph <- cph(SL~sex+V2+TotalProtein,x=T,y=T,surv=T,data=baseline) mysum.pro4 <- summary(mod.pro4.cph,TotalProtein=c(50,55))[c(2,6,4),c(4,6:7)] myanova.pro4 <- anova(mod.pro4.cph) mylab.pro4 <- "Baseline Protein (per 5 g increase)" mod.fat4.cph <- cph(SL~sex+V2+TotalFat,x=T,y=T,surv=T,data=baseline) mysum.fat4 <- summary(mod.fat4.cph,TotalFat=c(50,55))[c(2,6,4),c(4,6:7)] myanova.fat4 <- anova(mod.fat4.cph) mylab.fat4 <- "Baseline Fat (per 5 g increase)" mod.app4.cph <- cph(SL~sex+V2+as.factor(appetite),x=T,y=T,surv=T,data=baseline) mysum.app4 <- summary(mod.app4.cph,appetite=3)[c(2,4,6,8,10),c(4,6:7)] myanova.app4 <- anova(mod.app4.cph) mylab.app4 <- "Appetite" estimates.ene4 <- getTDest(mysum.ene4,myanova.ene4,mylab.ene4) estimates.car4 <- getTDest(mysum.car4,myanova.car4,mylab.car4) estimates.pro4 <- getTDest(mysum.pro4,myanova.pro4,mylab.pro4) estimates.fat4 <- getTDest(mysum.fat4,myanova.fat4,mylab.fat4) labels.app4 <- c("1$^{\\\\mathrm{st}}$ Principal Component","Male","Appetite","\\\\hspace{0.2in}None vs. Normal","\\\\hspace{0.2in}Little vs. Normal","\\\\hspace{0.2in}Hungry vs. Normal") estimates <- paste(sapply((mysum.app4[,1]),myformat3)," (",sapply((mysum.app4[,2]),myformat3),", ",sapply((mysum.app4[,3]),myformat3),")",sep="") pvals.app4 <- sapply(myanova.app4[c(2,1,3),3],myformat4) estimates.app4a <- data.frame(labels.app4,c(estimates[1:2],"",estimates[3:5]),c(pvals.app4[1:3],"","","")) estimates.app4 <- paste(paste(estimates.app4a[,1],estimates.app4a[,2],estimates.app4a[,3],sep=" & "),collapse=" \\\\\\\\ ") @ \subsection{Cox Models for Dietary Intake, Appetite and 90 Day Mortality} Cox regression models with baseline covariates are used to assess the relationship between time to death, baseline variables and baseline dietary intake or appetite. We use the first principal component of age, hemoglobin, CD4 count, BMI, and baseline phosphate marker to adjust for baseline status. Baseline energy, carbohydrate, protein, and fat intake and appetite have been identified as predictors of interest. \begin{center} \begin{threeparttable} \footnotesize \caption{Adjusted Hazard Ratios for Mortality in First 90 Days} \centering{ \begin{tabular}{llll} \hline & HR (95\% CI) & P-value \\ \hline \\ \Sexpr{estimates.ene3} \\ \hline \Sexpr{estimates.car3} \\ \hline \Sexpr{estimates.pro3} \\ \hline \Sexpr{estimates.fat3} \\ \hline \Sexpr{estimates.app3} \\ \\ \hline \\ \end{tabular} } \end{threeparttable} \end{center} \subsection{Cox Models for Dietary Intake, Appetite and 90 Day Mortality/LTFU} Cox regression models with baseline covariates are used to assess the relationship between time to death or LTFU, baseline variables and baseline dietary intake or appetite. \begin{center} \begin{threeparttable} \footnotesize \caption{Adjusted Hazard Ratios for Mortality/LTFU in First 90 Days} \centering{ \begin{tabular}{llll} \hline & HR (95\% CI) & P-value \\ \hline \\ \Sexpr{estimates.ene4} \\ \hline \Sexpr{estimates.car4} \\ \hline \Sexpr{estimates.pro4} \\ \hline \Sexpr{estimates.fat4} \\ \hline \Sexpr{estimates.app4} \\ \\ \hline \\ \end{tabular} } \end{threeparttable} \end{center} \end{document}