Hypothesis Testing

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")

Fundamentals

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 <- 65

If 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_fair

## [1] -7.054105

Alternatively, 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_biased

## [1] -2.484051

Next, we calculate the value of our test statistic for the observed data.

likelihood_ratio <- ln_l_if_biased - ln_l_if_fair
likelihood_ratio

## [1] 4.570054

Is 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 <- 1000000

First, 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)

## [1] 49 52 49 50 50 57

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)

## [1] -2.550679 -2.610102 -2.550679 -2.530876 -2.530876 -3.504272

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)

## [1] -2.530678 -2.530081 -2.530678 -2.530876 -2.530876 -2.521046

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_threshold

## [1] 2.013551

We 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')

Nucleotide model selection

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_aln

## 5 sequences with 895 character and 85 different site patterns.
## The states are a c g t

Imagine 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

## 
## Phylogenetic tree with 5 tips and 3 internal nodes.
## 
## Tip labels:
##   Human, Chimpanzee, Gorilla, Orangutan, Gibbon
## 
## Unrooted; includes branch lengths.

plot(brown_tree, type='unrooted')  # have a look at the tree

write.tree(brown_tree)  # to see the newick-formatted tree

## [1] "((Human:0.1,Chimpanzee:0.1):0.1,Gorilla:0.1,(Orangutan:0.1,Gibbon:0.1):0.1);"

Say 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')

## optimize edge weights:  -2933.577 --> -2843.877 
## optimize base frequencies:  -2843.877 --> -2838.466 
## optimize rate matrix:  -2838.466 --> -2669.197 
## optimize edge weights:  -2669.197 --> -2668.163 
## optimize base frequencies:  -2668.163 --> -2665.072 
## optimize rate matrix:  -2665.072 --> -2664.96 
## optimize edge weights:  -2664.96 --> -2664.93 
## optimize base frequencies:  -2664.93 --> -2664.926 
## optimize rate matrix:  -2664.926 --> -2664.926 
## optimize edge weights:  -2664.926 --> -2664.926 
## optimize base frequencies:  -2664.926 --> -2664.926 
## optimize rate matrix:  -2664.926 --> -2664.926 
## optimize edge weights:  -2664.926 --> -2664.926

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

## model: HKY 
## loglikelihood: -2664.926 
## unconstrained loglikelihood: -2476.965 
## 
## Rate matrix:
##         a       c       g       t
## a 0.00000 1.00000 9.30606 1.00000
## c 1.00000 0.00000 1.00000 9.30606
## g 9.30606 1.00000 0.00000 1.00000
## t 1.00000 9.30606 1.00000 0.00000
## 
## Base frequencies:  
##         a         c         g         t 
## 0.3059668 0.3233061 0.1130702 0.2576569

cat(write.tree(hky_ml$tree))  # tree with branch lengths

## ((Human:0.04156693165,Chimpanzee:0.05401620696):0.01741372836,Gorilla:0.05806931237,(Orangutan:0.1003618749,Gibbon:0.1394886842):0.05303070765);

hky_ml$df  # number of parameters in the model + branch lengths

## [1] 11

Then we optimise the model under the GTR model.

gtr_ml <- optim.pml(brown_pml, model='GTR')

## optimize edge weights:  -2933.577 --> -2843.877 
## optimize base frequencies:  -2843.877 --> -2838.466 
## optimize rate matrix:  -2838.466 --> -2660.085 
## optimize edge weights:  -2660.085 --> -2659.196 
## optimize base frequencies:  -2659.196 --> -2657.351 
## optimize rate matrix:  -2657.351 --> -2657.105 
## optimize edge weights:  -2657.105 --> -2657.099 
## optimize base frequencies:  -2657.099 --> -2657.062 
## optimize rate matrix:  -2657.062 --> -2657.055 
## optimize edge weights:  -2657.055 --> -2657.055 
## optimize base frequencies:  -2657.055 --> -2657.053 
## optimize rate matrix:  -2657.053 --> -2657.053 
## optimize edge weights:  -2657.053 --> -2657.053 
## optimize base frequencies:  -2657.053 --> -2657.053 
## optimize rate matrix:  -2657.053 --> -2657.053 
## optimize edge weights:  -2657.053 --> -2657.053 
## optimize base frequencies:  -2657.053 --> -2657.053 
## optimize rate matrix:  -2657.053 --> -2657.053 
## optimize edge weights:  -2657.053 --> -2657.053

