better.image <- function(x, y, z, zlim = range(z[is.finite(z)]), xlim = range(x), ylim = range(y), col = heat.colors(12), xaxs = "i", yaxs = "i", xlab, ylab, breaks, ...){ contains <- function(point1, segments1, point2, segments2) { if(sum(unlist(lapply(segments2, function(seg) !is.null(seg) && seg[2] == seg[4] && seg[2] <= point1[2] && seg[1] <= point1[1] && seg[3] >= point1[1]))) %% 2L) { return(1) } if(sum(unlist(lapply(segments1, function(seg) !is.null(seg) && seg[2] == seg[4] && seg[2] <= point2[2] && seg[1] <= point2[1] && seg[3] >= point2[1]))) %% 2L) { return(-1) } 0 } qsort <- function(list) { if(length(list) <= 1) return(list) pivot <- list[[1]] less <- NULL greater <- NULL center <- list[1] for(val in list[-1]) { resp <- contains(pivot\$path[1,], pivot\$segments, val\$path[1,], val\$segments) if(resp == 1) { less[[length(less)+1]] <- val } else if(resp == -1 || resp == 0) { greater[[length(greater)+1]] <- val } else { center[[length(center)+1]] <- val } } return(c(Recall(less), center, Recall(greater))) } get.index <- function(i, upper) ifelse(i < 1L | i > upper, NA, i) indx <- function(mat, modx=0L, mody=0L) get.index(row(mat) + modx, nrow(mat)) + (get.index(col(mat) + mody, ncol(mat)) -1L)*nrow(mat) dx <- 0.5 * diff(x) x <- c(x[1L] - dx[1L], x[-length(x)] + dx, x[length(x)] + dx[length(x) - 1]) dy <- 0.5 * diff(y) y <- c(y[1L] - dy[1L], y[-length(y)] + dy, y[length(y)] + dy[length(y) - 1]) edge <- matrix(ncol=ncol(z)*2L, nrow=nrow(z)*2L) row.indx <- (seq_len(nrow(edge))+2L) %/% 2L col.indx <- (seq_len(ncol(edge))+2L) %/% 2L stuff.indx <- (row(edge) +1L) %/% 2L + (col(edge) -1L) %/% 2L * nrow(z) plot(NA, NA, xlim = xlim, ylim = ylim, type = "n", xaxs = xaxs, yaxs = yaxs, xlab = xlab, ylab = ylab, ...) zi <- .C("bincode", as.double(z), length(z), as.double(breaks), length(breaks), code = integer(length(z)), (TRUE), (TRUE), nok = TRUE, NAOK = TRUE, DUP = FALSE, PACKAGE = "base")\$code poly <- NULL for(level in seq_len(length(breaks)-1)) { stuff <- !is.na(zi) & zi == level edge <- stuff[stuff.indx]*8L - ifelse(!is.na(indx(edge, -1L, +0L)) & stuff[stuff.indx[indx(edge, -1L, +0L)]], 1L, 0L) - ifelse(!is.na(indx(edge, +1L, +0L)) & stuff[stuff.indx[indx(edge, +1L, +0L)]], 1L, 0L) - ifelse(!is.na(indx(edge, +0L, -1L)) & stuff[stuff.indx[indx(edge, +0L, -1L)]], 1L, 0L) - ifelse(!is.na(indx(edge, +0L, +1L)) & stuff[stuff.indx[indx(edge, +0L, +1L)]], 1L, 0L) - ifelse(!is.na(indx(edge, -1L, -1L)) & stuff[stuff.indx[indx(edge, -1L, -1L)]], 1L, 0L) - ifelse(!is.na(indx(edge, +1L, -1L)) & stuff[stuff.indx[indx(edge, +1L, -1L)]], 1L, 0L) - ifelse(!is.na(indx(edge, +1L, +1L)) & stuff[stuff.indx[indx(edge, +1L, +1L)]], 1L, 0L) - ifelse(!is.na(indx(edge, -1L, +1L)) & stuff[stuff.indx[indx(edge, -1L, +1L)]], 1L, 0L) > 0L seg.max <- sum(edge) polySegments <- vector("list", seg.max) s.count <- 0L for(i in seq_len(nrow(edge))) { for(j in seq_len(ncol(edge))) { if(edge[i,j]) { path <- data.frame(x=integer(seg.max),y=integer(seg.max)) path[1,] <- c(i,j) count <- 1L x.end <- x.start <- i y.end <- y.start <- j newi <- i newj <- j edge[newi,newj] <- FALSE while(TRUE) { if(newj < ncol(edge) && edge[newi, newj + 1L]) { newj <- newj + 1L }else if(newi < nrow(edge) && edge[newi + 1L, newj]) { newi <- newi + 1L }else if(newj > 1L && edge[newi, newj - 1L]) { newj <- newj - 1L }else if(newi > 1L && edge[newi - 1L, newj]) { newi <- newi - 1L }else{ x.end <- path[1L,1L] y.end <- path[1L,2L] break } if(x.start == newi) { x.end <- newi y.end <- newj } else if(y.start == newj) { x.end <- newi y.end <- newj } else { count <- count+1L polySegments[[s.count + count - 1L]] <- if(x.start < x.end) c(x.start,y.start,x.end,y.end) else c(x.end,y.end,x.start,y.start) path[count,] <- c(x.end,y.end) x.start <- x.end x.end <- newi y.start <- y.end y.end <- newj } edge[newi,newj] <- FALSE } polySegments[[s.count + count]] <- if(x.start < x.end) c(x.start,y.start,x.end,y.end) else c(x.end,y.end,x.start,y.start) path <- path[seq_len(count),] x.point <- path[1,1] y.point <- path[1,2] # numcrossings <- sum(unlist(lapply(polySegments[seq_len(s.count+count)], function(seg) !is.null(seg) && seg[2] == seg[4] && seg[2] <= y.point && seg[1] <= x.point && seg[3] >= x.point))) %% 2L numcrossings <- 1 containsLevel <- stuff[stuff.indx[x.point +1L, y.point +1L]] if(numcrossings && containsLevel) { poly <- c(poly, list(list(level=level, segments=polySegments[seq.int(from=s.count+1,length.out=count)], path=path))) s.count <- s.count + count } else { polySegments[seq.int(from=s.count+1, length.out=count)] <- NULL } } } } } # poly <- qsort(poly) lapply(poly, function(p) polygon(x=x[row.indx[p\$path\$x]],y=y[col.indx[p\$path\$y]], border=NA, col=col[p\$level])) poly } layout.mat <- matrix(c(1,3,5,7,12,1,3,5,7,12,1,3,5,8,12,1,3,5,8,12,1,3,5,9,12,2,4,6,9,12,2,4,6,10,12,2,4,6,10,12,2,4,6,11,12,2,4,6,11,12),nrow=5) widths <- c(1,1,1,1,1,1,1,1,1,1) heights <- c(5,5,5,1,4) legend <- list("Animated VaxGen sensitivity analysis of SCE(t), the difference in the probability of initiating antiretroviral therapy by the given time (placebo minus vaccine) among nonwhites who would have been infected regardless of treatment assignment. Each", expression(paste("slide shows a sensitivity analysis of SCE(t), varying ", phi, ", ",beta[0],", and ",beta[1], ". The time after")), "HIV-infection, t, varies from 1 to 24. Contours represent the estimated SCE(t) at a given level of sensitivity parameters. Shaded regions correspond to those sensitivity parameters where the Wald-based 95% confidence intervals for SCE(t) do not contain 0. Estimates in the dark-shaded regions imply that among those who would have been infected regardless of treatment assignment, vaccination lowered the probability of starting antiretrovirals by t months of infection diagnosis (was beneficial), whereas estimates in the light-shaded regions imply that vaccination increased the probability of starting antiretrovirals within t months (was harmful).") mar1 <- c(0,0,1.1,0) mar2 <- c(2.6,2.6,1.1,1.1) mgp <- c(1.5, 0.5, 0) W <- 4.5 H <- W * sum(heights)/sum(widths) #W <- (H + ((-mar1[3]-mar1[1]) + 3*(-mar2[3]-mar2[1]))*0.2)*2/3 + ((mar1[4]+mar1[2])+2*(mar2[4]+mar2[2]))*0.2 val <- (sqrt(12)-2.5)/sqrt(12) rev <- data.frame(x=c(.5-2/3, .5-2/3,.5, .5,3, 3,5.5,5.5,3,3, .5,.5), y=c(2,-2,-2,-val, -2,-val, -2, 2,val, 2,val, 2)) rev\$x <- rev\$x - sum(range(rev\$x))/2 step <- data.frame(x=c(-2*sqrt(3),-0,-0,NA,2/3, 4/3, 4/3,2/3), y=c( 0,-2, 2,NA, -2,-2,2, 2)) step\$x <- step\$x - sum(range(step\$x,na.rm=TRUE))/2 pause <- data.frame(x=c( 0, 0, 2/3, 2/3, NA, 4/3, 4/3, 2, 2), y=c( 2,-2, -2, 2, NA, 2,-2, -2, 2)) pause\$x <- pause\$x - sum(range(pause\$x,na.rm=TRUE))/2 play <- data.frame(x=c(0,0,2*sqrt(3)), y=c(2,-2,0)) play\$x <- play\$x - sum(range(play\$x,na.rm=TRUE))/2 symbols <- list(list(poly=rev, pos=7, id="rewind"), list(poly=step, pos=8, id="rstep"), list(poly=pause, pos=9, id="pause"), list(poly=play, pos=9, id="play"), list(poly=-step, pos=10, id="fstep"), list(poly=-rev, pos=11, id="fforward")) ## x11(height=H, width=W) ## opar <- par(no.readonly=TRUE) ## par(oma=mar1) ## layout(layout.mat, widths=widths, heights=heights, respect=TRUE) ## par(mar=mar2, cex=1, mgp=mgp, cex.main=1) ## i <- 14 ## for(j in seq_along(phi)) { ## # plot.new() ## # plot.window(xlim=range(beta0),ylim=range(beta1)) ## pol<-better.image(beta0,beta1, reject1[i,,,j], breaks=c(-1.5,-0.5,0.5,1.5), ## col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="", axes=FALSE, asp=1, ann=FALSE, ## xlim=range(beta0), ylim=range(beta1)) ## # better.image(beta0,beta1, reject1[i,,,j], breaks=c(-1.5,-0.5,0.5,1.5), ## # col=c(gray(.9),gray(0),gray(.8)),xlab="",ylab="", axes=FALSE, asp=1, ann=FALSE, ## # xlim=range(beta0), ylim=range(beta1)) ## contour(beta0,beta1,ans\$SCE[i,,,j],add=T,lty=3,levels=c((0:20)/10-1),labcex=0.8, axes=FALSE) ## axis(side=1,at=pretty(beta0),label=format(pretty(beta0), nsmall=2,digits=2),cex.axis=0.792, lwd=0, lwd.tick=1) ## axis(side=2,at=pretty(beta1),label=format(pretty(beta1), nsmall=2,digits=2),cex.axis=0.792, lwd=0, lwd.tick=1) ## title(main=bquote(phi==.(format(phi[j], digits=3))), xlab=expression(beta[0]), ylab=expression(beta[1])) ## # mtext(side=1, line=2.4, expression(beta[0])) ## # mtext(side=2, line=2.4, expression(beta[1])) ## # mtext(side=3, line=0, bquote(phi==.(format(phi[j], digits=3)))) ## } ## mtext(side=3, line=0.1, outer=TRUE, ## bquote(Time==.(paste(t.point[i], " Months")))) ## par(mar=c(.25,.25,.25,.25)) ## for(sym in symbols[-4]) { ## plot.new() ## plot.window(xlim=c(-6,6), ylim=c(-3,3), asp=1) ## par(usr=c(-6,6,-3,3)) ## rect(-6,-3,6,3, col="grey") ## polygon(sym\$poly, col='black') ## } ## #dev.off() svg(file='vaxgen_axes_boxes.svg', height=H, width=W) par(oma=mar1) layout(ifelse(layout.mat == 12,7, ifelse(layout.mat%in%c(7:11), 0, layout.mat)), heights=heights, widths=widths, respect=TRUE) par(mar=mar2, cex=1, mgp=mgp, cex.main=1) axesValues <- c(0.05,.2,1,5,20) for(j in seq_along(phi)) { plot.new() plot.window(xlim=range(beta0),ylim=range(beta1), xaxs='i', yaxs='i', asp=1) par(cex.axis=0.6) box('plot') axis(side=1,at=log(axesValues)/12,label=format(axesValues, digits=1,trim=TRUE,drop0trailing=TRUE), line=NA) axis(side=2,at=log(axesValues)/12,label=format(axesValues, digits=1,trim=TRUE,drop0trailing=TRUE), line=NA) title(main=bquote(phi==.(format(phi[j], digits=3))), xlab=expression(OR[0] == exp(beta[0])), ylab=expression(OR[1] == exp(beta[1]))) } par(mar=c(.5,.5,.5,.5)) plot.new() plot.window(xlim=c(-60L,60L), ylim=c(-120L,120L)) par(cex=0.66, usr=c(-60L, 60L, -120L, 120L)) dy <- cumsum(c(118, - 3*par("cxy")[2]*par('cex'), -par("cxy")[2]*par('cex'))) text(-59, dy[1], legend[[1]], adj=c(0,1)) text(-59, dy[2], legend[[2]], adj=c(0,1)) text(-59, dy[3], legend[[3]], adj=c(0,1)) dev.off() for(sym in symbols) { svg(file=paste('vaxgen_button_',sym\$id,'.svg',sep=''), height=H,width=W) par(oma=mar1) layout(ifelse(layout.mat==sym\$pos,1,0), heights=heights, widths=widths, respect=TRUE) par(mar=mar2,cex=1, mgp=c(2.4,1,0), cex.main=1) # replicate(6, { plot.new(); plot.window(xlim=range(beta0),ylim=range(beta1));points(NA,NA) }) par(mar=c(.25,.25,.25,.25)) plot.new() plot.window(xlim=c(-6,6), ylim=c(-3,3), asp=1) par(usr=c(-6,6,-3,3)) rect(-6,-3,6,3, col="grey") polygon(sym\$poly, col='black') dev.off() } load('vaxgen_dataset.rda') vara1 <- sqrt(ans\$SCE.var) reject1 <- ifelse(ans\$SCE < qnorm(0.025)*vara1, -1, ifelse(ans\$SCE > qnorm(0.975)*vara1, 1, 0)) svg(file='vaxgen_animframe%03d.svg', height=H, width=W, onefile=FALSE) #x11(height=H,width=W) par(oma=mar1) layout(ifelse(layout.mat <=6, layout.mat, 0), widths=widths, heights=heights, respect=TRUE) par(mar=mar2, cex=1, mgp=mgp, cex.main=1) for(i in seq(length=dim(reject1)[1])) { for(j in seq_len(dim(reject1)[4])) { # image(beta0,beta1, reject1[i,,,j], breaks=c(-1.5,-0.5,0.5,1.5), # col=c(gray(.9),gray(1),gray(.8)),xlab="",ylab="", axes=FALSE, asp=1) pol<-better.image(beta0,beta1, reject1[i,,,j], breaks=c(-1.5,-0.5,0.5,1.5), col=c(gray(.85),gray(1),gray(.6)),xlab="",ylab="", axes=FALSE, asp=1, ann=FALSE, xlim=range(beta0), ylim=range(beta1)) contour(beta0,beta1,ans\$SCE[i,,,j],add=T,lty=3,levels=c((0:20)/10-1),labcex=0.8, axes=FALSE) cat(".") } mtext(side=3, line=0, outer=TRUE, bquote(Time==.(paste(t.point[i], " Months")))) cat("*") } cat('\n') dev.off()