########################################################################################### # Program Name: Chapter 9_Figures_Individual_PDF.R # # Purpose: Generating individual graphics in PDF format from # Chapter 9. Statistical Graphics in Clinical Oncology in the book # A Picture is Worth a Thousand Tables: Graphics in Life Sciences # # Author: Kye Gilder, PhD # NuVasive, Inc. # 7475 Lusk Blvd. # San Diego, CA 92121 # kgilder@nuvasive.com # # Version: 08AUG2012 # # R Version: 2.14.1 (2011-12-22) # # Notes: Comments and suggestions are welcomed. ########################################################################################### ########################################################################################### # Set-up ########################################################################################### rm(list=ls()) #PDF output path path <- "C:/Documents and Settings/kgilder/My Documents/ ... /Graphics/" ########################################################################################### # Packages ########################################################################################### library('rms') #survplot() function library('Hmisc') #event.chart(), errbar(), Ecdf() function library('survival') #ovarian and lung data library('PBSmodelling') #resetGraph() function library('graphics') #density() function library('lattice') library('grid') library('RColorBrewer') ########################################################################################### # Custom Functions ########################################################################################### capwords <- function(s, strict = FALSE) { cap <- function(s) paste(toupper(substring(s,1,1)), {s <- substring(s,2); if(strict) tolower(s) else s}, sep = "", collapse = " " ) sapply(strsplit(s, split = " "), cap, USE.NAMES = !is.null(names(s))) } ########################################################################################### # Custom Colors (RColorBrewer) ########################################################################################### mycolors <-brewer.pal(11,"PRGn")[c(2,10)] mycolors2 <-brewer.pal(4,"Dark2") mycolors3 <-c(brewer.pal(4,"Dark2")[2],brewer.pal(4,"Dark2")[1], brewer.pal(4,"Dark2")[3], brewer.pal(4,"Dark2")[4]) myredblue <-brewer.pal(11,"RdBu")[c(2,10)] myorred <-brewer.pal(9,"OrRd")[c(5:9)] ########################################################################################### # Figure 1 # Dose-Toxicity Profile (e.g., logistic, hyperbolic tangent, Lymann, Emax) ########################################################################################### p.tox <- function(a,b,x){round((1/(1+exp(-a-b*x))),3)} dose.starting <- 20 dose.actual <- c(20, 80, 140, 210, 280, 370, 500, 670, 900) doses <- dose.actual/dose.starting doses.label <- dose.starting*doses Scenario1 <- c(p.tox(-3.4, 0.4010, doses)[1:5], 0.974, 0.999, 0.999, 0.999) Scenario2 <- c(p.tox(-6.0, 0.2870, doses)[1:8], 0.999) Scenario3 <- c(p.tox(-3.0, 0.1247, doses)[1:8], 0.999) Scenario4 <- c(p.tox(-5.0, 0.1723, doses)[1:8], 0.999) Scenario5 <- c(p.tox(-6.0, 0.1585, doses)[1:8], 0.999) resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_01.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(1,1), mar=c(4.2,4.0,1.0,1.0)) plot(doses, seq(0, 1, length=length(doses)), type="n", ylim=c(0,1), axes=FALSE, xlab="Dose (mg)", ylab="Probability of Toxicity", main="") abline(h=1/3, col=8, lwd=2) points(doses, Scenario5, type="b", col=myorred[1], lwd=2.5, pch=16) #Scenario 5 points(doses, Scenario4, type="b", col=myorred[2], lwd=2.5, pch=16) #Scenario 4 points(doses, Scenario3, type="b", col=myorred[3], lwd=2.5, pch=16) #Scenario 3 points(doses, Scenario2, type="b", col=myorred[4], lwd=2.5, pch=16) #Scenario 2 points(doses, Scenario1, type="b", col=myorred[5], lwd=2.5, pch=16) #Scenario 1 axis(1, at=doses, labels=doses.label) axis(2, at=seq(0,1,0.1), tck=0.02, labels=FALSE, tick=TRUE) text(-2.5, seq(0,1,0.1), labels=as.character(seq(0,1,0.1)), srt=0, xpd=TRUE) axis(4, at=seq(0,1,0.1), tck=0.02, labels=FALSE, tick=TRUE) box() dev.off() ########################################################################################### # Figure 2 # Operating Characteristics of a Phase 1 'A+B' Dose Escallation Design ########################################################################################### #Dose-toxicity function p.tox <- function(a,b,x){round((1/(1+exp(-a-b*x))),3)} # Load two fuctions # (1) pmtd() is function for Operating Characteristics of A+B # (2) plot.pmtd() is for plotting of Operating Characteristics of A+B # Contact author for R program pmtd_Rbook.R source("C:\\Documents and Settings\\kgilder\\ ... \\pmtd_Rbook.r") dose.starting <- 20 dose.actual <- c(20, 80, 140, 210, 280, 370, 500, 670, 900) doses <- dose.actual/dose.starting doses.label <- dose.starting*doses Scenario <- c(p.tox(-6.0, 0.2870, doses)[1:8], 0.999) pmtd(doses.label, Scenario, design=c(3,3,1,1,1), deescalate=1) resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_02.pdf", sep=""), bg="white", onefile=TRUE) plot.pmtd(pmtd(doses.label, Scenario, design = c(3,3,1,1,1), deescalate=1), title="Scenario 1: '3+3' with De-escallation")#dev.off() dev.off() ########################################################################################### # Figure 3 # Operating Characteristics of a Phase 2 Predicted Probability Design vs. Simon 2-Stage 'Optimal' Design ########################################################################################### Simon <- data.frame(n=c(17,37), r=c(3,10)) PPD <- data.frame(n=c(10, 17, 21, 24, 27, 29, 31, 33, 34, 35, 36), r=c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)) resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_03.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(1,1), mar=c(6,4.0,1.0,1.0)) plot(10:38, seq(0,12,len=29), type="n", ylim=c(0,12), axes=FALSE, xlab="Patient", ylab="Rejection Region (Responders)") axis(1, at=10:38) axis(2, at=0:12, labels=FALSE) text(7.5, seq(0,12,1), as.character(seq(0,12,1)), srt=0, xpd=TRUE) axis(4, at=0:12, tck=0.02, labels=FALSE) segments(Simon$n,rep(0,2),Simon$n,Simon$r, col=myredblue[1], lwd=3) points(Simon$n,Simon$r, pch=16, col=myredblue[1], lwd=3) segments(PPD$n,rep(0,11),PPD$n,PPD$r, col=myredblue[2], lwd=3) points(PPD$n,PPD$r, pch=16, col=myredblue[2], lwd=3) segments(38,0,38,11, col="black", lwd=3) points(38,11, pch=16, col="black", lwd=3) points(9.5,9.5, pch=16, col="black", lwd=3) points(9.5,9.0, pch=16, col=myredblue[1], lwd=3) points(9.5,8.5, pch=16, col=myredblue[2], lwd=3) text(10, 12.0, substitute(paste(p[0]," = 0.20 (Null)")), adj=0) text(10, 11.5, substitute(paste(p[1]," = 0.40 (Alternative)")), adj=0) text(10, 11.0, "Type-1 error = 0.10", adj=0) text(10, 10.5, "Power = 90%", adj=0) #text(10, 10.0, "Max N = 37", adj=0) text(10, 09.5, substitute(paste("E[N | ", p[0],"] = 38 (Fixed sample size design)")) , adj=0) text(10, 09.0, substitute(paste("E[N | ", p[0],"] = 26.0 (Simon's optimal 2-stage design)")) , adj=0) text(10, 08.5, substitute(paste("E[N | ", p[0],"] = 27.7 (Predicted Probability design for maximize power)")) , adj=0) #text(15, 5., "E[N | p0] = 32 (PPD minimize N under p0 design)", adj=0, col=myredblue[2]) mtext(substitute(paste("Note: If number of ", responders <= Rejection, "Region, trial is terminated due to insufficient evidence of efficacy.")), side=1, line=4.25, at=7.5, cex=0.75, adj=0) box() dev.off() ########################################################################################### # Figure 4 # Triangle Test ########################################################################################### # Load function for reading SAS data files (i.e., .sas7bdat) #source("http://biostatmatt.com/R/sas7bdat.R") #Load simulation data #sim <- read.sas7bdat("C:/Documents and Settings/kgilder/My Documents/Working (Laptop)/Software/SAS/Programs/Clinical Programs/Triangle Test/simres.sas7bdat") #boundary <- read.sas7bdat("C:/Documents and Settings/kgilder/My Documents/Working (Laptop)/Software/SAS/Programs/Clinical Programs/Triangle Test/boundary.sas7bdat") #design <- read.sas7bdat("C:/Documents and Settings/kgilder/My Documents/Working (Laptop)/Software/SAS/Programs/Clinical Programs/Triangle Test/design.sas7bdat") sim <- read.csv("C:/Documents and Settings/kgilder/My Documents/Working (Laptop)/Software/SAS/Programs/Clinical Programs/Triangle Test/simres.csv", header=T, sep=",") boundary <- read.csv("C:/Documents and Settings/kgilder/My Documents/Working (Laptop)/Software/SAS/Programs/Clinical Programs/Triangle Test/boundary.csv", header=T, sep=",") design <- read.csv("C:/Documents and Settings/kgilder/My Documents/Working (Laptop)/Software/SAS/Programs/Clinical Programs/Triangle Test/design.csv", header=T, sep=",") attach(sim) attach(boundary) attach(design) names(sim) names(boundary) names(design) resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_04.pdf", sep=""), bg="white", onefile=TRUE) par(fig=c(0, 1, 0.25, 1), mar=c(4.0, 4.2, 1.0, 1.0), xpd=TRUE) plot(0:15, seq(-4,10,len=16), type="n", ylab="Z (Efficient Score)", xlab="", xlim=c(0,14.5), ylim=c(-4,10), axes=FALSE, xaxs="i") lines(x=c(7.65, 7.65), y=c(-6, 10.5), lty=4, col="black") segments(0, 0, VMAX, 0, col="darkgrey") segments(0, INTERCPT, VMAX, APEXZ, col="blue", lwd=3) segments(0, -INTERCPT, VMAX, APEXZ, col="red", lwd=3, lty=3) newy <- ((APEXZ-(-INTERCPT))/VMAX)*(VMAX/6)-INTERCPT segments((VMAX/6), newy, VMAX, APEXZ, col="red", lwd=3) points(VFINAL, ZFINAL, type="p", pch=16, cex=0.5) #axis(1, at=VFINAL, labels=as.character(TN)) axis(1, at=seq(0.9,14.4,.9), labels=FALSE, tick=TRUE) text(seq(0.9,14.4,.9),-5.3, labels=as.character(seq(10,160,10)), cex=0.9) axis(2, at=seq(-4,10,1), labels=FALSE, tick=TRUE) text(-0.6, seq(-4,10,1), labels=as.character(seq(-4,10,1)), srt=0) text(6.0, 10, "Fixed Sample Size", cex=0.75) text(9.0, -0.5, "Simulation Results", cex=0.90, adj=0) text(9.3, -1.0, paste("Simulation runs =", dim(sim)[1]), cex=0.75, adj=0) text(9.3, -1.5, paste("P(cross upper boundary) = ", round(table(bcrossed)[1]/dim(sim)[1],3)), cex=0.75, adj=0) text(9.3, -2.0, paste("P(cross lower boundary) = ", round(table(bcrossed)[2]/dim(sim)[1],3)), cex=0.75, adj=0) text(9.3, -2.5, paste("Mean sample size = ", round(mean(TN),2)), cex=0.75, adj=0) text(9.3, -3.0, paste("Median sample size = ", median(TN)), cex=0.75, adj=0) text(0.1, 3.6, "Superiority", cex=0.85, col="blue", adj=0) text(0.1, -2.4, "Inferiority", cex=0.85, col="red", adj=0) text(3.0, -0.5, "Futility", cex=0.75, col="red") box() par(fig=c(0, 1, 0, 0.3), new=TRUE) plot.density(density(sim[bcrossed==-1,]$VFINAL, kernel="gaussian"), xlab="", ylab="", main="", zero.line=FALSE, xlim=c(0,14.5), ylim=c(0,0.4), axes=FALSE, xaxs="i" ) lines(x=c(7.65, 7.65), y=c(-0.017, 1.37), lty=4, col="black") #axis(1, at=VFINAL, labels=as.character(TN)) axis(1, at=seq(0.9,14.4,.9), labels=FALSE, tick=TRUE) text(seq(0.9,14.4,.9), -0.1, labels=as.character(seq(10,160,10)), cex=0.9) axis(2, at=seq(0,0.4,0.1), labels=FALSE, tick=TRUE) text(-0.6, seq(0,0.4,0.1), labels=as.character(seq(0,0.4,0.1)), srt=0) mtext("Sample Size", side=1, line=-2, outer=TRUE) text(11, 0.27, "Distribution of Efficacy Boundary Crossing", cex=0.80) box() dev.off() ########################################################################################### # Figure 5 # 3-in-1 plot ########################################################################################### PFS <- read.csv("C:\\Documents and Settings\\kgilder\\My Documents\\Working (Laptop)\\Studies_BIIB\\Program.206.Volociximab\\Study.206-OC-202\\IVRS\\Bayesian Adaptation\\d_PFS.csv", header=T, sep=",") names(PFS) attach(PFS) plot3.fun1 <- function(Dmean, Dsd, X1mean, X1sd, X2mean, X2sd, old, new, daysfu, event, group, sensitivity){ ### Required inputs: ### 1. LD mean PFS ### 2. LD SD PFS ### 3. LD+q2w mean PFS ### 4. LD+q2w SD PFS ### 5. LD+qw mean PFS ### 6. LD+qw SD PFS ### 7. Vector, length 3, of previous N per group ### 8. Vector, length 3, of new N per group ### 9. vector, length is # of patients, of days follow-up per patient ### 10. vector, length is # of patients, 1 for progression, 0 for no progression ### 11. vector, length is # of patients, 1=Doxil, 2=+q2w, 3=+qw library(survival) N <- sum(old) kmfit <- survfit(Surv(exposure,censor)~group) #Log-Rank LR <- survdiff(Surv(exposure,censor)~group)$chisq pvalue <- (1-pchisq(LR,2)) #Stratified Log-Rank LR.S <- coxph(Surv(exposure, censor)~group + strata(sensitivity), method="breslow")$score pvalue.S <- (1-pchisq(LR.S,2)) par(mfrow=c(1,1), mar=c(4.0,4.0,2.0,1.0)) layout(matrix(c(1,2,1,3), 2, 2, byrow = TRUE), widths=c(3,2), heights=c(1,1)) #Kaplan-Meier Plot plot(kmfit, xlab="", ylab="", xlim=c(0,500), axes=F, lwd=2, col=mycolors2[1:3], main="Kaplan-Meier Plot") axis(1, at=seq(0,500,100), cex=0.75) mtext(side=1, line=2.2, "Progression-Free Survival (Days)", cex=0.85) axis(2, at=seq(0,1,0.1), labels=FALSE) text(-35, seq(0,1,0.1), as.character(seq(0,1,0.1)), srt=0, xpd=TRUE) mtext(side=2, line=2.5, "Probability", cex=0.85) #text(450, 0.80, paste("Log-rank statistic =", round(LR,dig=3)), cex=0.75, col="black") #text(450, 0.77, paste("p-value =", round(pvalue,dig=3)), cex=0.70, col="black") #text(450, 0.70, paste("Stratified Log-rank statistic =", round(LR.S,dig=3)), cex=0.75, col="black") #text(450, 0.67, paste("p-value =", round(pvalue.S,dig=4)), cex=0.70, col="black") legend("topright", border="black", title="Treatment", legend=c(paste("A (n=", old[1], ")", sep=""), paste("B (n=", old[2], ")", sep=""), paste("C (n=", old[3], ")", sep="")), lty=1, lwd=2, col=c(col=mycolors2[1], col=mycolors2[2], col=mycolors2[3]), inset = .05) box() #Posterior Distribution Plot x <- seq(0.0001, 800, length=1000000) B1 <- Dmean / Dsd^2 B2 <- X1mean / X1sd^2 B3 <- X2mean / X2sd^2 A1 <- Dmean^2 / Dsd^2 A2 <- X1mean^2 / X1sd^2 A3 <- X2mean^2 / X2sd^2 mx <- max(c(dgamma(x,A1,B1),dgamma(x,A2,B2),dgamma(x,A3,B3)))*1.1 plot(x,rep(0,length(x)), type="n", xlab="", ylab="", ylim=c(0,mx), axes=F, main="Posterior Distributions") axis(1, at=seq(0,800,100), cex=0.55) mtext(side=1, line=2.2, "Mean PFS (Days)", cex=0.85) lines(x,dgamma(x,A1,B1), lwd=2, col=mycolors2[1]) lines(x,dgamma(x,A2,B2), lwd=2, col=mycolors2[2]) lines(x,dgamma(x,A3,B3), lwd=2, col=mycolors2[3]) #Randomization Count Bar Chart hw <- 0.4 plot(0,0, type="n", xlim=c(0,4), ylim=c(0,80), axes=F, xlab="", ylab="", main="Randomization Allocation") axis(1, at=c(1,2,3), labels=c("A","B","C"), cex=0.75) mtext(side=1, line=2.2, "Treatment", cex=0.85) axis(2, at=seq(0,80,10), labels=FALSE) text(-0.65, seq(0,80,10), as.character(seq(0,80,10)), srt=0, xpd=TRUE) mtext(side=2, line=2.1, "Total Patients", cex=0.85) polygon(c(1-hw,1+hw,1+hw,1-hw),c(0,0,old[1],old[1]), density=-1,col=mycolors2[1]) polygon(c(2-hw,2+hw,2+hw,2-hw),c(0,0,old[2],old[2]), density=-1,col=mycolors2[2]) polygon(c(3-hw,3+hw,3+hw,3-hw),c(0,0,old[3],old[3]), density=-1,col=mycolors2[3]) polygon(c(1-hw,1+hw,1+hw,1-hw),c(old[1],old[1],old[1]+new[1],old[1]+new[1]), density=25, col=mycolors2[1]) polygon(c(2-hw,2+hw,2+hw,2-hw),c(old[2],old[2],old[2]+new[2],old[2]+new[2]), density=25, col=mycolors2[2]) polygon(c(3-hw,3+hw,3+hw,3-hw),c(old[3],old[3],old[3]+new[3],old[3]+new[3]), density=25, col=mycolors2[3]) text(1, old[1]+new[1]+8, new[1], cex=0.85) text(2, old[2]+new[2]+8, new[2], cex=0.85) text(3, old[3]+new[3]+8, new[3], cex=0.85) box(bty="l") } plot3.fun2 <- function(Dmean, Dsd, X1mean, X1sd, X2mean, X2sd, old, new, daysfu, event, group, sensitivity){ ### Required inputs: ### 1. LD mean PFS ### 2. LD SD PFS ### 3. LD+q2w mean PFS ### 4. LD+q2w SD PFS ### 5. LD+qw mean PFS ### 6. LD+qw SD PFS ### 7. Vector, length 3, of previous N per group ### 8. Vector, length 3, of new N per group ### 9. vector, length is # of patients, of days follow-up per patient ### 10. vector, length is # of patients, 1 for progression, 0 for no progression ### 11. vector, length is # of patients, 1=Doxil, 2=+q2w, 3=+qw library(survival) N <- sum(old) kmfit <- survfit(Surv(exposure,censor)~group) #Log-Rank LR <- survdiff(Surv(exposure,censor)~group)$chisq pvalue <- (1-pchisq(LR,2)) #Stratified Log-Rank LR.S <- coxph(Surv(exposure, censor)~group + strata(sensitivity), method="breslow")$score pvalue.S <- (1-pchisq(LR.S,2)) par(mfrow=c(1,1), mar=c(4.0,4.0,2.0,1.0)) layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE), widths=c(1,1), heights=c(3,2)) #Kaplan-Meier Plot plot(kmfit, xlab="", ylab="", xlim=c(0,500), axes=F, lwd=2, col=mycolors2[1:3], main="Kaplan-Meier Plot") axis(1, at=seq(0,500,100), cex=0.75) mtext(side=1, line=2.2, "Progression-Free Survival (Days)", cex=0.85) axis(2, at=seq(0,1,0.1), labels=FALSE) text(-20, seq(0,1,0.1), as.character(seq(0,1,0.1)), srt=0, xpd=TRUE) mtext(side=2, line=2.5, "Probability", cex=0.85) #text(450, 0.80, paste("Log-rank statistic =", round(LR,dig=3)), cex=0.75, col="black") text(325, 0.12, paste("p-value (logrank test) =", round(pvalue,dig=3)), cex=0.90, col="black", adj=0) #text(450, 0.70, paste("Stratified Log-rank statistic =", round(LR.S,dig=3)), cex=0.75, col="black") #text(450, 0.67, paste("p-value =", round(pvalue.S,dig=4)), cex=0.70, col="black") legend("topright", border="black", title="Treatment", legend=c(paste("A (n=", old[1], ")", sep=""), paste("B (n=", old[2], ")", sep=""), paste("C (n=", old[3], ")", sep="")), lty=1, lwd=2, col=c(col=mycolors2[1], col=mycolors2[2], col=mycolors2[3]), inset = .05) box() #Posterior Distribution Plot x <- seq(0.0001, 800, length=1000000) B1 <- Dmean / Dsd^2 B2 <- X1mean / X1sd^2 B3 <- X2mean / X2sd^2 A1 <- Dmean^2 / Dsd^2 A2 <- X1mean^2 / X1sd^2 A3 <- X2mean^2 / X2sd^2 mx <- max(c(dgamma(x,A1,B1),dgamma(x,A2,B2),dgamma(x,A3,B3)))*1.1 plot(x,rep(0,length(x)), type="n", xlab="", ylab="", ylim=c(0,mx), axes=F, main="Posterior Distributions") axis(1, at=seq(0,800,100), cex=0.55) mtext(side=1, line=2.2, "Mean PFS (Days)", cex=0.85) lines(x,dgamma(x,A1,B1), lwd=2, col=mycolors2[1]) lines(x,dgamma(x,A2,B2), lwd=2, col=mycolors2[2]) lines(x,dgamma(x,A3,B3), lwd=2, col=mycolors2[3]) #Randomization Count Bar Chart hw <- 0.4 plot(0,0, type="n", xlim=c(0,4), ylim=c(0,80), axes=F, xlab="", ylab="", main="Randomization Allocation") axis(1, at=c(1,2,3), labels=c("A","B","C"), cex=0.75) mtext(side=1, line=2.2, "Treatment", cex=0.85) axis(2, at=seq(0,80,10), labels=FALSE) text(-0.5, seq(0,80,10), as.character(seq(0,80,10)), srt=0, xpd=TRUE) mtext(side=2, line=2.1, "Total Patients", cex=0.85) polygon(c(1-hw,1+hw,1+hw,1-hw),c(0,0,old[1],old[1]), density=-1,col=mycolors2[1]) polygon(c(2-hw,2+hw,2+hw,2-hw),c(0,0,old[2],old[2]), density=-1,col=mycolors2[2]) polygon(c(3-hw,3+hw,3+hw,3-hw),c(0,0,old[3],old[3]), density=-1,col=mycolors2[3]) polygon(c(1-hw,1+hw,1+hw,1-hw),c(old[1],old[1],old[1]+new[1],old[1]+new[1]), density=25, col=mycolors2[1]) polygon(c(2-hw,2+hw,2+hw,2-hw),c(old[2],old[2],old[2]+new[2],old[2]+new[2]), density=25, col=mycolors2[2]) polygon(c(3-hw,3+hw,3+hw,3-hw),c(old[3],old[3],old[3]+new[3],old[3]+new[3]), density=25, col=mycolors2[3]) text(1, old[1]+new[1]+8, new[1], cex=0.85) text(2, old[2]+new[2]+8, new[2], cex=0.85) text(3, old[3]+new[3]+8, new[3], cex=0.85) box(bty="l") } resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_05.pdf", sep=""), bg="white", onefile=TRUE) plot3.fun2(424.517, 104.959, 398.877, 119.521, 428.632, 144.287, c(67,34,27), c(4,3,2), PFS$exposure, PFS$censor, PFS$treat, PFS$platresgrp) dev.off() #resetGraph(reset.mf=TRUE) #pdf(paste(path, "Gilder_Figure_05_extra.pdf", sep=""), bg="white", onefile=TRUE) #plot3.fun1(424.517, 104.959, 398.877, 119.521, 428.632, 144.287, c(67,34,27), c(4,3,2), PFS$exposure, PFS$censor, PFS$treat, PFS$platresgrp) #dev.off() ########################################################################################### # Figure 6 & Figure 7 # Waterfall Plots ########################################################################################### set.seed(9) n <- 60 CFB <- data.frame(subj=1:n, lesion=runif(n,min=-80,max=40), treat=rep(c(0,1),n/2)) CFB$status <- ifelse(CFB$lesion>=20, 0, ifelse(CFB$lesion>=-30, 1, 2)) #CFB$status <- factor((CFB$lesion <= -30) + (CFB$lesion <= 20), lab = c(0, 1, 2)) ss <- table(CFB$treat) for(i in 1:n) { CFB$status.mod[i] <- ifelse(CFB$status[i]==0, 0, ifelse(CFB$status[i]==1, rbinom(1,1,0.7), 2*rbinom(1,1,0.8)))+1 } CFBA <- CFB[CFB$treat==0,] CFBB <- CFB[CFB$treat==1,] #Figure Waterfall 1 resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_06.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(1,1), mar=c(4.2,4.0,1.0,1.0)) plot(rev(CFB[order(CFB$lesion),]$lesion), type="n", ylim=c(-80,80), xlab="Patients", ylab="Change in SLD from Baseline (%)", axes=FALSE) abline(h=c(-30,20), lty=2, lwd=2, col="lightgrey") points(rev(CFB[order(CFB$lesion),]$lesion), type="h", col=mycolors3[rev(CFB[order(CFB$lesion),]$status.mod)], lwd=4) abline(h=0, lwd=3) axis(1, at=seq(0,n,10)) axis(1, at=seq(5,n,10), tck=-0.01, labels=FALSE, tick=TRUE) axis(2, at=seq(-80,800,20), tck=-0.02, labels=FALSE, tick=TRUE) text(-4, seq(-80,80,20), labels=as.character(seq(-80,80,20)), srt=0, xpd=TRUE) axis(4, at=seq(-80,80,20), tck=0.02, labels=FALSE) text(30, 24, "20% increase", col="black", cex=0.9) text(30, -34, "30% decrease", col="black", cex=0.9) legend("topright", border="black", legend=c("PD", "SD", "PR", "CR"), lwd=5, col=mycolors3, inset=.05) box() dev.off() #Figure Waterfall 2 #resetGraph(reset.mf=TRUE) #pdf(paste(path, "Gilder_Figure_07_extra.pdf", sep=""), bg="white", onefile=TRUE) #plot(rev(CFB[order(CFB$lesion),]$lesion), type="n", ylim=c(-80,80), xlab="Patients", ylab="Max change in SLD from Baseline (%)", axes=FALSE) #abline(h=c(-30,20), lty=2, lwd=2, col="lightgrey") #points(rev(CFB[order(CFB$lesion),]$lesion), type="h", lwd=4, col=ifelse(rev(CFB[order(CFB$lesion),]$treat)==0,mycolors[1],mycolors[2])) #abline(h=0, lwd=3) #axis(1, at=seq(0,n,10)) #axis(1, at=seq(5,n,10), tck=-0.01, labels=FALSE, tick=TRUE) #axis(2, at=seq(-80,80,20), tck=-0.02, labels=FALSE, tick=TRUE) #text(-4, seq(-80,80,20), labels=as.character(seq(-80,80,20)), srt=0, xpd=TRUE) #axis(4, at=seq(-80,80,20), tck=0.02, labels=FALSE) #text(30, 24, "20% increase", col="black", cex=0.9) #text(30, -34, "30% decrease", col="black", cex=0.9) #legend("topright", border="black", title="Treatment", legend=c(paste("A (n=", ss[1], ")", sep=""), paste("B (n=", ss[2], ")", sep="")), lty=1, lwd=5, col=c(col=mycolors[1], col=mycolors[2]), inset=0.05) #box() #dev.off() #Figure Waterfall 3 resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_07.pdf", sep=""), bg="white", onefile=TRUE) par(mar=c(4.2,4.0,1.0,1.0), fig=c(0,1,0.45,1)) plot(rev(CFB[order(CFB$lesion),]$lesion), type="n", ylim=c(-80,80), xlab="", ylab="Max change in SLD from Baseline (%)", axes=FALSE) abline(h=c(-30,20), lty=2, lwd=2, col="lightgrey") points(rev(CFB[order(CFB$lesion),]$lesion), type="h", lwd=4, col=ifelse(rev(CFB[order(CFB$lesion),]$treat)==0,mycolors[1],mycolors[2])) abline(h=0, lwd=3) axis(1, at=seq(0,n,10)) axis(1, at=seq(5,n,10), tck=-0.01, labels=FALSE, tick=TRUE) axis(2, at=seq(-80,80,20), tck=-0.02, labels=FALSE, tick=TRUE) text(-4, seq(-80,80,20), labels=as.character(seq(-80,80,20)), srt=0, xpd=TRUE) axis(4, at=seq(-80,80,20), tck=0.02, labels=FALSE) text(30, 26, "20% increase", col="black", cex=0.9) text(30, -36, "30% decrease", col="black", cex=0.9) legend("topright", border="black", title="Treatment", legend=c(paste("A (n=", ss[1], ")", sep=""), paste("B (n=", ss[2], ")", sep="")), lty=1, lwd=5, col=c(col=mycolors[1], col=mycolors[2]), inset=0.05) box() par(mar=c(4.2,4.0,1.0,1.0), fig=c(0,0.5,0,0.5), new=TRUE) plot(rev(CFBA[order(CFBA$lesion),]$lesion), type="n", ylim=c(-80,80), xlab="Patients", ylab="Max change in SLD from Baseline (%)", axes=FALSE) abline(h=c(-30,20), lty=2, lwd=2, col="lightgrey") points(rev(CFBA[order(CFBA$lesion),]$lesion), type="h", col=mycolors3[rev(CFBA[order(CFBA$lesion),]$status.mod)], lwd=4) abline(h=0, lwd=3) axis(1, at=seq(0,ss[1],5)) axis(1, at=seq(5,ss[1],5), tck=-0.01, labels=FALSE, tick=TRUE) axis(2, at=seq(-80,80,20), tck=-0.02, labels=FALSE, tick=TRUE) text(-3, seq(-80,80,20), labels=as.character(seq(-80,80,20)), srt=0, xpd=TRUE) axis(4, at=seq(-80,80,20), tck=0.02, labels=FALSE) legend("topright", border="black", legend=c("PD", "SD", "PR", "CR"), lwd=5, col=mycolors3, cex=0.5, inset=.05) text(10, 80, "Treatment A", adj=0, cex=0.9) box() par(mar=c(4.2,4.0,1.0,1.0), fig=c(0.5,1,0,0.5), new=TRUE) plot(rev(CFBB[order(CFBB$lesion),]$lesion), type="n", ylim=c(-80,80), xlab="Patients", ylab="", axes=FALSE) abline(h=c(-30,20), lty=2, lwd=2, col="lightgrey") points(rev(CFBB[order(CFBB$lesion),]$lesion), type="h", col=mycolors3[rev(CFBB[order(CFBB$lesion),]$status.mod)], lwd=4) abline(h=0, lwd=3) axis(1, at=seq(0,ss[2],5)) axis(1, at=seq(5,ss[2],5), tck=-0.01, labels=FALSE, tick=TRUE) axis(2, at=seq(-80,80,20), tck=-0.02, labels=FALSE, tick=TRUE) text(-3, seq(-80,80,20), labels=as.character(seq(-80,80,20)), srt=0, xpd=TRUE) axis(4, at=seq(-80,80,20), tck=0.02, labels=FALSE) legend("topright", border="black", legend=c("PD", "SD", "PR", "CR"), lwd=5, col=mycolors3, cex=0.5, inset=.05) text(10, 80, "Treatment B", adj=0, cex=0.9) box() dev.off() ########################################################################################### # Figure 8 & Figure 9 # Forrest Plot ########################################################################################### #Figure Forest 1 (meta-analysis) # Contact author for R program pmtd_Rbook.R source("C:\\Documents and Settings\\kgilder\\My Documents\\ ... \\ForrestPlot.R") #Figure Forest 2 (subgroup analysis) #? ovarian names(ovarian) ovarian$age.f <- as.factor(ifelse(ovarian$age= 57 yrs")) attach(ovarian) #table(rx) #table(ecog.ps) #table(resid.ds) #summary(age) #describe(age) #table(age.f) #by(futime, rx, mean) #by(ovarian[ovarian$ecog.ps==1,]$futime, ovarian[ovarian$ecog.ps==1,]$rx, mean) #by(ovarian[ovarian$ecog.ps==2,]$futime, ovarian[ovarian$ecog.ps==2,]$rx, mean) #Figure Forest Plot (subgroup analysis) cphmfit <- coxph(Surv(futime,fustat)~rx, method="efron", data=ovarian) cphmfit summary(cphmfit) title1 <- paste("Primary Analysis (N = ", round(cphmfit$n,0),")", sep = "") result1 <- as.numeric(c(1, round(cphmfit$n,0), summary(cphmfit)$conf.int[c(1,3,4)])) cphmfit <- coxph(Surv(futime,fustat)~rx, method="efron", data=ovarian, subset=ecog.ps==1) title2 <- paste("ECOG 1 (n = ", round(cphmfit$n,0),")", sep = "") result2 <- c(2, round(cphmfit$n,0), summary(cphmfit)$conf.int[c(1,3,4)]) cphmfit <- coxph(Surv(futime,fustat)~rx, method="efron", data=ovarian, subset=ecog.ps==2) title3 <- paste("ECOG 2 (n = ", round(cphmfit$n,0),")", sep = "") result3 <- c(3, round(cphmfit$n,0), summary(cphmfit)$conf.int[c(1,3,4)]) cphmfit <- coxph(Surv(futime,fustat)~rx, method="efron", data=ovarian, subset=resid.ds==1) #yes title4 <- paste("Residual Disease (n = ", round(cphmfit$n,0),")", sep = "") result4 <- c(4, round(cphmfit$n,0), summary(cphmfit)$conf.int[c(1,3,4)]) cphmfit <- coxph(Surv(futime,fustat)~rx, method="efron", data=ovarian, subset=resid.ds==2) #no title5 <- paste("No Residual Disease (n = ", round(cphmfit$n,0),")", sep = "") result5 <- c(5, round(cphmfit$n,0), summary(cphmfit)$conf.int[c(1,3,4)]) cphmfit <- coxph(Surv(futime,fustat)~rx, method="efron", data=ovarian, subset=age.f=="< 57 yrs") #yes title6 <- paste("Age < 57 years (n = ", round(cphmfit$n,0),")", sep = "") result6 <- c(6, round(cphmfit$n,0), summary(cphmfit)$conf.int[c(1,3,4)]) cphmfit <- coxph(Surv(futime,fustat)~rx, method="efron", data=ovarian, subset=age.f==">= 57 yrs") #no title7 <- paste("Age >= 57 years (n = ", round(cphmfit$n,0),")", sep = "") result7 <- c(7, round(cphmfit$n,0), summary(cphmfit)$conf.int[c(1,3,4)]) forest <- cbind(data.frame(rbind(title1, title2, title3, title4, title5, title6, title7)), data.frame(rbind(result1, result2, result3, result4, result5, result6, result7))) names(forest) <- c("title", "index", "ss", "hr", "lb", "ub") forest$ub <- ifelse(forest$ub > 5, 5, forest$ub) attach(forest) forest$title <- ordered(forest$title, levels=c(title7, title6, title5, title4, title3, title2, title1)) resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_09X.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(1,1), mar=c(4.2,4.0,2.0,1.0)) plot(seq(-1,6,len=10), seq(0,10, len=10), type="n", xlab="", ylab="", axes=FALSE) Dotplot(title~Cbind(hr,lb,ub), xlim=c(-1,6), pch=16, cex=1.15, xlab="Hazard Ratio (95% Confidence Interval)", ylab="", main="", abline=list(v=1, col="darkgrey", lwd=2, lty=2), data=forest) text(c(1.26,2.85), 9.9, c("Favors Treatment 2", "Favors Treatment 1"), cex=0.75, xpd=TRUE) dev.off() ########################################################################################### # Figure 10 # Change from Baseline ########################################################################################### set.seed(15) mVAS <- c(7.7, 3.1) sdVAS <- c(1.1, 1.3) nVAS <- 40 d_sim <- data.frame(trtgrp=ordered(as.factor(c(rep("A",nVAS), rep("B",nVAS))),levels=c("A","B")), vas=c(rnorm(nVAS, mVAS[1], sdVAS[1]), rnorm(nVAS, mVAS[2], sdVAS[2]))) dim(d_sim) attach(d_sim) mVAS <- by(vas, trtgrp, mean) sdVAS <- by(vas, trtgrp, sd) LBVAS <- mVAS-qnorm(0.025, mean=0, sd=1, lower.tail = FALSE)*(sdVAS/sqrt(nVAS)) UBVAS <- mVAS+qnorm(0.025, mean=0, sd=1, lower.tail = FALSE)*(sdVAS/sqrt(nVAS)) #Figure Change from Baseline resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_10.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(2,2), oma=c(0,0.5,0,0), mar=c(4.5,4,2.5,2)) errbar(c(0.25,0.75), mVAS, UBVAS, LBVAS, xlim=c(0,1), ylim=c(0,10), pch=16, xlab="Treatment", ylab="Change from Baseline", axes=FALSE) abline(h=seq(0,10,1), col="lightgrey") errbar(c(0.25,0.75), mVAS, UBVAS, LBVAS, xlim=c(0,1), ylim=c(0,10), add=TRUE) points(c(0.25,0.75), mVAS, type="p", pch=16, col=c(mycolors[1],mycolors[2]), cex=1) axis(1, at=c(0.25,0.75), labels=levels(trtgrp)) axis(2, at=seq(0,10,1), labels=FALSE) text(-0.13, seq(0,10,1), as.character(seq(-5,5,1)), srt=0, xpd=TRUE) title("Dotplot with 95% Confidence Interval") box() boxplot(d_sim$vas ~ d_sim$trtgrp, xlab="Treatment", ylab="Change from Baseline", ylim=c(0,10), xlim=c(0,3), axes=FALSE, data=d_sim) abline(h=5, col="lightgrey") axis(1, at=c(1,2), labels=FALSE) text(c(1,2), -1.4, levels(trtgrp), xpd=TRUE) axis(2, at=seq(0,10,1), labels=FALSE) text(-0.4, seq(0,10,1), as.character(seq(-5,5,1)), srt=0, xpd=TRUE) title("Boxplot") box() plot(seq(0,10,1),seq(0,1,0.1), xlab="Change from Baseline", type="n", ylab="Cumulative Percent", xlim=c(0,10), ylim=c(0,1), main="", axes=FALSE) plot.ecdf(d_sim[trtgrp=="A",]$vas, xlab="", ylab="", lwd=1, col=mycolors[1], add=TRUE) plot.ecdf(d_sim[trtgrp=="B",]$vas, xlab="", ylab="", lwd=1, col=mycolors[2], add=TRUE) axis(1, at=seq(0,10,1), labels=as.character(seq(-5,5,1))) axis(2, at=seq(0,1,0.1), labels=FALSE) text(-1.4, seq(0,1,0.1), as.character(seq(0,1,0.1)), srt=0, xpd=TRUE) text(c(1.75,6.75),0.5, c("A","B"), adj=0, cex=0.9) title("Empirical Cumulative Distribution Plot") box() plot(density(d_sim$vas), type="n", xlab="Change from Baseline", ylab="Density", xlim=c(0,10), ylim=c(0,0.4), main="", axes=FALSE) lines(density(d_sim[trtgrp=="A",]$vas, kernel="gaussian"), lwd=2, col=mycolors[1]) lines(density(d_sim[trtgrp=="B",]$vas, kernel="gaussian"), lwd=2, col=mycolors[2]) axis(1, at=seq(0,10,1), labels=as.character(seq(-5,5,1))) axis(2, at=seq(0,0.4,0.1), labels=FALSE) text(-1.4, seq(0,0.4,0.1), as.character(seq(0,0.4,0.1)), srt=0, xpd=TRUE) text(c(3.0,7.5),c(0.25,0.39), c("A","B"), adj=0, cex=0.9) title("Density Plot") box() dev.off() #Ecdf(d_sim$vas, group=d_sim$trtgrp, xlab="", ylab="", label.curves=FALSE, lwd=1, col=c("red","blue"), main="", axes=FALSE) #abline(h=c(0,1), col="lightgray", lty=2) #axis(1, at=seq(0,10,1)) #axis(2, at=seq(0,1,0.1), labels=FALSE) #text(-0.4, seq(0,1,0.1), as.character(seq(0,1,0.1)), srt=0, xpd=TRUE) #mtext("Density plot", 3, line=1, cex=1) #box() ########################################################################################### # Figure 11 & Figure 12 # Event Charts ########################################################################################### n <- dim(ovarian)[1] ovarian$startdate <- 17200 ovarian$studyday <- floor(c(0,runif((n-1),min=0,max=(2*365.25)))) ovarian$randdate <- (ovarian$startdate + ovarian$studyday) ovarian$eventdate <- ifelse(ovarian$fustat==1,(ovarian$randdate + ovarian$futime),NA) ovarian$censdate <- ifelse(ovarian$fustat==0,(ovarian$randdate + ovarian$futime),NA) ovarian$survival <- ovarian$futime attach(ovarian) names(ovarian) ovarian.ec <- data.frame(randdate=as.numeric(randdate), eventdate=as.numeric(eventdate), censdate=as.numeric(censdate), survival=as.numeric(survival)) #Figure Event Chart 1 resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_11.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(1,1), mar=c(4.2,4.0,1.0,1.0)) event.chart(ovarian.ec, subset.c=c('randdate','eventdate','censdate'), sort.by = 'randdate', x.lab = 'Treatment Date', y.lab='Patients (Sorted by Treatment Date)', y.axis.custom.at=1:n, y.axis.custom.labels=c(1:n), y.axis='custom', titl='', point.pch=c(16,17,15), point.col=c("black",myredblue[1],myredblue[2]), point.cex=c(1,1,1), legend.plot=FALSE, legend.location='i', legend.cex=1, legend.point.text=c('Treatment','Death','Censored'), legend.point.at = list(c(17200, 17560), c(26, 22)), legend.bty='o') axis(4, at=seq(1,n,1), tck=0.02) legend(17200,27, border="black", legend=c('Treatment','Death','Censored'), pch=c(16,17,15), col=c("black",myredblue[1],myredblue[2])) dev.off() #Figure Event Chart 2 resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_12.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(1,1), mar=c(4.2,4.0,1.0,1.0)) event.chart(ovarian.ec, subset.c=c('randdate','eventdate','censdate'), sort.by = 'survival', x.lab = 'Survival Time (Days)', y.lab='Patients (Sorted by Survival Time)', y.axis.custom.at=1:n, y.axis.custom.labels=c(1:n), y.axis='custom', x.axis.custom.at=seq(0,1230,200), x.axis.custom.labels=seq(0,1230,200), x.axis = 'custom', titl='', point.pch=c(16,17,15), point.col=c("black",myredblue[1],myredblue[2]), point.cex=c(1,1,1), legend.plot=FALSE, legend.location='i', legend.cex=1, legend.point.text=c('Treatment','Death','Censored'), x.reference='randdate', x.julian=TRUE, legend.point.at = list(c(900, 1150), c(0, 4)), legend.bty='o') axis(4, at=seq(1,n,1), tck=0.02) legend(960,4, border="black", legend=c('Treatment','Death','Censored'), pch=c(16,17,15), col=c("black",myredblue[1],myredblue[2])) dev.off() ########################################################################################### # Figure 13 - 17 # Time-to-Event Analysis: Kaplain-Meier & Cox Proportional Hazard Model ########################################################################################### lung$inst.f <- as.factor(lung$inst) lung$age.f <- as.factor(ifelse(lung$age<65, "< 65 yrs", ">= 65 yrs")) lung$sex.f <- factor(lung$sex, levels=c(1,2), label=c("Male","Female")) lung$ph.ecog.group.f <- as.factor(ifelse(lung$ph.ecog==0, "0", ">= 1")) attach(lung) dim(lung) names(lung) table(ph.ecog) table(ph.ecog.group.f) table(sex.f) #Figure 13 kmfit1 <- survfit(Surv(time,status)~1, conf.int=0.95, data=lung) kmfit1 #plot(kmfit1, ylim=c(0,1), conf.int=TRUE, xlab="Survival Time (Days)", ylab="Survival Probability") #median <- kmfit1$time[min(which(kmfit1$surv<=0.5))] median <- summary(kmfit1)$table["median"] LB <- summary(kmfit1)$table["0.95LCL"] UB <- summary(kmfit1)$table["0.95UCL"] resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_13.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(1,1), mar=c(4.2,4.0,1.0,1.0)) plot(seq(0,1080,len=100), seq(0,1,len=100), type="n", xlab="Survival Time (Days)", ylab="Survival Probability", axes=FALSE) axis(1, at=seq(0, 1080, 180)) axis(2, at=seq(0,1,0.1), tck=0.02, labels=FALSE, tick=TRUE) text(-100, seq(0,1,0.1), labels=as.character(seq(0,1,0.1)), srt=0, xpd=TRUE) segments(0, 0.5, median ,0.5, lty=2, lwd=2, col="lightgrey") survplot(kmfit1, time.inc=180, conf="bands", conf.int=.95, col.fill="darkgrey", lwd=2, add=TRUE) text(median+100, 0.5, paste("Median Survival is ", median, " Days (95% CI: ", LB, " - ", UB,")", sep = ""), adj=0, cex=0.90) box() dev.off() #Figure 14 resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_14.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(1,1), mar=c(4.2,4.0,1.0,1.0)) plot(kmfit1, fun='event', ylim=c(0,1), conf.int=TRUE, xlab="Survival Time (Days)", ylab="Cumulative Incidence", axes=FALSE) axis(1, at=seq(0, 1080, 180)) axis(2, at=seq(0,1,0.1), tck=0.02, labels=FALSE, tick=TRUE) text(-50, seq(0,1,0.1), labels=as.character(seq(0,1,0.1)), srt=0, xpd=TRUE) box() dev.off() #Figure 15 kmfit2 <- survfit(Surv(time,status)~sex.f, conf.int=0.95, data=lung) LR <- survdiff(Surv(time,status)~sex.f, rho=0, data=lung)$chisq pvalue <- 1-pchisq(LR,1) resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_15.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(1,1), mar=c(4.2,4.0,1.0,1.0))#, oma=c(1,0.75,1,0.75)) plot(seq(0,1080,len=100), seq(-0.2,1,len=100), type="n", axes=FALSE, xlab="Survival Time (Days)", ylab="Survival Probability") abline(h=seq(0,1,0.1),lty=3,col="lightgrey") axis(1, at=seq(0, 1080, 180)) axis(1, at=seq(180, 1080, 180), lab=FALSE) axis(2, at=seq(0,1,0.1), tck=0.02, labels=FALSE, tick=TRUE) text(-100, seq(0,1,0.1), labels=as.character(seq(0,1,0.1)), srt=0, xpd=TRUE) axis(4, at=seq(0,1,0.1), tck=0.02, labels=FALSE, tick=TRUE) survplot(kmfit2, ylim=c(-0.2, 1), time.inc=180, conf="bars", conf.int=.95, lty=1, lwd=2, col=c(mycolors[1],mycolors[2]), levels.only=TRUE, label.curves=FALSE, n.risk=TRUE, sep.n.risk=0.04, cex.n.risk=0.80, adj.n.risk=0.5, add=TRUE) abline(h=0) text(c(50, 270), 0.75, c("Male", "Female"), col=c(mycolors[1],mycolors[2]), adj=0) text(700, 0.95, paste("p-value (logrank test) =", round(pvalue,dig=3)), cex=0.90, adj=0) text(-35, -0.10, "At risk n:", adj=0, cex=0.85) text(1000, c(-0.16, -0.2), c("Male", "Female"), col=c(mycolors[1],mycolors[2]), adj=0, cex=0.80) box() dev.off() #Figure 16 resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_16.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(1,1), mar=c(4.2,4.0,1.0,1.0)) plot(seq(0,1080,len=100), seq(-0.4,0.2,len=100), type="n", axes=FALSE, xlab="Survival Time (Days)", ylab="Difference in Survival Probability") axis(1, at=seq(0, 1080, 180)) axis(2, at=seq(-0.4,0.2,0.1), lab=FALSE) text(-100, seq(-0.4,0.2,0.1), labels=as.character(seq(-0.4,0.2,0.1)), srt=0, xpd=TRUE) text(200, c(0.015,-0.015), c("Favors Males","Favors Females"), cex=0.90, adj=0, col=myredblue[2]) text(700, 0.2, paste("p-value (logrank test) =", round(pvalue,dig=3)), cex=0.90, adj=0) survdiffplot(kmfit2, time.inc=180, lwd=2, conf="shaded", conf.int=0.95, add=TRUE) box() dev.off() #Figure 16 (extra) #resetGraph(reset.mf=TRUE) #pdf(paste(path, "Gilder_Figure_16_extra.pdf", sep=""), bg="white", onefile=TRUE) #par(mfrow=c(1,1), mar=c(4.2,4.0,1.0,1.0)) #plot(seq(0,1080,len=100), seq(-0.4,0.2,len=100), type="n", axes=FALSE, xlab="Survival Time (Days)", ylab="Difference in Survival Probability") #axis(1, at=seq(0, 1080, 180)) #axis(2, at=seq(-0.4,0.2,0.1), lab=FALSE) #text(-100, seq(-0.4,0.2,0.1), labels=as.character(seq(-0.4,0.2,0.1)), srt=0, xpd=TRUE) #text(200, c(0.015,-0.015), c("Favors Males","Favors Females"), cex=0.90, adj=0, col=myredblue[2]) #text(700, 0.2, paste("p-value (logrank test) =", round(pvalue,dig=3)), cex=0.90, adj=0) #survdiffplot(kmfit2, lwd=2, time.inc=180, conf="bands", conf.int=.95, add=TRUE) #box() #dev.off() #Figure 17 cphmfit <- coxph(Surv(time, status)~age + sex, method="efron", data=lung) cphmdiagnostic <- cox.zph(cphmfit, transform="km") cphmfit summary(cphmfit) cphmdiagnostic resetGraph(reset.mf=TRUE) pdf(paste(path, "Gilder_Figure_17.pdf", sep=""), bg="white", onefile=TRUE) par(mfrow=c(2,2), mar=c(4.2,4,2.5,1), oma=c(0,0.5,0,0)) survplot(kmfit2, time.inc=180, logt=TRUE, loglog=TRUE, xlim=c(2,7), title="", xlab="log(Survival Time (Days))", ylab="log(-log(Survival Probability))", ylim=c(-5,2), lty=1, lwd=2, col=c(mycolors[1],mycolors[2]), conf="none", label.curves=TRUE, levels.only=TRUE) title("A. Proportional Hazards") box() plot(cphmdiagnostic[2], ann=F, lwd=5, cex=0.75) abline(lm(cphmdiagnostic[2]$y ~ cphmdiagnostic[2]$x)$coefficients, lwd=2, lty=5, col=myredblue[2]) title(xlab="Survival Time (Days)", ylab="Beta(t) for Sex", main="B. Proportional Hazards") dfbeta <- residuals(cphmfit, type="dfbeta") plot(dfbeta[,2], ylab=capwords(names(coef(cphmfit))[2]), cex=0.75, main="C. Influential Observations") abline(h=0, lty=2) res <- residuals(cphmfit, type="martingale") X <- as.matrix(lung[,c("age", "sex")]) for (j in 1:1) { plot(X[,j], res, xlab=c("Age (Years)", "Sex")[j], ylab="Residuals", cex=0.75, main="D. Martingale-Residual Plot") abline(h=0, lty=2) lines(lowess(X[,j], res, iter=0)) } dev.off()