---
title: "Comparing causal models of binary traits using phylopath"
author: "Wouter van der Bijl"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Comparing causal models of binary traits using phylopath}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Introduction

This vignette gives a short example of how PPA can be applied to binary data sets using `phylopath`. A longer example with more explanation of the code can be found in the other vignette, "intro to phylopath".

### On phylogenetic logistic regression

There has been some discussion concerning how to best perform logistic regression with phylogenetic correction. This package uses `phylolm::phyloglm`, written by Lam Si Tung Ho, Robert Lachlan, Rachel Feldman and Cécile Ané. `phylopath`'s accuracy is directly dependent on the accuracy of that function, and if you don't trust `phyloglm` you should not trust binary models used in `phylo_path`.

`phylolm::phyloglm` performs checks for model convergence. In practice, these often fail. Especially in `phylopath`, where we often fit many models, it is likely that at least one model fails these checks and generates a warning. You can see the warnings using `show_warnings()`. You can try changing to the second method (`method = 'logistic_IG10'`), or passing parameters like `btol` and `log.alpha.bound` (see `?phylolm::phyloglm`). These parameters are then applied to all models. Another option is to use the `phylosem` package, which supports both Binomial and Poisson errors with a different implementation.

*The example below currently generates such warnings, and is only presented as an example.*

If you have useful opinions or information on these points, feel free to contact me.

### On the interpretation of coefficients

Coefficients from binary models are not the same quantity as coefficients from continuous ones, and it is worth being clear about this before looking at any output.

`phylopath` standardizes continuous variables before fitting, that is, it centers them and divides them by their standard deviation. So a path between two continuous variables is expressed in standard deviations: it is the familiar standardized regression coefficient. Binary variables are left as they are, so a path leading to a binary variable is a log odds ratio instead. In the example below every variable is binary, which means every coefficient is a log odds ratio between the two levels of a trait. Those are all the same kind of quantity, so they can be compared with each other directly, as we do below.

Where this becomes genuinely confusing is when a single causal model contains both binary *and* continuous traits. Such a model has coefficients in two different units, which is bad enough. On top of that, a binary predictor contributes a contrast between its two levels rather than a change per standard deviation, and for a trait where a proportion `p` of species has one of the levels, that contrast amounts to `1 / sqrt(p * (1 - p))` standard deviations, which is always more than two. Coefficients of this kind of mixed-type DAG therefore cannot be meaningfully compared with one another at all, and each has to be read on its own terms. `phylopath` will try to warn you when you fit such a model, and notes it on any output that shows its coefficients.

Importantly, this concerns the coefficients only. The comparison of the causal models themselves is unaffected, because the d-separation tests are performed on the unstandardized data: `CICc`, the model weights and the ranking of your models do not depend on the scale of the coefficients. See `?est_DAG` for more detail.

# Example analysis

### Data and hypotheses

This recreates the analysis from the following paper:

> Dey CJ, O'Connor CM, Wilkinson H, Shultz S, Balshine S & Fitzpatrick JL. 2017. Direct benefits and evolutionary transitions to complex societies. Nature Ecology & Evolution. 0137.

This is, to my knowledge, the first study to employ PPA on binary traits.

The study investigates the evolution of cooperative breeding in cichlids. In short (my summary), there has been intense debate about what factors drive species towards evolving systems of cooperative breeding. Many have argued (and provided evidence in birds and mammals) that cooperative breeding chiefly evolves from monogamous mating systems because helpers can gain indirect fitness benefits through kin selection. However, a non-exclusive alternative hypothesis is that direct benefits due to ecological factors may be important and provide direct benefits. Therefore, both hypotheses should be considered at the same time.

The data is included in this paper as `cichlids` and `cichlids_tree`

It contains five variables:

-   C: Cooperative breeding, cooperative or non-cooperative.
-   M: Mating system, monogamous or non-monogamous.
-   P: Parental care, maternal or biparental.
-   G: Social grouping, grouping or non-grouping.
-   D: Diet, piscivore or non-piscivore.

Under the indirect fitness hypothesis, monogamy is expected to be a major driver of cooperative breeding, while group living, biparental care and diet type may be important contributors towards a direct benefits scenario.

### Defining the causal models

Following the paper in question, we define 12 putative causal models.

```{r 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')
```

### Comparing the models

Now that we have our models, data and a tree, we can compare the models using `phylo_path`

```{r fit models}
(cichlids_results <- phylo_path(models, cichlids, cichlids_tree))
```

Note that three messages are printed. This is because there are missing values in our data set that are first being removed. Also, since the tree includes species for which data is missing, the tree had to be pruned. This is done automatically with a message to the user. You should check whether the amount of data removed is correct.

The third messages notes that warnings have been generated, see the **Important notes** at the start of this vignette. I will continue with these results here, but *models with warnings should not be trusted*.

`phylo_path` notes that indeed all variables are binary.

```{r 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)
```

We see that model F is the best supported model. This model notably does not include a link between monogamy and cooperative breeding, giving support to the direct benefits hypothesis.

### Evaluating variable importance

Now that we have selected F as our best model, we still have three factors that affect cooperative breeding: diet, social grouping and parental care. Which one is more important? For this we can fit the model and look at magnitude of the coefficients. In this case, since we want to use the best model we use the function `best()`. One can use `choice()` to choose any arbitrary model, or `average()` to average over several models.

```{r}
best_cichlids <- best(cichlids_results)
```

To see the individual coefficients and their standard errors, simply print `best_cichlids`:

```{r}
best_cichlids
```

Or plot those:

```{r, 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()
```

But we can also plot the final model:

```{r, 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)
```
