### R code from vignette source 'JohnRoCytokine.nw' rm(list = ls()) library(rms) library(Hmisc) library(foreign) library(ggplot2) library(grid) rawdata<-csv.get("linc_85_sub_cytokine_dataset_4jul12.csv",header=TRUE,lowernames=TRUE, na.strings=c("."," ", "#NULL!"),sep=",",stringsAsFactors =FALSE) describe(rawdata) ################################################### ### code chunk number 2: summary ################################################### x<-summary.formula(sex~age+race+sex+pi+bmi+cd4.count+hiv.viral.load+current.smoker+routine.asa.use+asa.or.nsaid.use +ln.hscrp+hscrp+ln.leptin+leptin+ln.il.6+il.6+ln.tnf.alpha+tnf.alpha+ln.il.1.beta+il.1.beta+ln.mcp.1 +mcp.1+ln.mip.1a+mip.1a+ln.cd14+cd14+tnf.a.receptor.1+ln.tnf.a.r1+ln.tnf.a.r2+tnf.a.receptor.2 +limbfatgrams+limbfatpct+trunkfat+trunkfatgrams+totalbodyfat+totalbodyfatgrams+totalbodymass+dexa.avail+age+pi, method='reverse', test=TRUE, overall=TRUE, data=rawdata) latex (x, size="small",prmsd=T,exclude1=0,npct='both',digits=2, longtable=TRUE,long=TRUE,landscape=TRUE, file="", caption="Patient Characteristics and Clinical Factors", where="!h") #anadata<-upData(anadata,drop=c("ec_surgtype","ec_comorblist","meds_cv_digind","lab_insulin_wk4")) #dd <- datadist(anadata, q.effect=c(0.25,0.75)) #options(datadist='dd') ################################################### ### code chunk number 3: JohnRoCytokine.nw:81-86 ################################################### anadata=rawdata dd <- datadist(anadata, q.effect=c(0.25,0.75)) options(datadist='dd') ################################################### ### code chunk number 4: JohnRoCytokine.nw:97-132 ################################################### c1 <- ggplot(anadata, aes(il.6, bmi)) + stat_smooth(method="loess")+ geom_point() c2 <- ggplot(anadata, aes(hscrp, bmi)) + stat_smooth(method="loess")+ geom_point() c3 <- ggplot(anadata, aes(tnf.a.receptor.1, bmi)) + stat_smooth(method="loess")+ geom_point() c4 <- ggplot(anadata, aes(il.6, limbfatgrams)) + stat_smooth() + geom_point() + scale_y_log10() c5 <- ggplot(anadata, aes(hscrp, limbfatgrams)) + stat_smooth() + geom_point() + scale_y_log10() c6 <- ggplot(anadata, aes(tnf.a.receptor.1, limbfatgrams)) + stat_smooth() + geom_point() c7 <- ggplot(anadata, aes(il.6, totalbodyfat)) + stat_smooth() + geom_point() c8 <- ggplot(anadata, aes(hscrp, totalbodyfat)) + stat_smooth() + geom_point() c9 <- ggplot(anadata, aes(tnf.a.receptor.1, totalbodyfat)) + stat_smooth() + geom_point() grid.newpage() pushViewport(viewport(layout = grid.layout(3, 3))) #print(c1, layout.pos.row=1, layout.pos.col=1) vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) print(c1, vp = vplayout(1, 1)) print(c2, vp = vplayout(1, 2)) print(c3, vp = vplayout(1, 3)) print(c4, vp = vplayout(2, 1)) print(c5, vp = vplayout(2, 2)) print(c6, vp = vplayout(2, 3)) print(c7, vp = vplayout(3, 1)) print(c8, vp = vplayout(3, 2)) print(c9, vp = vplayout(3, 3)) ################################################### ### code chunk number 5: JohnRoCytokine.nw:139-159 ################################################### ##independentall<- serumcyto<-c("il.6" ,"hscrp", "tnf.a.receptor.1","mip.1a" ,"cd14" ,"tnf.a.receptor.2","mcp.1") urinary<-c("leptin","trunkfatgrams","totalbodyfatgrams","limbfatgrams") pvalSpearman <-rho<- NULL for (i in 1:length(serumcyto)) { for (j in 1:length(urinary)) { spearman<-cor.test(anadata[[serumcyto[i]]],anadata[[urinary[j]]],method="spearman") pvalSpearman[4*(i-1)+j]<-round(spearman$p.value,3) rho[4*(i-1)+j]<-round(spearman$estimate,3) } } d1 <- expand.grid(x=serumcyto, y=urinary) corall<-data.frame(Metabolic=d1$y,urinary=d1$x,rho=rho,SpearmanPValue=pvalSpearman) corall<- cbind("Dexa"=as.character(d1$y), "Inflamation"=as.character(d1$x), '$\\rho$'=round(rho,3), 'P Value of Spearman Correlation'=ifelse(pvalSpearman<0.001, "<0.001", as.character(pvalSpearman))) ################################################### ### code chunk number 6: JohnRoCytokine.nw:162-164 ################################################### latex(corall ,cdec=c(3,3,3,3),col.just=rep('l',4), rowname=NULL, rowlabel='',caption='Spearman corrlation between Dexa variables and inflammation markers',file="", where="!h") ################################################### ### code chunk number 7: c11 ################################################### a=aregImpute(~age+sex+cd4.count+current.smoker+asa.or.nsaid.use+pi+ln.il.6+ln.hscrp+ln.mcp.1+ln.tnf.a.r2+ln.leptin+ln.il.6, data=anadata, n.impute=5) dependent<-c("ln.hscrp") independent<-c("ln.il.6" , "ln.tnf.a.r1","ln.mip.1a" ,"ln.cd14" ,"ln.tnf.a.r2","ln.mcp.1") OLSOutput=function(depen,indep){ form=as.formula(paste(depen,"~ ",indep,"+age+cd4.count+current.smoker+asa.or.nsaid.use+pi",sep="")) fit<-ols(form,data=anadata) } d1 <- expand.grid(x=independent, y=dependent) pval <- Beta <- SeBeta <-CIlow <-CIhigh<-NULL for (i in 1:1) { for (j in 1:6) { tempfit<-OLSOutput(depen=dependent[i],indep=independent[j]) pval[6*(i-1)+j]<-round(anova(tempfit)[1,5],3) Beta[6*(i-1)+j] <- exp(summary(tempfit)[1,"Effect"]) CIlow[6*(i-1)+j] <- exp(summary(tempfit)[1,"Lower 0.95"]) CIhigh[6*(i-1)+j] <- exp(summary(tempfit)[1,"Upper 0.95"]) } } #final1<-data.frame(dependent=d1$y, SerumCytokines=d1$x, c2=round(Beta,3),CILower=round(CIlow,3),CIHigh=round(CIhigh,3), pvalue=pval) b <- cbind("Model"=paste(d1$y, "$\\sim$", "c2*",d1$x,"+age+sex+cd4.count+current.smoker+asa.or.nsaid.use+asa+pi+cd_14", sep=" "), SerumCytokines=independent, 'c2 (95\\% CI)'=paste(round((Beta),3)," (",round((CIlow),3), ", ", round((CIhigh),3),")", sep=""), 'P Value of c2'=ifelse(pval<0.001, "<0.001", as.character(format(pval, nsmall=3)))) ################################################### ### code chunk number 8: JohnRoCytokine.nw:203-206 ################################################### latex(b ,cdec=c(3,3,3,3),col.just=rep('l',4), landscape=TRUE,rowname=NULL,n.rgroup=c(6),size="footnotesize", ctable=TRUE, rowlabel='',caption='The adjusted relationship between cytokines and hscrp',file="", where="htbp") ################################################### ### code chunk number 9: c11 ################################################### dependent<-c("ln.hscrp","ln.il.6" ,"ln.tnf.a.r1","ln.mip.1a" ,"ln.cd14" ,"ln.tnf.a.r2","ln.mcp.1") independent<-c("ln.leptin","trunkfatgrams","totalbodyfatgrams","limbfatgrams") OLSOutput=function(depen,indep){ form=as.formula(paste(depen,"~ ",indep,"+age+sex+cd4.count+current.smoker+asa.or.nsaid.use+pi",sep="")) fit<-ols(form,data=anadata) } d1 <- expand.grid(x=independent, y=dependent) pval <- Beta <- SeBeta <-CIlow <-CIhigh<-NULL for (i in 1:7) { for (j in 1:4) { tempfit<-OLSOutput(depen=dependent[i],indep=independent[j]) pval[4*(i-1)+j]<-round(anova(tempfit)[1,5],3) Beta[4*(i-1)+j] <- exp(summary(tempfit)[1,"Effect"]) CIlow[4*(i-1)+j] <- exp(summary(tempfit)[1,"Lower 0.95"]) CIhigh[4*(i-1)+j] <- exp(summary(tempfit)[1,"Upper 0.95"]) } } b.adiponectin <- cbind("Model"=paste(d1$y, "$\\sim$", "c2*",d1$x,"+age+sex+cd4.count+current.smoker+asa.or.nsaid.use+asa+pi+cd_14", sep=" "), SerumCytokines=independent, 'c2 (95\\% CI)'=paste(round((Beta),3)," (",round((CIlow),3), ", ", round((CIhigh),3),")", sep=""), 'P Value of c2'=ifelse(pval<0.001, "<0.001", as.character(format(pval, nsmall=3)))) ################################################### ### code chunk number 10: JohnRoCytokine.nw:239-243 ################################################### latex(b.adiponectin ,cdec=c(3,3,3,3),col.just=rep('l',4), landscape=TRUE,rowname=NULL,n.rgroup=c(4,4,4,4,4,4,4),size="footnotesize", ctable=TRUE, rowlabel='',caption='The adjusted relationship between Inflamation and DEXA',file="", where="htbp") ################################################### ### code chunk number 11: JohnRoCytokine.nw:253-283 ################################################### dependent<-c("ln.hscrp") pvalm<-NULL independent<-c("ln.leptin","trunkfatgrams","totalbodyfatgrams","limbfatgrams") OLSOutput=function(depen,indep){ form=as.formula(paste(depen,"~ ",indep,"+age+sex+cd4.count+current.smoker+asa.or.nsaid.use+pi+ln.mip.1a",sep="")) fit<-ols(form,data=anadata) } d1 <- expand.grid(x=independent, y=dependent) pval <- Beta <- SeBeta <-CIlow <-CIhigh<-NULL for (i in 1:1) { for (j in 1:4) { tempfit<-OLSOutput(depen=dependent[i],indep=independent[j]) pval[4*(i-1)+j]<-round(anova(tempfit)[1,5],3) pvalm[4*(i-1)+j]<-round(anova(tempfit)[8,5],3) Beta[4*(i-1)+j] <- exp(summary(tempfit)[1,"Effect"]) CIlow[4*(i-1)+j] <- exp(summary(tempfit)[1,"Lower 0.95"]) CIhigh[4*(i-1)+j] <- exp(summary(tempfit)[1,"Upper 0.95"]) } } b.adiponectin <- cbind("Model"=paste(d1$y, "$\\sim$", "c2*",d1$x,"+age+sex+cd4.count+current.smoker+asa.or.nsaid.use+asa+pi+cd_14+ln.mip.1a", sep=" "), SerumCytokines=independent, 'c2 (95\\% CI)'=paste(round((Beta),3)," (",round((CIlow),3), ", ", round((CIhigh),3),")", sep=""), 'P Value of c2'=ifelse(pval<0.001, "<0.001", as.character(format(pval, nsmall=3))), 'P Value of m'=ifelse(pval<0.001, "<0.001", as.character(format(pvalm, nsmall=3))) ) ################################################### ### code chunk number 12: JohnRoCytokine.nw:286-290 ################################################### latex(b.adiponectin ,cdec=c(3,3,3,3,3),col.just=rep('l',5), landscape=TRUE,rowname=NULL,n.rgroup=c(4),size="footnotesize", ctable=TRUE, rowlabel='',caption='Mediation effect of mip.1a',file="", where="htbp") ################################################### ### code chunk number 13: JohnRoCytokine.nw:299-329 ################################################### dependent<-c("ln.hscrp") pvalm<-NULL independent<-c("ln.leptin","trunkfatgrams","totalbodyfatgrams","limbfatgrams") OLSOutput=function(depen,indep){ form=as.formula(paste(depen,"~ ",indep,"+age+sex+cd4.count+current.smoker+asa.or.nsaid.use+pi+ln.tnf.a.r1",sep="")) fit<-ols(form,data=anadata) } d1 <- expand.grid(x=independent, y=dependent) pval <- Beta <- SeBeta <-CIlow <-CIhigh<-NULL for (i in 1:1) { for (j in 1:4) { tempfit<-OLSOutput(depen=dependent[i],indep=independent[j]) pval[4*(i-1)+j]<-round(anova(tempfit)[1,5],3) pvalm[4*(i-1)+j]<-round(anova(tempfit)[8,5],3) Beta[4*(i-1)+j] <- exp(summary(tempfit)[1,"Effect"]) CIlow[4*(i-1)+j] <- exp(summary(tempfit)[1,"Lower 0.95"]) CIhigh[4*(i-1)+j] <- exp(summary(tempfit)[1,"Upper 0.95"]) } } b.adiponectin <- cbind("Model"=paste(d1$y, "$\\sim$", "c2*",d1$x,"+age+sex+cd4.count+current.smoker+asa.or.nsaid.use+asa+pi+cd_14+ln.tnf.a.r1", sep=" "), SerumCytokines=independent, 'c2 (95\\% CI)'=paste(round((Beta),3)," (",round((CIlow),3), ", ", round((CIhigh),3),")", sep=""), 'P Value of c2'=ifelse(pval<0.001, "<0.001", as.character(format(pval, nsmall=3))), 'P Value of m'=ifelse(pval<0.001, "<0.001", as.character(format(pvalm, nsmall=3))) ) ################################################### ### code chunk number 14: JohnRoCytokine.nw:332-336 ################################################### latex(b.adiponectin ,cdec=c(3,3,3,3,3),col.just=rep('l',5), landscape=TRUE,rowname=NULL,n.rgroup=c(4),size="footnotesize", ctable=TRUE, rowlabel='',caption='Mediation effect of ln.tnf.a.r1',file="", where="htbp") ################################################### ### code chunk number 15: JohnRoCytokine.nw:338-340 ################################################### anova(ols(ln.hscrp~ln.leptin+age+sex+cd4.count+current.smoker+asa.or.nsaid.use+pi+ln.tnf.a.r1, data=rawdata)) summary(ols(ln.hscrp~ln.leptin+age+sex+cd4.count+current.smoker+asa.or.nsaid.use+pi+ln.tnf.a.r1, data=rawdata))