Pbw2025 update

_bibliography/references.bib

@@ -3558,7 +3558,7 @@ @article{Sanmartin2008
 @article{Sanmartin2010,
   title={Bayesian island biogeography in a continental setting: the Rand Flora case},
   author={Sanmart{\'i}n, Isabel and Anderson, Cajsa Lisa and Alarcon, Marisa and Ronquist, Fredrik and Aldasoro, Juan Jos{\'e}},
-  journal={Biology letters},
+  journal={Biology Letters},
   volume={6},
   number={5},
   pages={703--707},
@@ -3585,4 +3585,15 @@ @ARTICLE{Swiston2023
   year    =  2023,
 }
 
+@article{Swiston2025,
+  title={Testing relationships between multiple regional features and biogeographic processes of speciation, extinction, and dispersal},
+  author={Swiston, Sarah K and Landis, Michael J},
+  journal={Systematic Biology},
+  volume={74},
+  pages={282-300},
+  year={2025}
+}
+t:set nonu
+
+
 @Comment{jabref-meta: databaseType:bibtex;}

home/tutorials.md

@@ -196,16 +196,16 @@ Geographic State-dependent Speciation Extinction (GeoSSE), and Feature-Informed
 
 - Pre-requisites: **RevBayes fundamentals**, **Molecular phylogenetic analysis**, and **Background on state-dependent diversification-rate estimation**.
 
-Those interested in Dispersal-Extinction-Cladogenesis models should read:
-- [Introduction to DEC models]({{ base.url }}/tutorials/biogeo/biogeo_intro)
-- [Simple DEC]({{ base.url }}/tutorials/biogeo/biogeo_simple)
-- [Epoch DEC]({{ base.url }}/tutorials/biogeo/biogeo_epoch)
-- [Biogeographic Dating with DEC]({{ base.url }}/tutorials/biogeo/biogeo_dating)
-
-Those interested in GeoSSE/FIG models should read:
-- [Simple GeoSSE]({{ base.url }}/tutorials/geosse)
-- [MultiFIG]({{ base.url }}/tutorials/multifig)
-- [TimeFIG]({{ base.url }}/tutorials/timefig)
-- [Biogeographic Dating with TimeFIG]({{ base.url }}/tutorials/timefig_dating)
+These tutorials are for those interested in DEC models:
+- Provides a basic background for discrete biogeography models, with an emphasis on DEC: [Introduction to DEC models]({{ base.url }}/tutorials/biogeo/biogeo_intro)
+- A simple time-constant DEC analysis on a fixed phylogeny within a 4-region system: [Simple DEC]({{ base.url }}/tutorials/biogeo/biogeo_simple)
+- A paleogeography-informed DEC analysis on a fixed phylogeny within a 6-region system: [Epoch DEC]({{ base.url }}/tutorials/biogeo/biogeo_epoch)
+- A paleogeography-informed DEC analysis on a jointly inferred phylogeny (topology, divergence times) within a 6-region system: [Biogeographic Dating with DEC]({{ base.url }}/tutorials/biogeo/biogeo_dating)
+
+These tutorials are for those interested in GeoSSE/FIG models:
+- A simple GeoSSE analysis using a 2-region system: [Simple GeoSSE]({{ base.url }}/tutorials/geosse)
+- A time-constant FIG analysis in a 7-region system with multiple regional features:  [MultiFIG]({{ base.url }}/tutorials/multifig)
+- A paleogeography-informed FIG analysis on a fixed phylogeny in a 7-region system with multiple regional features:  [TimeFIG]({{ base.url }}/tutorials/timefig_simple)
+- A paleogeography-informed FIG analysis on a jointly inferred phylogeny (topology, divergence times) within a 7-region system with multiple regional features [Biogeographic Dating with TimeFIG]({{ base.url }}/tutorials/timefig_dating)
 
 {% endaside %}

tutorials/geosse/index.md

@@ -317,9 +317,9 @@ We also want MCMC to keep track of certain things while it runs. We want it to p
 ```
 mni = 1
 mn[mni++] = mnScreen(printgen=printgen)
-mn[mni++] = mnModel(printgen=printgen, filename=out_fp+"model.txt")
-mn[mni++] = mnJointConditionalAncestralState(glhbdsp=timetree, tree=timetree, printgen=printgen, filename=out_fp+"states.txt", withTips=true, withStartStates=true, type="NaturalNumbers")
-mn[mni++] = mnStochasticCharacterMap(glhbdsp=timetree, printgen=printgen, filename=out_fp+"stoch.txt")
+mn[mni++] = mnModel(printgen=printgen, filename=out_fp+"model.log")
+mn[mni++] = mnJointConditionalAncestralState(glhbdsp=timetree, tree=timetree, printgen=printgen, filename=out_fp+"states.log", withTips=true, withStartStates=true, type="NaturalNumbers")
+mn[mni++] = mnStochasticCharacterMap(glhbdsp=timetree, printgen=printgen, filename=out_fp+"stoch.log")
 ```
 
 Let's also store information for how the integer-valued ranges (0, 1, 2) relate to the regional presence-absence representation of ranges (A=10, B=01, AB=11).
@@ -338,16 +338,16 @@ Then we can start up the MCMC. It doesn't matter which model parameter you use t
 ```
 mdl = model(m_w)
 ch = mcmc(mv, mn, mdl)
-ch.burnin(n_burn, tuningInterval=10)
+ch.burnin(n_burn, tuningInterval=50)
 ch.run(n_gen)
 ```
 
 After the MCMC analysis has concluded, we can summarize the ancestral states we obtained, creating an ancestral state tree. This tree will be written to the file `ase.tre`. It may take a little while.
 
 ```
 f_burn = 0.2
-x_stoch = readAncestralStateTrace(file=out_fp+"stoch.txt")
-x_states = readAncestralStateTrace(file=out_fp+"states.txt")
+x_stoch = readAncestralStateTrace(file=out_fp+"stoch.log")
+x_states = readAncestralStateTrace(file=out_fp+"states.log")
 summarizeCharacterMaps(x_stoch,timetree,file=out_fp+"events.tsv",burnin=f_burn)
 state_tree = ancestralStateTree(tree=timetree,
                    ancestral_state_trace_vector=x_states,

tutorials/geosse/scripts/geosse.Rev

@@ -1,5 +1,6 @@
 # load TensorPhylo plugin
-loadPlugin("TensorPhylo")
+#loadPlugin("TensorPhylo")
+loadPlugin("TensorPhylo","/Users/mlandis/.local/lib/tensorphylo")
 
 # FILESYSTEM
 fp          = "./"
@@ -180,11 +181,11 @@ mv[mvi++] = mvSimplex(m_d_simplex, weight=m_d_simplex.size())
 # MCMC MONITORS
 mni = 1
 mn[mni++] = mnScreen(printgen=printgen)
-mn[mni++] = mnModel(printgen=printgen, filename=out_fp + "model.txt")
+mn[mni++] = mnModel(printgen=printgen, filename=out_fp + "model.log")
 mn[mni++] = mnJointConditionalAncestralState(glhbdsp=timetree, tree=timetree,
                                              printgen=printgen, withTips=true, withStartStates=true,
-                                             type="NaturalNumbers", filename=out_fp + "states.txt")
-mn[mni++] = mnStochasticCharacterMap(glhbdsp=timetree, printgen=printgen, filename=out_fp + "stoch.txt")
+                                             type="NaturalNumbers", filename=out_fp + "states.log")
+mn[mni++] = mnStochasticCharacterMap(glhbdsp=timetree, printgen=printgen, filename=out_fp + "stoch.log")
 
 # STATE LABELS
 write("index,range\n", file=lbl_fn)
@@ -196,13 +197,13 @@ for (i in 1:state_labels.size()) {
 # MCMC START
 mdl = model(m_w)
 ch = mcmc(mv, mn, mdl)
-ch.burnin(n_burn, tuningInterval=10)
+ch.burnin(n_burn, tuningInterval=50)
 ch.run(n_gen)
 
 # ANCESTRAL STATES
 f_burn = 0.2
-x_stoch = readAncestralStateTrace(file=out_fp + "stoch.txt")
-x_states = readAncestralStateTrace(file=out_fp + "states.txt")
+x_stoch = readAncestralStateTrace(file=out_fp + "stoch.log")
+x_states = readAncestralStateTrace(file=out_fp + "states.log")
 summarizeCharacterMaps(x_stoch,timetree,file=out_fp + "events.tsv",burnin=f_burn)
 state_tree = ancestralStateTree(tree=timetree,
                    ancestral_state_trace_vector=x_states,

tutorials/multifig/data/hawaii/qw_feature1.csv

@@ -1,2 +1,2 @@
 G,N,K,O,M,H,Z
-005,100,1431,1548,3155,10434,1500
+005,100,1431,1548,3155,4207,1500

tutorials/multifig/data/hawaii/qw_feature2.csv

@@ -1,2 +1,2 @@
 G,N,K,O,M,H,Z
-1.6094379124341,4.60517018598809,7.26612877955645,7.34471905414967,8.05674377497531,9.25282498358234,7.3132203870903
+1.6094379124341,4.60517018598809,7.26612877955645,7.34471905414967,8.05674377497531,8.34450508359,7.3132203870903

tutorials/multifig/index.md

@@ -15,7 +15,7 @@ prerequisites:
 
 In the previous examples, we used a GeoSSE model {% cite Goldberg2011 %} to investigate the evolution of the Hawaiian *Kadua*. The GeoSSE model allows us to estimate rates of within-region speciation, extinction, between-region speciation, and dispersal that differ among regions. Biologically, we expect that these different rates are informed by features of the regions where the species are evolving. For example, we might expect that species disperse at a lower rate between more distant islands, or go extinct at a higher rate on smaller islands.
 
-The FIG model {% cite Landis2022 %} and the Multiple Feature-Informed GeoSSE (MultiFIG) model {% cite Swiston2023 %} attempt to model this expectation. Rather than giving each region its own evolutionary rate parameters, FIG models use functions to link features of those regions to evolutionary rates. This allows us to test hypotheses about the importance of certain environmental features on evolutionary processes. It also has the benefit of reducing the number of parameters that need to be estimated. The number of parameters in the MultiFIG model is constant with respect to the number of regions, so we can investigate systems with more regions without suffering an explosion in the number of model parameters. In this tutorial, we will model the evolution and biogeography of *Kadua* using seven regions and eight regional features. Later tutorials will explore how to adapt FIG to allow for regional features, and their linked biogeographic rates, to change over time.
+The FIG model {% cite Landis2022 %} and the Multiple Feature-Informed GeoSSE (MultiFIG) model {% cite Swiston2025 %} attempt to model this expectation. Rather than giving each region its own evolutionary rate parameters, FIG models use functions to link features of those regions to evolutionary rates. This allows us to test hypotheses about the importance of certain environmental features on evolutionary processes. It also has the benefit of reducing the number of parameters that need to be estimated. The number of parameters in the MultiFIG model is constant with respect to the number of regions, so we can investigate systems with more regions without suffering an explosion in the number of model parameters. In this tutorial, we will model the evolution and biogeography of *Kadua* using seven regions and eight regional features. Later tutorials will explore how to adapt FIG to allow for regional features, and their linked biogeographic rates, to change over time.
 
 {% subsection The MultiFIG model %}
 
@@ -68,7 +68,7 @@ $$
 
 and behaves similarly to the quantitative effect variable. (In this example, we assume categorical variables take values 0 or 1 and treat them numerically to simplify notation. More complex relationships are written with more complex notation.)
 
-Each individual effect variable can be $<1$, $>1$, or $=1$, and so can the product of all effect variables. All relative rates are eventually rescaled by the same base rate (e.g. $\rho_w$). If the relative rates for regions $i$ and $j$ have the relationship $m_w(i) > m_w(j)$ then the absolute rates also follow $r_w(i) > r_w(j)$. Other $m$ functions behave in a similar manner. More details on the design of the $m$ functions are provided in {% cite Swiston2023 %}.
+Each individual effect variable can be $<1$, $>1$, or $=1$, and so can the product of all effect variables. All relative rates are eventually rescaled by the same base rate (e.g. $\rho_w$). If the relative rates for regions $i$ and $j$ have the relationship $m_w(i) > m_w(j)$ then the absolute rates also follow $r_w(i) > r_w(j)$. Other $m$ functions behave in a similar manner. More details on the design of the $m$ functions are provided in {% cite Swiston2025 %}.
 
 In this analysis, we are examining eight regional features. The first 4 are quantitative: maximum altitude (m), log maximum altitude (m), distance (km), and log distance (km). We include the log features because they will allow us to better understand the *shape* of the relationship between features and processes. For example, it may be that intermediate values of a particular feature are related to the highest rates of a particular process, so we would expect the feature strength parameter to be positive and the log-feature strength parameter to be negative. The other 4 features are categorical: age class (old/young), growth class (decay/growth), dispersal class (short/long), and relative age class (older/younger).
 
@@ -139,8 +139,9 @@ Eventually, we will be running an MCMC analysis at the end of the tutorial. This
 ```
 # MCMC variables
 num_proc  = 6
-num_gen   = 500             # set num_gen = 5000 for full analysis
-print_gen = 20
+num_gen   = 50             # set num_gen = 5000 for full analysis
+print_gen = 1
+save_gen  = 10
 moves     = VectorMoves()
 monitors  = VectorMonitors()
 ```
@@ -285,30 +286,30 @@ First, we will address the categorical feature effects for each process (w, e, d
 # initialize categorical feature effects, create moves, add monitor variables
 for (i in 1:feature_CW.size()) {
     sigma_w[i].setValue(0)
-    moves.append( mvScale(sigma_w[i], weight=2) )
-    moves.append( mvSlide(sigma_w[i], weight=2) )
-    moves.append( mvRJSwitch(sigma_w[i], weight=3) )
+    moves.append( mvScale(sigma_w[i], weight=1) )
+    moves.append( mvSlide(sigma_w[i], weight=1) )
+    moves.append( mvRJSwitch(sigma_w[i], weight=2) )
     use_sigma_w[i] := ifelse(sigma_w[i] == 0.0, 0, 1)
 }
 for (i in 1:feature_CW.size()) {
     sigma_e[i].setValue(0)
-    moves.append( mvScale(sigma_e[i], weight=2) )
-    moves.append( mvSlide(sigma_e[i], weight=2) )
-    moves.append( mvRJSwitch(sigma_e[i], weight=3) )
+    moves.append( mvScale(sigma_e[i], weight=1) )
+    moves.append( mvSlide(sigma_e[i], weight=1) )
+    moves.append( mvRJSwitch(sigma_e[i], weight=2) )
     use_sigma_e[i] := ifelse(sigma_e[i] == 0.0, 0, 1)
 }
 for (i in 1:feature_CB.size()) {
     sigma_d[i].setValue(0)
-    moves.append( mvScale(sigma_d[i], weight=2) )
-    moves.append( mvSlide(sigma_d[i], weight=2) )
-    moves.append( mvRJSwitch(sigma_d[i], weight=3) )
+    moves.append( mvScale(sigma_d[i], weight=1) )
+    moves.append( mvSlide(sigma_d[i], weight=1) )
+    moves.append( mvRJSwitch(sigma_d[i], weight=2) )
     use_sigma_d[i] := ifelse(sigma_d[i] == 0.0, 0, 1)
 }
 for (i in 1:feature_CB.size()) {
     sigma_b[i].setValue(0)
-    moves.append( mvScale(sigma_b[i], weight=2) )
-    moves.append( mvSlide(sigma_b[i], weight=2) )
-    moves.append( mvRJSwitch(sigma_b[i], weight=3) )
+    moves.append( mvScale(sigma_b[i], weight=1) )
+    moves.append( mvSlide(sigma_b[i], weight=1) )
+    moves.append( mvRJSwitch(sigma_b[i], weight=2) )
     use_sigma_b[i] := ifelse(sigma_b[i] == 0.0, 0, 1)
 }
 ```
@@ -319,30 +320,30 @@ Similarly, we will address the quantitative features for each process. These are
 # initialize quantitative feature effects, create moves, add monitor variables
 for (i in 1:feature_QW.size()) {
     phi_w[i].setValue(0)
-    moves.append( mvScale(phi_w[i], weight=2) )
-    moves.append( mvSlide(phi_w[i], weight=2) )
-    moves.append( mvRJSwitch(phi_w[i], weight=3) )
+    moves.append( mvScale(phi_w[i], weight=1) )
+    moves.append( mvSlide(phi_w[i], weight=1) )
+    moves.append( mvRJSwitch(phi_w[i], weight=2) )
     use_phi_w[i] := ifelse(phi_w[i] == 0.0, 0, 1)
 }
 for (i in 1:feature_QW.size()) {
     phi_e[i].setValue(0)
-    moves.append( mvScale(phi_e[i], weight=2) )
-    moves.append( mvSlide(phi_e[i], weight=2) )
-    moves.append( mvRJSwitch(phi_e[i], weight=3) )
+    moves.append( mvScale(phi_e[i], weight=1) )
+    moves.append( mvSlide(phi_e[i], weight=1) )
+    moves.append( mvRJSwitch(phi_e[i], weight=2) )
     use_phi_e[i] := ifelse(phi_e[i] == 0.0, 0, 1)
 }
 for (i in 1:feature_QB.size()) {
     phi_d[i].setValue(0)
-    moves.append( mvScale(phi_d[i], weight=2) )
-    moves.append( mvSlide(phi_d[i], weight=2) )
-    moves.append( mvRJSwitch(phi_d[i], weight=3) )
+    moves.append( mvScale(phi_d[i], weight=1) )
+    moves.append( mvSlide(phi_d[i], weight=1) )
+    moves.append( mvRJSwitch(phi_d[i], weight=2) )
     use_phi_d[i] := ifelse(phi_d[i] == 0.0, 0, 1)
 }
 for (i in 1:feature_QB.size()) {
     phi_b[i].setValue(0)
-    moves.append( mvScale(phi_b[i], weight=2) )
-    moves.append( mvSlide(phi_b[i], weight=2) )
-    moves.append( mvRJSwitch(phi_b[i], weight=3) )
+    moves.append( mvScale(phi_b[i], weight=1) )
+    moves.append( mvSlide(phi_b[i], weight=1) )
+    moves.append( mvRJSwitch(phi_b[i], weight=2) )
     use_phi_b[i] := ifelse(phi_b[i] == 0.0, 0, 1)
 }
 ```
@@ -513,20 +514,20 @@ Monitors are instructions for RevBayes to record MCMC output. Since we want RevB
 monitors.append( mnScreen(rho_d, rho_e, rho_w, rho_b, printgen=print_gen) )
 
 # file monitor for all simple model variables
-monitors.append( mnModel(printgen=print_gen, file=out_fn+".model.txt") )
+monitors.append( mnModel(printgen=save_gen, file=out_fn+".model.txt") )
 
 # file monitor for tree
-monitors.append( mnFile(timetree, printgen=print_gen, file=out_fn + ".tre") )
+monitors.append( mnFile(timetree, printgen=save_gen, file=out_fn + ".tre") )
 
 # monitor ancestral ranges at internal nodes
 monitors.append( mnJointConditionalAncestralState(
-    tree=timetree, glhbdsp=timetree, printgen=print_gen,
+    tree=timetree, glhbdsp=timetree, printgen=save_gen,
     filename=out_fn+".states.txt",
     withTips=true, withStartStates=true, type="NaturalNumbers") )
 
 # file monitor for biogeographic rates
 bg_mon_fn = out_fn + ".bg.txt"
-monitors.append( mnFile( filename = bg_mon_fn, printgen=print_gen,
+monitors.append( mnFile( filename = bg_mon_fn, printgen=save_gen,
                          rho_e, rho_w, rho_d, rho_b,
                          r_e, r_w,
                          r_d[1], r_d[2], r_d[3], r_d[4],