A Primer in MV-PMM

Introduction

This tutorial is associated with the article “A Primer in Multivariate Phylogenetic Mixed Modelling”

EcoEvoRxiv preprint DOI XXXXXXX

by

Ben Halliwell, Luke Yates & Barbara Holland

In this tutorial we will cover how to implement phylogenetic mixed models (PMM) in two popular R packages, MCMCglmm and brms. PMM is motivated by the realisation that, because closely related species tend to resemble each other phenotypically (a phenomenon known as phylogenetic signal), we cannot assume that species represent independent data points with respect to evolutionary hypotheses. Instead, we must test for non-independence and weight our analyses accordingly. A central assumption of PMM therefore, is that we can combine a molecular phylogeny with a model of evolution (nominally, Brownian Motion (BM)) to generate a covariance matrix that defines the expected phenotypic similarity among taxa. The strength of this signal is optimised from the data during model fit, scaling the magnitude of effects to produce a phylogenetically informed regression. For univariate (UV) cases, this approach only estimates phylogenetic signal in the response variable, while any predictor variables in the model are assumed not to display phylogenetic signal. This assumption is commonly violated in analyses of real world data sets, where predictor variables often represent species traits that are themselves subject to phylogenetic effects. This problem is effectively addressed by moving offending predictor variables to the LHS of the model equation in multivariate (MV) implementations. This way, phylogenetic signal is simultaneously evaluated for all relevant traits, allowing (co)variances between response variables to be partitioned across multiple levels in the model hierarchy (e.g. phylogenetic and residual). This partitioning is important because phenotypic covariances (i.e., covariance on the level of traits measured), Cov(y1,y2), represent the sum of covariances on hierarchically lower levels, which need not align in magnitude or even direction. Thus, for comparative analyses of species traits with sufficient replication, MV-PMM should be preferred, as it offers more informative outputs with respect to evolutionary hypotheses.

While frequentist implementations of the PMM are available (e.g., pglmm() in the “phyr” package), we advocate for a Bayesian approach due to 1) greater flexibility in fitting non-Gaussian response traits 2) superiority in estimating variance components associated with hierarchical random effects; 3) the ability to sample conditional multivariate distributions, avoiding intractable integration problems faced by frequentist techniques.

Throughout this tutorial, we purposefully explore a low dimensional example (tree with only 5 taxa) in order to show full workings of the matrix manipulations that underlie the models specified. This is for the sake of exposition only; analyses of the form presented will typically require much higher phylogenetic replication to produce reliable estimates of, for example, phylogenetic correlations (Housworth et al. 2004). Thus, while the mathematical workings presented follow a 5 taxa example, model outputs are based on analyses from simulations including 250 taxa.

Model 1 - Univariate Gaussian

Model Explanation

For the case of a single Gaussian trait, our response variable, \(y\), is modeled directly as a linear combination of the fixed effects \(\mu\), a phylogenetic random effect, \(b\), and residual variance, \(e\).

\[\begin{eqnarray} \mathbf{y} &=& \boldsymbol{\mu} + \mathbf{b} + \mathbf{e}\\[2mm] \end{eqnarray}\]

where \(\mu\) represents,

\[\begin{eqnarray} \boldsymbol{\mu} = \beta_0 + \beta_1\boldsymbol{X_1} + \ldots + \beta_{k}\boldsymbol{X_{k}} \end{eqnarray}\]

In order to capture the expected influence of co-ancestry on \(y\), (similar values among closely related taxa), the distribution of our random effect, \(b\), is multivariate normal, with mean 0 and (co)variance given by \(\sigma_{b}A\),

\[\begin{eqnarray} b &\sim& \mbox{MVN}(0, \sigma^2_{b}A)\\ \end{eqnarray}\]

Here, \(A\) is a fixed \(n\times n\) phylogenetic VCV matrix, where \(n\) is the number of taxa in the tree, and \(\sigma_{b}\), the phylogenetic variance, is a scalar to be estimated by the model. This model therefore takes as input a phylogeny (topology and branch lengths), which is assumed to be known without error. The \(A\) matrix is derived by taking the inverse of the distance matrix of the phylogeny. This gives the expected (co)variance among taxa in \(y\) assuming BM evolution.

