phylopath: Easy phylogenetic path analysis in R

Dataset

I will illustrate the use of the package by recreating a small part of the analysis by Gonzalez-Voyer et al. (2016). This study focused on the possible influence of brain size on the vulnerability to extinction in 474 mammalian species. Note that the goal of the analysis presented here is merely instructional; the original paper present a much more thorough analysis and should be used for biological inference.

The data used in the study is included in the package as red_list and red_list_tree. The data includes seven variables, listed in Table 1. Note that the species names are set as rownames and that these names match the tip labels of red_list_tree. This is how the package matches the observations to the phylogeny.

Defining the causal model set

I start out by defining various relationships common to all causal models. I assume that brain size is caused by body size (a result of allometry), gestation length is a causal parent of both litter size and weening age and that body size is a causal parent of population density, since these are all well-established relationships in the literature. I want to control for allometric effects of body size, and therefore include a direct effect of body size on status and an indirect effect through litter size. Additionally I also assume that the population density and life history variables all affect the vulnerability to extinction (which I will refer to as status), to limit the number of models that needs testing.

Since I am interested in testing for direct and indirect effects of brain size, I will vary those effects. Following the original authors, when considering indirect effects, brain size is a causal parent of litter size, gestation period and weaning age. When looking at direct effects, brain size is directly causally linked to status. This then leaves me with four causal hypotheses: a null model where brain size is irrelevant, a model with a direct effect, a model with indirect effects and a model with both.

I define these models using the define_model_set() function. I supply a list of formulas for each model, using c(). Formulas should be of the form child ∼ parent, or you can read the ∼ as “caused by,” and describe each path in your model. Multiple children of a single parent can be combined into a single formula: child ∼ parent1 + parent2. The paths that are shared between all models, can be included using the .common parameter. So I define our four models as follows:

library(phylopath)
m <- define_model_set(
 null = c(),
 direct = c(Status∼Br),
 indirect = c(L∼Br, G∼Br, W∼Br),
 both = c(Status∼Br, L∼Br, G∼Br, W∼Br),
 .common = c(Br∼B, P∼B, L∼B+G, W∼G, Status∼P+L+G+W+B)
)

It is easy to forget a path, or to make a typo. It is therefore good to make a quick plot to check. You can either plot a single model with, e.g., plot(m$direct), or plot all of them at once (Fig. 1A):

plot_model_set(m)

The nodes are laid out algorithmically. I mimic the lay-out used in the paper by manually defining the coordinates in a data.frame (Fig. 1B), which in this case looks much better:

positions <- data.frame(
 name = c('B', 'Br', 'P', 'L', 'G', 'W', 'Status'),
 x = c(2:3, c(1, 1.75, 3.25, 4), 2.5),
 y = c(3, 3, 2, 2, 2, 2, 1)
)
plot_model_set(m, manual_layout = positions, edge_width = 0.5)

Defining your model set is perhaps the most crucial part of PPA. Since the method is confirmative and not explorative, you want to strike a good balance between complexity and interpretability.

Evaluation of the hypotheses

p <- phylo_path(m, red_list, red_list_tree)

Printing the result gives us some basic information:

p
## A phylogenetic path analysis, on the variables:
## Continuous: G W B L P Status Br
## Binary:
##
## Evaluated for these models: null direct indirect both
##
## Containing 36 phylogenetic regressions, of which 18 unique

More importantly, asking for its summary and plotting it (Fig. 2) gives me the actual result of our comparison:

s <- summary(p)
s
##   model k q    C   p  CICc delta_CICc   l   w
## 1 indirect 8 20 19.205 0.258 61.059   0.000 1.000 0.696
## 2   both 7 21 18.671 0.178 62.715   1.656 0.437 0.304
## 3   null 11 17 247.625 0.000 282.967  221.908 0.000 0.000
## 4  direct 10 18 247.090 0.000 284.594  223.535 0.000 0.000
plot(s)