gtr_ml

## model: GTR 
## loglikelihood: -2657.053 
## unconstrained loglikelihood: -2476.965 
## 
## Rate matrix:
##           a         c         g         t
## a   0.00000  74.52946 461.57211  39.71261
## c  74.52946   0.00000  36.24128 453.33820
## g 461.57211  36.24128   0.00000   1.00000
## t  39.71261 453.33820   1.00000   0.00000
## 
## Base frequencies:  
##         a         c         g         t 
## 0.3016210 0.3197885 0.1164201 0.2621704

We 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

## [1] -2664.926

gtr_ml$logLik

## [1] -2657.053

Note, you can see inspect the maximum likelihood estimated tree with branch lengths using:

cat(write.tree(gtr_ml$tree))

## ((Human:0.04200909362,Chimpanzee:0.05352439774):0.01634363539,Gorilla:0.05875708429,(Orangutan:0.1001617194,Gibbon:0.1397753847):0.05265780496);

Part 1 (HKY vs GTR)

Question Calculate the likelihood ratio test statistic: 2 * (alternative log-likelihood - null log-likelihood)

## [1] 15.746

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.

## [1] 4

Question Calculate the p-value of the test statistic using pchisq(TEST_STATISTIC, df=DEGREES_OF_FREEDOM, lower.tail=F)

## [1] 0.00337963

Question Do we reject the null hypothesis?

## [1] TRUE

Part 2 (GTR vs GTR+G)

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)

## optimize edge weights:  -2893.676 --> -2828.074 
## optimize base frequencies:  -2828.074 --> -2823.131 
## optimize rate matrix:  -2823.131 --> -2638.025 
## optimize shape parameter:  -2638.025 --> -2634.99 
## optimize edge weights:  -2634.99 --> -2629.967 
## optimize base frequencies:  -2629.967 --> -2626.06 
## optimize rate matrix:  -2626.06 --> -2624.129 
## optimize shape parameter:  -2624.129 --> -2622.388 
## optimize edge weights:  -2622.388 --> -2621.091 
## optimize base frequencies:  -2621.091 --> -2620.789 
## optimize rate matrix:  -2620.789 --> -2620.068 
## optimize shape parameter:  -2620.068 --> -2619.494 
## optimize edge weights:  -2619.494 --> -2618.898 
## optimize base frequencies:  -2618.898 --> -2618.802 
## optimize rate matrix:  -2618.802 --> -2618.443 
## optimize shape parameter:  -2618.443 --> -2618.175 
## optimize edge weights:  -2618.175 --> -2617.855 
## optimize base frequencies:  -2617.855 --> -2617.815 
## optimize rate matrix:  -2617.815 --> -2617.608 
## optimize shape parameter:  -2617.608 --> -2617.457 
## optimize edge weights:  -2617.457 --> -2617.267 
## optimize base frequencies:  -2617.267 --> -2617.25 
## optimize rate matrix:  -2617.25 --> -2617.123 
## optimize shape parameter:  -2617.123 --> -2617.03 
## optimize edge weights:  -2617.03 --> -2616.899 
## optimize base frequencies:  -2616.899 --> -2616.891 
## optimize rate matrix:  -2616.891 --> -2616.799 
## optimize shape parameter:  -2616.799 --> -2616.746 
## optimize edge weights:  -2616.746 --> -2616.654 
## optimize base frequencies:  -2616.654 --> -2616.65 
## optimize rate matrix:  -2616.65 --> -2616.585 
## optimize shape parameter:  -2616.585 --> -2616.552 
## optimize edge weights:  -2616.552 --> -2616.492 
## optimize base frequencies:  -2616.492 --> -2616.489 
## optimize rate matrix:  -2616.489 --> -2616.446 
## optimize shape parameter:  -2616.446 --> -2616.422 
## optimize edge weights:  -2616.422 --> -2616.377 
## optimize base frequencies:  -2616.377 --> -2616.376 
## optimize rate matrix:  -2616.376 --> -2616.342 
## optimize shape parameter:  -2616.342 --> -2616.328 
## optimize edge weights:  -2616.328 --> -2616.296 
## optimize base frequencies:  -2616.296 --> -2616.295 
## optimize rate matrix:  -2616.295 --> -2616.271 
## optimize shape parameter:  -2616.271 --> -2616.261 
## optimize edge weights:  -2616.261 --> -2616.239

