#' @export
print.phylopath_summary <- function(x, ...) {
print.data.frame(x, digits = 3)
return(invisible(x))
}
#' @export
plot.phylopath_summary <- function(x, cut_off = 2, ...) {
x$model <- factor(x$model, rev(x$model))
ggplot2::ggplot(x, ggplot2::aes_(~w, ~model, fill = ~delta_CICc < cut_off, label = ~round(p, 3))) +
ggplot2::geom_col(col = 'black', alpha = 0.6) +
ggplot2::geom_text(hjust = "inward") +
ggplot2::scale_fill_manual(
values = c('TRUE' = 'firebrick', 'FALSE' = 'black'),
labels = c('TRUE' = paste('within', cut_off, 'CICc')),
breaks = c('TRUE')
) +
ggplot2::scale_x_continuous(position = 'top') +
ggplot2::coord_cartesian(expand = FALSE) +
ggplot2::guides(fill = ggplot2::guide_legend(title = NULL)) +
ggplot2::labs(y = "model weight", caption = "bar labels are p-values, signficance indicates rejection") +
ggplot2::theme(legend.position = 'bottom')
}
#' @export
print.phylopath <- function(x, ...) {
num_vars <- names(x$data)[sapply(x$data, is.numeric)]
bin_vars <- setdiff(names(x$data), num_vars)
cat('A phylogenetic path analysis, on the variables:\n')
cat('\tContinuous:\t', num_vars, '\n')
cat('\tBinary:\t\t', bin_vars, '\n')
cat('\n')
cat(' Evaluated for these models:', names(x$model_set), '\n')
cat('\n')
cat(' Containing', sum(purrr::map_dbl(x$d_sep, nrow)), 'phylogenetic regressions, of which',
length(unique(unlist(purrr::map(x$d_sep, 'd_sep')))), 'unique')
cat('\n')
}
#' Print out warnings from a phylopath analysis.
#'
#' Use this function after running `phylo_path()` to conveniently print any generated warnings
#' to the screen. You can either provide no arguments, which will only work if you run it directly
#' after the analysis, or you have to provide the `phylopath` object manually.
#'
#' @param phylopath A `phylopath` object of which the warnings should be printed.
#'
#' @export
show_warnings <- function(phylopath = NULL) {
if (is.null(phylopath)) phylopath <- .Last.value
if(!inherits(phylopath, 'phylopath'))
stop('This function expects a phylopath object.', call. = FALSE)
sink <- lapply(phylopath$warnings, warning, .call = FALSE)
return(invisible(NULL))
}
#' Plot a directed acyclic graph.
#'
#' @param x A `DAG`` object, usually created with the [define_model_set()] or [DAG()] function.
#' @param algorithm A layout algorithm from [igraph], see
#' [ggraph::create_layout()]. By default,
#' uses the Sugiyama layout algorithm, which is designed to minimize edge crossing in DAGs.
#' @param ... Not used.
#' @inheritParams plot_model_set
#'
#' @export
#'
#' @examples
#' d <- DAG(a ~ b + c + d)
#' plot(d)
#'
#' # Plot with manually defined positions:
#' ml <- data.frame(
#' name = c('a', 'b', 'c', 'd'),
#' x = c(1, 1, 2, 2),
#' y = c(1, 2, 1, 2)
#' )
#' plot(d, manual_layout = ml)
#'
plot.DAG <- function(x, labels = NULL, algorithm = 'sugiyama', manual_layout = NULL, text_size = 6,
box_x = 12, box_y = 8, edge_width = 1.5, curvature = 0.02, rotation = 0,
flip_x = FALSE, flip_y = FALSE,
arrow = grid::arrow(type = 'closed', 18, grid::unit(15, 'points')), ...) {
if (sum(x) == 0) {
stop('This DAG has no paths to plot.')
}
g <- igraph::graph_from_adjacency_matrix(x, weighted = TRUE)
l <- ggraph::create_layout(g, 'igraph', algorithm = algorithm)
if (!is.null(manual_layout)) {
l$x <- manual_layout$x[match(l$name, manual_layout$name)]
l$y <- manual_layout$y[match(l$name, manual_layout$name)]
}
l <- adjust_layout(l, rotation, flip_x, flip_y)
l <- combine_with_labels(l, labels)
ggraph::ggraph(l) +
ggraph::geom_edge_arc(
strength = curvature, arrow = arrow, edge_width = edge_width,
end_cap = ggraph::rectangle(box_x, box_y, 'mm'),
start_cap = ggraph::rectangle(box_x, box_y, 'mm')
) +
ggraph::geom_node_text(ggplot2::aes_(label = ~name), size = text_size) +
ggraph::theme_graph(base_family = 'sans')
}
#' Plot a directed acyclic graph with path coefficients.
#'
#' @param x An object of class `fitted_DAG`.
#' @param type How to express the weight of the path. Either `"width"`, or `"color"`.
#' @param algorithm A layout algorithm from \code{igraph}, see
#' [ggraph::create_layout()]. By default,
#' uses the Sugiyama layout algorithm, which is designed to minimize edge crossing in DAGs.
#' @param colors The end points of the continuous color scale. Keep in mind that red and green are
#' obvious colors to use, but are better to be avoided because of color blind users.
#' @param show.legend Whether a legend for the color scale should be shown.
#' @param width_const Deprecated.
#' @param ... Not used.
#' @inheritParams plot_model_set
#'
#' @export
#'
#' @examples
#' d <- DAG(LS ~ BM, NL ~ BM, DD ~ NL + LS)
#' d_fitted <- est_DAG(d, rhino, rhino_tree, 'lambda')
#' plot(d_fitted)
plot.fitted_DAG <- function(x, type = 'width', labels = NULL, algorithm = 'sugiyama',
manual_layout = NULL, text_size = 6, box_x = 12, box_y = 8,
edge_width = 1.25, curvature = 0.02, rotation = 0, flip_x = FALSE,
flip_y = FALSE,
arrow = grid::arrow(type = 'closed', 18, grid::unit(15, 'points')),
colors = c('firebrick', 'navy'), show.legend = TRUE,
width_const = NULL, ...) {
if (!is.null(width_const)) {
warning('width_const has been deprecated and is ignored.', call. = FALSE)
}
if (sum(x$coef) == 0) {
stop('This DAG has no paths to plot.')
}
type <- match.arg(type, c('width', 'color', 'colour'), FALSE)
if (type == 'colour') type <- 'color'
g <- igraph::graph_from_adjacency_matrix(x$coef, weighted = TRUE)
l <- ggraph::create_layout(g, 'igraph', algorithm = algorithm)
if (!is.null(manual_layout)) {
l$x <- manual_layout$x[match(l$name, manual_layout$name)]
l$y <- manual_layout$y[match(l$name, manual_layout$name)]
}
l <- adjust_layout(l, rotation, flip_x, flip_y)
l <- combine_with_labels(l, labels)
if (type == 'width') {
p <- ggplot2::ggplot(l) +
ggraph::geom_edge_arc(
ggplot2::aes_(width = ~abs(weight), color = ~weight < 0, label = ~round(weight, 2)),
strength = curvature, arrow = arrow, end_cap = ggraph::rectangle(box_x, box_y, 'mm'),
start_cap = ggraph::rectangle(box_x, box_y, 'mm'), show.legend = show.legend,
linejoin = c('bevel'), angle_calc = 'along', label_dodge = grid::unit(10, 'points')) +
ggraph::geom_node_text(ggplot2::aes_(label = ~name), size = text_size) +
ggraph::scale_edge_width_continuous(
limits = c(0, max(abs(igraph::E(g)$weight))), range = c(0, 2), guide = 'none'
) +
ggraph::scale_edge_color_manual(
name = NULL,
values = c('FALSE' = colors[2], 'TRUE' = colors[1]),
labels = c('positive', 'negative')
) +
ggraph::theme_graph(base_family = 'sans')
}
if (type == 'color') {
p <- ggplot2::ggplot(l) +
ggraph::geom_edge_arc(
ggplot2::aes_(colour = ~weight, label = ~round(weight, 2)),
edge_width = edge_width,
strength = curvature, arrow = arrow,
end_cap = ggraph::rectangle(box_x, box_y, 'mm'),
start_cap = ggraph::rectangle(box_x, box_y, 'mm'),
show.legend = show.legend,
linejoin = c('bevel'),
angle_calc = 'along',
label_dodge = grid::unit(10, 'points')
) +
ggraph::geom_node_text(ggplot2::aes_(label = ~name), size = text_size) +
ggraph::scale_edge_color_gradient2(
'standardized\npath coefficient',
low = colors[1], high = colors[2],
limits = c(-max(abs(igraph::E(g)$weight)),
max(abs(igraph::E(g)$weight))),
guide = ggraph::guide_edge_colorbar()
) +
ggraph::theme_graph(base_family = 'sans')
}
return(p)
}
#' Plot path coefficients and their confidence intervals or standard errors.
#'
#' @param fitted_DAG A fitted DAG, usually obtained by [best()], [average()] or [est_DAG()].
#' @param error_bar Whether to use confidence intervals (`"ci"`) or standard errors (`"se"`) as
#' error bars. Will force standard errors with a message if confidence intervals are not
#' available.
#' @param order_by By `"default"`, the paths are ordered as in the the model that is supplied.
#' Usually this is in the order that was established by `[phylo_path()]` for all combined graphs.
#' This can be change to `"causal"` to do a reordering based on the model at hand, or to
#' `"strength"` to order them by the standardized regression coefficient.
#' @param reverse_order If `TRUE`, the paths are plotted in reverse order.
#' Particularly useful in combination with [ggplot2::coord_flip()] to create
#' horizontal versions of the plot.
#' @param from Only show path coefficients from these nodes. Supply as a character vector.
#' @param to Only show path coefficients to these nodes. Supply as a character vector.
#'
#' @return A `ggplot` object.
#' @export
#'
#' @examples
#' d <- DAG(LS ~ BM, NL ~ BM, DD ~ NL + LS)
#' plot(d)
#' d_fitted <- est_DAG(d, rhino, rhino_tree, 'lambda')
#' plot(d_fitted)
#' coef_plot(d_fitted, error_bar = "se")
#' # to create a horizontal version, use this:
#' coef_plot(d_fitted, error_bar = "se", reverse_order = TRUE) + ggplot2::coord_flip()
coef_plot <- function(fitted_DAG, error_bar = 'ci', order_by = "default", from = NULL, to = NULL,
reverse_order = FALSE) {
stopifnot(inherits(fitted_DAG, 'fitted_DAG'))
error_bar <- match.arg(error_bar, c('ci', 'se'), several.ok = FALSE)
order_by <- match.arg(order_by, c('default', 'causal', 'strength'), FALSE)
if (error_bar == 'ci' & is.null(fitted_DAG$lower)) {
message(
'The fitted model does not contain confidence intervals, so showing standard errors instead. ',
'Fit the model with `boot` larger than 0 to get confidence intervals, or set `error_bar = "se"` ',
'to avoid this warning.'
)
error_bar <- 'se'
}
v <- colnames(fitted_DAG$coef)
df <- as.data.frame(fitted_DAG$coef)
df$from <- rownames(df)
df <- stats::reshape(df, varying = v, 'coef', direction = 'long')
df$to <- v[df$time]
if (error_bar == 'ci') {
df$lower <- c(fitted_DAG$lower)
df$upper <- c(fitted_DAG$upper)
} else {
df$lower <- c(fitted_DAG$coef - fitted_DAG$se)
df$upper <- c(fitted_DAG$coef + fitted_DAG$se)
}
df$path <- paste(df$from, df$to, sep = ' \U2192 ')
# Do the ordering of paths:
if (order_by == 'default') {
df <- df[order(match(df$from, v), match(df$to, v)), ]
}
if (order_by == 'causal') {
ordered_DAG <- fitted_DAG$coef > 0
ordered_DAG[, ] <- as.numeric(ordered_DAG)
order <- colnames(ggm::topSort(ordered_DAG))
df <- df[order(match(df$from, order), match(df$to, order)), ]
}
if (order_by == 'strength') {
df <- df[order(df$coef), ]
}
# Do the filtering of paths:
if (!is.null(from)) {
df <- df[df$from %in% from, ]
}
if (!is.null(to)) {
# make sure the variable is unambiguous
df <- df[df$to %in% to, ]
}
if (reverse_order) {
df$path <- factor(df$path, levels = rev(df$path))
} else {
df$path <- factor(df$path, levels = df$path)
}
df <- df[abs(df$coef) > .Machine$double.eps, ]
ggplot2::ggplot(df,
ggplot2::aes_(~path, ~coef, ymin = ~lower, ymax = ~upper)) +
ggplot2::geom_hline(yintercept = 0, size = 1, lty = 2) +
ggplot2::geom_pointrange(size = 0.75) +
ggplot2::xlab('') +
ggplot2::ylab(ifelse(error_bar == 'ci',
'standardized regression coefficient \U00B1 CI',
'standardized regression coefficient \U00B1 SE'))
}
#' Plot several causal hypothesis at once.
#'
#' @param model_set A list of `DAG` objects, usually created with [define_model_set()].
#' @param labels An optional set of labels to use for the nodes. This should be a named vector, of
#' the form `c(var1 = "label1", var2 = "label2")`.
#' If left at `NULL``, the variable names of the DAGs are used.
#' @param algorithm A layout algorithm from `igraph`, see
#' [ggraph::create_layout()]. By default, uses the Kamada-Kawai
#' layout algorithm. Another good option is `"sugiyama"`, which is
#' designed to minimize edge crossing in DAGs. However, it can often plot
#' nodes too close together.
#' @param manual_layout Alternatively, precisely define the layout yourself, by providing a
#' `data.frame` that at least has a column `name` with all variable names, and columns `x` and `y`
#' with positions to be plotted. Setting this parameter overrides `algorithm` but other changes,
#' such as `rotation` and `flip`s will still be applied.
#' @param text_size Size of the node label text.
#' @param box_x To avoid the arrows colliding with the nodes, specify the
#' rectangular dimensions of an invisible box around each node. If you have
#' long labels, you need to increase this.
#' @param box_y To avoid the arrows colliding with the nodes, specify the
#' rectangular dimensions of an invisible box around each node. If you have
#' multi-line labels, you need to increase this.
#' @param edge_width Width of the edges.
#' @param curvature Curvature of the edges. A slight curvature can look pretty.
#' @param rotation Supply the degrees you want to rotate the layout by. This is useful in order to
#' put rotate your upstream nodes towards the top if needed.
#' @param flip_x Whether to flip the node positions horizontally.
#' @param flip_y Whether to flip the node positions vertically.
#' @param nrow Number of rows to display the models on.
#' @param arrow A \code{grid::arrow} object, specifying the shape and size of the arrowheads.
#'
#' The order of facets is taken from the ordering of the list, with the facet
#' labels coming from the names of the list. If the list is unnamed, sequential
#' lettering is used.
#'
#' @return A `ggplot` object.
#' @export
#'
#' @examples
#' m <- list(one = DAG(a ~ b + c + d), two = DAG(a ~ b, b ~ c, d ~ d))
#' plot_model_set(m)
#' plot_model_set(m, algorithm = "sugiyama")
plot_model_set <- function(model_set, labels = NULL, algorithm = 'kk', manual_layout = NULL,
text_size = 5, box_x = 12, box_y = 10, edge_width = 1, curvature = 0.05,
rotation = 0, flip_x = FALSE, flip_y = FALSE, nrow = NULL,
arrow = grid::arrow(type = 'closed', 15, grid::unit(10, 'points'))) {
# Input checks
if (!is.list(model_set) | !all(purrr::map_lgl(model_set, ~inherits(., 'DAG')))) {
stop('model_set should be a list of DAG objects.')
}
if (is.null(names(model_set))) {
names(model_set) <- LETTERS[seq_along(model_set)]
}
var_names <- lapply(model_set, colnames)
if (length(model_set) > 1 &
(stats::var(lengths(model_set)) != 0 |
any(lengths(sapply(var_names[-1], setdiff, var_names[[1]])) != 0))) {
stop('All causal models need to include the same variables. Combined, your
models include the following variables:\n',
paste(sort(unique(unlist(var_names))), collapse = '\n'),
call. = FALSE)
}
# Build single complete graph
result <- igraph::make_empty_graph() + igraph::vertices(row.names(model_set[[1]]))
for (i in seq_along(model_set)) {
m <- model_set[[i]]
ind <- which(m == 1)
from <- ind %% nrow(model_set[[i]])
to <- (ind - from) / nrow(model_set[[i]]) + 1
result <- igraph::add_edges(result, c(rbind(rownames(m)[from], colnames(m)[to])),
attr = list(model = names(model_set)[[i]]))
}
igraph::edge.attributes(result)$model <- factor(igraph::E(result)$model,
names(model_set), names(model_set))
l <- ggraph::create_layout(result, 'igraph', algorithm = algorithm)
if (!is.null(manual_layout)) {
l$x <- manual_layout$x[match(l$name, manual_layout$name)]
l$y <- manual_layout$y[match(l$name, manual_layout$name)]
}
l <- adjust_layout(l, rotation, flip_x, flip_y)
l <- combine_with_labels(l, labels)
# Build plot.
ggraph::ggraph(l) +
ggraph::geom_edge_arc(strength = curvature, arrow = arrow, edge_width = edge_width,
end_cap = ggraph::rectangle(box_x, box_y, 'mm'),
start_cap = ggraph::rectangle(box_x, box_y, 'mm')) +
ggraph::geom_node_text(ggplot2::aes_(label = ~name), size = text_size) +
ggraph::facet_edges(~model, nrow = nrow) +
ggplot2::scale_x_continuous(expand = c(0.2, 0)) +
ggplot2::scale_y_continuous(expand = c(0.2, 0)) +
ggraph::theme_graph(foreground = 'grey80', base_family = 'sans')
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.