\[\begin{eqnarray} A &=& \begin{pmatrix}A_{11} & A_{12} & \ldots & A_{1n} \\ A_{21} & A_{22} & \ldots & A_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ A_{n1} & A_{n2} & \ldots & A_{nn} \\ \end{pmatrix} \end{eqnarray}\]

Our observation-level residual error, \(\mathbf{e}\), on the other hand, is assumed to be independent of species identity. This assumption is encoded by the identity matrix \(I\) (an \(n\times n\) matrix with 1 in all diagonal (variance) elements and 0 in all off-diagonal (covariance) elements) in the distributional statement:

\[\begin{eqnarray} \mathbf{e} &\sim& \mbox{MVN}(0, \sigma^2_{e}I)\\ \end{eqnarray}\]

Example

To cement these concepts, consider the following 5-taxon phylogeny with genetic distance, or dissimilarity, shown on the x axis:

t.toy2 <- read.tree(text='((1:1, 2:1):3, (3:3, (4:2, 5:2):1):1);')
plot(t.toy2, label.offset=0.1, edge.width=1.5)
axis(1, at = seq(0,4,by=1), labels = rev(seq(0,4,by=1)), line = 0.5)
title(xlab="distance")

For this phylogeny, the inverse distance matrix (elements expressing the sum of shared branch lengths between two taxa, \(i\) and \(j\)), \(A\), is:

\[\begin{eqnarray} A &=& \begin{pmatrix}A_{11} & A_{12} & A_{13} & A_{14} & A_{15} \\ A_{21} & A_{22} & A_{23} & A_{24} & A_{25} \\ A_{31} & A_{32} & A_{33} & A_{34} & A_{35} \\ A_{41} & A_{42} & A_{43} & A_{44} & A_{45} \\ A_{51} & A_{52} & A_{53} & A_{54} & A_{55} \\ \end{pmatrix} &=& \begin{pmatrix} 4 & 3 & 0 & 0 & 0 \\ 3 & 4 & 0 & 0 & 0 \\ 0 & 0 & 4 & 1 & 1 \\ 0 & 0 & 1 & 4 & 2 \\ 0 & 0 & 1 & 2 & 4 \\ \end{pmatrix} &=& \begin{pmatrix} &&&&\\&&&&\\&&A_{ij}&&\\&&&&\\&&&& \end{pmatrix} \\ \end{eqnarray}\]

In order to facilitate comparison between phylogenetic and residual contributions when partitioning the variance, \(A\) is standardised to the unit diagonal (i.e., re-scaled to a correlation matrix) prior to analyses, such that,

\[\begin{eqnarray} A_{cor} &=& \begin{pmatrix} 1 & 0.75 & 0 & 0 & 0 \\ 0.75 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0.25 & 0.25 \\ 0 & 0 & 0.25 & 1 & 0.5 \\ 0 & 0 & 0.25 & 0.5 & 1 \\ \end{pmatrix} \end{eqnarray}\]

The identity matrix \(I\) for this model is:

