Sections: Fundamentals / Nucleotide model selection / Alternatives to LRT / Testing phylogenetic trees / Model adequacy / Automatic model selection
You can download the files required for the practical from the wiki.
Set your working directory in R to the place you extracted the files for the practicals e.g.
setwd("~/Documents/2024/come/hypothesis_testing/hyptest-practical")Fair or biased coin example. Imagine we have a result of an experiment: the number of heads when tossing a coin 100 times.
number_of_coin_tosses <- 100
number_of_observed_heads <- 65If the coin is fair, the probability of getting heads is 1/2. This is our null model. We can calculate the log-likelihood of the observed data.
prob_of_heads_fair_coin <- 0.5
ln_l_if_fair <- log( dbinom(number_of_observed_heads, number_of_coin_tosses, prob_of_heads_fair_coin) )
ln_l_if_fairAlternatively, if the coin is biased towards heads, the probability of getting heads that maximises the likelihood is number of heads/tosses. This is the alternative model. We can calculate the likelihood of the observed data.
prob_of_heads_biased_coin <- number_of_observed_heads / number_of_coin_tosses
ln_l_if_biased <- log( dbinom(number_of_observed_heads, number_of_coin_tosses, prob_of_heads_biased_coin) )
ln_l_if_biasedNext, we calculate the value of our test statistic for the observed data.
likelihood_ratio <- ln_l_if_biased - ln_l_if_fair
likelihood_ratioIs this larger than expected? We can determine the distribution of the test statistic under the null model using simulation. We choose our level of significance to reject the null, and simulate the experiment many times.
alpha <- 0.05
number_of_null_experiments <- 1000000First, we run a “fair coin toss” (the null model) experiment. Each item in the array indicates the number of heads (we show the first few).
heads_from_null <- rbinom(number_of_null_experiments, number_of_coin_tosses, prob_of_heads_fair_coin)
head(heads_from_null)Next, calculate the log-likelihood of each simulation result as if coin were fair (null model)
ln_l_fair_from_null <- log( dbinom(heads_from_null, number_of_coin_tosses, prob_of_heads_fair_coin) )
head(ln_l_fair_from_null)And calculate the log-likelihood of each simulation result as if coin were biased (alternative model)
ml_prob_for_biased <- heads_from_null / number_of_coin_tosses
ln_l_biased_from_null <- log( dbinom(heads_from_null, number_of_coin_tosses, ml_prob_for_biased) )
head(ln_l_biased_from_null)Finally, we calculate the likelihood ratio for each of the simulations. The plot shows the distribution of the test statistics under the null model.
likelihood_ratio_from_null <- ln_l_biased_from_null - ln_l_fair_from_null
hist(likelihood_ratio_from_null, breaks=100)From this distribution, we can determine the test statistic cutoff (aka critical value) at the chosen significance level.
likelihood_ratio_from_null <- sort(likelihood_ratio_from_null)
t_threshold <- likelihood_ratio_from_null[number_of_null_experiments - (alpha * number_of_null_experiments)]
t_thresholdWe plot the distribution of the test statistic and show the threshold on the histogram as a vertical blue line. We see where the test statistic of the original, observed, result (added as red line). If the result it more extreme (larger) than the critical value, we reject the null model.
hist(likelihood_ratio_from_null, breaks=100)
abline(v=t_threshold, col='blue')
abline(v=likelihood_ratio, col='red')Load the ape and phangorn libraries.
library(ape)
library(phangorn)Load the nucleotide alignment.
brown_aln <- read.phyDat(file = 'brown.fasta', format = 'fasta', type='DNA')
brown_alnImagine we already have the topology. Load the unrooted tree.
example_tree <- "((Human:0.1, Chimpanzee:0.1):0.1, Gorilla:0.1, (Orangutan:0.1, Gibbon:0.1):0.1);"
brown_tree <- read.tree(text=example_tree)
brown_tree
plot(brown_tree, type='unrooted') # have a look at the tree
write.tree(brown_tree) # to see the newick-formatted treeSay we want to see whether our data supports a model with independent substitution rates. The null model is HKY (ts/tv ratio, base frequencies) & the alternative model is GTR (indepdendent substitution rates, base frequencies). Our significance threshold is set to 0.05.
We use two functions from the Phangorn library. pml
creates an object that calculates the likelihood of the data.
optim.pml takes that object and optimises the parameters of
the model (including branch lengths). First we optimise the model under
the HKY model.
alpha_signif <- 0.05
brown_pml <- pml(brown_tree, data = brown_aln)
hky_ml <- optim.pml(brown_pml, model='HKY')You can inspect the results and also get the number of parameters in the model (the number includes branch lengths).
hky_ml # substitution model maximum-likelihood, rate matrix, base frequencies
cat(write.tree(hky_ml$tree)) # tree with branch lengths
hky_ml$df # number of parameters in the model + branch lengthsThen we optimise the model under the GTR model.
gtr_ml <- optim.pml(brown_pml, model='GTR')
gtr_mlWe use the likelihood ratio test to determine whether there is evidence to reject the null hypothesis at significance of 0.05. The maximum log-likelihood of the optimised model is in the ‘logLik’ attribute.
hky_ml$logLik
gtr_ml$logLikNote, you can see inspect the maximum likelihood estimated tree with branch lengths using:
cat(write.tree(gtr_ml$tree))Question Calculate the likelihood ratio test
statistic:
2 * (alternative log-likelihood - null log-likelihood)
Question HKY is nested within the GTR model. How
many degrees of freedom should we use for the chi-squared distribution?
Hint: the total number of parameters for the HKY model is available
using hky_ml$df.
Question Calculate the p-value of the test statistic
using
pchisq(TEST_STATISTIC, df=DEGREES_OF_FREEDOM, lower.tail=F)
Question Do we reject the null hypothesis?
Lets do the same to see if data supports site-rate variation. We need
to add two arguments. k specifies the number of rate
categories for the gamma distribution. optGamma will
optimise the shape parameter of the gamma distribution.
gtrg_ml <- optim.pml(pml(brown_tree, data = brown_aln, k=4), model='GTR', optGamma=TRUE)Question Calculate the likelihood ratio. Is it greater than the critical value at 0.05 significance? Do we reject the null?
Extra credit If the null model is HKY+G, do we reject when alternative is GTR+G?
We want to test the molecular clock hypothesis. The null hypothesis specified the molecular clock i.e. root-to-tip distance is the same for all tips. The alternative hypothesis rejects the molecular clock. We fix the nucleotide model, and only enforce the constraint on branch lengths. This example uses GTR+G model. Load the rooted example tree.
tree.rooted <- read.tree(text='((((Human:0.1,Chimpanzee:0.1):0.1,Gorilla:0.2):0.1,Orangutan:0.3):0.1,Gibbon:0.4);')
tree.rooted
plot(tree.rooted)Now, optimise the branch lengths under the GTR+G gamma with molecular clock.
rooted_pml <- pml(tree.rooted, brown_aln, k=4)
gtrg_clock_ml <- optim.pml(rooted_pml, model='GTR', optGamma = TRUE, optRooted = TRUE)We already have the maximum likelihood estimation for the unrooted tree from the question above:
gtrg_mlQuestion Plot the maximum likelihood estimated tree for the data under the clock and GTR+G.
plot(gtrg_clock_ml$tree)Question If GTR+G (rooted & constrained tree) is the null model, and GTR+G (unrooted & unconstrained tree), what is the value of the test statistic?
Question How many degrees of freedom for the chi-squared table?
Question Do we reject the null hypothesis?
The Brown dataset is a concatenated alignment of protein coding sequences, with tRNA sequence in the middle. Given this information, say we decide to use a partitioned model to analyse the data. First, let’s build a table to specify the four partitions (first, second and third codon positions, and tRNA) for the alignment.
g1_parts <- c(rep(1:3, length.out=457))
trna_parts <- rep(4, length.out=198)
g2_parts <- rep(1:3, length.out=240)
partition_table <- table(attr(brown_aln, 'index'), c(g1_parts, trna_parts, g2_parts))Given the partitions, we supply a GTR model for each parti
gtr_ml <- optim.pml(pml(brown_tree, data = brown_aln), model='GTR')
gtr_partitioned <- pmlPart(edge ~ rate + Q + bf, gtr_ml, weight=partition_table)Question What is the log-likelihood of the GTR-partitioned model?
gtr_partitioned$logLik[1]We want to compare GTR+G and GTR-partitioned models, and use the one most suitable. They are not nested. Calculate the AIC using the formula below.
\[AIC=2k - 2\text{log}(L)\]
We can get the total number of parameters of each model:
df_gtrg <- gtrg_ml$df
df_gtr_partitioned <- attr(gtr_partitioned$logLik, "df")Question What is the AIC under the GTR+G model? What is the AIC under the GTR+partitioned model? Which model should we prefer?
aic_gtrg <- 2 * (df_gtrg) - 2 * (gtrg_ml$logLik)
aic_gtrg
aic_gtr_partitioned <- 2 * (df_gtr_partitioned) - 2 * (gtr_partitioned$logLik[1])
aic_gtr_partitionedWe will use the programs baseml and evolver
from the PAML
suite.
First, we analyse the Brown dataset using baseml. We
require three files:
brown.nuc - the alignment file (PHYLIP format)5s.unrooted.trees - the 15 possible unrooted trees for
five taxabaseml.ctl - the baseml control file, providing
arguments for running the programNote: some of the commands below should be run in the terminal, not in R (e.g. executing baseml, evolver and parsing the results file). These blocks have been identified with “# terminal”
Examine the control file. Any text after the asterisk are comments.
# terminal
cat baseml.ctlWe run baseml in the directory with our data files. We
analyse all trees under the GTR+G model and determine which is the
maximum likelihood tree. You’ll see this output on the console.
# terminal
baseml baseml.ctlThe results are written to a file results.out. Open it,
and scroll to the bottom of the file. There is a table with the
log-likelihood for each of the trees. The maximum likelihood tree is
marked with an asterisk (it should be tree 7). You can scroll up and get
the information about the analysis using tree 7.
Imagine we, based on some other information, believe that Tree 10 is in fact the true tree i.e. that despite tree 10 being the true tree, tree 7 is the maximum likelihood tree. But, in fact, the trees are both as good as each other.
We can use parametric bootstrap to test this hypothesis, similar to our procedure for testing the coin toss experiment. We simulate data under the null hypothesis (that tree 10 and the maximum likelihood tree are just as good as each other). Have a look at the model parameters for tree 10:
# terminal
grep -A 25 "TREE # 10" results.outWe examine the results file to calculate our test statistic delta, the log-likelihood difference between the maximum likelihood tree and tree 10.
Question What is the value of delta? Use the table at the end of the results file.
d <- -2616.074 - -2619.738We can use evolver to simulate data. We edit the
MCbase.dat file and specify the estimates for parameters
for Tree 10 analysis. Specifically, you need to change:
You file should look something like this:
# terminal
cat MCbase_Tree10.datWe can now run evolver to create replicate datasets
under Tree 10. The program will create a file mc.paml that
will contain 500 simulated alignments.
# terminal
evolver 5 MCbase_Tree10.datmc.paml contains 500 alignments that were simulated on
Tree 10. We now analyse all these data sets using the same process we
analysed the original data, running baseml with the same
options. Make a copy of the existing baseml.ctl control
file and modify it to analyse the tree 10 simulated data. Make the
following changes:
# terminal
cat baseml_tree10.ctl[Run this command in your terminal - it will take some time to
complete, and might crash RStudio!] Now run baseml giving
the new control file name as program option.
# terminal
baseml baseml_tree10.ctlThe results of analysing the simulated data are in
tree10.sim.out. We need to parse out all the results for
tree 10, specifically the line for tree 10 in the table of results at
the end of each dataset analysis. We can achieve this using grep.
# terminal
grep "^ 10" tree10.sim.out > tree10_loglik.txt
head tree10_loglik.txtEach line is the result of analysing a single dataset under tree 10. Note, sometimes tree 10 is the maximum likelihood tree (marked with asterisk) and other times some other tree had the maximum likelihood (no asterisk). The third column of the table is delta, the difference between the maximum likelihood tree and the tree 10 likelihood.
Read the table into R and look at the distribution of delta.
tree10_results <- read.table('tree10_loglik.txt', sep="")
tree10_deltas <- abs(tree10_results$V3)
hist(tree10_deltas)Finally, we can calculate the p-value of delta of the real data (3.664).
p = sum(d < tree10_deltas) / length(tree10_deltas)
pQuestion Do we reject the null hypothesis and say that tree 10 is not as good a tree as the maximum likelihood tree (tree 7)?
We can use the Goldman-Cox test to check for model adequacy.
We take the maximum likelihood estimates for the real data: tree 7 and the MLE parameters of the GTR+G model. The null hypothesis is that the GTR+G model could plausibly have generated the observed data. The alternative hypothesis is that an unconstrained multinomial model generated the data. The test statistic is the difference between the log-likelihoods of the two models.
In order to perform the test, we follow the same procedure as for the tests for topology.
Inspect the results.out file. The table at the end
indicates that tree 7 is the maximum likelihood tree (has an asterisk).
Note the maximum likelihood estimates for tree 7:
# terminal
grep -A 25 "TREE # 7" results.outAlso note the unconstrained log-likelihood:
# terminal
grep "unconstrained" results.outWe can calculate the value of the test statistic:
-2476.965133 - -2616.073763 = 139.1
d <- -2476.965133 - -2616.073763Use evolver to generate data sets simulated on tree 7
using the maximum likelihood estimated. You need to edit
MCbase.dat with the estimated tree, rates, alpha and base
frequencies parameters.
# terminal
cat MCbase_Tree7.datRun evolver as before.
# terminal
evolver 5 MCbase_Tree7.datThe simulated data sets are in mc.paml. We create a new
baseml control file to analyse the replicates. Note, I have
put tree 7 into it’s own file to make it faster (baseml will not analyse
all the possible topologies).
# terminal
cat baseml_tree7.ctlRun baseml:
# terminal
baseml baseml_tree7.ctlWe need to extract two numbers for each analysis: the maximum log-likelihood under the specified nucleotide model and the unconstrained log-likelihood.
You’ll need some of your UNIX skills to extract the numbers from the results file. Using grep:
# terminal
grep "^lnL" tree7.sim.out > tree7.sim.lnl.txt
grep "unconstrained" tree7.sim.out > tree7.sim.unconstrained.txtRead these files into R and calculate delta, the GTR+G log-likelihood subtracted from the unconstrained log-likelihood.
tree7_lnl <- read.table('tree7.sim.lnl.txt', sep='')$V5
tree7_unconstrained <- read.table('tree7.sim.unconstrained.txt', sep='')$V5
tree7_delta <- tree7_unconstrained - tree7_lnl
hist(tree7_delta)
abline(v=d)Calculate the p-value of the delta from the real data using the distribution of the test statistic from data simulated under the null.
p <- sum(d < tree7_delta) / length(tree7_delta)
p < alpha_signifQuestion What does this mean?
mT <- modelTest(brown_aln, brown_tree)
mT <- mT[order(mT$AIC),] # sort by AIC
mT