Question Calculate the likelihood ratio. Is it greater than the critical value at 0.05 significance? Do we reject the null?

## [1] 1.642616e-19

Extra credit If the null model is HKY+G, do we reject when alternative is GTR+G?

Part 3 (Molecular clock)

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

## 
## Phylogenetic tree with 5 tips and 4 internal nodes.
## 
## Tip labels:
##   Human, Chimpanzee, Gorilla, Orangutan, Gibbon
## 
## Rooted; includes branch lengths.

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)

## optimize edge weights:  -2982.882 --> -2828.954 
## optimize base frequencies:  -2828.954 --> -2823.991 
## optimize rate matrix:  -2823.991 --> -2638.367 
## optimize shape parameter:  -2638.367 --> -2635.145 
## optimize edge weights:  -2635.145 --> -2630.203 
## optimize base frequencies:  -2630.203 --> -2626.28 
## optimize rate matrix:  -2626.28 --> -2624.299 
## optimize shape parameter:  -2624.299 --> -2622.679 
## optimize edge weights:  -2622.679 --> -2621.524 
## optimize base frequencies:  -2621.524 --> -2621.227 
## optimize rate matrix:  -2621.227 --> -2620.586 
## optimize shape parameter:  -2620.586 --> -2620.037 
## optimize edge weights:  -2620.037 --> -2619.545 
## optimize base frequencies:  -2619.545 --> -2619.457 
## optimize rate matrix:  -2619.457 --> -2619.159 
## optimize shape parameter:  -2619.159 --> -2618.905 
## optimize edge weights:  -2618.905 --> -2618.636 
## optimize base frequencies:  -2618.636 --> -2618.6 
## optimize rate matrix:  -2618.6 --> -2618.424 
## optimize shape parameter:  -2618.424 --> -2618.302 
## optimize edge weights:  -2618.302 --> -2618.154 
## optimize base frequencies:  -2618.154 --> -2618.138 
## optimize rate matrix:  -2618.138 --> -2618.041 
## optimize shape parameter:  -2618.041 --> -2617.963 
## optimize edge weights:  -2617.963 --> -2617.866 
## optimize base frequencies:  -2617.866 --> -2617.858 
## optimize rate matrix:  -2617.858 --> -2617.79 
## optimize shape parameter:  -2617.79 --> -2617.75 
## optimize edge weights:  -2617.75 --> -2617.686 
## optimize base frequencies:  -2617.686 --> -2617.682 
## optimize rate matrix:  -2617.682 --> -2617.636 
## optimize shape parameter:  -2617.636 --> -2617.612 
## optimize edge weights:  -2617.612 --> -2617.569 
## optimize base frequencies:  -2617.569 --> -2617.566 
## optimize rate matrix:  -2617.566 --> -2617.534 
## optimize shape parameter:  -2617.534 --> -2617.519 
## optimize edge weights:  -2617.519 --> -2617.489 
## optimize base frequencies:  -2617.489 --> -2617.488 
## optimize rate matrix:  -2617.488 --> -2617.466 
## optimize shape parameter:  -2617.466 --> -2617.458 
## optimize edge weights:  -2617.458 --> -2617.435 
## optimize base frequencies:  -2617.435 --> -2617.434 
## optimize rate matrix:  -2617.434 --> -2617.417 
## optimize shape parameter:  -2617.417 --> -2617.412 
## optimize edge weights:  -2617.412 --> -2617.397

We already have the maximum likelihood estimation for the unrooted tree from the question above:

gtrg_ml

