R/cosinor-plots.R

In card: Cardiovascular Applications in Research Data

#' @title Graphical Assessment of Amplitude and Acrophase
#' @description This is a ggplot-styled graphical representation of the ellipse
#'   region generated by the cosinor analysis. It requires the same data used by
#'   cosinor model to be fit with the model [card::cosinor]. This includes
#'   the amplitude, acrophase,
#' @param object Requires a cosinor model to extract the correct statistics to
#'   generate the plot.
#' @param level Confidence level for ellipse
#' @param ... Additional parameters may be needed for extensibility
#' @examples
#' data("twins")
#' m <- cosinor(rDYX ~ hour, twins, tau = 24)
#' ggellipse(m)
#' @return Object of class `ggplot` to help identify confidence intervals
#' @import ggplot2
#' @export
ggellipse <- function(object, level = 0.95, ...) {

	if (object$type == "Population") {
		message("Ellipse may not be accurate for population-mean cosinor method.")
	}

	# Area
	area <- cosinor_area(object, level = level)$area
  gseq <- area[, "gseq"]
  bs1 <- area[, "bs1"]
  bs2 <- area[, "bs2"]

	# Confidence level
	a <- 1 - level

	# Parameters
	y <- object$model[,"y"]
	t <- object$model[,"t"]
	n <- length(t)
	p <- length(object$tau)

  # Create null variables
  mesor <- NULL
  for (i in 1:p) {
    assign(paste0("x", i), NULL)
    assign(paste0("z", i), NULL)
    assign(paste("amp", i), NULL)
    assign(paste("phi", i), NULL)
    assign(paste("beta", i), NULL)
    assign(paste("gamma", i), NULL)
  }
  x1 <- z1 <- beta1 <- gamma1 <- amp1 <- phi1 <- NULL

  for (i in 1:p) {
    assign(paste0("x", i), object$model[, paste0("x", i)])
    assign(paste0("z", i), object$model[, paste0("z", i)])
  }

	xmat <- object$xmat
	yhat <- object$fitted.values
	coefs <- object$coefficients
	names(coefs) <- object$coef_names
  for (i in 1:length(coefs)) {
    assign(names(coefs)[i], unname(coefs[i]))
	}

  # Necessary values for the plot
  theta_clock <- seq(0, 2 * pi, length.out = 24^2)
  clock <- cbind(2 * amp1 * cos(theta_clock), 2 * amp1 * sin(theta_clock))
  rad <- seq(0, 2 * pi - pi / 4, by = pi / 4)
  rad_clock <- cbind(2.2 * amp1 * cos(rad), 2.2 * amp1 * sin(rad))

  # GGplot the values
  ggplot() +
  	# Ellipse
  	geom_line(aes(x = gseq, y = bs1), col = "goldenrod", size = 1) +
  	geom_line(aes(x = gseq, y = bs2), col = "goldenrod", size = 1) +
  	# Line from origin to ellipse
  	geom_line(aes(x = c(0, gamma1),
  								y = c(0, beta1)),
  						lty = 1,
  						size = 1,
  						col = "black") +
  	# Line from ellipse to circumference
  	geom_line(aes(
  		x = c(gamma1, -2 * amp1 * sin(phi1)),
  		y = c(beta1, 2 * amp1 * cos(phi1)),
  		group = 0
  	),
  	size = 1,
  	col = "black",
  	lty = 3) +
  	# Axes
  	geom_line(aes(x = c(0, 0), y = c(-2 * amp1, 2 * amp1)),
  						lty = 5, col = "grey") +
  	geom_line(aes(y = c(0, 0), x = c(-2 * amp1, 2 * amp1)),
  						lty = 5, col = "grey") +
  	# "Clock" shape to help show degrees on unit circle
  	geom_path(aes(x = clock[, 1], y = clock[, 2]), col = "cornflowerblue") +
  	annotate(
  		geom = "text",
  		x = rad_clock[, 1],
  		y = rad_clock[, 2],
  		label = paste(rad * 180 / pi, "\u00B0")
  	) +
  	# Labels
  	xlab(expression(paste(gamma1))) +
  	ylab(expression(paste(beta1))) +
  	xlim(-2.5 * amp1, 2.5 * amp1) +
  	ylim(-2.5 * amp1, 2.5 * amp1) +
  	theme_minimal()
}

