rm(list=ls()) #------------------------# # Loading libraries # #------------------------# library(Hmisc) library(Design) library(lmtest) library(mice) #--------------------------------# # Functions # #--------------------------------# model.fun <- function(month,data){ if(month == 12){ main.mod <- ols(cd4_change_12 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3)*sex + race + pi_usage + age_fhaart,data=data) sec.mod <- ols(cd4_change_12 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3)*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart,data=data) no.int.mod <- ols(cd4_change_12 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3) + sex + race + pi_usage + age_fhaart,data=data) sec.no.int.mod <- ols(cd4_change_12 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3) + sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart,data=data) for.qq <- lm(cd4_change_12 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3)*sex + race + pi_usage + age_fhaart,data=data) for.qq2 <- lm(cd4_change_12 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3)*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart,data=data) } else { main.mod <- ols(cd4_change_24 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3)*sex + race + pi_usage + age_fhaart,data=data) sec.mod <- ols(cd4_change_24 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3)*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart,data=data) no.int.mod <- ols(cd4_change_24 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3) + sex + race + pi_usage + age_fhaart,data=data) sec.no.int.mod <- ols(cd4_change_24 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3) + sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart,data=data) for.qq <- lm(cd4_change_12 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3)*sex + race + pi_usage + age_fhaart,data=data) for.qq2 <- lm(cd4_change_24 ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3)*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart,data=data) } return(list(main.mod=main.mod, sec.mod=sec.mod, no.int.mod=no.int.mod, sec.no.int.mod=sec.no.int.mod, for.qq=for.qq, for.qq2=for.qq2)) } summary.fun.1 <- function(mod, sex, bmi.comp){ return(summary(mod, sex=sex, race="white",pi_usage="no", cd4_count_at_age_fhaart=c(200,300), age_fhaart=c(25,26), bmi_at_age_fhaart=bmi.comp)) } summary.fun <- function(mod, sex, bmi.comp){ return(summary(mod, sex=sex, race="white",pi_usage="no", cd4_count_at_age_fhaart=c(200,300), age_fhaart=c(25,26), bmi_at_age_fhaart=bmi.comp, loghiv1_rna_at_age_fhaart=c(3,4))) } impute.fun <- function(month=month, p.ht, r.ht, data){ ht.beta <- ht.beta.2 <- ht.beta.3 <- ht.beta.m <- ht.beta.2.m <- ht.beta.3.m <- NULL ht.se <- ht.se.2 <- ht.se.3 <- ht.se.m <- ht.se.2.m <- ht.se.3.m <- NULL beta.full <- beta.wo.bmi <- beta.wo.sex <- beta.wo.int <- beta.wo.rcs <- beta.wo.int.wo.bmi <- beta.wo.int.wo.rcs <- NULL lrtest.wo.bmi <- lrtest.wo.sex <- lrtest.wo.int <- lrtest.wo.rcs <- NULL ht.beta.wo.int <- ht.beta.wo.int.2 <- ht.beta.wo.int.3 <- NULL ht.se.wo.int <- ht.se.wo.int.2 <- ht.se.wo.int.3 <- NULL lrtest.wo.int.wo.bmi <- lrtest.wo.int.wo.rcs <- NULL new.bmi.df <- matrix(data=NA,nrow=dim(data)[1],ncol=10) vcov.full.mod <- vector("list",10) vcov.mod.wo.int <- vector("list",10) set.seed(1) m <- 10 for(i in 1:m){ ht.predicted <- p.ht + sample(r.ht,1) new.ht <- ifelse(is.na(data$height_meters),ht.predicted,data$height_meters) data$new.bmi <- ifelse(is.na(data$bmi_at_age_fhaart), data$baseline_weight_kg/(new.ht)^2, data$bmi_at_age_fhaart) new.bmi.df[,i] <- data$new.bmi ddist <<- datadist(data) options(datadist='ddist') if(month == 12){ y <- data$cd4_change_12 } else { y <- data$cd4_change_24 } mod.impute <- ols(y ~ cd4_count_at_age_fhaart + rcs(new.bmi,3)*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=data) mod.wo.int.ols <- ols(y ~ cd4_count_at_age_fhaart + rcs(new.bmi,3) + sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=data) # Collecting betas from various comparisons summary.fun <- function(mod, sex, bmi.comp){ return(summary(mod, sex=sex, race="white",pi_usage="no", cd4_count_at_age_fhaart=c(200,300), age_fhaart=c(25,26), new.bmi=bmi.comp, loghiv1_rna_at_age_fhaart=c(3,4))) } summ. <- summary.fun(mod.impute, sex="female",bmi.comp=c(25,20)) summ.2 <- summary.fun(mod.impute, sex="female",bmi.comp=c(25,30)) summ.3 <- summary.fun(mod.impute, sex="female", bmi.comp=c(25,40)) summ.m <- summary.fun(mod.impute, sex="male", bmi.comp=c(25,20)) summ.2.m <- summary.fun(mod.impute, sex="male", bmi.comp=c(25,30)) summ.3.m <- summary.fun(mod.impute, sex="male", bmi.comp=c(25,40)) summ.wo.int <- summary.fun(mod.wo.int.ols, sex="female", bmi.comp=c(25,20)) summ.wo.int.2 <- summary.fun(mod.wo.int.ols, sex="female", bmi.comp=c(25,30)) summ.wo.int.3 <- summary.fun(mod.wo.int.ols, sex="female", bmi.comp=c(25,40)) # For 1 df Wald tests and reporting overall beta's for BMI and sex that includes interactions, # etc, collecting beta's and SE's from summary function ht.beta <- c(ht.beta, summ.[,"Effect"]) ht.beta.2 <- c(ht.beta.2, summ.2[,"Effect"]) ht.beta.3 <- c(ht.beta.3, summ.3[,"Effect"]) ht.beta.m <- c(ht.beta.m, summ.m[,"Effect"]) ht.beta.2.m <- c(ht.beta.2.m,summ.2.m[,"Effect"]) ht.beta.3.m <- c(ht.beta.3.m,summ.3.m[,"Effect"]) ht.beta.wo.int <- c(ht.beta.wo.int, summ.wo.int[,"Effect"]) ht.beta.wo.int.2 <- c(ht.beta.wo.int.2, summ.wo.int.2[,"Effect"]) ht.beta.wo.int.3 <- c(ht.beta.wo.int.3, summ.wo.int.3[,"Effect"]) ht.se <- c(ht.se, summ.[,"S.E."]) ht.se.2 <- c(ht.se.2, summ.2[,"S.E."]) ht.se.3 <- c(ht.se.3, summ.3[,"S.E."]) ht.se.m <- c(ht.se.m, summ.m[,"S.E."]) ht.se.2.m <- c(ht.se.2.m, summ.2.m[,"S.E."]) ht.se.3.m <- c(ht.se.3.m, summ.3.m[,"S.E."]) ht.se.wo.int <- c(ht.se.wo.int, summ.wo.int[,"S.E."]) ht.se.wo.int.2 <- c(ht.se.wo.int.2, summ.wo.int.2[,"S.E."]) ht.se.wo.int.3 <- c(ht.se.wo.int.3, summ.wo.int.3[,"S.E."]) # For modified LR test method to test > 1 df, must collect unrestricted LR test statistics # Will do so using lrtest on lm rather than ols # Information for likelihood ratio tests full.mod <- lm(y ~ cd4_count_at_age_fhaart + rcs(new.bmi,3)*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=data) mod.wo.bmi <- lm(y ~ cd4_count_at_age_fhaart + sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=data, subset=!is.na(data$new.bmi)) mod.wo.sex <- lm(y ~ cd4_count_at_age_fhaart + rcs(new.bmi,3) + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=data) mod.wo.int <- lm(y ~ cd4_count_at_age_fhaart + rcs(new.bmi,3) + sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=data) mod.wo.rcs <- lm(y ~ cd4_count_at_age_fhaart + new.bmi*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=data) mod.wo.int.wo.bmi <- lm(y ~ cd4_count_at_age_fhaart + sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=data) mod.wo.int.wo.rcs <- lm(y ~ cd4_count_at_age_fhaart + new.bmi + sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=data) # Variance/Covariance matrix from full model for modified Wald test vcov.full.mod[[i]] <- vcov(full.mod) # Variance/Covariance matrix from model w/o int for modified Wald test vcov.mod.wo.int[[i]] <- vcov(mod.wo.int) # Collecting unrestricted LR test statistics lrtest.bmi <- lrtest(full.mod,mod.wo.bmi) lrtest.wo.bmi <- c(lrtest.wo.bmi, lrtest.bmi$Chisq[2]) lrtest.sex <- lrtest(full.mod,mod.wo.sex) lrtest.wo.sex <- c(lrtest.wo.sex, lrtest.sex$Chisq[2]) lrtest.int <- lrtest(full.mod,mod.wo.int) lrtest.wo.int <- c(lrtest.wo.int, lrtest.int$Chisq[2]) lrtest.rcs <- lrtest(full.mod,mod.wo.rcs) lrtest.wo.rcs <- c(lrtest.wo.rcs, lrtest.rcs$Chisq[2]) lrtest.no.int.bmi <- lrtest(mod.wo.int,mod.wo.int.wo.bmi) lrtest.wo.int.wo.bmi <- c(lrtest.wo.int.wo.bmi, lrtest.no.int.bmi$Chisq[2]) lrtest.no.int.rcs <- lrtest(mod.wo.int,mod.wo.int.wo.rcs) lrtest.wo.int.wo.rcs <- c(lrtest.wo.int.wo.rcs, lrtest.no.int.rcs$Chisq[2]) # Grabbing df of each of the LR tests k.bmi <- abs(lrtest.bmi$Df[2]) k.sex <- abs(lrtest.sex$Df[2]) k.int <- abs(lrtest.int$Df[2]) k.rcs <- abs(lrtest.rcs$Df[2]) k.no.int.bmi <- abs(lrtest.no.int.bmi$Df[2]) k.no.int.rcs <- abs(lrtest.no.int.rcs$Df[2]) beta.full <- c(beta.full, full.mod$coeff[-1]) beta.wo.bmi <- c(beta.wo.bmi, mod.wo.bmi$coeff[-1]) beta.wo.sex <- c(beta.wo.sex, mod.wo.sex$coeff[-1]) beta.wo.int <- c(beta.wo.int, mod.wo.int$coeff[-1]) beta.wo.rcs <- c(beta.wo.rcs, mod.wo.rcs$coeff[-1]) beta.wo.int.wo.bmi <- c(beta.wo.int.wo.bmi, mod.wo.int.wo.bmi$coeff[-1]) beta.wo.int.wo.rcs <- c(beta.wo.int.wo.rcs, mod.wo.int.wo.rcs$coeff[-1]) } beta.var.summary <- function(mod.beta, mod.se){ beta.overall <- tapply(as.numeric(mod.beta), INDEX=names(mod.beta), FUN=mean) btwn.var <- tapply(as.numeric(mod.beta), INDEX=names(mod.beta), FUN=var) mean.var <- tapply(as.numeric(mod.se)^2, INDEX=names(mod.beta), FUN=mean) var.overall <- mean.var + 11/10*btwn.var overall.lci <- beta.overall - qnorm(0.975)*sqrt(var.overall) overall.uci <- beta.overall + qnorm(0.975)*sqrt(var.overall) wald.st <- beta.overall/sqrt(var.overall) p.val <- (1-pchisq((wald.st)^2,1)) p.val <- ifelse(p.val < 0.01, format(round(p.val,3),nsmall=3),format(round(p.val,2),nsmall=2)) # Formatting beta.overall <- format(round(beta.overall,2),nsmall=2) overall.ci <- paste("(",format(round(overall.lci,2),nsmall=2),", ", format(round(overall.uci,2),nsmall=2),")",sep="") return(list(beta.overall=beta.overall, overall.ci=overall.ci, p.val=p.val)) } # Original model (comparing BMI 20 vs 25) orig.mod.info <- beta.var.summary(mod.beta=ht.beta, mod.se=ht.se) orig.mod.info.m <- beta.var.summary(mod.beta=ht.beta.m, mod.se=ht.se.m) # adjusting for 'male' # '.2' model (comparing BMI 25 vs 30) sec.mod.info <- beta.var.summary(mod.beta=ht.beta.2, mod.se=ht.se.2) sec.mod.info.m <- beta.var.summary(mod.beta=ht.beta.2.m, mod.se=ht.se.2.m) # adjusting for 'male' # '.3' model (comparing BMI 25 vs 40) third.mod.info <- beta.var.summary(mod.beta=ht.beta.3, mod.se=ht.se.3) third.mod.info.m <- beta.var.summary(mod.beta=ht.beta.3.m, mod.se=ht.se.3.m) # adjusting for 'male' # No interaction model (comparing BMI 25 vs 20) mod.wo.int.ols.info <- beta.var.summary(mod.beta=ht.beta.wo.int, mod.se=ht.se.wo.int) # No interaction model (comparing BMI 25 vs 30) mod.wo.int.ols.info.2 <- beta.var.summary(mod.beta=ht.beta.wo.int.2, mod.se=ht.se.wo.int.2) # No interaction model (comparing BMI 25 vs 40) mod.wo.int.ols.info.3 <- beta.var.summary(mod.beta=ht.beta.wo.int.3, mod.se=ht.se.wo.int.3) # Likelihood ratio beta's beta.full.overall <- tapply(as.numeric(beta.full), INDEX=names(beta.full), FUN=mean) beta.wo.bmi.overall <- tapply(as.numeric(beta.wo.bmi), INDEX=names(beta.wo.bmi), FUN=mean) beta.wo.sex.overall <- tapply(as.numeric(beta.wo.sex), INDEX=names(beta.wo.sex), FUN=mean) beta.wo.int.overall <- tapply(as.numeric(beta.wo.int), INDEX=names(beta.wo.int), FUN=mean) beta.wo.rcs.overall <- tapply(as.numeric(beta.wo.rcs), INDEX=names(beta.wo.rcs), FUN=mean) beta.wo.int.wo.bmi.overall <- tapply(as.numeric(beta.wo.int.wo.bmi), INDEX=names(beta.wo.int.wo.bmi), FUN=mean) beta.wo.int.wo.rcs.overall <- tapply(as.numeric(beta.wo.int.wo.rcs), INDEX=names(beta.wo.int.wo.rcs), FUN=mean) # Means of likelihood ratio test statistics lrtest.wo.bmi.mean <- mean(as.numeric(lrtest.wo.bmi)) lrtest.wo.sex.mean <- mean(as.numeric(lrtest.wo.sex)) lrtest.wo.int.mean <- mean(as.numeric(lrtest.wo.int)) lrtest.wo.rcs.mean <- mean(as.numeric(lrtest.wo.rcs)) lrtest.wo.int.wo.bmi.mean <- mean(as.numeric(lrtest.wo.int.wo.bmi)) lrtest.wo.int.wo.rcs.mean <- mean(as.numeric(lrtest.wo.int.wo.rcs)) # For each of m imputed data sets, calculate the likelihood for h1 with parameters constrained to # the mean. Repeat for the reduced model. lrtest.constrained.wo.bmi <- lrtest.constrained.wo.sex <- lrtest.constrained.wo.int <- NULL lrtest.constrained.wo.rcs <- NULL lrtest.constrained.wo.int.wo.bmi <- lrtest.constrained.wo.int.wo.rcs <- NULL for(i in 1:10){ data$new.bmi <- new.bmi.df[,i] rcs.bmi <- rcs(new.bmi.df[,i],3) data$linear.new.bmi <- rcs.bmi[,1] data$nonlinear.new.bmi <- rcs.bmi[,2] data$linear.new.bmi.w.sex <- with(data, linear.new.bmi*(as.numeric(sex)-1)) data$nonlinear.new.bmi.w.sex <- with(data, nonlinear.new.bmi*(as.numeric(sex)-1)) data$new.bmi.w.sex <- with(data, new.bmi*(as.numeric(sex)-1)) full.mod <- lm(y ~ offset(beta.full.overall[names(beta.full.overall)=="cd4_count_at_age_fhaart"]*cd4_count_at_age_fhaart) + offset(beta.full.overall[names(beta.full.overall)=="rcs(new.bmi, 3)new.bmi"]*linear.new.bmi) + offset(beta.full.overall[names(beta.full.overall)=="rcs(new.bmi, 3)new.bmi'"]*nonlinear.new.bmi) + offset(beta.full.overall[names(beta.full.overall)=="rcs(new.bmi, 3)new.bmi:sexmale"]*linear.new.bmi.w.sex) + offset(beta.full.overall[names(beta.full.overall)=="rcs(new.bmi, 3)new.bmi':sexmale"]*nonlinear.new.bmi.w.sex) + offset(beta.full.overall[names(beta.full.overall)=="sexmale"]*(as.numeric(sex)-1)) + offset(beta.full.overall[names(beta.full.overall)=="racenonwhite"]*(as.numeric(race)-1)) + offset(beta.full.overall[names(beta.full.overall)=="pi_usageyes"]*(as.numeric(pi_usage)-1)) + offset(beta.full.overall[names(beta.full.overall)=="age_fhaart"]*age_fhaart) + offset(beta.full.overall[names(beta.full.overall)=="loghiv1_rna_at_age_fhaart"]*loghiv1_rna_at_age_fhaart), data=data) mod.wo.bmi <- lm(y ~ offset(beta.wo.bmi.overall[names(beta.wo.bmi.overall)=="cd4_count_at_age_fhaart"]*cd4_count_at_age_fhaart) + offset(beta.wo.bmi.overall[names(beta.wo.bmi.overall)=="sexmale"]*(as.numeric(sex)-1)) + offset(beta.wo.bmi.overall[names(beta.wo.bmi.overall)=="racenonwhite"]*(as.numeric(race)-1)) + offset(beta.wo.bmi.overall[names(beta.wo.bmi.overall)=="pi_usageyes"]*(as.numeric(pi_usage)-1)) + offset(beta.wo.bmi.overall[names(beta.wo.bmi.overall)=="age_fhaart"]*age_fhaart) + offset(beta.wo.bmi.overall[names(beta.wo.bmi.overall)=="loghiv1_rna_at_age_fhaart"]*loghiv1_rna_at_age_fhaart), data=data) mod.wo.sex <- lm(y ~ offset(beta.wo.sex.overall[names(beta.wo.sex.overall)=="cd4_count_at_age_fhaart"]*cd4_count_at_age_fhaart) + offset(beta.wo.sex.overall[names(beta.wo.sex.overall)=="rcs(new.bmi, 3)new.bmi"]*linear.new.bmi) + offset(beta.wo.sex.overall[names(beta.wo.sex.overall)=="rcs(new.bmi, 3)new.bmi'"]*nonlinear.new.bmi) + offset(beta.wo.sex.overall[names(beta.wo.sex.overall)=="racenonwhite"]*(as.numeric(race)-1)) + offset(beta.wo.sex.overall[names(beta.wo.sex.overall)=="pi_usageyes"]*(as.numeric(pi_usage)-1)) + offset(beta.wo.sex.overall[names(beta.wo.sex.overall)=="age_fhaart"]*age_fhaart) + offset(beta.wo.sex.overall[names(beta.wo.sex.overall)=="loghiv1_rna_at_age_fhaart"]*loghiv1_rna_at_age_fhaart), data=data) mod.wo.int <- lm(y ~ offset(beta.wo.int.overall[names(beta.wo.int.overall)=="cd4_count_at_age_fhaart"]*cd4_count_at_age_fhaart) + offset(beta.wo.int.overall[names(beta.wo.int.overall)=="rcs(new.bmi, 3)new.bmi"]*linear.new.bmi) + offset(beta.wo.int.overall[names(beta.wo.int.overall)=="rcs(new.bmi, 3)new.bmi'"]*nonlinear.new.bmi) + offset(beta.wo.int.overall[names(beta.wo.int.overall)=="sexmale"]*(as.numeric(sex)-1)) + offset(beta.wo.int.overall[names(beta.wo.int.overall)=="racenonwhite"]*(as.numeric(race)-1)) + offset(beta.wo.int.overall[names(beta.wo.int.overall)=="pi_usageyes"]*(as.numeric(pi_usage)-1)) + offset(beta.wo.int.overall[names(beta.wo.int.overall)=="age_fhaart"]*age_fhaart) + offset(beta.wo.int.overall[names(beta.wo.int.overall)=="loghiv1_rna_at_age_fhaart"]*loghiv1_rna_at_age_fhaart), data=data) mod.wo.rcs <- lm(y ~ offset(beta.wo.rcs.overall[names(beta.wo.rcs.overall)=="cd4_count_at_age_fhaart"]*cd4_count_at_age_fhaart) + offset(beta.wo.rcs.overall[names(beta.wo.rcs.overall)=="new.bmi"]*new.bmi) + offset(beta.wo.rcs.overall[names(beta.wo.rcs.overall)=="new.bmi:sexmale"]*new.bmi.w.sex) + offset(beta.wo.rcs.overall[names(beta.wo.rcs.overall)=="sexmale"]*(as.numeric(sex)-1)) + offset(beta.wo.rcs.overall[names(beta.wo.rcs.overall)=="racenonwhite"]*(as.numeric(race)-1)) + offset(beta.wo.rcs.overall[names(beta.wo.rcs.overall)=="pi_usageyes"]*(as.numeric(pi_usage)-1)) + offset(beta.wo.rcs.overall[names(beta.wo.rcs.overall)=="age_fhaart"]*age_fhaart) + offset(beta.wo.rcs.overall[names(beta.wo.rcs.overall)=="loghiv1_rna_at_age_fhaart"]*loghiv1_rna_at_age_fhaart), data=data) mod.wo.int.wo.bmi <- lm(y ~ offset(beta.wo.int.wo.bmi.overall[names(beta.wo.int.wo.bmi.overall)=="cd4_count_at_age_fhaart"]* cd4_count_at_age_fhaart) + offset(beta.wo.int.wo.bmi.overall[names(beta.wo.int.wo.bmi.overall)=="sexmale"]*(as.numeric(sex)-1)) + offset(beta.wo.int.wo.bmi.overall[names(beta.wo.int.wo.bmi.overall)=="racenonwhite"]*(as.numeric(race)-1)) + offset(beta.wo.int.wo.bmi.overall[names(beta.wo.int.wo.bmi.overall)=="pi_usageyes"]*(as.numeric(pi_usage)-1)) + offset(beta.wo.int.wo.bmi.overall[names(beta.wo.int.wo.bmi.overall)=="age_fhaart"]*age_fhaart) + offset(beta.wo.int.wo.bmi.overall[names(beta.wo.int.wo.bmi.overall)=="loghiv1_rna_at_age_fhaart"]* loghiv1_rna_at_age_fhaart), data=data) mod.wo.int.wo.rcs <- lm(y ~ offset(beta.wo.int.wo.rcs.overall[names(beta.wo.int.wo.rcs.overall)=="cd4_count_at_age_fhaart"]* cd4_count_at_age_fhaart) + offset(beta.wo.int.wo.rcs.overall[names(beta.wo.int.wo.rcs.overall)=="new.bmi"]*new.bmi) + offset(beta.wo.int.wo.rcs.overall[names(beta.wo.int.wo.rcs.overall)=="sexmale"]*(as.numeric(sex)-1)) + offset(beta.wo.int.wo.rcs.overall[names(beta.wo.int.wo.rcs.overall)=="racenonwhite"]*(as.numeric(race)- 1)) + offset(beta.wo.int.wo.rcs.overall[names(beta.wo.int.wo.rcs.overall)=="pi_usageyes"]*(as.numeric(pi_usage )-1)) + offset(beta.wo.int.wo.rcs.overall[names(beta.wo.int.wo.rcs.overall)=="age_fhaart"]*age_fhaart) + offset(beta.wo.int.wo.rcs.overall[names(beta.wo.int.wo.rcs.overall)=="loghiv1_rna_at_age_fhaart"]* loghiv1_rna_at_age_fhaart), data=data) lrtest.constrained.wo.bmi <- c(lrtest.constrained.wo.bmi, lrtest(full.mod,mod.wo.bmi)[2,"Chisq"]) lrtest.constrained.wo.sex <- c(lrtest.constrained.wo.sex, lrtest(full.mod,mod.wo.sex)[2,"Chisq"]) lrtest.constrained.wo.int <- c(lrtest.constrained.wo.int, lrtest(full.mod,mod.wo.int)[2,"Chisq"]) lrtest.constrained.wo.rcs <- c(lrtest.constrained.wo.rcs, lrtest(full.mod,mod.wo.rcs)[2,"Chisq"]) lrtest.constrained.wo.int.wo.bmi <- c(lrtest.constrained.wo.int.wo.bmi, lrtest(mod.wo.int,mod.wo.int.wo.bmi)[2,"Chisq"]) lrtest.constrained.wo.int.wo.rcs <- c(lrtest.constrained.wo.int.wo.rcs, lrtest(mod.wo.int, mod.wo.int.wo.rcs)[2,"Chisq"]) } lrtest.constrained.wo.bmi.mean <- mean(lrtest.constrained.wo.bmi) lrtest.constrained.wo.sex.mean <- mean(lrtest.constrained.wo.sex) lrtest.constrained.wo.int.mean <- mean(lrtest.constrained.wo.int) lrtest.constrained.wo.rcs.mean <- mean(lrtest.constrained.wo.rcs) lrtest.constrained.wo.int.wo.bmi.mean <- mean(lrtest.constrained.wo.int.wo.bmi) lrtest.constrained.wo.int.wo.rcs.mean <- mean(lrtest.constrained.wo.int.wo.rcs) # Calculating the test statistic for the modified LR test D_l.wo.bmi <- lrtest.constrained.wo.bmi.mean/(k.bmi + ((m+1)/(m-1))* (lrtest.wo.bmi.mean - lrtest.constrained.wo.bmi.mean)) D_l.wo.sex <- lrtest.constrained.wo.sex.mean/(k.sex + ((m+1)/(m-1))* (lrtest.wo.sex.mean - lrtest.constrained.wo.sex.mean)) D_l.wo.int <- lrtest.constrained.wo.int.mean/(k.int + ((m+1)/(m-1))* (lrtest.wo.int.mean - lrtest.constrained.wo.int.mean)) D_l.wo.rcs <- lrtest.constrained.wo.rcs.mean/(k.rcs + ((m+1)/(m-1))* (lrtest.wo.rcs.mean - lrtest.constrained.wo.rcs.mean)) D_l.wo.int.wo.bmi <- lrtest.constrained.wo.int.wo.bmi.mean/(k.no.int.bmi + ((m+1)/(m-1))* (lrtest.wo.int.wo.bmi.mean - lrtest.constrained.wo.int.wo.bmi.mean)) D_l.wo.int.wo.rcs <- lrtest.constrained.wo.int.wo.rcs.mean/(k.no.int.rcs + ((m+1)/(m-1))* (lrtest.wo.int.wo.rcs.mean - lrtest.constrained.wo.int.wo.rcs.mean)) # Calculating degress of freedom nu.bmi <- k.bmi*(m-1) # k*(m-1) where k=df1 - df0 and m=# imputation interations nu.sex <- k.sex*(m-1) nu.int <- k.int*(m-1) nu.rcs <- k.rcs*(m-1) nu.wo.int.wo.bmi <- k.no.int.bmi*(m-1) nu.wo.int.wo.rcs <- k.no.int.rcs*(m-1) r_l.wo.bmi <- (m+1)/(k.bmi*(m-1))*(lrtest.wo.bmi.mean - lrtest.constrained.wo.bmi.mean) r_l.wo.sex <- (m+1)/(k.sex*(m-1))*(lrtest.wo.sex.mean - lrtest.constrained.wo.sex.mean) r_l.wo.int <- (m+1)/(k.int*(m-1))*(lrtest.wo.int.mean - lrtest.constrained.wo.int.mean) r_l.wo.rcs <- (m+1)/(k.rcs*(m-1))*(lrtest.wo.rcs.mean - lrtest.constrained.wo.rcs.mean) r_l.wo.int.wo.bmi <- (m+1)/(k.no.int.bmi*(m-1))* (lrtest.wo.int.wo.bmi.mean - lrtest.constrained.wo.int.wo.bmi.mean) r_l.wo.int.wo.rcs <- (m+1)/(k.no.int.rcs*(m-1))* (lrtest.wo.int.wo.rcs.mean - lrtest.constrained.wo.int.wo.rcs.mean) # Since nu > 4, use df = r + (nu-4){1+(1-2*nu^{-1})*r_l^{-1}}^2 df.wo.bmi <- 4 + (nu.bmi-4)*{1 + ((1 - 2/nu.bmi)/r_l.wo.bmi)}^2 df.wo.sex <- 4 + (nu.sex-4)*{1 + ((1 - 2/nu.sex)/r_l.wo.sex)}^2 df.wo.int <- 4 + (nu.int-4)*{1 + ((1 - 2/nu.int)/r_l.wo.int)}^2 df.wo.rcs <- 4 + (nu.rcs-4)*{1 + ((1 - 2/nu.rcs)/r_l.wo.rcs)}^2 df.wo.int.wo.bmi <- 4 + (nu.wo.int.wo.bmi-4)*{1 + ((1 - 2/nu.wo.int.wo.bmi)/r_l.wo.int.wo.bmi)}^2 df.wo.int.wo.rcs <- 4 + (nu.wo.int.wo.rcs-4)*{1 + ((1 - 2/nu.wo.int.wo.rcs)/r_l.wo.int.wo.rcs)}^2 p.wo.bmi <- 1- pf(q=D_l.wo.bmi, df1=k.bmi, df2=df.wo.bmi) p.wo.sex <- 1 - pf(q=D_l.wo.sex, df1=k.sex, df2=df.wo.sex) p.wo.int <- 1 - pf(q=D_l.wo.int, df1=k.int, df2=df.wo.int) p.wo.rcs <- 1 - pf(q=D_l.wo.rcs, df1=k.rcs, df2=df.wo.rcs) p.wo.int.wo.bmi <- 1 - pf(q=D_l.wo.int.wo.bmi, df1=k.no.int.bmi, df2=df.wo.int.wo.bmi) p.wo.int.wo.rcs <- 1 - pf(q=D_l.wo.int.wo.rcs, df1=k.no.int.rcs, df2=df.wo.int.wo.rcs) # Modified Wald test # First must calculate averaged variance/covariance from imputation # Formula for this is: # V = M^{-1}Sum(V_i) + (M+1/M)*B # where B is var/cov matrix from the matrix of parameter estimates. Each column is # one of the parameters, each row one of the iterations in the imputation # To create B matrix must take beta.full coefficients and rearrange using a tapply # and then a do.call. # Full model calculations junk <- cbind(as.data.frame(beta.full[names(beta.full)=="cd4_count_at_age_fhaart"]), as.data.frame(beta.full[names(beta.full)=="rcs(new.bmi, 3)new.bmi"]), as.data.frame(beta.full[names(beta.full)=="rcs(new.bmi, 3)new.bmi'"]), as.data.frame(beta.full[names(beta.full)=="sexmale"]), as.data.frame(beta.full[names(beta.full)=="racenonwhite"]), as.data.frame(beta.full[names(beta.full)=="pi_usageyes"]), as.data.frame(beta.full[names(beta.full)=="age_fhaart"]), as.data.frame(beta.full[names(beta.full)=="loghiv1_rna_at_age_fhaart"]), as.data.frame(beta.full[names(beta.full)=="rcs(new.bmi, 3)new.bmi:sexmale"]), as.data.frame(beta.full[names(beta.full)=="rcs(new.bmi, 3)new.bmi':sexmale"])) colnames(junk) <- c("cd4_count_at_age_fhaart", "rcs(new.bmi, 3)new.bmi", "rcs(new.bmi, 3)new.bmi'","sexmale", "racenonwhite","pi_usageyes", "age_fhaart","loghiv1_rna_at_age_fhaart", "rcs(new.bmi, 3)new.bmi:sexmale", "rcs(new.bmi, 3)new.bmi':sexmale") # rownames(junk) <- colnames(junk) B.matrix <- var(junk) matSums <- function(lis, na.rm=FALSE){ if(!is.list(lis) || !all(sapply(lis, is.matrix))){ stop("'lis' must be a list containing 2-dimensional arrays") } dims <- sapply(lis, dim) n <- dims[1, 1] p <- dims[2, 1] if(!all(n == dims[1, ]) || !all(p == dims[2, ])){ stop("the matrices must have the same dimensions") } mat <- matrix(unlist(lis), n * p, length(lis)) matrix(rowSums(mat, na.rm=na.rm), n, p) } #----------------------------------------------------------------------------------------------# #----------------------------------# # Loading data # #----------------------------------# load("bmi.rda") # in long format #----------------------------------# # Virologically suppressed # #----------------------------------# # 12 months bmi$cd4_change_12 <- with(bmi, cd4_count_12_months - cd4_count_at_age_fhaart) vl.suppressed.12.df <- subset(bmi, !duplicated(cfar_pid) & vl_suppression_at_12_months == 1) ddist.12 <- datadist(vl.suppressed.12.df) options(datadist = 'ddist.12') month.12.vl.supp <- model.fun(month=12,data=vl.suppressed.12.df) # 24 months bmi$cd4_change_24 <- with(bmi, cd4_count_24_months - cd4_count_at_age_fhaart) vl.suppressed.24.df <- subset(bmi, !duplicated(cfar_pid) & vl_suppression_at_24_months == 1) ddist.24 <- datadist(vl.suppressed.24.df) options(datadist = 'ddist.24') month.24.vl.supp <- model.fun(month=24,data=vl.suppressed.24.df) #----------------------------------# # All comers # #----------------------------------# ddist.all <- datadist(subset(bmi, !duplicated(cfar_pid))) options(datadist = 'ddist.all') # 12 months month.12.all <- model.fun(month=12,data=subset(bmi, !duplicated(cfar_pid))) # 24 months month.24.all <- model.fun(month=24,data=subset(bmi, !duplicated(cfar_pid))) #----------------------------------------------# # Assigning models in list to object # #----------------------------------------------# # Models with just virally suppressed subjects mod.12.vlsupp.wo.hiv <- month.12.vl.supp[[1]] mod.12.vlsupp.w.hiv <- month.12.vl.supp[[2]] mod.12.vlsupp.wo.hiv.no.int <- month.12.vl.supp[[3]] mod.12.vlsupp.w.hiv.no.int <- month.12.vl.supp[[4]] # Models with all comers # 12 months mod.12.all.wo.hiv <- month.12.all[[1]] mod.12.all.w.hiv <- month.12.all[[2]] mod.12.all.wo.hiv.no.int <- month.12.all[[3]] mod.12.all.w.hiv.no.int <- month.12.all[[4]] #---------------------------------# # Imputation for # # virally suppressed # #---------------------------------# # 12 months impute.vl.12 <- glm(height_meters ~ age_fhaart + sex + cd4_count_12_months + cd4_count_at_age_fhaart,data=vl.suppressed.12.df, family="gaussian") newdata.ht <- data.frame("age_fhaart"=vl.suppressed.12.df$age_fhaart, "sex"=vl.suppressed.12.df$sex, "cd4_count_12_months"=vl.suppressed.12.df$cd4_count_12_months, "cd4_count_at_age_fhaart"=vl.suppressed.12.df$cd4_count_at_age_fhaart) sub.p.ht <- predict(impute.vl.12, newdata=newdata.ht,na.action=na.omit) p.ht <- vector(mode="numeric",length=dim(vl.suppressed.12.df)[1]) p.ht[1:length(p.ht)] <- NA p.ht[as.numeric(names(sub.p.ht))] <- sub.p.ht ht.na <- length(which(is.na(vl.suppressed.12.df$height_meters))) r.ht <- residuals(impute.vl.12) impute.12.vl.fun <- impute.fun(month=12, p.ht, r.ht, data=vl.suppressed.12.df) #---------------------------------# # Imputation for # # all comers # #---------------------------------# # 12 months sub.bmi <- subset(bmi, !duplicated(cfar_pid)) impute.all.12 <- glm(height_meters ~ age_fhaart + sex + cd4_count_12_months + cd4_count_at_age_fhaart, data=sub.bmi,family="gaussian") newdata.ht <- data.frame("age_fhaart"=sub.bmi$age_fhaart, "sex"=sub.bmi$sex, "cd4_count_12_months"=sub.bmi$cd4_count_12_months, "cd4_count_at_age_fhaart"=sub.bmi$cd4_count_at_age_fhaart) sub.p.ht <- predict(impute.all.12, newdata=newdata.ht,na.action=na.omit) p.ht <- vector(mode="numeric",length=dim(sub.bmi)[1]) p.ht[1:length(p.ht)] <- NA p.ht[as.numeric(names(sub.p.ht))] <- sub.p.ht ht.na <- length(which(is.na(sub.bmi$height_meters))) r.ht <- residuals(impute.all.12) impute.12.all.fun <- impute.fun(month=12,p.ht, r.ht, data=sub.bmi) #---------------------------------# # Analysis in those # # with a weight change < 10% # #---------------------------------# bmi$pct_weight_change <- with(bmi, (weight_12_months - baseline_weight_kgs)/baseline_weight_kgs*100) bmi$excessive_weight_change <- ifelse(abs(bmi$pct_weight_change) > 10.0,1,0) # Subsetting bmi on those who have a percent weight change of less than 10% sub. <- subset(bmi, excessive_weight_change == 0 & !duplicated(cfar_pid)) ddist.12 <- datadist(sub.) options(datadist = 'ddist.12') mod.12.wt <- model.fun(month=12,data=sub.) #-----------------------------------# # Effect of BMI on reaching CD4 # # thresholds of 200/350 # #-----------------------------------# # Are baseline BMI and CD4 count associated with having a 12 month CD4 >= 200/350 in those who had a baseline CD4 < 200 sub. <- subset(bmi, cd4_count_at_age_fhaart < 200 & !duplicated(cfar_pid)) sub.$bmi.sq <- sub.$bmi_at_age_fhaart^2 sub.$cd4_count_12_months_200_ind <- ifelse(sub.$cd4_count_12_months >= 200 & !is.na(sub.$cd4_count_12_months),1, ifelse(!is.na(sub.$cd4_count_12_months),0,NA)) sub.$cd4_count_12_months_350_ind <- ifelse(sub.$cd4_count_12_months >= 350 & !is.na(sub.$cd4_count_12_months),1, ifelse(!is.na(sub.$cd4_count_12_months),0,NA)) # Fitting models lrm.gt.200 <- lrm(cd4_count_12_months_200_ind ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3)*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=sub.) lrm.gt.200.no.int <- lrm(cd4_count_12_months_200_ind ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3) + sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=sub.) lrm.gt.200.sq <- lrm(cd4_count_12_months_200_ind ~ cd4_count_at_age_fhaart + bmi_at_age_fhaart*sex + bmi.sq*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=sub.) lrm.gt.200.sq.2 <- glm(cd4_count_12_months_200_ind ~ cd4_count_at_age_fhaart + bmi_at_age_fhaart*sex + bmi.sq*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=sub., family="binomial") lrm.gt.350 <- lrm(cd4_count_12_months_350_ind ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3)*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=sub.) lrm.gt.350.no.int <- lrm(cd4_count_12_months_350_ind ~ cd4_count_at_age_fhaart + rcs(bmi_at_age_fhaart,3) + sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=sub.) lrm.gt.350.sq <- lrm(cd4_count_12_months_350_ind ~ cd4_count_at_age_fhaart + bmi_at_age_fhaart*sex + bmi.sq*sex + race + pi_usage + age_fhaart + loghiv1_rna_at_age_fhaart, data=sub.)