R/plot.R
In stc04003/reReg: Recurrent Event Regression
Defines functions
basebind
plotEvents.control
plotHaz
plotRate
plot.reReg
plotMCF
plotEvents
plot.Recur
Documented in
basebind
plotEvents
plotEvents.control
plotHaz
plotRate
plot.Recur
plot.reReg
globalVariables(c("time1", "time2", "group", "Y", "Y.upper", "Y.lower", "id", "event", "origin", "MCF", "bs"))
#' Produce Event Plot or Mean Cumulative Function Plot
#'
#' Plot the event plot or the mean cumulative function (MCF) from an \code{Recur} object.
#'
#' The argument \code{control} consists of options with argument defaults to a list with
#' the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is "Subject" for event plot and
#' "Cumulative mean" for MCF plot.}
#' \item{main}{customizable title, the default value is "Recurrent event plot"
#' when \code{mcf = FALSE} and
#' "Sample cumulative mean function plot" when \code{mcf = TRUE}.}
#' \item{terminal.name}{customizable label for terminal event,
#' the default value is "Terminal event".}
#' \item{recurrent.name}{customizable legend title for recurrent event,
#' the default value is "Recurrent events".}
#' \item{recurrent.types}{customizable label for recurrent event type,
#' the default value is \code{NULL}.}
#' \item{alpha}{between 0 and 1, controls the transparency of points.}
#' }
#' The \code{xlab}, \code{ylab} and \code{main} parameters can be specified
#' outside of the \code{control} list.
#'
#' @param x an object of class \code{Recur} returned by the \code{Recur()} function.
#' See \code{?Recur} for creating \code{Recur} objects.
#' @param event.result an optional character string that is passed to the
#' \code{plotEvents()} function as the \code{result} argument. See \code{\link{plotEvents}}.
#' This argument is used to specify whether the event plot is sorted by the subjects' terminal time.
#' The available options are
#' \describe{
#' \item{\code{increasing}}{sort the terminal time from in ascending order (default).
#' This places longer terminal times on top. }
#' \item{\code{decreasing}}{sort the terminal time from in descending order.
#' This places shorter terminal times on top. }
#' \item{\code{none}}{present the event plots as is, without sorting by the terminal times.}
#' }
#' @param event.calendarTime an optional logical value indicating whether to plot in calendar time.
#' When \code{event.calendarTime = FALSE} (default),
#' the event plot will have patient time on the x-axis.
#' @param control a list of control parameters. See \bold{Details}.
#' @param mcf.adjustRiskset an optional logical value that is passed to
#' the \code{mcf()} function as the \code{adjustRiskset} argument.
#' This argument indicates whether risk set size will be adjusted.
#' If \code{mcf.adjustRiskset = TRUE}, subjects leave the risk set after terminal times
#' as in the Nelson-Aalen estimator.
#' If \code{mcf.adjustRiskset = FALSE}, subjects remain in the risk set after terminal time.
#' @param mcf.conf.int an optional logical value that is passed to
#' the \code{mcf()} function as the \code{conf.int} argument. See \code{\link{mcf}} for details.
#' @param mcf an optional logical value indicating whether the mean cumulative function (MCF) will
#' be plotted instead of the event plot. When \code{mcf = TRUE},
#' the \code{mcf} is internally called. See \code{\link{mcf}} for details.
#' @param ... additional graphical parameters to be passed to methods.
#'
#' @seealso \code{\link{Recur}}, \code{\link{plotEvents}}, \code{\link{mcf}}
#'
#' @references Nelson, W. B. (1995) Confidence Limits for Recurrence Data-Applied to Cost
#' or Number of Product Repairs. \emph{Technometrics}, \bold{37}(2): 147--157.
#'
#' @keywords Plots
#' @export
#'
#' @return A \code{ggplot} object.
#' @example inst/examples/ex_plot_Recur.R
#' @method plot Recur
plot.Recur <- function(x, mcf = FALSE,
event.result = c("increasing", "decreasing", "asis"),
event.calendarTime = FALSE,
mcf.adjustRiskset = TRUE,
mcf.conf.int = FALSE,
control = list(), ...) {
event.result <- match.arg(event.result)
if (!is.Recur(x)) stop("Response must be a `Recur` object.")
ctrl <- plotEvents.control()
if (mcf & is.null(ctrl$ylab)) ctrl <- ctrl$ylab <- "Cumulative mean"
call <- match.call()
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) eval(call[[x]]))
}
if (!mcf) {
return(plotEvents(x, result = event.result, calendarTime = event.calendarTime, control = ctrl))
}
if (mcf) {
return(plot(mcf(x ~ 1, adjustRiskset = mcf.adjustRiskset), conf.int = mcf.conf.int))
## return(plotMCF(x, adjustRiskset = mcf.adjustRiskset, smooth = mcf.smooth, control = ctrl))
}
}
#' Produce Event Plots
#'
#' Plot the event plot for an \code{Recur} object.
#' The usage of the function is similar to that of \code{plot.Recur()} but with more flexible options.
#'
#' The argument \code{control} consists of options with argument defaults to a list with
#' the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is "Subject" for event plot and
#' "Cumulative mean" for MCF plot.}
#' \item{main}{customizable title, the default value is "Recurrent event plot"
#' when \code{mcf = FALSE} and
#' "Sample cumulative mean function plot" when \code{mcf = TRUE}.}
#' \item{terminal.name}{customizable label for terminal event,
#' the default value is "Terminal event".}
#' \item{recurrent.name}{customizable legend title for recurrent event,
#' the default value is "Recurrent events".}
#' \item{recurrent.types}{customizable label for recurrent event type,
#' the default value is \code{NULL}.}
#' \item{alpha}{between 0 and 1, controls the transparency of points.}
#' }
#' The \code{xlab}, \code{ylab} and \code{main} parameters can be specified
#' outside of the \code{control} list.
#'
#' @param formula a formula object, with the response on the left of a "~" operator,
#' and the predictors on the right.
#' The response must be a recurrent event survival object as returned by function \code{Recur()}.
#' @param data an optional data frame in which to interpret the variables occurring in
#' the "\code{formula}".
#' @param result an optional character string specifying whether the event plot is
#' sorted by the subjects' terminal time. The available options are
#' \describe{
#' \item{\code{increasing}}{sort the terminal time from in ascending order (default).
#' This places longer terminal times on top. }
#' \item{\code{decreasing}}{sort the terminal time from in descending order.
#' This places shorter terminal times on top. }
#' \item{\code{none}}{present the event plots as is, without sorting by the terminal times.}
#' }
#' @param calendarTime an optional logical value indicating whether to plot in calendar time.
#' When \code{calendarTime = FALSE} (default), the event plot will have patient time on the x-axis.
#' @param control a list of control parameters. See \bold{Details}.
#' @param ... graphical parameters to be passed to methods.
#' These include \code{xlab}, \code{ylab}, \code{main}, and more. See \bold{Details}.
#'
#' @seealso \code{\link{Recur}}, \code{\link{plot.Recur}}
#'
#' @keywords Plots
#' @export
#'
#' @return A \code{ggplot} object.
#' @importFrom ggplot2 geom_rect ggplot_build scale_y_continuous unit
#' @example inst/examples/ex_plot_event.R
plotEvents <- function(formula, data, result = c("increasing", "decreasing", "asis"),
calendarTime = FALSE, control = list(), ...) {
result <- match.arg(result)
ctrl <- plotEvents.control()
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
call <- match.call()
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) eval(call[[x]]))
}
nX <- 0
if (is.Recur(formula)) {
DF <- as.data.frame(formula@.Data)
isDate <- "Date" %in% formula@time_class
vNames <- NULL
} else {
if (missing(data)) obj <- eval(formula[[2]], parent.frame())
else obj <- eval(formula[[2]], data)
if (!is.Recur(obj)) stop("Response must be a `Recur` object.")
nX <- length(formula[[3]])
isDate <- "Date" %in% obj@time_class
if (formula[[3]] == 1) DF <- as.data.frame(obj@.Data)
if (formula[[3]] != 1 && nX == 1) {
if (missing(data)) DF <- data.frame(obj@.Data, eval(formula[[3]], parent.frame()))
if (!missing(data)) DF <- data.frame(obj@.Data, eval(formula[[3]], data))
colnames(DF) <- c(colnames(obj@.Data), paste0(formula[[3]], collapse = ""))
DF <- as.data.frame(DF)
}
if (formula[[3]] != 1 && nX > 1) {
DF <- as.data.frame(obj@.Data)
if (missing(data)) {
for (i in 2:nX) {
DF <- cbind(DF, eval(formula[[3]][[i]], parent.frame()))
}
} else {
for (i in 2:nX) {
DF <- cbind(DF, eval(formula[[3]][[i]], data))
}
}
vNames <- attr(terms(formula), "term.labels")
colnames(DF) <- c(colnames(obj@.Data), vNames)
}
vNames <- attr(terms(formula), "term.labels")
if (length(vNames) == 0) vNames <- NULL
}
## dat$status <- ifelse(is.na(dat$status), 0, dat$status)
## dat$Yi <- ifelse(is.na(dat$Yi), unlist(lapply(dat$tij, max)), dat$Yi)
newIDtime2 <- function(dat, result = "increasing") {
if (result == "asis") {
tmp <- table(dat$id)
dat$id <- rep(1:length(tmp), tmp[match(unique(dat$id), names(tmp))])
return(dat)
} else {
dat <- dat[order(dat$id),]
## if (all(dat$origin == 0))
if (!calendarTime)
tmp <- rank((dat$time2 - dat$origin)[dat$event == 0], ties.method = "first")
else
tmp <- rank(dat$time2[dat$event == 0], ties.method = "first")
}
if (result == "decreasing") tmp <- length(tmp) - tmp + 1
dat$id <- rep(tmp, table(dat$id))
return(dat)
}
if (nX == 0 || formula[[3]] == 1) {
DF <- newIDtime2(DF, result = result)
} else {
DF <- do.call(rbind,
lapply(split(DF, DF[, 7:ncol(DF)], drop = TRUE), newIDtime2, result = result))
rownames(DF) <- NULL
}
if (is.null(ctrl$cex)) sz <- 1 + 8 / (1 + exp(length(unique(DF$id)) / 30)) / max(1, nX)
else sz <- ctrl$cex
k <- length(unique(DF$event)) - 1 ## exclude event = 0
shp.val <- c(ctrl$terminal.shape, rep(ctrl$recurrent.shape, k))
if (is.null(ctrl$recurrent.color))
ctrl$recurrent.color <- hcl(h = seq(120, 360, length.out = k), l = 60, alpha = ctrl$alpha)
if (length(ctrl$recurrent.color) < k)
ctrl$recurrent.color <- alpha(rep(ctrl$recurrent.color, length.out = k), ctrl$alpha)
clr.val <- c(alpha(ctrl$terminal.color, ctrl$alpha), ctrl$recurrent.color)
if (k == 0) rec.lab <- NULL
else rec.lab <- paste("r", 1:k, sep = "")
if (k == 0) shp.lab <- ctrl$terminal.name
if (k == 1) shp.lab <- c(ctrl$terminal.name, ctrl$recurrent.name)
if (k > 1 & is.null(ctrl$recurrent.type)) {
shp.lab <- c(ctrl$terminal.name, paste(ctrl$recurrent.name, 1:k))
}
if (k > 1 & !is.null(ctrl$recurrent.type)) {
if (length(ctrl$recurrent.type) == k) {
shp.lab <- c(ctrl$terminal.name, ctrl$recurrent.type)
} else {
message('The length of "recurrent.type" mismatched, default names are used.')
shp.lab <- c(ctrl$terminal.name, paste(ctrl$recurrent.name, 1:k))
}
}
if (nX > 0) {
for (i in vNames) {
DF[,i] <- factor(DF[,i], labels = paste(i, "=", unique(DF[,i])))
}}
names(shp.val) <- names(clr.val) <- c("terminal", rec.lab)
## Bars
if (calendarTime) {
if (is.null(ctrl$width)) ctrl$width <- 0
gg <- ggplot(DF, aes(xmin = id - .45 - ctrl$width,
xmax = id + .45 + ctrl$width,
ymin = time1, ymax = time2)) +
geom_rect(fill = ctrl$bar.color) +
coord_flip()
} else gg <- ggplot(DF[DF$event == 0,], aes(id, time2 - origin)) +
geom_bar(stat = "identity", fill = ctrl$bar.color, width = ctrl$width) +
coord_flip()
## event dots
if (!calendarTime) DF$time2 <- DF$time2 - DF$origin
if (any(table(DF$id) > 0))
gg <- gg + geom_point(data = DF[DF$event > 0,],
aes(id, time2,
shape = factor(event, labels = rec.lab),
color = factor(event, labels = rec.lab),
stroke = ctrl$recurrent.stroke),
size = sz)
if (sum(DF$terminal, na.rm = TRUE) > 0)
gg <- gg + geom_point(data = DF[DF$terminal > 0,],
aes(id, time2, shape = "terminal", color = "terminal",
stroke = ctrl$terminal.stroke),
size = sz)
if (!is.null(ctrl$not.terminal.shape))
gg <- gg + geom_point(data = DF[DF$terminal == 0,], aes(id, time2),
shape = ctrl$not.terminal.shape,
color = ctrl$not.terminal.color,
stroke = ctrl$terminal.stroke,
size = sz)
if (nX > 0 && formula[[3]] != 1)
gg <- gg + facet_grid(as.formula(paste(formula[3], "~.", collapse = "")),
scales = "free", space = "free", switch = "both")
## Add theme and final touch ups
if (ctrl$main != "") gg <- gg + ggtitle(ctrl$main)
gg <- gg + scale_shape_manual(name = "", values = shp.val,
labels = shp.lab, breaks = c("terminal", rec.lab)) +
scale_color_manual(name = "", values = clr.val,
labels = shp.lab, breaks = c("terminal", rec.lab)) +
theme(panel.background = element_blank(),
axis.line = element_line(color = "black"),
axis.line.y = element_blank(),
axis.title.y = element_text(vjust = 0),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.title = element_text(size = 2 * ctrl$base_size),
strip.text = element_text(size = ctrl$base_size),
legend.text = element_text(size = 1.5 * ctrl$base_size),
legend.title = element_text(size = 1.5 * ctrl$base_size),
axis.text = element_text(size = ctrl$base_size),
axis.title = element_text(size = 1.5 * ctrl$base_size)) +
scale_x_continuous(expand = c(0, 1)) +
labs(x = ctrl$ylab, y = ctrl$xlab) +
guides(shape = guide_legend(override.aes = list(size = 2.7)))
## Handle legends
if (ctrl$legend.position %in% c("bottomright", "bottomleft","topleft", "topright")) {
## legend inside
if (ctrl$legend.position == "bottomright") legendAdj <- c(1, 0)
if (ctrl$legend.position == "bottomleft") legendAdj <- c(0, 0)
if (ctrl$legend.position == "topright") legendAdj <- c(1, 1)
if (ctrl$legend.position == "topleft") legendAdj <- c(0, 1)
gg <- gg + theme(
legend.position = legendAdj, legend.justification = legendAdj,
legend.key = element_rect(fill = "transparent", colour = "transparent"),
legend.background = element_rect(fill = alpha(1e-5)))
} else {
## legend outside
gg <- gg + theme(legend.position = ctrl$legend.position,
legend.key = element_rect(fill = "white", color = "white"))
}
if (isDate & calendarTime) {
xl <- ggplot_build(gg)$layout$panel_params[[1]]$x$breaks
xl <- xl[!is.na(xl)]
gg <- gg + scale_y_continuous(breaks = xl, labels = as.Date(xl, origin = "1970-01-01"))
}
gg
}
#' Produce Cumulative Sample Mean Function Plots
#'
#' Plot the mean cumulative function (MCF) for an \code{Recur} object.
#' The usage of the function is similar to that of \code{plot.Recur()}
#' but with more flexible options.
#'
#' When \code{adjustRiskset = TRUE}, the \code{plotMCF()} is equivalent to
#' the Nelson-Aalen estimator for the intensity function of the recurrent event process.
#' When \code{adjustRiskset = FALSE}, the \code{plotMCF()} does not adjust for the risk set and
#' assumes all subjects remain at risk after the last observed recurrent event.
#' This is known as the survivor rate function.
#' The argument \code{control} consists of options with argument defaults
#' to a list with the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is "Cumulative mean".}
#' \item{main}{customizable title, default value is "Sample cumulative mean function plot".}
#' }
#' The \code{xlab}, \code{ylab} and \code{main} parameters can also be
#' specified outside of the \code{control} list.
#'
#' @param formula a formula object, with the response on the left of a "~" operator,
#' and the predictors on the right.
#' The response must be a recurrent event survival object returned by the \code{Recur()} function.
#' @param data an optional data frame in which to interpret the variables occurring in
#' the "\code{formula}".
#' @param adjustRiskset an optional logical value that is passed to the \code{plotMCF()} function
#' as the \code{adjustRiskset} argument. See \code{\link{plotMCF}}.
#' This argument indicates whether risk set size will be adjusted. If \code{mcf.adjustRiskset = TRUE},
#' subjects leave the risk set after terminal times as in the Nelson-Aalen estimator.
#' If \code{mcf.adjustRiskset = FALSE}, subjects remain in the risk set after terminal time.
#' @param smooth an optional logical value indicating whether to add a smooth curve
#' obtained from a monotone increasing P-splines implemented in package \code{scam}.
#' This feature only works for data with one recurrent event type.
#' @param onePanel an optional logical value indicating whether the mean cumulative functions
#' will be plotted in the same panel.
#' This is only useful when there are multiple recurrent event types or
#' in the presence of (discrete) covariates.
#' @param control a list of control parameters.
#' @param ... graphical parameters to be passed to methods.
#' These include \code{xlab}, \code{ylab}, \code{main}, and more. See \bold{Details}.
#'
#' @seealso \code{\link{Recur}}, \code{\link{plot.Recur}}
#' @keywords Plots
#'
#' @return A \code{ggplot} object.
#'
#' @importFrom ggplot2 guides guide_legend
#' @importFrom scam scam
#'
#' @noRd
#'
#' @example inst/examples/ex_plot_MCF.R
plotMCF <- function(formula, data, adjustRiskset = TRUE, onePanel = FALSE,
smooth = FALSE, control = list(), ...) {
call <- match.call()
ctrl <- plotEvents.control()
if (is.null(ctrl$ylab)) ctrl$ylab <- "Cumulative mean"
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) eval(call[[x]]))
}
nX <- 0
if(is.Recur(formula)) {
DF <- as.data.frame(formula@.Data)
vNames <- NULL
} else {
if (missing(data)) obj <- eval(formula[[2]], parent.frame())
else obj <- eval(formula[[2]], data)
if (!is.Recur(obj)) stop("Response must be a `Recur` object.")
DF <- as.data.frame(obj@.Data)
nX <- length(formula[[3]])
if (formula[[3]] != 1 && nX == 1) {
if (missing(data)) DF <- cbind(obj@.Data, tmp = eval(formula[[3]], parent.frame()))
if (!missing(data)) DF <- cbind(obj@.Data, tmp = eval(formula[[3]], data))
colnames(DF) <- c(colnames(obj@.Data), paste0(formula[[3]], collapse = ""))
DF <- as.data.frame(DF)
}
if (formula[[3]] != 1 && nX > 1) {
DF <- as.data.frame(obj@.Data)
if (missing(data)) {
for (i in 2:nX) {
DF <- cbind(DF, eval(formula[[3]][[i]], parent.frame()))
}
} else {
for (i in 2:nX) {
DF <- cbind(DF, eval(formula[[3]][[i]], data))
}
}
colnames(DF) <- c(colnames(obj@.Data), attr(terms(formula), "term.labels"))
}
vNames <- attr(terms(formula), "term.labels")
if (length(vNames) == 0) vNames <- NULL
}
## dd <- subset(DF, select = c("event", vNames, "time2"))
dd <- DF[,c("event", vNames, "time2")]
dd <- dd[do.call(order, dd),]
nn <- unlist(aggregate(time2~., data = dd, table)$time2)
dd <- unique(dd)
dd$n <- as.integer(nn)
rownames(dd) <- NULL
k <- length(unique(dd$event)) - 1
if (!is.null(vNames)) { ## any covariates/stratifications?
dd2 <- DF[DF$event == 0, vNames, drop = FALSE]
dd2 <- dd2[do.call(order, dd2),,drop = FALSE]
nn <- c(t(table(dd2)))
dd2 <- unique(dd2)
dd2$n <- as.integer(nn) ## tmp2 in the 1st version
rownames(dd2) <- NULL
dd$GrpInd <- match(apply(dd[,vNames, drop = FALSE], 1, paste, collapse = ""),
apply(dd2[,vNames, drop = FALSE], 1, paste, collapse = ""))
tmp1 <- merge(dd, dd2, by = vNames)
## as.integer(eval(parse(text = paste(attr(terms(formula), "term.labels"), collapse = ":"))))
rec0 <- tmp1[tmp1$event == 0,]
dat0 <- do.call(rbind, lapply(split(tmp1, tmp1$GrpInd), function(x) {
x$adjustRiskset = apply(x, 1, function(y)
as.numeric(y['n.y']) - sum(rec0$n.x[as.numeric(y['time2'])> rec0$time2 &
rec0$GrpInd == as.numeric(y['GrpInd'])]))
return(x)}))
dat0$n.x <- dat0$n.x * (dat0$event > 0)
rec0$time2 <- rec0$n.x <- 0
rec0$adjustRiskset <- 1
dat0 <- unique(rbind(dat0, rec0))
rownames(dat0) <- NULL
## Number of recurrent types
if (k > 1) {
tmp <- dat0[dat0$event == 0,]
tmp <- tmp[rep(1:NROW(tmp), k),]
tmp$event <- rep(1:k, each = sum(dat0$event == 0))
dat0 <- rbind(dat0[dat0$event > 0,], tmp)
} else {
dat0$event <- dat0$event[1]
}
dat0 <- dat0[order(dat0$event, dat0$GrpInd, dat0$time2),]
if (adjustRiskset) {
dat0 <- do.call(rbind,
lapply(split(dat0, list(dat0$event, dat0$GrpInd)), function(x) {
x$mu <- x$n.x / x$adjustRiskset
x$MCF <- cumsum(x$mu)
return(x)}))
} else {
dat0 <- do.call(rbind,
lapply(split(dat0, list(dat0$event, dat0$GrpInd)), function(x) {
x$mu <- x$n.x / x$n.y
x$MCF <- cumsum(x$mu)
return(x)}))
}
} else { ## no covariates
dd$n.y <- length(unique(DF$id))
rec0 <- dd[dd$event == 0, ]
dd$adjustRiskset <- apply(dd, 1, function(x) x[4] - sum(rec0$n[x[2] > rec0$time2]))
dd$n <- dd$n * (dd$event > 0)
rec0$time2 <- rec0$n <- 0
rec0$adjustRiskset <- 1
dat0 <- unique(rbind(dd, rec0))
dat0 <- dat0[order(dat0$time2),]
dat0$mu <- dat0$n / (adjustRiskset * dat0$adjustRiskset + (!adjustRiskset) * dat0$n.y)
dat0$MCF <- cumsum(dat0$mu)
}
if (k == 1) {
dat0$event <- dat0$event[1]
dat0$event <- factor(dat0$event, labels = ctrl$recurrent.name)
}
if (k > 1) {
dat00 <- dat0[dat0$event == 0,]
dat0 <- dat0[dat0$event > 0,]
if (nrow(dat00) > 0) {
dat00 <- dat00[rep(1:nrow(dat00), k),]
dat00$event <- rep(1:k, each = nrow(dat00) / k)
dat0 <- rbind(dat0, dat00)
}
if (is.null(ctrl$recurrent.type))
dat0$event <- factor(dat0$event, labels = paste(ctrl$recurrent.name, 1:k))
if (!is.null(ctrl$recurrent.type)) {
if (length(ctrl$recurrent.type) == k) {
dat0$event <- factor(dat0$event, labels = ctrl$recurrent.type)
} else {
message('The length of "recurrent.type" mismatched, default names are used.')
dat0$event <- factor(dat0$event, labels = paste(ctrl$recurrent.name, 1:k))
}
}
}
if (nX > 0) {
for (i in vNames) {
dat0[,i] <- factor(dat0[,i], labels = paste(i, "=", unique(dat0[,i])))
}}
dat0 <- dat0[complete.cases(dat0),]
gg <- ggplot(data = dat0, aes(x = time2, y = MCF))
if (is.null(vNames) & k == 1) {
gg <- gg + geom_step(size = ctrl$lwd)
} else {
if (!onePanel & k == 1) gg <- gg + geom_step(size = ctrl$lwd)
if (!onePanel & k > 1)
gg <- gg + geom_step(aes(color = event), direction = "hv", size = ctrl$lwd) +
guides(color = guide_legend(title = ctrl$recurrent.name))
if (onePanel) {
dat0$GrpInd <- factor(dat0$GrpInd)
levels(dat0$GrpInd) <- apply(unique(dat0[,vNames, drop = FALSE]), 1, paste, collapse = ", ")
gg <- gg + geom_step(aes(color = dat0$GrpInd), direction = "hv", size = ctrl$lwd) +
guides(color = guide_legend(title = ""))
}
}
if (!onePanel && nX > 0 && formula[[3]] != 1)
gg <- gg + facet_grid(as.formula(paste(formula[3], "~.", collapse = "")),
scales = "free", space = "free", switch = "x")
if (!onePanel & k == 1)
gg <- gg + theme(legend.position = "none")
## if (onePanel & k == 1)
## gg <- gg + scale_color_discrete(name = "", labels = levels(interaction(vNames)))
## if (onePanel & k > 1) gg <- gg + scale_color_discrete(name = "")
if (smooth & k == 1 & !onePanel) {
if (is.null(dat0$GrpInd)) dat0$GrpInd <- 1
dat0 <- do.call(rbind, lapply(split(dat0, dat0$GrpInd), function(x){
x$bs <- scam(x$MCF ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values
return(x)}))
gg <- gg + geom_line(data = dat0, aes(time2, y = bs), color = 4, size = ctrl$lwd)
## geom_smooth(method = "scam", formula = y ~ s(x, k = 10, bs = "mpi"), size = ctrl$lwd, se = FALSE)
}
## gg <- gg + geom_smooth(method = "loess", size = ctrl$lwd, se = FALSE)
if (smooth & k > 1) message('Smoothing only works for data with one recurrent event type.')
if (ctrl$main != "") gg <- gg + ggtitle(ctrl$main)
gg <- gg + theme(axis.line = element_line(color = "black"),
plot.title = element_text(size = 2 * ctrl$base_size),
strip.text = element_text(size = ctrl$base_size),
legend.text = element_text(size = 1.5 * ctrl$base_size),
legend.title = element_text(size = 1.5 * ctrl$base_size),
axis.text = element_text(size = ctrl$base_size),
axis.title = element_text(size = 1.5 * ctrl$base_size)) +
labs(y = ctrl$ylab, x = ctrl$xlab)
## Handle legends
if (ctrl$legend.position %in% c("bottomright", "bottomleft","topleft", "topright")) {
## legend inside
if (ctrl$legend.position == "bottomright") legendAdj <- c(1, 0)
if (ctrl$legend.position == "bottomleft") legendAdj <- c(0, 0)
if (ctrl$legend.position == "topright") legendAdj <- c(1, 1)
if (ctrl$legend.position == "topleft") legendAdj <- c(0, 1)
gg <- gg + theme(
legend.position = legendAdj, legend.justification = legendAdj,
legend.key = element_rect(fill = "transparent", colour = "transparent"),
legend.background = element_rect(fill = alpha(1e-5)))
} else {
## legend outside
gg <- gg + theme(legend.position = ctrl$legend.position,
legend.key = element_rect(fill = "white", color = "white"))
}
gg
}
#' Plot the Baseline Cumulative Rate Function and the Baseline Cumulative Hazard Function
#'
#' Plot the baseline cumulative rate function and the baseline cumulative hazard function
#' (if applicable) for an \code{reReg} object.
#'
#' The argument \code{control} consists of options with argument defaults to a list
#' with the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is empty.}
#' \item{main}{customizable title, default value are "Baseline cumulative rate and
#' hazard function" when \code{baseline = "both"},
#' "Baseline cumulative rate function" when \code{baseline = "rate"},
#' and "Baseline cumulative hazard function" when \code{baseline = "hazard"}.}
#' }
## #' These arguments can also be passed down without specifying a \code{control} list.
#'
#' @param x an object of class \code{reReg}, returned by the \code{reReg} function.
#' @param smooth an optional logical value indicating whether to add a smooth curve
#' obtained from a monotone increasing P-splines implemented in package \code{scam}.
#' @param baseline a character string specifying which baseline function to plot.
#' \describe{
#' \item{\code{baseline = "both"}}{plot both the baseline cumulative rate and
#' the baseline cumulative hazard function (if applicable) in separate panels
#' within the same display (default).}
#' \item{\code{baseline = "rate"}}{plot the baseline cumulative rate function.}
#' \item{\code{baseline = "hazard"}}{plot the baseline cumulative hazard function.}
#' }
## #' @param rateType a character string specifying the type of rate function to be plotted.
## #' This argument is passed to the \code{plotRate()} function as the \code{type} argument.
## #' See \code{\link{plotRate}}.
## #' Options are "unrestricted", "scaled", "bounded".
#' @param newdata an optional data frame contains variables to include in the calculation
#' of the cumulative rate function.
#' If omitted, the baseline rate function will be plotted.
#' @param frailty an optional vector to specify the shared frailty for \code{newdata}.
#' If \code{newdata} is given and \code{frailty} is not specified, the
#' @param showName an optional logical value indicating whether to label the curves
#' when \code{newdata} is specified.
#' @param control a list of control parameters. See \bold{Details}.
#' @param ... additional graphical parameters to be passed to methods.
#'
#' @seealso \code{\link{reReg}}
#' @export
#' @keywords Plots
#'
#' @return A \code{ggplot} object.
#'
#' @importFrom ggplot2 geom_smooth geom_step ggplotGrob
#' @importFrom grid grid.draw
#' @example inst/examples/ex_plot_reReg.R
#' @method plot reReg
plot.reReg <- function(x,
baseline = c("both", "rate", "hazard"),
## type = c("unrestricted", "bounded", "scaled"),
smooth = FALSE, newdata = NULL, frailty = NULL, showName = FALSE,
control = list(), ...) {
baseline <- match.arg(baseline)
## type <- match.arg(type)
type <- "unrestricted"
if (x$typeRec %in% c("cox.GL", "cox.LWYY", "am.GL"))
stop("Baseline functions not yet available for this method.")
if (x$typeRec %in% c("cox.LWYY", "cox.HH"))
ctrl <- plotEvents.control(ylab = "Rate")
if (baseline %in% c("both", "rate"))
ctrl <- plotEvents.control(ylab = "Rate")
if (baseline == "hazard")
ctrl <- plotEvents.control(ylab = "Hazard")
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
call <- match.call()
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) call[[x]])
}
if (!is.reReg(x)) stop("Response must be a `reReg` class")
if (baseline == "rate")
return(plotRate(x, smooth = smooth, type = type,
newdata = newdata, frailty = frailty, showName = showName, control = ctrl))
if (baseline == "hazard")
return(plotHaz(x, smooth = smooth,
newdata = newdata, frailty = frailty, showName = showName, control = ctrl))
if (x$typeTem == ".") {
## cat(paste("Baseline cumulative hazard function is not available."))
## cat("\nOnly the baseline cumulative rate function is plotted.\n")
return(plotRate(x, smooth = smooth, type = type,
newdata = newdata, frailty = frailty, showName = showName, control = ctrl))
}
g1 <- plotRate(x, smooth = smooth, type = type,
newdata = newdata, frailty = frailty, showName = showName, control = ctrl)
g2 <- plotHaz(x, smooth = smooth, type = type,
newdata = newdata, frailty = frailty, showName = showName, control = ctrl)
if (ctrl$main != "") g1 <- g1 + ggtitle(ctrl$main)
g1 <- g1 + ylab("Rate") + xlab("") +
## facet_grid(factor(rep(1, nrow(x$DF)), levels = 1, labels = "Baseline cumulative rate")) +
theme(axis.title.x = element_blank(),
axis.ticks.x = element_blank(),
plot.margin = unit(c(0, 0, 0, 0), "cm"))
g2 <- g2 + ylab("Hazard") + xlab("Time") +
## facet_grid(factor(rep(1, nrow(x$DF)), levels = 1, labels = "Baseline cumulative hazard")) +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
g1 <- ggplotGrob(g1)
g2 <- ggplotGrob(g2)
grid.draw(rbind(g1, g2))
## gg + ggtitle(ctrl$main)
}
#' Plotting the Baseline Cumulative Rate Function for the Recurrent Event Process
#'
#' Plot the baseline cumulative rate function for an \code{reReg} object.
#'
#'
#' The \code{plotRate()} plots the estimated baseline cumulative rate function
#' depending on the identifiability assumption.
#' When \code{type = "unrestricted"} (default), the baseline cumulative rate function
#' is plotted under the assumption \eqn{E(Z) = 1}.
#' When \code{type = "scaled"}, the baseline cumulative rate function is plotted
#' under the assumption \eqn{\Lambda(\min(Y^\ast, \tau)) = 1}.
#' When \code{type = "bounded"}, the baseline cumulative rate function is plotted
#' under the assumption \eqn{\Lambda(\tau) = 1}.
#' See \code{?reReg} for the specification of the notations and underlying models.
#'
#' The argument \code{control} consists of options with argument defaults
#' to a list with the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is empty.}
#' \item{main}{customizable title, default value is "Baseline cumulative rate function".}
#' }
#' These arguments can also be specified outside of the \code{control} list.
#'
#' @param x an object of class \code{reReg}, usually returned by the \code{reReg} function.
#' @param type a character string specifying the type of rate function to be plotted.
#' Options are "unrestricted", "scaled", "bounded". See \bold{Details}.
#' @param smooth an optional logical value indicating whether to add a smooth curve
#' obtained from a monotone increasing P-splines implemented in package \code{scam}.
#' @param newdata an optional data frame contains variables to include in the calculation
#' of the cumulative rate function.
#' If omitted, the baseline rate function will be plotted.
#' @param frailty an optional vector to specify the shared frailty for \code{newdata}.
#' If \code{newdata} is given and \code{frailty} is not specified, the
#' @param showName an optional logical value indicating whether to label the curves
#' when \code{newdata} is specified.
#' @param control a list of control parameters.
#' @param ... graphical parameters to be passed to methods.
#' These include \code{xlab}, \code{ylab}, \code{main}, and more. See \bold{Details}.
#'
#' @seealso \code{\link{reReg}} \code{\link{plot.reReg}}
#' @export
#'
#' @return A \code{ggplot} object.
#'
#' @keywords Plots
#' @importFrom ggplot2 scale_x_continuous
#' @importFrom directlabels geom_dl
#'
#' @example inst/examples/ex_plot_rate.R
plotRate <- function(x, newdata = NULL, frailty = NULL, showName = FALSE,
type = c("unrestricted", "bounded", "scaled"),
smooth = FALSE, control = list(), ...) {
if (x$typeRec %in% c("cox.GL", "cox.LWYY", "am.GL"))
stop("Baseline cumulative rate function is not available")
if (is.null(frailty)) frailty <- exp(x$log.muZ)
if (length(frailty) > 1 & !is.null(newdata) && length(frailty) != nrow(newdata))
stop("newdata and frailty are different lengths")
## ctrl <- plot.reReg.control(main = "Baseline cumulative rate function")
ctrl <- plotEvents.control(ylab = "Rate")
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
call <- match.call()
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) call[[x]])
}
type <- match.arg(type)
if (x$typeRec == "nonparametric") {
type <- "bounded"
ctrl$ylab = "MCF Estimates"
}
if (!is.reReg(x)) stop("Response must be a `reReg` class")
dat <- x$DF[,"time2",drop = FALSE]
if (is.null(newdata)) {
if (type == "unrestricted") dat$Y <- x$Lam0(dat$time2) * exp(x$log.muZ)
if (type == "scaled") dat$Y <- x$Lam0(dat$time2) / x$Lam0(max(dat$time2))
if (type == "bounded") dat$Y <- x$Lam0(dat$time2)
if (!is.null(x$Lam0.upper)) {
if (type == "bounded") {
dat$Y.upper <- x$Lam0.upper(dat$time2)
dat$Y.lower <- x$Lam0.lower(dat$time2)
}
if (type == "unrestricted") {
dat$Y.upper <- x$Lam0.upper(dat$time2) * exp(x$log.muZ)
dat$Y.lower <- x$Lam0.lower(dat$time2) * exp(x$log.muZ)
}
if (type == "scaled") {
dat$Y.upper <- x$Lam0.upper(dat$time2) / x$Lam0(max(dat$time2))
dat$Y.lower <- x$Lam0.lower(dat$time2) / x$Lam0(max(dat$time2))
}
}
gg <- ggplot(data = dat, aes(x = time2, y = Y)) +
theme(axis.line = element_line(color = "black"))
if (smooth) {
dat$bs <- scam(dat$Y ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
gg <- gg + geom_line(aes(x = time2, y = dat$bs), color = 4)
if (!is.null(x$Lam0.upper)) {
dat$bs.upper <- scam(dat$Y.upper ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
dat$bs.lower <- scam(dat$Y.lower ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
gg <- gg + geom_line(aes(time2, y = dat$bs.upper), color = 4, lty = 2) +
geom_line(aes(time2, y = dat$bs.lower), color = 4, lty = 2)
}
} else {
gg <- gg + geom_step()
if (!is.null(x$Lam0.upper))
gg <- gg + geom_step(aes(x = time2, y = Y.upper), lty = 2) +
geom_step(aes(x = time2, y = Y.lower), lty = 2)
}
}
if (!is.null(newdata)) {
if (!is.null(x$xlevels)) {
for(i in which(names(newdata) %in% names(x$xlevels)))
newdata[,i] <- factor(newdata[,i], levels = x$xlevels[[i]])
newdata <- model.matrix(~., newdata)
}
X <- as.matrix(unique(newdata[,match(x$varNames, colnames(newdata)), drop = FALSE]))
if (ncol(X) != length(x$varNames))
stop(paste0("Variables ",
paste(setdiff(x$varNames, colnames(newdata)), collapse = ", "),
" are missing"))
p <- ncol(X)
exa1 <- exa2 <- 1
if (x$typeRec == "cox") exa2 <- exp(X %*% x$par1)
if (x$typeRec == "ar") exa1 <- exp(X %*% x$par1)
if (x$typeRec == "am") exa1 <- exa2 <- exp(X %*% x$par1)
if (x$typeRec == "sc") {
exa1 <- exp(X %*% x$par1)
exa2 <- exp(X %*% x$par2)
}
exa1 <- rep(drop(exa1), each = nrow(dat))
exa2 <- rep(drop(exa2), each = nrow(dat))
Y <- frailty * x$Lam0(dat$time2 * exa1) * exa2 / exa1
if (!is.null(x$Lam0.upper)) {
Y.upper <- frailty * x$Lam0.upper(dat$time2 * exa1) * exa2 / exa1
Y.lower <- frailty * x$Lam0.lower(dat$time2 * exa1) * exa2 / exa1
dat <- data.frame(time2 = dat$time2,
id = rep(rownames(X), each = nrow(dat)),
Y = Y, Y.upper = Y.upper, Y.lower = Y.lower)
} else
dat <- data.frame(time2 = dat$time2, id = rep(rownames(X), each = nrow(dat)), Y = Y)
gg <- ggplot(data = dat, aes(x = time2, y = Y, group = id)) +
theme(axis.line = element_line(color = "black"))
if (smooth) {
dat$bs <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
gg <- gg + geom_line(aes(x = time2, y = dat$bs, group = id), color = 4)
if (!is.null(x$Lam0.upper)) {
dat$bs.upper <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y.upper ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
dat$bs.lower <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y.lower ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
gg <- gg +
geom_line(aes(x = time2, y = dat$bs.upper, group = id), color = 4, lty = 2) +
geom_line(aes(x = time2, y = dat$bs.lower, group = id), color = 4, lty = 2)
}
} else {
gg <- gg + geom_step()
if (!is.null(x$Lam0.upper))
gg <- gg +
geom_step(aes(x = time2, y = Y.upper, group = id), lty = 2) +
geom_step(aes(x = time2, y = Y.lower, group = id), lty = 2)
}
if (showName)
gg <- gg + geom_dl(aes(label = paste(" Obs. =", id)), method = "right.polygons") +
scale_x_continuous(limits = c(0, max(dat$time2) * 1.1))
}
if (ctrl$main != "") gg <- gg + ggtitle(ctrl$main)
gg <- gg + labs(x = ctrl$xlab, y = ctrl$ylab) +
theme(plot.title = element_text(size = 2 * ctrl$base_size),
strip.text = element_text(size = ctrl$base_size),
legend.text = element_text(size = 1.5 * ctrl$base_size),
legend.title = element_text(size = 1.5 * ctrl$base_size),
axis.line = element_blank(),
axis.text = element_text(size = ctrl$base_size),
axis.title = element_text(size = 1.5 * ctrl$base_size))
attr(gg, "from") <- "reReg"
return(gg)
}
#' Plot the Baseline Cumulative Hazard Function for the Terminal Time
#'
#' Plot the baseline cumulative hazard function for an \code{reReg} object.
#' The 95\% confidence interval on the baseline cumulative rate function
#'
#' The argument \code{control} consists of options with argument
#' defaults to a list with the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is empty.}
#' \item{main}{customizable title, default value is "Baseline cumulative hazard function".}
#' }
#' These arguments can also be passed down without specifying a \code{control} list.
#'
#' @param x an object of class \code{reReg}, returned by the \code{reReg} function.
#' @param smooth an optional logical value indicating whether to add a smooth curve
#' obtained from a monotone increasing P-splines implemented in package \code{scam}.
#' @param newdata an optional data frame contains variables to include in the calculation
#' of the cumulative rate function.
#' If omitted, the baseline rate function will be plotted.
#' @param frailty an optional vector to specify the shared frailty for \code{newdata}.
#' If \code{newdata} is given and \code{frailty} is not specified, the
#' @param showName an optional logical value indicating whether to label the curves
#' when \code{newdata} is specified.
#' @param type a character string specifying the type of rate function to be plotted.
#' Options are "unrestricted", "scaled", "bounded". See \bold{Details}.
#' @param control a list of control parameters.
#' @param ... graphical parameters to be passed to methods.
#' These include \code{xlab}, \code{ylab}, \code{main}, and more. See \bold{Details}.
#'
#' @seealso \code{\link{reReg}} \code{\link{plot.reReg}}
#' @export
#'
#' @return A \code{ggplot} object.
#' @keywords Plots
#'
#' @example inst/examples/ex_plot_Haz.R
plotHaz <- function(x, newdata = NULL, frailty = NULL, showName = FALSE,
type = c("unrestricted", "bounded", "scaled"),
smooth = FALSE, control = list(), ...) {
if (x$typeRec %in% c("cox.GL", "cox.LWYY", "am.GL"))
stop("Baseline cumulative hazard function is not available.")
if (x$typeTem == ".") {
stop("Baseline cumulative hazard function is not available.")
}
if (is.null(frailty)) frailty <- exp(x$log.muZ)
if (length(frailty) > 1 & !is.null(newdata) && length(frailty) != nrow(newdata))
stop("newdata and frailty are different lengths")
## ctrl <- plot.reReg.control(main = "Baseline cumulative hazard function")
ctrl <- plotEvents.control(ylab = "Hazard")
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
call <- match.call()
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) call[[x]])
}
type <- match.arg(type)
if (!is.reReg(x)) stop("Response must be a `reReg` class.")
dat <- x$DF[, "time2", drop = FALSE]
if (is.null(newdata)) {
if (type == "unrestricted") dat$Y <- x$Haz0(dat$time2) * exp(x$log.muZ)
if (type == "scaled") dat$Y <- x$Haz0(dat$time2) / x$Haz0(max(dat$time2))
if (type == "bounded") dat$Y <- x$Haz0(dat$time2)
## dat$Y <- x$Haz0(dat$time2)
if (!is.null(x$Haz0.upper)) {
if (type == "bounded") {
dat$Y.upper <- x$Haz0.upper(dat$time2)
dat$Y.lower <- x$Haz0.lower(dat$time2)
}
if (type == "unrestricted") {
dat$Y.upper <- x$Haz0.upper(dat$time2) * exp(x$log.muZ)
dat$Y.lower <- x$Haz0.lower(dat$time2) * exp(x$log.muZ)
}
if (type == "scaled") {
dat$Y.upper <- x$Haz0.upper(dat$time2) / x$Haz0(max(dat$time2))
dat$Y.lower <- x$Haz0.lower(dat$time2) / x$Haz0(max(dat$time2))
}
}
gg <- ggplot(data = dat, aes(x = time2, y = Y)) +
theme(axis.line = element_line(color = "black"))
if (smooth) {
dat$bs <- scam(dat$Y ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
gg <- gg + geom_line(aes(time2, y = dat$bs), color = 4)
if (!is.null(x$Lam0.upper)) {
dat$bs.upper <- scam(dat$Y.upper ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
dat$bs.lower <- scam(dat$Y.lower ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
gg <- gg + geom_line(aes(time2, y = dat$bs.upper), color = 4, lty = 2) +
geom_line(aes(time2, y = dat$bs.lower), color = 4, lty = 2)
}
} else {
gg <- gg + geom_step()
if (!is.null(x$Lam0.upper))
gg <- gg + geom_step(aes(x = time2, y = Y.upper), lty = 2) +
geom_step(aes(x = time2, y = Y.lower), lty = 2)
}
}
if (!is.null(newdata)) {
X <- as.matrix(unique(newdata[,match(x$varNames, names(newdata)), drop = FALSE]))
if (ncol(X) != length(x$varNames))
stop(paste0("Variables ",
paste(setdiff(x$varNames, names(newdata)), collapse = ", "),
" are missing"))
p <- ncol(X)
exb1 <- exb2 <- 1
if (x$typeTem == "cox") exb2 <- exp(X %*% x$par3)
if (x$typeTem == "ar") exb1 <- exp(X %*% x$par3)
if (x$typeTem == "am") exb1 <- exb2 <- exp(X %*% x$par3)
if (x$typeTem == "sc") {
exb1 <- exp(X %*% x$par3)
exb2 <- exp(X %*% x$par4)
}
exb1 <- rep(drop(exb1), each = nrow(dat))
exb2 <- rep(drop(exb2), each = nrow(dat))
Y <- frailty * x$Haz0(dat$time2 * exb1) * exb2 / exb1
if (!is.null(x$Haz0.upper)) {
Y.upper <- frailty * x$Haz0.upper(dat$time2 * exb1) * exb2 / exb1
Y.lower <- frailty * x$Haz0.lower(dat$time2 * exb1) * exb2 / exb1
dat <- data.frame(time2 = dat$time2, id = rep(rownames(X), each = nrow(dat)),
Y = Y, Y.upper = Y.upper, Y.lower = Y.lower)
} else
dat <- data.frame(time2 = dat$time2, id = rep(rownames(X), each = nrow(dat)), Y = Y)
gg <- ggplot(data = dat, aes(x = time2, y = Y, group = id)) +
theme(axis.line = element_line(color = "black"))
if (smooth) {
dat$bs <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
gg <- gg + geom_line(aes(x = time2, y = dat$bs, group = id), color = 4)
if (!is.null(x$Haz0.upper)) {
dat$bs.upper <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y.upper ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
dat$bs.lower <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y.lower ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
gg <- gg +
geom_line(aes(x = time2, y = dat$bs.upper, group = id), color = 4, lty = 2) +
geom_line(aes(x = time2, y = dat$bs.lower, group = id), color = 4, lty = 2)
}
} else {
gg <- gg + geom_step()
if (!is.null(x$Haz0.upper))
gg <- gg +
geom_step(aes(x = time2, y = Y.upper, group = id), lty = 2) +
geom_step(aes(x = time2, y = Y.lower, group = id), lty = 2)
}
if (showName)
gg <- gg + geom_dl(aes(label = paste(" Obs. =", id)), method = "right.polygons") +
scale_x_continuous(limits = c(0, max(dat$time2) * 1.1))
}
if (ctrl$main != "") gg <- gg + ggtitle(ctrl$main)
gg <- gg + labs(x = ctrl$xlab, y = ctrl$ylab) +
theme(plot.title = element_text(size = 2 * ctrl$base_size),
strip.text = element_text(size = ctrl$base_size),
legend.text = element_text(size = 1.5 * ctrl$base_size),
legend.title = element_text(size = 1.5 * ctrl$base_size),
axis.line = element_blank(),
axis.text = element_text(size = ctrl$base_size),
axis.title = element_text(size = 1.5 * ctrl$base_size))
attr(gg, "from") <- "reReg"
return(gg)
}
#' Plot options for plotEvents
#'
#' This function provides the plotting options for the \code{plotEvents()} function.
#'
#' @param xlab a character string indicating the label for the x axis.
#' The default value is "Time".
#' @param ylab a character string indicating the label for the y axis.
#' The default value is "Subject".
#' @param main a character string indicating the title of the plot.
#' @param terminal.name a character string indicating the label for the terminal event
#' displayed in the legend. The default value is "Terminal event".
#' @param recurrent.name a character string indicating the label for the recurrent event
#' displayed in the legend. The default value is "Recurrent events".
#' @param recurrent.type a factor indicating the labels for the different recurrent event types.
#' This option is only available when there are more than one types of recurrent events.
#' The default value is "Recurrent events 1", "Recurrent events 2", ....
#' @param legend.position a character string specifies the position of the legend.
#' The available options are "none", "left", "right", "bottom", "top",
#' "bottomright", "bottomleft","topleft", "topright",
#' or a two-element numeric vector specifies the coordinate of the legend.
#' The legend is placed inside of the plotting region when
#' \code{legend.position} is specified as one of
#' "bottomright", "bottomleft","topleft", "topright".
#' Otherwise, the argument is passed to the \code{ggplot} theme environment.
#' @param base_size a numerical value to specify the base font size, given in pts.
#' This argument is passed to the \code{ggplot} theme environment.
#' The default value is 12.
#' @param cex a numerical value specifies the size of the points.
#' @param alpha a numerical value specifies the transparency of the points.
#' @param width a numerical value specifies the width of the event plot.
#' By \code{ggplot} default, set to 90\% of the resolution of the data.
#' @param bar.color a numerical value or a character string specifies
#' color for lines. Default to gray.
#' @param recurrent.color a numerical value or a character string
#' specifies color for recurrent events. Default to green.
#' @param terminal.color a numerical value or a character string
#' specifies color for terminal events. Default to red.
#' @param recurrent.shape a numerical value or a character string
#' specifies shape for recurrent events. Default to circle.
#' @param terminal.shape a numerical value or a character string
#' specifies shape for terminal events. Default to triangle.
#' @param recurrent.stroke a numerical value or a character string
#' specifies stroke for recurrent events. Default to circle.
#' @param terminal.stroke a numerical value or a character string
#' specifies stroke for terminal events. Default to triangle.
#' @param not.terminal.color a numerical value or a character string
#' specifies color for non-terminal events.
#' Non-terminal events are not plotted at default.
#' @param not.terminal.shape a numerical value or a character string
#' specifies shape for terminal events.
#' Non-terminal events are not plotted at default.
#'
#' @seealso \code{\link{plotEvents}}
#' @export
plotEvents.control <- function(xlab = NULL, ylab = NULL,
main = NULL,
terminal.name = NULL,
recurrent.name = NULL,
recurrent.type = NULL,
legend.position = NULL,
base_size = 12,
cex = NULL,
alpha = .7,
width = NULL,
bar.color = NULL,
recurrent.color = NULL,
recurrent.shape = NULL,
recurrent.stroke = NULL,
terminal.color = NULL,
terminal.shape = NULL,
terminal.stroke = NULL,
not.terminal.color = NULL,
not.terminal.shape = NULL) {
if (is.null(ylab)) ylab <- "Subject"
if (is.null(xlab)) xlab <- "Time"
if (is.null(main)) main <- ""
if (is.null(legend.position)) legend.position <- "bottomright"
## main <- "Recurrent event plot"
if (is.null(terminal.name)) terminal.name <- "Terminal event"
if (is.null(recurrent.name)) recurrent.name <- "Recurrent events"
if (is.null(bar.color)) bar.color <- "gray75"
if (is.null(terminal.color)) terminal.color <- "red"
if (is.null(terminal.shape)) terminal.shape <- 17
if (is.null(recurrent.shape)) recurrent.shape <- 19
list(xlab = xlab, ylab = ylab, main = main, cex = cex, width = width,
terminal.name = terminal.name, recurrent.name = recurrent.name,
recurrent.type = recurrent.type, alpha = alpha,
legend.position = legend.position, base_size = base_size,
bar.color = bar.color,
recurrent.color = recurrent.color,
recurrent.shape = recurrent.shape,
recurrent.stroke = recurrent.stroke,
not.terminal.color = not.terminal.color,
not.terminal.shape = not.terminal.shape,
terminal.color = terminal.color,
terminal.shape = terminal.shape,
terminal.stroke = terminal.stroke)
}
#' Function used to combine baseline functions in one plot
#'
#' Combine different plots into one.
#'
#' @param ... \code{ggplot} objects created by plotting \code{reReg} objects.
#' @param legend.title an optional character string to specify the legend title.
#' @param legend.labels an optional character string to specify the legend labels.
#' @param control a list of control parameters.
#'
#' @export
#' @keywords Plots
#'
#' @example inst/examples/ex_basebind.R
basebind <- function(..., legend.title, legend.labels, control = list()) {
gglst <- list(...)
if (any(sapply(gglst, function(x) attr(x, "from")) != "reReg"))
stop("Plots must be created from reReg objects")
if (missing(legend.title)) legend.title <- ""
ctrl <- plotEvents.control()
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
if (is.null(ctrl$ylab)) ctrl$ylab <- ""
nargs <- length(gglst)
if (missing(legend.labels)) legend.labels <- 1:nargs
if (length(legend.labels) != nargs) {
message('The length of "name" mismatched, default names are used.')
legend.labels <- 1:nargs
}
d <- do.call(rbind, lapply(gglst, function(x) x$data))
nobs <- sapply(gglst, function(x) nrow(x$data))
d$group <- factor(rep(1:nargs, nobs), labels = legend.labels)
gg <- ggplot(data = d, aes(x = time2, y = Y, color = group)) + geom_step()
if (!is.null(d$Y.upper))
gg <- gg + geom_step(aes(x = time2, y = Y.upper), lty = 2) +
geom_step(aes(x = time2, y = Y.lower), lty = 2)
gg <- gg + labs(x = gglst[[1]]$labels$x, y = gglst[[1]]$labels$y, color = legend.title) +
theme(plot.title = element_text(size = 2 * ctrl$base_size),
strip.text = element_text(size = ctrl$base_size),
legend.text = element_text(size = 1.5 * ctrl$base_size),
legend.title = element_text(size = 1.5 * ctrl$base_size),
axis.line = element_blank(),
axis.text = element_text(size = ctrl$base_size),
axis.title = element_text(size = 1.5 * ctrl$base_size))
## Handle legends
if (ctrl$legend.position %in% c("bottomright", "bottomleft","topleft", "topright")) {
## legend inside
if (ctrl$legend.position == "bottomright") legendAdj <- c(1, 0)
if (ctrl$legend.position == "bottomleft") legendAdj <- c(0, 0)
if (ctrl$legend.position == "topright") legendAdj <- c(1, 1)
if (ctrl$legend.position == "topleft") legendAdj <- c(0, 1)
gg <- gg + theme(
legend.position = legendAdj, legend.justification = legendAdj,
legend.key = element_rect(fill = "transparent", colour = "transparent"),
legend.background = element_rect(fill = alpha(1e-5)))
} else {
## legend outside
gg <- gg + theme(legend.position = ctrl$legend.position,
legend.key = element_rect(fill = "white", color = "white"))
}
gg
}
stc04003/reReg documentation built on Sept. 20, 2023, 7:30 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions basebind plotEvents.control plotHaz plotRate plot.reReg plotMCF plotEvents plot.Recur
Documented in basebind plotEvents plotEvents.control plotHaz plotRate plot.Recur plot.reReg
globalVariables(c("time1", "time2", "group", "Y", "Y.upper", "Y.lower", "id", "event", "origin", "MCF", "bs"))
#' Produce Event Plot or Mean Cumulative Function Plot
#'
#' Plot the event plot or the mean cumulative function (MCF) from an \code{Recur} object.
#'
#' The argument \code{control} consists of options with argument defaults to a list with
#' the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is "Subject" for event plot and
#' "Cumulative mean" for MCF plot.}
#' \item{main}{customizable title, the default value is "Recurrent event plot"
#' when \code{mcf = FALSE} and
#' "Sample cumulative mean function plot" when \code{mcf = TRUE}.}
#' \item{terminal.name}{customizable label for terminal event,
#' the default value is "Terminal event".}
#' \item{recurrent.name}{customizable legend title for recurrent event,
#' the default value is "Recurrent events".}
#' \item{recurrent.types}{customizable label for recurrent event type,
#' the default value is \code{NULL}.}
#' \item{alpha}{between 0 and 1, controls the transparency of points.}
#' }
#' The \code{xlab}, \code{ylab} and \code{main} parameters can be specified
#' outside of the \code{control} list.
#'
#' @param x an object of class \code{Recur} returned by the \code{Recur()} function.
#' See \code{?Recur} for creating \code{Recur} objects.
#' @param event.result an optional character string that is passed to the
#' \code{plotEvents()} function as the \code{result} argument. See \code{\link{plotEvents}}.
#' This argument is used to specify whether the event plot is sorted by the subjects' terminal time.
#' The available options are
#' \describe{
#' \item{\code{increasing}}{sort the terminal time from in ascending order (default).
#' This places longer terminal times on top. }
#' \item{\code{decreasing}}{sort the terminal time from in descending order.
#' This places shorter terminal times on top. }
#' \item{\code{none}}{present the event plots as is, without sorting by the terminal times.}
#' }
#' @param event.calendarTime an optional logical value indicating whether to plot in calendar time.
#' When \code{event.calendarTime = FALSE} (default),
#' the event plot will have patient time on the x-axis.
#' @param control a list of control parameters. See \bold{Details}.
#' @param mcf.adjustRiskset an optional logical value that is passed to
#' the \code{mcf()} function as the \code{adjustRiskset} argument.
#' This argument indicates whether risk set size will be adjusted.
#' If \code{mcf.adjustRiskset = TRUE}, subjects leave the risk set after terminal times
#' as in the Nelson-Aalen estimator.
#' If \code{mcf.adjustRiskset = FALSE}, subjects remain in the risk set after terminal time.
#' @param mcf.conf.int an optional logical value that is passed to
#' the \code{mcf()} function as the \code{conf.int} argument. See \code{\link{mcf}} for details.
#' @param mcf an optional logical value indicating whether the mean cumulative function (MCF) will
#' be plotted instead of the event plot. When \code{mcf = TRUE},
#' the \code{mcf} is internally called. See \code{\link{mcf}} for details.
#' @param ... additional graphical parameters to be passed to methods.
#'
#' @seealso \code{\link{Recur}}, \code{\link{plotEvents}}, \code{\link{mcf}}
#'
#' @references Nelson, W. B. (1995) Confidence Limits for Recurrence Data-Applied to Cost
#' or Number of Product Repairs. \emph{Technometrics}, \bold{37}(2): 147--157.
#'
#' @keywords Plots
#' @export
#'
#' @return A \code{ggplot} object.
#' @example inst/examples/ex_plot_Recur.R
#' @method plot Recur
plot.Recur <- function(x, mcf = FALSE,
event.result = c("increasing", "decreasing", "asis"),
event.calendarTime = FALSE,
mcf.adjustRiskset = TRUE,
mcf.conf.int = FALSE,
control = list(), ...) {
event.result <- match.arg(event.result)
if (!is.Recur(x)) stop("Response must be a `Recur` object.")
ctrl <- plotEvents.control()
if (mcf & is.null(ctrl$ylab)) ctrl <- ctrl$ylab <- "Cumulative mean"
call <- match.call()
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) eval(call[[x]]))
}
if (!mcf) {
return(plotEvents(x, result = event.result, calendarTime = event.calendarTime, control = ctrl))
}
if (mcf) {
return(plot(mcf(x ~ 1, adjustRiskset = mcf.adjustRiskset), conf.int = mcf.conf.int))
## return(plotMCF(x, adjustRiskset = mcf.adjustRiskset, smooth = mcf.smooth, control = ctrl))
}
}
#' Produce Event Plots
#'
#' Plot the event plot for an \code{Recur} object.
#' The usage of the function is similar to that of \code{plot.Recur()} but with more flexible options.
#'
#' The argument \code{control} consists of options with argument defaults to a list with
#' the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is "Subject" for event plot and
#' "Cumulative mean" for MCF plot.}
#' \item{main}{customizable title, the default value is "Recurrent event plot"
#' when \code{mcf = FALSE} and
#' "Sample cumulative mean function plot" when \code{mcf = TRUE}.}
#' \item{terminal.name}{customizable label for terminal event,
#' the default value is "Terminal event".}
#' \item{recurrent.name}{customizable legend title for recurrent event,
#' the default value is "Recurrent events".}
#' \item{recurrent.types}{customizable label for recurrent event type,
#' the default value is \code{NULL}.}
#' \item{alpha}{between 0 and 1, controls the transparency of points.}
#' }
#' The \code{xlab}, \code{ylab} and \code{main} parameters can be specified
#' outside of the \code{control} list.
#'
#' @param formula a formula object, with the response on the left of a "~" operator,
#' and the predictors on the right.
#' The response must be a recurrent event survival object as returned by function \code{Recur()}.
#' @param data an optional data frame in which to interpret the variables occurring in
#' the "\code{formula}".
#' @param result an optional character string specifying whether the event plot is
#' sorted by the subjects' terminal time. The available options are
#' \describe{
#' \item{\code{increasing}}{sort the terminal time from in ascending order (default).
#' This places longer terminal times on top. }
#' \item{\code{decreasing}}{sort the terminal time from in descending order.
#' This places shorter terminal times on top. }
#' \item{\code{none}}{present the event plots as is, without sorting by the terminal times.}
#' }
#' @param calendarTime an optional logical value indicating whether to plot in calendar time.
#' When \code{calendarTime = FALSE} (default), the event plot will have patient time on the x-axis.
#' @param control a list of control parameters. See \bold{Details}.
#' @param ... graphical parameters to be passed to methods.
#' These include \code{xlab}, \code{ylab}, \code{main}, and more. See \bold{Details}.
#'
#' @seealso \code{\link{Recur}}, \code{\link{plot.Recur}}
#'
#' @keywords Plots
#' @export
#'
#' @return A \code{ggplot} object.
#' @importFrom ggplot2 geom_rect ggplot_build scale_y_continuous unit
#' @example inst/examples/ex_plot_event.R
plotEvents <- function(formula, data, result = c("increasing", "decreasing", "asis"),
calendarTime = FALSE, control = list(), ...) {
result <- match.arg(result)
ctrl <- plotEvents.control()
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
call <- match.call()
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) eval(call[[x]]))
}
nX <- 0
if (is.Recur(formula)) {
DF <- as.data.frame(formula@.Data)
isDate <- "Date" %in% formula@time_class
vNames <- NULL
} else {
if (missing(data)) obj <- eval(formula[[2]], parent.frame())
else obj <- eval(formula[[2]], data)
if (!is.Recur(obj)) stop("Response must be a `Recur` object.")
nX <- length(formula[[3]])
isDate <- "Date" %in% obj@time_class
if (formula[[3]] == 1) DF <- as.data.frame(obj@.Data)
if (formula[[3]] != 1 && nX == 1) {
if (missing(data)) DF <- data.frame(obj@.Data, eval(formula[[3]], parent.frame()))
if (!missing(data)) DF <- data.frame(obj@.Data, eval(formula[[3]], data))
colnames(DF) <- c(colnames(obj@.Data), paste0(formula[[3]], collapse = ""))
DF <- as.data.frame(DF)
}
if (formula[[3]] != 1 && nX > 1) {
DF <- as.data.frame(obj@.Data)
if (missing(data)) {
for (i in 2:nX) {
DF <- cbind(DF, eval(formula[[3]][[i]], parent.frame()))
}
} else {
for (i in 2:nX) {
DF <- cbind(DF, eval(formula[[3]][[i]], data))
}
}
vNames <- attr(terms(formula), "term.labels")
colnames(DF) <- c(colnames(obj@.Data), vNames)
}
vNames <- attr(terms(formula), "term.labels")
if (length(vNames) == 0) vNames <- NULL
}
## dat$status <- ifelse(is.na(dat$status), 0, dat$status)
## dat$Yi <- ifelse(is.na(dat$Yi), unlist(lapply(dat$tij, max)), dat$Yi)
newIDtime2 <- function(dat, result = "increasing") {
if (result == "asis") {
tmp <- table(dat$id)
dat$id <- rep(1:length(tmp), tmp[match(unique(dat$id), names(tmp))])
return(dat)
} else {
dat <- dat[order(dat$id),]
## if (all(dat$origin == 0))
if (!calendarTime)
tmp <- rank((dat$time2 - dat$origin)[dat$event == 0], ties.method = "first")
else
tmp <- rank(dat$time2[dat$event == 0], ties.method = "first")
}
if (result == "decreasing") tmp <- length(tmp) - tmp + 1
dat$id <- rep(tmp, table(dat$id))
return(dat)
}
if (nX == 0 || formula[[3]] == 1) {
DF <- newIDtime2(DF, result = result)
} else {
DF <- do.call(rbind,
lapply(split(DF, DF[, 7:ncol(DF)], drop = TRUE), newIDtime2, result = result))
rownames(DF) <- NULL
}
if (is.null(ctrl$cex)) sz <- 1 + 8 / (1 + exp(length(unique(DF$id)) / 30)) / max(1, nX)
else sz <- ctrl$cex
k <- length(unique(DF$event)) - 1 ## exclude event = 0
shp.val <- c(ctrl$terminal.shape, rep(ctrl$recurrent.shape, k))
if (is.null(ctrl$recurrent.color))
ctrl$recurrent.color <- hcl(h = seq(120, 360, length.out = k), l = 60, alpha = ctrl$alpha)
if (length(ctrl$recurrent.color) < k)
ctrl$recurrent.color <- alpha(rep(ctrl$recurrent.color, length.out = k), ctrl$alpha)
clr.val <- c(alpha(ctrl$terminal.color, ctrl$alpha), ctrl$recurrent.color)
if (k == 0) rec.lab <- NULL
else rec.lab <- paste("r", 1:k, sep = "")
if (k == 0) shp.lab <- ctrl$terminal.name
if (k == 1) shp.lab <- c(ctrl$terminal.name, ctrl$recurrent.name)
if (k > 1 & is.null(ctrl$recurrent.type)) {
shp.lab <- c(ctrl$terminal.name, paste(ctrl$recurrent.name, 1:k))
}
if (k > 1 & !is.null(ctrl$recurrent.type)) {
if (length(ctrl$recurrent.type) == k) {
shp.lab <- c(ctrl$terminal.name, ctrl$recurrent.type)
} else {
message('The length of "recurrent.type" mismatched, default names are used.')
shp.lab <- c(ctrl$terminal.name, paste(ctrl$recurrent.name, 1:k))
}
}
if (nX > 0) {
for (i in vNames) {
DF[,i] <- factor(DF[,i], labels = paste(i, "=", unique(DF[,i])))
}}
names(shp.val) <- names(clr.val) <- c("terminal", rec.lab)
## Bars
if (calendarTime) {
if (is.null(ctrl$width)) ctrl$width <- 0
gg <- ggplot(DF, aes(xmin = id - .45 - ctrl$width,
xmax = id + .45 + ctrl$width,
ymin = time1, ymax = time2)) +
geom_rect(fill = ctrl$bar.color) +
coord_flip()
} else gg <- ggplot(DF[DF$event == 0,], aes(id, time2 - origin)) +
geom_bar(stat = "identity", fill = ctrl$bar.color, width = ctrl$width) +
coord_flip()
## event dots
if (!calendarTime) DF$time2 <- DF$time2 - DF$origin
if (any(table(DF$id) > 0))
gg <- gg + geom_point(data = DF[DF$event > 0,],
aes(id, time2,
shape = factor(event, labels = rec.lab),
color = factor(event, labels = rec.lab),
stroke = ctrl$recurrent.stroke),
size = sz)
if (sum(DF$terminal, na.rm = TRUE) > 0)
gg <- gg + geom_point(data = DF[DF$terminal > 0,],
aes(id, time2, shape = "terminal", color = "terminal",
stroke = ctrl$terminal.stroke),
size = sz)
if (!is.null(ctrl$not.terminal.shape))
gg <- gg + geom_point(data = DF[DF$terminal == 0,], aes(id, time2),
shape = ctrl$not.terminal.shape,
color = ctrl$not.terminal.color,
stroke = ctrl$terminal.stroke,
size = sz)
if (nX > 0 && formula[[3]] != 1)
gg <- gg + facet_grid(as.formula(paste(formula[3], "~.", collapse = "")),
scales = "free", space = "free", switch = "both")
## Add theme and final touch ups
if (ctrl$main != "") gg <- gg + ggtitle(ctrl$main)
gg <- gg + scale_shape_manual(name = "", values = shp.val,
labels = shp.lab, breaks = c("terminal", rec.lab)) +
scale_color_manual(name = "", values = clr.val,
labels = shp.lab, breaks = c("terminal", rec.lab)) +
theme(panel.background = element_blank(),
axis.line = element_line(color = "black"),
axis.line.y = element_blank(),
axis.title.y = element_text(vjust = 0),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.title = element_text(size = 2 * ctrl$base_size),
strip.text = element_text(size = ctrl$base_size),
legend.text = element_text(size = 1.5 * ctrl$base_size),
legend.title = element_text(size = 1.5 * ctrl$base_size),
axis.text = element_text(size = ctrl$base_size),
axis.title = element_text(size = 1.5 * ctrl$base_size)) +
scale_x_continuous(expand = c(0, 1)) +
labs(x = ctrl$ylab, y = ctrl$xlab) +
guides(shape = guide_legend(override.aes = list(size = 2.7)))
## Handle legends
if (ctrl$legend.position %in% c("bottomright", "bottomleft","topleft", "topright")) {
## legend inside
if (ctrl$legend.position == "bottomright") legendAdj <- c(1, 0)
if (ctrl$legend.position == "bottomleft") legendAdj <- c(0, 0)
if (ctrl$legend.position == "topright") legendAdj <- c(1, 1)
if (ctrl$legend.position == "topleft") legendAdj <- c(0, 1)
gg <- gg + theme(
legend.position = legendAdj, legend.justification = legendAdj,
legend.key = element_rect(fill = "transparent", colour = "transparent"),
legend.background = element_rect(fill = alpha(1e-5)))
} else {
## legend outside
gg <- gg + theme(legend.position = ctrl$legend.position,
legend.key = element_rect(fill = "white", color = "white"))
}
if (isDate & calendarTime) {
xl <- ggplot_build(gg)$layout$panel_params[[1]]$x$breaks
xl <- xl[!is.na(xl)]
gg <- gg + scale_y_continuous(breaks = xl, labels = as.Date(xl, origin = "1970-01-01"))
}
gg
}
#' Produce Cumulative Sample Mean Function Plots
#'
#' Plot the mean cumulative function (MCF) for an \code{Recur} object.
#' The usage of the function is similar to that of \code{plot.Recur()}
#' but with more flexible options.
#'
#' When \code{adjustRiskset = TRUE}, the \code{plotMCF()} is equivalent to
#' the Nelson-Aalen estimator for the intensity function of the recurrent event process.
#' When \code{adjustRiskset = FALSE}, the \code{plotMCF()} does not adjust for the risk set and
#' assumes all subjects remain at risk after the last observed recurrent event.
#' This is known as the survivor rate function.
#' The argument \code{control} consists of options with argument defaults
#' to a list with the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is "Cumulative mean".}
#' \item{main}{customizable title, default value is "Sample cumulative mean function plot".}
#' }
#' The \code{xlab}, \code{ylab} and \code{main} parameters can also be
#' specified outside of the \code{control} list.
#'
#' @param formula a formula object, with the response on the left of a "~" operator,
#' and the predictors on the right.
#' The response must be a recurrent event survival object returned by the \code{Recur()} function.
#' @param data an optional data frame in which to interpret the variables occurring in
#' the "\code{formula}".
#' @param adjustRiskset an optional logical value that is passed to the \code{plotMCF()} function
#' as the \code{adjustRiskset} argument. See \code{\link{plotMCF}}.
#' This argument indicates whether risk set size will be adjusted. If \code{mcf.adjustRiskset = TRUE},
#' subjects leave the risk set after terminal times as in the Nelson-Aalen estimator.
#' If \code{mcf.adjustRiskset = FALSE}, subjects remain in the risk set after terminal time.
#' @param smooth an optional logical value indicating whether to add a smooth curve
#' obtained from a monotone increasing P-splines implemented in package \code{scam}.
#' This feature only works for data with one recurrent event type.
#' @param onePanel an optional logical value indicating whether the mean cumulative functions
#' will be plotted in the same panel.
#' This is only useful when there are multiple recurrent event types or
#' in the presence of (discrete) covariates.
#' @param control a list of control parameters.
#' @param ... graphical parameters to be passed to methods.
#' These include \code{xlab}, \code{ylab}, \code{main}, and more. See \bold{Details}.
#'
#' @seealso \code{\link{Recur}}, \code{\link{plot.Recur}}
#' @keywords Plots
#'
#' @return A \code{ggplot} object.
#'
#' @importFrom ggplot2 guides guide_legend
#' @importFrom scam scam
#'
#' @noRd
#'
#' @example inst/examples/ex_plot_MCF.R
plotMCF <- function(formula, data, adjustRiskset = TRUE, onePanel = FALSE,
smooth = FALSE, control = list(), ...) {
call <- match.call()
ctrl <- plotEvents.control()
if (is.null(ctrl$ylab)) ctrl$ylab <- "Cumulative mean"
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) eval(call[[x]]))
}
nX <- 0
if(is.Recur(formula)) {
DF <- as.data.frame(formula@.Data)
vNames <- NULL
} else {
if (missing(data)) obj <- eval(formula[[2]], parent.frame())
else obj <- eval(formula[[2]], data)
if (!is.Recur(obj)) stop("Response must be a `Recur` object.")
DF <- as.data.frame(obj@.Data)
nX <- length(formula[[3]])
if (formula[[3]] != 1 && nX == 1) {
if (missing(data)) DF <- cbind(obj@.Data, tmp = eval(formula[[3]], parent.frame()))
if (!missing(data)) DF <- cbind(obj@.Data, tmp = eval(formula[[3]], data))
colnames(DF) <- c(colnames(obj@.Data), paste0(formula[[3]], collapse = ""))
DF <- as.data.frame(DF)
}
if (formula[[3]] != 1 && nX > 1) {
DF <- as.data.frame(obj@.Data)
if (missing(data)) {
for (i in 2:nX) {
DF <- cbind(DF, eval(formula[[3]][[i]], parent.frame()))
}
} else {
for (i in 2:nX) {
DF <- cbind(DF, eval(formula[[3]][[i]], data))
}
}
colnames(DF) <- c(colnames(obj@.Data), attr(terms(formula), "term.labels"))
}
vNames <- attr(terms(formula), "term.labels")
if (length(vNames) == 0) vNames <- NULL
}
## dd <- subset(DF, select = c("event", vNames, "time2"))
dd <- DF[,c("event", vNames, "time2")]
dd <- dd[do.call(order, dd),]
nn <- unlist(aggregate(time2~., data = dd, table)$time2)
dd <- unique(dd)
dd$n <- as.integer(nn)
rownames(dd) <- NULL
k <- length(unique(dd$event)) - 1
if (!is.null(vNames)) { ## any covariates/stratifications?
dd2 <- DF[DF$event == 0, vNames, drop = FALSE]
dd2 <- dd2[do.call(order, dd2),,drop = FALSE]
nn <- c(t(table(dd2)))
dd2 <- unique(dd2)
dd2$n <- as.integer(nn) ## tmp2 in the 1st version
rownames(dd2) <- NULL
dd$GrpInd <- match(apply(dd[,vNames, drop = FALSE], 1, paste, collapse = ""),
apply(dd2[,vNames, drop = FALSE], 1, paste, collapse = ""))
tmp1 <- merge(dd, dd2, by = vNames)
## as.integer(eval(parse(text = paste(attr(terms(formula), "term.labels"), collapse = ":"))))
rec0 <- tmp1[tmp1$event == 0,]
dat0 <- do.call(rbind, lapply(split(tmp1, tmp1$GrpInd), function(x) {
x$adjustRiskset = apply(x, 1, function(y)
as.numeric(y['n.y']) - sum(rec0$n.x[as.numeric(y['time2'])> rec0$time2 &
rec0$GrpInd == as.numeric(y['GrpInd'])]))
return(x)}))
dat0$n.x <- dat0$n.x * (dat0$event > 0)
rec0$time2 <- rec0$n.x <- 0
rec0$adjustRiskset <- 1
dat0 <- unique(rbind(dat0, rec0))
rownames(dat0) <- NULL
## Number of recurrent types
if (k > 1) {
tmp <- dat0[dat0$event == 0,]
tmp <- tmp[rep(1:NROW(tmp), k),]
tmp$event <- rep(1:k, each = sum(dat0$event == 0))
dat0 <- rbind(dat0[dat0$event > 0,], tmp)
} else {
dat0$event <- dat0$event[1]
}
dat0 <- dat0[order(dat0$event, dat0$GrpInd, dat0$time2),]
if (adjustRiskset) {
dat0 <- do.call(rbind,
lapply(split(dat0, list(dat0$event, dat0$GrpInd)), function(x) {
x$mu <- x$n.x / x$adjustRiskset
x$MCF <- cumsum(x$mu)
return(x)}))
} else {
dat0 <- do.call(rbind,
lapply(split(dat0, list(dat0$event, dat0$GrpInd)), function(x) {
x$mu <- x$n.x / x$n.y
x$MCF <- cumsum(x$mu)
return(x)}))
}
} else { ## no covariates
dd$n.y <- length(unique(DF$id))
rec0 <- dd[dd$event == 0, ]
dd$adjustRiskset <- apply(dd, 1, function(x) x[4] - sum(rec0$n[x[2] > rec0$time2]))
dd$n <- dd$n * (dd$event > 0)
rec0$time2 <- rec0$n <- 0
rec0$adjustRiskset <- 1
dat0 <- unique(rbind(dd, rec0))
dat0 <- dat0[order(dat0$time2),]
dat0$mu <- dat0$n / (adjustRiskset * dat0$adjustRiskset + (!adjustRiskset) * dat0$n.y)
dat0$MCF <- cumsum(dat0$mu)
}
if (k == 1) {
dat0$event <- dat0$event[1]
dat0$event <- factor(dat0$event, labels = ctrl$recurrent.name)
}
if (k > 1) {
dat00 <- dat0[dat0$event == 0,]
dat0 <- dat0[dat0$event > 0,]
if (nrow(dat00) > 0) {
dat00 <- dat00[rep(1:nrow(dat00), k),]
dat00$event <- rep(1:k, each = nrow(dat00) / k)
dat0 <- rbind(dat0, dat00)
}
if (is.null(ctrl$recurrent.type))
dat0$event <- factor(dat0$event, labels = paste(ctrl$recurrent.name, 1:k))
if (!is.null(ctrl$recurrent.type)) {
if (length(ctrl$recurrent.type) == k) {
dat0$event <- factor(dat0$event, labels = ctrl$recurrent.type)
} else {
message('The length of "recurrent.type" mismatched, default names are used.')
dat0$event <- factor(dat0$event, labels = paste(ctrl$recurrent.name, 1:k))
}
}
}
if (nX > 0) {
for (i in vNames) {
dat0[,i] <- factor(dat0[,i], labels = paste(i, "=", unique(dat0[,i])))
}}
dat0 <- dat0[complete.cases(dat0),]
gg <- ggplot(data = dat0, aes(x = time2, y = MCF))
if (is.null(vNames) & k == 1) {
gg <- gg + geom_step(size = ctrl$lwd)
} else {
if (!onePanel & k == 1) gg <- gg + geom_step(size = ctrl$lwd)
if (!onePanel & k > 1)
gg <- gg + geom_step(aes(color = event), direction = "hv", size = ctrl$lwd) +
guides(color = guide_legend(title = ctrl$recurrent.name))
if (onePanel) {
dat0$GrpInd <- factor(dat0$GrpInd)
levels(dat0$GrpInd) <- apply(unique(dat0[,vNames, drop = FALSE]), 1, paste, collapse = ", ")
gg <- gg + geom_step(aes(color = dat0$GrpInd), direction = "hv", size = ctrl$lwd) +
guides(color = guide_legend(title = ""))
}
}
if (!onePanel && nX > 0 && formula[[3]] != 1)
gg <- gg + facet_grid(as.formula(paste(formula[3], "~.", collapse = "")),
scales = "free", space = "free", switch = "x")
if (!onePanel & k == 1)
gg <- gg + theme(legend.position = "none")
## if (onePanel & k == 1)
## gg <- gg + scale_color_discrete(name = "", labels = levels(interaction(vNames)))
## if (onePanel & k > 1) gg <- gg + scale_color_discrete(name = "")
if (smooth & k == 1 & !onePanel) {
if (is.null(dat0$GrpInd)) dat0$GrpInd <- 1
dat0 <- do.call(rbind, lapply(split(dat0, dat0$GrpInd), function(x){
x$bs <- scam(x$MCF ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values
return(x)}))
gg <- gg + geom_line(data = dat0, aes(time2, y = bs), color = 4, size = ctrl$lwd)
## geom_smooth(method = "scam", formula = y ~ s(x, k = 10, bs = "mpi"), size = ctrl$lwd, se = FALSE)
}
## gg <- gg + geom_smooth(method = "loess", size = ctrl$lwd, se = FALSE)
if (smooth & k > 1) message('Smoothing only works for data with one recurrent event type.')
if (ctrl$main != "") gg <- gg + ggtitle(ctrl$main)
gg <- gg + theme(axis.line = element_line(color = "black"),
plot.title = element_text(size = 2 * ctrl$base_size),
strip.text = element_text(size = ctrl$base_size),
legend.text = element_text(size = 1.5 * ctrl$base_size),
legend.title = element_text(size = 1.5 * ctrl$base_size),
axis.text = element_text(size = ctrl$base_size),
axis.title = element_text(size = 1.5 * ctrl$base_size)) +
labs(y = ctrl$ylab, x = ctrl$xlab)
## Handle legends
if (ctrl$legend.position %in% c("bottomright", "bottomleft","topleft", "topright")) {
## legend inside
if (ctrl$legend.position == "bottomright") legendAdj <- c(1, 0)
if (ctrl$legend.position == "bottomleft") legendAdj <- c(0, 0)
if (ctrl$legend.position == "topright") legendAdj <- c(1, 1)
if (ctrl$legend.position == "topleft") legendAdj <- c(0, 1)
gg <- gg + theme(
legend.position = legendAdj, legend.justification = legendAdj,
legend.key = element_rect(fill = "transparent", colour = "transparent"),
legend.background = element_rect(fill = alpha(1e-5)))
} else {
## legend outside
gg <- gg + theme(legend.position = ctrl$legend.position,
legend.key = element_rect(fill = "white", color = "white"))
}
gg
}
#' Plot the Baseline Cumulative Rate Function and the Baseline Cumulative Hazard Function
#'
#' Plot the baseline cumulative rate function and the baseline cumulative hazard function
#' (if applicable) for an \code{reReg} object.
#'
#' The argument \code{control} consists of options with argument defaults to a list
#' with the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is empty.}
#' \item{main}{customizable title, default value are "Baseline cumulative rate and
#' hazard function" when \code{baseline = "both"},
#' "Baseline cumulative rate function" when \code{baseline = "rate"},
#' and "Baseline cumulative hazard function" when \code{baseline = "hazard"}.}
#' }
## #' These arguments can also be passed down without specifying a \code{control} list.
#'
#' @param x an object of class \code{reReg}, returned by the \code{reReg} function.
#' @param smooth an optional logical value indicating whether to add a smooth curve
#' obtained from a monotone increasing P-splines implemented in package \code{scam}.
#' @param baseline a character string specifying which baseline function to plot.
#' \describe{
#' \item{\code{baseline = "both"}}{plot both the baseline cumulative rate and
#' the baseline cumulative hazard function (if applicable) in separate panels
#' within the same display (default).}
#' \item{\code{baseline = "rate"}}{plot the baseline cumulative rate function.}
#' \item{\code{baseline = "hazard"}}{plot the baseline cumulative hazard function.}
#' }
## #' @param rateType a character string specifying the type of rate function to be plotted.
## #' This argument is passed to the \code{plotRate()} function as the \code{type} argument.
## #' See \code{\link{plotRate}}.
## #' Options are "unrestricted", "scaled", "bounded".
#' @param newdata an optional data frame contains variables to include in the calculation
#' of the cumulative rate function.
#' If omitted, the baseline rate function will be plotted.
#' @param frailty an optional vector to specify the shared frailty for \code{newdata}.
#' If \code{newdata} is given and \code{frailty} is not specified, the
#' @param showName an optional logical value indicating whether to label the curves
#' when \code{newdata} is specified.
#' @param control a list of control parameters. See \bold{Details}.
#' @param ... additional graphical parameters to be passed to methods.
#'
#' @seealso \code{\link{reReg}}
#' @export
#' @keywords Plots
#'
#' @return A \code{ggplot} object.
#'
#' @importFrom ggplot2 geom_smooth geom_step ggplotGrob
#' @importFrom grid grid.draw
#' @example inst/examples/ex_plot_reReg.R
#' @method plot reReg
plot.reReg <- function(x,
baseline = c("both", "rate", "hazard"),
## type = c("unrestricted", "bounded", "scaled"),
smooth = FALSE, newdata = NULL, frailty = NULL, showName = FALSE,
control = list(), ...) {
baseline <- match.arg(baseline)
## type <- match.arg(type)
type <- "unrestricted"
if (x$typeRec %in% c("cox.GL", "cox.LWYY", "am.GL"))
stop("Baseline functions not yet available for this method.")
if (x$typeRec %in% c("cox.LWYY", "cox.HH"))
ctrl <- plotEvents.control(ylab = "Rate")
if (baseline %in% c("both", "rate"))
ctrl <- plotEvents.control(ylab = "Rate")
if (baseline == "hazard")
ctrl <- plotEvents.control(ylab = "Hazard")
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
call <- match.call()
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) call[[x]])
}
if (!is.reReg(x)) stop("Response must be a `reReg` class")
if (baseline == "rate")
return(plotRate(x, smooth = smooth, type = type,
newdata = newdata, frailty = frailty, showName = showName, control = ctrl))
if (baseline == "hazard")
return(plotHaz(x, smooth = smooth,
newdata = newdata, frailty = frailty, showName = showName, control = ctrl))
if (x$typeTem == ".") {
## cat(paste("Baseline cumulative hazard function is not available."))
## cat("\nOnly the baseline cumulative rate function is plotted.\n")
return(plotRate(x, smooth = smooth, type = type,
newdata = newdata, frailty = frailty, showName = showName, control = ctrl))
}
g1 <- plotRate(x, smooth = smooth, type = type,
newdata = newdata, frailty = frailty, showName = showName, control = ctrl)
g2 <- plotHaz(x, smooth = smooth, type = type,
newdata = newdata, frailty = frailty, showName = showName, control = ctrl)
if (ctrl$main != "") g1 <- g1 + ggtitle(ctrl$main)
g1 <- g1 + ylab("Rate") + xlab("") +
## facet_grid(factor(rep(1, nrow(x$DF)), levels = 1, labels = "Baseline cumulative rate")) +
theme(axis.title.x = element_blank(),
axis.ticks.x = element_blank(),
plot.margin = unit(c(0, 0, 0, 0), "cm"))
g2 <- g2 + ylab("Hazard") + xlab("Time") +
## facet_grid(factor(rep(1, nrow(x$DF)), levels = 1, labels = "Baseline cumulative hazard")) +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
g1 <- ggplotGrob(g1)
g2 <- ggplotGrob(g2)
grid.draw(rbind(g1, g2))
## gg + ggtitle(ctrl$main)
}
#' Plotting the Baseline Cumulative Rate Function for the Recurrent Event Process
#'
#' Plot the baseline cumulative rate function for an \code{reReg} object.
#'
#'
#' The \code{plotRate()} plots the estimated baseline cumulative rate function
#' depending on the identifiability assumption.
#' When \code{type = "unrestricted"} (default), the baseline cumulative rate function
#' is plotted under the assumption \eqn{E(Z) = 1}.
#' When \code{type = "scaled"}, the baseline cumulative rate function is plotted
#' under the assumption \eqn{\Lambda(\min(Y^\ast, \tau)) = 1}.
#' When \code{type = "bounded"}, the baseline cumulative rate function is plotted
#' under the assumption \eqn{\Lambda(\tau) = 1}.
#' See \code{?reReg} for the specification of the notations and underlying models.
#'
#' The argument \code{control} consists of options with argument defaults
#' to a list with the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is empty.}
#' \item{main}{customizable title, default value is "Baseline cumulative rate function".}
#' }
#' These arguments can also be specified outside of the \code{control} list.
#'
#' @param x an object of class \code{reReg}, usually returned by the \code{reReg} function.
#' @param type a character string specifying the type of rate function to be plotted.
#' Options are "unrestricted", "scaled", "bounded". See \bold{Details}.
#' @param smooth an optional logical value indicating whether to add a smooth curve
#' obtained from a monotone increasing P-splines implemented in package \code{scam}.
#' @param newdata an optional data frame contains variables to include in the calculation
#' of the cumulative rate function.
#' If omitted, the baseline rate function will be plotted.
#' @param frailty an optional vector to specify the shared frailty for \code{newdata}.
#' If \code{newdata} is given and \code{frailty} is not specified, the
#' @param showName an optional logical value indicating whether to label the curves
#' when \code{newdata} is specified.
#' @param control a list of control parameters.
#' @param ... graphical parameters to be passed to methods.
#' These include \code{xlab}, \code{ylab}, \code{main}, and more. See \bold{Details}.
#'
#' @seealso \code{\link{reReg}} \code{\link{plot.reReg}}
#' @export
#'
#' @return A \code{ggplot} object.
#'
#' @keywords Plots
#' @importFrom ggplot2 scale_x_continuous
#' @importFrom directlabels geom_dl
#'
#' @example inst/examples/ex_plot_rate.R
plotRate <- function(x, newdata = NULL, frailty = NULL, showName = FALSE,
type = c("unrestricted", "bounded", "scaled"),
smooth = FALSE, control = list(), ...) {
if (x$typeRec %in% c("cox.GL", "cox.LWYY", "am.GL"))
stop("Baseline cumulative rate function is not available")
if (is.null(frailty)) frailty <- exp(x$log.muZ)
if (length(frailty) > 1 & !is.null(newdata) && length(frailty) != nrow(newdata))
stop("newdata and frailty are different lengths")
## ctrl <- plot.reReg.control(main = "Baseline cumulative rate function")
ctrl <- plotEvents.control(ylab = "Rate")
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
call <- match.call()
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) call[[x]])
}
type <- match.arg(type)
if (x$typeRec == "nonparametric") {
type <- "bounded"
ctrl$ylab = "MCF Estimates"
}
if (!is.reReg(x)) stop("Response must be a `reReg` class")
dat <- x$DF[,"time2",drop = FALSE]
if (is.null(newdata)) {
if (type == "unrestricted") dat$Y <- x$Lam0(dat$time2) * exp(x$log.muZ)
if (type == "scaled") dat$Y <- x$Lam0(dat$time2) / x$Lam0(max(dat$time2))
if (type == "bounded") dat$Y <- x$Lam0(dat$time2)
if (!is.null(x$Lam0.upper)) {
if (type == "bounded") {
dat$Y.upper <- x$Lam0.upper(dat$time2)
dat$Y.lower <- x$Lam0.lower(dat$time2)
}
if (type == "unrestricted") {
dat$Y.upper <- x$Lam0.upper(dat$time2) * exp(x$log.muZ)
dat$Y.lower <- x$Lam0.lower(dat$time2) * exp(x$log.muZ)
}
if (type == "scaled") {
dat$Y.upper <- x$Lam0.upper(dat$time2) / x$Lam0(max(dat$time2))
dat$Y.lower <- x$Lam0.lower(dat$time2) / x$Lam0(max(dat$time2))
}
}
gg <- ggplot(data = dat, aes(x = time2, y = Y)) +
theme(axis.line = element_line(color = "black"))
if (smooth) {
dat$bs <- scam(dat$Y ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
gg <- gg + geom_line(aes(x = time2, y = dat$bs), color = 4)
if (!is.null(x$Lam0.upper)) {
dat$bs.upper <- scam(dat$Y.upper ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
dat$bs.lower <- scam(dat$Y.lower ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
gg <- gg + geom_line(aes(time2, y = dat$bs.upper), color = 4, lty = 2) +
geom_line(aes(time2, y = dat$bs.lower), color = 4, lty = 2)
}
} else {
gg <- gg + geom_step()
if (!is.null(x$Lam0.upper))
gg <- gg + geom_step(aes(x = time2, y = Y.upper), lty = 2) +
geom_step(aes(x = time2, y = Y.lower), lty = 2)
}
}
if (!is.null(newdata)) {
if (!is.null(x$xlevels)) {
for(i in which(names(newdata) %in% names(x$xlevels)))
newdata[,i] <- factor(newdata[,i], levels = x$xlevels[[i]])
newdata <- model.matrix(~., newdata)
}
X <- as.matrix(unique(newdata[,match(x$varNames, colnames(newdata)), drop = FALSE]))
if (ncol(X) != length(x$varNames))
stop(paste0("Variables ",
paste(setdiff(x$varNames, colnames(newdata)), collapse = ", "),
" are missing"))
p <- ncol(X)
exa1 <- exa2 <- 1
if (x$typeRec == "cox") exa2 <- exp(X %*% x$par1)
if (x$typeRec == "ar") exa1 <- exp(X %*% x$par1)
if (x$typeRec == "am") exa1 <- exa2 <- exp(X %*% x$par1)
if (x$typeRec == "sc") {
exa1 <- exp(X %*% x$par1)
exa2 <- exp(X %*% x$par2)
}
exa1 <- rep(drop(exa1), each = nrow(dat))
exa2 <- rep(drop(exa2), each = nrow(dat))
Y <- frailty * x$Lam0(dat$time2 * exa1) * exa2 / exa1
if (!is.null(x$Lam0.upper)) {
Y.upper <- frailty * x$Lam0.upper(dat$time2 * exa1) * exa2 / exa1
Y.lower <- frailty * x$Lam0.lower(dat$time2 * exa1) * exa2 / exa1
dat <- data.frame(time2 = dat$time2,
id = rep(rownames(X), each = nrow(dat)),
Y = Y, Y.upper = Y.upper, Y.lower = Y.lower)
} else
dat <- data.frame(time2 = dat$time2, id = rep(rownames(X), each = nrow(dat)), Y = Y)
gg <- ggplot(data = dat, aes(x = time2, y = Y, group = id)) +
theme(axis.line = element_line(color = "black"))
if (smooth) {
dat$bs <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
gg <- gg + geom_line(aes(x = time2, y = dat$bs, group = id), color = 4)
if (!is.null(x$Lam0.upper)) {
dat$bs.upper <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y.upper ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
dat$bs.lower <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y.lower ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
gg <- gg +
geom_line(aes(x = time2, y = dat$bs.upper, group = id), color = 4, lty = 2) +
geom_line(aes(x = time2, y = dat$bs.lower, group = id), color = 4, lty = 2)
}
} else {
gg <- gg + geom_step()
if (!is.null(x$Lam0.upper))
gg <- gg +
geom_step(aes(x = time2, y = Y.upper, group = id), lty = 2) +
geom_step(aes(x = time2, y = Y.lower, group = id), lty = 2)
}
if (showName)
gg <- gg + geom_dl(aes(label = paste(" Obs. =", id)), method = "right.polygons") +
scale_x_continuous(limits = c(0, max(dat$time2) * 1.1))
}
if (ctrl$main != "") gg <- gg + ggtitle(ctrl$main)
gg <- gg + labs(x = ctrl$xlab, y = ctrl$ylab) +
theme(plot.title = element_text(size = 2 * ctrl$base_size),
strip.text = element_text(size = ctrl$base_size),
legend.text = element_text(size = 1.5 * ctrl$base_size),
legend.title = element_text(size = 1.5 * ctrl$base_size),
axis.line = element_blank(),
axis.text = element_text(size = ctrl$base_size),
axis.title = element_text(size = 1.5 * ctrl$base_size))
attr(gg, "from") <- "reReg"
return(gg)
}
#' Plot the Baseline Cumulative Hazard Function for the Terminal Time
#'
#' Plot the baseline cumulative hazard function for an \code{reReg} object.
#' The 95\% confidence interval on the baseline cumulative rate function
#'
#' The argument \code{control} consists of options with argument
#' defaults to a list with the following values:
#' \describe{
#' \item{xlab}{customizable x-label, default value is "Time".}
#' \item{ylab}{customizable y-label, default value is empty.}
#' \item{main}{customizable title, default value is "Baseline cumulative hazard function".}
#' }
#' These arguments can also be passed down without specifying a \code{control} list.
#'
#' @param x an object of class \code{reReg}, returned by the \code{reReg} function.
#' @param smooth an optional logical value indicating whether to add a smooth curve
#' obtained from a monotone increasing P-splines implemented in package \code{scam}.
#' @param newdata an optional data frame contains variables to include in the calculation
#' of the cumulative rate function.
#' If omitted, the baseline rate function will be plotted.
#' @param frailty an optional vector to specify the shared frailty for \code{newdata}.
#' If \code{newdata} is given and \code{frailty} is not specified, the
#' @param showName an optional logical value indicating whether to label the curves
#' when \code{newdata} is specified.
#' @param type a character string specifying the type of rate function to be plotted.
#' Options are "unrestricted", "scaled", "bounded". See \bold{Details}.
#' @param control a list of control parameters.
#' @param ... graphical parameters to be passed to methods.
#' These include \code{xlab}, \code{ylab}, \code{main}, and more. See \bold{Details}.
#'
#' @seealso \code{\link{reReg}} \code{\link{plot.reReg}}
#' @export
#'
#' @return A \code{ggplot} object.
#' @keywords Plots
#'
#' @example inst/examples/ex_plot_Haz.R
plotHaz <- function(x, newdata = NULL, frailty = NULL, showName = FALSE,
type = c("unrestricted", "bounded", "scaled"),
smooth = FALSE, control = list(), ...) {
if (x$typeRec %in% c("cox.GL", "cox.LWYY", "am.GL"))
stop("Baseline cumulative hazard function is not available.")
if (x$typeTem == ".") {
stop("Baseline cumulative hazard function is not available.")
}
if (is.null(frailty)) frailty <- exp(x$log.muZ)
if (length(frailty) > 1 & !is.null(newdata) && length(frailty) != nrow(newdata))
stop("newdata and frailty are different lengths")
## ctrl <- plot.reReg.control(main = "Baseline cumulative hazard function")
ctrl <- plotEvents.control(ylab = "Hazard")
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
call <- match.call()
namp <- names(match.call())
if (any(namp %in% names(ctrl))) {
namp <- namp[namp %in% names(ctrl)]
ctrl[namp] <- lapply(namp, function(x) call[[x]])
}
type <- match.arg(type)
if (!is.reReg(x)) stop("Response must be a `reReg` class.")
dat <- x$DF[, "time2", drop = FALSE]
if (is.null(newdata)) {
if (type == "unrestricted") dat$Y <- x$Haz0(dat$time2) * exp(x$log.muZ)
if (type == "scaled") dat$Y <- x$Haz0(dat$time2) / x$Haz0(max(dat$time2))
if (type == "bounded") dat$Y <- x$Haz0(dat$time2)
## dat$Y <- x$Haz0(dat$time2)
if (!is.null(x$Haz0.upper)) {
if (type == "bounded") {
dat$Y.upper <- x$Haz0.upper(dat$time2)
dat$Y.lower <- x$Haz0.lower(dat$time2)
}
if (type == "unrestricted") {
dat$Y.upper <- x$Haz0.upper(dat$time2) * exp(x$log.muZ)
dat$Y.lower <- x$Haz0.lower(dat$time2) * exp(x$log.muZ)
}
if (type == "scaled") {
dat$Y.upper <- x$Haz0.upper(dat$time2) / x$Haz0(max(dat$time2))
dat$Y.lower <- x$Haz0.lower(dat$time2) / x$Haz0(max(dat$time2))
}
}
gg <- ggplot(data = dat, aes(x = time2, y = Y)) +
theme(axis.line = element_line(color = "black"))
if (smooth) {
dat$bs <- scam(dat$Y ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
gg <- gg + geom_line(aes(time2, y = dat$bs), color = 4)
if (!is.null(x$Lam0.upper)) {
dat$bs.upper <- scam(dat$Y.upper ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
dat$bs.lower <- scam(dat$Y.lower ~ s(dat$time2, k = 10, bs = "mpi"))$fitted.values
gg <- gg + geom_line(aes(time2, y = dat$bs.upper), color = 4, lty = 2) +
geom_line(aes(time2, y = dat$bs.lower), color = 4, lty = 2)
}
} else {
gg <- gg + geom_step()
if (!is.null(x$Lam0.upper))
gg <- gg + geom_step(aes(x = time2, y = Y.upper), lty = 2) +
geom_step(aes(x = time2, y = Y.lower), lty = 2)
}
}
if (!is.null(newdata)) {
X <- as.matrix(unique(newdata[,match(x$varNames, names(newdata)), drop = FALSE]))
if (ncol(X) != length(x$varNames))
stop(paste0("Variables ",
paste(setdiff(x$varNames, names(newdata)), collapse = ", "),
" are missing"))
p <- ncol(X)
exb1 <- exb2 <- 1
if (x$typeTem == "cox") exb2 <- exp(X %*% x$par3)
if (x$typeTem == "ar") exb1 <- exp(X %*% x$par3)
if (x$typeTem == "am") exb1 <- exb2 <- exp(X %*% x$par3)
if (x$typeTem == "sc") {
exb1 <- exp(X %*% x$par3)
exb2 <- exp(X %*% x$par4)
}
exb1 <- rep(drop(exb1), each = nrow(dat))
exb2 <- rep(drop(exb2), each = nrow(dat))
Y <- frailty * x$Haz0(dat$time2 * exb1) * exb2 / exb1
if (!is.null(x$Haz0.upper)) {
Y.upper <- frailty * x$Haz0.upper(dat$time2 * exb1) * exb2 / exb1
Y.lower <- frailty * x$Haz0.lower(dat$time2 * exb1) * exb2 / exb1
dat <- data.frame(time2 = dat$time2, id = rep(rownames(X), each = nrow(dat)),
Y = Y, Y.upper = Y.upper, Y.lower = Y.lower)
} else
dat <- data.frame(time2 = dat$time2, id = rep(rownames(X), each = nrow(dat)), Y = Y)
gg <- ggplot(data = dat, aes(x = time2, y = Y, group = id)) +
theme(axis.line = element_line(color = "black"))
if (smooth) {
dat$bs <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
gg <- gg + geom_line(aes(x = time2, y = dat$bs, group = id), color = 4)
if (!is.null(x$Haz0.upper)) {
dat$bs.upper <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y.upper ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
dat$bs.lower <- unlist(lapply(split(dat, dat$id), function(x)
scam(x$Y.lower ~ s(x$time2, k = 10, bs = "mpi"))$fitted.values))
gg <- gg +
geom_line(aes(x = time2, y = dat$bs.upper, group = id), color = 4, lty = 2) +
geom_line(aes(x = time2, y = dat$bs.lower, group = id), color = 4, lty = 2)
}
} else {
gg <- gg + geom_step()
if (!is.null(x$Haz0.upper))
gg <- gg +
geom_step(aes(x = time2, y = Y.upper, group = id), lty = 2) +
geom_step(aes(x = time2, y = Y.lower, group = id), lty = 2)
}
if (showName)
gg <- gg + geom_dl(aes(label = paste(" Obs. =", id)), method = "right.polygons") +
scale_x_continuous(limits = c(0, max(dat$time2) * 1.1))
}
if (ctrl$main != "") gg <- gg + ggtitle(ctrl$main)
gg <- gg + labs(x = ctrl$xlab, y = ctrl$ylab) +
theme(plot.title = element_text(size = 2 * ctrl$base_size),
strip.text = element_text(size = ctrl$base_size),
legend.text = element_text(size = 1.5 * ctrl$base_size),
legend.title = element_text(size = 1.5 * ctrl$base_size),
axis.line = element_blank(),
axis.text = element_text(size = ctrl$base_size),
axis.title = element_text(size = 1.5 * ctrl$base_size))
attr(gg, "from") <- "reReg"
return(gg)
}
#' Plot options for plotEvents
#'
#' This function provides the plotting options for the \code{plotEvents()} function.
#'
#' @param xlab a character string indicating the label for the x axis.
#' The default value is "Time".
#' @param ylab a character string indicating the label for the y axis.
#' The default value is "Subject".
#' @param main a character string indicating the title of the plot.
#' @param terminal.name a character string indicating the label for the terminal event
#' displayed in the legend. The default value is "Terminal event".
#' @param recurrent.name a character string indicating the label for the recurrent event
#' displayed in the legend. The default value is "Recurrent events".
#' @param recurrent.type a factor indicating the labels for the different recurrent event types.
#' This option is only available when there are more than one types of recurrent events.
#' The default value is "Recurrent events 1", "Recurrent events 2", ....
#' @param legend.position a character string specifies the position of the legend.
#' The available options are "none", "left", "right", "bottom", "top",
#' "bottomright", "bottomleft","topleft", "topright",
#' or a two-element numeric vector specifies the coordinate of the legend.
#' The legend is placed inside of the plotting region when
#' \code{legend.position} is specified as one of
#' "bottomright", "bottomleft","topleft", "topright".
#' Otherwise, the argument is passed to the \code{ggplot} theme environment.
#' @param base_size a numerical value to specify the base font size, given in pts.
#' This argument is passed to the \code{ggplot} theme environment.
#' The default value is 12.
#' @param cex a numerical value specifies the size of the points.
#' @param alpha a numerical value specifies the transparency of the points.
#' @param width a numerical value specifies the width of the event plot.
#' By \code{ggplot} default, set to 90\% of the resolution of the data.
#' @param bar.color a numerical value or a character string specifies
#' color for lines. Default to gray.
#' @param recurrent.color a numerical value or a character string
#' specifies color for recurrent events. Default to green.
#' @param terminal.color a numerical value or a character string
#' specifies color for terminal events. Default to red.
#' @param recurrent.shape a numerical value or a character string
#' specifies shape for recurrent events. Default to circle.
#' @param terminal.shape a numerical value or a character string
#' specifies shape for terminal events. Default to triangle.
#' @param recurrent.stroke a numerical value or a character string
#' specifies stroke for recurrent events. Default to circle.
#' @param terminal.stroke a numerical value or a character string
#' specifies stroke for terminal events. Default to triangle.
#' @param not.terminal.color a numerical value or a character string
#' specifies color for non-terminal events.
#' Non-terminal events are not plotted at default.
#' @param not.terminal.shape a numerical value or a character string
#' specifies shape for terminal events.
#' Non-terminal events are not plotted at default.
#'
#' @seealso \code{\link{plotEvents}}
#' @export
plotEvents.control <- function(xlab = NULL, ylab = NULL,
main = NULL,
terminal.name = NULL,
recurrent.name = NULL,
recurrent.type = NULL,
legend.position = NULL,
base_size = 12,
cex = NULL,
alpha = .7,
width = NULL,
bar.color = NULL,
recurrent.color = NULL,
recurrent.shape = NULL,
recurrent.stroke = NULL,
terminal.color = NULL,
terminal.shape = NULL,
terminal.stroke = NULL,
not.terminal.color = NULL,
not.terminal.shape = NULL) {
if (is.null(ylab)) ylab <- "Subject"
if (is.null(xlab)) xlab <- "Time"
if (is.null(main)) main <- ""
if (is.null(legend.position)) legend.position <- "bottomright"
## main <- "Recurrent event plot"
if (is.null(terminal.name)) terminal.name <- "Terminal event"
if (is.null(recurrent.name)) recurrent.name <- "Recurrent events"
if (is.null(bar.color)) bar.color <- "gray75"
if (is.null(terminal.color)) terminal.color <- "red"
if (is.null(terminal.shape)) terminal.shape <- 17
if (is.null(recurrent.shape)) recurrent.shape <- 19
list(xlab = xlab, ylab = ylab, main = main, cex = cex, width = width,
terminal.name = terminal.name, recurrent.name = recurrent.name,
recurrent.type = recurrent.type, alpha = alpha,
legend.position = legend.position, base_size = base_size,
bar.color = bar.color,
recurrent.color = recurrent.color,
recurrent.shape = recurrent.shape,
recurrent.stroke = recurrent.stroke,
not.terminal.color = not.terminal.color,
not.terminal.shape = not.terminal.shape,
terminal.color = terminal.color,
terminal.shape = terminal.shape,
terminal.stroke = terminal.stroke)
}
#' Function used to combine baseline functions in one plot
#'
#' Combine different plots into one.
#'
#' @param ... \code{ggplot} objects created by plotting \code{reReg} objects.
#' @param legend.title an optional character string to specify the legend title.
#' @param legend.labels an optional character string to specify the legend labels.
#' @param control a list of control parameters.
#'
#' @export
#' @keywords Plots
#'
#' @example inst/examples/ex_basebind.R
basebind <- function(..., legend.title, legend.labels, control = list()) {
gglst <- list(...)
if (any(sapply(gglst, function(x) attr(x, "from")) != "reReg"))
stop("Plots must be created from reReg objects")
if (missing(legend.title)) legend.title <- ""
ctrl <- plotEvents.control()
namc <- names(control)
if (!all(namc %in% names(ctrl)))
stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])
ctrl[namc] <- control
if (is.null(ctrl$ylab)) ctrl$ylab <- ""
nargs <- length(gglst)
if (missing(legend.labels)) legend.labels <- 1:nargs
if (length(legend.labels) != nargs) {
message('The length of "name" mismatched, default names are used.')
legend.labels <- 1:nargs
}
d <- do.call(rbind, lapply(gglst, function(x) x$data))
nobs <- sapply(gglst, function(x) nrow(x$data))
d$group <- factor(rep(1:nargs, nobs), labels = legend.labels)
gg <- ggplot(data = d, aes(x = time2, y = Y, color = group)) + geom_step()
if (!is.null(d$Y.upper))
gg <- gg + geom_step(aes(x = time2, y = Y.upper), lty = 2) +
geom_step(aes(x = time2, y = Y.lower), lty = 2)
gg <- gg + labs(x = gglst[[1]]$labels$x, y = gglst[[1]]$labels$y, color = legend.title) +
theme(plot.title = element_text(size = 2 * ctrl$base_size),
strip.text = element_text(size = ctrl$base_size),
legend.text = element_text(size = 1.5 * ctrl$base_size),
legend.title = element_text(size = 1.5 * ctrl$base_size),
axis.line = element_blank(),
axis.text = element_text(size = ctrl$base_size),
axis.title = element_text(size = 1.5 * ctrl$base_size))
## Handle legends
if (ctrl$legend.position %in% c("bottomright", "bottomleft","topleft", "topright")) {
## legend inside
if (ctrl$legend.position == "bottomright") legendAdj <- c(1, 0)
if (ctrl$legend.position == "bottomleft") legendAdj <- c(0, 0)
if (ctrl$legend.position == "topright") legendAdj <- c(1, 1)
if (ctrl$legend.position == "topleft") legendAdj <- c(0, 1)
gg <- gg + theme(
legend.position = legendAdj, legend.justification = legendAdj,
legend.key = element_rect(fill = "transparent", colour = "transparent"),
legend.background = element_rect(fill = alpha(1e-5)))
} else {
## legend outside
gg <- gg + theme(legend.position = ctrl$legend.position,
legend.key = element_rect(fill = "white", color = "white"))
}
gg
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.