\documentclass[10pt]{article} \bibliographystyle{cj} \usepackage{longtable} % allows tables to break \usepackage[pdftex]{lscape} % allows tables to be landscaped %\usepackage[pdftex]{graphicx} \usepackage{pdfpages} \usepackage[letterpaper,margin=1in]{geometry} \usepackage{setspace,relsize} % needed for latex(describe()), \code \usepackage{moreverb} % for verbatimtabinput % \usepackage{tabularx} %\usepackage{colortbl} %\usepackage{placeins} %\usepackage[dotinlabels]{titletoc} %\usepackage{amsmath} %\usepackage{amsfonts} \usepackage{amssymb} \usepackage{fancyhdr} \usepackage{graphicx, subfigure} %\usepackage{vmargin} \usepackage{lscape} \usepackage{rotating} \usepackage{endnotes} \usepackage{color} \usepackage{enumerate, hyperref, pdflscape, xcolor} \hypersetup{ bookmarksnumbered=false, % true means bookmarks in left window are numbered bookmarksopen=false, % true means only level 1 is displayed. colorlinks=true, % false: boxed links; true: colored links linkcolor=blue, % color of internal links } \def\linkcol{blue} % \usepackage{tabularx} % required in the preamble %\setpapersize{USletter} % \usepackage{tabularx} % required in the preamble %\setpapersize{USletter} \begin{document} %Define all graphics to be 100\% of their width \setkeys{Gin}{width=1.0\textwidth} %default is 0.8 \center{The association of long-term CD4 count and BMI} \center{Report by: Cathy A. Jenkins} \vspace{0.1in} \tableofcontents \flushleft \section{Introduction} This report summarizes the results from the analysis investigating the relationship between BMI and long-term CD4 count using longitudinal mixed effects models. \section{Methods} Subjects were included in the analysis if they were $\ge$ 18, ART-naive (defined as no prior record of any antiretroviral agent exposure) and if they initiated combination ART (cART) treatment (defined as a combined regimen of $\ge$ 3 antiretroviral agents) after 1998. At least one CD4 count and BMI measurement subsequent to cART initiation had to be collected. \vspace{0.1in} The initial analysis included all comers and modeled the association of long-term CD4 count (square root transformed) with BMI adjusting for time since cART initiation, age at cART initiation, sex, race (white vs nonwhite), initial cART regimen class (NNRTI-based, PI-based, NRTI-only, and Other), year of cART start, baseline CD4 count and log10-transformed HIV1-RNA, and cohort. Interactions between BMI and time, sex, and race as well as between sex and race were included in the primary analysis. The significance of the interactions were evaluated to determine whether they should be included in the final analysis. All continuous covariates were fit with restricted cubic splines (6 knots). All records with a non-missing CD4 count were included in the analysis provided a BMI measurement within the window -180 days to +180 days of the CD4 count measurement could be identified. Repeat use of BMI measurements for different CD4 counts was also recorded. \vspace{0.1in} Several secondary analyses were also performed. These included: \begin{itemize} \item subsetting the cohort to those who were virologically suppressed $>$ 50\% of their follow-up time using the same method of determining time virologically suppressed that was used in previous analyses, \item creating a cohort of those virologically suppressed $>$ 50\% of their follow-up time simply by counting the number of suppressed and unsuppressed HIV1-RNA records and defining those with more suppressed records as being virologically suppressed $>$ 50\% of the time, \item censoring follow-up time for each subject at the viral load measurement prior to virologic failure as defined in CCASAnet analyses. In these analyses, someone was deemed as having virologically failed if they had one measurement $>$ 1000 or two consecutive ones that were detectable. In the first instance, the failure was said to have occurred when the viral load exceeded 1000; in the second instance, the failure was said to have occurred at the second of the two consecutive detectable measurements. \end{itemize} \vspace{0.1in} Because the data included repeated measures on subjects, the correlation structure needed to be estimated to ensure appropriate inference was made. Variograms were created and autocorrelation funtions were plotted to help decide which correlation structure was most appropriate. For completeness, several longitudinal methods were used along with several correlation structures. Longitudinal methods used included generalized least squares (gls), generalized estimating equations (gee), and linear mixed effects models (lme). Correlation structures applied to these models included compound symmetric in which every measure within an individual is equally correlated and an autoregressive correlation structure (AR1) in which repeated measures within an individual depend linearly on their own previous values. The final analysis reports results from the lme model using AR1 correlation. For completeness, results from the gls model with AR1 correlation were also reported. %------------------------------------------------------------------------% % BMI versus long-term CD4 count % %------------------------------------------------------------------------% <>= rm(list=ls()) options(width=180) opts_chunk$set(error=FALSE, warning=FALSE, message=FALSE) #----------------------------------------# # Libraries # #----------------------------------------# library(Hmisc,quietly=TRUE) library(rms,quietly=TRUE) library(nlme,quietly=TRUE) library(ggplot2,quietly=TRUE) library(lattice,quietly=TRUE) library(cmprsk,quietly=TRUE) library(sas7bdat) library(foreign) #----------------------------------------# # Loading data # #----------------------------------------# #BL edits: data <- read.csv('N1107excerpt_2014.csv', stringsAsFactors=FALSE, na.strings=c(NA,'')) # data <- read.csv('~/InfectiousDiseases/JohnKoethe/MSCIProject/NA-ACCORD data/N1107excerpt_2014.csv', stringsAsFactors=FALSE, na.strings=c(NA,'')) data_w_10_pct_chg <- read.csv('~/InfectiousDiseases/JohnKoethe/MSCIProject/NA-ACCORD data/Pregnant_women_as_10pctBMI_CH.csv', stringsAsFactors=FALSE, na.strings=c(NA,'')) #end BL edits names(data) <- tolower(names(data)) names(data_w_10_pct_chg) <- tolower(names(data_w_10_pct_chg)) # #BL edits: # ###################################################################### # non-excerpt data already has cohort # cohort.data <- read.csv('~/InfectiousDiseases/JohnKoethe/MSCIProject/NA-ACCORD data/n1107cohort.csv', stringsAsFactors=FALSE, na.strings=c(NA,'')) # # cohort.data <- read.csv('n1107cohort.csv', stringsAsFactors=FALSE, na.strings=c(NA,'')) # # names(cohort.data) <- tolower(names(cohort.data)) # # # Merging in cohort data # data$cohort <- cohort.data[match(data$newid, cohort.data$newid,nomatch=NA),'cohort'] # Creating time variables from date variables in additional data data_w_10_pct_chg$dt_art_stop <- as.Date(data_w_10_pct_chg$dt_art_stop, format='%m/%d/%Y') data_w_10_pct_chg$dt_art_start <- as.Date(data_w_10_pct_chg$dt_art_start, format='%m/%d/%Y') data_w_10_pct_chg$t_artstop <- with(data_w_10_pct_chg, as.numeric(dt_art_stop - dt_art_start)/365.25) data_w_10_pct_chg$dt_lastact <- as.Date(data_w_10_pct_chg$dt_lastact, format='%m/%d/%Y') data_w_10_pct_chg$t_lastact <- with(data_w_10_pct_chg, as.numeric(dt_lastact - dt_art_start)/365.25) data_w_10_pct_chg$dt_death <- as.Date(data_w_10_pct_chg$dt_death, format='%m/%d/%Y') data_w_10_pct_chg$t_dth <- with(data_w_10_pct_chg, as.numeric(dt_death - dt_art_start)/365.25) data_w_10_pct_chg$dt_cd4_bs <- as.Date(data_w_10_pct_chg$dt_cd4_bs, format='%m/%d/%Y') data_w_10_pct_chg$t_cd4_bs <- with(data_w_10_pct_chg, as.numeric(dt_cd4_bs - dt_art_start)/365.25) data_w_10_pct_chg$dt_vload_bs <- as.Date(data_w_10_pct_chg$dt_vload_bs, format='%m/%d/%Y') data_w_10_pct_chg$t_vload_bs <- with(data_w_10_pct_chg, as.numeric(dt_vload_bs - dt_art_start)/365.25) data_w_10_pct_chg$dt_bmi_bs <- as.Date(data_w_10_pct_chg$dt_bmi_bs, format='%m/%d/%Y') data_w_10_pct_chg$t_bmi_bs <- with(data_w_10_pct_chg, as.numeric(dt_bmi_bs - dt_art_start)/365.25) data_w_10_pct_chg$dt_lab <- as.Date(data_w_10_pct_chg$dt_lab, format='%m/%d/%Y') data_w_10_pct_chg$t_lab <- with(data_w_10_pct_chg, as.numeric(dt_lab - dt_art_start)/365.25) data_w_10_pct_chg <- upData(data_w_10_pct_chg, drop=c('dt_art_start','dt_art_stop','dt_lastact', 'dt_death','dt_cd4_bs','dt_vload_bs','dt_bmi_bs', 'dt_lab'), rename=c('naid'='newid')) data_total <- merge(data, data_w_10_pct_chg, all=TRUE) data_total_female <- subset(data_total, sex == 2) # end BL edits # ###################################################################### @ <>= # For assessing virologic failure for the > 90% suppressed cohort # Note: We want to censor at the age prior to the failure so must find the record just before it source('virologic_failure_function_adapted.R') source('model_fit_assessment.R') # This function is to do the necessary data management for John's descriptive graph based on the data of CD4 over time by BMI categories descriptive_graph_mgmt <- function(data){ # John wants a graph of CD4 over time with separate lines by BMI. # To accomplish, need to: # 1) Find CD4 counts at year 1, 2, 3, 4 # 2) Categorize the person's BMI at year 1, 2, 3, 4 # 3) Find the mean CD4 for a given BMI category/year group # 4) Plot the lines by BMI category cd4_by_year <- lapply(split(subset(data, select=c(cfar_pid,t_lab,bmi_closest,cd4_count)),data$cfar_pid),FUN=function(y){ # data_cat <- ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest']< 18.5,0, # ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest'] >= 18.5 & y[,'bmi_closest'] < 25.0,1, # ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest'] >= 25.0 & y[,'bmi_closest'] < 30.0,2, # ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest'] >= 30.0 & y[,'bmi_closest'] < 40.0,3, # ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest'] >= 40.0,4,NA))))) # data_cat.factor <- factor(data_cat, levels=c(0,1,2,3,4), labels=c('Underweight','Normal','Overweight','Obese','Morbidly obese')) data_cat <- ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest']< 18.5,0, ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest'] >= 18.5 & y[,'bmi_closest'] <= 22.0,1, ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest'] > 22.0 & y[,'bmi_closest'] < 25.0,2, ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest'] >= 25.0 & y[,'bmi_closest'] < 30.0,3, ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest'] >= 30.0 & y[,'bmi_closest'] < 40.0,4, ifelse(!is.na(y[,'bmi_closest']) & y[,'bmi_closest'] >= 40.0,5,NA)))))) data_cat.factor <- factor(data_cat, levels=c(0,1,2,3,4,5), labels=c('< 18.5','18,5-22.0','22.1-24.9','25.0-29.9','30.0-39.9', '>= 40.0')) diff1 <- y[,'t_lab'] - 1 diff2 <- y[,'t_lab'] - 2 diff3 <- y[,'t_lab'] - 3 diff4 <- y[,'t_lab'] - 4 diff5 <- y[,'t_lab'] - 5 diff1 <- ifelse(abs(diff1) > 0.5,NA,diff1) diff2 <- ifelse(abs(diff2) > 0.5,NA,diff2) diff3 <- ifelse(abs(diff3) > 0.5,NA,diff3) diff4 <- ifelse(abs(diff4) > 0.5,NA,diff4) diff5 <- ifelse(abs(diff5) > 0.5,NA,diff5) which.1 <- which(abs(diff1) == abs(min(diff1,na.rm=TRUE)))[1] which.2 <- which(abs(diff2) == abs(min(diff2,na.rm=TRUE)))[1] which.3 <- which(abs(diff3) == abs(min(diff3,na.rm=TRUE)))[1] which.4 <- which(abs(diff4) == abs(min(diff4,na.rm=TRUE)))[1] which.5 <- which(abs(diff5) == abs(min(diff5,na.rm=TRUE)))[1] bmi_year1 <- bmi_year2 <- bmi_year3 <- bmi_year4 <- NA cd4_year1 <- cd4_year2 <- cd4_year3 <- cd4_year4 <- NA bmi_cat_year1 <- bmi_cat_year2 <- bmi_cat_year3 <- bmi_cat_year4 <- NA if(length(which.1) > 0){ bmi_year1 <- y[which.1,'bmi_closest'] cd4_year1 <- y[which.1,'cd4_count'] bmi_cat_year1 <- data_cat.factor[which.1] } if(length(which.2) > 0){ bmi_year2 <- y[which.2,'bmi_closest'] cd4_year2 <- y[which.2,'cd4_count'] bmi_cat_year2 <- data_cat.factor[which.2] } if(length(which.3) > 0){ bmi_year3 <- y[which.3,'bmi_closest'] cd4_year3 <- y[which.3,'cd4_count'] bmi_cat_year3 <- data_cat.factor[which.3] } if(length(which.4) > 0){ bmi_year4 <- y[which.4,'bmi_closest'] cd4_year4 <- y[which.4,'cd4_count'] bmi_cat_year4 <- data_cat.factor[which.4] } if(length(which.5) > 0){ bmi_year5 <- y[which.5,'bmi_closest'] cd4_year5 <- y[which.5,'cd4_count'] bmi_cat_year5 <- data_cat.factor[which.5] } df <- data.frame('cfar_pid'=y[1,'cfar_pid'], 'bmi_year1'=bmi_year1, 'bmi_year2'=bmi_year2, 'bmi_year3'=bmi_year3, 'bmi_year4'=bmi_year4, 'bmi_year5'=bmi_year5, 'bmi_cat_year1'=bmi_cat_year1, 'bmi_cat_year2'=bmi_cat_year2, 'bmi_cat_year3'=bmi_cat_year3, 'bmi_cat_year4'=bmi_cat_year4, 'bmi_cat_year5'=bmi_cat_year5, 'cd4_year1'=cd4_year1, 'cd4_year2'=cd4_year2, 'cd4_year3'=cd4_year3, 'cd4_year4'=cd4_year4, 'cd4_year5'=cd4_year5) return(df) }) cd4_by_year_df <- do.call('rbind',cd4_by_year) data$bmi_at_age_fhaart_special <- ifelse(data$bmi_at_age_fhaart < 18.5,0, ifelse(data$bmi_at_age_fhaart >= 18.5 & data$bmi_at_age_fhaart <= 22.0,1, ifelse(data$bmi_at_age_fhaart > 22.0 & data$bmi_at_age_fhaart < 25.0,2, ifelse(data$bmi_at_age_fhaart >= 25.0 & data$bmi_at_age_fhaart < 30.0,3, ifelse(data$bmi_at_age_fhaart >= 30.0 & data$bmi_at_age_fhaart < 40.0,4, ifelse(data$bmi_at_age_fhaart >= 40.0,5,NA)))))) data$bmi_at_age_fhaart_special.factor <- factor(data$bmi_at_age_fhaart_special, levels=c(0,1,2,3,4,5), labels=c('< 18.5','18.5-22.0','22.1-24.9','25.0-29.9','30.0-39.9','>= 40.0')) cd4_by_year_combined0 <- tapply(data$cd4_count_at_age_fhaart[!duplicated(data$cfar_pid)], INDEX=data$bmi_at_age_fhaart_special.factor[!duplicated(data$cfar_pid)], FUN=function(y){mean(y,na.rm=TRUE)}) cd4_by_year_combined1 <- tapply(cd4_by_year_df$cd4_year1, INDEX=cd4_by_year_df$bmi_cat_year1,FUN=function(y){ mean_cd4 <- mean(y,na.rm=TRUE) }) cd4_by_year_combined2 <- tapply(cd4_by_year_df$cd4_year2, INDEX=cd4_by_year_df$bmi_cat_year2,FUN=function(y){ mean_cd4 <- mean(y,na.rm=TRUE) }) cd4_by_year_combined3 <- tapply(cd4_by_year_df$cd4_year3, INDEX=cd4_by_year_df$bmi_cat_year3,FUN=function(y){ mean_cd4 <- mean(y,na.rm=TRUE) }) cd4_by_year_combined4 <- tapply(cd4_by_year_df$cd4_year4, INDEX=cd4_by_year_df$bmi_cat_year4,FUN=function(y){ mean_cd4 <- mean(y,na.rm=TRUE) }) cd4_by_year_combined5 <- tapply(cd4_by_year_df$cd4_year5, INDEX=cd4_by_year_df$bmi_cat_year5,FUN=function(y){ mean_cd5 <- mean(y,na.rm=TRUE) }) cd4_by_year_combinednew1 <- data.frame('Time'=rep(c(0,1,2,3,4,5),each=6), 'CD4_count'=c(cd4_by_year_combined0,cd4_by_year_combined1, cd4_by_year_combined2,cd4_by_year_combined3,cd4_by_year_combined4,cd4_by_year_combined5), # 'BMI'=rep(c('Underweight','Normal','Overweight','Obese','Morbidly obese'),times=5)) 'BMI'=rep(c('< 18.5','18.5-22.0','22.1-24.9','25.0-29.9','30.0-39.9','>= 40.0'),times=6)) cd4_by_year_combinednew1$BMI <- factor(cd4_by_year_combinednew1$BMI, # levels=c('Underweight','Normal','Overweight','Obese','Morbidly obese')) levels=c('< 18.5','18.5-22.0','22.1-24.9','25.0-29.9','30.0-39.9','>= 40.0')) return(cd4_by_year_combinednew1) } # end of function cd4_by_bmi_figure <- function(mod_predict,predict_data, fig_limits, plot_method, subset_method){ if(subset_method == 'overall'){ rge_wm <- fig_limits[['bmi_rge_wm']] rge_nwm <- fig_limits[['bmi_rge_nwm']] rge_wf <- fig_limits[['bmi_rge_wf']] rge_nwf <- fig_limits[['bmi_rge_nwf']] } else if(subset_method == 'lme'){ rge_wm <- fig_limits[['tmp_rge_wm']] rge_nwm <- fig_limits[['tmp_rge_nwm']] rge_wf <- fig_limits[['tmp_rge_wf']] rge_nwf <- fig_limits[['bmi_rge_nwf']] } else if(subset_method == 'wo_outliers'){ rge_wm <- fig_limits[['tmp_w_residuals_rge_wm']] rge_nwm <- fig_limits[['tmp_w_residuals_rge_nwm']] rge_wf <- fig_limits[['tmp_w_residuals_rge_wf']] rge_nwf <- fig_limits[['tmp_w_residuals_rge_nwf']] } overall_x_rge <- fig_limits[['overall_x_rge']] if(plot_method == 'ggplot2'){ predict_data$predicted_cd4 <- mod_predict predict_data$Group <- paste(predict_data$race.factor, predict_data$sex.factor,sep=' ') predict_data$Group <- factor(predict_data$Group, levels=c('white male','nonwhite male','white female','nonwhite female'), labels=c('White male','Nonwhite male','White female','Nonwhite female')) predict_data$sex.factor <- factor(predict_data$sex.factor, levels=c('male','female')) predict_data$race.factor <- factor(predict_data$race.factor, levels=c('white','nonwhite')) ## Not sure how to only plot each line on the range of BMIs in the data... ggplot(subset(predict_data, bmi_closest <= 50.0), aes(x=bmi_closest, y=predicted_cd4, colour=Group)) + geom_line(data=predict_data[which(predict_data$sex.factor == 'male' & predict_data$race.factor == 'white' & predict_data$bmi_closest >= rge_wm[1] & predict_data$bmi_closest <= min(c(rge_wm[2],50))),]) + geom_line(data=predict_data[which(predict_data$sex.factor == 'male' & predict_data$race.factor == 'nonwhite' & predict_data$bmi_closest >= rge_nwm[1] & predict_data$bmi_closest <= min(c(rge_nwm[2],50))),]) + geom_line(data=predict_data[which(predict_data$sex.factor == 'female' & predict_data$race.factor == 'white' & predict_data$bmi_closest >= rge_wf[1] & predict_data$bmi_closest <= min(c(rge_wf[2],50))),]) + geom_line(data=predict_data[which(predict_data$sex.factor == 'female' & predict_data$race.factor == 'nonwhite' & predict_data$bmi_closest >= rge_nwf[1] & predict_data$bmi_closest <= min(c(rge_nwf[2],50))),]) + scale_colour_manual(values=c('green','red','blue','black'), name='Group', breaks=levels(predict_data$Group)) + ylim(0,1200) + xlab('BMI') + ylab('CD4 count') } else if(plot_method == 'base'){ # White males plot(mod_predict[which(predict_data$sex.factor == 'male' & predict_data$race.factor == 'white' & predict_data$bmi_closest >= rge_wm[1] & predict_data$bmi_closest <= min(c(rge_wm[2],50)))] ~ predict_data$bmi_closest[which(predict_data$sex.factor == 'male' & predict_data$race.factor == 'white' & predict_data$bmi_closest >= rge_wm[1] & predict_data$bmi_closest <= min(c(rge_wm[2],50)))], type='l',ylim=rge,ylab='CD4 count',xlab='BMI',xlim=c(10,50)) # Non-white males lines(mod_predict[which(predict_data$sex.factor == 'male' & predict_data$race.factor == 'nonwhite' & predict_data$bmi_closest >= rge_nwm[1] & predict_data$bmi_closest <= min(c(rge_nwm[2],50)))] ~ predict_data$bmi_closest[which(predict_data$sex.factor == 'male' & predict_data$race.factor == 'nonwhite' & predict_data$bmi_closest >= rge_nwm[1] & predict_data$bmi_closest <= min(c(rge_nwm[2],50)))],col='red') # White females lines(mod_predict[which(predict_data$sex.factor == 'female' & predict_data$race.factor == 'white' & predict_data$bmi_closest >= rge_wf[1] & predict_data$bmi_closest <= min(c(rge_wf[2],50)))] ~ predict_data$bmi_closest[which(predict_data$sex.factor == 'female' & predict_data$race.factor == 'white' & predict_data$bmi_closest >= rge_wf[1] & predict_data$bmi_closest <= min(c(rge_wf[2],50)))],col='blue') # Non-white females lines(mod_predict[which(predict_data$sex.factor == 'female' & predict_data$race.factor == 'nonwhite' & predict_data$bmi_closest >= rge_nwf[1] & predict_data$bmi_closest <= min(c(rge_nwf[2],50)))] ~ predict_data$bmi_closest[which(predict_data$sex.factor == 'female' & predict_data$race.factor == 'nonwhite' & predict_data$bmi_closest >= rge_nwf[1] & predict_data$bmi_closest <= min(c(rge_nwf[2],50)))],col='green') legend('topleft',fill=c('black','red','blue','green'),legend=c('White male','Non-white male','White female','Non-white female')) } }# end of function cd4_by_time_model_based_figure <- function(data, mod, method,main){ med_age <- median(data$age_fhaart[which(!duplicated(data$cfar_pid) & abs(data$residuals) <= 6)],na.rm=TRUE) med_cd4 <- median(data$cd4_count_at_age_fhaart_sqrt[which(!duplicated(data$cfar_pid) & abs(data$residuals) <= 6)],na.rm=TRUE) med_vl <- median(data$loghiv1_rna_at_age_fhaart[which(!duplicated(data$cfar_pid) & abs(data$residuals) <= 6)],na.rm=TRUE) mode_of_cohort_tbl <- sort(table(data$cohort),decreasing=TRUE) mode_of_cohort <- names(mode_of_cohort_tbl)[1] cd4_by_time_model_based <- data.frame('t_lab'=rep(seq(0,5,length.out=100),times=5), 'bmi_closest'=rep(c(18.5,22.0,25.0,30.0,40.0),each=100), 'sex.factor'=rep('male',times=500), 'race.factor'=rep('nonwhite',times=500), 'age_fhaart'=rep(med_age,times=500), 'regimen_class_in_fhaart.factor'=rep('NNRTI-based',times=500), 'year_of_first_haart_start'=rep(median(data$year_of_first_haart_start[which(!duplicated(data$cfar_pid) & abs(data$residuals) <= 6)]),times=500), 'cd4_count_at_age_fhaart_sqrt'=rep(med_cd4,times=500), 'loghiv1_rna_at_age_fhaart'=rep(med_vl,times=500), 'cohort'=rep(mode_of_cohort,times=500)) # 'cohort'=rep(sort(unique(data$cohort))[1],times=500)) mod_predict <- predict(mod, newdata=cd4_by_time_model_based,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed mod_predict <- mod_predict^2 if(method == 'base'){ rge_lme <- range(mod_predict) # x.rge_lme <- range(bmi$bmi_closest,na.rm=TRUE) # x.rge_lme <- c(x.rge_lme[1],quantile(bmi$bmi_closest,prob=0.90,na.rm=TRUE)) # BMI of 18.5 plot(mod_predict[which(cd4_by_time_model_based$bmi_closest == 18.5)] ~ cd4_by_time_model_based$t_lab[which(cd4_by_time_model_based$bmi_closest == 18.5)],type='l',ylim=c(0,800), ylab='CD4 count',xlab='Time since ART initiation (years)',xlim=c(0,5)) # BMI of 22.0 lines(mod_predict[which(cd4_by_time_model_based$bmi_closest == 22.0)] ~ cd4_by_time_model_based$t_lab[which(cd4_by_time_model_based$bmi_closest == 22.0)],col='red') # BMI of 25.0 lines(mod_predict[which(cd4_by_time_model_based$bmi_closest == 25.0)] ~ cd4_by_time_model_based$t_lab[which(cd4_by_time_model_based$bmi_closest == 25.0)],col='green') # BMI of 30.0 lines(mod_predict[which(cd4_by_time_model_based$bmi_closest == 30.0)] ~ cd4_by_time_model_based$t_lab[which(cd4_by_time_model_based$bmi_closest == 30.0)],col='blue') # BMI of 40.0 lines(mod_predict[which(cd4_by_time_model_based$bmi_closest == 40.0)] ~ cd4_by_time_model_based$t_lab[which(cd4_by_time_model_based$bmi_closest == 40.0)],col='purple') # legend('topleft',legend=c('18.5','22.0','25.0','30.0','40.0'),fill=c('black','red','green','blue','purple'),cex=0.75) legend('topleft',legend=c('40.0','30.0','25.0','22.0','18.5'),fill=c('purple','blue','green','red','black'),cex=0.75) } else if(method == 'ggplot2') { cd4_by_time_model_based$predicted_cd4 <- mod_predict cd4_by_time_model_based$bmi_closest <- factor(cd4_by_time_model_based$bmi_closest, levels=c(18.5,22,25,30,40)) ggplot(cd4_by_time_model_based, aes(x=t_lab, y=predicted_cd4,group=bmi_closest, colour=bmi_closest)) + geom_line() + ylab('Predicted CD4') + xlab('Time since ART initiation (years)') + scale_colour_manual(values=c('black','red','green','blue','purple'), name='BMI') } } # end of function predicted_cd4_function <- function(data, mod){ med_age <- median(data$age_fhaart[which(!duplicated(data$cfar_pid) & abs(data$residuals) <= 6)],na.rm=TRUE) med_cd4 <- median(data$cd4_count_at_age_fhaart_sqrt[which(!duplicated(data$cfar_pid) & abs(data$residuals) <= 6)],na.rm=TRUE) med_vl <- median(data$loghiv1_rna_at_age_fhaart[which(!duplicated(data$cfar_pid) & abs(data$residuals) <= 6)],na.rm=TRUE) mode_of_cohort_tbl <- sort(table(data$cohort),decreasing=TRUE) mode_of_cohort <- names(mode_of_cohort_tbl)[1] cd4_by_time_model_based <- data.frame('t_lab'=rep(seq(0,5,length.out=100),times=5), 'bmi_closest'=rep(c(18.5,22.0,25.0,30.0,40.0),each=100), 'sex.factor'=rep('male',times=500), 'race.factor'=rep('nonwhite',times=500), 'age_fhaart'=rep(med_age,times=500), 'regimen_class_in_fhaart.factor'=rep('NNRTI-based',times=500), 'year_of_first_haart_start'=rep(median(data$year_of_first_haart_start[which(!duplicated(data$cfar_pid) & abs(data$residuals) <= 6)]),times=500), 'cd4_count_at_age_fhaart_sqrt'=rep(med_cd4,times=500), 'loghiv1_rna_at_age_fhaart'=rep(med_vl,times=500), 'cohort'=rep(mode_of_cohort,times=500)) # 'cohort'=rep(sort(unique(data$cohort))[1],times=500)) mod_predict <- predict(mod, newdata=cd4_by_time_model_based,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed mod_predict <- mod_predict^2 # Calculating predicted CD4 values for years 0,1,3, and 5 for each of the BMI cut-offs: year0_18.5 <- mod_predict[which(cd4_by_time_model_based$bmi_closest == 18.5 & cd4_by_time_model_based$t_lab == 0.0)][1] year0_22 <- mod_predict[which(cd4_by_time_model_based$bmi_closest == 22.0 & cd4_by_time_model_based$t_lab == 0.0)][1] year0_25 <- mod_predict[which(cd4_by_time_model_based$bmi_closest == 25.0 & cd4_by_time_model_based$t_lab == 0.0)][1] year0_30 <- mod_predict[which(cd4_by_time_model_based$bmi_closest == 30.0 & cd4_by_time_model_based$t_lab == 0.0)][1] year0_40 <- mod_predict[which(cd4_by_time_model_based$bmi_closest == 40.0 & cd4_by_time_model_based$t_lab == 0.0)][1] year1_18.5 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 18.5 & round(cd4_by_time_model_based$t_lab,1) == 1.0))][1] year1_22 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 22.0 & round(cd4_by_time_model_based$t_lab,1) == 1.0))][1] year1_25 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 25.0 & round(cd4_by_time_model_based$t_lab,1) == 1.0))][1] year1_30 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 30.0 & round(cd4_by_time_model_based$t_lab,1) == 1.0))][1] year1_40 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 40.0 & round(cd4_by_time_model_based$t_lab,1) == 1.0))][1] year3_18.5 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 18.5 & round(cd4_by_time_model_based$t_lab,1) == 3.0))][1] year3_22 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 22.0 & round(cd4_by_time_model_based$t_lab,1) == 3.0))][1] year3_25 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 25.0 & round(cd4_by_time_model_based$t_lab,1) == 3.0))][1] year3_30 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 30.0 & round(cd4_by_time_model_based$t_lab,1) == 3.0))][1] year3_40 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 40.0 & round(cd4_by_time_model_based$t_lab,1) == 3.0))][1] year5_18.5 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 18.5 & round(cd4_by_time_model_based$t_lab,1) == 5.0))][1] year5_22 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 22.0 & round(cd4_by_time_model_based$t_lab,1) == 5.0))][1] year5_25 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 25.0 & round(cd4_by_time_model_based$t_lab,1) == 5.0))][1] year5_30 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 30.0 & round(cd4_by_time_model_based$t_lab,1) == 5.0))][1] year5_40 <- mod_predict[max(which(cd4_by_time_model_based$bmi_closest == 40.0 & round(cd4_by_time_model_based$t_lab,1) == 5.0))][1] return(list(year0_18.5=year0_18.5, year0_22=year0_22, year0_25=year0_25, year0_30=year0_30, year0_40=year0_40, year1_18.5=year1_18.5, year1_22=year1_22, year1_25=year1_25, year1_30=year1_30, year1_40=year1_40, year3_18.5=year3_18.5, year3_22=year3_22, year3_25=year3_25, year3_30=year3_30, year3_40=year3_40, year5_18.5=year5_18.5, year5_22=year5_22, year5_25=year5_25, year5_30=year5_30, year5_40=year5_40)) } # end of function graph_limits <- function(bmi){ # Overall bmi_rge_wm <- range(bmi$bmi_closest[which(bmi$sex.factor == 'male' & bmi$race.factor == 'white')],na.rm=TRUE) bmi_rge_wf <- range(bmi$bmi_closest[which(bmi$sex.factor == 'female' & bmi$race.factor == 'white')],na.rm=TRUE) bmi_rge_nwm <- range(bmi$bmi_closest[which(bmi$sex.factor == 'male' & bmi$race.factor == 'nonwhite')],na.rm=TRUE) bmi_rge_nwf <- range(bmi$bmi_closest[which(bmi$sex.factor == 'female' & bmi$race.factor == 'nonwhite')],na.rm=TRUE) bmi_rge_bm <- range(bmi$bmi_closest[which(bmi$sex.factor == 'male' & bmi$race_for_reviewers.factor == 'black')],na.rm=TRUE) bmi_rge_bf <- range(bmi$bmi_closest[which(bmi$sex.factor == 'female' & bmi$race_for_reviewers.factor == 'black')],na.rm=TRUE) bmi_rge_hm <- range(bmi$bmi_closest[which(bmi$sex.factor == 'male' & bmi$race_for_reviewers.factor == 'hispanic/other/unknown')],na.rm=TRUE) bmi_rge_hf <- range(bmi$bmi_closest[which(bmi$sex.factor == 'female' & bmi$race_for_reviewers.factor == 'hispanic/other/unknown')],na.rm=TRUE) # tmp used for lme tmp_rge_wm <- range(tmp$bmi_closest[which(tmp$sex.factor == 'male' & tmp$race.factor == 'white')],na.rm=TRUE) tmp_rge_wf <- range(tmp$bmi_closest[which(tmp$sex.factor == 'female' & tmp$race.factor == 'white')],na.rm=TRUE) tmp_rge_nwm <- range(tmp$bmi_closest[which(tmp$sex.factor == 'male' & tmp$race.factor == 'nonwhite')],na.rm=TRUE) tmp_rge_nwf <- range(tmp$bmi_closest[which(tmp$sex.factor == 'female' & tmp$race.factor == 'nonwhite')],na.rm=TRUE) tmp_rge_bm <- range(tmp_review$bmi_closest[which(tmp_review$sex.factor == 'male' & tmp_review$race_for_reviewers.factor == 'black')],na.rm=TRUE) tmp_rge_bf <- range(tmp_review$bmi_closest[which(tmp_review$sex.factor == 'female' & tmp_review$race_for_reviewers.factor == 'black')],na.rm=TRUE) tmp_rge_hm <- range(tmp_review$bmi_closest[which(tmp_review$sex.factor == 'male' & tmp_review$race_for_reviewers.factor == 'hispanic/other/unknown')],na.rm=TRUE) tmp_rge_hf <- range(tmp_review$bmi_closest[which(tmp_review$sex.factor == 'female' & tmp_review$race_for_reviewers.factor == 'hispanic/other/unknown')],na.rm=TRUE) # lme excluding records with extreme residuals tmp_w_residuals_rge_wm <- range(tmp$bmi_closest[which(tmp$sex.factor == 'male' & tmp$race.factor == 'white' & abs(tmp$residuals) <= 6)],na.rm=TRUE) tmp_w_residuals_rge_wf <- range(tmp$bmi_closest[which(tmp$sex.factor == 'female' & tmp$race.factor == 'white' & abs(tmp$residuals) <= 6)],na.rm=TRUE) tmp_w_residuals_rge_nwm <- range(tmp$bmi_closest[which(tmp$sex.factor == 'male' & tmp$race.factor == 'nonwhite' & abs(tmp$residuals) <= 6)],na.rm=TRUE) tmp_w_residuals_rge_nwf <- range(tmp$bmi_closest[which(tmp$sex.factor == 'female' & tmp$race.factor == 'nonwhite' & abs(tmp$residuals) <= 6)],na.rm=TRUE) tmp_w_residuals_rge_bm <- range(tmp_review$bmi_closest[which(tmp_review$sex.factor == 'male' & tmp_review$race_for_reviewers.factor == 'black' & abs(tmp_review$residuals) <= 6)],na.rm=TRUE) tmp_w_residuals_rge_bf <- range(tmp_review$bmi_closest[which(tmp_review$sex.factor == 'female' & tmp_review$race_for_reviewers.factor == 'black' & abs(tmp_review$residuals) <= 6)],na.rm=TRUE) tmp_w_residuals_rge_hm <- range(tmp_review$bmi_closest[which(tmp_review$sex.factor == 'male' & tmp_review$race_for_reviewers.factor == 'hispanic/other/unknown' & abs(tmp_review$residuals) <= 6)],na.rm=TRUE) tmp_w_residuals_rge_hf <- range(tmp_review$bmi_closest[which(tmp_review$sex.factor == 'female' & tmp_review$race_for_reviewers.factor == 'hispanic/other/unknown' & abs(tmp_review$residuals) <= 6)],na.rm=TRUE) overall_x_rge <- range(bmi$bmi_closest,na.rm=TRUE) return(list(bmi_rge_wm=bmi_rge_wm, bmi_rge_nwm = bmi_rge_nwm, bmi_rge_wf = bmi_rge_wf, bmi_rge_nwf = bmi_rge_nwf, bmi_rge_bm=bmi_rge_bm, bmi_rge_bf=bmi_rge_bf, bmi_rge_hm=bmi_rge_hm, bmi_rge_hf=bmi_rge_hf, tmp_rge_wm = tmp_rge_wm, tmp_rge_nwm = tmp_rge_nwm, tmp_rge_wf = tmp_rge_wf, tmp_rge_nwf = tmp_rge_nwf, tmp_rge_bm=tmp_rge_bm, tmp_rge_bf=tmp_rge_bf, tmp_rge_hm=tmp_rge_hm, tmp_rge_hf=tmp_rge_hf, tmp_w_residuals_rge_wm = tmp_w_residuals_rge_wm, tmp_w_residuals_rge_nwm = tmp_w_residuals_rge_nwm, tmp_w_residuals_rge_wf = tmp_w_residuals_rge_wf, tmp_w_residuals_rge_nwf = tmp_w_residuals_rge_nwf, tmp_w_residuals_rge_bm=tmp_w_residuals_rge_bm, tmp_w_residuals_rge_bf=tmp_w_residuals_rge_bf, tmp_w_residuals_rge_hm=tmp_w_residuals_rge_hm, tmp_w_residuals_rge_hf=tmp_w_residuals_rge_hf, overall_x_rge = overall_x_rge)) } # end of function data_mgmt_function <- function(data){ data$cohort <- factor(data$cohort) data <- upData(data, rename=c('newid'='cfar_pid', 'age_art'='age_fhaart', 'cd4_bs'='cd4_count_at_age_fhaart', 'vload_bs'='hiv1_rna_at_age_fhaart', 'cd4'='cd4_count', 'vload'='hiv1_rna', 'bmi_bs'='bmi_at_age_fhaart', 'ht_m_bs'='baseline_height_meters', 'wt_kg_bs'='baseline_weight_kgs', 'year_art'='year_of_first_haart_start', 'hepc_at_art'='has_hcv_diagnosis', 'pi_in_art'='pi_usage', 'idu'='evidence_of_idu', 'aids_at_art'='prior_ade')) # Creating age variable data$age <- data$age_fhaart + data$t_lab # Creating age_at_death variable data$age_at_death <- data$age_fhaart + data$t_dth # Creating age_at_last_visit (need to check with Bryan Lau whether t_lastact would correspond to this) data$age_at_last_visit <- data$age_fhaart + data$t_lastact # Modifying sex variable from 1,2 variable to male/female data$sex.factor <- factor(data$sex, levels=c(2,1), labels=c('female','male')) # Creating factor for has_hcv_diagnosis data$has_hcv_diagnosis.factor <- factor(data$has_hcv_diagnosis, levels=c(0,1), labels=c('No','Yes')) # Creating factor for art_naive data$art_naive.factor <- factor(data$art_naive, levels=c(0,1), labels=c('No','Yes')) # Creating factor for IDU data$evidence_of_idu.factor <- factor(data$evidence_of_idu, levels=c(0,1), labels=c('No','Yes')) # Creating factor variable for year_of_first_haart_start data$year_of_first_haart_start_categorized <- ifelse(data$year_of_first_haart_start %in% c(1998,1999,2000),0, ifelse(data$year_of_first_haart_start %in% c(2001,2002,2003),1, ifelse(data$year_of_first_haart_start %in% c(2004,2005,2006),2, ifelse(data$year_of_first_haart_start %in% c(2007,2008,2009,2010),3,NA)))) data$year_of_first_haart_start_categorized.factor <- factor(data$year_of_first_haart_start_categorized, levels=c(0,1,2,3), labels=c('1998-2000','2001-2003','2004-2006','2007-2010')) # Creating factor for PI usage data$pi_usage.factor <- factor(data$pi_usage, levels=c(0,1), labels=c('No','Yes')) # Creating factor for regimen class in ART data$regimen_class_in_art.factor <- factor(data$regimen_class_in_art, levels=c(1,2,3,4), labels=c('PI-based','NNRTI-based','NRTI-only','Other')) data <- data[order(data$cfar_pid, data$t_lab),] data$regimen_class_in_fhaart <- NA split(data$regimen_class_in_fhaart, data$cfar_pid) <- lapply(split(data$regimen_class_in_art, data$cfar_pid),FUN=function(y){ return(y[1]) }) data$regimen_class_in_fhaart.factor <- factor(data$regimen_class_in_fhaart, levels=c(1,2,3,4), labels=c('PI-based','NNRTI-based','NRTI-only','Other')) # Creating factor for prior ADE (prior to HAART start) data$prior_ade.factor <- factor(data$prior_ade, levels=c(0,1), labels=c('No','Yes')) # Put some kind of check in to make sure calculated correctly? data$bmi_check <- round(data$baseline_weight_kgs/(data$ht_m)^2,2) # Creating dichotomous race variable -- white vs nonwhite data$race.orig <- data$race n.race.orig <- as.numeric(table(data$race.orig)[which(names(table(data$race.orig))==5)]) n.race.orig <- ifelse(length(n.race.orig) == 0,0,n.race.orig) data$race.orig.factor <- factor(data$race.orig, levels=c(1,2,3,4,5), labels=c('White','Black','Hispanic','Other','Unknown')) data$race_for_reviewers <- ifelse(data$race %in% c(3,4,5),3,data$race) data$race_for_reviewers.factor <- factor(data$race_for_reviewers, levels=c(1,2,3), labels=c('white','black','hispanic/other/unknown')) data$race <- ifelse(!is.na(data$race) & data$race %in% c(2,3,4),1, ifelse(!is.na(data$race) & data$race == 1,0,NA)) data$race.factor <- factor(data$race, levels=c(0,1),labels=c('white','nonwhite')) # Treating level 5 - unknown/missing - as its own level data$race.3lvl <- ifelse(!is.na(data$race.orig) & data$race.orig %in% c(2,3,4),1, ifelse(data$race.orig == 5,2, ifelse(!is.na(data$race.orig),0,NA))) data$race.3lvl <- factor(data$race.3lvl, levels=c(0,1,2),labels=c('white','nonwhite','unknown/missing')) # Creating loghiv1_rna_at_age_fhaart data$loghiv1_rna_at_age_fhaart <- log10(data$hiv1_rna_at_age_fhaart) # Creating factor variables data$pi_usage.factor <- factor(data$pi_usage, levels=c(0,1), labels=c('No','Yes')) data$prior_ade.factor <- factor(data$prior_ade, levels=c(0,1), labels=c('No','Yes')) data$has_hcv_diagnosis.factor <- factor(data$has_hcv_diagnosis, levels=c(0,1), labels=c('No','Yes')) # Creating factor for death variable data$death.factor <- factor(data$death, levels=c(0,1), labels=c('Alive','Dead')) # There are duplicate records for id/date with wt differing but cd4 the same (2); but also have # duplicate records of id/date where wt differs data$pid.time <- paste(data$cfar_pid, data$t_lab,sep='.') data <- data[order(data$cfar_pid, data$t_lab),] data$bmi_new <- NA split(data$bmi_new, data$pid.time) <- lapply(split(data$bmi, data$pid.time),FUN=function(y){ if(length(y) > 1 & all(!is.na(y))){ mean(y,na.rm=TRUE) } else if(length(y) > 1 & all(is.na(y))) { NA } else { y }}) data$cd4_count_new <- NA split(data$cd4_count_new, data$pid.time) <- lapply(split(data$cd4_count, data$pid.time),FUN=function(y){ if(length(y) > 1 & all(!is.na(y))){ mean(y,na.rm=TRUE) } else if(length(y) > 1 & all(is.na(y))){ NA } else { y }}) data <- subset(data, !duplicated(pid.time)) data <- upData(data, drop=c('cd4_count','bmi')) data <- upData(data, rename=c('cd4_count_new'='cd4_count', 'bmi_new'='bmi')) bmi <- data #--------------------------------------------# # Categorizing BMI # #--------------------------------------------# bmi$year_of_first_haart_start.factor <- factor(bmi$year_of_first_haart_start, levels=c('1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009','2010')) bmi$bmi_at_age_fhaart_categorized_smaller <- ifelse(!is.na(bmi$bmi_at_age_fhaart) & bmi$bmi_at_age_fhaart < 25.0,0, ifelse(!is.na(bmi$bmi_at_age_fhaart) & bmi$bmi_at_age_fhaart >= 25.0 & bmi$bmi_at_age_fhaart < 30.0,1, ifelse(!is.na(bmi$bmi_at_age_fhaart) & bmi$bmi_at_age_fhaart >= 30.0,2,NA))) bmi$bmi_at_age_fhaart_categorized_smaller.factor <- factor(bmi$bmi_at_age_fhaart_categorized_smaller, levels=c(0,1,2), labels=c('< 25','[25-30)','>= 30')) bmi$bmi_at_age_fhaart_categorized_larger <- ifelse(!is.na(bmi$bmi_at_age_fhaart) & bmi$bmi_at_age_fhaart < 18.5,0, ifelse(!is.na(bmi$bmi_at_age_fhaart) & bmi$bmi_at_age_fhaart >= 18.5 & bmi$bmi_at_age_fhaart < 25.0,1, ifelse(!is.na(bmi$bmi_at_age_fhaart) & bmi$bmi_at_age_fhaart >= 25.0 & bmi$bmi_at_age_fhaart < 30.0,2, ifelse(!is.na(bmi$bmi_at_age_fhaart) & bmi$bmi_at_age_fhaart >= 30.0 & bmi$bmi_at_age_fhaart < 40.0,3, ifelse(!is.na(bmi$bmi_at_age_fhaart) & bmi$bmi_at_age_fhaart > 40.0,4,NA))))) bmi$bmi_at_age_fhaart_categorized_larger.factor <- factor(bmi$bmi_at_age_fhaart_categorized_larger, levels=c(0,1,2,3,4), labels=c('< 18.5','[18.5-25.0)','[25.0-30.0)','[30.0-40.0)','>= 40')) # Matching a BMI to a non-missing CD4 count that is within +/- 6 months of the CD4 count bmi$age_bmi <- ifelse(!is.na(bmi$bmi),bmi$age,NA) bmi$age_cd4 <- ifelse(!is.na(bmi$cd4_count),bmi$age,NA) bmi$bmi_closest <- NA bmi$age_bmi_closest <- NA for(i in 1:nrow(bmi)){ if(!is.na(bmi$cd4_count[i]) & !is.na(bmi$bmi[i])){ bmi$bmi_closest[i] <- bmi$bmi[i] bmi$age_bmi_closest[i] <- bmi$age[i] } else if(!is.na(bmi$cd4_count[i])){ id <- bmi$cfar_pid[i] current_age <- bmi$age_cd4[i] tmp <- subset(bmi, cfar_pid == id & !is.na(bmi) & abs(current_age - age_bmi) <= 0.5,select=c(cfar_pid,age,age_bmi,bmi)) if(nrow(tmp) > 0){ tmp$age_diff <- with(tmp, abs(current_age - age_bmi)) which.row <- which(tmp$age_diff == min(tmp$age_diff))[1] bmi$bmi_closest[i] <- tmp[which.row,'bmi'] bmi$age_bmi_closest[i] <- tmp$age[which.row] } } } # Need to subset bmi so that there are no non-missing CD4 counts (or missing BMI values) num_cd4_counts_total <- length(which(!is.na(bmi$cd4_count))) num_bmi_total <- length(which(!is.na(bmi$bmi))) pre_num_cd4_per_person <- tapply(bmi$cd4_count, bmi$cfar_pid, FUN=function(y){length(which(!is.na(y)))}) pre_num_bmi_per_person <- tapply(bmi$bmi, bmi$cfar_pid, FUN=function(y){length(which(!is.na(y)))}) bmi <- subset(bmi, !is.na(cd4_count) | !is.na(bmi_closest)) num_bmi_total2 <- length(unique(bmi$id_age_bmi_closest)) cd4_counts_per_person2 <- tapply(bmi$cd4_count, INDEX=bmi$cfar_pid, FUN=function(y){length(which(!is.na(y)))}) desc_cd4_counts_per_person <- summary(cd4_counts_per_person2) bmi_counts_per_person2 <- tapply(bmi$age_bmi_closest, INDEX=bmi$cfar_pid, FUN=function(y){length(unique(y))}) desc_bmi_per_person <- summary(bmi_counts_per_person2) # Will describe distribution of repeated bmi_closest values using age_bmi_closest bmi$id_age_bmi_closest <- paste(bmi$cfar_pid, round(bmi$age_bmi_closest,5), sep='.') nr <- tapply(bmi$id_age_bmi_closest, INDEX=bmi$cfar_pid, FUN=function(y){length(unique(y))}) summ_nr <- summary(nr) tk <- tapply(bmi$age_bmi_closest, INDEX=bmi$cfar_pid, FUN=function(y){length(which(duplicated(y)))/length(y)*100}) desc_tk <- summary(tk) bmi$cd4_count_sqrt <- sqrt(bmi$cd4_count) bmi$cd4_count_at_age_fhaart_sqrt <- as.numeric(sqrt(bmi$cd4_count_at_age_fhaart)) label(bmi$age_fhaart) <- 'Age at ART initiation' label(bmi$sex.factor) <- 'Sex' label(bmi$race.factor) <- 'Race' label(bmi$evidence_of_idu.factor) <- 'Evidence of IDU' label(bmi$year_of_first_haart_start) <- 'Year of ART initiation' label(bmi$year_of_first_haart_start_categorized.factor) <- 'Year of ART initiation (categorized)' label(bmi$art_naive.factor) <- 'ART naive' label(bmi$pi_usage.factor) <- 'PI in initial ART' label(bmi$regimen_class_in_fhaart.factor) <- 'Regimen class of initial ART' label(bmi$prior_ade.factor) <- 'Prior ADE' label(bmi$has_hcv_diagnosis.factor) <- 'HCV diagnosis' label(bmi$death.factor) <- 'Status' label(bmi$cd4_count_at_age_fhaart) <- 'CD4 count at ART initiation' label(bmi$loghiv1_rna_at_age_fhaart) <- 'HIV1-RNA at ART initiaion (log10-transformed)' label(bmi$bmi_at_age_fhaart) <- 'BMI at ART initiation' label(bmi$bmi_at_age_fhaart_categorized_larger.factor) <- 'BMI at ART initiation (categorized)' label(bmi$cd4_count_sqrt) <- 'CD4 count (square root transformed)' label(bmi$cd4_count_at_age_fhaart_sqrt) <- 'CD4 count at ART initiation (square root transformed)' label(bmi$race_for_reviewers.factor) <- 'Race (3-level)' return(list(bmi=bmi, n.race.orig=n.race.orig, num_cd4_counts_total=num_cd4_counts_total, num_bmi_total=num_bmi_total, pre_num_cd4_per_person=pre_num_cd4_per_person, pre_num_bmi_per_person=pre_num_bmi_per_person, num_bmi_total2=num_bmi_total2, cd4_counts_per_person2=cd4_counts_per_person2, desc_cd4_counts_per_person=desc_cd4_counts_per_person, bmi_counts_per_person2=bmi_counts_per_person2, desc_bmi_per_person=desc_bmi_per_person, nr=nr, summ_nr=summ_nr, tk=tk, desc_tk=desc_tk)) } # end of function @ <>= initial_data <- data_mgmt_function(data) bmi <- initial_data[['bmi']] n.race.orig <- initial_data[['n.race.orig']] num_cd4_counts_total <- initial_data[['num_cd4_counts_total']] num_bmi_total <- initial_data[['num_bmi_total']] pre_num_cd4_per_person <- initial_data[['pre_num_cd4_per_person']] pre_num_bmi_per_person <- initial_data[['pre_num_bmi_per_person']] num_bmi_total2 <- initial_data[['num_bmi_total2']] cd4_counts_per_person2 <- initial_data[['cd4_counts_per_person2']] desc_cd4_counts_per_person <- initial_data[['desc_cd4_counts_per_person']] bmi_counts_per_person2 <- initial_data[['bmi_counts_per_person2']] desc_bmi_per_person <- initial_data[['desc_bmi_per_person']] nr <- initial_data[['nr']] summ_nr <- initial_data[['summ_nr']] tk <- initial_data[['tk']] desc_tk <- initial_data[['desc_tk']] sens_data <- data_mgmt_function(data_total) bmi_sens <- sens_data[['bmi']] n.race.orig_sens <- sens_data[['n.race.orig']] num_cd4_counts_total_sens <- sens_data[['num_cd4_counts_total']] num_bmi_total_sens <- sens_data[['num_bmi_total']] pre_num_cd4_per_person_sens <- sens_data[['pre_num_cd4_per_person']] pre_num_bmi_per_person_sens <- sens_data[['pre_num_bmi_per_person']] num_bmi_total2_sens <- sens_data[['num_bmi_total2']] cd4_counts_per_person2_sens <- sens_data[['cd4_counts_per_person2']] desc_cd4_counts_per_person_sens <- sens_data[['desc_cd4_counts_per_person']] bmi_counts_per_person2_sens <- sens_data[['bmi_counts_per_person2']] desc_bmi_per_person_sens <- sens_data[['desc_bmi_per_person']] nr_sens <- sens_data[['nr']] summ_nr_sens <- sens_data[['summ_nr']] tk_sens <- sens_data[['tk']] desc_tk_sens <- sens_data[['desc_tk']] sens_data_female <- data_mgmt_function(data_total_female) bmi_sens_female <- sens_data_female[['bmi']] n.race.orig_sens_female <- sens_data_female[['n.race.orig']] num_cd4_counts_total_sens_female <- sens_data_female[['num_cd4_counts_total']] num_bmi_total_sens_female <- sens_data_female[['num_bmi_total']] pre_num_cd4_per_person_sens_female <- sens_data_female[['pre_num_cd4_per_person']] pre_num_bmi_per_person_sens_female <- sens_data_female[['pre_num_bmi_per_person']] num_bmi_total2_sens_female <- sens_data_female[['num_bmi_total2']] cd4_counts_per_person2_sens_female <- sens_data_female[['cd4_counts_per_person2']] desc_cd4_counts_per_person_sens_female <- sens_data_female[['desc_cd4_counts_per_person']] bmi_counts_per_person2_sens_female <- sens_data_female[['bmi_counts_per_person2']] desc_bmi_per_person_sens_female <- sens_data_female[['desc_bmi_per_person']] nr_sens_female <- sens_data_female[['nr']] summ_nr_sens_female <- sens_data_female[['summ_nr']] tk_sens_female <- sens_data_female[['tk']] desc_tk_sens_female <- sens_data_female[['desc_tk']] @ <>= # They want table 2 by men and women # Need 3 year BMI (in aim2_amendment_b.nw use the last BMI in data truncated at 3.5 years) tmp_3yr <- subset(bmi, t_lab <= 3.5, select=c('cfar_pid','bmi','t_lab')) tmp_3yr$three_yr_bmi <- NA split(tmp_3yr$three_yr_bmi, tmp_3yr$cfar_pid) <- lapply(split(tmp_3yr$bmi, tmp_3yr$cfar_pid),FUN=function(y){ lgth <- length(y) return(y[lgth]) }) bmi$three_yr_bmi <- tmp_3yr[match(bmi$cfar_pid, tmp_3yr$cfar_pid,nomatch=NA),'three_yr_bmi'] bmi$three_yr_bmi_categorized_larger <- ifelse(!is.na(bmi$three_yr_bmi) & bmi$three_yr_bmi < 18.5,0, ifelse(!is.na(bmi$three_yr_bmi) & bmi$three_yr_bmi >= 18.5 & bmi$three_yr_bmi < 25,1, ifelse(!is.na(bmi$three_yr_bmi) & bmi$three_yr_bmi >= 25 & bmi$three_yr_bmi < 30,2, ifelse(!is.na(bmi$three_yr_bmi) & bmi$three_yr_bmi >= 30 & bmi$three_yr_bmi < 40,3, ifelse(!is.na(bmi$three_yr_bmi) & bmi$three_yr_bmi >= 40,4,NA))))) bmi$three_yr_bmi_categorized_larger.factor <- factor(bmi$three_yr_bmi_categorized_larger, levels=c(0,1,2,3,4),labels=c('< 18.5','[18.5,25)','[25,30)','[30,40)','>= 40')) # Need 1 year BMI (in aim2_amendment_b.nw use the last BMI in data truncated at 1.5 years) tmp_1yr <- subset(bmi, t_lab <= 1.5, select=c('cfar_pid','bmi','t_lab')) tmp_1yr$one_yr_bmi <- NA split(tmp_1yr$one_yr_bmi, tmp_1yr$cfar_pid) <- lapply(split(tmp_1yr$bmi, tmp_1yr$cfar_pid),FUN=function(y){ lgth <- length(y) return(y[lgth]) }) bmi$one_yr_bmi <- tmp_1yr[match(bmi$cfar_pid, tmp_1yr$cfar_pid, nomatch=NA),'one_yr_bmi'] bmi$one_yr_bmi_categorized_larger <- ifelse(!is.na(bmi$one_yr_bmi) & bmi$one_yr_bmi < 18.5,0, ifelse(!is.na(bmi$one_yr_bmi) & bmi$one_yr_bmi >= 18.5 & bmi$one_yr_bmi < 25,1, ifelse(!is.na(bmi$one_yr_bmi) & bmi$one_yr_bmi >= 25 & bmi$one_yr_bmi < 30,2, ifelse(!is.na(bmi$one_yr_bmi) & bmi$one_yr_bmi >= 30 & bmi$one_yr_bmi < 40,3, ifelse(!is.na(bmi$one_yr_bmi) & bmi$one_yr_bmi >= 40,4,NA))))) bmi$one_yr_bmi_categorized_larger.factor <- factor(bmi$one_yr_bmi_categorized_larger, levels=c(0,1,2,3,4),labels=c('< 18.5','[18.5,25)','[25,30)','[30,40)','>= 40')) # Creating tables for men # 3 year versus baseline tbl_3_yr_m <- table(bmi$bmi_at_age_fhaart_categorized_larger.factor[which(bmi$sex.factor == 'male' & !duplicated(bmi$cfar_pid))], bmi$three_yr_bmi_categorized_larger.factor[which(bmi$sex.factor == 'male' & !duplicated(bmi$cfar_pid))]) tbl_3_yr_f <- table(bmi$bmi_at_age_fhaart_categorized_larger.factor[which(bmi$sex.factor == 'female' & !duplicated(bmi$cfar_pid))], bmi$three_yr_bmi_categorized_larger.factor[which(bmi$sex.factor == 'female' & !duplicated(bmi$cfar_pid))]) # 1 year versus baseline tbl_1_yr_m <- table(bmi$bmi_at_age_fhaart_categorized_larger.factor[which(bmi$sex.factor == 'male' & !duplicated(bmi$cfar_pid))], bmi$one_yr_bmi_categorized_larger.factor[which(bmi$sex.factor == 'male' & !duplicated(bmi$cfar_pid))]) tbl_1_yr_f <- table(bmi$bmi_at_age_fhaart_categorized_larger.factor[which(bmi$sex.factor == 'female' & !duplicated(bmi$cfar_pid))], bmi$one_yr_bmi_categorized_larger.factor[which(bmi$sex.factor == 'female' & !duplicated(bmi$cfar_pid))]) @ \section{Results} The data contain \Sexpr{nrow(bmi)} records representing \Sexpr{length(unique(bmi$cfar_pid))} subjects. There were \Sexpr{num_cd4_counts_total} non-missing CD4 counts in the data, and \Sexpr{num_bmi_total} non-missing BMI measurements in the data. For a given CD4 count, a BMI measurement was matched to it provided it fell within -180 to +180 days of the CD4 measurement. Therefore, in some instances, a single BMI measurement was matched to multiple CD4 measurements. Once a BMI measurement was matched to a given CD4 count, the distribution of the frequency of both measurements within a subject were calculated. Median (Interquartile Range [IQR]) CD4 counts per person were \Sexpr{desc_cd4_counts_per_person[3]} (\Sexpr{desc_cd4_counts_per_person[2]}, \Sexpr{desc_cd4_counts_per_person[5]}); there were \Sexpr{desc_bmi_per_person[3]} (\Sexpr{desc_bmi_per_person[2]}, \Sexpr{desc_bmi_per_person[5]}) median (IQR) BMI measurements per person. The median (IQR) percentage of records for a given person with BMI measurements matched to multiple CD4 counts is \Sexpr{desc_tk[3]} (\Sexpr{desc_tk[2]},\Sexpr{desc_tk[5]}). \vspace{0.1in} Given the diagnostics run on the subset of data used for coding, a linear mixed effect (lme) model with AR1 correlation was used as the primary model of interest. A generalized least squares (gls) model with AR1 correlation was also fit for completeness. \vspace{0.1in} Table 1 summarizes descriptive statistics on all relevant data by categorized BMI values and overall. Assuming the models with all interactions will be retained in the final analysis, Figures \ref{fig:cd4_by_bmi_lme_alt} and \ref{fig:cd4_by_bmi_lme_alt_wo_outliers} illustrate the relationship between time-varying CD4 and BMI using the lme model results stratified by sex and race. Similarly, Figures \ref{fig:cd4_by_bmi_gls} and \ref{fig:cd4_by_bmi_gls_alt} illustrate the relationship between time-varying CD4 and BMI using the gls model results stratified by sex and race. Figure \ref{fig:describe_ggplot} is a descriptive graph for CD4 over time by BMI categories based solely on the data; Figure \ref{fig:model_base} is a more model-based approach to Figure \ref{fig:describe_ggplot}. Analyses on each of the cohorts specified repeat the tables and graphs presented in the overall analysis. <>= tbl1 <- summary(bmi_at_age_fhaart_categorized_larger.factor ~ age_fhaart + sex.factor + race.factor + race_for_reviewers.factor + evidence_of_idu.factor + year_of_first_haart_start + year_of_first_haart_start_categorized.factor + pi_usage.factor + regimen_class_in_fhaart.factor + prior_ade.factor + has_hcv_diagnosis.factor + death.factor + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + bmi_at_age_fhaart + cohort, data=bmi, method='reverse',overall=TRUE,subset=!duplicated(cfar_pid),test=TRUE,lines.page=45) latex(tbl1, file='', long=TRUE,exclude1=FALSE,longtable=TRUE,caption='Descriptive statistics by categorized BMI at first HAART',digits=2,size='tiny',landscape=TRUE,prtest='P',lines.page=40) @ <>= tbl1_supp <- summary(year_of_first_haart_start_categorized.factor ~ sex.factor + race.factor + race_for_reviewers.factor + age_fhaart + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + pi_usage.factor + has_hcv_diagnosis.factor, data=bmi,subset=!duplicated(cfar_pid), method='reverse',overall=TRUE,test=TRUE) tbl1_supp$labels <- c('Sex','Race','Baseline age','Baseline CD4 count','Baseline viral load (log-transformed)','PI usage','Prior ADE','HCV diagnosis') latex(tbl1_supp, file='', exclude1=FALSE, long=TRUE, digits=2, prtest='P',landscape=TRUE,caption='Descriptive statistics on cohort by year of ART initiation') @ \clearpage <>= latex(tbl_1_yr_f, file='', caption='Baseline BMI versus BMI at 1 year in women.',rowname=rownames(tbl_1_yr_f),rgroup='Baseline BMI',n.rgroup=c(5),cgroup='One year BMI (women)',n.cgroup=c(5),where='!h',rowlabel='', insert.bottom=paste('Number of women with a non-missing BMI at 1 year = ',length(unique(bmi$cfar_pid[which(bmi$sex.factor == 'female' & !is.na(bmi$one_yr_bmi))])),'.',sep=''),size='footnotesize') latex(tbl_1_yr_m, file='', caption='Baseline BMI versus BMI at 1 year in men.',rowname=rownames(tbl_1_yr_m),rgroup='Baseline BMI',n.rgroup=c(5),cgroup='One year BMI (men)',n.cgroup=c(5),where='!h',rowlabel='', insert.bottom=paste('Number of men with a non-missing BMI at 1 year = ',length(unique(bmi$cfar_pid[which(bmi$sex.factor == 'male' & !is.na(bmi$one_yr_bmi))])),'.',sep=''),size='footnotesize') latex(tbl_3_yr_f, file='', caption='Baseline BMI versus BMI at 3 years in women.',rowname=rownames(tbl_3_yr_f),rgroup='Baseline BMI',n.rgroup=c(5),cgroup='Three year BMI (women)',n.cgroup=c(5),where='!h',rowlabel='', insert.bottom=paste('Number of women with a non-missing BMI at 3 years = ',length(unique(bmi$cfar_pid[which(bmi$sex.factor == 'female' & !is.na(bmi$three_yr_bmi))])),'.',sep=''),size='footnotesize') latex(tbl_3_yr_m, file='', caption='Baseline BMI versus BMI at 3 years in men.',rowname=rownames(tbl_3_yr_m),rgroup='Baseline BMI',n.rgroup=c(5),cgroup='Three year BMI (men)',n.cgroup=c(5),where='!h',rowlabel='', insert.bottom=paste('Number of men with a non-missing BMI at 3 years = ',length(unique(bmi$cfar_pid[which(bmi$sex.factor == 'male' & !is.na(bmi$three_yr_bmi))])),'.',sep=''),size='footnotesize') @ \clearpage <>= ddist <- datadist(bmi) options(datadist = 'ddist') #------------------------# # Formulas # #------------------------# fmla_wo_int <- as.formula('cd4_count_sqrt ~ rcs(bmi_closest,6) + rcs(t_lab,6) + sex.factor + race.factor + rcs(age_fhaart, 6) + regimen_class_in_fhaart.factor + rcs(year_of_first_haart_start,6) + rcs(cd4_count_at_age_fhaart,6) + rcs(loghiv1_rna_at_age_fhaart,6) + cohort') fmla_w_int <- as.formula('cd4_count_sqrt ~ rcs(bmi_closest,c(19.16,22.08,24.28,26.39,29.24,36.48))* (rcs(t_lab,c(0.099,0.74,1.66,2.88,4.70,8.24)) + sex.factor + race.factor) + sex.factor*race.factor + rcs(age_fhaart,c(24.87,34.13,38.44,42.37,46.87,54.73)) + regimen_class_in_fhaart.factor + rcs(year_of_first_haart_start,c(1998,2000,2001,2003,2004,2008)) + rcs(cd4_count_at_age_fhaart_sqrt,c(2.65,7.14,12.41,15.62,17.86,23.66)) + rcs(loghiv1_rna_at_age_fhaart,c(2.83,4.23,4.69,5.00,5.47,5.88)) + cohort') fmla_w_int_for_reviewers <- as.formula('cd4_count_sqrt ~ rcs(bmi_closest,c(19.16,22.08,24.28,26.39,29.24,36.48))* (rcs(t_lab,c(0.099,0.74,1.66,2.88,4.70,8.24)) + sex.factor + race_for_reviewers.factor) + sex.factor*race_for_reviewers.factor + rcs(age_fhaart,c(24.87,34.13,38.44,42.37,46.87,54.73)) + regimen_class_in_fhaart.factor + rcs(year_of_first_haart_start,c(1998,2000,2001,2003,2004,2008)) + rcs(cd4_count_at_age_fhaart_sqrt,c(2.65,7.14,12.41,15.62,17.86,23.66)) + rcs(loghiv1_rna_at_age_fhaart,c(2.83,4.23,4.69,5.00,5.47,5.88)) + cohort') fmla_w_int_wo_outliers <- as.formula('cd4_count_sqrt ~ rcs(bmi_closest,c(19.16,22.08,24.28,26.39,29.24,36.48))* (rcs(t_lab,c(0.099,0.74,1.66,2.88,4.70,8.24)) + sex.factor + race.factor) + sex.factor*race.factor + rcs(age_fhaart,c(24.87,34.13,38.44,42.37,46.87,54.73)) + regimen_class_in_fhaart.factor + rcs(year_of_first_haart_start,c(1998,2000,2001,2003,2004,2008)) + rcs(cd4_count_at_age_fhaart_sqrt,c(3,7.14,12.41,15.62,17.86,20)) + rcs(loghiv1_rna_at_age_fhaart,c(2.83,4.23,4.69,5.00,5.47,5.88)) + cohort') fmla_w_int_wo_outliers_for_reviewers <- as.formula('cd4_count_sqrt ~ rcs(bmi_closest,c(19.16,22.08,24.28,26.39,29.24,36.48))* (rcs(t_lab,c(0.099,0.74,1.66,2.88,4.70,8.24)) + sex.factor + race_for_reviewers.factor) + sex.factor*race_for_reviewers.factor + rcs(age_fhaart,c(24.87,34.13,38.44,42.37,46.87,54.73)) + regimen_class_in_fhaart.factor + rcs(year_of_first_haart_start,c(1998,2000,2001,2003,2004,2008)) + rcs(cd4_count_at_age_fhaart_sqrt,c(3,7.14,12.41,15.62,17.86,20)) + rcs(loghiv1_rna_at_age_fhaart,c(2.83,4.23,4.69,5.00,5.47,5.88)) + cohort') @ <>= #-----------------------------------------# # Generalized least squares # #-----------------------------------------# # With AR1 correlated data ddist <- datadist(bmi) options(datadist = 'ddist') gls_mod_wo_int_ar1 <- Gls(fmla_wo_int, data=bmi, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) gls_mod_w_int_ar1 <- Gls(fmla_w_int, data=bmi, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) # Sensitivity models ddist <- datadist(bmi_sens) options(datadist = 'ddist') gls_mod_wo_int_ar1_sens <- Gls(fmla_wo_int, data=bmi_sens, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) gls_mod_w_int_ar1_sens <- Gls(fmla_w_int, data=bmi_sens, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) #-----------------------------------------# # Mixed effects linear model # #-----------------------------------------# # AR1 correlation in lme ddist <- datadist(bmi) options(datadist = 'ddist') tmp <- subset(bmi, select=c('cfar_pid','cd4_count_sqrt','bmi_closest','t_lab','sex.factor','race.factor','age_fhaart','regimen_class_in_fhaart.factor', 'year_of_first_haart_start','cd4_count_at_age_fhaart_sqrt','loghiv1_rna_at_age_fhaart','cohort')) lme_mod_w_int_ar1 <- lme(fmla_w_int, random=~1|cfar_pid, data=tmp, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) # wo_high_residuals # The predicted values does not include those records with missing values so need to mark them to remove them tmp$missing <- apply(tmp, 1, FUN=function(y){ifelse(any(is.na(y)),1,0)}) tmp <- subset(tmp, missing == 0) tmp$residuals <- lme_mod_w_int_ar1[['residuals']][,'fixed'] lme_mod_w_int_ar1_wo_outliers <- lme(fmla_w_int_wo_outliers, random=~1|cfar_pid, data=tmp, subset=abs(residuals) <= 6, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) # Sensitivity models # AR1 correlation in lme ddist <- datadist(bmi_sens) options(datadist = 'ddist') tmp_sens <- subset(bmi_sens, select=c('cfar_pid','cd4_count_sqrt','bmi_closest','t_lab','sex.factor','race.factor','age_fhaart','regimen_class_in_fhaart.factor', 'year_of_first_haart_start','cd4_count_at_age_fhaart_sqrt','loghiv1_rna_at_age_fhaart','cohort')) lme_mod_w_int_ar1_sens <- lme(fmla_w_int, random=~1|cfar_pid, data=tmp_sens, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) # wo_high_residuals # The predicted values does not include those records with missing values so need to mark them to remove them tmp_sens$missing <- apply(tmp_sens, 1, FUN=function(y){ifelse(any(is.na(y)),1,0)}) tmp_sens <- subset(tmp_sens, missing == 0) tmp_sens$residuals <- lme_mod_w_int_ar1_sens[['residuals']][,'fixed'] lme_mod_w_int_ar1_wo_outliers_sens <- lme(fmla_w_int_wo_outliers, random=~1|cfar_pid, data=tmp_sens, subset=abs(residuals) <= 6, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) # -------------------------------------------- # # For reviewers # AR1 correlation in lme ddist <- datadist(bmi) options(datadist = 'ddist') tmp_review <- subset(bmi, select=c('cfar_pid','cd4_count_sqrt','bmi_closest','t_lab','sex.factor','race_for_reviewers.factor','age_fhaart','regimen_class_in_fhaart.factor', 'year_of_first_haart_start','cd4_count_at_age_fhaart_sqrt','loghiv1_rna_at_age_fhaart','cohort')) lme_mod_w_int_ar1_review <- lme(fmla_w_int_for_reviewers, random=~1|cfar_pid, data=tmp_review, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) tmp_review$missing <- apply(tmp_review, 1, FUN=function(y){ifelse(any(is.na(y)),1,0)}) tmp_review <- subset(tmp_review, missing == 0) tmp_review$residuals <- lme_mod_w_int_ar1_review[['residuals']][,'fixed'] lme_mod_w_int_ar1_wo_outliers_review <- lme(fmla_w_int_wo_outliers_for_reviewers, random=~1|cfar_pid, data=tmp_review, subset=abs(residuals) <= 6, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) @ <>= # In primary models main_mod_assessment <- model_fit_assessment(gls_mod=gls_mod_w_int_ar1, lme_mod=lme_mod_w_int_ar1, data=bmi, fmla=fmla_w_int) cor_df <- main_mod_assessment[['cor_df']] dt <- main_mod_assessment[['dt']] dr <- main_mod_assessment[['dr']] dt_wo_outliers <- main_mod_assessment[['dt_wo_outliers']] dr_wo_outliers <- main_mod_assessment[['dr_wo_outliers']] lm_residuals_wo_outliers <- main_mod_assessment[['lm_residuals_wo_outliers']] # In sensitivity main_mod_assessment_sens <- model_fit_assessment(gls_mod=gls_mod_w_int_ar1_sens, lme_mod=lme_mod_w_int_ar1_sens, data=bmi_sens, fmla=fmla_w_int) cor_df_sens <- main_mod_assessment_sens[['cor_df']] dt_sens <- main_mod_assessment_sens[['dt']] dr_sens <- main_mod_assessment_sens[['dr']] dt_wo_outliers_sens <- main_mod_assessment_sens[['dt_wo_outliers']] dr_wo_outliers_sens <- main_mod_assessment_sens[['dr_wo_outliers']] lm_residuals_wo_outliers_sens <- main_mod_assessment_sens[['lm_residuals_wo_outliers']] @ \subsection{Assessing appropriate correlation structure} <>= # Table to evaluate the minimum AIC for each of the models. cor_df @ \clearpage %\begin{figure} %\caption{Variogram for assessing correlation structure} %<>= %#pdf('Variogram_wo_points_and_w_lowess_line_capped.pdf') %# plot(dt_wo_outliers, dr_wo_outliers,main='',xlab='dt',ylab='dr',type='n',ylim=c(0,30)) %# lines(lowess(x=dt_wo_outliers,y=dr_wo_outliers),col='red') %# abline(h=var(lm_residuals_wo_outliers),col='blue') %#dev.off() %@ %\end{figure} %\begin{figure} %\caption{Autocorrelation function graphs for the lme and gls models} %<>= %# # Autocorrelation function graphs %# acf1 <- plot(ACF(lme_mod_w_int_ar1, maxLag=10,resType='n'),alpha=0.01,main='lme model',cex.main=0.65) %# acf2 <- plot(ACF(gls_mod_w_int_ar1,resType="normalized"),alpha=0.05,main='gls model',cex.main=0.65) %# %# print(acf1, position=c(0,0.5,1,1.0),more=TRUE) %# print(acf2, position=c(0,0,1,0.45)) %@ %\end{figure} %\begin{figure} %\caption{QQ plots for all data and for data in which residuals whose absolute value exceeds 6 are removed} %<>= %# par(mfrow=c(2,2),mar=c(1,1,1,1),oma=c(1,1,1,1)) %# # QQ plot %# #pdf('test_qq.pdf') %# qq1 <- qqnorm(lme_mod_w_int_ar1, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1),cex.lab=0.65) %# #dev.off() %# %# qq2 <- qqnorm(lme_mod_w_int_ar1_wo_outliers, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1),xlab='Normalized residuals (outliers removed)',cex.lab=0.65) %# %# print(qq1, position=c(0,0.5,1,1.0),more=TRUE) %# print(qq2, position=c(0,0,1,0.45)) %@ %\end{figure} \clearpage \subsection{Comparing model output} \tiny <>= # ---------------------------------------- # # Comparing the different model outputs: # # ---------------------------------------- # # Gls model using AR1 correlation form=~1|ID #summary(gls_mod_w_int_ar1) gls_coeff <- gls_mod_w_int_ar1$coefficients # lme model using AR1 correlation form=~1|ID #summary(lme_mod_w_int_ar1) lme_coeff <- lme_mod_w_int_ar1$coefficients$fixed cbind(format(round(lme_coeff,5),scientific=FALSE),format(round(gls_coeff,5),scientific=FALSE)) @ \clearpage \subsection{Evaluating models with interaction terms} <>= # ----------------------------- # # Evaluating significance of # # interaction terms # # ----------------------------- # # ANOVA table for gls models anova(gls_mod_w_int_ar1, tol=1e-13) # ANOVA table for lme models anova(lme_mod_w_int_ar1) @ <>= # Bryan wants the individual lines on the graphs to only span the BMIs found for that race/sex group main_mod_limits <- graph_limits(bmi) bmi_rge_wm <- main_mod_limits[['bmi_rge_wm']] bmi_rge_nwm <- main_mod_limits[['bmi_rge_nwm']] bmi_rge_wf <- main_mod_limits[['bmi_rge_wf']] bmi_rge_nwf <- main_mod_limits[['bmi_rge_nwf']] bmi_rge_bm <- main_mod_limits[['bmi_rge_bm']] bmi_rge_bf <- main_mod_limits[['bmi_rge_bf']] bmi_rge_hm <- main_mod_limits[['bmi_rge_hm']] bmi_rge_hf <- main_mod_limits[['bmi_rge_hf']] tmp_rge_wm <- main_mod_limits[['tmp_rge_wm']] tmp_rge_nwm <- main_mod_limits[['tmp_rge_nwm']] tmp_rge_wf <- main_mod_limits[['tmp_rge_wf']] tmp_rge_nwf <- main_mod_limits[['tmp_rge_nwf']] tmp_rge_bm <- main_mod_limits[['tmp_rge_bm']] tmp_rge_bf <- main_mod_limits[['tmp_rge_bf']] tmp_rge_hm <- main_mod_limits[['tmp_rge_hm']] tmp_rge_hf <- main_mod_limits[['tmp_rge_hf']] tmp_w_residuals_rge_wm <- main_mod_limits[['tmp_w_residuals_rge_wm']] tmp_w_residuals_rge_nwm <- main_mod_limits[['tmp_w_residuals_rge_nwm']] tmp_w_residuals_rge_wf <- main_mod_limits[['tmp_w_residuals_rge_wf']] tmp_w_residuals_rge_nwf <- main_mod_limits[['tmp_w_residuals_rge_nwf']] tmp_w_residuals_rge_bm <- main_mod_limits[['tmp_w_residuals_rge_bm']] tmp_w_residuals_rge_bf <- main_mod_limits[['tmp_w_residuals_rge_bf']] tmp_w_residuals_rge_hm <- main_mod_limits[['tmp_w_residuals_rge_hm']] tmp_w_residuals_rge_hf <- main_mod_limits[['tmp_w_residuals_rge_hf']] overall_x_rge <- main_mod_limits[['overall_x_rge']] @ <>= # From Gls gls_predict <- Predict(gls_mod_w_int_ar1,bmi_closest,sex.factor,race.factor) # Transforming CD4 back from square root transformed gls_predict$yhat <- gls_predict$yhat^2 gls_predict$lower <- gls_predict$lower^2 gls_predict$upper <- gls_predict$upper^2 rge <- range(c(gls_predict$lower,gls_predict$upper)) x.rge <- range(c(gls_predict$bmi_closest)) x.rge <- c(x.rge[1],quantile(gls_predict$bmi_closest,prob=0.90)) x.rge_gls_alt <- x.rge rge_gls_alt <- rge # ------------------------------ # # From lme # Getting predicted CD4 count values (square root transformed) cohort_mode_tbl <- sort(table(bmi$cohort),decreasing=TRUE) cohort_mode <- names(cohort_mode_tbl)[1] lme_predict_data <- data.frame('bmi_closest'=rep(seq(min(bmi$bmi_closest,na.rm=TRUE),max(bmi$bmi_closest,na.rm=TRUE),length.out=100),times=4), 't_lab'=rep(median(bmi$t_lab),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100),rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(bmi$age_fhaart[!duplicated(bmi$cfar_pid)]),times=400), # 'regimen_class_in_fhaart.factor'=rep(levels(bmi$regimen_class_in_fhaart.factor)[1],times=400), 'regimen_class_in_fhaart.factor'=rep('NNRTI-based',times=400), 'year_of_first_haart_start'=rep(median(bmi$year_of_first_haart_start),times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(bmi$cd4_count_at_age_fhaart_sqrt[!duplicated(bmi$cfar_pid)]),times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(bmi$loghiv1_rna_at_age_fhaart[!duplicated(bmi$cfar_pid)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode,times=400)) # 'cohort'=rep(levels(bmi$cohort)[1],times=400)) lme_predict <- predict(lme_mod_w_int_ar1, newdata=lme_predict_data,na.action=na.omit,level=0) ## Transforming CD4 back from square root transformed lme_predict <- lme_predict^2 # rge_lme <- range(lme_predict) # x.rge_lme <- range(bmi$bmi_closest,na.rm=TRUE) # x.rge_lme <- c(x.rge_lme[1],quantile(bmi$bmi_closest,prob=0.90,na.rm=TRUE)) @ % \begin{figure} % \caption{Predicted CD4 count by BMI using lme}\label{fig:cd4_by_bmi_lme_alt} % <>= % # # Getting predicted CD4 count values (square root transformed) % # lme_predict_data <- data.frame('bmi_closest'=rep(seq(min(bmi$bmi_closest,na.rm=TRUE),max(bmi$bmi_closest,na.rm=TRUE),length.out=100),times=4), % # 't_lab'=rep(median(bmi$t_lab),times=400), % # 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), % # 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100),rep('white',times=100),rep('nonwhite',times=100)), % # 'age_fhaart'=rep(median(bmi$age_fhaart[!duplicated(bmi$cfar_pid)]),times=400), % # 'regimen_class_in_art.factor'=rep(levels(bmi$regimen_class_in_art.factor)[1],times=400), % # 'year_of_first_haart_start'=rep(median(bmi$year_of_first_haart_start),times=400), % # 'cd4_count_at_age_fhaart_sqrt'=rep(median(bmi$cd4_count_at_age_fhaart_sqrt[!duplicated(bmi$cfar_pid)]),times=400), % # 'loghiv1_rna_at_age_fhaart'=rep(median(bmi$loghiv1_rna_at_age_fhaart[!duplicated(bmi$cfar_pid)],na.rm=TRUE),times=400), % # 'cohort'=rep(levels(bmi$cohort)[1],times=400)) % # % # cjb_lme <- predict(lme_mod_w_int_ar1, newdata=lme_predict_data,na.action=na.omit,level=0) % # % # % # ## Transforming CD4 back from square root transformed % # cjb_lme <- cjb_lme^2 % # % # rge_lme <- range(cjb_lme) % # x.rge_lme <- range(bmi$bmi_closest,na.rm=TRUE) % # x.rge_lme <- c(x.rge_lme[1],quantile(bmi$bmi_closest,prob=0.90,na.rm=TRUE)) % % # White males % plot(1:400,1:400,type='n',xlab='BMI',ylab='CD4 count',xlim=overall_x_rge,ylim=rge) % lines(cjb_lme[which(lme_predict_data$sex.factor == 'male' & lme_predict_data$race.factor == 'white' & lme_predict_data$bmi_closest >= bmi_rge_wm[1] & lme_predict_data$bmi_closest <= bmi_rge_wm[2])] ~ % lme_predict_data$bmi_closest[which(lme_predict_data$sex.factor == 'male' & lme_predict_data$race.factor == 'white' & lme_predict_data$bmi_closest >= bmi_rge_wm[1] & lme_predict_data$bmi_closest <= bmi_rge_wm[2])], % ylim=rge,ylab='CD4 count',xlab='BMI',xlim=bmi_rge_wm) % % # Non-white males % lines(cjb_lme[which(lme_predict_data$sex.factor == 'male' & lme_predict_data$race.factor == 'nonwhite' & lme_predict_data$bmi_closest >= bmi_rge_nwm[1] & lme_predict_data$bmi_closest <= bmi_rge_nwm[2])] ~ % lme_predict_data$bmi_closest[which(lme_predict_data$sex.factor == 'male' & lme_predict_data$race.factor == 'nonwhite' & lme_predict_data$bmi_closest >= bmi_rge_nwm[1] & lme_predict_data$bmi_closest <= bmi_rge_nwm[2])], % col='red') % % # White females % lines(cjb_lme[which(lme_predict_data$sex.factor == 'female' & lme_predict_data$race.factor == 'white' & lme_predict_data$bmi_closest >= bmi_rge_wf[1] & lme_predict_data$bmi_closest <= bmi_rge_wf[2])] ~ % lme_predict_data$bmi_closest[which(lme_predict_data$sex.factor == 'female' & lme_predict_data$race.factor == 'white' & lme_predict_data$bmi_closest >= bmi_rge_wf[1] & lme_predict_data$bmi_closest <= bmi_rge_wf[2])], % col='blue') % % # Non-white females % lines(cjb_lme[which(lme_predict_data$sex.factor == 'female' & lme_predict_data$race.factor == 'nonwhite' & lme_predict_data$bmi_closest >= bmi_rge_nwf[1] & lme_predict_data$bmi_closest <= bmi_rge_nwf[2])] ~ % lme_predict_data$bmi_closest[which(lme_predict_data$sex.factor == 'female' & lme_predict_data$race.factor == 'nonwhite' & lme_predict_data$bmi_closest >= bmi_rge_nwf[1] & lme_predict_data$bmi_closest <= bmi_rge_nwf[2])], % col='green') % % legend('topleft',fill=c('black','red','blue','green'),legend=c('White male','Non-white male','White female','Non-white female')) % @ % \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI using lme (ggplot version)}\label{fig:cd4_by_bmi_lme_alt} <>= cd4_by_bmi_figure(lme_predict,lme_predict_data, main_mod_limits, plot_method='ggplot2', subset_method='overall') @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI using lme and omitting those with abs(residuals) $>$ 6 (base graphics version)}\label{fig:cd4_by_bmi_lme_alt_wo_outliers} <>= cohort_mode_tbl <- sort(table(bmi$cohort),decreasing=TRUE) cohort_mode <- names(cohort_mode_tbl)[1] lme_predict_data_wo_outliers <- data.frame('bmi_closest'=rep(seq(min(tmp$bmi_closest[which(abs(tmp$residuals) <= 6)],na.rm=TRUE),max(tmp$bmi_closest[which(abs(tmp$residuals) <= 6)],na.rm=TRUE),length.out=100),times=4), 't_lab'=rep(median(tmp$t_lab[which(abs(tmp$residuals) <= 6)]),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100),rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'= rep(median(tmp$age_fhaart[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)]),times=400), # 'regimen_class_in_fhaart.factor'=rep(sort(unique(tmp$regimen_class_in_fhaart.factor[which(abs(tmp$residuals) <= 6)]))[1],times=400), 'regimen_class_in_fhaart.factor'=rep(sort(unique(tmp$regimen_class_in_fhaart.factor[which(abs(tmp$residuals) <= 6)]))[2],times=400), 'year_of_first_haart_start'=rep(median(tmp$year_of_first_haart_start[which(abs(tmp$residuals) <= 6)]),times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp$cd4_count_at_age_fhaart_sqrt[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)]),times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp$loghiv1_rna_at_age_fhaart[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode,times=400)) # 'cohort'=rep(sort(unique(tmp$cohort[which(abs(tmp$residuals) <= 6)]))[1],times=400)) lme_predict_wo_outliers <- predict(lme_mod_w_int_ar1_wo_outliers, newdata=lme_predict_data_wo_outliers,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed lme_predict_wo_outliers <- lme_predict_wo_outliers^2 cd4_by_bmi_figure(mod_predict=lme_predict_wo_outliers, predict_data=lme_predict_data_wo_outliers, fig_limits=main_mod_limits, plot_method='base', subset_method='wo_outliers') @ \end{figure} \clearpage \begin{figure} \caption{Predicted CD4 count by BMI in men using lme and omitting those with abs(residuals) $>$ 6 (base graphics version) and using 3-level race requested by reviewers.}\label{fig:cd4_by_bmi_lme_alt_wo_outliers_review_m} <>= cohort_mode_tbl_review <- sort(table(bmi$cohort),decreasing=TRUE) cohort_mode_review <- names(cohort_mode_tbl_review)[1] lme_predict_data_wo_outliers_review <- data.frame('bmi_closest'=rep(seq(min(tmp_review$bmi_closest[which(abs(tmp_review$residuals) <= 6)],na.rm=TRUE), max(tmp_review$bmi_closest[which(abs(tmp_review$residuals) <= 6)],na.rm=TRUE),length.out=80),times=6), 't_lab'=rep(median(tmp_review$t_lab[which(abs(tmp_review$residuals) <= 6)]),times=480), 'sex.factor'=c(rep('female',times=240),rep('male',times=240)), 'race_for_reviewers.factor'=c(rep('white',times=80),rep('black',times=80),rep('hispanic/other/unknown',times=80), rep('white',times=80),rep('black',times=80),rep('hispanic/other/unknown',times=80)), 'age_fhaart'= rep(median(tmp_review$age_fhaart[which(!duplicated(tmp_review$cfar_pid) & abs(tmp_review$residuals) <= 6)]),times=480), 'regimen_class_in_fhaart.factor'=rep(sort(unique(tmp_review$regimen_class_in_fhaart.factor[which(abs(tmp_review$residuals) <= 6)]))[2],times=480), 'year_of_first_haart_start'=rep(median(tmp_review$year_of_first_haart_start[which(abs(tmp_review$residuals) <= 6)]),times=480), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp_review$cd4_count_at_age_fhaart_sqrt[which(!duplicated(tmp_review$cfar_pid) & abs(tmp_review$residuals) <= 6)]),times=480), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp_review$loghiv1_rna_at_age_fhaart[which(!duplicated(tmp_review$cfar_pid) & abs(tmp_review$residuals) <= 6)],na.rm=TRUE),times=480), 'cohort'=rep(cohort_mode_review,times=480)) lme_predict_wo_outliers_review <- predict(lme_mod_w_int_ar1_wo_outliers_review, newdata=lme_predict_data_wo_outliers_review,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed lme_predict_wo_outliers_review <- lme_predict_wo_outliers_review^2 cd4_by_bmi_figure_for_reviewers <- function(mod_predict,predict_data, fig_limits, plot_method, subset_method,sex){ rge_wm <- fig_limits[['tmp_w_residuals_rge_wm']] rge_bm <- fig_limits[['tmp_w_residuals_rge_bm']] rge_hm <- fig_limits[['tmp_w_residuals_rge_hm']] rge_wf <- fig_limits[['tmp_w_residuals_rge_wf']] rge_bf <- fig_limits[['tmp_w_residuals_rge_bf']] rge_hf <- fig_limits[['tmp_w_residuals_rge_hf']] overall_x_rge <- fig_limits[['overall_x_rge']] # White plot(mod_predict[which(predict_data$sex.factor == sex & predict_data$race_for_reviewers.factor == 'white' & predict_data$bmi_closest >= rge_wm[1] & predict_data$bmi_closest <= min(c(rge_wm[2],50)))] ~ predict_data$bmi_closest[which(predict_data$sex.factor == sex & predict_data$race_for_reviewers.factor == 'white' & predict_data$bmi_closest >= rge_wm[1] & predict_data$bmi_closest <= min(c(rge_wm[2],50)))], type='l',ylim=range(mod_predict),ylab='CD4 count',xlab='BMI',xlim=c(10,50),main=capitalize(sex)) # Black lines(mod_predict[which(predict_data$sex.factor == sex & predict_data$race_for_reviewers.factor == 'black' & predict_data$bmi_closest >= rge_bm[1] & predict_data$bmi_closest <= min(c(rge_bm[2],50)))] ~ predict_data$bmi_closest[which(predict_data$sex.factor == sex & predict_data$race_for_reviewers.factor == 'black' & predict_data$bmi_closest >= rge_bm[1] & predict_data$bmi_closest <= min(c(rge_bm[2],50)))],col='red') # Hispanic/Other/Unknown lines(mod_predict[which(predict_data$sex.factor == sex & predict_data$race_for_reviewers.factor == 'hispanic/other/unknown' & predict_data$bmi_closest >= rge_hm[1] & predict_data$bmi_closest <= min(c(rge_hm[2],50)))] ~ predict_data$bmi_closest[which(predict_data$sex.factor == sex & predict_data$race_for_reviewers.factor == 'hispanic/other/unknown' & predict_data$bmi_closest >= rge_hm[1] & predict_data$bmi_closest <= min(c(rge_hm[2],50)))],col='green') legend('topleft',fill=c('black','red','green'),legend=c('White','Black','Hispanic/other/unknown'),cex=0.90,bty='n') }# end of function cd4_by_bmi_figure_for_reviewers(mod_predict=lme_predict_wo_outliers_review, predict_data=lme_predict_data_wo_outliers_review, fig_limits=main_mod_limits,sex='male') @ \end{figure} \clearpage \begin{figure} \caption{Predicted CD4 count by BMI in women using lme and omitting those with abs(residuals) $>$ 6 (base graphics version) and using 3-level race requested by reviewers.}\label{fig:cd4_by_bmi_lme_alt_wo_outliers_review_f} <>= cd4_by_bmi_figure_for_reviewers(mod_predict=lme_predict_wo_outliers_review, predict_data=lme_predict_data_wo_outliers_review, fig_limits=main_mod_limits,sex='female') @ \end{figure} \clearpage \begin{figure} \caption{Predicted CD4 count by BMI from gls. Predicted values (and their 95\% CIs) for the designated race and sex combination are plotted in purple; the black lines represent the predicted values for the other sex/race combinations for reference.}\label{fig:cd4_by_bmi_gls} <>= ##pdf('test_text.pdf') par(mfrow=c(2,2),oma=c(1,3,1,1),mar=c(4,5,1,1),cex.axis=0.90) #rge <- range(c(cjb$lower,cjb$upper)) #x.rge <- range(c(cjb$bmi_closest)) #x.rge <- c(x.rge[1],quantile(cjb$bmi_closest,prob=0.90)) rge <- range(c(gls_predict$lower,gls_predict$upper)) x.rge <- range(c(gls_predict$bmi_closest)) x.rge <- c(x.rge[1],quantile(gls_predict$bmi_closest,prob=0.90)) x.rge[2] <- min(c(x.rge[2],50)) # White males plot(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')], type='l',ylim=rge,xlab='',ylab='CD4 count',main='White',xlim=x.rge) polygon(c(gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')], rev(gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')])), c(gls_predict$lower[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')],rev(gls_predict$upper[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')])), col=rgb(0,0,0.255,0.35),border=FALSE) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')],col='blue') lines(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')],col='black',lty=2) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')],col='black',lty=3) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')],col='black',lty=4) # text(x=-1.75,y=(rge[1]+(rge[2]-rge[1])/2),'Male',xpd=NA,srt=90,cex=1.1,font=2) # text(x=-(x.rge[1] + (x.rge[2]-x.rge[1])/2),y=(rge[1]+(rge[2]-rge[1])/2),'Male',xpd=NA,srt=90,cex=1.1,font=2) text(x=-10,y=(rge[1]+(rge[2]-rge[1])/2),'Male',xpd=NA,srt=90,cex=1.1,font=2) # mtext(side=3,outer=TRUE,text='CD4 count',line=-1,font=2) # Non-white males plot(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')], type='l',ylim=rge,xlab='',ylab='',main='Non-white',xlim=x.rge) polygon(c(gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')], rev(gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')])), c(gls_predict$lower[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')],rev(gls_predict$upper[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')])), col=rgb(0,0,0.255,0.35),border=FALSE) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')],col='blue') lines(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')],col='black',lty=2) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')],col='black',lty=3) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')],col='black',lty=4) # White females plot(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')], type='l',ylim=rge,xlab='BMI', ylab='CD4 count',main='',xlim=x.rge) polygon(c(gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')], rev(gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')])), c(gls_predict$lower[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')],rev(gls_predict$upper[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')])), col=rgb(0.255,0.20,0.147,0.25),border=FALSE) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')],col=6) # text(x=-1.75,y=(rge[1]+(rge[2]-rge[1])/2),'Female',xpd=NA,srt=90,cex=1.1,font=2) # text(x=-(x.rge[1] + (x.rge[2]-x.rge[1])/2),y=(rge[1]+(rge[2]-rge[1])/2),'Female',xpd=NA,srt=90,cex=1.1,font=2) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')],col='black',lty=2) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')],col='black',lty=3) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')],col='black',lty=4) text(x=-10,y=(rge[1]+(rge[2]-rge[1])/2),'Female',xpd=NA,srt=90,cex=1.1,font=2) # Non-white females plot(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')], type='l',ylim=rge,xlab='BMI',ylab='',main='',xlim=x.rge) polygon(c(gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')], rev(gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')])), c(gls_predict$lower[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')],rev(gls_predict$upper[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')])), col=rgb(0.255,0.20,0.147,0.25),border=FALSE) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')],col=6) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')],col='black',lty=2) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite')],col='black',lty=3) lines(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white')] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite')],col='black',lty=4) @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI using gls}\label{fig:cd4_by_bmi_gls_alt} <>= # White males plot(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white' & gls_predict$bmi_closest >= bmi_rge_wm[1] & gls_predict$bmi_closest <= min(c(bmi_rge_wm[2],50)))] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'white' & gls_predict$bmi_closest >= bmi_rge_wm[1] & gls_predict$bmi_closest <= min(c(bmi_rge_wm[2],50)))], type='l',ylim=rge,ylab='CD4 count',xlab='BMI',xlim=c(overall_x_rge[1],min(c(overall_x_rge[2],50)))) # Non-white males lines(gls_predict$yhat[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite' & gls_predict$bmi_closest >= bmi_rge_nwm[1] & gls_predict$bmi_closest <= min(c(bmi_rge_nwm[2],50)))] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'male' & gls_predict$race.factor == 'nonwhite' & gls_predict$bmi_closest >= bmi_rge_nwm[1] & gls_predict$bmi_closest <= min(c(bmi_rge_nwm[2],50)))], col='red') # White females lines(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white' & gls_predict$bmi_closest >= bmi_rge_wf[1] & gls_predict$bmi_closest <= min(c(bmi_rge_wf[2],50)))] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'white' & gls_predict$bmi_closest >= bmi_rge_wf[1] & gls_predict$bmi_closest <= min(c(bmi_rge_wf[2],50)))], col='blue') # Non-white females lines(gls_predict$yhat[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite' & gls_predict$bmi_closest >= bmi_rge_nwf[1] & gls_predict$bmi_closest <= min(c(bmi_rge_nwf[2],50)))] ~ gls_predict$bmi_closest[which(gls_predict$sex.factor == 'female' & gls_predict$race.factor == 'nonwhite' & gls_predict$bmi_closest >= bmi_rge_nwf[1] & gls_predict$bmi_closest <= min(c(bmi_rge_nwf[2],50)))], col='green') legend('topleft',fill=c('black','red','blue','green'),legend=c('White male','Non-white male','White female','Non-white female')) @ \end{figure} <>= cd4_by_year_combinednew1 <- descriptive_graph_mgmt(bmi) @ \begin{figure} \caption{Descriptive graph of CD4 count by year since ART initiation and BMI category. Yearly CD4 counts and BMI measurements were defined as the closest CD4 count within +/- 180 days. A subject could switch BMI categories from year to year. Points at a given year and BMI category are the average of all CD4 counts within that year and BMI category.}\label{fig:describe_ggplot} <>= ggplot(cd4_by_year_combinednew1, aes(x=Time,y=CD4_count, group=BMI, colour=BMI)) + geom_line() + xlab('Time since ART initiation (years)') + ylab('CD4 count') + scale_colour_manual(values=c('black','red','green','blue','yellow','purple')) @ \end{figure} \begin{figure} \caption{Model-based graph of CD4 counts by year since ART initiation and BMI cut points corresponding to common BMI categories. Predicted CD4 values were calculated over 4 years for each of the cut points listed in the legend.}\label{fig:model_base} <>= cd4_by_time_model_based_figure(data=tmp, mod=lme_mod_w_int_ar1_wo_outliers, method='base') @ \end{figure} \clearpage <>= all_comers_predicted_cd4 <- predicted_cd4_function(data=tmp, mod=lme_mod_w_int_ar1_wo_outliers) predicted_df <- data.frame('Year'=c(0,1,3,5), '18.5'=round(c(all_comers_predicted_cd4[['year0_18.5']], all_comers_predicted_cd4[['year1_18.5']], all_comers_predicted_cd4[['year3_18.5']], all_comers_predicted_cd4[['year5_18.5']])), '22'=round(c(all_comers_predicted_cd4[['year0_22']], all_comers_predicted_cd4[['year1_22']], all_comers_predicted_cd4[['year3_22']], all_comers_predicted_cd4[['year5_22']])), '25'=round(c(all_comers_predicted_cd4[['year0_25']], all_comers_predicted_cd4[['year1_25']], all_comers_predicted_cd4[['year3_25']], all_comers_predicted_cd4[['year5_25']])), '30'=round(c(all_comers_predicted_cd4[['year0_30']], all_comers_predicted_cd4[['year1_30']], all_comers_predicted_cd4[['year3_30']], all_comers_predicted_cd4[['year5_30']])), '40'=round(c(all_comers_predicted_cd4[['year0_40']], all_comers_predicted_cd4[['year1_40']], all_comers_predicted_cd4[['year3_40']], all_comers_predicted_cd4[['year5_40']])),check.names=FALSE) latex(predicted_df, file='', rowname=NULL, caption='Predicted CD4 values by year and BMI cut-off. These should correspond to the previous graph',where='!h') @ \begin{figure} \caption{Alternate version of above in ggplot2}\label{fig:model_ggplot} <>= ## Alternate version in ggplot2 cd4_by_time_model_based_figure(data=tmp, mod=lme_mod_w_int_ar1_wo_outliers, method='ggplot2') @ \end{figure} \clearpage \subsection{Virologically suppressed cohort using original definition} % ---------------------------------------------------------------------------------- % % Virally suppressed cohort % % ---------------------------------------------------------------------------------- % <>= #---------------------------------------# # Virologic suppression: # # Defined as whether subject spent # # more than 50% of the follow-up time # #---------------------------------------# # Original cohort bmi$vl_lt_400 <- ifelse(!is.na(bmi$hiv1_rna) & bmi$hiv1_rna < 400,1, ifelse(!is.na(bmi$hiv1_rna),0,NA)) vl_suppression_50 <- lapply(split(subset(bmi,select=c(cfar_pid,vl_lt_400,age)),bmi$cfar_pid),FUN=function(y){ visit_interval <- c(diff(y[,"age"]),NA) sum.vl.supp.1 <- sum(visit_interval[y[,"vl_lt_400"] == 1],na.rm=TRUE) sum.vl.supp.0 <- sum(visit_interval[y[,"vl_lt_400"] == 0],na.rm=TRUE) if(sum.vl.supp.1 > sum.vl.supp.0){ vl_suppression_50 <- as.numeric(TRUE) } else { vl_suppression_50 <- as.numeric(FALSE) } df <- data.frame('cfar_pid'=y[1,'cfar_pid'], 'sum.vl.supp.1'=sum.vl.supp.1, 'sum.vl.supp.0'=sum.vl.supp.0, 'vl_suppression_50'=vl_suppression_50) return(df) }) vl_suppression_50.df <- do.call("rbind",vl_suppression_50) vl_suppression_50.df <- as.data.frame(vl_suppression_50.df) bmi$vl_suppression_50 <- vl_suppression_50.df[match(bmi$cfar_pid, vl_suppression_50.df$cfar_pid,nomatch=NA),'vl_suppression_50'] # Creating factor for vl_suppression_at_12_months bmi$vl_suppression_50.factor <- factor(bmi$vl_suppression_50, levels=c(0,1), labels=c('Suppressed <= 50 % of the year', 'Suppressed > 50 % of the year')) vl_suppression_50_alt <- lapply(split(subset(bmi, select=c(cfar_pid, vl_lt_400)),bmi$cfar_pid),FUN=function(y){ suppressed <- sum(y[,'vl_lt_400'],na.rm=TRUE) unsuppressed <- as.integer(sum(y[,'vl_lt_400'] == 0,na.rm=TRUE)) vl_suppression_50_alt <- ifelse(suppressed > unsuppressed,1, ifelse(!is.na(suppressed) & !is.na(unsuppressed),0,NA)) df <- data.frame('cfar_pid'=y[1,'cfar_pid'], 'vl_suppression_50_alt'=vl_suppression_50_alt) return(df) }) vl_suppression_50_alt.df <- do.call('rbind',vl_suppression_50_alt) bmi$vl_suppression_50_alt <- vl_suppression_50_alt.df[match(bmi$cfar_pid, vl_suppression_50_alt.df$cfar_pid,nomatch=NA),'vl_suppression_50_alt'] # Creating factor for vl_suppression_at_12_months bmi$vl_suppression_50_alt.factor <- factor(bmi$vl_suppression_50_alt, levels=c(0,1), labels=c('Suppressed <= 50 % of the year', 'Suppressed > 50 % of the year')) # Sensitivity cohort bmi_sens$vl_lt_400 <- ifelse(!is.na(bmi_sens$hiv1_rna) & bmi_sens$hiv1_rna < 400,1, ifelse(!is.na(bmi_sens$hiv1_rna),0,NA)) vl_suppression_50_sens <- lapply(split(subset(bmi_sens,select=c(cfar_pid,vl_lt_400,age)),bmi_sens$cfar_pid),FUN=function(y){ visit_interval <- c(diff(y[,"age"]),NA) sum.vl.supp.1 <- sum(visit_interval[y[,"vl_lt_400"] == 1],na.rm=TRUE) sum.vl.supp.0 <- sum(visit_interval[y[,"vl_lt_400"] == 0],na.rm=TRUE) if(sum.vl.supp.1 > sum.vl.supp.0){ vl_suppression_50 <- as.numeric(TRUE) } else { vl_suppression_50 <- as.numeric(FALSE) } df <- data.frame('cfar_pid'=y[1,'cfar_pid'], 'sum.vl.supp.1'=sum.vl.supp.1, 'sum.vl.supp.0'=sum.vl.supp.0, 'vl_suppression_50'=vl_suppression_50) return(df) }) vl_suppression_50.df_sens <- do.call("rbind",vl_suppression_50_sens) vl_suppression_50.df_sens <- as.data.frame(vl_suppression_50.df_sens) bmi_sens$vl_suppression_50 <- vl_suppression_50.df_sens[match(bmi_sens$cfar_pid, vl_suppression_50.df_sens$cfar_pid,nomatch=NA),'vl_suppression_50'] # Creating factor for vl_suppression_at_12_months bmi_sens$vl_suppression_50.factor <- factor(bmi_sens$vl_suppression_50, levels=c(0,1), labels=c('Suppressed <= 50 % of the year','Suppressed > 50 % of the year')) vl_suppression_50_alt_sens <- lapply(split(subset(bmi_sens, select=c(cfar_pid, vl_lt_400)),bmi_sens$cfar_pid),FUN=function(y){ suppressed <- sum(y[,'vl_lt_400'],na.rm=TRUE) unsuppressed <- as.integer(sum(y[,'vl_lt_400'] == 0,na.rm=TRUE)) vl_suppression_50_alt <- ifelse(suppressed > unsuppressed,1, ifelse(!is.na(suppressed) & !is.na(unsuppressed),0,NA)) df <- data.frame('cfar_pid'=y[1,'cfar_pid'], 'vl_suppression_50_alt'=vl_suppression_50_alt) return(df) }) vl_suppression_50_alt.df_sens <- do.call('rbind',vl_suppression_50_alt_sens) bmi_sens$vl_suppression_50_alt <- vl_suppression_50_alt.df_sens[match(bmi_sens$cfar_pid, vl_suppression_50_alt.df_sens$cfar_pid,nomatch=NA),'vl_suppression_50_alt'] # Creating factor for vl_suppression_at_12_months bmi_sens$vl_suppression_50_alt.factor <- factor(bmi_sens$vl_suppression_50_alt, levels=c(0,1),labels=c('Suppressed <= 50 % of the year','Suppressed > 50 % of the year')) @ <>= tbl1_suppressed <- summary(bmi_at_age_fhaart_categorized_larger.factor ~ age_fhaart + sex.factor + race.factor + race_for_reviewers.factor + evidence_of_idu.factor + year_of_first_haart_start + pi_usage.factor + regimen_class_in_fhaart.factor + prior_ade.factor + has_hcv_diagnosis.factor + death.factor + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + bmi_at_age_fhaart + cohort, data=bmi, method='reverse',overall=TRUE,subset=!duplicated(cfar_pid) & vl_suppression_50 == 1, test=TRUE) latex(tbl1_suppressed, file='', long=TRUE,exclude1=FALSE,longtable=TRUE,caption='Descriptive statistics by categorized BMI at first HAART in those virally suppressed $>$ 50 percent of their follow up time', digits=2, size='tiny', landscape=TRUE,where='!h',prtest='P') @ \clearpage <>= #-----------------------------------------# # Generalized least squares # #-----------------------------------------# # With AR1 correlated data gls_mod_wo_int_suppressed_ar1 <- Gls(fmla_wo_int, data=bmi, subset=vl_suppression_50 == 1, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) gls_mod_w_int_suppressed_ar1 <- Gls(fmla_w_int, data=bmi, subset=vl_suppression_50 == 1, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) # In data for sensitivity analysis gls_mod_wo_int_suppressed_ar1_sens <- Gls(fmla_wo_int, data=bmi_sens, subset=vl_suppression_50 == 1, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) gls_mod_w_int_suppressed_ar1_sens <- Gls(fmla_w_int, data=bmi_sens, subset=vl_suppression_50 == 1, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) #-----------------------------------------# # Mixed effects linear model # #-----------------------------------------# # AR1 correlation in lme tmp <- subset(bmi, vl_suppression_50 == 1, select=c('cfar_pid','vl_suppression_50','cd4_count_sqrt','bmi_closest','t_lab','sex.factor','race.factor', 'age_fhaart','regimen_class_in_fhaart.factor','year_of_first_haart_start','cd4_count_at_age_fhaart_sqrt', 'loghiv1_rna_at_age_fhaart','cohort')) lme_mod_w_int_suppressed_ar1 <- lme(fmla_w_int, random=~1|cfar_pid, data=tmp, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) # In data for sensitivity analysis tmp_sens <- subset(bmi_sens, vl_suppression_50 == 1, select=c('cfar_pid','vl_suppression_50','cd4_count_sqrt','bmi_closest','t_lab','sex.factor','race.factor', 'age_fhaart','regimen_class_in_fhaart.factor','year_of_first_haart_start','cd4_count_at_age_fhaart_sqrt', 'loghiv1_rna_at_age_fhaart','cohort')) lme_mod_w_int_suppressed_ar1_sens <- lme(fmla_w_int, random=~1|cfar_pid, data=tmp_sens, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) @ <>= # The predicted values does not include those records with missing values so need to mark them to remove them tmp$missing <- apply(tmp, 1, FUN=function(y){ifelse(any(is.na(y)),1,0)}) tmp <- subset(tmp, missing == 0) tmp$residuals <- lme_mod_w_int_suppressed_ar1[['residuals']][,'fixed'] lme_mod_w_int_suppressed_ar1_wo_outliers <- lme(fmla_w_int_wo_outliers, random=~1|cfar_pid, data=tmp, subset=abs(tmp$residuals) <= 6, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) # In data for sensitivity analysis tmp_sens$missing <- apply(tmp_sens, 1, FUN=function(y){ifelse(any(is.na(y)),1,0)}) tmp_sens <- subset(tmp_sens, missing == 0) tmp_sens$residuals <- lme_mod_w_int_suppressed_ar1_sens[['residuals']][,'fixed'] lme_mod_w_int_suppressed_ar1_wo_outliers_sens <- lme(fmla_w_int_wo_outliers, random=~1|cfar_pid, data=tmp_sens, subset=abs(tmp_sens$residuals) <= 6, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) @ <>= suppressed_mod_assessment <- model_fit_assessment(gls_mod=gls_mod_w_int_suppressed_ar1, lme_mod=lme_mod_w_int_suppressed_ar1, data=bmi[which(bmi$vl_suppression_50 == 1),], fmla=fmla_w_int) cor_df_suppressed <- suppressed_mod_assessment[['cor_df']] dt_suppressed <- suppressed_mod_assessment[['dt']] dr_suppressed <- suppressed_mod_assessment[['dr']] dt_wo_outliers_suppressed <- suppressed_mod_assessment[['dt_wo_outliers']] dr_wo_outliers_suppressed <- suppressed_mod_assessment[['dr_wo_outliers']] lm_residuals_wo_outliers_suppressed <- suppressed_mod_assessment[['lm_residuals_wo_outliers']] # In data for sensivitiy analysis suppressed_mod_assessment_sens <- model_fit_assessment(gls_mod=gls_mod_w_int_suppressed_ar1_sens, lme_mod=lme_mod_w_int_suppressed_ar1_sens, data=bmi_sens[which(bmi$vl_suppression_50 == 1),], fmla=fmla_w_int) cor_df_suppressed_sens <- suppressed_mod_assessment_sens[['cor_df']] dt_suppressed_sens <- suppressed_mod_assessment_sens[['dt']] dr_suppressed_sens <- suppressed_mod_assessment_sens[['dr']] dt_wo_outliers_suppressed_sens <- suppressed_mod_assessment_sens[['dt_wo_outliers']] dr_wo_outliers_suppressed_sens <- suppressed_mod_assessment_sens[['dr_wo_outliers']] lm_residuals_wo_outliers_suppressed_sens <- suppressed_mod_assessment_sens[['lm_residuals_wo_outliers']] @ \subsubsection{Assessing appropriate correlation structure in the virally suppressed cohort} \normalsize <>= # Table to evaluate the minimum AIC for each of the models from the virally suppressed cohort. cor_df_suppressed @ \clearpage \begin{figure} \caption{Variogram for assessing correlation structure in the virally suppressed cohort} <>= #pdf('Variogram_wo_points_and_w_lowess_line_capped.pdf') plot(dt_wo_outliers_suppressed, dr_wo_outliers_suppressed,main='',xlab='dt',ylab='dr',type='n',ylim=c(0,30)) lines(lowess(x=dt_wo_outliers_suppressed,y=dr_wo_outliers_suppressed),col='red') abline(h=var(lm_residuals_wo_outliers_suppressed),col='blue') #dev.off() @ \end{figure} \begin{figure} \caption{Autocorrelation function graphs for the lme and gls models in the virally suppressed cohort} <>= # Autocorrelation function graphs acf1_suppressed <- plot(ACF(lme_mod_w_int_suppressed_ar1, maxLag=10,resType='n'),alpha=0.01,main='lme model',cex.main=0.65) acf2_suppressed <- plot(ACF(gls_mod_w_int_suppressed_ar1,resType="normalized"),alpha=0.05,main='gls model',cex.main=0.65) print(acf1_suppressed, position=c(0,0.5,1,1.0),more=TRUE) print(acf2_suppressed, position=c(0,0,1,0.45)) @ \end{figure} \begin{figure} \caption{QQ plots for all data in the virally suppressed cohort and for data with abs(residuals) $>$ 6 removed.} <>= par(mfrow=c(2,2),mar=c(1,1,1,1),oma=c(1,1,1,1)) # QQ plot #pdf('test_qq.pdf') qq1_suppressed <- qqnorm(lme_mod_w_int_suppressed_ar1, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1),cex.lab=0.65) #dev.off() qq2_suppressed <- qqnorm(lme_mod_w_int_suppressed_ar1_wo_outliers, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1), xlab='Normalized residuals (outliers removed)',cex.lab=0.65) print(qq1_suppressed, position=c(0,0.5,1,1.0),more=TRUE) print(qq2_suppressed, position=c(0,0,1,0.45)) @ \end{figure} \clearpage \subsubsection{Comparing model output in virally suppressed cohort} \tiny <>= # ---------------------------------------- # # Comparing the different model outputs: # # ---------------------------------------- # # Gls model using AR1 correlation form=~1|ID #summary(gls_mod_w_int_suppressed_ar1) gls_coeff_suppressed <- gls_mod_w_int_suppressed_ar1$coefficients # lme model using AR1 correlation form=~1|ID #summary(lme_mod_w_int_suppressed_ar1) lme_coeff_suppressed <- lme_mod_w_int_suppressed_ar1$coefficients$fixed cbind(format(round(lme_coeff_suppressed,5),scientific=FALSE),format(round(gls_coeff_suppressed,5),scientific=FALSE)) @ \clearpage \subsubsection{Evaluating models with interaction terms in virally suppressed cohort} <>= # ----------------------------- # # Evaluating significance of # # interaction terms # # ----------------------------- # # ANOVA table for gls models anova(gls_mod_w_int_suppressed_ar1, tol=1e-13) # ANOVA table for lme models anova(lme_mod_w_int_suppressed_ar1) @ <>= # Bryan wants the individual lines on the graphs to only span the BMIs found for that race/sex group suppressed_mod_limits <- graph_limits(bmi[which(bmi$vl_suppression_50 == 1),]) bmi_rge_wm_suppressed <- suppressed_mod_limits[['bmi_rge_wm']] bmi_rge_wf_suppressed <- suppressed_mod_limits[['bmi_rge_wf']] bmi_rge_nwm_suppressed <- suppressed_mod_limits[['bmi_rge_nwm']] bmi_rge_nwf_suppressed <- suppressed_mod_limits[['bmi_rge_nwf']] tmp_rge_wm_suppressed <- suppressed_mod_limits[['tmp_rge_wm']] tmp_rge_wf_suppressed <- suppressed_mod_limits[['tmp_rge_wf']] tmp_rge_nwm_suppressed <- suppressed_mod_limits[['tmp_rge_nwm']] tmp_rge_nwf_suppressed <- suppressed_mod_limits[['tmp_rge_nwf']] tmp_w_residuals_rge_wm_suppressed <- suppressed_mod_limits[['tmp_w_residuals_rge_wm']] tmp_w_residuals_rge_wf_suppressed <- suppressed_mod_limits[['tmp_w_residuals_rge_wf']] tmp_w_residuals_rge_nwm_suppressed <- suppressed_mod_limits[['tmp_w_residuals_rge_nwm']] tmp_w_residuals_rge_nwf_suppressed <- suppressed_mod_limits[['tmp_w_residuals_rge_nwf']] overall_x_rge_suppressed <- suppressed_mod_limits[['overall_x_rge']] @ <>= # From Gls gls_predict_suppressed <- Predict(gls_mod_w_int_suppressed_ar1,bmi_closest,sex.factor,race.factor) # Transforming CD4 back from square root transformed gls_predict_suppressed$yhat <- gls_predict_suppressed$yhat^2 gls_predict_suppressed$lower <- gls_predict_suppressed$lower^2 gls_predict_suppressed$upper <- gls_predict_suppressed$upper^2 rge_suppressed <- range(c(gls_predict_suppressed$lower,gls_predict_suppressed$upper)) x.rge_suppressed <- range(c(gls_predict_suppressed$bmi_closest)) x.rge_suppressed <- c(x.rge_suppressed[1],quantile(gls_predict_suppressed$bmi_closest,prob=0.90)) x.rge_gls_alt_suppressed <- x.rge_suppressed rge_gls_alt_suppressed <- rge_suppressed # ------------------------------ # # From lme # Getting predicted CD4 count values (square root transformed) cohort_mode_tbl <- sort(table(bmi$cohort),decreasing=TRUE) cohort_mode <- names(cohort_mode_tbl)[1] lme_predict_data_suppressed <- data.frame('bmi_closest'=rep(seq(min(bmi$bmi_closest[which(bmi$vl_suppression_50 == 1)],na.rm=TRUE), max(bmi$bmi_closest[which(bmi$vl_suppression_50 == 1)],na.rm=TRUE), length.out=100),times=4), 't_lab'=rep(median(bmi$t_lab[which(bmi$vl_suppression_50 == 1)]),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(bmi$age_fhaart[which(!duplicated(bmi$cfar_pid) & bmi$vl_suppression_50 == 1)]), times=400), 'regimen_class_in_fhaart.factor'=rep(levels(bmi$regimen_class_in_fhaart.factor[which(bmi$vl_suppression_50 == 1)])[2], times=400), # 'regimen_class_in_fhaart.factor'=rep(levels(bmi$regimen_class_in_fhaart.factor[which(bmi$vl_suppression_50 == 1)])[1], # times=400), 'year_of_first_haart_start'=rep(median(bmi$year_of_first_haart_start[which(bmi$vl_suppression_50 == 1)]), times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(bmi$cd4_count_at_age_fhaart_sqrt[which(!duplicated(bmi$cfar_pid) & bmi$vl_suppression_50 == 1)]), times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(bmi$loghiv1_rna_at_age_fhaart[which(!duplicated(bmi$cfar_pid) & bmi$vl_suppression_50 == 1)], na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode,times=400)) # 'cohort'=rep(levels(bmi$cohort[which(bmi$vl_suppression_50 == 1)])[1],times=400)) lme_predict_suppressed <- predict(lme_mod_w_int_suppressed_ar1, newdata=lme_predict_data_suppressed,na.action=na.omit,level=0) ## Transforming CD4 back from square root transformed lme_predict_suppressed <- lme_predict_suppressed^2 @ \begin{figure} \caption{Predicted CD4 count by BMI in the virally suppressed cohort using lme}\label{fig:cd4_by_bmi_lme_alt_suppressed} <>= cd4_by_bmi_figure(mod_predict=lme_predict_suppressed, predict_data=lme_predict_data_suppressed, fig_limits= suppressed_mod_limits, plot_method='base',subset_method='lme') @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI in the virally suppressed cohort using lme and omitting those with abs(residuals) $>$ 6}\label{fig:cd4_by_bmi_lme_alt_wo_outliers_suppressed} <>= # Getting predicted CD4 count values (square root transformed) cohort_mode_tbl <- sort(table(tmp$cohort),decreasing=TRUE) cohort_mode <- names(cohort_mode_tbl)[1] lme_predict_data_wo_outliers_suppressed <- data.frame('bmi_closest'=rep(seq(min(tmp$bmi_closest[which(abs(tmp$residuals) <= 6)],na.rm=TRUE), max(tmp$bmi_closest[which(abs(tmp$residuals) <= 6)],na.rm=TRUE), length.out=100),times=4), 't_lab'=rep(median(tmp$t_lab),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(tmp$age_fhaart[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)]), times=400), 'regimen_class_in_fhaart.factor'= rep('NNRTI-based', times=400), # 'regimen_class_in_fhaart.factor'= rep(sort(unique(tmp$regimen_class_in_fhaart.factor[which(abs(tmp$residuals) <= 6)]))[1], times=400), 'year_of_first_haart_start'=rep(median(tmp$year_of_first_haart_start[which(abs(tmp$residuals) <= 6)]),times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp$cd4_count_at_age_fhaart_sqrt[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)]),times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp$loghiv1_rna_at_age_fhaart[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode,times=400)) # 'cohort'=rep(sort(unique(tmp$cohort[which(abs(tmp$residuals) <= 6)]))[1],times=400)) lme_predict_wo_outliers_suppressed <- predict(lme_mod_w_int_suppressed_ar1_wo_outliers, newdata=lme_predict_data_wo_outliers_suppressed,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed lme_predict_wo_outliers_suppressed <- lme_predict_wo_outliers_suppressed^2 cd4_by_bmi_figure(mod_predict=lme_predict_wo_outliers_suppressed, predict_data=lme_predict_data_wo_outliers_suppressed, fig_limits=main_mod_limits, plot_method='base', subset_method='wo_outliers') @ \end{figure} <>= cd4_by_year_combinednew1_suppressed <- descriptive_graph_mgmt(subset(bmi, vl_suppression_50 == 1)) @ \begin{figure} \caption{Descriptive graph of CD4 count by year since ART initiation and BMI category. Yearly CD4 counts and BMI measurements were defined as the closest CD4 count within +/- 180 days. A subject could switch BMI categories from year to year. Points at a given year and BMI category are the average of all CD4 counts within that year and BMI category. This graph is on the cohort of subjects who were virally suppressed at least 50\% of the time.} <>= ggplot(cd4_by_year_combinednew1_suppressed, aes(x=Time,y=CD4_count, group=BMI, colour=BMI)) + geom_line() + xlab('Time since ART initiation (years)') + ylab('CD4 count') + scale_colour_manual(values=c('black','red','green','blue','yellow','purple')) @ \end{figure} \clearpage \begin{figure} \caption{Model-based graph of CD4 count by year since ART initiation and by BMI cut points corresponding to common BMI categories. Predicted CD4 values were calculated over 4 years for each of the cut points listed in the legend. This is on the cohort of subjects who were virally suppressed at least 50\% of the time and with outliers removed.} <>= cd4_by_time_model_based_figure(data=tmp, mod=lme_mod_w_int_suppressed_ar1_wo_outliers, method='base') @ \end{figure} \begin{figure} \caption{Alternate version of above in ggplot2 in those suppressed $>$ 50 percent of the time} <>= ## Alternate version in ggplot2 cd4_by_time_model_based_figure(data=tmp, mod=lme_mod_w_int_suppressed_ar1_wo_outliers, method='ggplot2') @ \end{figure} \clearpage \subsection{An alternate method for determining those who were virologically suppressed more than 50 \% of the time.} \normalsize Results are presented in this section in which the cohort of those who were virologically suppressed more than 50\% of the time was determined by simply counting the number of records in which a suppressed viral load measurement occurred and comparing it to the number of records in which a detectable viral load measurement occurred. Those who had more measurements in which viral load was suppressed were considered suppressed more than 50\% of their follow-up time. <>= tbl1_suppressed_alt <- summary(bmi_at_age_fhaart_categorized_larger.factor ~ age_fhaart + sex.factor + race.factor + race_for_reviewers.factor + evidence_of_idu.factor + year_of_first_haart_start + pi_usage.factor + regimen_class_in_fhaart.factor + prior_ade.factor + has_hcv_diagnosis.factor + death.factor + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + bmi_at_age_fhaart + cohort, data=bmi, method='reverse',overall=TRUE,subset=!duplicated(cfar_pid) & vl_suppression_50_alt == 1, test=TRUE) latex(tbl1_suppressed_alt, file='', long=TRUE,exclude1=FALSE,longtable=TRUE,caption='Descriptive statistics by categorized BMI at first HAART in those virally suppressed $>$ 50 percent of their follow up time using the alternate definition', digits=2, size='tiny', landscape=TRUE,where='!h',prtest='P') @ \clearpage <>= #-----------------------------------------# # Generalized least squares # #-----------------------------------------# # With AR1 correlated data gls_mod_wo_int_suppressed_ar1_alt <- Gls(fmla_wo_int, data=bmi, subset=vl_suppression_50_alt == 1, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) gls_mod_w_int_suppressed_ar1_alt <- Gls(fmla_w_int, data=bmi, subset=vl_suppression_50_alt == 1, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) # For sensitivity analysis gls_mod_wo_int_suppressed_ar1_alt_sens <- Gls(fmla_wo_int, data=bmi_sens, subset=vl_suppression_50_alt == 1,corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) gls_mod_w_int_suppressed_ar1_alt_sens <- Gls(fmla_w_int, data=bmi_sens, subset=vl_suppression_50_alt == 1,corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) #-----------------------------------------# # Mixed effects linear model # #-----------------------------------------# # AR1 correlation in lme tmp <- subset(bmi, vl_suppression_50_alt == 1, select=c('cfar_pid','vl_suppression_50_alt','cd4_count_sqrt','bmi_closest','t_lab','sex.factor','race.factor', 'age_fhaart','regimen_class_in_fhaart.factor','year_of_first_haart_start','cd4_count_at_age_fhaart_sqrt', 'loghiv1_rna_at_age_fhaart','cohort')) lme_mod_w_int_suppressed_ar1_alt <- lme(fmla_w_int, random=~1|cfar_pid, data=tmp, subset=vl_suppression_50_alt == 1, na.action=na.omit, correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) # For sensitivity analysis tmp_sens <- subset(bmi_sens, vl_suppression_50_alt == 1, select=c('cfar_pid','vl_suppression_50_alt','cd4_count_sqrt','bmi_closest','t_lab','sex.factor','race.factor', 'age_fhaart','regimen_class_in_fhaart.factor','year_of_first_haart_start','cd4_count_at_age_fhaart_sqrt', 'loghiv1_rna_at_age_fhaart','cohort')) lme_mod_w_int_suppressed_ar1_alt_sens <- lme(fmla_w_int, random=~1|cfar_pid, data=tmp_sens, subset=vl_suppression_50_alt == 1, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) @ <>= # The predicted values does not include those records with missing values so need to mark them to remove them tmp$missing <- apply(tmp, 1, FUN=function(y){ifelse(any(is.na(y)),1,0)}) tmp <- subset(tmp, missing == 0) # tmp_w_residuals_suppressed_alt <- subset(tmp, vl_suppression_50_alt == 1) tmp$residuals <- lme_mod_w_int_suppressed_ar1_alt[['residuals']][,'fixed'] lme_mod_w_int_suppressed_ar1_wo_outliers_alt <- lme(fmla_w_int_wo_outliers, random=~1|cfar_pid, data=tmp, subset=abs(residuals) <= 6, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) # For sensitivity analysis tmp_sens$missing <- apply(tmp_sens, 1, FUN=function(y){ifelse(any(is.na(y)),1,0)}) tmp_sens <- subset(tmp_sens, missing == 0) # tmp_w_residuals_suppressed_alt <- subset(tmp, vl_suppression_50_alt == 1) tmp_sens$residuals <- lme_mod_w_int_suppressed_ar1_alt_sens[['residuals']][,'fixed'] lme_mod_w_int_suppressed_ar1_wo_outliers_alt_sens <- lme(fmla_w_int_wo_outliers, random=~1|cfar_pid, data=tmp_sens, subset=abs(residuals) <= 6,na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) @ <>= suppressed_mod_assessment_alt <- model_fit_assessment(gls_mod=gls_mod_w_int_suppressed_ar1_alt, lme_mod=lme_mod_w_int_suppressed_ar1_alt, data=bmi[which(bmi$vl_suppression_50_alt == 1),], fmla=fmla_w_int) cor_df_suppressed_alt <- suppressed_mod_assessment_alt[['cor_df']] dt_suppressed_alt <- suppressed_mod_assessment_alt[['dt_wo_outliers']] dr_suppressed_alt <- suppressed_mod_assessment_alt[['dr_wo_outliers']] lm_residuals_wo_outliers_suppressed_alt <- suppressed_mod_assessment_alt[['lm_residuals_wo_outliers']] # For sensitivity analysis suppressed_mod_assessment_alt_sens <- model_fit_assessment(gls_mod=gls_mod_w_int_suppressed_ar1_alt_sens, lme_mod=lme_mod_w_int_suppressed_ar1_alt_sens, data=bmi_sens[which(bmi_sens$vl_suppression_50_alt == 1),], fmla=fmla_w_int) cor_df_suppressed_alt_sens <- suppressed_mod_assessment_alt_sens[['cor_df']] dt_suppressed_alt_sens <- suppressed_mod_assessment_alt_sens[['dt_wo_outliers']] dr_suppressed_alt_sens <- suppressed_mod_assessment_alt_sens[['dr_wo_outliers']] lm_residuals_wo_outliers_suppressed_alt_sens <- suppressed_mod_assessment_alt_sens[['lm_residuals_wo_outliers']] @ \subsubsection{Assessing appropriate correlation structure in the virally suppressed cohort using the alternate definition for viral suppression} <>= # Table to evaluate the minimum AIC for each of the models from the virally suppressed cohort. cor_df_suppressed_alt @ \clearpage \begin{figure} \caption{Variogram for assessing correlation structure in the virally suppressed cohort using the alternate definition for viral suppression} <>= #pdf('Variogram_wo_points_and_w_lowess_line_capped.pdf') plot(dt_suppressed_alt, dr_suppressed_alt,main='',xlab='dt',ylab='dr',type='n',ylim=c(0,30)) lines(lowess(x=dt_suppressed_alt,y=dr_suppressed_alt),col='red') abline(h=var(lm_residuals_wo_outliers_suppressed_alt),col='blue') #dev.off() @ \end{figure} \begin{figure} \caption{Autocorrelation function graphs for the lme and gls models in the virally suppressed cohort using the alternate definition for viral suppression} <>= # Autocorrelation function graphs acf1_suppressed_alt <- plot(ACF(lme_mod_w_int_suppressed_ar1_alt, maxLag=10,resType='n'),alpha=0.01,main='lme model',cex.main=0.65) acf2_suppressed_alt <- plot(ACF(gls_mod_w_int_suppressed_ar1_alt,resType="normalized"),alpha=0.05,main='gls model',cex.main=0.65) print(acf1_suppressed_alt, position=c(0,0.5,1,1.0),more=TRUE) print(acf2_suppressed_alt, position=c(0,0,1,0.45)) @ \end{figure} \begin{figure} \caption{QQ plots for all data and when outliers were removed in the virally suppressed cohort using the alternate definition for viral suppression} <>= par(mfrow=c(2,2),mar=c(1,1,1,1),oma=c(1,1,1,1)) # QQ plot #pdf('test_qq.pdf') qq1_suppressed_alt <- qqnorm(lme_mod_w_int_suppressed_ar1_alt, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1),cex.lab=0.65) #dev.off() qq2_suppressed_alt <- qqnorm(lme_mod_w_int_suppressed_ar1_wo_outliers_alt, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1), xlab='Normalized residuals (outliers removed)',cex.lab=0.65) print(qq1_suppressed_alt, position=c(0,0.5,1,1.0),more=TRUE) print(qq2_suppressed_alt, position=c(0,0,1,0.45)) @ \end{figure} \clearpage \subsubsection{Comparing model output in virally suppressed cohort using the alternate definition for viral suppression} \tiny <>= # ---------------------------------------- # # Comparing the different model outputs: # # ---------------------------------------- # # Gls model using AR1 correlation form=~1|ID #summary(gls_mod_w_int_suppressed_ar1) gls_coeff_suppressed_alt <- gls_mod_w_int_suppressed_ar1_alt$coefficients # lme model using AR1 correlation form=~1|ID #summary(lme_mod_w_int_suppressed_ar1) lme_coeff_suppressed_alt <- lme_mod_w_int_suppressed_ar1_alt$coefficients$fixed cbind(format(round(lme_coeff_suppressed_alt,5),scientific=FALSE),format(round(gls_coeff_suppressed_alt,5),scientific=FALSE)) @ \clearpage \subsubsection{Evaluating models with interaction terms in virally suppressed cohort using the alternate definition for viral suppression} <>= # ----------------------------- # # Evaluating significance of # # interaction terms # # ----------------------------- # # ANOVA table for gls models anova(gls_mod_w_int_suppressed_ar1_alt, tol=1e-13) # ANOVA table for lme models anova(lme_mod_w_int_suppressed_ar1_alt) @ <>= # Bryan wants the individual lines on the graphs to only span the BMIs found for that race/sex group suppressed_mod_limits_alt <- graph_limits(bmi[which(bmi$vl_suppression_50_alt == 1),]) bmi_rge_wm_suppressed_alt <- suppressed_mod_limits_alt[['bmi_rge_wm']] bmi_rge_wf_suppressed_alt <- suppressed_mod_limits_alt[['bmi_rge_wf']] bmi_rge_nwm_suppressed_alt <- suppressed_mod_limits_alt[['bmi_rge_nwm']] bmi_rge_nwf_suppressed_alt <- suppressed_mod_limits_alt[['bmi_rge_nwf']] tmp_rge_wm_suppressed_alt <- suppressed_mod_limits_alt[['tmp_rge_wm']] tmp_rge_wf_suppressed_alt <- suppressed_mod_limits_alt[['tmp_rge_wf']] tmp_rge_nwm_suppressed_alt <- suppressed_mod_limits_alt[['tmp_rge_nwm']] tmp_rge_nwf_suppressed_alt <- suppressed_mod_limits_alt[['tmp_rge_nwf']] tmp_w_residuals_rge_wm_suppressed_alt <- suppressed_mod_limits_alt[['tmp_w_residuals_rge_wm']] tmp_w_residuals_rge_wf_suppressed_alt <- suppressed_mod_limits_alt[['tmp_w_residuals_rge_wf']] tmp_w_residuals_rge_nwm_suppressed_alt <- suppressed_mod_limits_alt[['tmp_w_residuals_rge_nwm']] tmp_w_residuals_rge_nwf_suppressed_alt <- suppressed_mod_limits_alt[['tmp_w_residuals_rge_nwf']] overall_x_rge_suppressed_alt <- suppressed_mod_limits_alt[['overall_x_rge']] @ <>= # From Gls gls_predict_suppressed_alt <- Predict(gls_mod_w_int_suppressed_ar1_alt,bmi_closest,sex.factor,race.factor) # Transforming CD4 back from square root transformed gls_predict_suppressed_alt$yhat <- gls_predict_suppressed_alt$yhat^2 gls_predict_suppressed_alt$lower <- gls_predict_suppressed_alt$lower^2 gls_predict_suppressed_alt$upper <- gls_predict_suppressed_alt$upper^2 # ------------------------------ # # From lme # Getting predicted CD4 count values (square root transformed) cohort_mode_tbl <- sort(table(tmp$cohort),decreasing=TRUE) cohort_mode <- names(cohort_mode_tbl)[1] lme_predict_data_suppressed_alt <- data.frame('bmi_closest'=rep(seq(min(tmp$bmi_closest,na.rm=TRUE), max(tmp$bmi_closest,na.rm=TRUE), length.out=100),times=4), 't_lab'=rep(median(tmp$t_lab),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(tmp$age_fhaart[!duplicated(tmp$cfar_pid)]), times=400), 'regimen_class_in_fhaart.factor'=rep('NNRTI-based', times=400), # 'regimen_class_in_fhaart.factor'=rep(levels(tmp$regimen_class_in_fhaart.factor)[1], times=400), 'year_of_first_haart_start'=rep(median(tmp$year_of_first_haart_start), times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp$cd4_count_at_age_fhaart_sqrt[!duplicated(tmp$cfar_pid)]), times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp$loghiv1_rna_at_age_fhaart[!duplicated(tmp$cfar_pid)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode,times=400)) # 'cohort'=rep(levels(tmp$cohort)[1],times=400)) lme_predict_suppressed_alt <- predict(lme_mod_w_int_suppressed_ar1_alt, newdata=lme_predict_data_suppressed_alt,na.action=na.omit,level=0) ## Transforming CD4 back from square root transformed lme_predict_suppressed_alt <- lme_predict_suppressed_alt^2 @ \begin{figure} \caption{Predicted CD4 count by BMI in the virally suppressed cohort using the alternate suppression definition and using lme}\label{fig:cd4_by_bmi_lme_alt_suppressed_alt} <>= cd4_by_bmi_figure(mod_predict=lme_predict_suppressed_alt, predict_data=lme_predict_data_suppressed_alt, fig_limits= suppressed_mod_limits_alt, plot_method='base', subset_method='lme') @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI in the virally suppressed cohort using the alternate definition for suppression and using lme and omitting those with abs(residuals) $>$ 6}\label{fig:cd4_by_bmi_lme_alt_wo_outliers_suppressed_alt} <>= cohort_mode_tbl <- sort(table(tmp$cohort),decreasing=TRUE) cohort_mode <- names(cohort_mode_tbl)[1] lme_predict_data_wo_outliers_suppressed_alt <- data.frame('bmi_closest'=rep(seq(min(tmp$bmi_closest[which(abs(tmp$residuals) <= 6)],na.rm=TRUE), max(tmp$bmi_closest[which(abs(tmp$residuals) <= 6)],na.rm=TRUE),length.out=100),times=4), 't_lab'=rep(median(tmp$t_lab[which(abs(tmp$residuals) <= 6)]),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(tmp$age_fhaart[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)]), times=400), 'regimen_class_in_fhaart.factor'=rep('NNRTI-based', times=400), # 'regimen_class_in_fhaart.factor'=rep(sort(unique(tmp$regimen_class_in_fhaart.factor[which(abs(tmp$residuals) <= 6)]))[1], times=400), 'year_of_first_haart_start'=rep(median(tmp$year_of_first_haart_start[which(abs(tmp$residuals) <= 6)]),times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp$cd4_count_at_age_fhaart_sqrt[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)]),times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp$loghiv1_rna_at_age_fhaart[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode,times=400)) # 'cohort'=rep(sort(unique(tmp$cohort[which(abs(tmp$residuals) <= 6)]))[1],times=400)) lme_predict_wo_outliers_suppressed_alt <- predict(lme_mod_w_int_suppressed_ar1_wo_outliers_alt, newdata=lme_predict_data_wo_outliers_suppressed_alt,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed lme_predict_wo_outliers_suppressed_alt <- lme_predict_wo_outliers_suppressed_alt^2 cd4_by_bmi_figure(mod_predict=lme_predict_wo_outliers_suppressed_alt, predict_data=lme_predict_data_wo_outliers_suppressed_alt, fig_limits=main_mod_limits, plot_method='base', subset_method='wo_outliers') @ \end{figure} <>= cd4_by_year_combinednew1_suppressed_alt <- descriptive_graph_mgmt(subset(bmi, vl_suppression_50_alt == 1)) @ \begin{figure} \caption{Descriptive graph of CD4 count by year since ART initiation and BMI category. Yearly CD4 counts and BMI measurements were defined as the closest CD4 count within +/- 180 days. A subject could switch BMI categories from year to year. Points at a given year and BMI category are the average of all CD4 counts within that year and BMI category. This graph is on the cohort of subjects who were virally suppressed at least 50\% of the time using the alternate definition for viral suppression.} <>= ggplot(cd4_by_year_combinednew1_suppressed_alt, aes(x=Time,y=CD4_count, group=BMI, colour=BMI)) + geom_line() + xlab('Time since ART initiation (years)') + ylab('CD4 count') + scale_colour_manual(values=c('black','red','green','blue','yellow','purple')) @ \end{figure} \clearpage \begin{figure} \caption{Model-based graph of CD4 count by year since ART initiation and BMI cut points corresponding to common BMI categories. Predicted CD4 values were calculated over 4 years for each of the cut points listed in the legend. This is on the cohort of subjects who were virally suppressed at least 50\% of the time using the alternate definition for viral suppression and outliers have been removed.} <>= cd4_by_time_model_based_figure(data=tmp, mod=lme_mod_w_int_suppressed_ar1_wo_outliers_alt, method='base') @ \end{figure} \begin{figure} \caption{Alternate version of above in ggplot2 in those suppressed $>$ 50 percent of the time using the alternate definition for viral suppression and outliers removed} <>= ## Alternate version in ggplot2 cd4_by_time_model_based_figure(data=tmp, mod=lme_mod_w_int_suppressed_ar1_wo_outliers_alt, method='ggplot2') @ \end{figure} \clearpage \subsubsection{Censoring follow-up when a subject becomes unsuppressed} \normalsize This analysis censors subjects who meet the CCASAnet definition for virologic failure. Subjects are censored at the first record having a non-missing viral load measurement prior to the failure. Recall that the age of failure was the age at the 2nd of two consecutive detectable viral loads or the age at which one measurement exceeded 1000. <>= # Total follow-up time for someone in the data bmi$max_age <- NA split(bmi$max_age, bmi$cfar_pid) <- lapply(split(bmi$age, bmi$cfar_pid), FUN=function(y){ max(y) }) bmi$total_followup_yrs <- bmi$max_age - bmi$age_fhaart # And then using the CCASAnet definition of virologic failure (2 consecutive readings > 400 or 1 reading > 1000) vl_failure_90 <- virologic_failure_function_adapted(rna=bmi[,c('cfar_pid','age_fhaart','age','hiv1_rna')],detectable_level=400) # Without the 1000 ignores the 1 reading > 1000 criterion bmi$vl_failure <- vl_failure_90[match(bmi$cfar_pid, vl_failure_90$cfar_pid, nomatch=NA),'vl_failure'] bmi$vl_failure <- ifelse(is.na(bmi$vl_failure), 0,bmi$vl_failure) bmi$age_vl_failure <- vl_failure_90[match(bmi$cfar_pid, vl_failure_90$cfar_pid, nomatch=NA),'age_vl_failure'] bmi$age_censor <- vl_failure_90[match(bmi$cfar_pid, vl_failure_90$cfar_pid, nomatch=NA),'age_censor'] # With the 1000 includes the 1 reading > 1000 criterion bmi$vl_failure.1000 <- vl_failure_90[match(bmi$cfar_pid, vl_failure_90$cfar_pid, nomatch=NA),'vl_failure.1000'] bmi$vl_failure.1000 <- ifelse(is.na(bmi$vl_failure.1000), 0,bmi$vl_failure.1000) bmi$age_vl_failure.1000 <- vl_failure_90[match(bmi$cfar_pid, vl_failure_90$cfar_pid, nomatch=NA),'age_vl_failure.1000'] bmi$age_censor.1000 <- vl_failure_90[match(bmi$cfar_pid, vl_failure_90$cfar_pid, nomatch=NA),'age_censor.1000'] bmi$cum_inc_followup_yrs_alt.1000 <- bmi$age_vl_failure.1000 - bmi$age_fhaart bmi$cum_inc_followup_yrs_alt.1000 <- ifelse(is.na(bmi$cum_inc_followup_yrs_alt.1000),bmi$max_age - bmi$age_fhaart,bmi$cum_inc_followup_yrs_alt.1000) bmi$cum_inc_censor.1000 <- ifelse(bmi$vl_failure.1000 == 1 & bmi$age == bmi$age_vl_failure.1000,1,0) cum_inc_unsuppressed_alt.1000 <- cuminc(ftime=bmi$cum_inc_followup_yrs_alt.1000[which(!duplicated(bmi$cfar_pid) & (is.na(bmi$age_censor.1000) | bmi$age_censor.1000 != -99))], fstatus=bmi$vl_failure.1000[which(!duplicated(bmi$cfar_pid) & (is.na(bmi$age_censor.1000) | bmi$age_censor.1000 != -99))], cencode='0',rho=0) which_50_suppressed <- max(which(cum_inc_unsuppressed_alt.1000[[1]]$est <= 0.50)) time_to_50 <- cum_inc_unsuppressed_alt.1000[[1]]$time[which_50_suppressed] # Total follow-up time for someone in the data for sensitivity analysis bmi_sens$max_age <- NA split(bmi_sens$max_age, bmi_sens$cfar_pid) <- lapply(split(bmi_sens$age, bmi_sens$cfar_pid), FUN=function(y){ max(y) }) bmi_sens$total_followup_yrs <- bmi_sens$max_age - bmi_sens$age_fhaart # And then using the CCASAnet definition of virologic failure (2 consecutive readings > 400 or 1 reading > 1000) vl_failure_90_sens <- virologic_failure_function_adapted(rna=bmi_sens[,c('cfar_pid','age_fhaart','age','hiv1_rna')],detectable_level=400) # Without the 1000 ignores the 1 reading > 1000 criterion bmi_sens$vl_failure <- vl_failure_90_sens[match(bmi_sens$cfar_pid, vl_failure_90_sens$cfar_pid, nomatch=NA),'vl_failure'] bmi_sens$vl_failure <- ifelse(is.na(bmi_sens$vl_failure), 0,bmi_sens$vl_failure) bmi_sens$age_vl_failure <- vl_failure_90_sens[match(bmi_sens$cfar_pid, vl_failure_90_sens$cfar_pid, nomatch=NA),'age_vl_failure'] bmi_sens$age_censor <- vl_failure_90_sens[match(bmi_sens$cfar_pid, vl_failure_90_sens$cfar_pid, nomatch=NA),'age_censor'] # With the 1000 includes the 1 reading > 1000 criterion bmi_sens$vl_failure.1000 <- vl_failure_90_sens[match(bmi_sens$cfar_pid, vl_failure_90_sens$cfar_pid, nomatch=NA),'vl_failure.1000'] bmi_sens$vl_failure.1000 <- ifelse(is.na(bmi_sens$vl_failure.1000), 0,bmi_sens$vl_failure.1000) bmi_sens$age_vl_failure.1000 <- vl_failure_90_sens[match(bmi_sens$cfar_pid, vl_failure_90_sens$cfar_pid, nomatch=NA),'age_vl_failure.1000'] bmi_sens$age_censor.1000 <- vl_failure_90_sens[match(bmi_sens$cfar_pid, vl_failure_90_sens$cfar_pid, nomatch=NA),'age_censor.1000'] bmi_sens$cum_inc_followup_yrs_alt.1000 <- bmi_sens$age_vl_failure.1000 - bmi_sens$age_fhaart bmi_sens$cum_inc_followup_yrs_alt.1000 <- ifelse(is.na(bmi_sens$cum_inc_followup_yrs_alt.1000),bmi_sens$max_age - bmi_sens$age_fhaart,bmi_sens$cum_inc_followup_yrs_alt.1000) bmi_sens$cum_inc_censor.1000 <- ifelse(bmi_sens$vl_failure.1000 == 1 & bmi_sens$age == bmi_sens$age_vl_failure.1000,1,0) cum_inc_unsuppressed_alt.1000_sens <- cuminc(ftime=bmi_sens$cum_inc_followup_yrs_alt.1000[which(!duplicated(bmi_sens$cfar_pid) & (is.na(bmi_sens$age_censor.1000) | bmi_sens$age_censor.1000 != -99))], fstatus=bmi_sens$vl_failure.1000[which(!duplicated(bmi_sens$cfar_pid) & (is.na(bmi_sens$age_censor.1000) | bmi_sens$age_censor.1000 != -99))], cencode='0',rho=0) which_50_suppressed_sens <- max(which(cum_inc_unsuppressed_alt.1000_sens[[1]]$est <= 0.50)) time_to_50_sens <- cum_inc_unsuppressed_alt.1000_sens[[1]]$time[which_50_suppressed_sens] @ \begin{figure} \caption{Cumulative incidence function for virologic failure using the CCASAnet definition} <>= # pdf('Cumulative_incidence_of_detectable_vl_alternate2.pdf') par(mfrow=c(1,1),mar=c(4,5,1,1),oma=c(1,2,1,1)) plot(cum_inc_unsuppressed_alt.1000[[1]]$est[which(cum_inc_unsuppressed_alt.1000[[1]]$time <= 10.0)] ~ cum_inc_unsuppressed_alt.1000[[1]]$time[which(cum_inc_unsuppressed_alt.1000[[1]]$time <= 10.0)], type='l',axes=FALSE, xlab='Years since ART initiation',ylab='Cumulative incidence of virologic failure\nusing CCASAnet definition') axis(1,at=c(0,2,4,6,8,10)) axis(2) box(bty='l') text(x=7,y=0.02,paste(round(time_to_50,2),' years until 50% of subjects have failed.',sep=''),cex=0.65) @ \end{figure} <>= # Creating cohort in which follow-up is truncated at their time maximum follow-up or at their age of censoring. Also need to remove # those subjects who failed at their first record (age_censor.1000 == -99) which means they would have had to have been censored at # the record prior to their first record. ids_to_remove <- unique(bmi$cfar_pid[which(bmi$age_censor.1000 == -99)]) bmi.1000 <- subset(bmi, cfar_pid %nin% ids_to_remove & (is.na(age_censor.1000) | age <= age_censor.1000)) # Now need to define max age once censored to get total follow-up time bmi.1000$max_age.1000 <- NA split(bmi.1000$max_age.1000, bmi.1000$cfar_pid) <- lapply(split(bmi.1000$age, bmi.1000$cfar_pid), FUN=function(y){ max(y) }) bmi.1000$total_followup_yrs.1000 <- with(bmi.1000, max_age.1000 - age_fhaart) label(bmi.1000$total_followup_yrs.1000) <- 'Total follow-up (yrs)' tbl_90 <- summary(bmi_at_age_fhaart_categorized_larger.factor ~ age_fhaart + sex.factor + race.factor + evidence_of_idu.factor + year_of_first_haart_start + pi_usage.factor + regimen_class_in_fhaart.factor + prior_ade.factor + has_hcv_diagnosis.factor + death.factor + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + bmi_at_age_fhaart + total_followup_yrs.1000 + cohort, data=bmi.1000, method='reverse',overall=TRUE,subset=!duplicated(cfar_pid),test=TRUE) latex(tbl_90, file='', long=TRUE,exclude1=FALSE,longtable=TRUE,caption='Descriptive statistics by categorized BMI at first HAART in cohort followed only until their viral load becomes detectable using the CCASAnet definition of virologic failure', digits=2, size='tiny', landscape=TRUE,prtest='P') # In sensitivity analysis data ids_to_remove_sens <- unique(bmi_sens$cfar_pid[which(bmi_sens$age_censor.1000 == -99)]) bmi.1000_sens <- subset(bmi_sens, cfar_pid %nin% ids_to_remove_sens & (is.na(age_censor.1000) | age <= age_censor.1000)) # Now need to define max age once censored to get total follow-up time bmi.1000_sens$max_age.1000 <- NA split(bmi.1000_sens$max_age.1000, bmi.1000_sens$cfar_pid) <- lapply(split(bmi.1000_sens$age, bmi.1000_sens$cfar_pid), FUN=function(y){ max(y) }) bmi.1000_sens$total_followup_yrs.1000 <- with(bmi.1000_sens, max_age.1000 - age_fhaart) label(bmi.1000_sens$total_followup_yrs.1000) <- 'Total follow-up (yrs)' @ <>= #-----------------------------------------# # Generalized least squares # #-----------------------------------------# # With AR1 correlated data gls_mod_wo_int_ar1.1000 <- Gls(fmla_wo_int, data=bmi.1000, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) gls_mod_w_int_ar1.1000 <- Gls(fmla_w_int, data=bmi.1000, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) # For sensitivity analysis gls_mod_wo_int_ar1.1000_sens <- Gls(fmla_wo_int, data=bmi.1000_sens, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) gls_mod_w_int_ar1.1000_sens <- Gls(fmla_w_int, data=bmi.1000_sens, corr=corAR1(form = ~ 1 | cfar_pid),na.action=na.omit) #-----------------------------------------# # Mixed effects linear model # #-----------------------------------------# # AR1 correlation in lme tmp <- subset(bmi.1000, select=c('cfar_pid','cd4_count_sqrt','bmi_closest','t_lab','sex.factor','race.factor', 'age_fhaart','regimen_class_in_fhaart.factor','year_of_first_haart_start','cd4_count_at_age_fhaart_sqrt', 'loghiv1_rna_at_age_fhaart','cohort')) lme_mod_w_int_ar1.1000 <- lme(fmla_w_int, random=~1|cfar_pid, data=tmp, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) # AR1 correlation in lme for sensitivity analysis tmp_sens <- subset(bmi.1000_sens, select=c('cfar_pid','cd4_count_sqrt','bmi_closest','t_lab','sex.factor','race.factor', 'age_fhaart','regimen_class_in_fhaart.factor','year_of_first_haart_start','cd4_count_at_age_fhaart_sqrt', 'loghiv1_rna_at_age_fhaart','cohort')) lme_mod_w_int_ar1.1000_sens <- lme(fmla_w_int, random=~1|cfar_pid, data=tmp_sens, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) @ <>= # The predicted values does not include those records with missing values so need to mark them to remove them tmp$missing <- apply(tmp, 1, FUN=function(y){ifelse(any(is.na(y)),1,0)}) tmp <- subset(tmp, missing == 0) tmp$residuals <- lme_mod_w_int_ar1.1000[['residuals']][,'fixed'] lme_mod_w_int_ar1_wo_outliers.1000 <- lme(fmla_w_int_wo_outliers, random=~1|cfar_pid, data=tmp, subset=abs(tmp$residuals) <= 6, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) # For sensitivity analysis tmp_sens$missing <- apply(tmp_sens, 1, FUN=function(y){ifelse(any(is.na(y)),1,0)}) tmp_sens <- subset(tmp_sens, missing == 0) tmp_sens$residuals <- lme_mod_w_int_ar1.1000_sens[['residuals']][,'fixed'] lme_mod_w_int_ar1_wo_outliers.1000_sens <- lme(fmla_w_int_wo_outliers, random=~1|cfar_pid, data=tmp_sens, subset=abs(tmp_sens$residuals) <= 6, na.action=na.omit,correlation=corAR1(form=~1|cfar_pid),keep.data=TRUE) @ <>= mod_assessment.1000 <- model_fit_assessment(gls_mod=gls_mod_w_int_ar1.1000, lme_mod=lme_mod_w_int_ar1.1000, data=bmi.1000, fmla=fmla_w_int) cor_df_1000 <- mod_assessment.1000[['cor_df']] dt_1000 <- mod_assessment.1000[['dt']] dr_1000 <- mod_assessment.1000[['dr']] dt_wo_outliers_1000 <- mod_assessment.1000[['dt_wo_outliers']] dr_wo_outliers_1000 <- mod_assessment.1000[['dr_wo_outliers']] lm_residuals_wo_outliers_1000 <- mod_assessment.1000[['lm_residuals_wo_outliers']] # For sensitivity analysis mod_assessment.1000_sens <- model_fit_assessment(gls_mod=gls_mod_w_int_ar1.1000_sens, lme_mod=lme_mod_w_int_ar1.1000_sens, data=bmi.1000_sens, fmla=fmla_w_int) cor_df_1000_sens <- mod_assessment.1000_sens[['cor_df']] dt_1000_sens <- mod_assessment.1000_sens[['dt']] dr_1000_sens <- mod_assessment.1000_sens[['dr']] dt_wo_outliers_1000_sens <- mod_assessment.1000_sens[['dt_wo_outliers']] dr_wo_outliers_1000_sens <- mod_assessment.1000_sens[['dr_wo_outliers']] lm_residuals_wo_outliers_1000_sens <- mod_assessment.1000_sens[['lm_residuals_wo_outliers']] @ \subsubsection{Assessing appropriate correlation structure in the cohort censored using CCASAnet definitions for virologic failure} \normalsize <>= # Table to evaluate the minimum AIC for each of the models from the cohort censored using CCASAnet definitions for virologic failure cor_df_1000 @ \clearpage \begin{figure} \caption{Variogram for assessing correlation structure in the cohort censored using CCASAnet definitions for virologic failure} <>= #pdf('Variogram_wo_points_and_w_lowess_line_capped.pdf') plot(dt_wo_outliers_1000, dr_wo_outliers_1000,main='',xlab='dt',ylab='dr',type='n',ylim=c(0,30)) lines(lowess(x=dt_wo_outliers_1000,y=dr_wo_outliers_1000),col='red') abline(h=var(lm_residuals_wo_outliers_1000),col='blue') #dev.off() @ \end{figure} \begin{figure} \caption{Autocorrelation function graphs for the lme and gls models in the cohort censored using CCASAnet definitions for virologic failure} <>= # Autocorrelation function graphs acf1_1000 <- plot(ACF(lme_mod_w_int_ar1.1000, maxLag=10,resType='n'),alpha=0.01,main='lme model',cex.main=0.65) acf2_1000 <- plot(ACF(gls_mod_w_int_ar1.1000,resType="normalized"),alpha=0.05,main='gls model',cex.main=0.65) print(acf1_1000, position=c(0,0.5,1,1.0),more=TRUE) print(acf2_1000, position=c(0,0,1,0.45)) @ \end{figure} \begin{figure} \caption{QQ plots for all data in the cohort censored using CCASAnet definitions for virologic failure and for data with abs(residuals) $>$ 6 removed.} <>= par(mfrow=c(2,2),mar=c(1,1,1,1),oma=c(1,1,1,1)) # QQ plot #pdf('test_qq.pdf') qq1_1000 <- qqnorm(lme_mod_w_int_ar1.1000, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1),cex.lab=0.65) #dev.off() qq2_1000 <- qqnorm(lme_mod_w_int_ar1_wo_outliers.1000, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1), xlab='Normalized residuals (outliers removed)',cex.lab=0.65) print(qq1_1000, position=c(0,0.5,1,1.0),more=TRUE) print(qq2_1000, position=c(0,0,1,0.45)) @ \end{figure} \clearpage \subsubsection{Comparing model output in cohort censored using CCASAnet definitions for virologic failure} \tiny <>= # ---------------------------------------- # # Comparing the different model outputs: # # ---------------------------------------- # # Gls model using AR1 correlation form=~1|ID #summary(gls_mod_w_int_ar1.1000) gls_coeff_1000 <- gls_mod_w_int_ar1.1000$coefficients # lme model using AR1 correlation form=~1|ID #summary(lme_mod_w_int_ar1.1000) lme_coeff_1000 <- lme_mod_w_int_ar1.1000$coefficients$fixed cbind(format(round(lme_coeff_1000,5),scientific=FALSE),format(round(gls_coeff_1000,5),scientific=FALSE)) @ \clearpage \subsubsection{Evaluating models with interaction terms in the cohort censored using CCASAnet definitions for virologic failure} <>= # ----------------------------- # # Evaluating significance of # # interaction terms # # ----------------------------- # # ANOVA table for gls models anova(gls_mod_w_int_ar1.1000, tol=1e-13) # ANOVA table for lme models anova(lme_mod_w_int_ar1.1000) @ <>= # Bryan wants the individual lines on the graphs to only span the BMIs found for that race/sex group mod_limits_1000 <- graph_limits(bmi.1000) bmi_rge_wm_1000 <- mod_limits_1000[['bmi_rge_wm']] bmi_rge_wf_1000 <- mod_limits_1000[['bmi_rge_wf']] bmi_rge_nwm_1000 <- mod_limits_1000[['bmi_rge_nwm']] bmi_rge_nwf_1000 <- mod_limits_1000[['bmi_rge_nwf']] tmp_rge_wm_1000 <- mod_limits_1000[['tmp_rge_wm']] tmp_rge_wf_1000 <- mod_limits_1000[['tmp_rge_wf']] tmp_rge_nwm_1000 <- mod_limits_1000[['tmp_rge_nwm']] tmp_rge_nwf_1000 <- mod_limits_1000[['tmp_rge_nwf']] tmp_w_residuals_rge_wm_1000 <- mod_limits_1000[['tmp_w_residuals_rge_wm']] tmp_w_residuals_rge_wf_1000 <- mod_limits_1000[['tmp_w_residuals_rge_wf']] tmp_w_residuals_rge_nwm_1000 <- mod_limits_1000[['tmp_w_residuals_rge_nwm']] tmp_w_residuals_rge_nwf_1000 <- mod_limits_1000[['tmp_w_residuals_rge_nwf']] overall_x_rge_1000 <- mod_limits_1000[['overall_x_rge']] @ <>= # From Gls gls_predict_1000 <- Predict(gls_mod_w_int_ar1.1000,bmi_closest,sex.factor,race.factor) # Transforming CD4 back from square root transformed gls_predict_1000$yhat <- gls_predict_1000$yhat^2 gls_predict_1000$lower <- gls_predict_1000$lower^2 gls_predict_1000$upper <- gls_predict_1000$upper^2 # ------------------------------ # # From lme # Getting predicted CD4 count values (square root transformed) cohort_mode_tbl <- sort(table(bmi.1000$cohort),decreasing=TRUE) cohort_mode <- names(cohort_mode_tbl)[1] lme_predict_data_1000 <- data.frame('bmi_closest'=rep(seq(min(bmi.1000$bmi_closest,na.rm=TRUE), max(bmi.1000$bmi_closest,na.rm=TRUE), length.out=100),times=4), 't_lab'=rep(median(bmi.1000$t_lab),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(bmi.1000$age_fhaart[which(!duplicated(bmi.1000$cfar_pid))]), times=400), 'regimen_class_in_fhaart.factor'=rep('NNRTI-based', times=400), # 'regimen_class_in_fhaart.factor'=rep(levels(bmi.1000$regimen_class_in_fhaart.factor)[1], # times=400), 'year_of_first_haart_start'=rep(median(bmi.1000$year_of_first_haart_start[which(!duplicated(bmi.1000$cfar_pid))]), times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(bmi.1000$cd4_count_at_age_fhaart_sqrt[which(!duplicated(bmi.1000$cfar_pid))]), times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(bmi.1000$loghiv1_rna_at_age_fhaart[which(!duplicated(bmi.1000$cfar_pid))], na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode,times=400)) # 'cohort'=rep(levels(bmi.1000$cohort)[1],times=400)) lme_predict_1000 <- predict(lme_mod_w_int_ar1.1000, newdata=lme_predict_data_1000,na.action=na.omit,level=0) ## Transforming CD4 back from square root transformed lme_predict_1000 <- lme_predict_1000^2 @ \begin{figure} \caption{Predicted CD4 count by BMI in the cohort censored using CCASAnet definitions for virologic failure and using lme}\label{fig:cd4_by_bmi_lme_alt_1000} <>= cd4_by_bmi_figure(mod_predict=lme_predict_1000, predict_data=lme_predict_data_1000, fig_limits=mod_limits_1000, plot_method='base', subset_method='overall') @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI in the cohort censored using CCASAnet definitions for virologic failure and using lme and omitting those with abs(residuals) $>$ 6}\label{fig:cd4_by_bmi_lme_alt_wo_outliers_1000} <>= cohort_mode_tbl <- sort(table(tmp$cohort),decreasing=TRUE) cohort_mode <- names(cohort_mode_tbl)[1] lme_predict_data_wo_outliers_1000 <- data.frame('bmi_closest'=rep(seq(min(tmp$bmi_closest[which(abs(tmp$residuals) <= 6)],na.rm=TRUE), max(tmp$bmi_closest[which(abs(tmp$residuals) <= 6)],na.rm=TRUE), length.out=100),times=4), 't_lab'=rep(median(tmp$t_lab),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(tmp$age_fhaart[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)]), times=400), 'regimen_class_in_fhaart.factor'= rep('NNRTI-based', times=400), # 'regimen_class_in_fhaart.factor'= rep(sort(unique(tmp$regimen_class_in_fhaart.factor[which(abs(tmp$residuals) <= 6)]))[1], times=400), 'year_of_first_haart_start'=rep(median(tmp$year_of_first_haart_start[which(abs(tmp$residuals) <= 6)]),times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp$cd4_count_at_age_fhaart_sqrt[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)]),times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp$loghiv1_rna_at_age_fhaart[which(!duplicated(tmp$cfar_pid) & abs(tmp$residuals) <= 6)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode,times=400)) # 'cohort'=rep(sort(unique(tmp$cohort[which(abs(tmp$residuals) <= 6)]))[1],times=400)) lme_predict_wo_outliers_1000 <- predict(lme_mod_w_int_ar1_wo_outliers.1000, newdata=lme_predict_data_wo_outliers_1000,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed lme_predict_wo_outliers_1000 <- lme_predict_wo_outliers_1000^2 cd4_by_bmi_figure(mod_predict=lme_predict_wo_outliers_1000, predict_data=lme_predict_data_wo_outliers_1000, fig_limits=main_mod_limits, plot_method='base', subset_method='wo_outliers') @ \end{figure} <>= cd4_by_year_combinednew1_1000 <- descriptive_graph_mgmt(bmi.1000) @ \begin{figure} \caption{Descriptive graph of CD4 count by year since ART initiation and BMI category. Yearly CD4 counts and BMI measurements were defined as the closest CD4 count within +/- 180 days. A subject could switch BMI categories from year to year. Points at a given year and BMI category are the average of all CD4 counts within that year and BMI category. This graph is on the cohort of subjects who were censored using the CCASAnet definition of virologic failure.} <>= ggplot(cd4_by_year_combinednew1_1000, aes(x=Time,y=CD4_count, group=BMI, colour=BMI)) + geom_line() + xlab('Time since ART initiation (years)') + ylab('CD4 count') + scale_colour_manual(values=c('black','red','green','blue','yellow','purple')) @ \end{figure} \clearpage \begin{figure} \caption{Model-based graph of CD4 count by year since ART initiation and BMI cut points corresponding to common BMI categories. Predicted CD4 values were calculated over 4 years for each of the cut points listed in the legend. This is on the cohort of subjects who were censored using the CCASAnet definition of virologic failure and with outliers removed.} <>= cd4_by_time_model_based_figure(data=tmp, mod=lme_mod_w_int_ar1_wo_outliers.1000, method='base') @ \end{figure} \begin{figure} \caption{Alternate version of above in ggplot2 in cohort censored using CCASAnet definition of virologic failure} <>= ## Alternate version in ggplot2 cd4_by_time_model_based_figure(data=tmp, mod=lme_mod_w_int_ar1_wo_outliers.1000, method='ggplot2') @ \end{figure} \clearpage \section{Sensitivity analysis including women in the analysis who had more than 10\% weight change in a 6 month period.} \normalsize The data contain \Sexpr{nrow(bmi_sens)} records representing \Sexpr{length(unique(bmi_sens$cfar_pid))} subjects. There were \Sexpr{num_cd4_counts_total_sens} non-missing CD4 counts in the data, and \Sexpr{num_bmi_total_sens} non-missing BMI measurements in the data. For a given CD4 count, a BMI measurement was matched to it provided it fell within -180 to +180 days of the CD4 measurement. Therefore, in some instances, a single BMI measurement was matched to multiple CD4 measurements. Once a BMI measurement was matched to a given CD4 count, the distribution of the frequency of both measurements within a subject were calculated. Median (Interquartile Range [IQR]) CD4 counts per person were \Sexpr{desc_cd4_counts_per_person_sens[3]} (\Sexpr{desc_cd4_counts_per_person_sens[2]}, \Sexpr{desc_cd4_counts_per_person_sens[5]}); there were \Sexpr{desc_bmi_per_person_sens[3]} (\Sexpr{desc_bmi_per_person_sens[2]}, \Sexpr{desc_bmi_per_person_sens[5]}) median (IQR) BMI measurements per person. The median (IQR) percentage of records for a given person with BMI measurements matched to multiple CD4 counts is \Sexpr{desc_tk_sens[3]} (\Sexpr{desc_tk_sens[2]},\Sexpr{desc_tk_sens[5]}). \vspace{0.1in} Given the diagnostics run on the subset of data used for coding, a linear mixed effect (lme) model with AR1 correlation was used as the primary model of interest. A generalized least squares (gls) model with AR1 correlation was also fit for completeness. \vspace{0.1in} Table 11 summarizes descriptive statistics on all relevant data in the sensivitiy analysis by categorized BMI values and overall. Assuming the models with all interactions will be retained in the final analysis, Figures \ref{fig:cd4_by_bmi_lme_alt_sens} and \ref{fig:cd4_by_bmi_lme_alt_wo_outliers_sens} illustrate the relationship between time-varying CD4 and BMI using the lme model results stratified by sex and race. Similarly, Figures \ref{fig:cd4_by_bmi_gls_sens} and \ref{fig:cd4_by_bmi_gls_alt_sens} illustrate the relationship between time-varying CD4 and BMI using the gls model results stratified by sex and race. Figure \ref{fig:describe_ggplot_sens} is a descriptive graph for CD4 over time by BMI categories based solely on the data; Figure \ref{fig:model_base_sens} is a more model-based approach to Figure \ref{fig:describe_ggplot_sens}. Analyses on each of the cohorts specified repeat the tables and graphs presented in the overall analysis. <>= tbl1_sens <- summary(bmi_at_age_fhaart_categorized_larger.factor ~ age_fhaart + sex.factor + race.factor + race_for_reviewers.factor + evidence_of_idu.factor + year_of_first_haart_start + year_of_first_haart_start_categorized.factor + pi_usage.factor + regimen_class_in_fhaart.factor + prior_ade.factor + has_hcv_diagnosis.factor + death.factor + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + bmi_at_age_fhaart + cohort, data=bmi_sens, method='reverse',overall=TRUE,subset=!duplicated(cfar_pid),test=TRUE) latex(tbl1_sens, file='', long=TRUE,exclude1=FALSE,longtable=TRUE,caption='Descriptive statistics by categorized BMI at first HAART in data from the sensivitiy analysis.',digits=2,size='tiny',landscape=TRUE,prtest='P',lines.page=40) @ <>= tbl1_supp_sens <- summary(year_of_first_haart_start_categorized.factor ~ sex.factor + race.factor + race_for_reviewers.factor + age_fhaart + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + pi_usage.factor + has_hcv_diagnosis.factor, data=bmi_sens,subset=!duplicated(cfar_pid), method='reverse',overall=TRUE,test=TRUE) tbl1_supp_sens$labels <- c('Sex','Race','Baseline age','Baseline CD4 count','Baseline viral load (log-transformed)','PI usage','Prior ADE','HCV diagnosis') latex(tbl1_supp_sens, file='', exclude1=FALSE, long=TRUE, digits=2, prtest='P',landscape=TRUE,caption='Descriptive statistics on cohort by year of ART initiation in the data from the sensitivity analysis.') @ <>= bmi_sens_female$bmi_change.factor <- factor(bmi_sens_female$bmi_change, levels=c(0,1), labels=c('Original data','Additional data')) tbl1_sens_female <- summary(bmi_change.factor ~ age_fhaart + race.factor + race_for_reviewers.factor + bmi_at_age_fhaart_categorized_larger.factor + evidence_of_idu.factor + year_of_first_haart_start + year_of_first_haart_start_categorized.factor + pi_usage.factor + regimen_class_in_fhaart.factor + prior_ade.factor + has_hcv_diagnosis.factor + death.factor + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + bmi_at_age_fhaart + cohort, data=bmi_sens_female, method='reverse',overall=TRUE,subset=!duplicated(cfar_pid),test=TRUE) latex(tbl1_sens_female, file='', exclude1=FALSE, long=TRUE, digits=2, prtest='P',caption='Descriptive statistics for women by whether they were included in the original analysis or not.',where='!h',size='footnotesize') @ \clearpage \subsection{Assessing appropriate correlation structure} <>= # Table to evaluate the minimum AIC for each of the models. cor_df_sens @ \clearpage \begin{figure} \caption{Variogram for assessing correlation structure in sensitivity models} <>= #pdf('Variogram_wo_points_and_w_lowess_line_capped.pdf') plot(dt_wo_outliers_sens, dr_wo_outliers_sens,main='',xlab='dt',ylab='dr',type='n',ylim=c(0,30)) lines(lowess(x=dt_wo_outliers_sens,y=dr_wo_outliers_sens),col='red') abline(h=var(lm_residuals_wo_outliers_sens),col='blue') #dev.off() @ \end{figure} \begin{figure} \caption{Autocorrelation function graphs for the lme and gls sensitivity models} <>= # Autocorrelation function graphs acf1_sens <- plot(ACF(lme_mod_w_int_ar1_sens, maxLag=10,resType='n'),alpha=0.01,main='lme model',cex.main=0.65) acf2_sens <- plot(ACF(gls_mod_w_int_ar1_sens,resType="normalized"),alpha=0.05,main='gls model',cex.main=0.65) print(acf1_sens, position=c(0,0.5,1,1.0),more=TRUE) print(acf2_sens, position=c(0,0,1,0.45)) @ \end{figure} \begin{figure} \caption{QQ plots for all data and for data in which residuals whose absolute value exceeds 6 are removed in sensitivity models} <>= par(mfrow=c(2,2),mar=c(1,1,1,1),oma=c(1,1,1,1)) # QQ plot #pdf('test_qq.pdf') qq1_sens <- qqnorm(lme_mod_w_int_ar1_sens, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1),cex.lab=0.65) #dev.off() qq2_sens <- qqnorm(lme_mod_w_int_ar1_wo_outliers_sens, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1),xlab='Normalized residuals (outliers removed)',cex.lab=0.65) print(qq1_sens, position=c(0,0.5,1,1.0),more=TRUE) print(qq2_sens, position=c(0,0,1,0.45)) @ \end{figure} \clearpage \subsection{Comparing sensitivity model output} \tiny <>= # ---------------------------------------- # # Comparing the different model outputs: # # ---------------------------------------- # # Gls model using AR1 correlation form=~1|ID #summary(gls_mod_w_int_ar1) gls_coeff_sens <- gls_mod_w_int_ar1_sens$coefficients lme_coeff_sens <- lme_mod_w_int_ar1_sens$coefficients$fixed cbind(format(round(lme_coeff_sens,5),scientific=FALSE),format(round(gls_coeff_sens,5),scientific=FALSE)) @ \clearpage \subsection{Evaluating sensitivity models with interaction terms} <>= # ----------------------------- # # Evaluating significance of # # interaction terms # # ----------------------------- # # ANOVA table for gls models anova(gls_mod_w_int_ar1_sens, tol=1e-13) # ANOVA table for lme models anova(lme_mod_w_int_ar1_sens) @ <>= # Bryan wants the individual lines on the graphs to only span the BMIs found for that race/sex group main_mod_limits_sens <- graph_limits(bmi_sens) bmi_rge_wm_sens <- main_mod_limits_sens[['bmi_rge_wm']] bmi_rge_nwm_sens <- main_mod_limits_sens[['bmi_rge_nwm']] bmi_rge_wf_sens <- main_mod_limits_sens[['bmi_rge_wf']] bmi_rge_nwf_sens <- main_mod_limits_sens[['bmi_rge_nwf']] tmp_rge_wm_sens <- main_mod_limits_sens[['tmp_rge_wm']] tmp_rge_nwm_sens <- main_mod_limits_sens[['tmp_rge_nwm']] tmp_rge_wf_sens <- main_mod_limits_sens[['tmp_rge_wf']] tmp_rge_nwf_sens <- main_mod_limits_sens[['tmp_rge_nwf']] tmp_w_residuals_rge_wm_sens <- main_mod_limits_sens[['tmp_w_residuals_rge_wm']] tmp_w_residuals_rge_nwm_sens <- main_mod_limits_sens[['tmp_w_residuals_rge_nwm']] tmp_w_residuals_rge_wf_sens <- main_mod_limits_sens[['tmp_w_residuals_rge_wf']] tmp_w_residuals_rge_nwf_sens <- main_mod_limits_sens[['tmp_w_residuals_rge_nwf']] overall_x_rge_sens <- main_mod_limits_sens[['overall_x_rge']] @ <>= # From Gls gls_predict_sens <- Predict(gls_mod_w_int_ar1_sens,bmi_closest,sex.factor,race.factor) # Transforming CD4 back from square root transformed gls_predict_sens$yhat <- gls_predict_sens$yhat^2 gls_predict_sens$lower <- gls_predict_sens$lower^2 gls_predict_sens$upper <- gls_predict_sens$upper^2 rge_sens <- range(c(gls_predict_sens$lower,gls_predict_sens$upper)) x.rge_sens <- range(c(gls_predict_sens$bmi_closest)) x.rge_sens <- c(x.rge_sens[1],quantile(gls_predict_sens$bmi_closest,prob=0.90)) x.rge_gls_alt_sens <- x.rge_sens rge_gls_alt_sens <- rge_sens # ------------------------------ # # From lme # Getting predicted CD4 count values (square root transformed) cohort_mode_tbl_sens <- sort(table(bmi_sens$cohort),decreasing=TRUE) cohort_mode_sens <- names(cohort_mode_tbl_sens)[1] lme_predict_data_sens <- data.frame('bmi_closest'=rep(seq(min(bmi_sens$bmi_closest,na.rm=TRUE),max(bmi_sens$bmi_closest,na.rm=TRUE),length.out=100),times=4), 't_lab'=rep(median(bmi_sens$t_lab),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100),rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(bmi_sens$age_fhaart[!duplicated(bmi_sens$cfar_pid)]),times=400), 'regimen_class_in_fhaart.factor'=rep('NNRTI-based',times=400), 'year_of_first_haart_start'=rep(median(bmi_sens$year_of_first_haart_start),times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(bmi_sens$cd4_count_at_age_fhaart_sqrt[!duplicated(bmi_sens$cfar_pid)]),times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(bmi_sens$loghiv1_rna_at_age_fhaart[!duplicated(bmi_sens$cfar_pid)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode_sens,times=400)) lme_predict_sens <- predict(lme_mod_w_int_ar1_sens, newdata=lme_predict_data_sens,na.action=na.omit,level=0) ## Transforming CD4 back from square root transformed lme_predict_sens <- lme_predict_sens^2 @ \begin{figure} \caption{Predicted CD4 count by BMI in sensitivyt model using lme (ggplot version)}\label{fig:cd4_by_bmi_lme_alt_sens} <>= cd4_by_bmi_figure(lme_predict_sens,lme_predict_data_sens, main_mod_limits_sens, plot_method='ggplot2', subset_method='overall') @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI in sensitivity model using lme and omitting those with abs(residuals) $>$ 6 (base graphics version)}\label{fig:cd4_by_bmi_lme_alt_wo_outliers_sens} <>= cohort_mode_tbl_sens <- sort(table(bmi_sens$cohort),decreasing=TRUE) cohort_mode_sens <- names(cohort_mode_tbl_sens)[1] lme_predict_data_wo_outliers_sens <- data.frame('bmi_closest'=rep(seq(min(tmp_sens$bmi_closest[which(abs(tmp_sens$residuals) <= 6)],na.rm=TRUE), max(tmp_sens$bmi_closest[which(abs(tmp_sens$residuals) <= 6)],na.rm=TRUE),length.out=100),times=4), 't_lab'=rep(median(tmp_sens$t_lab[which(abs(tmp_sens$residuals) <= 6)]),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100),rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'= rep(median(tmp_sens$age_fhaart[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)]),times=400), 'regimen_class_in_fhaart.factor'=rep(sort(unique(tmp_sens$regimen_class_in_fhaart.factor[which(abs(tmp_sens$residuals) <= 6)]))[2],times=400), 'year_of_first_haart_start'=rep(median(tmp_sens$year_of_first_haart_start[which(abs(tmp_sens$residuals) <= 6)]),times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp_sens$cd4_count_at_age_fhaart_sqrt[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)]),times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp_sens$loghiv1_rna_at_age_fhaart[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode_sens,times=400)) lme_predict_wo_outliers_sens <- predict(lme_mod_w_int_ar1_wo_outliers_sens, newdata=lme_predict_data_wo_outliers_sens,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed lme_predict_wo_outliers_sens <- lme_predict_wo_outliers_sens^2 cd4_by_bmi_figure(mod_predict=lme_predict_wo_outliers_sens, predict_data=lme_predict_data_wo_outliers_sens, fig_limits=main_mod_limits_sens, plot_method='base', subset_method='wo_outliers') @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI from gls sensitivity model. Predicted values (and their 95\% CIs) for the designated race and sex combination are plotted in purple; the black lines represent the predicted values for the other sex/race combinations for reference.}\label{fig:cd4_by_bmi_gls_sens} <>= ##pdf('test_text.pdf') par(mfrow=c(2,2),oma=c(1,3,1,1),mar=c(4,5,1,1),cex.axis=0.90) rge_sens <- range(c(gls_predict_sens$lower,gls_predict_sens$upper)) x.rge_sens <- range(c(gls_predict_sens$bmi_closest)) x.rge_sens <- c(x.rge_sens[1],quantile(gls_predict_sens$bmi_closest,prob=0.90)) x.rge_sens[2] <- min(c(x.rge_sens[2],50)) # White males plot(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')], type='l', ylim=rge_sens,xlab='',ylab='CD4 count',main='White',xlim=x.rge_sens) polygon(c(gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')], rev(gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')])), c(gls_predict_sens$lower[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')], rev(gls_predict_sens$upper[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')])),col=rgb(0,0,0.255,0.35),border=FALSE) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')],col='blue') lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')],col='black',lty=2) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')],col='black',lty=3) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')],col='black',lty=4) text(x=-10,y=(rge_sens[1]+(rge_sens[2]-rge_sens[1])/2),'Male',xpd=NA,srt=90,cex=1.1,font=2) # Non-white males plot(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')], type='l',ylim=rge_sens,xlab='',ylab='',main='Non-white',xlim=x.rge_sens) polygon(c(gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')], rev(gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')])),c(gls_predict_sens$lower[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')],rev(gls_predict_sens$upper[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')])),col=rgb(0,0,0.255,0.35),border=FALSE) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')],col='blue') lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')],col='black',lty=2) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')],col='black',lty=3) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')],col='black',lty=4) # White females plot(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')], type='l',ylim=rge_sens,xlab='BMI',ylab='CD4 count',main='',xlim=x.rge_sens) polygon(c(gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')], rev(gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')])),c(gls_predict_sens$lower[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')], rev(gls_predict_sens$upper[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')])), col=rgb(0.255,0.20,0.147,0.25),border=FALSE) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')],col=6) # text(x=-1.75,y=(rge[1]+(rge[2]-rge[1])/2),'Female',xpd=NA,srt=90,cex=1.1,font=2) # text(x=-(x.rge[1] + (x.rge[2]-x.rge[1])/2),y=(rge[1]+(rge[2]-rge[1])/2),'Female',xpd=NA,srt=90,cex=1.1,font=2) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')],col='black',lty=2) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')],col='black',lty=3) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')],col='black',lty=4) text(x=-10,y=(rge_sens[1]+(rge_sens[2]-rge_sens[1])/2),'Female',xpd=NA,srt=90,cex=1.1,font=2) # Non-white females plot(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')], type='l',ylim=rge_sens,xlab='BMI',ylab='',main='',xlim=x.rge_sens) polygon(c(gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')], rev(gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')])),c(gls_predict_sens$lower[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')],rev(gls_predict_sens$upper[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')])), col=rgb(0.255,0.20,0.147,0.25),border=FALSE) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')],col=6) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')],col='black',lty=2) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite')],col='black',lty=3) lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white')] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite')],col='black',lty=4) @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI using sensitivity model gls}\label{fig:cd4_by_bmi_gls_alt_sens} <>= # White males plot(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white' & gls_predict_sens$bmi_closest >= bmi_rge_wm_sens[1] & gls_predict_sens$bmi_closest <= min(c(bmi_rge_wm_sens[2],50)))] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'white' & gls_predict_sens$bmi_closest >= bmi_rge_wm_sens[1] & gls_predict_sens$bmi_closest <= min(c(bmi_rge_wm_sens[2],50)))], type='l',ylim=rge,ylab='CD4 count',xlab='BMI',xlim=c(overall_x_rge[1],min(c(overall_x_rge[2],50)))) # Non-white males lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite' & gls_predict_sens$bmi_closest >= bmi_rge_nwm_sens[1] & gls_predict_sens$bmi_closest <= min(c(bmi_rge_nwm_sens[2],50)))] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'male' & gls_predict_sens$race.factor == 'nonwhite' & gls_predict_sens$bmi_closest >= bmi_rge_nwm_sens[1] & gls_predict_sens$bmi_closest <= min(c(bmi_rge_nwm_sens[2],50)))],col='red') # White females lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white' & gls_predict_sens$bmi_closest >= bmi_rge_wf_sens[1] & gls_predict_sens$bmi_closest <= min(c(bmi_rge_wf_sens[2],50)))] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'white' & gls_predict_sens$bmi_closest >= bmi_rge_wf_sens[1] & gls_predict_sens$bmi_closest <= min(c(bmi_rge_wf_sens[2],50)))], col='blue') # Non-white females lines(gls_predict_sens$yhat[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite' & gls_predict_sens$bmi_closest >= bmi_rge_nwf_sens[1] & gls_predict_sens$bmi_closest <= min(c(bmi_rge_nwf_sens[2],50)))] ~ gls_predict_sens$bmi_closest[which(gls_predict_sens$sex.factor == 'female' & gls_predict_sens$race.factor == 'nonwhite' & gls_predict_sens$bmi_closest >= bmi_rge_nwf_sens[1] & gls_predict_sens$bmi_closest <= min(c(bmi_rge_nwf_sens[2],50)))],col='green') legend('topleft',fill=c('black','red','blue','green'),legend=c('White male','Non-white male','White female','Non-white female')) @ \end{figure} <>= cd4_by_year_combinednew1_sens <- descriptive_graph_mgmt(bmi_sens) @ \begin{figure} \caption{Descriptive graph of CD4 count by year since ART initiation and BMI category. Yearly CD4 counts and BMI measurements were defined as the closest CD4 count within +/- 180 days. A subject could switch BMI categories from year to year. Points at a given year and BMI category are the average of all CD4 counts within that year and BMI category. This graph is using data from the sensitivity analysis.}\label{fig:describe_ggplot_sens} <>= ggplot(cd4_by_year_combinednew1_sens, aes(x=Time,y=CD4_count, group=BMI, colour=BMI)) + geom_line() + xlab('Time since ART initiation (years)') + ylab('CD4 count') + scale_colour_manual(values=c('black','red','green','blue','yellow','purple')) @ \end{figure} \begin{figure} \caption{Model-based graph of CD4 counts by year since ART initiation and BMI cut points corresponding to common BMI categories. Predicted CD4 values were calculated over 4 years for each of the cut points listed in the legend. This graph uses data from the sensitivity analysis.}\label{fig:model_base_sens} <>= cd4_by_time_model_based_figure(data=tmp_sens, mod=lme_mod_w_int_ar1_wo_outliers_sens, method='base') @ \end{figure} \clearpage <>= all_comers_predicted_cd4_sens <- predicted_cd4_function(data=tmp_sens, mod=lme_mod_w_int_ar1_wo_outliers_sens) predicted_df_sens <- data.frame('Year'=c(0,1,3,5), '18.5'=round(c(all_comers_predicted_cd4_sens[['year0_18.5']], all_comers_predicted_cd4_sens[['year1_18.5']], all_comers_predicted_cd4_sens[['year3_18.5']], all_comers_predicted_cd4_sens[['year5_18.5']])), '22'=round(c(all_comers_predicted_cd4_sens[['year0_22']], all_comers_predicted_cd4_sens[['year1_22']], all_comers_predicted_cd4_sens[['year3_22']], all_comers_predicted_cd4_sens[['year5_22']])), '25'=round(c(all_comers_predicted_cd4_sens[['year0_25']], all_comers_predicted_cd4_sens[['year1_25']], all_comers_predicted_cd4_sens[['year3_25']], all_comers_predicted_cd4_sens[['year5_25']])), '30'=round(c(all_comers_predicted_cd4_sens[['year0_30']], all_comers_predicted_cd4_sens[['year1_30']], all_comers_predicted_cd4_sens[['year3_30']], all_comers_predicted_cd4_sens[['year5_30']])), '40'=round(c(all_comers_predicted_cd4_sens[['year0_40']], all_comers_predicted_cd4_sens[['year1_40']], all_comers_predicted_cd4_sens[['year3_40']], all_comers_predicted_cd4_sens[['year5_40']])),check.names=FALSE) latex(predicted_df_sens, file='', rowname=NULL, caption='Predicted CD4 values by year and BMI cut-off. These should correspond to the previous graph using data from the sensitivity analysis.',where='!h') @ \begin{figure} \caption{Alternate version of above in ggplot2 using data from the sensitivity analysis}\label{fig:model_ggplot_sens} <>= ## Alternate version in ggplot2 cd4_by_time_model_based_figure(data=tmp_sens, mod=lme_mod_w_int_ar1_wo_outliers_sens, method='ggplot2') @ \end{figure} \clearpage \subsection{Virologically suppressed cohort using original definition in data from the sensitivity analysis.} % ------------------------------------------------------- % % Virologically suppressed % % ------------------------------------------------------- % <>= tbl1_suppressed_sens <- summary(bmi_at_age_fhaart_categorized_larger.factor ~ age_fhaart + sex.factor + race.factor + race_for_reviewers.factor + evidence_of_idu.factor + year_of_first_haart_start + pi_usage.factor + regimen_class_in_fhaart.factor + prior_ade.factor + has_hcv_diagnosis.factor + death.factor + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + bmi_at_age_fhaart + cohort, data=bmi_sens, method='reverse',overall=TRUE,subset=!duplicated(cfar_pid) & vl_suppression_50 == 1, test=TRUE) latex(tbl1_suppressed_sens, file='', long=TRUE,exclude1=FALSE,longtable=TRUE,caption='Descriptive statistics by categorized BMI at first HAART in those virally suppressed $>$ 50 percent of their follow up time in data from the sensitivity analysis.', digits=2, size='tiny', landscape=TRUE,where='!h',prtest='P') @ \subsubsection{Assessing appropriate correlation structure in the virally suppressed cohort in data for the sensitivity analysis.} \normalsize <>= # Table to evaluate the minimum AIC for each of the models from the virally suppressed cohort. cor_df_suppressed_sens @ \clearpage \begin{figure} \caption{Variogram for assessing correlation structure in the virally suppressed cohort in data for the sensitivity analysis.} <>= #pdf('Variogram_wo_points_and_w_lowess_line_capped.pdf') plot(dt_wo_outliers_suppressed_sens, dr_wo_outliers_suppressed_sens,main='',xlab='dt',ylab='dr',type='n',ylim=c(0,30)) lines(lowess(x=dt_wo_outliers_suppressed_sens,y=dr_wo_outliers_suppressed_sens),col='red') abline(h=var(lm_residuals_wo_outliers_suppressed_sens),col='blue') #dev.off() @ \end{figure} \begin{figure} \caption{Autocorrelation function graphs for the lme and gls models in the virally suppressed cohort in data for the sensitivity analysis} <>= # Autocorrelation function graphs acf1_suppressed_sens <- plot(ACF(lme_mod_w_int_suppressed_ar1_sens, maxLag=10,resType='n'),alpha=0.01,main='lme model',cex.main=0.65) acf2_suppressed_sens <- plot(ACF(gls_mod_w_int_suppressed_ar1_sens,resType="normalized"),alpha=0.05,main='gls model',cex.main=0.65) print(acf1_suppressed_sens, position=c(0,0.5,1,1.0),more=TRUE) print(acf2_suppressed_sens, position=c(0,0,1,0.45)) @ \end{figure} \begin{figure} \caption{QQ plots for all data in the virally suppressed cohort and for data with abs(residuals) $>$ 6 removed in data for the sensitivity analysis.} <>= par(mfrow=c(2,2),mar=c(1,1,1,1),oma=c(1,1,1,1)) # QQ plot #pdf('test_qq.pdf') qq1_suppressed_sens <- qqnorm(lme_mod_w_int_suppressed_ar1_sens, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1),cex.lab=0.65) #dev.off() qq2_suppressed_sens <- qqnorm(lme_mod_w_int_suppressed_ar1_wo_outliers_sens, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1), xlab='Normalized residuals (outliers removed)',cex.lab=0.65) print(qq1_suppressed_sens, position=c(0,0.5,1,1.0),more=TRUE) print(qq2_suppressed_sens, position=c(0,0,1,0.45)) @ \end{figure} \clearpage \subsubsection{Comparing model output in virally suppressed cohort in data for the sensitivity analysis.} \tiny <>= # ---------------------------------------- # # Comparing the different model outputs: # # ---------------------------------------- # # Gls model using AR1 correlation form=~1|ID #summary(gls_mod_w_int_suppressed_ar1) gls_coeff_suppressed_sens <- gls_mod_w_int_suppressed_ar1_sens$coefficients lme_coeff_suppressed_sens <- lme_mod_w_int_suppressed_ar1_sens$coefficients$fixed cbind(format(round(lme_coeff_suppressed_sens,5),scientific=FALSE),format(round(gls_coeff_suppressed_sens,5),scientific=FALSE)) @ \clearpage \subsubsection{Evaluating models with interaction terms in virally suppressed cohort in data for the sensitivity analysis.} <>= # ----------------------------- # # Evaluating significance of # # interaction terms # # ----------------------------- # # ANOVA table for gls models anova(gls_mod_w_int_suppressed_ar1_sens, tol=1e-13) # ANOVA table for lme models anova(lme_mod_w_int_suppressed_ar1_sens) @ <>= # Bryan wants the individual lines on the graphs to only span the BMIs found for that race/sex group suppressed_mod_limits_sens <- graph_limits(bmi_sens[which(bmi_sens$vl_suppression_50 == 1),]) bmi_rge_wm_suppressed_sens <- suppressed_mod_limits_sens[['bmi_rge_wm']] bmi_rge_wf_suppressed_sens <- suppressed_mod_limits_sens[['bmi_rge_wf']] bmi_rge_nwm_suppressed_sens <- suppressed_mod_limits_sens[['bmi_rge_nwm']] bmi_rge_nwf_suppressed_sens <- suppressed_mod_limits_sens[['bmi_rge_nwf']] tmp_rge_wm_suppressed_sens <- suppressed_mod_limits_sens[['tmp_rge_wm']] tmp_rge_wf_suppressed_sens <- suppressed_mod_limits_sens[['tmp_rge_wf']] tmp_rge_nwm_suppressed_sens <- suppressed_mod_limits_sens[['tmp_rge_nwm']] tmp_rge_nwf_suppressed_sens <- suppressed_mod_limits_sens[['tmp_rge_nwf']] tmp_w_residuals_rge_wm_suppressed_sens <- suppressed_mod_limits_sens[['tmp_w_residuals_rge_wm']] tmp_w_residuals_rge_wf_suppressed_sens <- suppressed_mod_limits_sens[['tmp_w_residuals_rge_wf']] tmp_w_residuals_rge_nwm_suppressed_sens <- suppressed_mod_limits_sens[['tmp_w_residuals_rge_nwm']] tmp_w_residuals_rge_nwf_suppressed_sens <- suppressed_mod_limits_sens[['tmp_w_residuals_rge_nwf']] overall_x_rge_suppressed_sens <- suppressed_mod_limits_sens[['overall_x_rge']] @ <>= # From Gls gls_predict_suppressed_sens <- Predict(gls_mod_w_int_suppressed_ar1_sens,bmi_closest,sex.factor,race.factor) # Transforming CD4 back from square root transformed gls_predict_suppressed_sens$yhat <- gls_predict_suppressed_sens$yhat^2 gls_predict_suppressed_sens$lower <- gls_predict_suppressed_sens$lower^2 gls_predict_suppressed_sens$upper <- gls_predict_suppressed_sens$upper^2 rge_suppressed_sens <- range(c(gls_predict_suppressed_sens$lower,gls_predict_suppressed_sens$upper)) x.rge_suppressed_sens <- range(c(gls_predict_suppressed_sens$bmi_closest)) x.rge_suppressed_sens <- c(x.rge_suppressed_sens[1],quantile(gls_predict_suppressed_sens$bmi_closest,prob=0.90)) x.rge_gls_alt_suppressed_sens <- x.rge_suppressed_sens rge_gls_alt_suppressed_sens <- rge_suppressed_sens # ------------------------------ # # From lme # Getting predicted CD4 count values (square root transformed) cohort_mode_tbl_sens <- sort(table(bmi_sens$cohort),decreasing=TRUE) cohort_mode_sens <- names(cohort_mode_tbl_sens)[1] lme_predict_data_suppressed_sens <- data.frame('bmi_closest'=rep(seq(min(bmi_sens$bmi_closest[which(bmi_sens$vl_suppression_50 == 1)],na.rm=TRUE), max(bmi_sens$bmi_closest[which(bmi_sens$vl_suppression_50 == 1)],na.rm=TRUE), length.out=100),times=4), 't_lab'=rep(median(bmi_sens$t_lab[which(bmi_sens$vl_suppression_50 == 1)]),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(bmi_sens$age_fhaart[which(!duplicated(bmi_sens$cfar_pid) & bmi_sens$vl_suppression_50 == 1)]), times=400), 'regimen_class_in_fhaart.factor'=rep(levels(bmi_sens$regimen_class_in_fhaart.factor[which(bmi_sens$vl_suppression_50 == 1)])[2], times=400), 'year_of_first_haart_start'=rep(median(bmi_sens$year_of_first_haart_start[which(bmi_sens$vl_suppression_50 == 1)]), times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(bmi_sens$cd4_count_at_age_fhaart_sqrt[which(!duplicated(bmi_sens$cfar_pid) & bmi_sens$vl_suppression_50 == 1)]), times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(bmi_sens$loghiv1_rna_at_age_fhaart[which(!duplicated(bmi_sens$cfar_pid) & bmi_sens$vl_suppression_50 == 1)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode_sens,times=400)) lme_predict_suppressed_sens <- predict(lme_mod_w_int_suppressed_ar1_sens, newdata=lme_predict_data_suppressed_sens,na.action=na.omit,level=0) ## Transforming CD4 back from square root transformed lme_predict_suppressed_sens <- lme_predict_suppressed_sens^2 @ \begin{figure} \caption{Predicted CD4 count by BMI in the virally suppressed cohort using lme in data for sensitivity analysis.}\label{fig:cd4_by_bmi_lme_alt_suppressed_sens} <>= cd4_by_bmi_figure(mod_predict=lme_predict_suppressed_sens, predict_data=lme_predict_data_suppressed_sens, fig_limits= suppressed_mod_limits_sens, plot_method='base',subset_method='lme') @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI in the virally suppressed cohort using lme and omitting those with abs(residuals) $>$ 6 in data for sensitivity analysis.}\label{fig:cd4_by_bmi_lme_alt_wo_outliers_suppressed_sens} <>= # Getting predicted CD4 count values (square root transformed) cohort_mode_tbl_sens <- sort(table(tmp_sens$cohort),decreasing=TRUE) cohort_mode_sens <- names(cohort_mode_tbl_sens)[1] lme_predict_data_wo_outliers_suppressed_sens <- data.frame('bmi_closest'=rep(seq(min(tmp_sens$bmi_closest[which(abs(tmp_sens$residuals) <= 6)],na.rm=TRUE), max(tmp_sens$bmi_closest[which(abs(tmp_sens$residuals) <= 6)],na.rm=TRUE),length.out=100),times=4), 't_lab'=rep(median(tmp_sens$t_lab),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(tmp_sens$age_fhaart[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)]), times=400), 'regimen_class_in_fhaart.factor'= rep('NNRTI-based', times=400), 'year_of_first_haart_start'=rep(median(tmp_sens$year_of_first_haart_start[which(abs(tmp_sens$residuals) <= 6)]),times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp_sens$cd4_count_at_age_fhaart_sqrt[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)]),times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp_sens$loghiv1_rna_at_age_fhaart[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode_sens,times=400)) lme_predict_wo_outliers_suppressed_sens <- predict(lme_mod_w_int_suppressed_ar1_wo_outliers_sens, newdata=lme_predict_data_wo_outliers_suppressed_sens,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed lme_predict_wo_outliers_suppressed_sens <- lme_predict_wo_outliers_suppressed_sens^2 cd4_by_bmi_figure(mod_predict=lme_predict_wo_outliers_suppressed_sens, predict_data=lme_predict_data_wo_outliers_suppressed_sens, fig_limits=main_mod_limits_sens, plot_method='base', subset_method='wo_outliers') @ \end{figure} <>= cd4_by_year_combinednew1_suppressed_sens <- descriptive_graph_mgmt(subset(bmi_sens, vl_suppression_50 == 1)) @ \begin{figure} \caption{Descriptive graph of CD4 count by year since ART initiation and BMI category. Yearly CD4 counts and BMI measurements were defined as the closest CD4 count within +/- 180 days. A subject could switch BMI categories from year to year. Points at a given year and BMI category are the average of all CD4 counts within that year and BMI category. This graph is on the cohort of subjects who were virally suppressed at least 50\% of the time and using data for the sensitivity analysis.} <>= ggplot(cd4_by_year_combinednew1_suppressed_sens, aes(x=Time,y=CD4_count, group=BMI, colour=BMI)) + geom_line() + xlab('Time since ART initiation (years)') + ylab('CD4 count') + scale_colour_manual(values=c('black','red','green','blue','yellow','purple')) @ \end{figure} \clearpage \begin{figure} \caption{Model-based graph of CD4 count by year since ART initiation and by BMI cut points corresponding to common BMI categories. Predicted CD4 values were calculated over 4 years for each of the cut points listed in the legend. This is on the cohort of subjects who were virally suppressed at least 50\% of the time and with outliers removed. This graph uses data from the sensitivity analysis.} <>= cd4_by_time_model_based_figure(data=tmp_sens, mod=lme_mod_w_int_suppressed_ar1_wo_outliers_sens, method='base') @ \end{figure} \begin{figure} \caption{Alternate version of above in ggplot2 in those suppressed $>$ 50 percent of the time. This graph uses data from the sensitivity analysis.} <>= ## Alternate version in ggplot2 cd4_by_time_model_based_figure(data=tmp_sens, mod=lme_mod_w_int_suppressed_ar1_wo_outliers_sens, method='ggplot2') @ \end{figure} \subsection{An alternate method for determining those who were virologically suppressed more than 50 \% of the time using data for the sensitivity analysis.} \normalsize Results are presented in this section in which the cohort of those who were virologically suppressed more than 50\% of the time was determined by simply counting the number of records in which a suppressed viral load measurement occurred and comparing it to the number of records in which a detectable viral load measurement occurred. Those who had more measurements in which viral load was suppressed were considered suppressed more than 50\% of their follow-up time. <>= tbl1_suppressed_alt_sens <- summary(bmi_at_age_fhaart_categorized_larger.factor ~ age_fhaart + sex.factor + race.factor + race_for_reviewers.factor + evidence_of_idu.factor + year_of_first_haart_start + pi_usage.factor + regimen_class_in_fhaart.factor + prior_ade.factor + has_hcv_diagnosis.factor + death.factor + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + bmi_at_age_fhaart + cohort, data=bmi_sens, method='reverse',overall=TRUE,subset=!duplicated(cfar_pid) & vl_suppression_50_alt == 1, test=TRUE) latex(tbl1_suppressed_alt_sens, file='', long=TRUE,exclude1=FALSE,longtable=TRUE,caption='Descriptive statistics by categorized BMI at first HAART in those virally suppressed $>$ 50 percent of their follow up time using the alternate definition and data for the sensitivity analysis.', digits=2, size='tiny', landscape=TRUE,where='!h',prtest='P') @ \clearpage \subsubsection{Assessing appropriate correlation structure in the virally suppressed cohort using the alternate definition for viral suppression for the sensitivity analysis.} <>= # Table to evaluate the minimum AIC for each of the models from the virally suppressed cohort. cor_df_suppressed_alt_sens @ \clearpage \begin{figure} \caption{Variogram for assessing correlation structure in the virally suppressed cohort using the alternate definition for viral suppression for the sensitivity analysis.} <>= #pdf('Variogram_wo_points_and_w_lowess_line_capped.pdf') plot(dt_suppressed_alt_sens, dr_suppressed_alt_sens,main='',xlab='dt',ylab='dr',type='n',ylim=c(0,30)) lines(lowess(x=dt_suppressed_alt_sens,y=dr_suppressed_alt_sens),col='red') abline(h=var(lm_residuals_wo_outliers_suppressed_alt_sens),col='blue') #dev.off() @ \end{figure} \begin{figure} \caption{Autocorrelation function graphs for the lme and gls models in the virally suppressed cohort using the alternate definition for viral suppression for the sensitivity analysis.} <>= # Autocorrelation function graphs acf1_suppressed_alt_sens <- plot(ACF(lme_mod_w_int_suppressed_ar1_alt_sens, maxLag=10,resType='n'),alpha=0.01,main='lme model',cex.main=0.65) acf2_suppressed_alt_sens <- plot(ACF(gls_mod_w_int_suppressed_ar1_alt_sens,resType="normalized"),alpha=0.05,main='gls model',cex.main=0.65) print(acf1_suppressed_alt_sens, position=c(0,0.5,1,1.0),more=TRUE) print(acf2_suppressed_alt_sens, position=c(0,0,1,0.45)) @ \end{figure} \begin{figure} \caption{QQ plots for all data and when outliers were removed in the virally suppressed cohort using the alternate definition for viral suppression for the sensitivity analysis.} <>= par(mfrow=c(2,2),mar=c(1,1,1,1),oma=c(1,1,1,1)) # QQ plot #pdf('test_qq.pdf') qq1_suppressed_alt_sens <- qqnorm(lme_mod_w_int_suppressed_ar1_alt_sens, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1),cex.lab=0.65) #dev.off() qq2_suppressed_alt_sens <- qqnorm(lme_mod_w_int_suppressed_ar1_wo_outliers_alt_sens, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1), xlab='Normalized residuals (outliers removed)',cex.lab=0.65) print(qq1_suppressed_alt_sens, position=c(0,0.5,1,1.0),more=TRUE) print(qq2_suppressed_alt_sens, position=c(0,0,1,0.45)) @ \end{figure} \clearpage \subsubsection{Comparing model output in virally suppressed cohort using the alternate definition for viral suppression for the sensitivity analysis.} \tiny <>= # ---------------------------------------- # # Comparing the different model outputs: # # ---------------------------------------- # # Gls model using AR1 correlation form=~1|ID #summary(gls_mod_w_int_suppressed_ar1) gls_coeff_suppressed_alt_sens <- gls_mod_w_int_suppressed_ar1_alt_sens$coefficients # lme model using AR1 correlation form=~1|ID #summary(lme_mod_w_int_suppressed_ar1) lme_coeff_suppressed_alt_sens <- lme_mod_w_int_suppressed_ar1_alt_sens$coefficients$fixed cbind(format(round(lme_coeff_suppressed_alt_sens,5),scientific=FALSE),format(round(gls_coeff_suppressed_alt_sens,5),scientific=FALSE)) @ \clearpage \subsubsection{Evaluating models with interaction terms in virally suppressed cohort using the alternate definition for viral suppression for the sensitivity analysis.} <>= # ----------------------------- # # Evaluating significance of # # interaction terms # # ----------------------------- # # ANOVA table for gls models anova(gls_mod_w_int_suppressed_ar1_alt_sens, tol=1e-13) # ANOVA table for lme models anova(lme_mod_w_int_suppressed_ar1_alt_sens) @ <>= # Bryan wants the individual lines on the graphs to only span the BMIs found for that race/sex group suppressed_mod_limits_alt_sens <- graph_limits(bmi_sens[which(bmi_sens$vl_suppression_50_alt == 1),]) bmi_rge_wm_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['bmi_rge_wm']] bmi_rge_wf_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['bmi_rge_wf']] bmi_rge_nwm_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['bmi_rge_nwm']] bmi_rge_nwf_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['bmi_rge_nwf']] tmp_rge_wm_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['tmp_rge_wm']] tmp_rge_wf_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['tmp_rge_wf']] tmp_rge_nwm_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['tmp_rge_nwm']] tmp_rge_nwf_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['tmp_rge_nwf']] tmp_w_residuals_rge_wm_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['tmp_w_residuals_rge_wm']] tmp_w_residuals_rge_wf_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['tmp_w_residuals_rge_wf']] tmp_w_residuals_rge_nwm_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['tmp_w_residuals_rge_nwm']] tmp_w_residuals_rge_nwf_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['tmp_w_residuals_rge_nwf']] overall_x_rge_suppressed_alt_sens <- suppressed_mod_limits_alt_sens[['overall_x_rge']] @ <>= # From Gls gls_predict_suppressed_alt_sens <- Predict(gls_mod_w_int_suppressed_ar1_alt_sens,bmi_closest,sex.factor,race.factor) # Transforming CD4 back from square root transformed gls_predict_suppressed_alt_sens$yhat <- gls_predict_suppressed_alt_sens$yhat^2 gls_predict_suppressed_alt_sens$lower <- gls_predict_suppressed_alt_sens$lower^2 gls_predict_suppressed_alt_sens$upper <- gls_predict_suppressed_alt_sens$upper^2 # ------------------------------ # # From lme # Getting predicted CD4 count values (square root transformed) cohort_mode_tbl_sens <- sort(table(tmp_sens$cohort),decreasing=TRUE) cohort_mode_sens <- names(cohort_mode_tbl_sens)[1] lme_predict_data_suppressed_alt_sens <- data.frame('bmi_closest'=rep(seq(min(tmp_sens$bmi_closest,na.rm=TRUE), max(tmp_sens$bmi_closest,na.rm=TRUE), length.out=100),times=4), 't_lab'=rep(median(tmp_sens$t_lab),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(tmp_sens$age_fhaart[!duplicated(tmp_sens$cfar_pid)]), times=400), 'regimen_class_in_fhaart.factor'=rep('NNRTI-based', times=400), 'year_of_first_haart_start'=rep(median(tmp_sens$year_of_first_haart_start), times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp_sens$cd4_count_at_age_fhaart_sqrt[!duplicated(tmp_sens$cfar_pid)]), times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp_sens$loghiv1_rna_at_age_fhaart[!duplicated(tmp_sens$cfar_pid)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode_sens,times=400)) lme_predict_suppressed_alt_sens <- predict(lme_mod_w_int_suppressed_ar1_alt_sens, newdata=lme_predict_data_suppressed_alt_sens,na.action=na.omit,level=0) ## Transforming CD4 back from square root transformed lme_predict_suppressed_alt_sens <- lme_predict_suppressed_alt_sens^2 @ \begin{figure} \caption{Predicted CD4 count by BMI in the virally suppressed cohort using the alternate suppression definition and using lme}\label{fig:cd4_by_bmi_lme_alt_suppressed_alt using data for the sensitivity analysis.} <>= cd4_by_bmi_figure(mod_predict=lme_predict_suppressed_alt_sens, predict_data=lme_predict_data_suppressed_alt_sens, fig_limits= suppressed_mod_limits_alt_sens, plot_method='base', subset_method='lme') @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI in the virally suppressed cohort using the alternate definition for suppression and using lme and omitting those with abs(residuals) $>$ 6. This graph uses data for the sensitivity analysis}\label{fig:cd4_by_bmi_lme_alt_wo_outliers_suppressed_alt_sens} <>= cohort_mode_tbl_sens <- sort(table(tmp_sens$cohort),decreasing=TRUE) cohort_mode_sens <- names(cohort_mode_tbl_sens)[1] lme_predict_data_wo_outliers_suppressed_alt_sens <- data.frame('bmi_closest'=rep(seq(min(tmp_sens$bmi_closest[which(abs(tmp_sens$residuals) <= 6)],na.rm=TRUE), max(tmp_sens$bmi_closest[which(abs(tmp_sens$residuals) <= 6)],na.rm=TRUE),length.out=100),times=4), 't_lab'=rep(median(tmp_sens$t_lab[which(abs(tmp_sens$residuals) <= 6)]),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(tmp_sens$age_fhaart[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)]), times=400), 'regimen_class_in_fhaart.factor'=rep('NNRTI-based', times=400), 'year_of_first_haart_start'=rep(median(tmp_sens$year_of_first_haart_start[which(abs(tmp_sens$residuals) <= 6)]),times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp_sens$cd4_count_at_age_fhaart_sqrt[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)]),times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp_sens$loghiv1_rna_at_age_fhaart[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode_sens,times=400)) lme_predict_wo_outliers_suppressed_alt_sens <- predict(lme_mod_w_int_suppressed_ar1_wo_outliers_alt_sens, newdata=lme_predict_data_wo_outliers_suppressed_alt_sens,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed lme_predict_wo_outliers_suppressed_alt_sens <- lme_predict_wo_outliers_suppressed_alt_sens^2 cd4_by_bmi_figure(mod_predict=lme_predict_wo_outliers_suppressed_alt_sens, predict_data=lme_predict_data_wo_outliers_suppressed_alt_sens, fig_limits=main_mod_limits_sens, plot_method='base', subset_method='wo_outliers') @ \end{figure} <>= cd4_by_year_combinednew1_suppressed_alt_sens <- descriptive_graph_mgmt(subset(bmi_sens, vl_suppression_50_alt == 1)) @ \begin{figure} \caption{Descriptive graph of CD4 count by year since ART initiation and BMI category. Yearly CD4 counts and BMI measurements were defined as the closest CD4 count within +/- 180 days. A subject could switch BMI categories from year to year. Points at a given year and BMI category are the average of all CD4 counts within that year and BMI category. This graph is on the cohort of subjects who were virally suppressed at least 50\% of the time using the alternate definition for viral suppression and data for the sensitivity analysis.} <>= ggplot(cd4_by_year_combinednew1_suppressed_alt_sens, aes(x=Time,y=CD4_count, group=BMI, colour=BMI)) + geom_line() + xlab('Time since ART initiation (years)') + ylab('CD4 count') + scale_colour_manual(values=c('black','red','green','blue','yellow','purple')) @ \end{figure} \clearpage \begin{figure} \caption{Model-based graph of CD4 count by year since ART initiation and BMI cut points corresponding to common BMI categories. Predicted CD4 values were calculated over 4 years for each of the cut points listed in the legend. This is on the cohort of subjects who were virally suppressed at least 50\% of the time using the alternate definition for viral suppression and outliers have been removed and using data for the sensitivity analysis.} <>= cd4_by_time_model_based_figure(data=tmp_sens, mod=lme_mod_w_int_suppressed_ar1_wo_outliers_alt_sens, method='base') @ \end{figure} \begin{figure} \caption{Alternate version of above in ggplot2 in those suppressed $>$ 50 percent of the time using the alternate definition for viral suppression and outliers removed and using data for the sensitivity analysis.} <>= ## Alternate version in ggplot2 cd4_by_time_model_based_figure(data=tmp_sens, mod=lme_mod_w_int_suppressed_ar1_wo_outliers_alt_sens, method='ggplot2') @ \end{figure} \clearpage \subsubsection{Censoring follow-up when a subject becomes unsuppressed in data for sensitivity analysis.} \normalsize This analysis censors subjects who meet the CCASAnet definition for virologic failure. Subjects are censored at the first record having a non-missing viral load measurement prior to the failure. Recall that the age of failure was the age at the 2nd of two consecutive detectable viral loads or the age at which one measurement exceeded 1000. \begin{figure} \caption{Cumulative incidence function for virologic failure using the CCASAnet definition in the data for the sensitivity analysis.} <>= # pdf('Cumulative_incidence_of_detectable_vl_alternate2.pdf') par(mfrow=c(1,1),mar=c(4,5,1,1),oma=c(1,2,1,1)) plot(cum_inc_unsuppressed_alt.1000_sens[[1]]$est[which(cum_inc_unsuppressed_alt.1000_sens[[1]]$time <= 10.0)] ~ cum_inc_unsuppressed_alt.1000_sens[[1]]$time[which(cum_inc_unsuppressed_alt.1000_sens[[1]]$time <= 10.0)], type='l',axes=FALSE, xlab='Years since ART initiation',ylab='Cumulative incidence of virologic failure\nusing CCASAnet definition') axis(1,at=c(0,2,4,6,8,10)) axis(2) box(bty='l') text(x=7,y=0.02,paste(round(time_to_50_sens,2),' years until 50% of subjects have failed.',sep=''),cex=0.65) @ \end{figure} <>= tbl_90_sens <- summary(bmi_at_age_fhaart_categorized_larger.factor ~ age_fhaart + sex.factor + race.factor + evidence_of_idu.factor + year_of_first_haart_start + pi_usage.factor + regimen_class_in_fhaart.factor + prior_ade.factor + has_hcv_diagnosis.factor + death.factor + cd4_count_at_age_fhaart + loghiv1_rna_at_age_fhaart + bmi_at_age_fhaart + total_followup_yrs.1000 + cohort, data=bmi.1000_sens, method='reverse',overall=TRUE,subset=!duplicated(cfar_pid),test=TRUE) latex(tbl_90_sens, file='', long=TRUE,exclude1=FALSE,longtable=TRUE,caption='Descriptive statistics by categorized BMI at first HAART in cohort followed only until their viral load becomes detectable using the CCASAnet definition of virologic failure and data for sensitivity analysis.', digits=2, size='tiny', landscape=TRUE,prtest='P') @ \subsubsection{Assessing appropriate correlation structure in the cohort censored using CCASAnet definitions for virologic failure using data for the sensitivity analysis} \normalsize <>= # Table to evaluate the minimum AIC for each of the models from the cohort censored using CCASAnet definitions for virologic failure cor_df_1000_sens @ \clearpage \begin{figure} \caption{Variogram for assessing correlation structure in the cohort censored using CCASAnet definitions for virologic failure using data for the sensitivity analysis.} <>= #pdf('Variogram_wo_points_and_w_lowess_line_capped.pdf') plot(dt_wo_outliers_1000_sens, dr_wo_outliers_1000_sens,main='',xlab='dt',ylab='dr',type='n',ylim=c(0,30)) lines(lowess(x=dt_wo_outliers_1000_sens,y=dr_wo_outliers_1000_sens),col='red') abline(h=var(lm_residuals_wo_outliers_1000_sens),col='blue') #dev.off() @ \end{figure} \begin{figure} \caption{Autocorrelation function graphs for the lme and gls models in the cohort censored using CCASAnet definitions for virologic failure using data for the sensitivity analysis.} <>= # Autocorrelation function graphs acf1_1000_sens <- plot(ACF(lme_mod_w_int_ar1.1000_sens, maxLag=10,resType='n'),alpha=0.01,main='lme model',cex.main=0.65) acf2_1000_sens <- plot(ACF(gls_mod_w_int_ar1.1000_sens,resType="normalized"),alpha=0.05,main='gls model',cex.main=0.65) print(acf1_1000_sens, position=c(0,0.5,1,1.0),more=TRUE) print(acf2_1000_sens, position=c(0,0,1,0.45)) @ \end{figure} \begin{figure} \caption{QQ plots for all data in the cohort censored using CCASAnet definitions for virologic failure and for data with abs(residuals) $>$ 6 removed using data for the sensitivity analysis.} <>= par(mfrow=c(2,2),mar=c(1,1,1,1),oma=c(1,1,1,1)) # QQ plot #pdf('test_qq.pdf') qq1_1000_sens <- qqnorm(lme_mod_w_int_ar1.1000_sens, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1),cex.lab=0.65) #dev.off() qq2_1000_sens <- qqnorm(lme_mod_w_int_ar1_wo_outliers.1000_sens, ~resid(., type='n'),main='',cex.main=0.65,pch=18,abline=c(0,1), xlab='Normalized residuals (outliers removed)',cex.lab=0.65) print(qq1_1000_sens, position=c(0,0.5,1,1.0),more=TRUE) print(qq2_1000_sens, position=c(0,0,1,0.45)) @ \end{figure} \clearpage \subsubsection{Comparing model output in cohort censored using CCASAnet definitions for virologic failure using data for the sensitivity analysis.} \tiny <>= # ---------------------------------------- # # Comparing the different model outputs: # # ---------------------------------------- # # Gls model using AR1 correlation form=~1|ID #summary(gls_mod_w_int_ar1.1000) gls_coeff_1000_sens <- gls_mod_w_int_ar1.1000_sens$coefficients lme_coeff_1000_sens <- lme_mod_w_int_ar1.1000_sens$coefficients$fixed cbind(format(round(lme_coeff_1000_sens,5),scientific=FALSE),format(round(gls_coeff_1000_sens,5),scientific=FALSE)) @ \clearpage \subsubsection{Evaluating models with interaction terms in the cohort censored using CCASAnet definitions for virologic failure using data for the sensitivity analysis.} <>= # ----------------------------- # # Evaluating significance of # # interaction terms # # ----------------------------- # # ANOVA table for gls models anova(gls_mod_w_int_ar1.1000_sens, tol=1e-13) # ANOVA table for lme models anova(lme_mod_w_int_ar1.1000_sens) @ <>= # Bryan wants the individual lines on the graphs to only span the BMIs found for that race/sex group mod_limits_1000_sens <- graph_limits(bmi.1000_sens) bmi_rge_wm_1000_sens <- mod_limits_1000_sens[['bmi_rge_wm']] bmi_rge_wf_1000_sens <- mod_limits_1000_sens[['bmi_rge_wf']] bmi_rge_nwm_1000_sens <- mod_limits_1000_sens[['bmi_rge_nwm']] bmi_rge_nwf_1000_sens <- mod_limits_1000_sens[['bmi_rge_nwf']] tmp_rge_wm_1000_sens <- mod_limits_1000_sens[['tmp_rge_wm']] tmp_rge_wf_1000_sens <- mod_limits_1000_sens[['tmp_rge_wf']] tmp_rge_nwm_1000_sens <- mod_limits_1000_sens[['tmp_rge_nwm']] tmp_rge_nwf_1000_sens <- mod_limits_1000_sens[['tmp_rge_nwf']] tmp_w_residuals_rge_wm_1000_sens <- mod_limits_1000_sens[['tmp_w_residuals_rge_wm']] tmp_w_residuals_rge_wf_1000_sens <- mod_limits_1000_sens[['tmp_w_residuals_rge_wf']] tmp_w_residuals_rge_nwm_1000_sens <- mod_limits_1000_sens[['tmp_w_residuals_rge_nwm']] tmp_w_residuals_rge_nwf_1000_sens <- mod_limits_1000_sens[['tmp_w_residuals_rge_nwf']] overall_x_rge_1000_sens <- mod_limits_1000_sens[['overall_x_rge']] @ <>= # From Gls gls_predict_1000_sens <- Predict(gls_mod_w_int_ar1.1000_sens,bmi_closest,sex.factor,race.factor) # Transforming CD4 back from square root transformed gls_predict_1000_sens$yhat <- gls_predict_1000_sens$yhat^2 gls_predict_1000_sens$lower <- gls_predict_1000_sens$lower^2 gls_predict_1000_sens$upper <- gls_predict_1000_sens$upper^2 # ------------------------------ # # From lme # Getting predicted CD4 count values (square root transformed) cohort_mode_tbl_sens <- sort(table(bmi.1000_sens$cohort),decreasing=TRUE) cohort_mode_sens <- names(cohort_mode_tbl_sens)[1] lme_predict_data_1000_sens <- data.frame('bmi_closest'=rep(seq(min(bmi.1000_sens$bmi_closest,na.rm=TRUE), max(bmi.1000_sens$bmi_closest,na.rm=TRUE), length.out=100),times=4), 't_lab'=rep(median(bmi.1000_sens$t_lab),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(bmi.1000_sens$age_fhaart[which(!duplicated(bmi.1000_sens$cfar_pid))]), times=400), 'regimen_class_in_fhaart.factor'=rep('NNRTI-based', times=400), 'year_of_first_haart_start'=rep(median(bmi.1000_sens$year_of_first_haart_start[which(!duplicated(bmi.1000_sens$cfar_pid))]), times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(bmi.1000_sens$cd4_count_at_age_fhaart_sqrt[which(!duplicated(bmi.1000_sens$cfar_pid))]), times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(bmi.1000_sens$loghiv1_rna_at_age_fhaart[which(!duplicated(bmi.1000_sens$cfar_pid))],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode_sens,times=400)) lme_predict_1000_sens <- predict(lme_mod_w_int_ar1.1000_sens, newdata=lme_predict_data_1000_sens,na.action=na.omit,level=0) ## Transforming CD4 back from square root transformed lme_predict_1000_sens <- lme_predict_1000_sens^2 @ \begin{figure} \caption{Predicted CD4 count by BMI in the cohort censored using CCASAnet definitions for virologic failure and using lme and data for sensitivity analysis.}\label{fig:cd4_by_bmi_lme_alt_1000_sens} <>= cd4_by_bmi_figure(mod_predict=lme_predict_1000_sens, predict_data=lme_predict_data_1000_sens, fig_limits=mod_limits_1000_sens, plot_method='base', subset_method='overall') @ \end{figure} \begin{figure} \caption{Predicted CD4 count by BMI in the cohort censored using CCASAnet definitions for virologic failure and using lme and omitting those with abs(residuals) $>$ 6 and using data for sensitivity analysis.}\label{fig:cd4_by_bmi_lme_alt_wo_outliers_1000_sens} <>= cohort_mode_tbl_sens <- sort(table(tmp_sens$cohort),decreasing=TRUE) cohort_mode_sens <- names(cohort_mode_tbl_sens)[1] lme_predict_data_wo_outliers_1000_sens <- data.frame('bmi_closest'=rep(seq(min(tmp_sens$bmi_closest[which(abs(tmp_sens$residuals) <= 6)],na.rm=TRUE), max(tmp_sens$bmi_closest[which(abs(tmp_sens$residuals) <= 6)],na.rm=TRUE),length.out=100),times=4), 't_lab'=rep(median(tmp_sens$t_lab),times=400), 'sex.factor'=c(rep('female',times=200),rep('male',times=200)), 'race.factor'=c(rep('white',times=100),rep('nonwhite',times=100), rep('white',times=100),rep('nonwhite',times=100)), 'age_fhaart'=rep(median(tmp_sens$age_fhaart[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)]), times=400), 'regimen_class_in_fhaart.factor'= rep('NNRTI-based', times=400), 'year_of_first_haart_start'=rep(median(tmp_sens$year_of_first_haart_start[which(abs(tmp_sens$residuals) <= 6)]),times=400), 'cd4_count_at_age_fhaart_sqrt'=rep(median(tmp_sens$cd4_count_at_age_fhaart_sqrt[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)]),times=400), 'loghiv1_rna_at_age_fhaart'=rep(median(tmp_sens$loghiv1_rna_at_age_fhaart[which(!duplicated(tmp_sens$cfar_pid) & abs(tmp_sens$residuals) <= 6)],na.rm=TRUE),times=400), 'cohort'=rep(cohort_mode_sens,times=400)) lme_predict_wo_outliers_1000_sens <- predict(lme_mod_w_int_ar1_wo_outliers.1000_sens, newdata=lme_predict_data_wo_outliers_1000_sens,na.action=na.omit,level=0) # Transforming CD4 back from square root transformed lme_predict_wo_outliers_1000_sens <- lme_predict_wo_outliers_1000_sens^2 cd4_by_bmi_figure(mod_predict=lme_predict_wo_outliers_1000_sens, predict_data=lme_predict_data_wo_outliers_1000_sens, fig_limits=main_mod_limits_sens, plot_method='base', subset_method='wo_outliers') @ \end{figure} <>= cd4_by_year_combinednew1_1000_sens <- descriptive_graph_mgmt(bmi.1000_sens) @ \begin{figure} \caption{Descriptive graph of CD4 count by year since ART initiation and BMI category. Yearly CD4 counts and BMI measurements were defined as the closest CD4 count within +/- 180 days. A subject could switch BMI categories from year to year. Points at a given year and BMI category are the average of all CD4 counts within that year and BMI category. This graph is on the cohort of subjects who were censored using the CCASAnet definition of virologic failure and uses data for the sensitivity analysis.} <>= ggplot(cd4_by_year_combinednew1_1000_sens, aes(x=Time,y=CD4_count, group=BMI, colour=BMI)) + geom_line() + xlab('Time since ART initiation (years)') + ylab('CD4 count') + scale_colour_manual(values=c('black','red','green','blue','yellow','purple')) @ \end{figure} \clearpage \begin{figure} \caption{Model-based graph of CD4 count by year since ART initiation and BMI cut points corresponding to common BMI categories. Predicted CD4 values were calculated over 4 years for each of the cut points listed in the legend. This is on the cohort of subjects who were censored using the CCASAnet definition of virologic failure and with outliers removed and uses data for the sensitivity analysis.} <>= cd4_by_time_model_based_figure(data=tmp_sens, mod=lme_mod_w_int_ar1_wo_outliers.1000_sens, method='base') @ \end{figure} \begin{figure} \caption{Alternate version of above in ggplot2 in cohort censored using CCASAnet definition of virologic failure and uses data for the sensitivity analysis.} <>= ## Alternate version in ggplot2 cd4_by_time_model_based_figure(data=tmp_sens, mod=lme_mod_w_int_ar1_wo_outliers.1000_sens, method='ggplot2') @ \end{figure} <>= ##Run your model #model <- lme(weight ~ Time + Diet, BodyWeight, ~ 1 | Rat) #summary(model) ##Predict the values ##predict.lme is a pain because you have to specify which rat ##you are interested in, but we don't want that ##manually predicting things instead #times <- seq.int(0, 65, 0.1) #mcf <- model$coefficients$fixed #predicted <- # mcf["(Intercept)"] + # rep.int(mcf["Time"] * times, nlevels(BodyWeight$Diet)) + # rep(c(0, mcf["Diet2"], mcf["Diet3"]), each = length(times)) #prediction_data <- data.frame( # weight = predicted, # Time = rep.int(times, nlevels(BodyWeight$Diet)), # Diet = rep(levels(BodyWeight$Diet), each = length(times)) #) ##Draw the plot (using ggplot2) #pdf('test_lme_predict.pdf') #(p <- ggplot(BodyWeight, aes(Time, weight, colour = Diet)) + # #geom_point() + # geom_line(data = prediction_data) #) #dev.off() #predframe <- with(weeds,expand.grid(month=levels(month), # weed.cover=seq(min(weed.cover),max(weed.cover),25)) #predframe$biomass <- predict.lme(model,level=0,newdata=predframe) #then I would be tempted to use ggplot: #library(ggplot2) #ggplot(weeds,aes(x=weed.cover,y=biomass,colour=month))+ # geom_point()+ # geom_line(data=predframe) #predframe <- lme_mod_w_int_ar1[['data']] #predframe$cd4_count_sqrt <- predict(lme_mod_w_int_ar1, level=0,newdata=predframe,na.action=na.exclude) #jk <- predict(lme_mod_w_int_ar1, level=0,newdata=bmi,na.action=na.omit) @ <>= # # First trying a definition of any viral load > 400 # t2_first_suppression <- lapply(split(subset(bmi, select=c(cfar_pid, age, age_fhaart, vl_lt_400)), bmi$cfar_pid), FUN=function(y){ # diff. <- 365.25*as.numeric(y[,'age'] - y[,'age_fhaart']) # min_time <- min(diff.[which(y[,'vl_lt_400'] == 1)],na.rm=TRUE) # min_time <- ifelse(min_time == Inf,NA,min_time) # first_non_na <- min(diff.[which(!is.na(y[,'vl_lt_400']))],na.rm=TRUE) # first_non_na <- ifelse(first_non_na == Inf, NA,first_non_na) # min_time_gt_180 <- min(diff.[which(y[,'vl_lt_400'] == 0 & diff. > 180)],na.rm=TRUE) # min_time_gt_180 <- ifelse(min_time_gt_180 == Inf,NA,min_time_gt_180) # if(!is.na(min_time_gt_180)){ # age_unsuppressed_gt_180 <- y[which(diff. == min_time_gt_180),'age'] # } else { # age_unsuppressed_gt_180 <- NA # } # # df <- data.frame('cfar_pid'=y[1,'cfar_pid'], # 'min_days_to_first_suppression'=min_time, # 'min_days_to_first_non_na_vl'=first_non_na, # 'min_days_to_unsuppressed_gt_180'=min_time_gt_180, # 'age_unsuppressed_gt_180'=age_unsuppressed_gt_180) # return(df) }) #t2_first_suppression.df <- do.call('rbind',t2_first_suppression) #bmi$age_unsuppressed_gt_180 <- t2_first_suppression.df[match(bmi$cfar_pid, t2_first_suppression.df$cfar_pid, nomatch=NA),'age_unsuppressed_gt_180'] #six_months <- length(which(is.na(t2_first_suppression.df$min_days_to_unsuppressed_gt_180) | t2_first_suppression.df$min_days_to_unsuppressed_gt_180 > 180)) #twelve_months <- length(which(is.na(t2_first_suppression.df$min_days_to_unsuppressed_gt_180) | t2_first_suppression.df$min_days_to_unsuppressed_gt_180 >= 365.25)) #eighteen_months <- length(which(is.na(t2_first_suppression.df$min_days_to_unsuppressed_gt_180) | t2_first_suppression.df$min_days_to_unsuppressed_gt_180 >= 1.5*365.25)) #twentyfour_months <- length(which(is.na(t2_first_suppression.df$min_days_to_unsuppressed_gt_180) | t2_first_suppression.df$min_days_to_unsuppressed_gt_180 >= 2*365.25)) #t2_first_suppression.df$total_followup_yrs <- bmi[match(t2_first_suppression.df$cfar_pid, bmi$cfar_pid, nomatch=NA),'total_followup_yrs'] #t2_first_suppression.df$total_followup_days <- 365.25*t2_first_suppression.df$total_followup_yrs #t2_first_suppression.df$time_to_unsuppressed_days <- ifelse(is.na(t2_first_suppression.df$min_days_to_unsuppressed_gt_180),t2_first_suppression.df$total_followup_days, # t2_first_suppression.df$min_days_to_unsuppressed_gt_180) #t2_first_suppression.df$time_to_unsuppressed_yrs <- ifelse(is.na(t2_first_suppression.df$min_days_to_unsuppressed_gt_180),t2_first_suppression.df$total_followup_days/365.25, # t2_first_suppression.df$min_days_to_unsuppressed_gt_180/365.25) #t2_first_suppression.df$censor <- ifelse(is.na(t2_first_suppression.df$min_days_to_unsuppressed_gt_180),0,1) #cum_inc_unsuppressed <- cuminc(ftime=t2_first_suppression.df$time_to_unsuppressed_yrs, # fstatus=t2_first_suppression.df$censor, # cencode='0',rho=0) # # pdf('Cumulative_incidence_of_detectable_vl.pdf') # par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1)) # plot(cum_inc_unsuppressed[[1]]$est ~ cum_inc_unsuppressed[[1]]$time, type='l',axes=FALSE, # xlab='Years since ART initiation',ylab='Cumulative incidence of having a detectable viral load') # axis(1) # axis(2) # box(bty='l') # # dev.off() @ \end{document}