################################################################################ # Standalone script for: # # General Semiparametric Shared Frailty Model: # Estimation and Simulation with frailtySurv # # Authors: # John V. Monaco, Malka Gorfine, Li Hsu # # **Note** # This script is very computationally demanding: each section takes about a day. # # Some results are saved as RData files and then loaded when the # script is re-executed. This includes the case studies and # simulation results. # # To run as a background script and log the results: # # $ nohup Rscript code.R > log 2> log.err < /dev/null & # ################################################################################ library(grid) library(gtable) library(ggplot2) library(gridExtra) library(latex2exp) library(plyr) library(reshape2) library(survival) library(frailtypack) library(frailtySurv) SEED <- 2015 # Seed used for data gen and bootstrap functions CORES <- 0 # number of cores to use, set to 0 for all R.Version()$version.string packageDescription("frailtySurv", fields = "Version") # Helper function to load/save results to avoid re-running long simulations load_or_run <- function(name, fn, args) { fname <- paste(name, ".RData", sep = "") if (file.exists(fname)) { cat("Loading", name, "...\n") load(fname) } else { cat("Running", name, "...\n") assign(name, do.call(fn, args)) save(list = name, file = fname) } get(name) } ################################################################################ # Section 1: Introduction ################################################################################ # Figure 1 pdf(NULL) # To avoid generating Rplots.pdf palette <- c("#009E73", "#E69F00", "#0072B2", "#D55E00", "#CC79A7") names(palette) <- c("gamma", "lognormal", "invgauss", "pvf", "posstab") tau1 <- 0.3 tau2 <- 0.45 tau3 <- 0.70 tau1.theta.posstab <- 1 - tau1 tau1.theta.gamma <- frailtySurv:::theta.given.tau(tau1, "gamma") tau1.theta.pvf <- frailtySurv:::theta.given.tau(tau1, "pvf") tau1.theta.lognormal <- frailtySurv:::theta.given.tau(tau1, "lognormal") tau1.theta.invgauss <- frailtySurv:::theta.given.tau(tau1, "invgauss") tau2.theta.posstab <- 1 - tau2 tau2.theta.gamma <- frailtySurv:::theta.given.tau(tau2, "gamma") tau2.theta.lognormal <- frailtySurv:::theta.given.tau(tau2, "lognormal") tau2.theta.invgauss <- frailtySurv:::theta.given.tau(tau2, "invgauss") tau3.theta.posstab <- 1 - tau3 tau3.theta.gamma <- frailtySurv:::theta.given.tau(tau3, "gamma") tau3.theta.lognormal <- frailtySurv:::theta.given.tau(tau3, "lognormal") x <- seq(0, 3, length.out = 1000) x.pvf <- x[x > 0.01] x.ps1 <- seq(0.15, 3, length.out = 1000) x.ps2 <- seq(0.05, 3, length.out = 1000) x.ig <- c(seq(0, 0.01, length.out = 1000), seq(0.01, 3, length.out = 1000)) tau1.densities <- rbind( data.frame( x = x, y = frailtySurv:::dgamma_c(x, tau1.theta.gamma), distr = "gamma" ), data.frame( x = x, y = frailtySurv:::dlognormal_c(x, tau1.theta.lognormal), distr = "lognormal" ), data.frame( x = x, y = frailtySurv:::dinvgauss_c(x, tau1.theta.invgauss), distr = "invgauss" ), data.frame( x = x.pvf, y = frailtySurv:::dpvf_c(x.pvf, tau1.theta.pvf), distr = "pvf" ), data.frame( x = x.ps1, y = frailtySurv:::dposstab_c(x.ps1, tau1.theta.posstab), distr = "posstab" ) ) tau2.densities <- rbind( data.frame( x = x, y = frailtySurv:::dgamma_c(x, tau2.theta.gamma), distr = "gamma" ), data.frame( x = x, y = frailtySurv:::dlognormal_c(x, tau2.theta.lognormal), distr = "lognormal" ), data.frame( x = x.ig, y = frailtySurv:::dinvgauss_c(x.ig, tau2.theta.invgauss), distr = "invgauss" ), data.frame( x = x.ps2, y = frailtySurv:::dposstab_c(x.ps2, tau2.theta.posstab), distr = "posstab" ) ) tau3.densities <- rbind( data.frame( x = x, y = frailtySurv:::dgamma_c(x, tau3.theta.gamma), distr = "gamma" ), data.frame( x = x, y = frailtySurv:::dlognormal_c(x, tau3.theta.lognormal), distr = "lognormal" ), data.frame( x = x, y = frailtySurv:::dposstab_c(x, tau3.theta.posstab), distr = "posstab" ) ) colscale <- scale_colour_manual(values = palette, breaks = c("gamma", "lognormal", "invgauss", "pvf", "ps"), labels = c("Gamma", "LN", "IG", "PVF", "PS")) p1 <- ggplot(tau1.densities, aes( x = x, y = y, color = distr, group = distr )) + geom_line() + scale_colour_manual(values = palette, breaks = c("gamma", "lognormal", "invgauss", "pvf", "posstab"), labels = c(sprintf("Gamma(%.2f)", tau1.theta.gamma), sprintf("LN(%.2f)", tau1.theta.lognormal), sprintf("IG(%.2f)", tau1.theta.invgauss), sprintf("PVF(%.2f)", tau1.theta.pvf), sprintf("PS(%.2f)", tau1.theta.posstab))) + xlab(TeX("$\\omega$")) + ylab("Density") + ylim(0, 1.5) + theme( plot.margin = unit(c(0,-1, 0, 0), "lines"), plot.background = element_blank(), legend.position = c(1, 1), legend.justification = c(1, 1), legend.title = element_blank() ) + ggtitle(TeX(sprintf("$ \\kappa = %.2f $", tau1))) + theme( plot.title = element_text( family = "Times", size = 12, hjust = 0.5 ), legend.position = c(0.94, 0.94) ) p2 <- ggplot(tau2.densities, aes( x = x, y = y, color = distr, group = distr )) + geom_line() + scale_colour_manual(values = palette, breaks = c("gamma", "lognormal", "invgauss", "posstab"), labels = c(sprintf("Gamma(%.2f)", tau2.theta.gamma), sprintf("LN(%.2f)", tau2.theta.lognormal), sprintf("IG(%.2f)", tau2.theta.invgauss), sprintf("PS(%.2f)", tau2.theta.posstab))) + xlab(TeX("$\\omega$")) + ylim(0, 1.5) + theme( plot.margin = unit(c(0,-0.5, 0,-0.5), "lines"), axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank(), plot.background = element_blank(), legend.position = c(1, 1), legend.justification = c(1, 1), legend.title = element_blank() ) + ggtitle(TeX(sprintf("$ \\kappa = %.2f $", tau2))) + theme( plot.title = element_text( family = "Times", size = 12, hjust = 0.5 ), legend.position = c(0.94, 0.94) ) p3 <- ggplot(tau3.densities, aes( x = x, y = y, color = distr, group = distr )) + geom_line() + scale_colour_manual(values = palette, breaks = c("gamma", "lognormal", "posstab"), labels = c(sprintf("Gamma(%.2f)", tau3.theta.gamma), sprintf("LN(%.2f)", tau3.theta.lognormal), sprintf("PS(%.2f)", tau3.theta.posstab))) + xlab(TeX("$\\omega$")) + ylim(0, 1.5) + theme( plot.margin = unit(c(0, 0, 0,-1), "lines"), axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank(), plot.background = element_blank(), legend.position = c(1, 1), legend.justification = c(1, 1), legend.title = element_blank() ) + ggtitle(TeX(sprintf("$ \\kappa = %.2f $", tau3))) + theme( plot.title = element_text( family = "Times", size = 12, hjust = 0.5 ), legend.position = c(0.94, 0.94) ) gt1 <- ggplot_gtable(ggplot_build(p1)) gt2 <- ggplot_gtable(ggplot_build(p2)) gt3 <- ggplot_gtable(ggplot_build(p3)) newWidth = unit.pmax(gt1$widths[2:3], gt2$widths[2:3], gt3$widths[2:3]) gt1$widths[2:3] = as.list(newWidth) gt2$widths[2:3] = as.list(newWidth) gt3$widths[2:3] = as.list(newWidth) pdf("figure1.pdf", width = 7, height = 3) print(grid.arrange(gt1, gt2, gt3, ncol = 3)) dev.off() ################################################################################ # Section 2: Data generation ################################################################################ # Specify the baseline hazard set.seed(SEED) dat <- genfrail( N = 300, K = 2, beta = c(log(2), log(3)), frailty = "gamma", theta = 2, lambda_0 = function(t, tau = 4.6, C = 0.01) (tau * (C * t) ^ tau) / t ) head(dat, 3) # Specify the cumulative baseline hazard set.seed(SEED) dat.cbh <- genfrail( N = 300, K = 2, beta = c(log(2), log(3)), frailty = "gamma", theta = 2, Lambda_0 = function(t, tau = 4.6, C = 0.01) (C * t) ^ tau ) head(dat.cbh, 3) # Specify the inverse cumulative baseline hazard set.seed(SEED) dat.inv <- genfrail( N = 300, K = 2, beta = c(log(2), log(3)), frailty = "gamma", theta = 2, Lambda_0_inv = function(t, tau = 4.6, C = 0.01) (t ^ (1 / tau)) / C ) head(dat.inv, 3) # PVF frailty set.seed(SEED) dat.pvf <- genfrail( N = 300, K = 2, beta = c(log(2), log(3)), frailty = "pvf", theta = 0.3, censor.rate = 0.40, Lambda_0_inv = function(t, tau = 4.6, C = 0.01) (t ^ (1 / tau)) / C ) summary(dat.pvf) ################################################################################ # Section 3: Model estimation ################################################################################ if (file.exists("fit-example.RData")) { load("fit-example.RData") } else { # Fit a model by maximizing the log-likelihood fit <- fitfrail( Surv(time, status) ~ Z1 + Z2 + cluster(family), dat, frailty = "gamma", fitmethod = "loglik" ) # Fit a model by solving the score equations fit.score <- fitfrail( Surv(time, status) ~ Z1 + Z2 + cluster(family), dat, frailty = "gamma", fitmethod = "score" ) COV.est <- vcov(fit) set.seed(SEED) COV.boot <- vcov(fit, boot = TRUE, B = 500) save(fit, COV.est, COV.boot, fit.score, file = "fit-example.RData") } fit fit.score # Summarize the survival curve head(summary(fit), 3) tail(summary(fit), 3) # Summarize the baseline hazard at specific time points summary(fit, type = "cumhaz", Lambda.times = c(20, 50, 80, 110)) # Estimated variance sqrt(diag(COV.est)) # Weighted-bootstrap variance sqrt(diag(COV.boot))[1:8] ################################################################################ # Section 4: Simulation ################################################################################ if (file.exists("sim1.example.RData")) { load("sim1.example.RData") } else { # Note that 1000 reps are used for the simulation in the main text set.seed(SEED) sim <- simfrail( 1000, genfrail.args = alist( beta = c(log(2), log(3)), frailty = "gamma", censor.rate = 0.30, N = 300, K = 2, theta = 2, covar.distr = "uniform", covar.param = c(0, 1), Lambda_0 = function(t, tau = 4.6, C = 0.01) (C * t) ^ tau ), fitfrail.args = alist( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), frailty = "gamma" ), Lambda.times = 1:120, cores = CORES ) # Using coxph set.seed(SEED) sim.coxph <- frailtySurv:::simcoxph( 1000, genfrail.args = alist( beta = c(log(2), log(3)), frailty = "gamma", censor.rate = 0.30, N = 300, K = 2, theta = 2, covar.distr = "uniform", covar.param = c(0, 1), Lambda_0 = function(t, tau = 4.6, C = 0.01) (C * t) ^ tau ), coxph.args = alist(formula = Surv(time, status) ~ Z1 + Z2 + frailty(family)), Lambda.times = 1:120 ) save(sim, sim.coxph, file = "sim1.example.RData") } summary(sim) summary(sim.coxph) # Correlation between coefficients and frailty parameter sapply(names(sim)[grepl("^hat.beta|^hat.theta", names(sim))], function(name) cor(sim[[name]], sim.coxph[[name]])) # Mean correlation between cumulative baseline hazard mean(sapply(names(sim)[grepl("^hat.Lambda", names(sim))], function(name) cor(sim[[name]], sim.coxph[[name]])), na.rm = TRUE) ################################################################################ # Section 5: Case study ################################################################################ if (file.exists("drs.RData")) { load("drs.RData") } else { # Diabetic Retinopathy Study data("drs", package = "frailtySurv") fit.drs <- fitfrail(Surv(time, status) ~ treated + cluster(subject_id), drs, frailty = "gamma") COV.drs <- vcov(fit.drs) ## Warning: computationally expensive set.seed(SEED) COV.boot.drs <- vcov(fit.drs, boot = TRUE) save(fit.drs, COV.drs, COV.boot.drs, file = "drs.RData") } fit.drs sqrt(diag(COV.drs)) # Parameter estimation trace pdf("figure2.pdf", width = 8, height = 3) plot(fit.drs, type = "trace") dev.off() set.seed(SEED) # Seed the weights generated by the bootstrap procedure pdf("figure3.pdf", width = 8, height = 3) plot(fit.drs, type = "cumhaz", CI = 0.95) dev.off() # Compare to coxph estimates library("survival") coxph(Surv(time, status) ~ treated + frailty(subject_id), drs) if (file.exists("hdfail.RData")) { load("hdfail.RData") } else { # Hard drive failure data("hdfail", package = "frailtySurv") hdfail.sub <- subset(hdfail, grepl("WDC", model)) fit.hd <- fitfrail( Surv(time, status) ~ temp + rer + rsc + psc + cluster(model), hdfail.sub, frailty = "gamma", fitmethod = "score" ) ## Warning: computationally expensive weighted bootstrap procedure set.seed(SEED) COV <- vcov(fit.hd, boot = TRUE, cores=1) save(fit.hd, COV, file = "hdfail.RData") } fit.hd set.seed(SEED) COV <- vcov(fit.hd, boot = TRUE, cores=1) # Result should be cached # Coefficients standard errors and pvalues se <- sqrt(diag(COV)[c("temp", "rsc", "rer", "psc", "theta.1")]) se pvalues <- pnorm(abs(c(fit.hd$beta, fit.hd$theta)) / se, lower.tail = FALSE) * 2 pvalues # This should run quickly since the bootstraped variance was already calculated # above and cached in the fit.hd object set.seed(SEED) pdf("figure4.pdf", width = 8, height = 3) plot(fit.hd, type = "cumhaz", CI = 0.95, end = 365 * 6) dev.off() ################################################################################ # Appendix B: Simulation results ################################################################################ # Default simulation parameters unless otherwise noted REPS <- 1000 # Number of repetitions for each simulation LAMBDA.TIMES <- 1:120 # where to evaluate the baseline hazard N <- 300 # Number of clusters K <- 2 # default cluster sizes BETA <- c(log(2), log(3)) THETA <- 2 THETA.PVF <- 0.3 FRAILTY <- "gamma" CENSOR.PARAM <- c(130, 15) CENSOR.RATE <- 0.30 COVAR.DISTR <- "uniform" # Simple baseline hazard lambda_0 <- function(t, tau = 4.6, C = 0.01) { (tau * (C * t) ^ tau) / t } Lambda_0 <- function(t, tau = 4.6, C = 0.01) { (C * t) ^ tau } # Increasing oscillating baseline hazard lambda_0.osc <- function(t, tau = 4.6, C = 0.01, A = 2, f = 0.1) { A ^ sin(f * pi * t) * (tau * (C * t) ^ tau) / t } Lambda_0.osc <- Vectorize(function(t, tau = 4.6, C = 0.01, A = 2, f = 0.1) { ifelse(t > 0, integrate(lambda_0.osc, 0, t, subdivisions = 1000L)$value, 0) }) # Load previous simulation results or save new results and print summary/figure simfrail.output <- function(name, figname, seed = SEED, ...) { if (file.exists(paste(name, ".RData", sep = ""))) { cat("Loading", name, "...\n") load(paste(name, ".RData", sep = "")) } else { cat("Running", name, "...\n") set.seed(seed) assign(name, simfrail(...)) save(list = name, file = paste(name, ".RData", sep = "")) } sim <- get(name) print(summary(sim)) pdf(paste(figname, ".pdf", sep = ""), width = 8, height = 3) print(plot(sim, "cumhaz")) dev.off() invisible() } # Section B.1 simfrail.output( "sim1", "figure5", reps = REPS, genfrail.args = alist( beta = BETA, frailty = "gamma", censor.rate = CENSOR.RATE, N = N, K = K, theta = THETA, covar.distr = COVAR.DISTR, Lambda_0 = Lambda_0 ), fitfrail.args = alist( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), frailty = "gamma" ), Lambda.times = LAMBDA.TIMES, cores = CORES ) # Extended version of sim1, increasing N to show decreasing residuals if (file.exists("sim1.ext.RData")) { cat("Loading sim1.ext...\n") load("sim1.ext.RData") } else { cat("Running sim1.ext...\n") sim1.ext <- frailtySurv:::simfrail.enum( reps = 100, # Only need enough to estimate standard error seed = SEED, genfrail.args = alist( beta = BETA, frailty = "gamma", censor.rate = CENSOR.RATE, K = K, theta = THETA, covar.distr = COVAR.DISTR, Lambda_0 = Lambda_0 ), fitfrail.args = alist( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), frailty = "gamma" ), Lambda.times = LAMBDA.TIMES, param.name = "N", param.values = c(25, 50, 250, 500), cores = CORES ) save(sim1.ext, file = "sim1.ext.RData") } pdf("figure6.pdf", width = 8, height = 3) plot(sim1.ext, "residuals", Lambda.times = c(30, 60, 90)) # Section B.2 simfrail.output( "sim2", "figure7", reps = REPS, genfrail.args = alist( beta = BETA, frailty = "gamma", censor.rate = CENSOR.RATE, N = 100, K = 6, theta = THETA, covar.distr = COVAR.DISTR, Lambda_0 = Lambda_0 ), fitfrail.args = alist( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), frailty = "gamma" ), Lambda.times = LAMBDA.TIMES, cores = CORES ) # Section B.3 simfrail.output( "sim3", "figure8", reps = REPS, genfrail.args = alist( beta = BETA, frailty = "gamma", censor.rate = CENSOR.RATE, N = N, K = K, theta = THETA, covar.distr = COVAR.DISTR, Lambda_0 = Lambda_0, round.base = 10 ), fitfrail.args = alist( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), frailty = "gamma" ), Lambda.times = LAMBDA.TIMES, cores = CORES ) # Section B.4 simfrail.output( "sim4", "figure9", reps = REPS, genfrail.args = alist( beta = BETA, frailty = "gamma", censor.rate = CENSOR.RATE, N = N, K = K, theta = THETA, covar.distr = COVAR.DISTR, Lambda_0 = Lambda_0.osc ), fitfrail.args = alist( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), frailty = "gamma" ), Lambda.times = LAMBDA.TIMES, cores = CORES ) # Section B.5 simfrail.output( "sim5", "figure10", reps = REPS, genfrail.args = alist( beta = BETA, frailty = "pvf", censor.rate = CENSOR.RATE, N = N, K = K, theta = THETA.PVF, covar.distr = COVAR.DISTR, Lambda_0 = Lambda_0 ), fitfrail.args = alist( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), frailty = "pvf" ), Lambda.times = LAMBDA.TIMES, cores = CORES ) # Section B.6 simfrail.output( "sim6", "figure11", reps = REPS, genfrail.args = alist( beta = BETA, frailty = "pvf", censor.rate = CENSOR.RATE, N = N, K = "poisson", K.param = c(2, 0), theta = THETA.PVF, covar.distr = COVAR.DISTR, Lambda_0 = Lambda_0 ), fitfrail.args = alist( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), frailty = "pvf" ), Lambda.times = LAMBDA.TIMES, cores = CORES ) # Section B.7 simfrail.output( "sim7", "figure12", reps = REPS, genfrail.args = alist( beta = BETA, frailty = "lognormal", censor.rate = CENSOR.RATE, N = N, K = K, theta = THETA, covar.distr = COVAR.DISTR, Lambda_0 = Lambda_0 ), fitfrail.args = alist( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), frailty = "lognormal" ), Lambda.times = LAMBDA.TIMES, cores = CORES ) # Section B.8 simfrail.output( "sim8", "figure13", reps = REPS, genfrail.args = alist( beta = BETA, frailty = "invgauss", censor.rate = CENSOR.RATE, N = N, K = K, theta = THETA, covar.distr = COVAR.DISTR, Lambda_0 = Lambda_0 ), fitfrail.args = alist( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), frailty = "invgauss" ), Lambda.times = LAMBDA.TIMES, cores = CORES ) ################################################################################ # Appendix C: Performance analysis ################################################################################ FITMETHOD <- "score" REPS <- 100 N <- 300 K <- 2 BETA <- c(log(2), log(3)) KENDALLS.TAU <- 0.3 FRAILTY.ALL <- c("gamma", "pvf", "lognormal", "invgauss") FRAILTY.INT <- c("lognormal", "invgauss") lambda_0 <- function(t, c = 0.01, d = 4.6) { (d * (c * t) ^ d) / t } Lambda_0 <- function(t, c = 0.01, d = 4.6) { (c * t) ^ d } Lambda_0_inv <- function(t, c = 0.01, d = 4.6) { (t ^ (1 / d)) / c } TOLS <- c(1, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9) # Function to compare the timings of core frailtySurv functions compare.frailtySurv.timings <- function(reps, N) { params <- expand.grid( seed = sample(1:1e7, reps, replace = TRUE), frailty = FRAILTY.ALL, N = N, stringsAsFactors = FALSE ) fn <- function(p) { set.seed(p[["seed"]]) genfrail.bh.time <- system.time(dat.bh <- do.call( "genfrail", list( N = as.integer(p[["N"]]), K = K, beta = BETA, frailty = p[["frailty"]], theta = frailtySurv:::theta.given.tau(KENDALLS.TAU, p[["frailty"]]), lambda_0 = lambda_0 ) ))[["elapsed"]] set.seed(p[["seed"]]) genfrail.cbh.time <- system.time(dat.cbh <- do.call( "genfrail", list( N = as.integer(p[["N"]]), K = K, beta = BETA, frailty = p[["frailty"]], theta = frailtySurv:::theta.given.tau(KENDALLS.TAU, p[["frailty"]]), Lambda_0 = Lambda_0 ) ))[["elapsed"]] set.seed(p[["seed"]]) genfrail.icbh.time <- system.time(dat.icbh <- do.call( "genfrail", list( N = as.integer(p[["N"]]), K = K, beta = BETA, frailty = p[["frailty"]], theta = frailtySurv:::theta.given.tau(KENDALLS.TAU, p[["frailty"]]), Lambda_0_inv = Lambda_0_inv ) ))[["elapsed"]] set.seed(p[["seed"]]) fitfrail.loglik.time <- system.time(fit.loglik <- do.call( "fitfrail", list( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), dat = dat.icbh, frailty = p[["frailty"]], fitmethod = "loglik" ) ))[["elapsed"]] set.seed(p[["seed"]]) fitfrail.score.time <- system.time(fit.score <- do.call( "fitfrail", list( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), dat = dat.icbh, frailty = p[["frailty"]], fitmethod = "score" ) ))[["elapsed"]] set.seed(p[["seed"]]) vcov.time <- system.time(v <- do.call("vcov", list(fit.score)))[["elapsed"]] c( seed = p[["seed"]], N = as.integer(p[["N"]]), frailty = p[["frailty"]], genfrail.bh.time = genfrail.bh.time, genfrail.cbh.time = genfrail.cbh.time, genfrail.icbh.time = genfrail.icbh.time, fitfrail.loglik.time = fitfrail.loglik.time, fitfrail.score.time = fitfrail.score.time, vcov.time = vcov.time ) } timings <- parallel::mclapply( split(params, 1:nrow(params)), fn, mc.preschedule = FALSE, mc.set.seed = TRUE, mc.cores = parallel::detectCores(), mc.silent = FALSE ) data.frame(t(simplify2array(timings))) } # function to compare the timings of three different shared frailty packages compare.package.timings <- function(reps, N) { params <- expand.grid( seed = sample(1:1e7, reps, replace = TRUE), frailty = c("gamma", "lognormal"), N = N, stringsAsFactors = FALSE ) fn <- function(p) { if (p[["frailty"]] == "gamma") { coxph.frailty <- "gamma" frailtyPenal.frailty <- "Gamma" } else if (p[["frailty"]] == "lognormal") { coxph.frailty <- "gaussian" frailtyPenal.frailty <- "LogN" } set.seed(p[["seed"]]) dat <- genfrail( N = as.integer(p[["N"]]), K = K, beta = BETA, frailty = p[["frailty"]], theta = frailtySurv:::theta.given.tau(KENDALLS.TAU, p[["frailty"]]), Lambda_0_inv = Lambda_0_inv ) set.seed(p[["seed"]]) fitfrail.time <- system.time(fit.loglik <- do.call( "fitfrail", list( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), dat = dat, frailty = p[["frailty"]], fitmethod = "score" ) ))[["elapsed"]] set.seed(p[["seed"]]) coxph.time <- system.time(fit.score <- do.call("coxph", list( formula = Surv(time, status) ~ Z1 + Z2 + frailty(family, coxph.frailty), data = dat )))[["elapsed"]] set.seed(p[["seed"]]) frailtyPenal.time <- system.time(fit.score <- do.call( "frailtyPenal", list( formula = Surv(time, status) ~ Z1 + Z2 + cluster(family), data = dat, n.knots = 10, kappa = 2, RandDist = frailtyPenal.frailty ) ))[["elapsed"]] c( seed = p[["seed"]], N = as.integer(p[["N"]]), frailty = p[["frailty"]], fitfrail.time = fitfrail.time, coxph.time = coxph.time, frailtyPenal.time = frailtyPenal.time ) } timings <- parallel::mclapply( split(params, 1:nrow(params)), fn, mc.preschedule = FALSE, mc.set.seed = TRUE, mc.cores = parallel::detectCores(), mc.silent = FALSE ) data.frame(t(simplify2array(timings))) } # function to compare the timings of three different shared frailty packages compare.fitfrail.accuracy <- function(reps, frailty, fitfrail.param, fitfrail.param.values) { params <- expand.grid( seed = sample(1:1e7, reps, replace = TRUE), frailty = frailty, tol = fitfrail.param.values, stringsAsFactors = FALSE ) fn <- function(p) { set.seed(p[["seed"]]) dat <- genfrail( N = N, K = K, beta = BETA[1], frailty = p[["frailty"]], theta = frailtySurv:::theta.given.tau(KENDALLS.TAU, p[["frailty"]]), Lambda_0_inv = Lambda_0_inv ) set.seed(p[["seed"]]) fitfrail.args <- list( formula = Surv(time, status) ~ Z1 + cluster(family), dat = dat, frailty = p[["frailty"]], fitmethod = FITMETHOD ) fitfrail.args[[fitfrail.param]] <- as.numeric(p[["tol"]]) if (fitfrail.param == "abstol") { fitfrail.args[["reltol"]] <- 0 } else if (fitfrail.param == "reltol") { fitfrail.args[["abstol"]] <- 0 } else if (fitfrail.param == "int.abstol") { fitfrail.args[["int.reltol"]] <- 0 } else if (fitfrail.param == "int.reltol") { fitfrail.args[["int.abstol"]] <- 0 } fitfrail.time <- system.time(fit <- do.call("fitfrail", fitfrail.args))[["elapsed"]] fitfrail.beta.1.res <- (fit$beta[1] - attr(dat, "beta")[1]) fitfrail.theta.res <- (fit$theta - attr(dat, "theta")) c( seed = p[["seed"]], frailty = p[["frailty"]], tol = p[["tol"]], fitfrail.time = fitfrail.time, fitfrail.beta.1.res = fitfrail.beta.1.res, fitfrail.theta.res = fitfrail.theta.res ) } timings <- parallel::mclapply( split(params, 1:nrow(params)), fn, mc.preschedule = FALSE, mc.set.seed = TRUE, mc.cores = parallel::detectCores(), mc.silent = FALSE ) data.frame(t(simplify2array(timings))) } # Section C.1 set.seed(SEED) frailtySurv.timings <- load_or_run("frailtySurv.timings", compare.frailtySurv.timings, list(reps = REPS, N = seq(50, 200, 10))) # Figure 14 pdf(NULL) indx <- sapply(frailtySurv.timings, is.factor) indx[["frailty"]] <- FALSE frailtySurv.timings[indx] <- lapply(frailtySurv.timings[indx], function(x) as.numeric(as.character(x))) xcenter <- min(frailtySurv.timings$N) + (max(frailtySurv.timings$N) - min(frailtySurv.timings$N)) / 2 ymax <- max(frailtySurv.timings$genfrail.bh.time) + sd(frailtySurv.timings$genfrail.bh.time) colscale <- scale_colour_manual(values = palette, breaks = c("gamma", "lognormal", "invgauss", "pvf"), labels = c("Gamma", "LN", "IG", "PVF")) p1 <- ggplot(frailtySurv.timings, aes(x = N, y = genfrail.bh.time, color = frailty)) + colscale + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylab("Runtime (s)") + ylim(0, ymax) + theme( plot.margin = unit(c(0, 0, 0.2, 0), "lines"), plot.background = element_blank(), legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank() ) + annotate( "text", x = xcenter, y = ymax, label = "Baseline hazard", hjust = 0.5, vjust = 1 ) ymax <- max(frailtySurv.timings$genfrail.cbh.time) + sd(frailtySurv.timings$genfrail.cbh.time) p2 <- ggplot(frailtySurv.timings, aes(x = N, y = genfrail.cbh.time, color = frailty)) + colscale + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylab("Runtime (s)") + ylim(0, ymax) + theme( plot.margin = unit(c(0.2, 0, 0.2, 0), "lines"), plot.background = element_blank(), legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank() ) + annotate( "text", x = xcenter, y = ymax, label = "Cumulative baseline hazard", hjust = 0.5, vjust = 1 ) ymax <- max(frailtySurv.timings$genfrail.icbh.time) + sd(frailtySurv.timings$genfrail.icbh.time) p3 <- ggplot(frailtySurv.timings, aes(x = N, y = genfrail.icbh.time, color = frailty)) + colscale + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + xlab("N") + ylab("Runtime (s)") + ylim(0, ymax) + theme( legend.position = "bottom", legend.justification = c(0.5, 1), legend.direction = "horizontal", legend.title = element_blank(), plot.margin = unit(c(0.2, 0, 0, 0), "lines"), plot.background = element_blank() ) + annotate( "text", x = xcenter, y = ymax, label = "Inverse cumulative baseline hazard", hjust = 0.5, vjust = 1 ) gt1 <- ggplotGrob(p1) gt2 <- ggplotGrob(p2) gt3 <- ggplotGrob(p3) pdf("figure14.pdf", width = 7, height = 6) print(grid.arrange(gtable_rbind(gt1, gt2, gt3))) dev.off() # Figure 15 pdf(NULL) xcenter <- min(frailtySurv.timings$N) + (max(frailtySurv.timings$N) - min(frailtySurv.timings$N)) / 2 ymax <- 75 colscale <- scale_colour_manual(values = palette, breaks = c("gamma", "lognormal", "invgauss", "pvf"), labels = c("Gamma", "LN", "IG", "PVF")) p1 <- ggplot(frailtySurv.timings, aes(x = N, y = fitfrail.loglik.time, color = frailty)) + colscale + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylab("Runtime (s)") + ylim(0, ymax) + theme(plot.margin = unit(c(0, 0, 0.2, 0), "lines"), plot.background = element_blank()) + theme( plot.title = element_text(family = "Times", size = 12), legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank() ) + annotate( "text", x = xcenter, y = ymax, label = "Loglikelihood", hjust = 0.5, vjust = 1 ) ymax <- 75 p2 <- ggplot(frailtySurv.timings, aes(x = N, y = fitfrail.score.time, color = frailty)) + colscale + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + xlab("N") + ylab("Runtime (s)") + ylim(0, ymax) + theme( plot.margin = unit(c(0.2, 0, 0, 0), "lines"), plot.background = element_blank(), legend.position = "bottom", legend.justification = c(0.5, 1), legend.direction = "horizontal", legend.title = element_blank() ) + annotate( "text", x = xcenter, y = ymax, label = "Score", hjust = 0.5, vjust = 1 ) gt1 <- ggplotGrob(p1) gt2 <- ggplotGrob(p2) pdf("figure15.pdf", width = 7, height = 4) print(grid.arrange(gtable_rbind(gt1, gt2))) dev.off() # Figure 16 pdf(NULL) xcenter <- min(frailtySurv.timings$N) + (max(frailtySurv.timings$N) - min(frailtySurv.timings$N)) / 2 ymax <- max(frailtySurv.timings$vcov.time) + sd(frailtySurv.timings$vcov.time) colscale <- scale_colour_manual(values = palette, breaks = c("gamma", "lognormal", "invgauss", "pvf"), labels = c("Gamma", "LN", "IG", "PVF")) p1 <- ggplot(frailtySurv.timings, aes(x = N, y = vcov.time, color = frailty)) + colscale + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + xlab("N") + ylab("Runtime (s)") + ylim(0, ymax) + theme( plot.margin = unit(c(0.2, 0, 0, 0), "lines"), plot.background = element_blank(), legend.position = "bottom", legend.justification = c(0.5, 1), legend.direction = "horizontal", legend.title = element_blank() ) gt1 <- ggplotGrob(p1) pdf("figure16.pdf", width = 7, height = 2.5) print(grid.arrange(gtable_rbind(gt1))) dev.off() # Section C.2 set.seed(SEED) package.timings <- load_or_run("package.timings", compare.package.timings, list(reps = REPS, N = seq(50, 200, 10))) # Figure 17 pdf(NULL) indx <- sapply(package.timings, is.factor) indx[["frailty"]] <- FALSE package.timings[indx] <- lapply(package.timings[indx], function(x) as.numeric(as.character(x))) xcenter <- min(package.timings$N) + (max(package.timings$N) - min(package.timings$N)) / 2 ymax <- 75 colscale <- scale_colour_manual(values = palette, breaks = c("gamma", "lognormal"), labels = c("Gamma", "LN")) p1 <- ggplot(package.timings, aes(x = N, y = fitfrail.time, color = frailty)) + colscale + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylab("Runtime (s)") + ylim(0, ymax) + theme(plot.margin = unit(c(0, 0, 0.2, 0), "lines"), plot.background = element_blank()) + theme( plot.title = element_text(family = "Times", size = 12), legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank() ) + annotate( "text", x = xcenter, y = ymax, label = "fitfrail", hjust = 0.5, vjust = 1 ) ymax <- 0.2 p2 <- ggplot(package.timings, aes(x = N, y = coxph.time, color = frailty)) + colscale + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylab("Runtime (s)") + ylim(0, ymax) + theme(plot.margin = unit(c(0.2, 0, 0.2, 0), "lines"), plot.background = element_blank()) + theme( legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank() ) + annotate( "text", x = xcenter, y = ymax, label = "coxph", hjust = 0.5, vjust = 1 ) ymax <- 20 p3 <- ggplot(package.timings, aes(x = N, y = frailtyPenal.time, color = frailty)) + colscale + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + xlab("N") + ylab("Runtime (s)") + ylim(0, ymax) + theme( legend.position = "bottom", legend.justification = c(0.5, 1), legend.direction = "horizontal", legend.title = element_blank(), plot.margin = unit(c(0.2, 0, 0, 0), "lines"), plot.background = element_blank() ) + annotate( "text", x = xcenter, y = ymax, label = "frailtyPenal", hjust = 0.5, vjust = 1 ) gt1 <- ggplotGrob(p1) gt2 <- ggplotGrob(p2) gt3 <- ggplotGrob(p3) pdf("figure17.pdf", width = 7, height = 6) print(grid.arrange(gtable_rbind(gt1, gt2, gt3))) dev.off() # Section C.3 set.seed(SEED) fitfrail.accuracy.abstol <- load_or_run( "fitfrail.accuracy.abstol", compare.fitfrail.accuracy, list( reps = REPS, frailty = FRAILTY.ALL, fitfrail.param = "abstol", fitfrail.param.values = TOLS ) ) set.seed(SEED) fitfrail.accuracy.reltol <- load_or_run( "fitfrail.accuracy.reltol", compare.fitfrail.accuracy, list( reps = REPS, frailty = FRAILTY.ALL, fitfrail.param = "reltol", fitfrail.param.values = TOLS ) ) set.seed(SEED) fitfrail.accuracy.int.abstol <- load_or_run( "fitfrail.accuracy.int.abstol", compare.fitfrail.accuracy, list( reps = REPS, frailty = FRAILTY.INT, fitfrail.param = "int.abstol", fitfrail.param.values = TOLS ) ) set.seed(SEED) fitfrail.accuracy.int.reltol <- load_or_run( "fitfrail.accuracy.int.reltol", compare.fitfrail.accuracy, list( reps = REPS, frailty = FRAILTY.INT, fitfrail.param = "int.reltol", fitfrail.param.values = TOLS ) ) # Figure 18 pdf(NULL) beta.ylim.1 <- 0.2 theta.ylim.1 <- 1.1 beta.ylim.2 <- 0.1 theta.ylim.2 <- 0.25 indx <- sapply(fitfrail.accuracy.abstol, is.factor) indx[["frailty"]] <- FALSE fitfrail.accuracy.abstol[indx] <- lapply(fitfrail.accuracy.abstol[indx], function(x) as.numeric(as.character(x))) indx <- sapply(fitfrail.accuracy.reltol, is.factor) indx[["frailty"]] <- FALSE fitfrail.accuracy.reltol[indx] <- lapply(fitfrail.accuracy.reltol[indx], function(x) as.numeric(as.character(x))) colscale <- scale_colour_manual(values = palette, breaks = c("gamma", "lognormal", "invgauss", "pvf"), labels = c("Gamma", "LN", "IG", "PVF")) BREAKS <- unique(fitfrail.accuracy.abstol$tol) LABELS <- parse(text = paste("10^", log10(BREAKS), sep = "")) p1 <- ggplot( fitfrail.accuracy.abstol, aes( x = fitfrail.accuracy.abstol$tol, y = fitfrail.accuracy.abstol$fitfrail.time, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylab("Runtime (s)") + ylim(0, 150) + theme( plot.margin = unit(c(0, 0, 0.2, 0), "lines"), plot.background = element_blank(), legend.position = "top", legend.justification = c(0.5, 1), legend.direction = "horizontal", legend.title = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank() ) p2 <- ggplot( fitfrail.accuracy.abstol, aes( x = fitfrail.accuracy.abstol$tol, y = fitfrail.accuracy.abstol$fitfrail.beta.1.res, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylab(TeX(sprintf("$ \\widehat{ \\beta } - \\beta $"))) + ylim(-beta.ylim.1, beta.ylim.1) + theme( plot.margin = unit(c(0.2, 0, 0.2, 0), "lines"), plot.background = element_blank(), legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank() ) p3 <- ggplot( fitfrail.accuracy.reltol, aes( x = fitfrail.accuracy.abstol$tol, y = fitfrail.accuracy.abstol$fitfrail.theta.res, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS, labels = LABELS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + xlab("abstol") + ylab(TeX(sprintf("$ \\widehat{ \\theta } - \\theta $"))) + ylim(-theta.ylim.1, theta.ylim.1) + theme( plot.margin = unit(c(0.2, 0, 0, 0), "lines"), plot.background = element_blank(), legend.position = "none", axis.title.x = element_text(family = "mono") ) p4 <- ggplot( fitfrail.accuracy.reltol, aes( x = fitfrail.accuracy.reltol$tol, y = fitfrail.accuracy.reltol$fitfrail.time, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylim(0, 150) + theme( plot.margin = unit(c(0, 0, 0.2, 0), "lines"), plot.background = element_blank(), legend.position = "top", legend.justification = c(0.5, 1), legend.direction = "horizontal", legend.title = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank() ) p5 <- ggplot( fitfrail.accuracy.reltol, aes( x = fitfrail.accuracy.reltol$tol, y = fitfrail.accuracy.reltol$fitfrail.beta.1.res, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylim(-beta.ylim.1, beta.ylim.1) + theme( plot.margin = unit(c(0.2, 0, 0.2, 0), "lines"), plot.background = element_blank(), legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank() ) p6 <- ggplot( fitfrail.accuracy.reltol, aes( x = fitfrail.accuracy.reltol$tol, y = fitfrail.accuracy.reltol$fitfrail.theta.res, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS, labels = LABELS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + xlab("reltol") + ylim(-theta.ylim.1, theta.ylim.1) + theme( plot.margin = unit(c(0.2, 0, 0, 0), "lines"), plot.background = element_blank(), legend.position = "none", axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank(), axis.title.x = element_text(family = "mono") ) mylegend <- gtable_filter(ggplot_gtable(ggplot_build(p1)), "guide-box") gt1 <- ggplotGrob(p1 + theme(legend.position = "none")) gt2 <- ggplotGrob(p2) gt3 <- ggplotGrob(p3) gt4 <- ggplotGrob(p4 + theme(legend.position = "none")) gt5 <- ggplotGrob(p5) gt6 <- ggplotGrob(p6) pdf("figure18.pdf", width = 7, height = 6) print(grid.arrange( gtable_cbind(gtable_rbind(gt1, gt2, gt3), gtable_rbind(gt4, gt5, gt6)), mylegend, heights = unit.c(unit(1, "npc") - mylegend$height, mylegend$height) )) dev.off() # Figure 19 pdf(NULL) indx <- sapply(fitfrail.accuracy.int.abstol, is.factor) indx[["frailty"]] <- FALSE fitfrail.accuracy.int.abstol[indx] <- lapply(fitfrail.accuracy.int.abstol[indx], function(x) as.numeric(as.character(x))) indx <- sapply(fitfrail.accuracy.int.reltol, is.factor) indx[["frailty"]] <- FALSE fitfrail.accuracy.int.reltol[indx] <- lapply(fitfrail.accuracy.int.reltol[indx], function(x) as.numeric(as.character(x))) FRAILTY.LABELS <- c("LN", "IG") colscale <- scale_colour_manual(values = palette, breaks = c("gamma", "lognormal", "invgauss", "pvf"), labels = c("Gamma", "LN", "IG", "PVF")) BREAKS <- unique(fitfrail.accuracy.int.abstol$tol) LABELS <- parse(text = paste("10^", log10(BREAKS), sep = "")) p1 <- ggplot( fitfrail.accuracy.int.abstol, aes( x = fitfrail.accuracy.int.abstol$tol, y = fitfrail.accuracy.int.abstol$fitfrail.time, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + xlab("int.abstol") + ylab("Runtime (s)") + ylim(0, 700) + theme( plot.margin = unit(c(0, 0, 0.2, 0), "lines"), plot.background = element_blank(), legend.position = "top", legend.justification = c(0.5, 1), legend.direction = "horizontal", legend.title = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank() ) p2 <- ggplot( fitfrail.accuracy.int.abstol, aes( x = fitfrail.accuracy.int.abstol$tol, y = fitfrail.accuracy.int.abstol$fitfrail.beta.1.res, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylab(TeX(sprintf("$ \\widehat{ \\beta } - \\beta $"))) + ylim(-beta.ylim.2, beta.ylim.2) + theme( plot.margin = unit(c(0.2, 0, 0.2, 0), "lines"), plot.background = element_blank(), legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank() ) p3 <- ggplot( fitfrail.accuracy.int.reltol, aes( x = fitfrail.accuracy.int.abstol$tol, y = fitfrail.accuracy.int.abstol$fitfrail.theta.res, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS, labels = LABELS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + xlab("int.abstol") + ylab(TeX(sprintf("$ \\widehat{ \\theta } - \\theta $"))) + ylim(-theta.ylim.2, theta.ylim.2) + theme( plot.margin = unit(c(0.2, 0, 0, 0), "lines"), plot.background = element_blank(), legend.position = "none", axis.title.x = element_text(family = "mono") ) p4 <- ggplot( fitfrail.accuracy.int.reltol, aes( x = fitfrail.accuracy.int.reltol$tol, y = fitfrail.accuracy.int.reltol$fitfrail.time, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylim(0, 700) + theme( plot.margin = unit(c(0, 0, 0.2, 0), "lines"), plot.background = element_blank(), legend.position = "top", legend.justification = c(0.5, 1), legend.direction = "horizontal", legend.title = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank() ) p5 <- ggplot( fitfrail.accuracy.int.reltol, aes( x = fitfrail.accuracy.int.reltol$tol, y = fitfrail.accuracy.int.reltol$fitfrail.beta.1.res, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + ylim(-beta.ylim.2, beta.ylim.2) + theme( plot.margin = unit(c(0.2, 0, 0.2, 0), "lines"), plot.background = element_blank(), legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank() ) p6 <- ggplot( fitfrail.accuracy.int.reltol, aes( x = fitfrail.accuracy.int.reltol$tol, y = fitfrail.accuracy.int.reltol$fitfrail.theta.res, color = frailty ) ) + colscale + scale_x_log10(breaks = BREAKS, labels = LABELS) + stat_summary(fun.data = mean_cl_boot, geom = "smooth") + xlab("int.reltol") + ylim(-theta.ylim.2, theta.ylim.2) + theme( plot.margin = unit(c(0.2, 0, 0, 0), "lines"), plot.background = element_blank(), legend.position = "none", axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank(), axis.title.x = element_text(family = "mono") ) mylegend <- gtable_filter(ggplot_gtable(ggplot_build(p1)), "guide-box") gt1 <- ggplotGrob(p1 + theme(legend.position = "none")) gt2 <- ggplotGrob(p2) gt3 <- ggplotGrob(p3) gt4 <- ggplotGrob(p4 + theme(legend.position = "none")) gt5 <- ggplotGrob(p5) gt6 <- ggplotGrob(p6) pdf("figure19.pdf", width = 7, height = 6) print(grid.arrange( gtable_cbind(gtable_rbind(gt1, gt2, gt3), gtable_rbind(gt4, gt5, gt6)), mylegend, heights = unit.c(unit(1, "npc") - mylegend$height, mylegend$height) )) dev.off()