### R code from vignette source 'LocalControl-jss-2020.Rnw' ################################################### ### code chunk number 1: preliminaries ################################################### options(prompt = "R> ", continue = "+ ", width = 70, useFancyQuotes = FALSE) ################################################### ### code chunk number 2: init_libs_vars_calcs ################################################### # toggle long computations for fast pdf building doCalcs = F # load packages if(!require("data.table")){ install.packages("data.table") library("data.table") } if(!require("colorspace")){ install.packages("colorspace") library("colorspace") } if(!require("RColorBrewer")){ install.packages("RColorBrewer") library("RColorBrewer") } if(!require("gridExtra")){ install.packages("gridExtra") library("gridExtra") } if(!require("ggplot2")){ install.packages("ggplot2") library("ggplot2") } if(!require("rpart")){ install.packages("rpart") library("rpart") } if(!require("rpart.plot")){ install.packages("rpart.plot") library("rpart.plot") } if(!require("LocalControl")){ install.packages("LocalControl") library("LocalControl") } if(!require("xtable")){ install.packages("xtable") library("xtable") } # frequently used variables data(lindner) all7Vars <- c("stent", "height", "female", "diabetic", "acutemi", "ejecfrac", "ves1proc") numClusters <- c(1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50) t0Red = "#E41A1C" t0Fill = rgb(red = 228, green = 26, blue = 28, alpha = 127, maxColorValue = 255) t1Blue = "#377EB8" t1Fill = rgb(red = 55, green = 126, blue = 184, alpha = 127, maxColorValue = 255) t2Green = "#4DAF4A" if (doCalcs) { linResults = LocalControlClassic( data = lindner, clusterVars = all7Vars, treatmentColName = "abcix", outcomeColName = "cardbill", clusterCounts = numClusters) linRes = LocalControl(data = lindner, clusterVars = all7Vars, treatmentColName = "abcix", outcomeColName = "cardbill", treatmentCode = 1) linSults = LocalControl(data = lindner, treatmentColName = "abcix", outcomeColName = "cardbill", treatmentCode = 1, clusterVars = all7Vars, radDecayRate = .95, radMinFract = .01 ) sdata = lindner linSummary = summary(linSults) linSummary$male_stent = 0 linSummary$female_stent = 0 linSummary$male_without = 0 linSummary$female_without = 0 linSummary$maleAverage = 0 linSummary$femaleAverage = 0 linSummary$medianLTD = 0 stentMales = which(sdata$female == 0 & sdata$stent == 1) stentFemales = which(sdata$female == 1 & sdata$stent == 1) nosMales = which(sdata$female == 0 & sdata$stent == 0) nosFemales = which(sdata$female == 1 & sdata$stent == 0) allMen = which(sdata$female == 0) allWomen = which(sdata$female == 1) for (i in 1:nrow(linSummary)) { linSummary$male_stent[i] = mean(linSults$ltds[stentMales,i], na.rm = T) linSummary$female_stent[i] = mean(linSults$ltds[stentFemales,i], na.rm = T) linSummary$male_without[i] = mean(linSults$ltds[nosMales,i], na.rm = T) linSummary$female_without[i] = mean(linSults$ltds[nosFemales,i], na.rm = T) linSummary$maleAverage[i] = mean(linSults$ltds[allMen,i], na.rm = T) linSummary$femaleAverage[i] = mean(linSults$ltds[allWomen,i], na.rm = T) linSummary$medianLTD[i] = median(linSults$ltds[,i], na.rm = T) } linRad33 = linSults$ltds$rad_33 saveRDS(linRes, 'calcs/linRes.rds') saveRDS(linSummary, 'calcs/linSummary.rds') saveRDS(linRad33, 'calcs/linRad33.rds') lindnerResults = list() ltdds = list() for (cc in numClusters) { objName = paste0("UPSnnltd", cc) nnObj = linResults[[objName]] ltdds[[objName]] = nnObj$pbindif[nnObj$nnhbindf[,1]] } lindnerResults[["summary"]] = linResults[["summary"]] lindnerResults[["ltdds"]] = ltdds lindnerResults[["yname"]] = linResults[["UPSaccum.pars"]][,"yvar"] saveRDS(lindnerResults, 'calcs/lindnerResults.rds') } else{ lindnerResults = readRDS('calcs/lindnerResults.rds') linSummary = readRDS('calcs/linSummary.rds') linRad33 = readRDS('calcs/linRad33.rds') linRes = readRDS('calcs/linRes.rds') } set.seed(253748) N = 10000 weight = c(rnorm(N/2, 75, 15), rnorm(N/2, 75, 15)) dosage = weight + c(rnorm(N/2, 0, 15), rnorm(N/2, 0, 5)) trmt = c(rep(1, N/2), rep(0, N/2)) ADR = abs(weight - dosage) noise1 = rnorm(n = N, 0, 1) noise2 = rnorm(n = N, 0, 1) simData = data.frame(weight, trmt, dosage, ADR, noise1, noise2) if (doCalcs) { xSults = LocalControl(data = simData, treatmentColName = "trmt", treatmentCode = 1, outcomeColName = "ADR", clusterVars = c("weight", "dosage"), radMinFract = .01, radDecayRate = 0.95, numThreads = 4) xSultsT1 = xSults$outcomes$T1[,"rad_91"] xSultsT0 = xSults$outcomes$T0[,"rad_91"] saveRDS(xSultsT1, 'calcs/xSultsT1.rds') saveRDS(xSultsT0, 'calcs/xSultsT0.rds') } else{ xSultsT1 = readRDS( 'calcs/xSultsT1.rds') xSultsT0 = readRDS( 'calcs/xSultsT0.rds') } if (doCalcs) { noisyVars = c("weight", "dosage", "noise1", "noise2") noisySults = LocalControl( simData, treatmentCode = 1, treatmentColName = "trmt", outcomeColName = "ADR", clusterVars = noisyVars, radMinFract = .01, radDecayRate = 0.95 , numThreads = 4) fixedRads = summary(noisySults)$radius varCombinations = expand.grid(0:1, 0:1, 0:1, 0:1) ltext = apply(X = varCombinations, MARGIN = 1, FUN = function(x) paste0(x, collapse = '')) ltdVecs = list() ltdVecs[[1]] = rep(summary(noisySults)$ltd[1], nrow(summary(noisySults))) for (i in 2:16) { varSS = noisyVars[which(varCombinations[i,] == 1)] scaleFactor = sqrt(length(varSS)) / sqrt(length(noisyVars)) scaleRads = fixedRads * scaleFactor sults = LocalControl( simData, treatmentCode = 1, treatmentColName = "trmt", outcomeColName = "ADR", clusterVars = varSS, radiusLevels = scaleRads, numThreads = 4) ltdVecs[[i]] = summary(sults)$ltd } ltdFrame = data.frame(ltdVecs) names(ltdFrame) = ltext saveRDS(object = varCombinations, file = "calcs/ffvc.rds") resultSummary = summary(sults) saveRDS(object = resultSummary, file = "calcs/resultSummary.rds") saveRDS(object = ltdVecs, file = "calcs/ffltdv.rds") saveRDS(object = ltdFrame, file = "calcs/ffltdf.rds") } else{ varCombinations = readRDS(file = "calcs/ffvc.rds") resultSummary = readRDS(file = "calcs/resultSummary.rds") ltdVecs = readRDS(file = "calcs/ffltdv.rds") ltdFrame = readRDS(file = "calcs/ffltdf.rds") } if (doCalcs) { linCI = LocalControlNearestNeighborsConfidence(data = lindner, clusterVars = all7Vars, treatmentColName = "abcix", outcomeColName = "cardbill", treatmentCode = 1, nBootstrap = 250) saveRDS(object = linCI, file = "calcs/linCI.rds") } else{ linCI = readRDS(file = "calcs/linCI.rds") } # generate survival simulation weibullSim = function(N, lambda, rho, betaage, betabmi) { bmi = rnorm(N, mean = 26, sd = 4) age = runif(N) * 47 + 18 pbmi = (bmi - min(bmi)) / (max(bmi) - min(bmi)) * 0.8 + 0.1 page = (age - min(age)) / (max(age) - min(age)) * 0.8 + 0.1 drug = 1 - rbinom(N, 1, (pbmi + page) / 2) et = exp(bmi * betabmi + age * betaage) Tlat = (-log(runif(n = N)) / (lambda * et))^(1 / rho) C = runif(N) * 30 time = pmin(Tlat, C) status = as.numeric(Tlat <= C) data.frame(id = 1:N, drug, age, bmi, time, status) } set.seed(24563568) survSimData = weibullSim(10000, 1e-10, 2.6, log(1.2), log(1.45)) if (doCalcs) { sSults = LocalControl(survSimData, treatmentColName = "drug", timeColName = "time", outcomeType = "survival", outcomeColName = "status",treatmentCode = 1, clusterVars = c("age", "bmi"), radMinFract = .01) sSimFTs = sSults$Failtimes sSimT1R1 = sSults$KM$T1$rad_1 sSimT0R1 = sSults$KM$T0$rad_1 sSimT1R21 = sSults$KM$T1$rad_21 sSimT0R21 = sSults$KM$T0$rad_21 saveRDS(object = sSimFTs, file = "calcs/sSimFTs.rds") saveRDS(object = sSimT1R1, file = "calcs/sSimT1R1.rds") saveRDS(object = sSimT0R1, file = "calcs/sSimT0R1.rds") saveRDS(object = sSimT1R21, file = "calcs/sSimT1R21.rds") saveRDS(object = sSimT0R21, file = "calcs/sSimT0R21.rds") } else{ sSimFTs = readRDS(file = "calcs/sSimFTs.rds") sSimT1R1 = readRDS(file = "calcs/sSimT1R1.rds") sSimT0R1 = readRDS(file = "calcs/sSimT0R1.rds") sSimT1R21 = readRDS(file = "calcs/sSimT1R21.rds") sSimT0R21 = readRDS(file = "calcs/sSimT0R21.rds") } data(framingham) shadeVect = c("Nonsmoker" = "#80000055", "Smoker" = "#0000FF55") colVect = c("Smoker" = "red", "Nonsmoker" = "blue") ltyVect = c("Hypertension" = "dashed", "Death" = "solid") if (doCalcs) { FHSResults = LocalControl( data = framingham, outcomeType = "survival", treatmentColName = "cursmoke", treatmentCode = 1, timeColName = "time_outcome", outcomeColName = "outcome", clusterVars = c("female", "totchol", "age", "bmi", "BPVar", "heartrte", "glucose")) FHSConfidence = LocalControlCompetingRisksConfidence(LCCompRisk = FHSResults, confLevel = "90%") getLines2Plot <- function(sult, conf, rad, fc, tr) { outFrame = data.frame(TIMES = sult$Failtimes, UPPER = conf[[tr]][[fc]]$UPPER_CONF[,rad], LOWER = conf[[tr]][[fc]]$LOWER_CONF[,rad], CIF = sult$CIF[[fc]][[tr]][,rad]) outFrame } FHSPlotLines = list() for (radLim in c("rad_1","rad_11")) { plotLines = list() t1lines2plot = getLines2Plot(FHSResults, FHSConfidence, radLim, "Failcode_1", "T1") t1lines2plot$CLASS = "Smoker" t1lines2plot$FAIL = "Hypertension" t1lines2plot$NG = paste(t1lines2plot$CLASS, t1lines2plot$FAIL) plotLines[["t1HT"]] = t1lines2plot t0lines2plot = getLines2Plot(FHSResults, FHSConfidence, radLim, "Failcode_1", "T0") t0lines2plot$CLASS = "Nonsmoker" t0lines2plot$FAIL = "Hypertension" t0lines2plot$NG = paste(t0lines2plot$CLASS, t0lines2plot$FAIL) plotLines[["t0HT"]] = t0lines2plot t1DeathLines = getLines2Plot(FHSResults, FHSConfidence, radLim, "Failcode_2", "T1") t1DeathLines$CLASS = "Smoker" t1DeathLines$FAIL = "Death" t1DeathLines$NG = paste(t1DeathLines$CLASS, t1DeathLines$FAIL) plotLines[["t1Death"]] = t1DeathLines t0DeathLines = getLines2Plot(FHSResults, FHSConfidence, radLim, "Failcode_2", "T0") t0DeathLines$CLASS = "Nonsmoker" t0DeathLines$FAIL = "Death" t0DeathLines$NG = paste(t0DeathLines$CLASS, t0DeathLines$FAIL) plotLines[["t0Death"]] = t0DeathLines FHSPlotLines[[radLim]] = plotLines } FHSSummary = summary(FHSResults) saveRDS(object = FHSSummary, file = "calcs/FHSSummary.rds") saveRDS(object = FHSPlotLines, file = "calcs/FHSPlotLines.rds") } else{ FHSPlotLines = readRDS(file = "calcs/FHSPlotLines.rds") FHSSummary = readRDS(file = "calcs/FHSSummary.rds") } avgDif <- function(uncorrected, corrected) { return(sum(uncorrected - corrected, na.rm = TRUE)/ length(which(!is.na(corrected)))) } setColAlpha = function(col, intensity) { iint = 1 - as.double(intensity) cc = col2rgb(col)/255 nc = rgb(red = cc[1], green = cc[2], blue = cc[3], alpha = iint) nc } if (doCalcs) { lind = data.table(lindner) lindavgs = lind[, .(mean(cardbill)), by = .(abcix)] lindLC = LocalControlClassic( data = lindner, clusterVars = c("stent", "female", "acutemi"), treatmentColName = "abcix", outcomeColName = "cardbill", clusterCounts = c(30)) lind$c = lindLC$UPSnnltd30$nnhbindf localOutcomes = data.table(lindLC$UPSnnltd30$binmean) localOutcomes$LTDobs = localOutcomes$`abcix = 1` - localOutcomes$`abcix = 0` setkey(lind, c) setkey(localOutcomes, BIN) lind = lind[localOutcomes, nomatch = NA] set.seed(1772) randomData = lind # copy for randomizing randomLTDs = data.table(cluster = 1:30, artltd = rep(0.0, 30)) #aggregate df setkey(randomLTDs, cluster) for (i in 1:100) { randomData[, "c"] = sample(x = lind$c, replace = F) # random 2000 clusters w/ same sizes randomAverages = randomData[, .(mean(cardbill)), by = .(abcix, c)] # avg outcome per trt in cluster randomAverages[abcix == 0, "V1"] = -randomAverages[ abcix == 0, "V1"] # flip sign on t0s before sum rltds = randomAverages[, .(sum(V1)), by = .(c)] # (t1 - t0) per cluster setkey(rltds, c) randomLTDs = randomLTDs[ rltds, nomatch = NA] } ccc = lind setkey(ccc, c) # keys are used to join data.tables setkey(randomLTDs, cluster) # next line joins randomized cluster values back with main table based on cluster value. ccc = ccc[randomLTDs, nomatch = NA] # set artltd value for each patient based on original clustering ccc[ which(is.na(ccc$LTDobs)), "artltd"] = NA # remove same number of noninformative clusters/patients # above uses means (11 values for 100 resamplings) # this uses all 1100, rather than the average randomLTDs[which(is.na(localOutcomes$LTDobs)), 3:ncol(randomLTDs)] = NA allrands = as.vector(as.matrix(randomLTDs[, 3:ncol(randomLTDs)])) saveRDS(object = allrands, file = "calcs/allRands.rds") saveRDS(object = ccc, file = "calcs/ccc.rds") saveRDS(object = localOutcomes, file = "calcs/localoutcomes.rds") } else{ ccc = readRDS(file = "calcs/ccc.rds") allrands = readRDS(file = "calcs/allRands.rds") localOutcomes = readRDS(file = "calcs/localoutcomes.rds") } if (doCalcs) { dflindner = lindner dflindner$c = lindLC$UPSnnltd30$nnhbindf$HclusBin dflindner$t = dflindner$abcix dflindner$y = dflindner$cardbill # # [2] Compute observed Local Treatment Differences (within clusters) = # Local Differences in Mean y-outcomes between Treatment Groups. # dfLTD <- do.call( rbind, lapply( split(dflindner, dflindner$c), function(x) { n1 = sum(x$t) n0 = length(x$t) - n1 if (n1 == 0 || n0 == 0) LTD = NA else LTD = sum(x$y * x$t)/n1 - sum(x$y * (1 - x$t))/n0 data.frame(c = x$c[1], LTD = LTD, w = length(x$c)) } ) ) # Create data.frame to hold full replicates of Local Treatment Effect-Size Estimates (many within-cluster ties) dfLreps = as.data.frame(cbind(dflindner$c, dflindner$y)) colnames(dfLreps) = c("c","LTD") # Here, Local Treatment Effect-Size Estimates are LTDs... for (i in 1:length(dfLreps[,1])) { dfLreps$LTD[i] = dfLTD$LTD[dflindner$c[i]] # insert observed LTD values } dfLTDobs = dfLreps # Save Observed LTD Distribution from lindner data... #hist(dfLTDobs[,2]) ############################################# # Calculate random permutation (NULL) distribution of LTDs for N full Replications of given Clusters / Subgroups... N = 100 # Number of Random Permutation Replications in LC Confirm Phase Calculations... nobs = nrow(lindner) set.seed(12345) # Set seed for Monte Carlo pseudo random sequence... for (i in 1:N) { crand <- as.vector(rnorm(nobs)) # next vector of Random numbers... srand <- dflindner$c[order(crand)] # corresponding vector of Permuted cluster numbers pdf <- as.data.frame(cbind(srand, dflindner$y, dflindner$t)) # pdf = PERMUTED, 3 variable data.frame colnames(pdf) <- c("c","y","t") pdfc <- do.call( rbind, lapply( split(pdf, pdf$c), function(x) { n1 = sum(x$t) n0 = length(x$t) - n1 if (n1 == 0 || n0 == 0) LTD = NA else LTD = sum(x$y * x$t)/n1 - sum(x$y * (1 - x$t))/n0 data.frame(c = x$c[1], LTD = LTD, w = length(x$c)) } ) ) colnames(pdfc) <- c("c","LTD","w") for (j in 1:length(dfLreps[,1])) { dfLreps$LTD[j] = pdfc$LTD[dflindner$c[j]] # insert permuted LTD values } if ( i == 1 ) dfconfirm = dfLreps else dfconfirm = rbind(dfconfirm, dfLreps) } #str(dfconfirm) # full random permutation distribution of NULL LTD estimates... permLTDs = dfconfirm$LTD ksobserved = ks.test(dfLTDobs$LTD, dfconfirm$LTD) # Warning message: # In ks.test(dfLTDobs$LTD, dfconfirm$LTD) : # p-value will be approximate in the presence of ties #length(dfLTDobs$LTD) # 996 #length(dfconfirm$LTD) # 99600 #ksobserved # Observed KS D = 0.42208 #ksobserved$p.value # = 0 ##################################################################### NKS = 1000 set.seed(54321) ksDperm <- as.vector(rnorm(NKS)) for (i in 1:NKS) { crand <- as.vector(rnorm(nobs)) srand <- dflindner$c[order(crand)] pdf <- as.data.frame(cbind(srand, dflindner$y, dflindner$t)) colnames(pdf) <- c("c","y","t") pdfc <- do.call( rbind, lapply( split(pdf, pdf$c), function(x) { n1 = sum(x$t) n0 = length(x$t) - n1 if(n1 == 0 || n0 == 0) LTD = NA else LTD = sum(x$y * x$t)/n1 - sum(x$y * (1 - x$t))/n0 data.frame(c = x$c[1], LTD = LTD, w = length(x$c)) } ) ) colnames(pdfc) <- c("c","LTD","w") for (j in 1:length(dfLreps[,1])) { dfLreps$LTD[j] = pdfc$LTD[dflindner$c[j]] # insert permuted LTD values } ksobserved <- ks.test(dfLreps$LTD, dfconfirm$LTD) ksDperm[i] <- ksobserved$statistic } #ksDperm[order(ksDperm)] # Max of 100 NULL KS D stats = 0.4122992 } ################################################### ### code chunk number 3: echo_install_lc (eval = FALSE) ################################################### ## install.packages("LocalControl") ## library("LocalControl") ## data("lindner") ################################################### ### code chunk number 4: echo_lcc_covars (eval = FALSE) ################################################### ## all7Vars = c("stent", "height", "female", "diabetic", ## "acutemi", "ejecfrac", "ves1proc") ## numClusters = c(1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50) ################################################### ### code chunk number 5: echo_lcc_execution (eval = FALSE) ################################################### ## linResults = LocalControlClassic( data = lindner, ## clusterVars = all7Vars, treatmentColName = "abcix", ## outcomeColName = "cardbill", clusterCounts = numClusters) ################################################### ### code chunk number 6: echo_lindner_plot (eval = FALSE) ################################################### ## UPSLTDdist(linResults, ylim = c(-2500, 5000)) ################################################### ### code chunk number 7: figure_2 ################################################### ltdds = lindnerResults[["ltdds"]] infPatFrac = lindnerResults[["summary"]]$Informative_fraction actClusterCounts = lindnerResults[["summary"]]$Clusters_created ltdm = sapply(ltdds, function(x) mean(x,na.rm=T)) boxes = lapply(ltdds, FUN = function(x) boxplot(x, plot=FALSE)) bStats = sapply(boxes, FUN = function(x) x$stats) yname = lindnerResults[["yname"]] globMean = ltdm[1] par(mar = c(5,5,2,5)) plot(x = actClusterCounts, y = ltdm, ylim = c(-3000,5000), lwd = 3, type = "l", col = "black", xlab = "Number of Clusters", ylab = paste0(yname," LTD Distribution"), main = paste0(yname," Local Treatment Differences vs. Number of Clusters")) grid() lines(actClusterCounts, bStats[2,], col = "black", lty = 3, lwd = 2)#lower hinge lines(actClusterCounts, bStats[4,], col = "black", lty = 3, lwd = 2)#upper hinge abline(h = globMean, col = "purple", lty = 2, lwd = 2) par(new = T) plot( x = actClusterCounts, y = infPatFrac, type = "l", lwd = 3, lty = 3, axes=F, xlab=NA, ylab=NA, cex=1.2, col = "green", ylim = c(0,1)) axis(side = 4) mtext(side = 4, line = 3, 'Fraction of Patients Informative') legend( "bottomright", col = c("purple", "black", "black", "green"), lty = c(2, 1, 3, 3), lwd = c(2, 3, 2, 2), bg = "white", legend = c( "Global Treatment Difference", "Mean of the LTD distribution", "Lower and upper quartile", "Fraction of Patients Informative")) ################################################### ### code chunk number 8: figure_2_6 ################################################### par(mfrow = c(1,1), mar = c(4.1, 4.1, .2, .2)) plxmin = min(allrands, na.rm =T) plxmax = max(allrands, na.rm =T) plot(ecdf(allrands), verticals=TRUE, do.points=FALSE, ann=FALSE, col="green3", lwd=2, pch=46, xlim=c(plxmin,plxmax), xaxt="n") grid() lines(ecdf(localOutcomes$LTDobs), verticals=TRUE, do.points=FALSE, ann=FALSE, col="red", pch=46) axis(1, at=axTicks(1), labels=sprintf(ifelse(axTicks(1) >= 0, "$%s","-$%s"), abs(axTicks(1))), las=1) abline(v=0, lty="dashed") legend("topleft", legend = c("Observed", "Artifical"), col = c("red", "green3"), lty = 1, bg = "white") title(main=" ", ylab="Cumulative Probability", xlab="Local Treatment Differences") ################################################### ### code chunk number 9: figure_3 ################################################### xSim = simData[order(simData$ADR),] xSim$colIntensity = seq(1 , 240 , length.out = nrow(xSim))/255 xSim$pc = ifelse(xSim$trmt == 1, t1Blue, t0Red) xSim$pca = apply(X = xSim, MARGIN = 1, FUN = function(x) {setColAlpha(col = x["pc"], intensity = x["colIntensity"])}) par(mar = c(2.1,2.5,1.1,1), xaxs = "i",yaxs = "r") plot( NA, xlab="", ylab="", xaxt = "n", yaxt = "n", xlim = c(35,115), ylim = c(35,115), pch = 16, lwd=7, cex.lab = 1.75, cex.axis = 1.4) title(ylab = "Dosage (40-110 mcg/kg)", line=0.3, cex.lab = 1.5) title(xlab = "Weight (40-110 kg)", line=0.3, cex.lab = 1.5) grid() box() points(xSim$weight, xSim$dosage, col = xSim$pca, pch = 16) legend("bottomright", legend = c("Treatment 0", "Treatment 1"), col = c(t0Red, t1Blue), pch = 16) ################################################### ### code chunk number 10: table_3 (eval = FALSE) ################################################### ## t1t0 = xSim ## t1s = xSim[which(xSim$trmt == 1),] ## t0s = xSim[which(xSim$trmt == 0),] ## ## npats = c(nrow(t1t0), nrow(t0s), nrow(t1s), NA) ## wmean = c(mean(t1t0$weight), mean(t0s$weight), mean(t1s$weight), t.test(t1s$weight, t0s$weight)$p.value) ## wvars = c(sqrt(var(t1t0$weight)), sqrt(var(t0s$weight)), sqrt(var(t1s$weight)), var.test(t1s$weight, t0s$weight)$p.value) ## dmean = c(mean(t1t0$dosage), mean(t0s$dosage), mean(t1s$dosage), t.test(t1s$dosage, t0s$dosage)$p.value) ## dvars = c(sqrt(var(t1t0$dosage)), sqrt(var(t0s$dosage)), sqrt(var(t1s$dosage)), var.test(t1s$dosage, t0s$dosage)$p.value) ## amean = c(mean(t1t0$ADR), mean(t0s$ADR), mean(t1s$ADR), t.test(t1s$ADR, t0s$ADR)$p.value) ## avars = c(sqrt(var(t1t0$ADR)), sqrt(var(t0s$ADR)), sqrt(var(t1s$ADR)), var.test(t1s$ADR, t0s$ADR)$p.value) ## ## tabOut = t(data.frame(n = npats, weightmean = wmean, weightvar = wvars, dosagemean = dmean, dosagevar = dvars, adrmean = amean, adrvar = avars)) ################################################### ### code chunk number 11: echo_sim_generation (eval = FALSE) ################################################### ## set.seed(253748) ## N = 10000 ## weight = c(rnorm(N/2, 75, 15), rnorm(N/2, 75, 15)) ## dosage = weight + c(rnorm(N/2, 0, 15), rnorm(N/2, 0, 5)) ## trmt = c(rep(1, N/2), rep(0, N/2)) ## ADR = abs(weight - dosage) ## noise1 = rnorm(n = N, 0, 1) ## noise2 = rnorm(n = N, 0, 1) ## xSim = data.frame(weight, trmt, dosage, ADR, noise1, noise2) ################################################### ### code chunk number 12: figure_4 ################################################### t1UVals = xSim[which(xSim$trmt == 1),"ADR"] t0UVals = xSim[which(xSim$trmt == 0),"ADR"] ulegtxt = c("Treatment 0", "Treatment 1" , paste0("T0 mean outcome = ", round(mean(t0UVals, na.rm = T),digits = 2)), paste0("T1 mean outcome = ", round(mean(t1UVals, na.rm = T),digits = 3))) t1CVals = xSultsT1 t0CVals = xSultsT0 clegtxt = c("Treatment 0", "Treatment 1" , paste0("T0 mean outcome = ", round(mean(t0CVals, na.rm = T),digits = 2)), paste0("T1 mean outcome = ", round(mean(t1CVals, na.rm = T),digits = 3))) par(mar = c(3.1,2.1,1.1,1.1), mfrow = c(1,2)) hist(t1UVals,breaks = seq(0, 60, by = 2), ylim = c(0,2500), xlim = c(0,30), yaxt = "n", panel.first = {grid(); box()}, xaxs="i", yaxs="i", xlab="", col = t1Fill, lwd = 1, main = "", ylab = "",mgp=c(1.2,0.5,0)) hist(t0UVals, col = t0Fill, lwd = 5, add = T,breaks = seq(0,60, by = 2)) legend("topright", legend = ulegtxt, ncol = 1,cex = 1.1,x.intersp = c(-1, -1, 1, 1), fill = c(t0Fill, t1Fill, NA, NA), bty = "n", border = c("black", "black", NA, NA), lty = c(NA,NA,3,3), lwd = c(NA,NA,5,5), col = c(NA,NA,t0Red,t1Blue)) abline(v = mean(t1UVals, na.rm = TRUE), col = t1Blue, lty = 3, lwd = 5) abline(v = mean(t0UVals, na.rm = TRUE), col = t0Red, lty = 3, lwd = 5) title(ylab = "Frequency", line=0.3, cex.lab = 1.5) title(xlab = "ADR", line = 2, cex.lab = 1.5) hist(t1CVals,breaks = seq(0, 60, by = 2), ylim = c(0,2500), xlim = c(0,30), yaxt = "n", panel.first = {grid(); box()}, xaxs="i", yaxs="i", xlab="", col = t1Fill, lwd = 1, main = "", ylab = "",mgp=c(1.2,0.5,0)) hist(t0CVals, col = t0Fill, lwd = 5, add = T,breaks = seq(0,60, by = 2)) legend("topright", legend = clegtxt, ncol = 1,cex = 1.1,x.intersp = c(-1, -1, 1, 1), fill = c(t0Fill, t1Fill, NA, NA), bty = "n", border = c("black", "black", NA, NA), lty = c(NA,NA,3,3), lwd = c(NA,NA,5,5), col = c(NA,NA,t0Red,t1Blue)) abline(v = mean(t1CVals, na.rm = TRUE), col = t1Blue, lty = 3, lwd = 5) abline(v = mean(t0CVals, na.rm = TRUE), col = t0Red, lty = 3, lwd = 5) title(ylab = "Frequency", line=0.3, cex.lab = 1.5) title(xlab = "ADR", line = 2, cex.lab = 1.5) ################################################### ### code chunk number 13: echo_nnsim_execution (eval = FALSE) ################################################### ## xSults = LocalControl(data = xSim, ## treatmentColName = "trmt", ## treatmentCode = 1, ## outcomeColName = "ADR", ## clusterVars = c("weight", "dosage"), ## radMinFract = .01, ## radDecayRate = 0.95, ## numThreads = 4) ################################################### ### code chunk number 14: table_4 ################################################### #ltdFrame = readRDS('calcs/ffltds.rds') noisyVars = c("weight", "dosage", "noise1", "noise2") difs = numeric() difs[1] = 0 for(i in 2:ncol(ltdFrame)){ difs[i] = avgDif(ltdFrame[1:92,1], ltdFrame[1:92,i]) } outmat = expand.grid(c(-1,1),c(-1,1),c(-1,1),c(-1,1)) outmat = data.frame(outmat) names(outmat) = noisyVars outmat$difs = difs dm = matrix(data = 0, nrow = 16, ncol = 6) dm[,6] = 2 print(xtable(outmat, digits = dm), include.rownames = F, sanitize.text.function = function(x){x}, floating=FALSE) ################################################### ### code chunk number 15: echo_ff_step1 (eval = FALSE) ################################################### ## noisyVars = c("weight", "dosage", "noise1", "noise2") ## noisySults = LocalControl( xSim, ## treatmentColName = "trmt", ## outcomeColName = "ADR", ## clusterVars = noisyVars, ## radMinFract = .01, ## radDecayRate = 0.95 ) ## fixedRads = summary(noisySults)$limits ################################################### ### code chunk number 16: echo_ff_step2 (eval = FALSE) ################################################### ## varCombinations = expand.grid(0:1, 0:1, 0:1, 0:1) ## ltext = apply(X = varCombinations, MARGIN = 1, ## FUN = function(x) paste0(x, collapse = '')) ## ltdVecs = list() ## ltdVecs[[1]] = rep(summary(noisySults)$ltd[1], nrow(summary(noisySults))) ## for(i in 2:16){ ## varSS = noisyVars[which(varCombinations[i, ] == 1)] ## scaleFactor = sqrt(length(varSS)) / sqrt(length(noisyVars)) ## scaleRads = fixedRads * scaleFactor ## sults = LocalControl( xSim, ## treatmentColName = "trmt", ## outcomeColName = "ADR", ## clusterVars = varSS, ## radiusLevels = scaleRads ) ## ltdVecs[[i]] = summary(sults)$ltd ## } ## ltdFrame = data.frame(ltdVecs) ## names(ltdFrame) = ltext ################################################### ### code chunk number 17: echo_ff_step3 (eval = FALSE) ################################################### ## avgDif = function(uncorrected, corrected){ ## return(sum(uncorrected - corrected, na.rm = TRUE)/ ## length(which(!is.na(corrected)))) ## } ## difs = numeric() ## difs[1] = 0 ## for(i in 2:ncol(ltdFrame)){ ## difs[i] = avgDif(ltdFrame[1:92, 1], ltdFrame[1:92, i]) ## } ## outmat = data.frame(expand.grid(c(-1,1), c(-1,1), c(-1,1), c(-1,1))) ## names(outmat) = noisyVars ## outmat$difs = difs ################################################### ### code chunk number 18: figure_5 ################################################### ltext = apply(varCombinations, 1, function(x) paste0(x, collapse = "")) ltys = apply(varCombinations[,3:4], 1, function(x) ifelse(sum(x) == 1, 5, ifelse(sum(x) == 2, 2, 1))) cols = apply(varCombinations, 1, function(x) ifelse(x[1] == 1, ifelse(x[2] == 1, "purple", "red"), ifelse(x[2] == 1, "blue", "green"))) cols[1] = "black" prads = resultSummary$pct_radius[1:(nrow(resultSummary) - 1)] xl = c(1,min(prads)) par(mar=c(4.1,4.1,2.1,0.5)) plot(NA, xlim = xl, ylim = c(0,9), xlab = "Fraction of maximum radius", ylab = "ADR treatment difference", log = "x", panel.first = grid(equilogs=FALSE),main = "Full-factorial simulation analysis") for(i in 16:1){ ltdz = ltdVecs[[i]] lines(x = prads, y = ltdz[1:(length(ltdz)-1)], xlim = xl, ylim = c(0, 9), col = cols[i], lty = ltys[i], pch = 16, lwd = 3) } lt = c("Uncorrected", "Weight and noise only", "Dosage and noise only", "Both weight and dosage", "Noise only") box() legend("bottomleft",col = cols, lwd = 2, lty = ltys, legend = ltext) legend(x = .702, y = 1.66,cex = 1, legend = lt, col = c("black", "red", "blue", "purple", "green"), lwd = NA, pch = 15) ################################################### ### code chunk number 19: echo_nnlc_regression ################################################### model <- formula("difs ~ (weight + dosage + noise1 + noise2)^4") fit <- glm( difs ~ 1, data = outmat, family = gaussian) fit.AIC <- step( fit, model, k = 2, trace = 0, direction = "both") regTable = summary.glm( fit.AIC )$coef ################################################### ### code chunk number 20: table_5 ################################################### dm = matrix(data = 2, nrow = 4, ncol = 5) dm[,5] = -2 print(xtable(regTable, digits = dm), include.rownames = T, floating=FALSE) ################################################### ### code chunk number 21: echo_nnlc_regression (eval = FALSE) ################################################### ## model = formula("difs ~ (weight + dosage + noise1 + noise2)^4") ## fit = glm( difs ~ 1, data = outmat, family = gaussian) ## fit.AIC = step( fit, model, direction = "both", k = 2, trace = 0) ## regTable = summary.glm( fit.AIC )$coef ################################################### ### code chunk number 22: echo_lind_nnlc_execution (eval = FALSE) ################################################### ## linRes = LocalControl(data = lindner, ## clusterVars = all7Vars, treatmentColName= "abcix", ## outcomeColName = "cardbill", treatmentCode = 1) ## linCI = LocalControlNearestNeighborsConfidence(data = lindner, ## clusterVars = all7Vars, treatmentColName = "abcix", ## outcomeColName = "cardbill", treatmentCode = 1, ## nBootstrap = 100) ## plot(linRes, nnConfidence = linCI) ################################################### ### code chunk number 23: figure_6 ################################################### plot(linRes, nnConfidence = linCI, main = "LocalControl confidence intervals") ################################################### ### code chunk number 24: table_7 ################################################### tas = survSimData[which(survSimData$drug == 1),] tbs = survSimData[which(survSimData$drug == 0),] ot = data.frame(matrix(data = 0, nrow = 3, ncol = 4), row.names = c("n (patients)", "age (years)", "BMI ($\\frac{kg}{m^2}$)")) names(ot) = c("A + B", "A", "B", "$p$~value") ot[1,] = as.numeric(c(nrow(survSimData), nrow(tas),nrow(tbs), NA)) ageTest = t.test(tas$age, tbs$age) ot[2,] = as.numeric(c(mean(survSimData$age), mean(tas$age), mean(tbs$age), as.numeric(ageTest$p.value))) bmiTest = t.test(tas$bmi, tbs$bmi) ot[3,] = as.numeric(c(mean(survSimData$bmi), mean(tas$bmi), mean(tbs$bmi), as.numeric(bmiTest$p.value))) ctext = "Survival simulation cohort summary. A hypothetical hypertension Treatment A (blue) is prescribed more frequently to younger, healthier patients with a low body mass index (BMI), Treatment B (red) is prescribed to older patients with a higher body mass index. Significant treatment biases exist for age and BMI." dm = matrix(data = 2, nrow = 3, ncol = 5) dm[,5] = -2 dm[1,] = 0 print(xtable(ot, caption = ctext, label = "tab:ssimData", digits = dm, align = "lrrrr"), sanitize.text.function = function(x){x}, floating=FALSE, hline.after=NULL, add.to.row=list(pos=list(-1, 0, nrow(ot)), command = c('\\hline\n','\\hline\n','\\hline\n'))) ################################################### ### code chunk number 25: echo_survival_sim_generation (eval = FALSE) ################################################### ## weibullSim = function(N, lambda, rho, betaage, betabmi) { ## bmi = rnorm(N, mean = 26, sd = 4) ## age = runif(N) * 47 + 18 ## pbmi = (bmi - min(bmi)) / (max(bmi) - min(bmi)) * 0.8 + 0.1 ## page = (age - min(age)) / (max(age) - min(age)) * 0.8 + 0.1 ## drug = 1 - rbinom(N, 1, (pbmi + page) / 2) ## et = exp(bmi * betabmi + age * betaage) ## Tlat = (-log(runif(n=N)) / (lambda * et))^(1 / rho) ## C = runif(N) * 30 ## time = pmin(Tlat, C) ## status = as.numeric(Tlat <= C) ## data.frame(id = 1:N, drug, age, bmi, time, status) ## } ## survSimData = weibullSim(10000, 1e-10, 2.6, log(1.2), log(1.45)) ################################################### ### code chunk number 26: echo_call_cardsim (eval = FALSE) ################################################### ## results = LocalControl( data = cardSim, ## outcomeType = "survival", ## treatmentColName = "drug", ## timeColName = "time", ## outcomeColName = "status", ## clusterVars = c("age", "bmi")) ################################################### ### code chunk number 27: figure_7 ################################################### sSim = survSimData maxtime = max(sSim$time) sSim$pc = ifelse(sSim$drug == 1, t1Blue, t0Red) sSim$ci = 1 - sSim$time/maxtime sSim$pca = apply(X = sSim, MARGIN = 1, FUN = function(x) {setColAlpha(col = x["pc"], intensity = x["ci"])}) sSim = sSim[sample(nrow(sSim), nrow(sSim)),] # Main plot setup par(mar = c(2.6,2.6,1,1), mgp = c(1.6, 0.5, 0)) plot(NA,xlim = c(0,30), ylim = c(0,1),xaxs="i", yaxs="i", xlab = "Time in years",cex.lab = 1, cex.axis = 1, cex = 1, ylab = "Survival probability") grid() # Adding survival curves lines(x = sSimFTs, y = sSimT1R1, col = t1Blue, lty = 2, lwd = 2) lines(x = sSimFTs, y = sSimT0R1, col = t0Red, lty = 2, lwd = 2) lines(x = sSimFTs, y = sSimT1R21, col = t1Blue, lwd = 2) lines(x = sSimFTs, y = sSimT0R21, col = t0Red, lwd = 2) box() legend("bottomleft", legend = c("Treatment A (uncorrected)", "Treatment A (corrected)", "Treatment B (uncorrected)", "Treatment B (corrected)"), lty = c(2, 1, 2, 1), lwd = c(2,2,2,2), pt.cex = 1, cex = 1,col = c(t1Blue, t1Blue, t0Red, t0Red)) # Subplot setup # Calculate margin widths in NDC top_margin = grconvertY(1, "lines", "ndc") - grconvertY(0, "lines", "ndc") right_margin = grconvertX(1, "lines", "ndc") - grconvertX(0, "lines", "ndc") # Define normalized dimensions for the subplot based on main plot dimensions u <- par("usr") v <- c( grconvertX(u[1:2], "user", "ndc"), grconvertY(u[3:4], "user", "ndc") ) insetW = 0.5 * (v[2] - v[1]) insetH = 0.5 * (v[4] - v[3]) # Adjust subplot_position to be at upper right subplot_position <- c(v[2] - insetW, v[2], v[4] - insetH, v[4]) # Set parameters for the subplot including margins to show labels par(fig=subplot_position, new=TRUE, mar=c(3, 3, 0, 0)) plot(x = sSim$age, y = sSim$bmi, pch=20, col=sSim$pca, xlab="Age (18-65)", ylab="BMI (10-40)", cex.lab=0.7, cex.axis=0.7, cex=0.7, mgp=c(2, 0.7, 0), axes=FALSE) axis(1, las=1, tck=-0.02, cex.axis=0.7) # Ensure ticks are outward axis(2, las=2, tck=-0.02, cex.axis=0.7) # Ensure ticks are outward box() # Restore default parameters par(fig=c(0, 1, 0, 1), mar=c(5, 5, 4, 4) + 0.1, new=FALSE) ################################################### ### code chunk number 28: table_8 ################################################### data(framingham) smokers = framingham[which(framingham$cursmoke == 1), ] nonsmokers = framingham[which(framingham$cursmoke == 0), ] framVars = c("female", "totchol", "age", "bmi", "BPVar", "heartrte", "glucose") ot = data.frame(matrix(data = 0, nrow = 8, ncol = 4), row.names = c("n (patients)", "female", "totchol ($\\frac{mg}{dL}$)", "age (years)", "BMI ($\\frac{kg}{m^2}$)", "BPVar (mm Hg)", "heartrte (bpm)", "glucose ($\\frac{mg}{dL}$)")) names(ot) = c("All patients", "Smokers", "Nonsmokers", "$p$~value") ot[1,] = as.numeric(c(nrow(framingham), nrow(smokers),nrow(nonsmokers), NA)) for(i in 1:length(framVars)){ cv = framVars[i] ot[i+1,] = as.numeric(c(mean(framingham[, cv]), mean(smokers[, cv]), mean(nonsmokers[, cv]), as.numeric(t.test(smokers[, cv],nonsmokers[, cv])$p.value))) } ctext = 'Framingham Heart Study cohort biases. Patients with preexisting cardiovascular conditions are dropped from the study. Fisher\'s exact test is used for the comparison of the female binary covariate. For the remaining continuous covariates, a t-test is used to compare the two groups. Smoking "treatment" bias significantly affects sex, age, BMI, blood pressure, and heart rate.' dm = matrix(data = 2, nrow = 8, ncol = 5) dm[,5] = -2 dm[1,] = 0 print(xtable(ot, caption = ctext, label = "tab:framData", digits = dm), sanitize.text.function = function(x){x}, table.placement="ht") ################################################### ### code chunk number 29: echo_fram_ana (eval = FALSE) ################################################### ## data("framingham") ## ## framVars = c("female","totchol","age","bmi","BPVar","heartrte","glucose") ## FHSResults = LocalControl( data = framingham, outcomeType = "survival", ## treatmentColName = "cursmoke", treatmentCode = 1, ## timeColName = "time_outcome", outcomeColName = "outcome", ## clusterVars = framVars) ## print(summary(FHSResults)) ################################################### ### code chunk number 30: table_9 ################################################### dm = matrix(data = 2, nrow = 11, ncol = 6) dm[1,] = 0 dm[,6] = -2 print(xtable(FHSSummary[1:11,], digits = dm, caption = "Framingham Local Control summary. Each row corresponds to one radius of correction. The values in the first column are the radius lengths in normalized units. The second column contains the fraction of observations who are informative at the given radius. The pct\\_radius column is the size of the radius as a fraction of the maximum distance between any two observations. The last two columns contain the results from the hypothesis tests comparing the hypertension CIFs for the two treatment groups (as described in Section~\\ref{sec:pepe}).", label = 'tab:framSumm'), include.rownames = T,NA.string = "-", table.placement="hb") ################################################### ### code chunk number 31: figure_8 ################################################### plotz = list() for(radLim in c("rad_1","rad_11")){ t1lines2plot = FHSPlotLines[[radLim]][["t1HT"]] t0lines2plot = FHSPlotLines[[radLim]][["t0HT"]] t1DeathLines = FHSPlotLines[[radLim]][["t1Death"]] t0DeathLines = FHSPlotLines[[radLim]][["t0Death"]] ptitle = ifelse(radLim=="rad_1", "Uncorrected cumulative incidence", "Local Control corrected") plotz[[radLim]]<-ggplot() + theme_bw(base_size = 12) + theme(plot.title = element_text(hjust = 0.5)) + ggtitle(ptitle) + geom_line(data = t1lines2plot, aes(x = TIMES, y=CIF, colour=NG, linetype = NG)) + geom_ribbon(data = t1lines2plot,aes(x = TIMES,ymin=LOWER, ymax=UPPER, fill = NG), alpha=0.2) + geom_line(data = t0lines2plot, aes(x = TIMES, y=CIF, colour=NG, linetype = NG)) + geom_ribbon(data = t0lines2plot,aes(x = TIMES,ymin=LOWER, ymax=UPPER, fill = NG), alpha=0.2) + geom_line(data = t0DeathLines, aes(x = TIMES, y=CIF, colour=NG, linetype = NG)) + geom_ribbon(data = t0DeathLines,aes(x = TIMES,ymin=LOWER, ymax=UPPER, fill = NG), alpha=0.2) + geom_line(data = t1DeathLines, aes(x = TIMES, y=CIF, colour=NG, linetype = NG)) + geom_ribbon(data = t1DeathLines,aes(x = TIMES,ymin=LOWER, ymax=UPPER, fill = NG), alpha=0.2) + scale_fill_manual(values = c( "#377EB8", "#377EB8","#E41A1C","#E41A1C"), name = "Smoker status and outcome", labels = c("Nonsmoker Death", "Nonsmoker Hypertension","Smoker Death", "Smoker Hypertension")) + scale_colour_manual( values = c( "#377EB8", "#377EB8","#E41A1C","#E41A1C"), name = "Smoker status and outcome", labels = c("Nonsmoker Death", "Nonsmoker Hypertension","Smoker Death", "Smoker Hypertension"))+ scale_linetype_manual(name = "Smoker status and outcome", labels = c("Nonsmoker Death", "Nonsmoker Hypertension","Smoker Death", "Smoker Hypertension"), values = c(1, 2, 1, 2)) + scale_x_continuous(limits =c(0,24),expand = c(0, 0), breaks = seq(0, 24, 2), name = "Time (years)") + scale_y_continuous(limits =c(0,.7), expand = c(0, 0), breaks = seq(0, .7, .1), name = "Cumulative incidence of death and hypertension") + theme(legend.justification=c(0,1), legend.position=c(0,1), legend.title=element_blank(), legend.direction='vertical', legend.box='horizontal', legend.background = element_rect(colour = 'transparent', fill = 'transparent')) } grid.arrange(plotz$rad_1, plotz$rad_11, ncol=1) ################################################### ### code chunk number 32: figure_9 ################################################### #sults = linSults lindner$rad_33 = linRad33 # sults$ltds$rad_33 lindner$stent = ifelse(lindner$stent == 1, "1", "0") lindner$female = ifelse(lindner$female == 1, "1", "0") fit <- rpart(formula = rad_33 ~ female + stent, cost = c(0.5, 2), data = lindner, control=rpart.control(maxdepth=2,cp=-1), method = "anova") boxlabs = c("Entire population", "All male", "Male without stent", "Male with stent","All female", "Female with stent", "Female without stent") boxlabs = paste0(boxlabs, "\ncb = ") prp(fit, main=NA, type=4, fallen=T, branch=.3, round=0, leaf.round=0, clip.right.labs=F, under.cex=1, border.col = c("purple", "blue","blue","blue", "#FF0099FF", "#FF0099FF","#FF0099FF"), box.palette=0, prefix = boxlabs, branch.col = c("#FF0099FF", "blue", "blue", "blue", "#FF0099FF", "#FF0099FF"), branch.lwd=3, branch.lty = rep(1,6), # branch.lty = c(2,1,2,1,1,1), extra=1, under=F, lt=" < ", ge=" >= ") par(mar = c(0,0,0,0)) botrowbot = -0.01 botrowtop = 0.165 sidePoints = 5 tbPoints = 8 pcex = 0.7 #MEN WITHOUT STENT BOX box1LX = -0.015 box1RX = 0.1875 points(x = rep(box1LX, sidePoints), y = seq(from = botrowbot, to = botrowtop, length.out = sidePoints) , col = "blue", pch = 1, cex = pcex) #left points(x = rep(box1RX, sidePoints), y = seq(from = botrowbot, to = botrowtop, length.out = sidePoints) , col = "blue", pch = 1, cex = pcex) #right points(x = seq(from = box1LX, to = box1RX, length.out = tbPoints), y = rep(botrowtop, tbPoints) , col = "blue", pch = 1, cex = pcex) #top points(x = seq(from = box1LX, to = box1RX, length.out = tbPoints), y = rep(botrowbot, tbPoints) , col = "blue", pch = 1, cex = pcex) #bot box2LX = 0.273 box2RX = 0.4423 points(x = rep(box2LX, sidePoints), y = seq(from = botrowbot, to = botrowtop, length.out = sidePoints) , col = "blue", pch = 16, cex = pcex) #left points(x = rep(box2RX, sidePoints), y = seq(from = botrowbot, to = botrowtop, length.out = sidePoints) , col = "blue", pch = 16, cex = pcex) #right points(x = seq(from = box2LX, to = box2RX, length.out = tbPoints), y = rep(botrowtop, tbPoints) , col = "blue", pch = 16, cex = pcex) #top points(x = seq(from = box2LX, to = box2RX, length.out = tbPoints), y = rep(botrowbot, tbPoints) , col = "blue", pch = 16, cex = pcex) #bot box3LX = 0.53 box3RX = .73 points(x = rep(box3LX, sidePoints), y = seq(from = botrowbot, to = botrowtop, length.out = sidePoints) , col = "#FF0099FF", pch = 16, cex = pcex) #left points(x = rep(box3RX, sidePoints), y = seq(from = botrowbot, to = botrowtop, length.out = sidePoints) , col = "#FF0099FF", pch = 16, cex = pcex) #right points(x = seq(from = box3LX, to = box3RX, length.out = tbPoints), y = rep(botrowtop, tbPoints) , col = "#FF0099FF", pch = 16, cex = pcex) #top points(x = seq(from = box3LX, to = box3RX, length.out = tbPoints), y = rep(botrowbot, tbPoints) , col = "#FF0099FF", pch = 16, cex = pcex) #bot box4LX = 0.789 box4RX = 1.02 points(x = rep(box4LX, sidePoints), y = seq(from = botrowbot, to = botrowtop, length.out = sidePoints) , col = "#FF0099FF", pch = 1, cex = pcex) #left points(x = rep(box4RX, sidePoints), y = seq(from = botrowbot, to = botrowtop, length.out = sidePoints) , col = "#FF0099FF", pch = 1, cex = pcex) #right points(x = seq(from = box4LX, to = box4RX, length.out = tbPoints), y = rep(botrowtop, tbPoints) , col = "#FF0099FF", pch = 1, cex = pcex) #top points(x = seq(from = box4LX, to = box4RX, length.out = tbPoints), y = rep(botrowbot, tbPoints) , col = "#FF0099FF", pch = 1, cex = pcex) #bot ################################################### ### code chunk number 33: figure_10 ################################################### pMatches = linSummary[nrow(linSummary),] everyone = linSummary[1:(nrow(linSummary)-1),] logXlimits = c(1,min(everyone$pct_radius)) layMatrix = matrix(1, nrow = 20, ncol = 20, byrow = TRUE) layMatrix[1:5,1:18] = 3 layMatrix[6:20,19:20] = 2 layMatrix[1:5,19:20] = 4 layout(layMatrix) #Bottomleft - Main plot par(mar = c(4.1,4.1,0,0)) { plot(x = everyone$pct_radius,y = everyone$ltd, xlim = logXlimits,type = "l", xlab = "Fraction of maximum radius", ylim = c(-16001, 8000), col = "purple", ylab = "Cardbill difference (Abciximab - control)", log = "x", panel.first=grid(equilogs=FALSE), lwd= 2) lines(x = everyone$pct_radius, y = everyone$maleAverage, col = "blue", lwd= 2) lines(x = everyone$pct_radius, y = everyone$femaleAverage, col = "#FF0099FF", lwd= 2) lines(x = everyone$pct_radius,y = everyone$male_without, col = "blue", pch = 1, type = "o") lines(x = everyone$pct_radius,y = everyone$female_without, col = "#FF0099FF", pch = 1, type = "o") lines(x = everyone$pct_radius,y = everyone$male_stent, col = "blue", pch = 16, type = "o") lines(x = everyone$pct_radius,y = everyone$female_stent, col = "#FF0099FF", pch = 16, type = "o") abline(v=linSummary["rad_33","pct_radius"],col="forestgreen", lwd = 3, lty = 3) legend("bottomleft", cex = 1.1, legend = c( "Entire population", "Women", "Men", "Women with stents", "Men with stents", "Women without stents", "Men without stents", "Fraction of informative data","Tree-radius"), col = c("purple", rep(c("#FF0099FF", "blue" ),3), "black","forestgreen"), lwd = c(2,2,2,rep(1,4),2,3), lty = c(rep(1,8),3), pch = c(NA,NA,NA,16,16,1,1,NA,NA)) } #Bottomright - Perfect matches par(mar = c(4.1, 0.5,0,2)) { plot(NA, xaxt = "n", yaxt = "n", ylim = c(-16001, 8000) ,xlim=c(-1,1), xlab = NA, ylab=NA, panel.first = grid(nx = 2, ny = NULL)) axis(side = 1, at = 0) mtext(side = 1, text = "Perfect\n matches", cex = 0.6, line = 3) lines(x = c(-1/2,1/2),y = c(pMatches$ltd,pMatches$ltd), col = "purple", lwd = 2) lines(x = c(-1/2,1/2),y = c(pMatches$maleAverage,pMatches$maleAverage), col = "blue", lwd = 2) lines(x = c(-1/2,1/2),y = c(pMatches$femaleAverage,pMatches$femaleAverage), col = "#FF0099FF", lwd = 2) points(x = c(-1/2,0,1/2),y = rep(pMatches$male_without,3), col = "blue", pch = 1, type = "o") points(x = c(-1/2,0,1/2),y = rep(pMatches$female_without,3), col = "#FF0099FF", pch = 1, type = "o") points(x = c(-1/2,0,1/2),y = rep(pMatches$male_stent,3), col = "blue", pch = 16, type = "o") points(x = c(-1/2,0,1/2),y = rep(pMatches$female_stent,3), col = "#FF0099FF", pch = 16, type = "o") } #Topleft - pct info for main plot par(mar = c(0.5, 4.1, 2, 0)) { plot(x = everyone$pct_radius, y = everyone$pct_data, xaxt = "n",yaxt="n",type = "l", lwd = 3, ylab =NA , xlab=NA, xlim = logXlimits, ylim = c(0,1), log="x", panel.first=grid(ny = 5,nx=NULL,equilogs = F)) abline(v=linSummary["rad_33","pct_radius"],col="forestgreen", lwd = 3, lty = 3) mtext(side = 2, text = "Fraction\ninformative data", cex = 0.6, line = 2) axis(side = 2, at = c(0,0.2,0.4,0.6,0.8,1), labels = c(0,".2",".4",".6",".8","1"),cex.axis = .9) } #Topright - pct info for perfect matches par(mar = c(0.5, 0.5, 2, 2)) { plot(x=c(-1/2,1/2), y= c(pMatches$pct_data, pMatches$pct_data), type = "l", xlim = c(-1,1), ylim = c(0,1), xaxt = "n", lwd = 3, yaxt = "n", xlab = NA, ylab=NA, panel.first=grid(ny = 5,nx=2)) }