The summary reports the results table as used by Gonzalez-Voyer & von Hardenberg (2014). Specifically, it reports the model name, the number of independence claims made by the model (k), the number of parameters (q), the C statistic and the accompanying p-value. A significant p-value would indicate that the available evidence rejects the model. It also reports model selection information: the C-statistic information criterion corrected for small sample sizes (CICc), the difference in CICc with the top model (delta_CICc) and finally the associated relative likelihoods (l) and CICc weights (w).

In this example, there is strong support for the indirect pathway. The addition of the direct path in the both model did lead to a small improvement (the C-statistic is lower) but not enough to put it ahead of the indirect model.

Choosing a final model

So what is the best causal model? Firstly, the null and direct models are not supported since they have significant p-values and should therefore be discarded. The indirect pathway is certainly important, but what about the direct pathway? There are several philosophies of dealing with this issue. In this particular case the two top-ranked causal models are directly nested, they share all the same paths except for one. One can think of this like nested regression models. Typically, the extra path should lower the CICc by at least some margin, often two. In this case it does not and I elect to choose the top ranked model (see Arnold, 2010 for a discussion on AIC and uninformative parameters).

After I have found my final model, I can estimate the relative importance of each of the paths. To estimate the paths in the highest ranked model, use the best function:

b <- best(p)
plot(b, manual_layout = positions)

This will return both the standardized regression coefficients, as well as their standard errors. The resulting plot is shown in Fig. 3. In order to get confidence intervals as well, you need to take bootstrap replicates using the boot argument: e.g. b_ci <- best(p, boot = 500), which uses the bootstrap methods of the phylolm package (see “Implementation Notes”). This is disabled by default because it is slow. Using plot will give a visualization of the causal model. You can fit any arbitrary causal model that you evaluated with choice, so in this case choice(p, "both") would give the second ranked model.

A second way to look at a fitted model is to more directly look at the standardized coefficients and errors of the paths using coef_plot. This can be used it to quickly compare the importance of the different variables that affect Status. Although I have modeled five effects on status, they are not necessarily all important and certainly litter size and body size have small effects (Fig. 4A).

coef_plot(b, error_bar = "se", order_by = "strength", to = "Status") + ggplot2::coord_flip()

Model averaging

In many cases it may not be obvious or correct to choose one model. While in this case the two top competing models were nested, they do not have to be. In cases like these, it may be useful to perform model averaging instead, as discussed and used in the original paper (von Hardenberg & Gonzalez-Voyer, 2013). phylopath makes model averaging easy, and you can quickly average over a selection of the top models, or all models considered. One should take care to not include models with significant C-statistics in the averaging, as these models are not supported. Models are weighted by their likelihood, and these weights can be found in the original summary table in the w column. One needs to choose how to deal with paths that do not occur in all models. One can average path coefficients only between those models that include that path. This is often called conditional averaging and was used by von Hardenberg & Gonzalez-Voyer (2013) and is the default behavior in phylopath. Alternatively, one can consider missing paths to have a coefficient of zero and average over all models, which is often called full averaging. The latter results in shrinkage, where the path coefficients that do not occur in all models will shrink toward zero.

In this case, I could choose to average the two competing models. I use full averaging, as I would like uncertain paths to experience shrinkage, and re-evaluate the strength of the coefficients toward Status (Fig. 4B):

avg <- average(p, avg_method = "full")
coef_plot(avg, error_bar = "se", order_by = "strength", to = "Status") + ggplot2::coord_flip()

The average function selects the competing models, estimates the standardized path coefficients and then averages them. Note that only the two top models have been averaged, since by default the cut_off is set to two CICc. You can average over all models in the set by using cut_off = Inf (but should only do so when all C-statistics are non-significant, see above).

Analysis conclusion

A clear rejection of the null model indicates that brain size is related to the vulnerability to extinction of mammals, where large-brained animals high a higher vulnerability. This effect is mediated through life history, where the weaning and gestation periods are more important than litter size. There is no strong evidence in support of a direct effect of brain size on vulnerability to extinction that is independent of life history. The original analysis came to the same conclusion.