#' @title ggplot of cosinor model
#' @description ggplot of cosinor model that can visualize a variety of cosinor
#'   model subtypes, including single-component, multiple-component, individual,
#'   and population cosinor models, built using [card::cosinor]. For single
#'   component cosinor, the following values are plotted:
#'
#'   * M = midline estimating statistic of rhythm
#'
#'   * A = amplitude
#'
#'   * P = phi or acrophase (shift from 0 to peak)
#'
#'   If using a multiple-component cosinor, the terms are different. If the
#'   periods or frequencies resonate or are harmonic, then the following are
#'   calculated. If the periods are not harmonic, the values are just
#'   descriptors of the curve.
#'
#'   * M = midline estimating statistic of rhythm
#'
#'   * Ag = global amplitude, which is the distance between peak and trough
#'   (this is the same value as the amplitude from single component)
#'
#'   * Po = orthophase (the equivalent of the acrophase in a single component),
#'   the lag time to peak value
#'
#'   * Pb = bathyphase, the lag time to trough value
#'
#' @param object Model of class `cosinor`. If instead of a single cosinor model,
#'   multiple objects are to be plotted, can provide a list of cosinor models.
#'   Plotting multiple models simultaneously is preferred if the outcome
#'   variable is similar in scale.
#' @param labels Logical value if annotations should be placed on plot, default
#'   = TRUE. The labels depend on the type of plot. The labels are attempted to
#'   be placed "smartly" using the [ggrepel::geom_label_repel()] function.
#' @param ... For extensibility. This function will use different
#'   implementations based on the type of model (single or multiple component).
#'   Attributes of the object will be passed down, or calculated on the fly.
#' @return Object of class `ggplot` that can be layered
#' @examples
#' data(triplets)
#' m1 <- cosinor(rDYX ~ hour, twins, tau = 24)
#' m2 <- cosinor(rDYX ~ hour, twins, tau = c(24, 12))
#' ggcosinor(m1, labels = FALSE)
#' ggcosinor(m2)
#' ggcosinor(list(single = m1, multiple = m2))
#' @import ggplot2
#' @family cosinor
#' @export
ggcosinor <- function(object, labels = TRUE, ...) {

	# Check to make sure just single object
	if ("cosinor" %in% class(object)) {
		# Model basic data
		tau <- object$tau
		p <- length(tau) # Single or multicomponent

	  # Create null variables
	  mesor <- NULL
	  for (i in 1:p) {
	    assign(paste0("x", i), NULL)
	    assign(paste0("z", i), NULL)
	    assign(paste("amp", i), NULL)
	    assign(paste("phi", i), NULL)
	    assign(paste("beta", i), NULL)
	    assign(paste("gamma", i), NULL)
	  }
	  .fitted <- .resid <- x <- y <- term <- value <- NULL
	  orthophase <- bathyphase <- trough <- ampGlobal <- NULL

		type <- dplyr::case_when(
			p == 1 & object$type == "Individual" ~ "scomp",
			p > 1 & object$type == "Individual" ~ "mcomp",
		)
			if (p == 1 & object$type == "Individual") {
				"single"
			} else if (p > 1 & object$type == "Individual") {
				"multiple"
			} else if (object$type == "Population") {
				stop("`ggcosinor` does not currently support plotting population cosinor models.", call. = FALSE)
			}
	} else if (is.vector(object) & "cosinor" %in% class(object[[1]])) {
		gg <- ggmulticosinor(object, labels)
		return(gg)
	} else {
		stop("Cannot determine if model is cosinor object.", call. = FALSE)
	}

	# Identify which function to call
	aug <- augment(object)

	# Coefficients (include amplitudes)
	coefs <- object$coefficients
	names(coefs) <- object$coef_names
  for (i in 1:length(coefs)) {
    assign(names(coefs)[i], unname(coefs[i]))
  }

	# Acrophases
	for (i in 1:p) {
		assign(paste0("acro", i), abs((get(paste0("phi", i)) * tau[i]) / (2 * pi)))
	}

	## Common geom mappings
	# Amplitude
	amp <- list()
	for (i in 1:p) {
		amp[[i]] <-
			geom_segment(
			aes(
				x = !! get(paste0("acro", i)),
				xend = !! get(paste0("acro", i)),
				y = !! mesor,
				yend = !! (mesor + get(paste0("amp", i)))
			),
			linetype = "twodash", lineend = "butt", linejoin = "mitre",
		)
	}

	# Acrophases
	acro <- list()
	for (i in 1:p) {
		acro[[i]] <-
			geom_segment(
			aes(
				x = 0,
				xend = !! get(paste0("acro", i)),
				y = !! (mesor + get(paste0("amp", i))),
				yend = !! (mesor + get(paste0("amp", i)))
			),
			linetype = "twodash", lineend = "butt", linejoin = "mitre",
		)
	}

	## Data points / labels

	# Features that may be needed
	features <- cosinor_features(object)
	orthophase <- features$orthophase
	bathyphase <- features$bathyphase
	peak <- features$peak
	trough <- features$trough
	zero <- min(aug$t)
	mid <- mean(aug$t)

	glabs <-
		dplyr::bind_cols(
		term = c("M", paste0("A", 1:p), paste0("P", 1:p), "Po", "Pb"),
		value = c(
			mesor,
			unlist(mget(paste0("amp", 1:p))),
			unlist(mget(paste0("acro", 1:p))),
			orthophase,
			bathyphase
		),
		x = c(
			zero,
			unlist(mget(paste0("acro", 1:p))),
			(zero + unlist(mget(paste0("acro", 1:p))))/2,
			orthophase,
			bathyphase
		),
		y = c(
			mesor,
			mesor + unlist(mget(paste0("amp", 1:p)))/2,
			unlist(mget(paste0("amp", 1:p))) + mesor,
			peak,
			trough
		)
	)

	# Basic plot
	g <- ggplot() +
		stat_smooth(
			aes(x = t, y = .fitted), data = aug,
			method = "gam", color = "black", size = 1.2
		) +
		geom_hline(
			yintercept = mesor,
			color = "grey"
		) +
		geom_vline(
			xintercept = orthophase,
			color = "grey"
		)


	# Eventually add in residuals
	gres <- list(
		geom_segment(
			aes(x = t, xend = t, y = y, yend = .fitted), data = aug,
			alpha = 0.3
		),
		geom_point(
			aes(x = t, y = y, colour = abs(.resid), size = abs(.resid)), data = aug
		),
		scale_color_viridis_c(option = "magma")
	)

	# Switch to correct function
	switch(
		type,
		scomp = {

			# For single component, just peak and trough values
			glabs$term[glabs$term == "Po"] <- "Peak"
			glabs$term[glabs$term == "Pb"] <- "Trough"
			glabs$value[glabs$term == "Peak"] <- glabs$y[glabs$term == "Peak"]
			glabs$value[glabs$term == "Trough"] <- glabs$y[glabs$term == "Trough"]

			# Labels if needed
			gl <- unlist(list(

				# Amplitude and Acrophase
				amp,
				acro,

				# Peak and trough
				geom_point(
					aes(x = orthophase, y = peak),
					shape = 18, size = 5
				),
				geom_point(
					aes(x = bathyphase, y = trough),
					shape = 18, size = 5
				),

				# All labeling points from "glabs" created above
				geom_point(aes(x = x, y = y), glabs, alpha = 0),

				# Repelling labels
				ggrepel::geom_label_repel(
					aes(
						x = x, y = y, label = paste0(term, " = ", round(value, 3))
					),
					data = glabs,
					label.size = NA, label.r = 0.25, label.padding = 0.25,
					force = 20,
					nudge_y =
						ifelse(
							glabs$term == "P1", 0.01 * mesor,
							ifelse(glabs$term %in% c("M", "Trough"), -0.01 * mesor, 0)
						),
					nudge_x =
						ifelse(
							glabs$term %in% c("M", "A1"), 0.10 * mid,
							ifelse(glabs$term == "Peak", 0.20 * mid, 0)
						),
					segment.color = "transparent"
				)
			))

			if(labels) {g <- g + gl} else {g <- g}
		},

		mcomp = {

			gl <- list(

				# Mesor
				geom_text(
					aes(x = 2*zero, y = 1.01*mesor),
					label = paste0("M = ", round(mesor, 3))
				),

				# Orthophase
				geom_segment(
					aes(x = zero, xend = orthophase, y = peak, yend = peak),
					linetype = "dotted", size = 0.5
				),
				geom_text(
					aes(x = (orthophase - zero) / 2, y = peak + 0.01*mesor),
					label = paste0("Po = ", round(orthophase, 3))
				),

				# Bathyphase
				geom_vline(xintercept = bathyphase, color = "grey"),
				geom_segment(
					aes(x = zero, xend = bathyphase, y = trough, yend = trough),
					linetype = "dotted", size = 0.5
				),
				geom_text(
					aes(x = (bathyphase + zero) / 2, y = trough - 0.01*mesor),
					label = paste0("Pb = ", round(bathyphase, 3))
				),

				# Global Amplitude
				geom_segment(
					aes(
						x = orthophase, xend = (orthophase + bathyphase)/2,
						y = peak, yend = peak
					),
					linetype = "twodash", size = 0.5
				),
				geom_segment(
					aes(
						x = bathyphase, xend = (orthophase + bathyphase)/2,
						y = trough, yend = trough
					),
					linetype = "twodash", size = 0.5
				),
				geom_segment(
					aes(
						x = (orthophase + bathyphase)/2, xend = (orthophase + bathyphase)/2,
						y = mesor, yend = trough
					),
					linetype = "twodash", size = 0.5,
					arrow = arrow(type = "closed", length = unit(0.02, "npc"))
				),
				geom_segment(
					aes(
						x = (orthophase + bathyphase)/2, xend = (orthophase + bathyphase)/2,
						y = mesor, yend = peak
					),
					linetype = "twodash", size = 0.5,
					arrow = arrow(type = "closed", length = unit(0.02, "npc"))
				),
				geom_text(
					aes(x = (bathyphase + orthophase)/2 + max(aug$t)/4, y = 1.01*mesor),
					label = paste0("Ag = ", round((peak - trough) / 2), 3)
				)
			)

			# Return plot
			if(labels) {g <- g + gl} else {g <- g}
		}
	)

	# Dress up output
	gg <-
		g +
		theme_minimal() +
		theme(
			legend.position = "none",
			panel.grid.major = element_blank(),
			panel.grid.minor = element_blank()
		)

	# Return
	return(gg)

}

