rm(list=ls()) load("df12.rda") d<-df2 table(d$male[!duplicated(d$pid)]) n<-length(d$male[!duplicated(d$pid)]) sum(d$male[!duplicated(d$pid)])/n b<-d[which(!duplicated(d$pid)),] sum(b$male)/n sum(b$black)/n sum(b$idu)/n summary(b$baseline.age) par(mfrow=c(1,1)) ds<-d[which(!is.na(d$y.cd4.ade)),] bs<-b[which(!is.na(b$y.cd4.ade)),] hist(ds$x,nclass=20) hist(b$y.cd4.ade,nclass=20) w.tot<-d$w*d$v1*d$v2*d$z hist(w.tot) q.w<-quantile(w.tot,c(.025,.975)) ## old was .05, .95 w.t<-ds$w*ds$v1*ds$v2*ds$z w.t<-w.tot w.t<-ifelse(w.tq.w[2],q.w[2],w.t)) library(Design) y.cd4.ade<-d$y.cd4.ade y.complete.ade<-d$y.complete.ade x<-d$x w.t<-as.numeric(w.t) ddist <- datadist(y.cd4.ade,y.complete.ade,x,w.t) options(datadist='ddist') mod.cd4.ade <- ols(y.cd4.ade ~ rcs(x,6),weights=w.t) stuff<-plot(mod.cd4.ade) dev.off() mod2<- ols(y.complete.ade ~ rcs(x,6), weights=w.t) stuff2<-plot(mod2) dev.off() pred1<-predict(mod.cd4.ade) max.x.1<-unique(x[pred1==max(pred1,na.rm=TRUE)&!is.na(pred1)]) pred2<-predict(mod2) max.x.2<-unique(x[pred2==max(pred2,na.rm=TRUE)&!is.na(pred2)]) test <- predict(mod.cd4.ade) test.max.cd4.ade <- x[test==max(test,na.rm=TRUE) & !is.na(test)][1] library(survey) library(hexbin) #### function from Thomas Lumley svyplot<-function (formula, design, style = c("bubble", "hex", "grayhex", "subsample", "transparent"), sample.size = 500, subset = NULL, legend = 1, inches = 0.05, amount = NULL, basecol = "black", alpha = c(0, 0.8), xbins = 30, ...) { style <- match.arg(style) if (style %in% c("hex", "grayhex") && !require(hexbin)) { stop(style, " plots require the hexbin package (from Bioconductor)") } subset <- substitute(subset) subset <- with(design$variables, subset) if (length(subset) > 0) design <- design[subset, ] W <- weights(design, "sampling") mf <- model.frame(formula, design$variables, na.action = na.pass) Y <- model.response(mf) X <- mf[, attr(attr(mf, "terms"), "term.labels")] switch(style, bubble = { symbols(X, Y, circles = sqrt(W), inches = inches, ...) }, hex = { if (exists("hcell")) { rval <- hexbin(X, Y, xbins = xbins) cell <- hcell(X, Y)$cell rval$cnt <- tapply(W, cell, sum) rval$xcm <- tapply(1:length(X), cell, function(ii) weighted.mean(X[ii], W[ii])) rval$ycm <- tapply(1:length(Y), cell, function(ii) weighted.mean(Y[ii], W[ii])) plot(rval, legend = legend, style = "centroids", ...) } else { rval <- hexbin(X, Y, ID = TRUE, xbins = xbins) cell <- rval@cID rval@count <- as.vector(tapply(W, cell, sum)) rval@xcm <- as.vector(tapply(1:length(X), cell, function(ii) weighted.mean(X[ii], W[ii]))) rval@ycm <- as.vector(tapply(1:length(Y), cell, function(ii) weighted.mean(Y[ii], W[ii]))) gplot.hexbin(rval, legend = legend, style = "centroids", ...) } }, grayhex = { if (exists("hcell")) { rval <- hexbin(X, Y, xbins = xbins) cell <- hcell(X, Y)$cell rval$cnt <- tapply(W, cell, sum) plot(rval, legend = legend, ...) } else { rval <- hexbin(X, Y, ID = TRUE, xbins = xbins) cell <- rval@cID rval@count <- as.vector(tapply(W, cell, sum)) gplot.hexbin(rval, legend = legend, ...) } }, subsample = { index <- sample(length(X), sample.size, replace = TRUE, prob = W) if (is.numeric(X)) xs <- jitter(X[index], factor = 3, amount = amount$x) else xs <- X[index] if (is.numeric(Y)) ys <- jitter(Y[index], factor = 3, amount = amount$y) else ys <- Y[index] plot(xs, ys, ...) }, transparent = { transcol <- function(base, opacity) { rgbs <- col2rgb(base)/255 rgb(rgbs[1, ], rgbs[2, ], rgbs[3, ], alpha = opacity) } if (is.function(basecol)) basecol <- basecol(model.frame(design)) w <- weights(design) maxw <- max(w) minw <- 0 alphas <- (alpha[1] * (maxw - w) + alpha[2] * (w - minw))/(maxw - minw) plot(X, Y, col = transcol(basecol, alphas), ...) }) } d$w.t<-w.t svystuff<-svydesign(id=~pid, data=d, weights=~w.t) postscript("figure2a.eps",height=6,width=6,horiz=FALSE) hist(ds$x,breaks=seq(200,750, by=25),axes=FALSE,ylab="Frequency",xlab="CD4 Rule",cex.lab=1.3,main="") axis(1) axis(2,at=c(0,2000,4000,6000,8000,10000),label=c("0","","4000","","8000","")) mtext("A. Histogram of HAART Initiation Rules",cex=1.5,line=1.5) dev.off() postscript("figure2b.eps",height=6,width=6,horiz=FALSE) hist(w.t,nclass=20,xlab="weights",cex.axis=1,main="",cex.lab=1.3) mtext("B. Histogram of Inverse Probability Weights",cex=1.5,line=1.5) dev.off() postscript("figure2c.eps",height=6,width=6,horiz=FALSE) t.ynew<-function(x) { (x-12-12)*100/18 } hist(t.ynew(b$y.cd4.ade),nclass=20,xlab="CD4 count-based Utility",cex.axis=1,cex.lab=1.3,main="",axes=FALSE) axis(2) axis(1, at=c(0,200,400,600,800,1000,1200),c(0,200,400,600,800,1000,1200)) mtext("C. Histogram of Utility 1",cex=1.5,line=1.5) dev.off() postscript("figure2e.eps",height=6,width=6,horiz=FALSE) hist(bs$y.complete.ade,nclass=20,xlab="Quality of Life Utility",cex.axis=1,xlim=c(0,1),cex.lab=1.3,main="") mtext("E. Histogram of Utility 2",cex=1.5,line=1.5) dev.off() postscript("figure2d.eps",height=6,width=6,horiz=FALSE) pstuff<-stuff$x.xbeta the.plot<-svyplot(t.ynew(d$y.cd4.ade)~d$x, design=svystuff ,style="grayhex",legend=0,xlab="CD4 Rule", ylab="CD4 count-based Utility",xbins=50,yaxt="n", main=paste("D. Utility 1 maximized at CD4 Rule = ", round(max.x.1),sep="") ) grid.lines(x=unit(stuff$x$x,"native"),y=unit(t.ynew(stuff$x$y.cd4.ade),"native"), gp=gpar(col="white", lwd=2),vp=the.plot$plot.vp@hexVp.on) grid.points(x=unit(max.x.1,"native"),y=unit(t.ynew(max(pred1,na.rm=TRUE)),"native"), gp=gpar(cex=2,lwd=2,col="white"),vp=the.plot$plot.vp@hexVp.on) grid.yaxis(label=FALSE, at=c(0,200,400,600,800,1000,1200),vp=the.plot$plot.vp@hexVp.off) grid.text(c(0,200,400,600,800,1000,1200), y= unit(c(0,200,400,600,800,1000,1200), "native"), x=unit(-1, "lines"), gp = gpar(fontsize = 12), rot = 90,vp=the.plot$plot.vp@hexVp.off) dev.off() postscript("figure2f.eps",height=6,width=6,horiz=FALSE) pstuff<-stuff2$x.xbeta the.plot2<-svyplot(y.complete.ade~x, design=svystuff ,style="grayhex",legend=0,xlab="CD4 Rule", ylab="Quality of Life Utility",xbins=50,yaxt="n", main=paste("F. Utility 2 maximized at CD4 Rule = ", round(max.x.2),sep="") ) grid.lines(x=unit(stuff2$x$x,"native"),y=unit(stuff2$x$y.complete.ade,"native"), gp=gpar(col="black", lwd=2),vp=the.plot2$plot.vp@hexVp.on) grid.points(x=unit(max.x.2,"native"),y=unit(max(pred2,na.rm=TRUE),"native"), gp=gpar(cex=2,lwd=2,col="black"),vp=the.plot2$plot.vp@hexVp.on) grid.yaxis(label=FALSE, at=c(0,.2,.4,.6,.8,1),vp=the.plot2$plot.vp@hexVp.off) grid.text(c(0,.2,.4,.6,.8,1), y= unit(c(0,.2,.4,.6,.8,1), "native"), x=unit(-1, "lines"), gp = gpar(fontsize = 12), rot = 90,vp=the.plot2$plot.vp@hexVp.off) dev.off()