Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions tutorials/fbd_simple/tests.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
FBD_tutorial.Rev
40 changes: 20 additions & 20 deletions tutorials/host_rep/host_rep.md
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
Expand Down
2 changes: 2 additions & 0 deletions tutorials/intro/tests.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
linear_regression_generative.Rev
linear_regression.Rev
Empty file modified tutorials/intro_posterior_prediction/index.md
100755 → 100644
Empty file.
Empty file.
Empty file.
Empty file modified tutorials/intro_posterior_prediction/pps_TwoNormals_results.txt
100755 → 100644
Empty file.
Empty file modified tutorials/intro_posterior_prediction/twoNormals_means.log
100755 → 100644
Empty file.
Empty file modified tutorials/intro_posterior_prediction/twoNormals_posterior.log
100755 → 100644
Empty file.
Empty file modified tutorials/intro_posterior_prediction/twoNormals_sds.log
100755 → 100644
Empty file.
1 change: 1 addition & 0 deletions tutorials/mcmc/tests.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Binomial_MCMC.Rev
2 changes: 1 addition & 1 deletion tutorials/mcmc_troubleshooting/Priors.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
7 changes: 4 additions & 3 deletions tutorials/mcmc_troubleshooting/scripts/age_calibrations.Rev
Original file line number Diff line number Diff line change
@@ -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 )
Expand Down Expand Up @@ -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)
Expand Down
1 change: 1 addition & 0 deletions tutorials/mcmc_troubleshooting/tests.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
age_calibrations.Rev