## model: GTR+G(4) 
## loglikelihood: -2616.239 
## unconstrained loglikelihood: -2476.965 
## Discrete gamma model
## Number of rate categories: 4 
## Shape parameter: 0.2167888 
## 
## Rate matrix:
##           a          c          g         t
## a    0.0000  273.44358 4804.69523  147.8975
## c  273.4436    0.00000   88.01125 4045.8167
## g 4804.6952   88.01125    0.00000    1.0000
## t  147.8975 4045.81669    1.00000    0.0000
## 
## Base frequencies:  
##          a          c          g          t 
## 0.32149694 0.34442252 0.09696143 0.23711911

Question 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?

## [1] FALSE

Alternatives to LRT

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')

## optimize edge weights:  -2933.577 --> -2843.877 
## optimize base frequencies:  -2843.877 --> -2838.466 
## optimize rate matrix:  -2838.466 --> -2660.085 
## optimize edge weights:  -2660.085 --> -2659.196 
## optimize base frequencies:  -2659.196 --> -2657.351 
## optimize rate matrix:  -2657.351 --> -2657.105 
## optimize edge weights:  -2657.105 --> -2657.099 
## optimize base frequencies:  -2657.099 --> -2657.062 
## optimize rate matrix:  -2657.062 --> -2657.055 
## optimize edge weights:  -2657.055 --> -2657.055 
## optimize base frequencies:  -2657.055 --> -2657.053 
## optimize rate matrix:  -2657.053 --> -2657.053 
## optimize edge weights:  -2657.053 --> -2657.053 
## optimize base frequencies:  -2657.053 --> -2657.053 
## optimize rate matrix:  -2657.053 --> -2657.053 
## optimize edge weights:  -2657.053 --> -2657.053 
## optimize base frequencies:  -2657.053 --> -2657.053 
## optimize rate matrix:  -2657.053 --> -2657.053 
## optimize edge weights:  -2657.053 --> -2657.053

gtr_partitioned <- pmlPart(edge ~ rate + Q + bf, gtr_ml, weight=partition_table)

## loglik: -2657.053 --> -2452.759 
## loglik: -2452.759 --> -2452.092 
## loglik: -2452.092 --> -2452.011 
## loglik: -2452.011 --> -2451.996 
## loglik: -2451.996 --> -2451.993 
## loglik: -2451.993 --> -2451.991 
## loglik: -2451.991 --> -2451.991 
## loglik: -2451.991 --> -2451.991 
## loglik: -2451.991 --> -2451.991

Question What is the log-likelihood of the GTR-partitioned model?

gtr_partitioned$logLik[1]

## [1] -2451.991

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

## [1] 5264.477

aic_gtr_partitioned <- 2 * (df_gtr_partitioned) - 2 * (gtr_partitioned$logLik[1])
aic_gtr_partitioned

## [1] 4987.982

Testing phylogenetic trees

We will use the programs baseml and evolver from the PAML suite.

First, we analyse the Brown dataset using baseml. We require three files:

1. brown.nuc - the alignment file (PHYLIP format)
2. 5s.unrooted.trees - the 15 possible unrooted trees for five taxa
3. baseml.ctl - the baseml control file, providing arguments for running the program

Note: 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”

Analysing the real data

Examine the control file. Any text after the asterisk are comments.

# terminal
cat baseml.ctl

## seqfile = brown.nuc
## outfile = results.out
## treefile = 5s.unrooted.trees
## model = 7  * GTR
## fix_alpha = 0 * 0: estimate alpha
## alpha = 0.2
## nhomo = 1

We 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.ctl

The 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.out

