## ----eval=FALSE---------------------------------------------------------------
# library(ape)
# my_tree <- read.tree('my_tree.tre') # For Newick format trees
# my_tree <- read.nexus('my_tree.nex') # For NEXUS format trees

## ----eval=FALSE---------------------------------------------------------------
# rownames(my_data) <- my_data$species_name

## ----eval=FALSE---------------------------------------------------------------
# my_tree$tip.label # Check the tip labels of your tree
# rownames(my_data) <- gsub(' ', '_', my_data$species_name_with_spaces)

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(dev = "png", fig.height = 5, fig.width = 5, dpi = 300, out.width = "450px")

## -----------------------------------------------------------------------------
library(phylopath)

models <- define_model_set(
  one   = c(RS ~ DD),
  two   = c(DD ~ NL, RS ~ LS + DD),
  three = c(RS ~ NL),
  four  = c(RS ~ BM + NL),
  five  = c(RS ~ BM + NL + DD),
  six   = c(NL ~ RS, RS ~ BM),
  seven = c(NL ~ RS, RS ~ LS + BM),
  eight = c(NL ~ RS),
  nine  = c(NL ~ RS, RS ~ LS),
  .common = c(LS ~ BM, NL ~ BM, DD ~ NL)
)

## -----------------------------------------------------------------------------
models$one

## ----fig.height = 5, fig.width = 5, dpi = 300, fig.alt = "A causal diagram of model one, with arrows from body mass to litter size and nose length, from nose length to dry days, and from dry days to range size."----
plot(models$one)

## ----fig.height=8, fig.width=8, out.width = "600px", fig.alt = "A grid of nine causal diagrams, one panel per candidate model, sharing the same five variables but differing in which arrows are present."----
plot_model_set(models)

## -----------------------------------------------------------------------------
result <- phylo_path(models, data = rhino, tree = rhino_tree, model = 'lambda')

## -----------------------------------------------------------------------------
result

## -----------------------------------------------------------------------------
(s <- summary(result))

## ----fig.alt = "A dot and line plot of CICc against model, with the models ordered from best to worst supported and a dashed line marking the cut off two CICc units above the best model."----
plot(s)

## -----------------------------------------------------------------------------
(best_model <- best(result))

## ----warning = FALSE, fig.width = 6, fig.alt = "The best supported causal model, with each arrow labelled by its standardized path coefficient and drawn with a width proportional to the strength of the effect. All arrows are green, since all coefficients are positive."----
plot(best_model)

## ----fig.width = 7, fig.alt = "The conditionally averaged causal model, labelled with averaged path coefficients. It contains arrows in both directions between nose length and range size, so it is cyclical rather than a DAG."----
average_model <- average(result)
plot(average_model, algorithm = 'mds', curvature = 0.1) # increase the curvature to avoid overlapping edges

## ----fig.width = 7, fig.alt = "The fully averaged causal model. It has the same arrows as the conditionally averaged model, but the coefficients for paths present in only some models are shrunk towards zero."----
average_model_full <- average(result, avg_method = "full")
plot(average_model_full, algorithm = 'mds', curvature = 0.1)

## ----fig.alt = "A point and error bar plot of the averaged path coefficients, one path per position along the horizontal axis, with a dashed horizontal line at zero."----
coef_plot(average_model)

## ----fig.height=3.5, fig.alt = "The same coefficient plot for the fully averaged model, drawn horizontally in black and white. Several intervals now overlap zero, reflecting the shrinkage of weakly supported paths."----
coef_plot(average_model_full, reverse_order = TRUE) +
  ggplot2::coord_flip() +
  ggplot2::theme_bw()

## -----------------------------------------------------------------------------
result$d_sep$one