## ----define_models, fig.align='center', fig.width=10, fig.height=8, out.height="600px", fig.dpi = 600, fig.alt = "A grid of twelve causal diagrams, one panel per candidate model, each relating cooperative breeding, mating system, parental care, social grouping and diet with directed arrows."----
library(phylopath)

models <- define_model_set(
  A = c(C~M+D),
  B = c(C~D),
  C = c(C~D, P~M),
  D = c(C~D, M~P, G~P),
  E = c(C~D, P~M, G~P),
  F = c(C~D, P~M+G),
  G = c(C~D, M~P, P~G),
  H = c(C~D, M~P),
  I = c(C~D, M~M, G~P),
  J = c(M~P, G~D),
  K = c(P~M, G~D),
  L = c(C~M+D, P~M+G),
  .common = c(C~P+G)
)

plot_model_set(models, algorithm = 'kk')

## ----fit models---------------------------------------------------------------
(cichlids_results <- phylo_path(models, cichlids, cichlids_tree))

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

## -----------------------------------------------------------------------------
best_cichlids <- best(cichlids_results)

## -----------------------------------------------------------------------------
best_cichlids

## ----fig.alt = "A horizontal point and error bar plot of the path coefficients of the best supported model, with standard errors and a dashed line at zero."----
coef_plot(best_cichlids, error_bar = "se", reverse_order = TRUE) + ggplot2::coord_flip()

## ----fig.align='center', fig.width=8, fig.height=4, out.width="600px", fig.dpi = 300, fig.alt = "The fitted best supported model, with arrows labelled by their coefficients and drawn with a width proportional to the strength of the effect."----
plot(best_cichlids)