## TREE # 10:  ((1, 3), 2, (4, 5));  MP score: 358.00
## lnL(ntime:  7  np: 16):  -2619.737562    +6.134185
##    6..7     7..1     7..3     6..2     6..8     8..4     8..5  
##  0.006681 0.059413 0.107933 0.069709 0.152878 0.349886 0.547262 0.862012 0.022762 0.000010 0.052939 0.010879 0.234885 0.344774 0.323997 0.189846
## 
## tree length =   1.29376
## 
## ((Human, Gorilla), Chimpanzee, (Orangutan, Gibbon));
## 
## ((Human: 0.059413, Gorilla: 0.107933): 0.006681, Chimpanzee: 0.069709, (Orangutan: 0.349886, Gibbon: 0.547262): 0.152878);
## 
## Detailed output identifying parameters
## rate (kappa or abcde) under REV:  0.86201  0.02276  0.00001  0.05294  0.01088
## Base frequencies:
##   0.23489  0.34477  0.32400  0.09634
## Rate matrix Q, Average Ts/Tv =  12.6163
##    -1.396743    1.362918    0.033820    0.000004
##     0.928520   -1.011984    0.078657    0.004807
##     0.024518    0.083701   -0.550045    0.441826
##     0.000011    0.017201    1.485810   -1.503021
## 
## alpha (gamma, K=4) =  0.18985
## rate:   0.00037  0.02813  0.35664  3.61486
## freq:   0.25000  0.25000  0.25000  0.25000

We 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.738

Simulating data using parametric bootstrap

We 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:

1. The tree with branch lengths
2. The rate parameters
3. The alpha shape parameter
4. The base frequencies

You file should look something like this:

# terminal
cat MCbase_Tree10.dat

##  0          * 0: paml format (mc.paml); 1:paup format (mc.nex)
## 23456721   * random number seed (odd number)
## 
## 5 895 500  * <# seqs>  <# nucleotide sites>  <# replicates>
## -1         * <tree length, use -1 if tree below has absolute branch lengths>
## ((Human: 0.058193, Gorilla: 0.100083): 0.004490, Chimpanzee: 0.068743, (Orangutan: 0.278122, Gibbon: 0.442771): 0.120082);
## 
## 7          * model: 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85, 5:T92, 6:TN93, 7:REV
## 0.87875  0.02635  0.00001  0.07245  0.01731
## 0.21924 4                                  * alpha and # of categories for discrete gamma
## 0.25318  0.32894  0.31196  0.10592
## 
## ==================================================
## The rest of this data file are notes, ignored by the program evolver.
## Change the values of the parameters, but do not delete them.
## evolver simulates nucleotide sequences under the REV+Gamma model
## and its simpler forms.
## 
## More notes:
##   Parameter kappa or rate parameters in the substituton model: 
##     For TN93, two kappa values are required, while for REV, 5 values
##     (a,b,c,d,e) are required (see Yang 1994 for the definition of these
##     parameters).  
##     The kappa parameter is defined differently under HKY85 (when k=1 means
##     no transition bias) and under F84 (when k=0 means no bias).
##     JC69 and F81 are considered species cases of HKY85, so use 1 for kappa
##     for those two models.  Notation is from my two papers in JME in 1994.
##   
##   Use equal base frequencies (0.25) for JC69 and K80.
##   Use 0 for alpha to have the same rate for all sites.
##   Use 0 for <#categories for discrete gamma> to use the continuous gamma
## 
## =========!! Check screen output carefully !! =====

We 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.dat

Analyse the simulated data

mc.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:

1. Change seqfile to mc.paml
2. Add a option ndata to specify multiple datasets
3. Change outfile to save output to another file

# terminal
cat baseml_tree10.ctl

## seqfile = mc.txt
## ndata = 500
## noisy = 1
## verbose = 0
## outfile = tree10.sim.out
## treefile = 5s.unrooted.trees
## model = 7  * GTR
## fix_alpha = 0 * 0: estimate alpha
## alpha = 0.2
## nhomo = 1

[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.ctl

The 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.

• Pipe the grep output to a file.
• Have a look at the output

# terminal
grep "^    10" tree10.sim.out > tree10_loglik.txt
head tree10_loglik.txt

##     10    -2691.897    -0.000     0.000   0.276     0.726    0.162
##     10    -2574.089    -0.110     0.492   0.412     0.800    0.401
##     10*   -2652.183     0.000     0.000  -1.000    -1.000    0.661
##     10    -2428.675    -0.616     1.373   0.327     0.751    0.208
##     10*   -2631.700     0.000     0.000  -1.000    -1.000    0.482
##     10*   -2600.009     0.000     0.000  -1.000    -1.000    0.421
##     10    -2508.045    -0.001     0.054   0.490     0.823    0.212
##     10    -2669.740    -0.000     0.000   0.121     0.796    0.065
##     10*   -2608.340     0.000     0.000  -1.000    -1.000    0.440
##     10    -2594.436    -0.296     0.861   0.365     0.719    0.203

Each 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)
p

