From b559238c5f3911b6a38ff0e1fac432ed77c28ecf Mon Sep 17 00:00:00 2001 From: Benjamin Redelings Date: Thu, 14 Aug 2025 14:57:02 -0400 Subject: [PATCH 1/2] More tests. --- tutorials/fbd_simple/tests.txt | 1 + tutorials/host_rep/host_rep.md | 40 +++++++++---------- tutorials/intro/tests.txt | 2 + tutorials/intro_posterior_prediction/index.md | 0 .../pps_SingleNormal_results.txt | 0 .../pps_TwoNormals_fullMCMC_results.txt | 0 .../pps_TwoNormals_results.txt | 0 .../twoNormals_means.log | 0 .../twoNormals_posterior.log | 0 .../twoNormals_sds.log | 0 tutorials/mcmc/tests.txt | 1 + tutorials/mcmc_troubleshooting/Priors.md | 2 +- .../scripts/age_calibrations.Rev | 7 ++-- tutorials/mcmc_troubleshooting/tests.txt | 1 + 14 files changed, 30 insertions(+), 24 deletions(-) create mode 100644 tutorials/fbd_simple/tests.txt create mode 100644 tutorials/intro/tests.txt mode change 100755 => 100644 tutorials/intro_posterior_prediction/index.md mode change 100755 => 100644 tutorials/intro_posterior_prediction/pps_SingleNormal_results.txt mode change 100755 => 100644 tutorials/intro_posterior_prediction/pps_TwoNormals_fullMCMC_results.txt mode change 100755 => 100644 tutorials/intro_posterior_prediction/pps_TwoNormals_results.txt mode change 100755 => 100644 tutorials/intro_posterior_prediction/twoNormals_means.log mode change 100755 => 100644 tutorials/intro_posterior_prediction/twoNormals_posterior.log mode change 100755 => 100644 tutorials/intro_posterior_prediction/twoNormals_sds.log create mode 100644 tutorials/mcmc/tests.txt create mode 100644 tutorials/mcmc_troubleshooting/tests.txt diff --git a/tutorials/fbd_simple/tests.txt b/tutorials/fbd_simple/tests.txt new file mode 100644 index 000000000..d0dadcdcd --- /dev/null +++ b/tutorials/fbd_simple/tests.txt @@ -0,0 +1 @@ +FBD_tutorial.Rev diff --git a/tutorials/host_rep/host_rep.md b/tutorials/host_rep/host_rep.md index e517e71ab..09f9e264f 100644 --- a/tutorials/host_rep/host_rep.md +++ b/tutorials/host_rep/host_rep.md @@ -119,60 +119,60 @@ In this tutorial we are going to perform the analysis with two different host tr Now, let's begin. First, create file management variables for input - +{% snippet scripts/run_nymphalini.Rev %} phy_host_fn = "data/angio_25tips_time.phy" phy_parasite_fn = "data/Nymphalini.phy" dat_parasite_fn = "data/interaction_matrix.nex" - +{% endsnippet %} then read in our character data - +{% snippet scripts/run_nymphalini.Rev %} dat_parasite <- readDiscreteCharacterData(dat_parasite_fn) - +{% endsnippet %} For this tutorial we'll assume we know the host and parasite phylogenies without error. Note that our host repertoire inference method uses a root branch length to estimate the stationary probabilities at the root node. Our parasite tree file (`Nymphalini.phy`) is modified to have a branch length assigned to the root node. If you provide a tree without a root branch length, the software will consider it to be the same length as the tree height. - +{% snippet scripts/run_nymphalini.Rev %} phy_parasite <- readTrees(phy_parasite_fn)[1] - +{% endsnippet %} Here is where you can change the host tree to `angio_25tips_bl1.phy` when you repeat the analysis - +{% snippet scripts/run_nymphalini.Rev %} phy_host <- readTrees(phy_host_fn)[1] - +{% endsnippet %} Retrieve dataset dimensions - +{% snippet scripts/run_nymphalini.Rev %} n_host_tips <- phy_host.ntips() n_host_branches <- 2 * n_host_tips - 2 n_parasite_branches <- 2 * phy_parasite.ntips() - 2 n_sites <- dat_parasite.nchar() - +{% endsnippet %} Add more information to the name of output files - +{% snippet scripts/run_nymphalini.Rev %} out_str = "out.time" out_fn = "output/" + out_str - +{% endsnippet %} We need to create vectors of *moves* and *monitors* to later inform how our Markov chain Monte Carlo (MCMC) analysis needs to propose and sample new model parameters and host repertoire histories. Also, we use monitors to record the information we want to use - +{% snippet scripts/run_nymphalini.Rev %} moves = VectorMoves() monitors = VectorMonitors() - +{% endsnippet %} Next, we'll build the transition rate matrix for the model. In this example, the rate matrix requires four rates: two gain rates (0->1 and 1->2) and two loss rates (1->0 and 2->1). First, create a vector containing all transition rates and assign it a move - +{% snippet scripts/run_nymphalini.Rev %} switch_rates_pos ~ dnDirichlet( [1,1,1,1] ) moves.append( mvSimplex(switch_rates_pos, alpha=10, weight=2, tune=false) ) - +{% endsnippet %} We'll now create a set of deterministic nodes to help us map our simplex of transition rates onto specific host gain and loss events - +{% snippet scripts/run_nymphalini.Rev %} switch_rate_0_to_1 := switch_rates_pos[1] switch_rate_0_to_2 := 0. switch_rate_1_to_0 := switch_rates_pos[2] switch_rate_1_to_2 := switch_rates_pos[3] switch_rate_2_to_0 := 0. switch_rate_2_to_1 := switch_rates_pos[4] - +{% endsnippet %} Next, we assemble our named rate variables into a vector - +{% snippet scripts/run_nymphalini.Rev %} switch_rates := v( switch_rate_0_to_1, switch_rate_0_to_2, switch_rate_1_to_0, switch_rate_1_to_2, switch_rate_2_to_0, switch_rate_2_to_1 ) - +{% endsnippet %} We then construct a rate matrix for three states (0: non-host, 1: potential host, 2: actual host) using our vector of named rates. We found that the MCMC mixes better when the Q matrix is not rescaled such that the expected number of events per unit time per character is 1 (`rescaled=false`). This might not be true for every data set and you can always change it to `rescaled=true`. Q_char := fnFreeK( transition_rates=switch_rates, rescaled=false ) diff --git a/tutorials/intro/tests.txt b/tutorials/intro/tests.txt new file mode 100644 index 000000000..dcae8298e --- /dev/null +++ b/tutorials/intro/tests.txt @@ -0,0 +1,2 @@ +linear_regression_generative.Rev +linear_regression.Rev diff --git a/tutorials/intro_posterior_prediction/index.md b/tutorials/intro_posterior_prediction/index.md old mode 100755 new mode 100644 diff --git a/tutorials/intro_posterior_prediction/pps_SingleNormal_results.txt b/tutorials/intro_posterior_prediction/pps_SingleNormal_results.txt old mode 100755 new mode 100644 diff --git a/tutorials/intro_posterior_prediction/pps_TwoNormals_fullMCMC_results.txt b/tutorials/intro_posterior_prediction/pps_TwoNormals_fullMCMC_results.txt old mode 100755 new mode 100644 diff --git a/tutorials/intro_posterior_prediction/pps_TwoNormals_results.txt b/tutorials/intro_posterior_prediction/pps_TwoNormals_results.txt old mode 100755 new mode 100644 diff --git a/tutorials/intro_posterior_prediction/twoNormals_means.log b/tutorials/intro_posterior_prediction/twoNormals_means.log old mode 100755 new mode 100644 diff --git a/tutorials/intro_posterior_prediction/twoNormals_posterior.log b/tutorials/intro_posterior_prediction/twoNormals_posterior.log old mode 100755 new mode 100644 diff --git a/tutorials/intro_posterior_prediction/twoNormals_sds.log b/tutorials/intro_posterior_prediction/twoNormals_sds.log old mode 100755 new mode 100644 diff --git a/tutorials/mcmc/tests.txt b/tutorials/mcmc/tests.txt new file mode 100644 index 000000000..56a5d9528 --- /dev/null +++ b/tutorials/mcmc/tests.txt @@ -0,0 +1 @@ +Binomial_MCMC.Rev diff --git a/tutorials/mcmc_troubleshooting/Priors.md b/tutorials/mcmc_troubleshooting/Priors.md index 56ca95904..1f30dd99a 100644 --- a/tutorials/mcmc_troubleshooting/Priors.md +++ b/tutorials/mcmc_troubleshooting/Priors.md @@ -99,7 +99,7 @@ This interaction needs to be taken into account in order to correctly interpret Each of these clades has an associated prior distribution on the age of their most recent common ancestor (MRCA): {{ calib_script | snippet:"block#","7" }} -We also have a prior set on the origin time of the tree: +We also have a prior set on the origin time of the trees: {{ calib_script | snippet:"block#","4" }} When we run the analysis, we obtain the posterior distribution shown in {% ref fig_age_posterior %} for the _Ursavus_ clade that we have set. The first impression we might get from this result is that our inference strongly supports younger ages for this clade, because our original prior was a uniform distribution on the full range [25.0 ; 36.0] My, whereas the posterior has much higher densities for values in the first half of the range. However, this is not taking into account the interactions between priors, which may lead to an **effective** prior on the clade age which is not uniform at all. diff --git a/tutorials/mcmc_troubleshooting/scripts/age_calibrations.Rev b/tutorials/mcmc_troubleshooting/scripts/age_calibrations.Rev index 32ed5c65f..ac64ad906 100644 --- a/tutorials/mcmc_troubleshooting/scripts/age_calibrations.Rev +++ b/tutorials/mcmc_troubleshooting/scripts/age_calibrations.Rev @@ -1,11 +1,12 @@ #### Interacting priors - age calibrations # Read the full list of taxa (including all fossils and extant species # -taxa <- readTaxonData("/users/biodiv/barido/Downloads/FBD/data/bears_taxa.tsv") +seed(1) +taxa <- readTaxonData("data/bears_taxa.tsv") # Import the morphological character matrix # # this file contains only the taxa for which morphological characters are available # -morpho <- readDiscreteCharacterData("/users/biodiv/barido/Downloads/FBD/data/bears_morphology.nex") +morpho <- readDiscreteCharacterData("data/bears_morphology.nex") # Add the missing taxa to each data partition # morpho.addMissingTaxa( taxa ) @@ -74,7 +75,7 @@ phyMorpho.clamp(morpho) ### MCMC mymodel = model(fbd_tree) monitors = VectorMonitors() -monitors.append( mnModel(filename="/users/biodiv/barido/Downloads/FBD/output/bears_orig.log", printgen=1000) ) +monitors.append( mnModel(filename="output/bears_orig.log", printgen=1000) ) mymcmc = mcmc(mymodel, monitors, moves) mymcmc.run(generations=1000000) diff --git a/tutorials/mcmc_troubleshooting/tests.txt b/tutorials/mcmc_troubleshooting/tests.txt new file mode 100644 index 000000000..a0a23d18e --- /dev/null +++ b/tutorials/mcmc_troubleshooting/tests.txt @@ -0,0 +1 @@ +age_calibrations.txt From ca6a6e5c913083cfec670513378e9813bbe82c71 Mon Sep 17 00:00:00 2001 From: Benjamin Redelings Date: Thu, 14 Aug 2025 15:01:23 -0400 Subject: [PATCH 2/2] Fix typo. --- tutorials/mcmc_troubleshooting/tests.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/mcmc_troubleshooting/tests.txt b/tutorials/mcmc_troubleshooting/tests.txt index a0a23d18e..d29e82ecb 100644 --- a/tutorials/mcmc_troubleshooting/tests.txt +++ b/tutorials/mcmc_troubleshooting/tests.txt @@ -1 +1 @@ -age_calibrations.txt +age_calibrations.Rev