\[\begin{eqnarray} I &=& \begin{pmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{pmatrix} \end{eqnarray}\]

Data Simulation and Model Specification

Implementing these models in R is fairly straight forward in both MCMCglmm and brms. But first, let’s simulate some data. To simulate data under this model, we simply construct our variable, \(y\), from the linear combination of predictors specified by the model, by supplying vectors of length \(n\) for each of the quantities \(\mu\), \(b\), and \(e\). For Gaussian variables this can be done directly, while for non-Gaussian variables we must make predictions on the link scale (see section 3). In this example, we will simulate from an intercept only model (\(\mu\) = \(\beta_0\)).

## SIMULATION

# simulate tree
n = 250 # number of species
t.toy <- geiger::sim.bdtree(b=1, d=0, stop="taxa", n=n, extinct=FALSE) # N.B. may be best not to use pure birth trees (Adams and Collyer 2018)
t.toy <- multi2di(t.toy) # ensure tree is fully bifurcating

# create VCV matrix from tree, scale to correlation matrix with corr = T
A.mat <- ape::vcv.phylo(t.toy, corr = T)

# phylogenetic VCV matrix
A = A.mat

# identity matrix
I = diag(n)

# fixed effects (intercept only)
u = 0

# phylogenetic variance
sigma_phy = 0.5

# residual variance
sigma_res = 0.1

# simulate phylogenetic random effects as one draw from a MVN
b = MASS::mvrnorm(1,rep(0,n),sigma_phy*A)

# simulate residuals
e = MASS::mvrnorm(1,rep(0,n),sigma_res*I)

# construct response from the linear predictor
y = u + b + e

# generate df
animal <- t.toy$tip.label # "animal" is a reserved term in MCMCglmm for taxon ID
d.toy <- data.frame(animal,y)

## PLOTS TO EXPLORE THIS SIMULATED DATA?

With our tree and trait data simulated, we can now fit the model shown in ():

## MCMCglmm 
p.m.1 <- list(G = list(G1 = list(V = 1, nu = 0.002)), 
              R = list(V = 1, nu = 0.002))
m.1 <- MCMCglmm(y ~ 1,
                random = ~animal,
                pedigree = t.toy, 
                family = c("gaussian"), 
                data = d.toy, 
                prior = p.m.1,
                nitt = 210000, burnin = 10000, thin = 200,
                pr = TRUE, verbose = FALSE)


## BRMS
p.b.1 <- list()
b.1 <- brm(y ~ 1 + (1|gr(animal, cov = A)), 
            data = d.toy,
            prior = p.b.1
            family = gaussian(), 
            data2 = list(A = A.mat),
            cores = 4, chains = 4, iter = 6000, thin = 3)

Model Syntax (MCMCglmm vs. brms)

The code above fits the same model in MCMCglmm and brms. In MCMCglmm, animal is a reserved term used to identify individuals/taxa in a quantitative genetic/phylogenetic analysis. The random effect specification ~animal instructs MCMCglmm to fit a random effect at the individual/taxon level with covariance structure supplied via the pedigree argument. For pedigree, MCMCglmm accepts a phylo object directly and performs the necessary steps to derive the inverse distance matrix \(A\) under the hood. The remaining code specifies the sampler settings, whether the posterior distribution of random effect should be saved and whether to print sampler progress to the console.

In brms, the phylogenetic random effect can be encoded in the style of lme4 (1|gr(animal, cov = A)), where (1|x) specifies random intercepts and gr(animal, cov = A) encodes the phylogenetic group structure with covariance \(A\). It can also be encoded by a call to cov_ranef, e.g. ... + (1|animal), cov_ranef = list(animal = A), ....

It is not necessary to specify a prior to fit this model in either MCMCglmm or brms; if no prior is specified, default priors are used automatically. However, we choose to explicitly specify (current default) priors to ensure repeatability, should the default priors of either package change in the future. Also, because a notable difference between these packages is that MCMCglmm uses Gibbs sampling, where-as brms uses Hamiltonian sampling. These different computational approaches place different requirements on prior specification. In practice, this can mean a different level of knowledge and coding is required of the user to effectively tune models and avoid convergence issues. MCMCglmm uses inverse Wishart priors for (co)variance structures G and R (see MCMCglmm course notes for details), where-as as brms uses….

Diagnostics

Basic diagnostics are easily achieved for each package

# visual check for convergence in MCMC chains
plot(m.1$VCV)
b.1 %>% plot

# calculate autocorrelation in MCMC draws
autocorr(m.1$VCV)
# brms equiv?

Conveniently, brms reports \(\widehat{R}\) statistics, or potential scale reduction factors, for each parameter in the model summary to assess MCMC chain convergence. It is also very straightforward to run posterior predictive checks in brms.

# posterior predictive checks
b.1 %>% pp_check(nsamples = 100)

Inference

We can compare model outputs to confirm that we are indeed fitting the same model.

summary(m.1) # MCMCglmm
summary(b.1) # brms

Intercept estimates are directly comparable between MCMCglmm and brms. However, MCMCglmm reports variance components where-as brms reports these quantities as standard deviations. As our simulations also specify standard deviations, we will re-scale estimates from MCMCglmm by calculating posterior means from the sqrt() of MCMC samples.

## PARAMETER ESTIMATES ##

# intercept
summary(m.1)$solutions
summary(b.1)[["fixed"]]

# phylogenetic variance (sigma_phy) = sqrt(phy_var)?
sqrt(m.1$VCV[,"animal"])
summary(b.1)[["random"]]

# residual variance (sigma_res)
sqrt(summary(m.1)$Rcovariances)
summary(b.1)[["spec_pars"]]

An estimate of phylogenetic signal, the amount of variation in a trait explained by shared ancestry between species, in \(y\) can be calculated as

\[\begin{eqnarray} h^2 = \sigma^2_{b}/(\sigma^2_{b}+\sigma^2_{e}) \end{eqnarray}\]

# MCMCglmm
mean(m.1$VCV[,"animal"]/(m.1$VCV[,"animal"]+m.1$VCV[,"units"]))

# brms (clean this up)
mean(b.1$fit@sim$samples[[1]]$sd_animal__g1_Intercept/(b.1$fit@sim$samples[[1]]$sd_animal__g1_Intercept+b.1$fit@sim$samples[[1]]$sigma_g1))

We can also confirm that \(h^2\) from our PMM is estimating the same quantity as the MLE of lambda from an equivalent pgls model (see Halliwell et al. 2022 for details)

## compare to MLE of lambda from PGLS fit
comp_dat <- comparative.data(t.toy, d.toy, animal, vcv=TRUE)
mod <- pgls(y ~ 1, data = comp_dat, lambda = "ML")
summary(mod)$param["lambda"]

Model 2 - Multivariate Gaussian

Model Explanation

In a multivariate (MV) implementation, our response variables are modelled jointly, allowing both phylogenetic and residual (co)variances to be estimated. With this design, the phylogenetic random effects must also be realised jointly as a multivariate normal, with variance given by \(\Omega\), the Kronecker product of a trait-level correlation matrix \(\Sigma^{\mathrm{b}}\) and the phylogenetic VCV matrix, \(A\). With \(m\) traits (response variables) and \(n\) taxa in the phylogeny, \(\Omega\) is an \(mn\times mn\) VCV matrix containing phylogenetic (co)variances for all traits as well as between all traits. Notice that the two taxon-level variance structures we encountered in the univariate case above (\(A\) and \(I\)) remain, but we must now also consider trait-level covariance operating at both the phylogenetic and residual level. This is achieved via the Kronecker operation, e.g. \(\Omega = \Sigma^{\mathrm{b}}\otimes A\). A MV-PMM with two Gaussian response variables takes the form,

\[\begin{eqnarray} \begin{pmatrix}\mathbf{y}_1 \\ \mathbf{y}_2 \end{pmatrix} &=& \begin{pmatrix}\boldsymbol{\mu}_1 + \mathbf{b}_1 + \mathbf{e}_1 \\ \boldsymbol{\mu}_2 + \mathbf{b}_2 + \mathbf{e}_2 \end{pmatrix}\\ \end{eqnarray}\]

where fixed effects are expressed as a linear combination of predictors,

\[\begin{eqnarray} \begin{pmatrix}\boldsymbol{\mu}_1 \\ \boldsymbol{\mu}_2 \end{pmatrix} &=& \begin{pmatrix}\beta_{0,1}\mathbf{1} + \beta_{1,1}\boldsymbol{X_{1,1}} + \ldots + \beta_{k,1}\boldsymbol{X_{k,1}} \\ \beta_{0,2}\mathbf{1} + \beta_{1,2}\boldsymbol{X_{1,2}} + \ldots + \beta_{k,2}\boldsymbol{X_{k,2}} \end{pmatrix}\\ \end{eqnarray}\]

and the phylogenetic random effects and residuals are drawn from multivariate normal distributions,

\[\begin{eqnarray} \begin{pmatrix}\mathbf{b}_1 \\ \mathbf{b}_2 \end{pmatrix} &\sim& \mbox{MVN}(0, \Omega)\\ \begin{pmatrix}\mathbf{e}_1 \\ \mathbf{e}_2 \end{pmatrix} &\sim& \mbox{MVN}(0, \Sigma)\\[2mm] \Omega &=& \Sigma^{\mathrm{b}}\otimes A \\ \Sigma &=& \Sigma^{\mathrm{e}}\otimes I \end{eqnarray}\]

Here, \(\Sigma^{\mathrm{b}}\) is an \(m\times m\) phylogenetic trait-level VCV matrix and \(m\) is the number of response traits in the model,

\[\begin{eqnarray} \Sigma^{\mathrm{b}} = \begin{pmatrix}\Sigma^{\mathrm{b}}_{11} & \Sigma^{\mathrm{b}}_{12} & \ldots & \Sigma^{\mathrm{b}}_{1m} \\ \Sigma^{\mathrm{b}}_{21} & \Sigma^{\mathrm{b}}_{22} & \ldots & \Sigma^{\mathrm{b}}_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ \Sigma^{\mathrm{b}}_{m1} & \Sigma^{\mathrm{b}}_{m2} & \ldots & \Sigma^{\mathrm{b}}_{mm} \\ \end{pmatrix} &=& \begin{pmatrix}(\sigma^{\mathrm{b}}_1)^2 & \rho^{\mathrm{b}}_{12}\sigma^{\mathrm{b}}_1 \sigma^{\mathrm{b}}_2 & \ldots & \rho^{\mathrm{b}}_{1m}\sigma^{\mathrm{b}}_1 \sigma^{\mathrm{b}}_m \\ \rho^{\mathrm{b}}_{21}\sigma^{\mathrm{b}}_2 \sigma^{\mathrm{b}}_1 & (\sigma^{\mathrm{b}}_2)^2 & \ldots & \rho^{\mathrm{b}}_{2m}\sigma^{\mathrm{b}}_2 \sigma^{\mathrm{b}}_m \\ \vdots & \vdots & \ddots & \vdots \\ \rho^{\mathrm{b}}_{m1}\sigma^{\mathrm{b}}_m \sigma^{\mathrm{b}}_1 & \rho^{\mathrm{b}}_{m2}\sigma^{\mathrm{b}}_m\sigma^{\mathrm{b}}_2 & \ldots & (\sigma^{\mathrm{b}}_m)^2 \\ \end{pmatrix}\\[2mm] \end{eqnarray}\]

such that elements of \(\Sigma^{\mathrm{b}}\) have the general form,

\[\begin{eqnarray}\ \Sigma^{\mathrm{b}}_{ij} &=& \rho^{\mathrm{b}}_{ij}\sigma^{\mathrm{b}}_i \sigma^{\mathrm{b}}_j \\ \end{eqnarray}\]

where \(\rho_{ij}\) is the correlation between traits \(i\) and \(j\) and \(\sigma_i\) is the standard deviation of trait \(i\). Thus, returning to our 5-taxon example, a MV model with two Gaussian response traits (i.e., \(m = 2\)) will have VCV matrix for the phylogenetic random effects, \(\Omega\), of the form:

\[\begin{eqnarray} \Omega &=& \Sigma^{\mathrm{b}} \otimes A &=& \begin{pmatrix} \Sigma^{\mathrm{b}}_{11}A & \Sigma^{\mathrm{b}}_{12}A \\ \Sigma^{\mathrm{b}}_{21}A & \Sigma^{\mathrm{b}}_{22}A \\ \end{pmatrix} &=& \begin{pmatrix} \Sigma^{\mathrm{b}}_{11} \left(\begin{smallmatrix} 4 & 3 & 0 & 0 & 0 \\ 3 & 4 & 0 & 0 & 0 \\ 0 & 0 & 4 & 1 & 1 \\ 0 & 0 & 1 & 4 & 2 \\ 0 & 0 & 1 & 2 & 4 \\ \end{smallmatrix}\right) & \Sigma^{\mathrm{b}}_{12}\begin{pmatrix} &&\\&A_{ij}&\\&& \end{pmatrix} \\ \Sigma^{\mathrm{b}}_{21}\begin{pmatrix} &&\\&A_{ij}&\\&& \end{pmatrix} & \Sigma^{\mathrm{b}}_{22}\begin{pmatrix} &&\\&A_{ij}&\\&& \end{pmatrix} \\ \end{pmatrix} \end{eqnarray}\]

Similarly, the VCV matrix for the residual errors, \(\Sigma\), is:

\[\begin{eqnarray} \Sigma &=& \Sigma^{\mathrm{e}} \otimes I &=& \begin{pmatrix} \Sigma^{\mathrm{e}}_{11}I & \Sigma^{\mathrm{e}}_{12}I \\ \Sigma^{\mathrm{e}}_{21}I & \Sigma^{\mathrm{e}}_{22}I \\ \end{pmatrix} &=& \begin{pmatrix} \Sigma^{\mathrm{e}}_{11} \left(\begin{smallmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{smallmatrix}\right) & \Sigma^{\mathrm{e}}_{12}\begin{pmatrix} &&\\&I_{ij}&\\&& \end{pmatrix} \\ \Sigma^{\mathrm{e}}_{21}\begin{pmatrix} &&\\&I_{ij}&\\&& \end{pmatrix} & \Sigma^{\mathrm{e}}_{22}\begin{pmatrix} &&\\&I_{ij}&\\&& \end{pmatrix} \\ \end{pmatrix} \end{eqnarray}\]

N.B. This example is purely for illustrative purposes. The parameters of this model are not estimable with such low phylogenetic replication (n = 5 taxa). However, the low dimension of \(A\) and \(I\) allow us to explicitly show the working of matrix manipulations.

Data Simulation and Model Specification

With multiple response variables, we must now consider covariance between responses at both the phylogenetic and residual level. To incorporate this into our data simulations, we can specify correlation matrices which define the (co)variance structure of our random effects and residual errors.

## SIMULATE DATA

# number of traits
k = 2

# fixed effects (intercepts for each trait)
u = c(1,2)

## B is the phylogenetic trait-level VCV matrix (\Sigma^{\mathrm{b}}). B specifies the phylogenetic variance in each trait (diagonals) as well as the phylogenetic covariance between traits (off-diagonals). Each unique element in B is to be estimated by the model. To construct B, we must first define its (co)variance components.
sig.B <- c(b11 = 0.4, b22 = 0.6) # sqrt of the phylogenetic variance for each trait
b12_rho = 0.5 # phylogenetic correlation coefficient
Bcor <- matrix(c(c(1,b12_rho), # phylogenetic correlation matrix
                 c(b12_rho,1)),k,k, byrow = T) 

# Construct B as point-wise product
# N.B. Kronecker used here just for ease of matrix formatting. Do not confuse with Kronecker operation described in the text for specifying the covariance of random effects a and residuals e (below)
B <- matrix(kronecker(sig.B, sig.B),k,k)*Bcor 

## C is the residual trait-level VCV matrix (\Sigma^{\mathrm{e}}).
sig.C <- c(c11 = 0.1, c22 = 0.1) # sqrt of the residual variance for each trait
c12_rho = 0 # residual correlation coefficient
Ccor <- matrix(c(c(1,c12_rho), # residual correlation matrix
                 c(c12_rho,1)),k,k, byrow = T) 
C <- matrix(kronecker(sig.C, sig.C),k,k)*Ccor

# simulate phylogenetic random effects
a = mvrnorm(1,rep(0,n*k),kronecker(B,A)) # In the Kronecker, trait-level covariance captured by B, taxon-level covariance captured by A. 

# extract random effects for each trait from the vector, a.
a1 <- a[1:n]
a2 <- a[1:n + n]

# simulate residuals (on the link scale) and extract
e = mvrnorm(1,rep(0,n*k),kronecker(C,I))
e1 <- e[1:n]
e2 <- e[1:n + n]

# construct response traits from each linear predictor
y1 = u[1] + a1 + e1
y2 = u[2] + a2 + e2

# generate df
species <- t.toy$tip.label
d.toy <- data.frame(species,y1,y2)
d.toy$animal <- d.toy$species
d.toy$obs <- 1:nrow(d.toy)

A MV model with two Gaussian response traits, such as the example provided above, can be specified in MCMCglmm and brms as follows:

## MCMCglmm 
p.m.2 <- list(G = list(G1 = list(V = diag(2), nu = 1.002)), 
              R = list(V = diag(2), nu = 1.002))
m.2 <- MCMCglmm(cbind(y1, y2) ~ trait-1,
                random = ~us(trait):animal,
                rcov = ~us(trait):units,
                pedigree = t.toy,
                family = c("gaussian","gaussian"), 
                data = d.toy, prior=p.m.2, 
                nitt=210000, burnin=10000, thin=200,
                pr = TRUE, verbose = FALSE)

## BRMS
p.b.2 <- ()
b.2 <- brm(mvbind(y1, y2) ~ (1|a|gr(animal, cov = A)),
            data = d.toy,
            prior = p.b.2,
            family = gaussian(), 
            data2 = list(A = A.mat),
            cores = 4, chains = 4, iter = 6000, thin = 3)

Model Syntax

In MCMCglmm, we use cbind() to specify multiple response variables. The reserved term trait is used to specify fixed effects for all responses onm the LHS. In our example, trait-1 suppresses the global intercept, instead fitting separate intercepts for each trait \(y1\) and \(y2\). With multiple response variables, it is now possible to consider covariance between responses at different levels of the model hierarchy. MCMCglmm offers several options for defining the structure of covariance matrices to be estimated (see Hadfield XXX). The argument ~us(trait) specifies an unstructured (all elements estimated) covariance matrix with dimension equal to the number of response traits considered (here 2 x 2). The suffix :animal specifies this covariance matrix is on the level of the grouping factor animal, where-as :units encodes a residual covariance matrix.

In brms, we use mvbind() to specify multiple response variables. An index (here |a|, but any unique identifier is accepted) is used within the random effect specification to instruct brms to estimate the correlation between \(y1\) and \(y2\) at the animal group level. If omitted, only the correlation between \(y1\) and \(y2\) on the residual level will be estimated.

Inference

brms reports group level (here, phylogenetic) and residual correlations in the summary output. For MCMCglmm, correlations can be calculated form the reported (co)variances. In general, estimates of the correlation between traits can be calculated as

\[\begin{eqnarray} \rho_{12} = \frac{\Sigma_{12}}{\sqrt{(\Sigma_{11} \times \Sigma_{22})}} \end{eqnarray}\]

Applying the above calculation to the phylogenetic and residual (co)variance estimates from our MCMCglmm model fit yields phylogenetic and residual correlations, respectively.

## MCMCglmm
# phylogenetic correlation between y1 and y2
mean(m.2$VCV[,"traity1:traity2.animal"]/sqrt(m.2$VCV[,"traity1:traity1.animal"]*m.2$VCV[,"traity2:traity2.animal"]))
# residual correlation between y1 and y2
mean(m.2$VCV[,"traity1:traity2.units"]/sqrt(m.2$VCV[,"traity1:traity1.units"]*m.2$VCV[,"traity2:traity2.units"]))

## brms
# phylogenetic correlation = cor(y1_Intercept,y2_Intercept)
# residual correlation = rescor(y1,y2)
summary(b.2)

We can see close agreement between the estimates from MCMCglmm and brms (i.e., within Monte Carlo error), with both packages effectivly able to recover the generating simulation parameters.

Model 3 - Multivariate Non-Gaussian (Gaussian, Bernoulli)

Model Explanation

In previous simulations, we saw that Gaussian responses can be drawn directly from the MV structure, i.e., modeled via the identity link function. For non-Gaussian variables, (co)variances must be modeled on the link scale, e.g., for Poisson regression the mean is modeled with the log link while for binomial data, the probability of success is modeled with the logit link. To demonstrate this, we will again simulate and fit data under a bivariate model, this time taking a binomial response for \(y2\).

\[\begin{eqnarray} \begin{pmatrix}\mathbf{y_1} \\ \mathrm{logit}(\mathbf{p_2}) \end{pmatrix} &=& \begin{pmatrix}\boldsymbol{\mu}_1 + \mathbf{b}_1 + \mathbf{e}_1 \\ \boldsymbol{\mu}_2 + \mathbf{b}_2 + \mathbf{e}_2 \\ \end{pmatrix}\\[2mm] \\ \mathbf{y_2} &\sim& \mathrm{binomial}(\mathbf{p_2}) \end{eqnarray}\]

Data Simulation and Model Specification

Using our vectors of random effects and residuals simulated previously (on the link scale), we can now simply use the inverse link function to realise y2 as a binomial variable.

## SIMULATE DATA

# construct response traits from each linear predictor
y1 = u[1] + a1 + e1 # gaussian
y2 = rbinom(n,1,plogis(u[2] + a2 + e2)) # binomial

# generate df
species <- t.toy$tip.label
d.toy <- data.frame(species,y1,y2)
d.toy$animal <- d.toy$species
d.toy$obs <- 1:nrow(d.toy)

Specification of this model changes very little for MCMCglmm; the second argument to family is simply specified as “categorical.” Prior specification requires more careful attention. For binomial variables, the residual variance is not identifiable (REF). This is handled in MCMCglmm by fixing the residual variance of the binomial response to a nominal value (e.g., V = 1) in the prior specification. Fixing at higher values of V can improve mixing of the chain, but may also lead to numerical problems (see Hadfeild XXX). While it is theoretically possible to estimate residual correlations with a fixed residual variance (i.e., fixing the width of the error variance for one trait does not prevent correlation between joint multivariate draws), we point the user toward some practical concerns identified in the literature (REF).

For brms, it is useful to specify separate formulae for each response variable when considering different error families. Unlike MCMCglmm, brms does not model additive overdispersion by default for non-gaussian traits. Thus, in order to model residual (co)variances, it is necessary to specify an additive overdispersion term in the form an observation level random effect, (1|q|obs). However, this introduces a non-identifiability issue for the Gaussian error term, which cannot be silenced with default brms coding. Here, for simplicity, we have just constrained the Gaussian error term to be small (sigma = 0.1). A more appropriate (albeit technically demanding) solution is to edit the underlying stan code to prevent estimation of the redundant Gaussian error term (shown here - link).

# MCMCglmm

# NOTE - NEED TO USE DIFFERENT PRIORS FOR NON-GAUSSIAN RESPONSES
p3 <- list(G = list(G1 = list(V = diag(4), nu = 3.002)),
           R = list(V = diag(4), nu = 3.002))

m.3 <- MCMCglmm(cbind(y1, y2) ~ trait-1,
                random = ~us(trait):animal,
                rcov = ~us(trait):units,
                pedigree=t.toy,
                family = c("gaussian", "categorical"), 
                data = d.toy, prior=p3, 
                nitt=210000, burnin=10000, thin=200,
                pr=TRUE,verbose = FALSE)

# BRMS
bf_y1 <- bf(y1 ~ 1 + (1|p|gr(animal, cov = A)) + (1|q|obs), sigma = 0.1) + gaussian()
bf_y2 <- bf(y2 ~ 1 + (1|p|gr(animal, cov = A)) + (1|q|obs)) + bernoulli()

b.3 <- brm(bf_y1 + bf_y2 + set_rescor(FALSE),
          data = d.toy, 
          family = gaussian(), 
          data2 = list(A = A.mat),
          cores=4,
          chains=4, iter = 2000, thin = 1
)

Inference

Because the residual variance for binomial traits is fixed in MCMCglmm, it may be desirable to re-scale the phylogenetic correlations relative to the magnitude of this fixed residual variance (see Hadfield XXXX).

Another consequence of fixing the residual variance is that \(h^2\) is no longer an appropriate estimate of phylogenetic signal. For our example, the amount of variation in our binomial trait \(y2\) explained by shared ancestry may instead be estimated as the intraclass correlation coefficient

\[\begin{eqnarray} ICC = \Sigma^{\mathrm{b}}_{11}/(\Sigma^{\mathrm{b}}_{11}+V+\pi^{2/3}) \end{eqnarray}\]

where \(V\) is the level of residual variance (additive over-dispersion) fixed in the prior specification (here, V = 1) and \(\pi^{2/3}\) is the distribution specific variance term for the binomial (see Hadfield XXX; Nakagawa and Schielzeth XXX)

Model Validation

Blah blah

NOTES

For a reference use at symbol Royle (2004) , see the bib_phmm.bib file for standard bibtex.

References

Royle, J. Andrew. 2004. “N-Mixture Models for Estimating Population Size from Spatially Replicated Counts.” Biometrics 60 (1): 108–15. https://doi.org/10.1111/j.0006-341X.2004.00142.x.