## [1] 0.012

Question Do we reject the null hypothesis and say that tree 10 is not as good a tree as the maximum likelihood tree (tree 7)?

Model adequacy

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.

1. Calculate the test statistic for the real data
2. Generate simulated data sets using the maximum likelihood estimates.
3. Analyse each simulated data set with the same model and topology as the MLE.
4. Calculate delta for each simulated data set: the difference between the log-likelihood of the simulated data and the log-likelihood of the unconstrained model.
5. See whether the test statistic from the real data is greater than 5% (say) of the collected deltas from the simulated data.

Step 1

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.out

## TREE #  7:  ((1, 2), 3, (4, 5));  MP score: 357.00
## lnL(ntime:  7  np: 16):  -2616.073763    +9.797983
##    6..7     7..1     7..2     6..3     6..8     8..4     8..5  
##  0.035490 0.057986 0.074612 0.074353 0.131399 0.350171 0.544628 0.867230 0.026478 0.000010 0.049286 0.011098 0.235002 0.345871 0.322998 0.192488
## 
## tree length =   1.26864
## 
## ((Human, Chimpanzee), Gorilla, (Orangutan, Gibbon));
## 
## ((Human: 0.057986, Chimpanzee: 0.074612): 0.035490, Gorilla: 0.074353, (Orangutan: 0.350171, Gibbon: 0.544628): 0.131399);
## 
## Detailed output identifying parameters
## rate (kappa or abcde) under REV:  0.86723  0.02648  0.00001  0.04929  0.01110
## Base frequencies:
##   0.23500  0.34587  0.32300  0.09613
## Rate matrix Q, Average Ts/Tv =  12.8774
##    -1.409676    1.370593    0.039079    0.000004
##     0.931250   -1.008866    0.072741    0.004875
##     0.028432    0.077892   -0.545578    0.439253
##     0.000011    0.017539    1.475914   -1.493464
## 
## alpha (gamma, K=4) =  0.19249
## rate:   0.00041  0.02956  0.36377  3.60627
## freq:   0.25000  0.25000  0.25000  0.25000

Also note the unconstrained log-likelihood:

# terminal
grep "unconstrained" results.out

## ln Lmax (unconstrained) = -2476.965133

We can calculate the value of the test statistic: -2476.965133 - -2616.073763 = 139.1

d <- -2476.965133 - -2616.073763

Step 2

Use 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.dat

##  0          * 0: paml format (mc.paml); 1:paup format (mc.nex)
## 23456721   * random number seed (odd number)
## 
## 5 895 500  * <# seqs>  <# nucleotide sites>  <# replicates>
## -1         * <tree length, use -1 if tree below has absolute branch lengths>
## ((Human: 0.057986, Chimpanzee: 0.074612): 0.035490, Gorilla: 0.074353, (Orangutan: 0.350170, Gibbon: 0.544628): 0.131399);
## 
## 7          * model: 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85, 5:T92, 6:TN93, 7:REV
## 0.86723  0.02648  0.00001  0.04929  0.01110
## 0.19249 4                                  * alpha and # of categories for discrete gamma
## 0.23500  0.34587  0.32300  0.09613
## 
## ==================================================
## The rest of this data file are notes, ignored by the program evolver.
## Change the values of the parameters, but do not delete them.
## evolver simulates nucleotide sequences under the REV+Gamma model
## and its simpler forms.
## 
## More notes:
##   Parameter kappa or rate parameters in the substituton model: 
##     For TN93, two kappa values are required, while for REV, 5 values
##     (a,b,c,d,e) are required (see Yang 1994 for the definition of these
##     parameters).  
##     The kappa parameter is defined differently under HKY85 (when k=1 means
##     no transition bias) and under F84 (when k=0 means no bias).
##     JC69 and F81 are considered species cases of HKY85, so use 1 for kappa
##     for those two models.  Notation is from my two papers in JME in 1994.
##   
##   Use equal base frequencies (0.25) for JC69 and K80.
##   Use 0 for alpha to have the same rate for all sites.
##   Use 0 for <#categories for discrete gamma> to use the continuous gamma
## 
## =========!! Check screen output carefully !! =====

