David Wells *12/07/2020*

Simiulate data for a phylogenetic generalised linear mixed model (GLMM) and fit that model to assess accuracy of estiamtes. A phylogenetic GLMM is a comparative method which accounts for similarity between related species. It is very similar to an "animal model" which uses a pedigree rather than a phylogeny to account for similarity between related individuals.

The phylogenetic GLMM fits a random effect which takes a separate value for each species but which follows a correlation structure based on the phylogeny. I.e. the random intercept is more closely correlated for more closely related species.

A lot of this borrows from chapter 6 in the MCMCglmm course notes (which is incomplete).

In [28]:

```
packages <- c("MCMCglmm", "dplyr", "ggplot2", "ggtree", "diosR", "ape", "MASS")
lapply(packages, require, character.only=T)
sessionInfo()
options(repr.plot.width=6, repr.plot.height=4)
```

Simulate a random tree. The structure of the tree can be expressed as a relatedness matrix $A$. $A$ is a symetric matrix where each element $A_{ij}$ represents the branch length from root to split of species $i$ and $j$. This means that the diagonal is the distance of a taxa from the root. You can verrify this by comparing the table and tree below. However, this matrix is often scaled so that the distance from root to tip is 1. This is only possible for ultrametric trees so in this case we do not scale it. Plotted next to the tree in colour it is easy to see that this $A$ matrix ensure high covariance amongst closely related species.

In [105]:

```
n <- 5
tree<- rtree(n)
plot(tree)
edgelabels(round(tree$edge.length,3), bg="white", frame="none", col="dimgrey")
```

In [106]:

```
Aphy <- vcv.phylo(tree, cor=F)
round(Aphy, 2)
```

In [107]:

```
gheatmap(
ggtree(tree) + geom_tiplab() + ggtitle("Tree and Aphy"),
Aphy
) +
scale_fill_gradientn(colours=c("#185d75", "#49ebbe"))
# scale_fill_gradientn(colours=c("darkturquoise", "orange"))
```

For our actual analysis we simulate a much larger tree. Then we use this $A$ matrix `Aphy`

to enforce the desired corelation structure on to the random intercepts. The variance of this random effect is specified as `phyV`

and the residual variance (of the overall model) as `resV`

. The product of these variables gives the covariance structure that we sample the phylogenetic random effects `phyfx`

from. The dataset `sim`

is simulated with x as a fixed effect (for simplicity we use `x`

as is which is equivalent to a coefficient of 1) and the random intercept drawn from `phyfx`

.

In [4]:

```
n <- 1000
phyV <- 7
resV <- 2
tree<- rtree(n)
Aphy <- vcv.phylo(tree, cor=F)
phyfx <- mvrnorm(1, mu=rep(0, n), Sigma=Aphy * phyV)
```

In [5]:

```
sim <- data.frame(
x = runif(n),
phylo = paste0("t",1:n),
y = NA
)
sim$y <- rnorm(n, sim$x + phyfx[as.character(sim$phylo)], sd=sqrt(resV))
```

Next we set the prior for the random effect `G`

and the residual variance `R`

. Then we're ready to run the model, again setting `scale = F`

to allow a non-ultrametric tree. The random effect is fit as `phylo`

which indicates the column of the dataset which matches rows up with nodes in the tree. We leave the thinning interval and number of iterations as defaults here but strongly recommend that you assess and adjust your model to ensure good mixing and convergence.

Finally we plot the distribution of the estimated phylogenetic variance and the residual variance with the true value used for simulations in red.

In [6]:

```
prior1 <- list(
G=list(
G1=list(V=1, nu=0.002)
),
R=list(V=1, nu=0.002)
)
m1 <- MCMCglmm(y~x,
random=~phylo,
ginverse=list(phylo=inverseA(tree, scale=F)$Ainv),
data=sim,
prior=prior1,
verbose=F)
```

In [7]:

```
densplot(m1$VCV[,'phylo'], main="Phylogenetic variance")
abline(v=phyV, col="red")
densplot(m1$VCV[,'units'], main="Residual variance")
abline(v=resV, col="red")
```