# Copyright 2024 Observational Health Data Sciences and Informatics
#
# This file is part of SelfControlledCaseSeries
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#' Plot the age ranges spanned by each observation period.
#'
#' @template StudyPopulation
#' @param maxPersons The maximum number of persons to plot. If there are more than this number of persons
#' a random sample will be taken to avoid visual clutter.
#' @param fileName Name of the file where the plot should be saved, for example 'plot.png'.
#' See the function [ggplot2::ggsave()] for supported file formats.
#' @param title Optional: the main title for the plot
#'
#' @details
#' Plots a line per patient from their age at observation start to their age at observation end.
#'
#' @return
#' A ggplot object. Use the [ggplot2::ggsave()] function to save to file in a different
#' format.
#'
#' @export
plotAgeSpans <- function(studyPopulation,
maxPersons = 10000,
title = NULL,
fileName = NULL) {
errorMessages <- checkmate::makeAssertCollection()
checkmate::assertList(studyPopulation, min.len = 1, add = errorMessages)
checkmate::assertInt(maxPersons, lower = 1, add = errorMessages)
checkmate::assertCharacter(title, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::assertCharacter(fileName, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::reportAssertions(collection = errorMessages)
cases <- studyPopulation$cases %>%
transmute(startAge = .data$ageAtObsStart + .data$startDay,
endAge = .data$ageAtObsStart + .data$endDay) %>%
arrange(.data$startAge, .data$endAge) %>%
mutate(rank = row_number())
ageLabels <- floor(min(cases$startAge) / 365.25):ceiling(max(cases$endAge) / 365.25)
if (length(ageLabels) > 10) {
ageLabels <- 10 * (floor(min(cases$startAge) / 3652.5):floor(max(cases$endAge) / 3652.5))
}
ageBreaks <- ageLabels * 365.25
if (nrow(cases) > maxPersons) {
warning("There are ", nrow(cases), " cases. Random sampling ", maxPersons, " cases.")
cases <- cases[sample.int(nrow(cases), maxPersons), ]
}
theme <- ggplot2::element_text(colour = "#000000", size = 12)
themeRA <- ggplot2::element_text(colour = "#000000", size = 12, hjust = 1)
plot <- ggplot2::ggplot(cases, ggplot2::aes(x = .data$startAge, xmin = .data$startAge, xmax = .data$endAge, y = rank)) +
ggplot2::geom_vline(xintercept = ageBreaks, colour = "#AAAAAA", lty = 1, linewidth = 0.2) +
ggplot2::geom_errorbarh(color = rgb(0, 0, 0.8), alpha = 0.8) +
ggplot2::scale_x_continuous("Age (years)", breaks = ageBreaks, labels = ageLabels) +
ggplot2::scale_y_continuous("Case rank") +
ggplot2::theme(
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_rect(fill = "#FAFAFA", colour = NA),
panel.grid.major = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
axis.text.y = themeRA,
axis.text.x = theme,
strip.text.x = theme,
strip.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5),
legend.title = ggplot2::element_blank(),
legend.position = "top"
)
if (!is.null(title)) {
plot <- plot + ggplot2::ggtitle(title)
}
# fileName <- "S:/temp/plot.png"
if (!is.null(fileName)) {
ggplot2::ggsave(fileName, plot, width = 7, height = 5, dpi = 400)
}
return(plot)
}
computeTimeToObsEnd <- function(studyPopulation) {
outcomes <- studyPopulation$outcomes %>%
group_by(.data$caseId) %>%
summarise(outcomeDay = min(.data$outcomeDay), .groups = "drop_last") %>%
inner_join(studyPopulation$cases, by = join_by("caseId")) %>%
transmute(
daysFromEvent = .data$endDay - .data$outcomeDay,
censoring = case_when(
.data$noninformativeEndCensor == 1 ~ "Uncensored",
TRUE ~ "Censored"
)
)
return(outcomes)
}
#' Plot time from event to observation end for censored and uncensored time.
#'
#' @template StudyPopulation
#' @param fileName Name of the file where the plot should be saved, for example 'plot.png'.
#' See the function [ggplot2::ggsave()] for supported file formats.
#' @param title Optional: the main title for the plot
#'
#' @details
#' This plot shows whether there is a difference in time between (first) event and the observation period end for periods that are '
#' censored' and those that are 'uncensored'. By 'censored' we mean periods that end before we would normally expect. Here, we define
#' periods to be uncensored if they end at either the study end date (if specified), database end date (i.e. the date after which no
#' data is captured in the database), or maximum age (if specified). All other periods are assumed to be censored.
#'
#' As proposed by Farrington et al., by comparing the two plots, we can gain some insight into whether the censoring is dependent
#' on the occurrence of the event.
#'
#' @references
#' Farrington P, Whitaker H, Ghebremichael Weldeselassie Y (2018), Self-controlled case series studies: A modelling guide with R, Taylor &
#' Francis
#'
#' @return
#' A ggplot object. Use the [ggplot2::ggsave()] function to save to file in a different
#' format.
#'
#' @export
plotEventObservationDependence <- function(studyPopulation,
title = NULL,
fileName = NULL) {
errorMessages <- checkmate::makeAssertCollection()
checkmate::assertList(studyPopulation, min.len = 1, add = errorMessages)
checkmate::assertCharacter(title, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::assertCharacter(fileName, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::reportAssertions(collection = errorMessages)
outcomes <- computeTimeToObsEnd(studyPopulation)
ageLabels <- 0:ceiling(max(outcomes$daysFromEvent) / 365.25)
ageBreaks <- ageLabels * 365.25
theme <- ggplot2::element_text(colour = "#000000", size = 12)
themeRA <- ggplot2::element_text(colour = "#000000", size = 12, hjust = 1)
plot <- ggplot2::ggplot(outcomes, ggplot2::aes(x = .data$daysFromEvent)) +
ggplot2::geom_vline(xintercept = ageBreaks, colour = "#AAAAAA", lty = 1, linewidth = 0.2) +
ggplot2::geom_histogram(binwidth = 30.5, fill = rgb(0, 0, 0.8), alpha = 0.8) +
ggplot2::scale_x_continuous("Years from event", breaks = ageBreaks, labels = ageLabels) +
ggplot2::scale_y_continuous("Frequency") +
ggplot2::facet_grid(censoring ~ ., scales = "free_y") +
ggplot2::theme(
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_rect(fill = "#FAFAFA", colour = NA),
panel.grid.major = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
axis.text.y = themeRA,
axis.text.x = theme,
strip.text.y = theme,
strip.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5),
legend.title = ggplot2::element_blank(),
legend.position = "top"
)
if (!is.null(title)) {
plot <- plot + ggplot2::ggtitle(title)
}
if (!is.null(fileName)) {
ggplot2::ggsave(fileName, plot, width = 7, height = 5, dpi = 400)
}
return(plot)
}
computeTimeToEvent <- function(studyPopulation,
sccsData,
exposureEraId) {
cases <- studyPopulation$cases %>%
select("caseId", "startDay", "endDay")
exposures <- sccsData$eras %>%
filter(.data$eraId == exposureEraId & .data$eraType == "rx") %>%
group_by(.data$caseId) %>%
inner_join(cases,
by = join_by("caseId", "eraStartDay" >= "startDay", "eraStartDay" < "endDay"),
copy = TRUE) %>%
collect()
if (nrow(exposures) == 0) {
warning("No exposures found with era ID ", exposureEraId)
tibble(
number = 1,
start = 1.0,
end = 1.0,
observed = 1,
eventsExposed = 1.0,
eventsUnexposed = 1.0
) %>%
filter(.data$number == -1) %>%
return()
}
firstExposures <- exposures %>%
group_by(.data$caseId, .data$startDay, .data$endDay) %>%
summarise(
eraStartDay = min(.data$eraStartDay, na.rm = TRUE),
eraEndDay = min(.data$eraEndDay, na.rm = TRUE),
.groups = "drop_last"
)
outcomes <- studyPopulation$outcomes %>%
inner_join(firstExposures, by = join_by("caseId")) %>%
mutate(delta = .data$outcomeDay - .data$eraStartDay) %>%
select("caseId", "outcomeDay", "delta")
exposedoutcomes <- exposures %>%
inner_join(outcomes, by = join_by("caseId"), relationship = "many-to-many") %>%
filter(
.data$outcomeDay >= .data$eraStartDay,
.data$outcomeDay <= .data$eraEndDay
) %>%
select("caseId", "delta") %>%
mutate(exposed = 1)
outcomes <- outcomes %>%
left_join(exposedoutcomes, by = join_by("caseId", "delta"), relationship = "many-to-many") %>%
mutate(exposed = coalesce(.data$exposed, 0))
weeks <- dplyr::tibble(number = -26:25) %>%
mutate(
start = .data$number * 7,
end = .data$number * 7 + 7
)
events <- weeks %>%
cross_join(select(outcomes, "delta", "exposed")) %>%
filter(.data$delta >= .data$start, .data$delta < .data$end) %>%
group_by(.data$number, .data$start, .data$end) %>%
summarise(
eventsExposed = sum(.data$exposed),
eventsUnexposed = n() - sum(.data$exposed),
.groups = "drop"
)
observed <- weeks %>%
cross_join(transmute(firstExposures, startDelta = .data$startDay -.data$eraStartDay , endDelta = .data$endDay - .data$eraStartDay)) %>%
filter(.data$endDelta >= .data$start, .data$startDelta < .data$end) %>%
group_by(.data$number, .data$start, .data$end) %>%
summarise(
observed = n(),
.groups = "drop"
)
result <- observed %>%
left_join(events, by = join_by("number", "start", "end")) %>%
mutate(
eventsExposed = if_else(is.na(.data$eventsExposed), 0, .data$eventsExposed),
eventsUnexposed = if_else(is.na(.data$eventsUnexposed), 0, .data$eventsUnexposed)
)
return(result)
}
#' Plot information centered around the start of exposure
#'
#' @param exposureEraId The exposure to create the era data for. If not specified it is
#' assumed to be the one exposure for which the data was loaded from
#' the database.
#' @template StudyPopulation
#' @template SccsData
#' @param highlightExposedEvents Highlight events that occurred during the exposure era using a different color?
#' @param fileName Name of the file where the plot should be saved, for example 'plot.png'.
#' See the function [ggplot2::ggsave()] for supported file formats.
#' @param title Optional: the main title for the plot
#'
#' @details
#' This plot shows the number of events and the number of subjects under observation in week-sized intervals relative to the start
#' of the first exposure.
#'
#' @return
#' A ggplot object. Use the [ggplot2::ggsave()] function to save to file in a different
#' format.
#'
#' @export
plotExposureCentered <- function(studyPopulation,
sccsData,
exposureEraId = NULL,
highlightExposedEvents = TRUE,
title = NULL,
fileName = NULL) {
errorMessages <- checkmate::makeAssertCollection()
checkmate::assertList(studyPopulation, min.len = 1, add = errorMessages)
checkmate::assertClass(sccsData, "SccsData", add = errorMessages)
checkmate::assertInt(exposureEraId, null.ok = TRUE, add = errorMessages)
checkmate::assertLogical(highlightExposedEvents, len = 1, add = errorMessages)
checkmate::assertCharacter(title, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::assertCharacter(fileName, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::reportAssertions(collection = errorMessages)
if (is.null(exposureEraId)) {
exposureEraId <- attr(sccsData, "metaData")$exposureIds
if (length(exposureEraId) != 1) {
stop("No exposure ID specified, but multiple exposures found")
}
}
data <- computeTimeToEvent(studyPopulation, sccsData, exposureEraId)
if (nrow(data) == 0) {
return(NULL)
}
if (highlightExposedEvents) {
events <- data %>%
transmute(.data$start,
.data$end,
type = "Events",
count1 = .data$eventsUnexposed,
count2 = .data$eventsExposed
)
} else {
events <- data %>%
transmute(.data$start,
.data$end,
type = "Events",
count1 = .data$eventsUnexposed + .data$eventsExposed,
count2 = NA
)
}
observed <- data %>%
transmute(.data$start,
.data$end,
type = "Subjects under observation",
count1 = .data$observed,
count2 = NA
)
data <- bind_rows(events, observed)
breaks <- seq(-150, 150, 30)
theme <- ggplot2::element_text(colour = "#000000", size = 12)
themeRA <- ggplot2::element_text(colour = "#000000", size = 12, hjust = 1)
plot <- ggplot2::ggplot(data, ggplot2::aes(x = .data$start, xmin = .data$start, xmax = .data$end, ymax = .data$count1, ymin = 0)) +
ggplot2::geom_rect(fill = rgb(0, 0, 0.8), alpha = 0.8) +
ggplot2::geom_rect(ggplot2::aes(ymax = .data$count1 + .data$count2, ymin = .data$count1), fill = rgb(0.8, 0, 0), alpha = 0.8) +
ggplot2::geom_vline(xintercept = 0, colour = "#000000", lty = 1, linewidth = 1) +
ggplot2::scale_x_continuous("Days since first exposure start", breaks = breaks, labels = breaks) +
ggplot2::scale_y_continuous("Count") +
ggplot2::facet_grid(type ~ ., scales = "free_y") +
ggplot2::theme(
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_rect(fill = "#FAFAFA", colour = NA),
panel.grid.major = ggplot2::element_line(colour = "#AAAAAA"),
axis.ticks = ggplot2::element_blank(),
axis.text.y = themeRA,
axis.text.x = theme,
strip.text.y = theme,
strip.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5),
legend.title = ggplot2::element_blank(),
legend.position = "top"
)
if (!is.null(title)) {
plot <- plot + ggplot2::ggtitle(title)
}
if (!is.null(fileName)) {
ggplot2::ggsave(fileName, plot, width = 7, height = 5, dpi = 400)
}
return(plot)
}
#' Plot the ratio of observed to expected events over calendar time.
#'
#' @details
#' Plot the ratio of observed to expected events over calendar time. The expected count expected
#' rate considers which persons were observed during that month, and if specified in the model,
#' the adjustment for season and calendar time.
#'
#' @template StudyPopulation
#' @param sccsModel Optional: A fitted SCCS model as created using [fitSccsModel()]. If the
#' model contains splines for seasonality and or calendar time a panel will
#' be added with outcome counts adjusted for these splines.
#' @param fileName Name of the file where the plot should be saved, for example 'plot.png'.
#' See the function [ggplot2::ggsave()] for supported file formats.
#' @param title Optional: the main title for the plot
#'
#' @return
#' A ggplot object. Use the [ggplot2::ggsave()] function to save to file in a different
#' format.
#'
#' @export
plotEventToCalendarTime <- function(studyPopulation,
sccsModel = NULL,
title = NULL,
fileName = NULL) {
errorMessages <- checkmate::makeAssertCollection()
checkmate::assertList(studyPopulation, min.len = 1, add = errorMessages)
checkmate::assertClass(sccsModel, "SccsModel", null.ok = TRUE, add = errorMessages)
checkmate::assertCharacter(title, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::assertCharacter(fileName, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::reportAssertions(collection = errorMessages)
data <- computeOutcomeRatePerMonth(studyPopulation, sccsModel)
plotData <- data %>%
select("month", "monthStartDate", "monthEndDate", value = "ratio") %>%
mutate(type = "Assuming constant rate")
levels <- c("Assuming constant rate")
if (!is.null(sccsModel) && (hasCalendarTimeEffect(sccsModel) || hasSeasonality(sccsModel))) {
types <- c("cal. time", "season")[c(hasCalendarTimeEffect(sccsModel), hasSeasonality(sccsModel))]
type <- paste("Adj. for", paste(types, collapse = " and "))
plotData <- bind_rows(
plotData,
select(data, "month", "monthStartDate", "monthEndDate", value = "adjustedRatio") %>%
mutate(type = !!type),
)
levels <- c(levels, type)
}
plotData$type <- factor(plotData$type, levels = rev(levels))
theme <- ggplot2::element_text(colour = "#000000", size = 12)
themeRA <- ggplot2::element_text(colour = "#000000", size = 12, hjust = 1)
plot <- ggplot2::ggplot(plotData, ggplot2::aes(xmin = .data$monthStartDate, xmax = .data$monthEndDate + 1)) +
ggplot2::geom_rect(ggplot2::aes(ymax = .data$value), ymin = 0, fill = rgb(0, 0, 0.8), alpha = 0.8, linewidth = 0) +
ggplot2::scale_x_date("Calendar time") +
ggplot2::scale_y_continuous("Observed / expected", limits = c(0, NA)) +
ggplot2::facet_grid(.data$type ~ ., scales = "free_y") +
ggplot2::theme(
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_rect(fill = "#FAFAFA", colour = NA),
panel.grid.major = ggplot2::element_line(colour = "#AAAAAA"),
axis.ticks = ggplot2::element_blank(),
axis.text.y = themeRA,
axis.text.x = theme,
strip.text.y = theme,
strip.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5),
legend.title = ggplot2::element_blank(),
legend.position = "top"
)
# plot
if (!is.null(title)) {
plot <- plot + ggplot2::ggtitle(title)
}
if (!is.null(fileName)) {
ggplot2::ggsave(fileName, plot, width = 7, height = 1 + (2 * length(levels)), dpi = 400)
}
return(plot)
}
#' Plot the age effect
#'
#' @details
#' Plot the spline curve of the age effect.
#'
#' @param sccsModel An object of type \code{sccsModel} as created using the
#' \code{\link{fitSccsModel}} function.
#' @param rrLim The limits on the incidence rate ratio scale in the plot.
#' @param fileName Name of the file where the plot should be saved, for example 'plot.png'. See the
#' function \code{ggsave} in the ggplot2 package for supported file formats.
#' @param title Optional: the main title for the plot
#'
#' @return
#' A Ggplot object. Use the ggsave function to save to file.
#'
#' @export
plotAgeEffect <- function(sccsModel,
rrLim = c(0.1, 10),
title = NULL,
fileName = NULL) {
errorMessages <- checkmate::makeAssertCollection()
checkmate::assertClass(sccsModel, "SccsModel", add = errorMessages)
checkmate::assertNumeric(rrLim, len = 2, add = errorMessages)
checkmate::assertCharacter(title, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::assertCharacter(fileName, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::reportAssertions(collection = errorMessages)
if (!hasAgeEffect(sccsModel)) {
stop("The model does not contain an age effect.")
}
estimates <- sccsModel$estimates
splineCoefs <- estimates[estimates$covariateId >= 100 & estimates$covariateId < 200, "logRr"]
ageKnots <- sccsModel$metaData$age$ageKnots
age <- seq(min(ageKnots), max(ageKnots), length.out = 100)
ageDesignMatrix <- splines::bs(age,
knots = ageKnots[2:(length(ageKnots) - 1)],
Boundary.knots = ageKnots[c(1, length(ageKnots))],
degree = 2
)
logRr <- apply(ageDesignMatrix %*% splineCoefs, 1, sum)
logRr <- logRr - mean(logRr)
rr <- exp(logRr)
data <- data.frame(age = age, rr = rr)
breaks <- c(0.1, 0.25, 0.5, 1, 2, 4, 6, 8, 10)
ageLabels <- floor(min(ageKnots) / 365.25):floor(max(ageKnots) / 365.25)
if (length(ageLabels) > 10) {
ageLabels <- 10 * (floor(min(ageKnots) / 3652.5):floor(max(ageKnots) / 3652.5))
}
ageBreaks <- ageLabels * 365.25
theme <- ggplot2::element_text(colour = "#000000", size = 12)
themeRA <- ggplot2::element_text(colour = "#000000", size = 12, hjust = 1)
plot <- ggplot2::ggplot(data, ggplot2::aes(x = .data$age, y = .data$rr)) +
ggplot2::geom_hline(yintercept = breaks, colour = "#AAAAAA", lty = 1, linewidth = 0.2) +
ggplot2::geom_line(color = rgb(0, 0, 0.8), alpha = 0.8, linewidth = 1) +
ggplot2::scale_x_continuous("Age", breaks = ageBreaks, labels = ageLabels) +
ggplot2::scale_y_continuous("Relative risk",
limits = rrLim,
trans = "log10",
breaks = breaks,
labels = breaks
) +
ggplot2::theme(
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_rect(fill = "#FAFAFA", colour = NA),
panel.grid.major = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
axis.text.y = themeRA,
axis.text.x = theme,
strip.text.x = theme,
strip.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5),
legend.title = ggplot2::element_blank(),
legend.position = "top"
)
if (!is.null(title)) {
plot <- plot + ggplot2::ggtitle(title)
}
if (!is.null(fileName)) {
ggplot2::ggsave(fileName, plot, width = 7, height = 5, dpi = 400)
}
return(plot)
}
#' Plot the seasonality effect
#'
#' @details
#' Plot the spline curve of the seasonality effect.
#'
#' @param sccsModel An object of type \code{sccsModel} as created using the
#' \code{\link{fitSccsModel}} function.
#' @param rrLim The limits on the incidence rate ratio scale in the plot.
#' @param fileName Name of the file where the plot should be saved, for example 'plot.png'. See the
#' function \code{ggsave} in the ggplot2 package for supported file formats.
#' @param title Optional: the main title for the plot
#'
#' @return
#' A Ggplot object. Use the ggsave function to save to file.
#'
#' @export
plotSeasonality <- function(sccsModel,
rrLim = c(0.1, 10),
title = NULL,
fileName = NULL) {
errorMessages <- checkmate::makeAssertCollection()
checkmate::assertClass(sccsModel, "SccsModel", add = errorMessages)
checkmate::assertNumeric(rrLim, len = 2, add = errorMessages)
checkmate::assertCharacter(title, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::assertCharacter(fileName, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::reportAssertions(collection = errorMessages)
if (!hasSeasonality(sccsModel)) {
stop("The model does not contain seasonality.")
}
estimates <- sccsModel$estimates
splineCoefs <- estimates[estimates$covariateId >= 200 & estimates$covariateId < 300, "logRr"]
seasonKnots <- sccsModel$metaData$seasonality$seasonKnots
season <- sort(unique(c(seq(min(seasonKnots), max(seasonKnots), length.out = 100),
seasonKnots)))
seasonDesignMatrix <- cyclicSplineDesign(season, seasonKnots)
logRr <- apply(seasonDesignMatrix %*% splineCoefs, 1, sum)
logRr <- logRr - mean(logRr)
rr <- exp(logRr)
data <- tibble(season = season, rr = rr)
knotData <- data %>%
filter(.data$season %in% seasonKnots)
breaks <- c(0.1, 0.25, 0.5, 1, 2, 4, 6, 8, 10)
seasonBreaks <- 1:12
theme <- ggplot2::element_text(colour = "#000000", size = 12)
themeRA <- ggplot2::element_text(colour = "#000000", size = 12, hjust = 1)
plot <- ggplot2::ggplot(data, ggplot2::aes(x = season, y = rr)) +
ggplot2::geom_hline(yintercept = breaks, colour = "#AAAAAA", lty = 1, linewidth = 0.2) +
ggplot2::geom_line(color = rgb(0, 0, 0.8), alpha = 0.8, linewidth = 1) +
ggplot2::geom_point(data = knotData) +
ggplot2::scale_x_continuous("Month", breaks = seasonBreaks, labels = seasonBreaks) +
ggplot2::scale_y_continuous("Relative risk",
limits = rrLim,
trans = "log10",
breaks = breaks,
labels = breaks
) +
ggplot2::theme(
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_rect(fill = "#FAFAFA", colour = NA),
panel.grid.major = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
axis.text.y = themeRA,
axis.text.x = theme,
strip.text.x = theme,
strip.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5),
legend.title = ggplot2::element_blank(),
legend.position = "top"
)
if (!is.null(title)) {
plot <- plot + ggplot2::ggtitle(title)
}
if (!is.null(fileName)) {
ggplot2::ggsave(fileName, plot, width = 7, height = 5, dpi = 400)
}
return(plot)
}
#' Plot the calendar time ranges spanned by each observation period.
#'
#' @template StudyPopulation
#' @param maxPersons The maximum number of persons to plot. If there are more than this number of persons
#' a random sample will be taken to avoid visual clutter.
#' @param fileName Name of the file where the plot should be saved, for example 'plot.png'.
#' See the function [ggplot2::ggsave()] for supported file formats.
#' @param title Optional: the main title for the plot
#'
#' @details
#' Plots a line per patient from their observation start to their observation end.
#'
#' @return
#' A ggplot object. Use the [ggplot2::ggsave()] function to save to file in a different
#' format.
#'
#' @export
plotCalendarTimeSpans <- function(studyPopulation,
maxPersons = 10000,
title = NULL,
fileName = NULL) {
errorMessages <- checkmate::makeAssertCollection()
checkmate::assertList(studyPopulation, add = errorMessages)
checkmate::assertInt(maxPersons, lower = 1, add = errorMessages)
checkmate::assertCharacter(title, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::assertCharacter(fileName, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::reportAssertions(collection = errorMessages)
cases <- studyPopulation$cases %>%
mutate(startDate = .data$observationPeriodStartDate + .data$startDay,
endDate = .data$observationPeriodStartDate + .data$endDay) %>%
select("startDate", "endDate") %>%
arrange(.data$startDate, .data$endDate) %>%
mutate(rank = row_number())
if (nrow(cases) > maxPersons) {
warning("There are ", nrow(cases), " cases. Random sampling ", maxPersons, " cases.")
cases <- cases[sample.int(nrow(cases), maxPersons), ]
}
theme <- ggplot2::element_text(colour = "#000000", size = 12)
themeRA <- ggplot2::element_text(colour = "#000000", size = 12, hjust = 1)
plot <- ggplot2::ggplot(cases, ggplot2::aes(x = .data$startDate, xmin = .data$startDate, xmax = .data$endDate, y = rank)) +
ggplot2::geom_errorbarh(color = rgb(0, 0, 0.8)) +
ggplot2::scale_x_date("Calendar time") +
ggplot2::scale_y_continuous("Case rank") +
ggplot2::theme(
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_rect(fill = "#FAFAFA", colour = NA),
panel.grid.major.x = ggplot2::element_line(colour = "#AAAAAA", linewidth = 0.2),
panel.grid.major.y = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
axis.text.y = themeRA,
axis.text.x = theme,
strip.text.x = theme,
strip.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5),
legend.title = ggplot2::element_blank(),
legend.position = "top"
)
if (!is.null(title)) {
plot <- plot + ggplot2::ggtitle(title)
}
# fileName <- "S:/temp/plot.png"
if (!is.null(fileName)) {
ggplot2::ggsave(fileName, plot, width = 7, height = 5, dpi = 400)
}
return(plot)
}
#' Plot the calendar time effect
#'
#' @details
#' Plot the spline curve of the calendar time effect.
#'
#' @param sccsModel An object of type \code{sccsModel} as created using the
#' \code{\link{fitSccsModel}} function.
#' @param rrLim The limits on the incidence rate ratio scale in the plot.
#' @param fileName Name of the file where the plot should be saved, for example 'plot.png'. See the
#' function \code{ggsave} in the ggplot2 package for supported file formats.
#' @param title Optional: the main title for the plot
#'
#' @return
#' A Ggplot object. Use the ggsave function to save to file.
#'
#' @export
plotCalendarTimeEffect <- function(sccsModel,
rrLim = c(0.1, 10),
title = NULL,
fileName = NULL) {
errorMessages <- checkmate::makeAssertCollection()
checkmate::assertClass(sccsModel, "SccsModel", add = errorMessages)
checkmate::assertNumeric(rrLim, len = 2, add = errorMessages)
checkmate::assertCharacter(title, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::assertCharacter(fileName, len = 1, null.ok = TRUE, add = errorMessages)
checkmate::reportAssertions(collection = errorMessages)
if (!hasCalendarTimeEffect(sccsModel)) {
stop("The model does not contain a calendar time effect.")
}
estimates <- sccsModel$estimates
splineCoefs <- estimates[estimates$covariateId >= 300 & estimates$covariateId < 400, "logRr"]
calendarTimeKnotsInPeriods <- sccsModel$metaData$calendarTime$calendarTimeKnotsInPeriods
allKnots <- do.call(c, calendarTimeKnotsInPeriods)
x <- sort(unique(c(seq(min(allKnots), max(allKnots), length.out = 200), allKnots)))
designMatrix <- createMultiSegmentDesignMatrix(x = x, knotsPerSegment = calendarTimeKnotsInPeriods)
logRr <- apply(designMatrix %*% splineCoefs, 1, sum)
# logRr <- logRr - mean(logRr)
rr <- exp(logRr)
data <- tibble(date = convertMonthToStartDate(x) + 14,
calendarTime = x,
rr = rr,
segment = 0)
for (i in seq_along(calendarTimeKnotsInPeriods)) {
knots <- calendarTimeKnotsInPeriods[[i]]
data <- data %>%
mutate(segment = if_else(.data$calendarTime >= min(knots) & .data$calendarTime <= max(knots),
i,
.data$segment))
}
data <- data %>%
filter(.data$segment > 0)
knotData <- data %>%
filter(.data$calendarTime %in% allKnots) %>%
arrange(.data$calendarTime)
breaks <- c(0.1, 0.25, 0.5, 1, 2, 4, 6, 8, 10)
theme <- ggplot2::element_text(colour = "#000000", size = 12)
themeRA <- ggplot2::element_text(colour = "#000000", size = 12, hjust = 1)
plot <- ggplot2::ggplot(data, ggplot2::aes(x = date, y = rr)) +
ggplot2::geom_line(ggplot2::aes(group = .data$segment), color = rgb(0, 0, 0.8), alpha = 0.8, linewidth = 1) +
ggplot2::geom_point(data = knotData) +
ggplot2::scale_x_date("Calendar Time") +
ggplot2::scale_y_continuous("Relative risk",
limits = rrLim,
trans = "log10",
breaks = breaks,
labels = breaks
) +
ggplot2::theme(
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_rect(fill = "#FAFAFA", colour = NA),
panel.grid.major = ggplot2::element_line(colour = "#AAAAAA", linewidth = 0.2),
axis.ticks = ggplot2::element_blank(),
axis.text.y = themeRA,
axis.text.x = theme,
strip.text.x = theme,
strip.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5),
legend.title = ggplot2::element_blank(),
legend.position = "top"
)
# plot
if (!is.null(title)) {
plot <- plot + ggplot2::ggtitle(title)
}
if (!is.null(fileName)) {
ggplot2::ggsave(fileName, plot, width = 7, height = 5, dpi = 400)
}
return(plot)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.