Run evolver as before.

# terminal
evolver 5 MCbase_Tree7.dat

Step 3

The 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.ctl

## seqfile = mc.txt
## ndata = 500
## noisy = 1
## verbose = 0
## outfile = tree7.sim.out
## treefile = 5s.unrooted.trees
## model = 7  * GTR
## fix_alpha = 0 * 0: estimate alpha
## alpha = 0.2
## nhomo = 1

Run baseml:

# terminal
baseml baseml_tree7.ctl

Step 4

We 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:

1. Extract and save the lines beginning with “lnL”
2. Extract and save the lines containing “unconstrained”

# terminal
grep "^lnL" tree7.sim.out > tree7.sim.lnl.txt
grep "unconstrained" tree7.sim.out > tree7.sim.unconstrained.txt

Read 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)

Step 5

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_signif

## [1] FALSE

Question What does this mean?

Automatic model selection

mT <- modelTest(brown_aln, brown_tree)

## Model        df  logLik   AIC      BIC
##           JC 7 -2914.115 5842.23 5875.808 
##         JC+I 8 -2901.647 5819.294 5857.669 
##      JC+G(4) 8 -2902.185 5820.369 5858.744 
##    JC+G(4)+I 9 -2902.162 5822.325 5865.496 
##          F81 10 -2838.461 5696.923 5744.891 
##        F81+I 11 -2821.84 5665.68 5718.445 
##     F81+G(4) 11 -2822.169 5666.337 5719.102 
##   F81+G(4)+I 12 -2822.156 5668.313 5725.875 
##          K80 8 -2748.411 5512.822 5551.197 
##        K80+I 9 -2727.199 5472.399 5515.57 
##     K80+G(4) 9 -2726.494 5470.988 5514.159 
##   K80+G(4)+I 10 -2726.503 5473.005 5520.973 
##          HKY 11 -2664.926 5351.851 5404.616 
##        HKY+I 12 -2619.036 5262.072 5319.634 
##     HKY+G(4) 12 -2618.111 5260.222 5317.784 
##   HKY+G(4)+I 13 -2617.936 5261.872 5324.23 
##         TrNe 9 -2737.837 5493.675 5536.846 
##       TrNe+I 10 -2716.408 5452.817 5500.785 
##    TrNe+G(4) 10 -2715.902 5451.804 5499.772 
##  TrNe+G(4)+I 11 -2715.919 5453.837 5506.602 
##          TrN 12 -2664.908 5353.816 5411.378 
##        TrN+I 13 -2616.792 5259.584 5321.943 
##     TrN+G(4) 13 -2617.896 5261.791 5324.15 
##   TrN+G(4)+I 14 -2616.233 5260.466 5327.622 
##         TPM1 9 -2746.1 5510.2 5553.371 
##       TPM1+I 10 -2724.635 5469.27 5517.239 
##    TPM1+G(4) 10 -2723.912 5467.824 5515.792 
##  TPM1+G(4)+I 11 -2723.918 5469.837 5522.602 
##          K81 9 -2746.1 5510.2 5553.371 
##        K81+I 10 -2724.635 5469.27 5517.239 
##     K81+G(4) 10 -2723.912 5467.824 5515.792 
##   K81+G(4)+I 11 -2723.918 5469.837 5522.602 
##        TPM1u 12 -2663.726 5351.452 5409.014 
##      TPM1u+I 13 -2618.521 5263.042 5325.401 
##   TPM1u+G(4) 13 -2617.702 5261.405 5323.763 
## TPM1u+G(4)+I 14 -2617.52 5263.039 5330.195 
##         TPM2 9 -2735.73 5489.459 5532.631 
##       TPM2+I 10 -2714.854 5449.708 5497.676 
##    TPM2+G(4) 10 -2714.14 5448.281 5496.249 
##  TPM2+G(4)+I 11 -2714.163 5450.325 5503.09 
##        TPM2u 12 -2661.025 5346.049 5403.611 
##      TPM2u+I 13 -2617.104 5260.208 5322.567 
##   TPM2u+G(4) 13 -2616.87 5259.739 5322.098 
## TPM2u+G(4)+I 14 -2616.435 5260.869 5328.025 
##         TPM3 9 -2744.343 5506.686 5549.858 
##       TPM3+I 10 -2724.5 5469 5516.968 
##    TPM3+G(4) 10 -2723.71 5467.42 5515.389 
##  TPM3+G(4)+I 11 -2723.718 5469.435 5522.2 
##        TPM3u 12 -2661.738 5347.476 5405.038 
##      TPM3u+I 13 -2618.824 5263.648 5326.006 
##   TPM3u+G(4) 13 -2617.64 5261.281 5323.639 
## TPM3u+G(4)+I 14 -2617.607 5263.214 5330.37 
##        TIM1e 10 -2735.55 5491.099 5539.068 
##      TIM1e+I 11 -2713.978 5449.957 5502.722 
##   TIM1e+G(4) 11 -2713.463 5448.925 5501.69 
## TIM1e+G(4)+I 12 -2713.468 5450.936 5508.498 
##         TIM1 13 -2663.706 5353.413 5415.772 
##       TIM1+I 14 -2616.21 5260.419 5327.575 
##    TIM1+G(4) 14 -2617.512 5263.023 5330.179 
##  TIM1+G(4)+I 15 -2615.736 5261.473 5333.425 
##        TIM2e 10 -2725.071 5470.142 5518.11 
##      TIM2e+I 11 -2704.522 5431.044 5483.809 
##   TIM2e+G(4) 11 -2704.285 5430.57 5483.335 
## TIM2e+G(4)+I 12 -2704.29 5432.581 5490.143 
##         TIM2 13 -2661.021 5348.041 5410.4 
##       TIM2+I 14 -2614.988 5257.977 5325.132 
##    TIM2+G(4) 14 -2616.75 5261.5 5328.656 
##  TIM2+G(4)+I 15 -2614.798 5259.597 5331.549 
##        TIM3e 10 -2733.797 5487.593 5535.561 
##      TIM3e+I 11 -2713.786 5449.572 5502.337 
##   TIM3e+G(4) 11 -2713.06 5448.12 5500.885 
## TIM3e+G(4)+I 12 -2713.081 5450.161 5507.723 
##         TIM3 13 -2661.715 5349.431 5411.79 
##       TIM3+I 14 -2616.413 5260.826 5327.981 
##    TIM3+G(4) 14 -2617.383 5262.766 5329.922 
##  TIM3+G(4)+I 15 -2615.844 5261.689 5333.641 
##         TVMe 11 -2730.753 5483.506 5536.271 
##       TVMe+I 12 -2711.22 5446.44 5504.002 
##    TVMe+G(4) 12 -2710.586 5445.171 5502.733 
##  TVMe+G(4)+I 13 -2710.61 5447.22 5509.579 
##          TVM 14 -2657.059 5342.118 5409.274 
##        TVM+I 15 -2617.188 5264.376 5336.328 
##     TVM+G(4) 15 -2616.376 5262.753 5334.705 
##   TVM+G(4)+I 16 -2616.286 5264.572 5341.321 
##          SYM 12 -2720.114 5464.229 5521.791 
##        SYM+I 13 -2700.846 5427.693 5490.051 
##     SYM+G(4) 13 -2700.519 5427.037 5489.396 
##   SYM+G(4)+I 14 -2700.534 5429.067 5496.223 
##          GTR 15 -2657.053 5344.105 5416.058 
##        GTR+I 16 -2614.639 5261.279 5338.028 
##     GTR+G(4) 16 -2616.223 5264.447 5341.196 
##   GTR+G(4)+I 17 -2614.439 5262.877 5344.423

mT <- mT[order(mT$AIC),]  # sort by AIC
mT