#' @noRd
#' @family cosinor
ggmulticosinor <- function(object, labels, ...) {

	# Number of cosinor objects
	n <- length(object)
	p <- length(object$tau)

	# Null variables
  mesor <- NULL
  for (i in 1:p) {
    assign(paste0("x", i), NULL)
    assign(paste0("z", i), NULL)
    assign(paste("amp", i), NULL)
    assign(paste("phi", i), NULL)
    assign(paste("beta", i), NULL)
    assign(paste("gamma", i), NULL)
  }
  .fitted <- .resid <- orthophase <- bathyphase <- trough <- ampGlobal <- .id <- peak <- NULL

	# Can be character vector or NULL
	objNames <- if(is.null(names(object))) {
		c(1:n)
	} else {
		names(object)
	}

	# Augmented data of all types
	aug <- list()
	features <- list()
	mesor <- list()
	for(i in 1:n) {
		aug[[i]] <-
			augment(object[[i]])[c("y", "t", ".fitted", ".resid")]
		features[[i]] <- cosinor_features(object[[i]])
		mesor[[i]] <- object[[i]]$coefficients[1]
	}

	# Merge together
	aug <- dplyr::bind_rows(aug, .id = ".id")
	features <-
		dplyr::bind_rows(features, .id = ".id") |>
		tibble::add_column(mesor = unlist(mesor))

	# Plot
	g <- ggplot() +
		stat_smooth(
			aes(x = t, y = .fitted, colour = .id), data = aug,
			method = "gam", size = 1.2
		) +
		geom_hline(
			aes(yintercept = mesor, colour = .id), data = features, alpha = 0.5
		) +
		geom_vline(
			aes(xintercept = orthophase, colour = .id), data = features, alpha = 0.5
		) +
		geom_vline(
			aes(xintercept = bathyphase, colour = .id), data = features, alpha = 0.5
		) +
		scale_color_viridis_d(
			option = "plasma", end = 0.8,
			name = "Objects",
			labels = objNames
		) +
		theme_minimal()

	# Labels
	gl <- list(

		# Mesor
		ggrepel::geom_label_repel(
			aes(
				x = min(aug$t), y = mesor, colour = .id,
				label = paste0("MESOR = ", round(mesor, 2))
			),
			data = features,
			show.legend = FALSE,
			label.size = NA, label.r = 0.25, label.padding = 0.25,
			force = 10,
			segment.color = "transparent"
		),

		# Peaks
		geom_point(
			aes(x = orthophase, y = peak, colour = .id), data = features,
			shape = 19, alpha = 0
		),
		ggrepel::geom_label_repel(
			aes(
				x = orthophase, y = peak, colour = .id,
				label = paste0("Peak = ", round(peak, 2))
			),
			data = features,
			show.legend = FALSE,
			label.size = NA, label.r = 0.25, label.padding = 0.25,
			force = 50,
			segment.color = "transparent"
		),

		# Troughs
		geom_point(
			aes(x = bathyphase, y = trough, colour = .id), data = features,
			shape = 19, alpha = 0
		),
		ggrepel::geom_label_repel(
			aes(
				x = bathyphase, y = trough, colour = .id,
				label = paste0("Trough = ", round(trough, 2))
			),
			data = features,
			show.legend = FALSE,
			label.size = NA, label.r = 0.25, label.padding = 0.25,
			force = 50,
			segment.color = "transparent"
		),

		# Amplitude
		ggrepel::geom_label_repel(
			aes(
				x = (orthophase + bathyphase) / 2, y = mesor + ampGlobal/2, colour = .id,
				label = paste0("Amplitude = ", round(ampGlobal, 2))
			),
			data = features,
			show.legend = FALSE,
			label.size = NA, label.r = 0.25, label.padding = 0.25,
			force = 1,
			segment.color = "transparent"
		)
	)

	if (labels) return(g + gl) else return(g)

}