################################################### ### chunk number 1: ################################################### ## Load Hmisc, chron packages library(Hmisc) library(Design) library(Cairo) ## Create variables to describe versions of R, packages R.version <- gsub(',.*','',R.version.string) Hmisc.version <- packageDescription('Hmisc')$Version chron.version <- packageDescription('chron')$Version ################################################### ### chunk number 2: ################################################### analyze <- read.csv("abc_oneobs.csv", sep=",", header=TRUE, na.strings="NA") ## Add labels, take only subset abc.analyze <- upData(subset(analyze, select=c(id, group.pp, vent.free.death.28, vent.free.death.14)), labels=c(vent.free.death.28="Ventilator-free days during first 28 days of study", vent.free.death.14="Ventilator-free days during first 14 days of study", group.pp="Study arm")) ## Create summary of vent-free days sum.vfd <- summary.formula(group.pp ~ vent.free.death.28, data=abc.analyze, method='reverse', test=TRUE) sum.vfd.sat <- summary.formula(group.pp ~ vent.free.death.28, data=subset(abc.analyze, group.pp=='SBT+SAT'), method='reverse') sum.vfd.sbt <- summary.formula(group.pp ~ vent.free.death.28, data=subset(abc.analyze, group.pp=='SBT'), method='reverse') ################################################### ### chunk number 3: ################################################### file <- 'msmtg_18oct07_histograms.pdf' CairoPDF(file, width=8, height=10) par(mfrow=c(3, 1)) hist(abc.analyze$vent.free.death.28, plot=TRUE, main="Histogram of Ventilator-Free Days, All Patients", col="light blue", xlab="Vent-free days", xlim=c(0, 30), breaks=10) hist(subset(abc.analyze, group.pp=='SBT+SAT')$vent.free.death.28, plot=TRUE, xlim=c(0, 30), ylim=c(0, 120), breaks=10, main="Histogram of Ventilator-Free Days, Intervention Group", col="light blue", xlab="Vent-free days") hist(subset(abc.analyze, group.pp=='SBT')$vent.free.death.28, plot=TRUE, xlim=c(0, 30), ylim=c(0, 120), breaks=10, main="Histogram of Ventilator-Free Days, Control Group", col="light blue", xlab="Vent-free days") par(mfrow=c(1, 1)) invisible(dev.off()) cat('\\includegraphics{',file,'}\n',sep='') ################################################### ### chunk number 4: ################################################### latex(sum.vfd, file='', where='!htbp', prmsd=TRUE, digits=1) ################################################### ### chunk number 5: ################################################### bootdif <- function(y, g) { ## Treatment group g <- as.factor(g) ## Use the smean.cl.boot function to bootstrap means for # variable y for each treatment group (A and B); this code # uses 2000 samples, but can easily be changed a <- attr(smean.cl.boot(y[g==levels(g)[1]], B=2000, reps=TRUE), 'reps') b <- attr(smean.cl.boot(y[g==levels(g)[2]], B=2000, reps=TRUE), 'reps') ## Calculate the observed mean difference between groups meandif <- diff(tapply(y, g, mean, na.rm=TRUE)) ## Calculate the 2.5 and 97.5 percentiles of the differences # in bootstrapped means a.b <- quantile(b-a, c(.025,.975)) ## Prepare object to return res <- c(meandif, a.b) names(res) <- c('Mean','.025','.975') res } ################################################### ### chunk number 6: ################################################### head(subset(abc.analyze, select=c(id, group.pp, vent.free.death.28))) ################################################### ### chunk number 7: ################################################### ## VFDs g <- abc.analyze$group.pp y <- abc.analyze$vent.free.death.28 set.seed(1) ci.vfds <- bootdif(y, g) ci.vfds ################################################### ### chunk number 8: ################################################### meandiff.vfds <- round(ci.vfds[1], 1) lcl.vfds <- round(ci.vfds[2], 1) ucl.vfds <- round(ci.vfds[3], 1) ################################################### ### chunk number 9: ################################################### ## VFDs during first 14 days of study g <- abc.analyze$group.pp y <- abc.analyze$vent.free.death.14 set.seed(1) ci.vfds.14 <- bootdif(y, g) meandiff.vfds.14 <- round(ci.vfds.14[1], 1) lcl.vfds.14 <- round(ci.vfds.14[2], 1) ucl.vfds.14 <- round(ci.vfds.14[3], 1) boot.cis <- data.frame(outcome = c("Ventilator-Free Days During First 28 Study Days", "Ventilator-Free Days During First 14 Study Days"), numbers=c( paste(meandiff.vfds, ' (', lcl.vfds, ', ', ucl.vfds, ')', sep=''), paste(meandiff.vfds.14, ' (', lcl.vfds.14, ', ', ucl.vfds.14, ')', sep=''))) ################################################### ### chunk number 10: ################################################### latex(boot.cis, file='', rowname=NULL, where='!htbp', col.just=c('l', 'c'), caption="Difference in Means (Intervention - Usual Care) and Bootstrapped 95\\% CIs", colheads=c("Outcome", "Difference in Means (CI)"))