In USCbiostats/phylogenetic: Statistical Inference and Prediction of Annotations in Phylogenetic Trees

knitr::opts_chunk$set(fig.width = 6, fig.height = 6, fig.align = "center")

This is a simple example of simulation, inference, and prediction with the aphylo R package.

Setup

library(aphylo)

# Parameter estimates
psi  <- c(.05, .025)
mu_d <- c(.2, .1)
mu_s <- c(.1, .05)
Pi   <- .5

Simulation

set.seed(223)
x <- raphylo(n = 200, psi = psi, mu_d = mu_d, mu_s = mu_s, Pi = Pi)
plot(x)

The simulation function generates an aphylo class object which is simply a wrapper containing:

• a 0/1 integer matrix (annotations),
• a phylo tree (from the ape package), and
• some information about the tree and annotations.

If needed, we can export the data as follows:

# Edgelist describing parent->offspring relations
write.csv(x$tree, file = "tree.tree", row.names = FALSE)

# Tip annotations
ann <- with(x, rbind(tip.annotation, node.annotation))
write.csv(ann, file = "annotations.csv", row.names = FALSE)

# Event types
events <- with(x, cbind(c(tip.type*NA, node.type)))
rownames(events) <- 1:nrow(events)
write.csv(events, file = "events.csv", row.names = FALSE)

Inference

To fit the data, we can use MCMC as follows:

ans <- aphylo_mcmc(x ~ psi + mu_d + mu_s + Pi)
ans

For goodness-of-fit analysis, we have a couple of tools. We can compare the predicted values with the observed values:

plot(ans)

We can also take a look at the surface of the posterior function

plot_logLik(ans)

And we can also take a look at the prediction scores

ps <- prediction_score(ans)
ps       # Printing
plot(ps) # and plotting

USCbiostats/phylogenetic documentation built on Oct. 28, 2023, 7:23 a.m.