Identical likelihoods across alternative dispersal hypotheses in BioGeoBEARS

[Emanuel Seculi]

Dear BioGeoBEARS community,

I am currently working on the historical biogeography of Ceratosauria (a clade of non-avian theropods, with a particular focus on Abelisauridae) using BioGeoBEARS, and I am encountering an issue that I have not been able to resolve despite testing multiple approaches.

When evaluating alternative dispersal hypotheses (e.g., origins in Europe, South America, or Africa), I consistently obtain identical likelihoods, parameter estimates (d, e, j), and AICc values across all models (DEC, DEC+J, DIVALIKE, DIVALIKE+J, BAYAREALIKE, BAYAREALIKE+J) (see the attached image). This occurs even when I substantially modify dispersal multiplier matrices and area constraints.

Given this, I would like to ask:

-Could this behavior be driven by the j parameter dominating the inference (particularly in DEC+J and DIVALIKE+J)?
-Are likelihoods generally insensitive to dispersal multiplier matrices in time-homogeneous (non-stratified) analyses?
-Is it expected that non-stratified models may fail to distinguish strongly directional dispersal scenarios?
-Could overly restrictive settings (e.g., zeros in dispersal multipliers or in areas_allowed) effectively collapse the likelihood space?
-Would implementing time-stratified analyses (i.e., multiple dispersal matrices across time slices) be necessary to meaningfully test these alternative hypotheses?

I have attached the R scripts and input data used in these analyses.

Any guidance, suggestions, or relevant references would be greatly appreciated.

[Ivan Magalhães]

Hi Emanuel,

Quickly inspecting your code, it seems that your matrices have been incorrectly fed to the bgb_run objects. It seems that BioGeoBEARS is not actually using your "areas allowed" matrices, and I think the same goes for the "dispersal matrices". 

This could be because your geography file uses areas A–F: 33 6 (A B C D E F)

And your matrices use A, B, C, AS, AF, MA.

It is better to use single-letter areas, because if you plan to do biogeographic stochastic mapping, it only handles single-letter areas.

We can know BioGeoBEARS is not using your matrices because your areas allowed matrix disallows area B:

$areas_allowed
   A B C AS AF MA
A  0 0 1  1  1  1
B  0 0 0  0  0  0
C  1 0 1  1  1  1
AS 1 0 1  1  1  1
AF 1 0 1  1  1  1
MA 1 0 1  1  1  0

Because you have a single stratum, this means that no species can ever occupy area B. You do have species in this area in your geography file. So if BioGeoBEARS were using this matrix, the run would be aborted immediately, because there would be no possible estimate compatible with species occupying B, and a matrix disallowing area B. This indicates the matrix is not being used.

I think the syntax of how you are feeding the matrix to the object might be incorrect. You have something like:

$lists_of_matrices

$areas_allowed
   A B C AS AF MA
A  0 0 1  1  1  1
B  0 0 0  0  0  0
C  1 0 1  1  1  1
AS 1 0 1  1  1  1
AF 1 0 1  1  1  1
MA 1 0 1  1  1  0

A correct object that I set in my own analysis looks like this:

$list_of_areas_allowed_mats
$list_of_areas_allowed_mats[[1]]  

  M F U N S
1 1 1 1 0 0
2 1 1 1 0 0
3 1 1 1 0 0
4 0 0 0 0 0
5 0 0 0 0 0

(notice this disallows areas N and S in this stratum; they are allowed in other strata) 

My analysis was set using BioGeoBEARS native functions:

BioGeoBEARS_run_object$areas_allowed_fn = areas_allowed_fn

BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object)

I suppose you can feed the matrix in a similar way you are trying to do, but you have to get the syntax right. Same for the dispersal multiplier matrices. To be safe, perhaps use the "native" syntax above.

Additionally (and perhaps more importantly), you are setting DIVA-like and BAYAREA-like incorrectly. 

  bgb_run_DIV <- bgb_run
  bgb_run_DIV$BioGeoBEARS_model_object@params_table["e","type"] <- "fixed"
  bgb_run_DIV$BioGeoBEARS_model_object@params_table["j","type"] <- "fixed"

This results only in a modified version of DEC were e is fixed (in your case, to the starting value, which I think is 0.1) and j is fixed to 0. It essentially is a DEC with fixed extinction (probably not what you want).

For correctly setting up a DIVA-like model, you need to modify the cladogenetic events by allowing wide vicariance, disallowing subset sympatry, and adjusting the linkage between the parameters.

This comes from Nick's template scripts, so you have an idea:

# Remove subset-sympatry
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","init"] = 0.0
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["s","est"] = 0.0
# Adjust linkage between parameters
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ysv","type"] = "2-j"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ys","type"] = "ysv*1/2"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["y","type"] = "ysv*1/2"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","type"] = "ysv*1/2"
# Allow classic, widespread vicariance; all events equiprobable
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","type"] = "fixed"
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","init"] = 0.5
BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01v","est"] = 0.5

You should check Nick's example to see how DIVA-like and BAYAREA-like models are constructed: http://phylo.wikidot.com/biogeobears

Best,

Ivan

--
You received this message because you are subscribed to the Google Groups "BioGeoBEARS" group.
To unsubscribe from this group and stop receiving emails from it, send an email to biogeobears...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/biogeobears/9f2e130c-f126-4e3c-b936-ca8062c0be48n%40googlegroups.com.

[Emanuel Seculi]

Dear Ivan,

Thank you very much for your help, it was extremely useful.

You were absolutely right: there were two main issues in my analysis. First, a mismatch between the geographic file and the areas allowed matrix, and second, my implementation did not follow the standard BioGeoBEARS workflow. 

I also really appreciate your clarification regarding the parametrization of DIVA-LIKE and BAYAREA-LIKE models. I now understand that my previous settings were effectively reducing them to DEC-, which was a conceptual mistake on my side.

After correcting these issues, I am now obtaining different results across my alternative origin hypotheses. 

I would like to ask one additional question, more from a biological perspective. This clade of non-avian theropods spans a time interval characterized by major paleogeographic changes, including the breakup of Pangea into Gondwana and Laurasia, and the subsequent fragmentation of Gondwana into the present-day continents. In this context, would you recommend implementing time-stratified analyses and dispersal matrices to explicitly model these processes? Or do you think that using  areas allowed, as in my current approach, can provide a sufficiently realistic approximation of the historical geographic scenario?

Thanks again Ivan for your time and help , I really appreciate it.

Best regards, 

Emanuel

To view this discussion visit https://groups.google.com/d/msgid/biogeobears/CAAsMhYrRJxEZHH-sGYm1DFK4gkH12TmwnouOeNv3F%2BPP58z40Q%40mail.gmail.com.

--

Licenciado en Ciencias Biológicas 

Instituto Patagónico de Geología y Paleontología (IPGP)-CENPAT-CONICET

[Ivan Magalhães]

Hi Emanuel

Glad it helped.

There is no way to know a priori if a time-stratified analysis with different dispersal matrices will improve your estimates, but if you think it makes sense from a biological perspective (i.e., if you think that dispersal may be hindered once the continents are apart), I think it is worth exploring.

Also, I think it does not make sense to use an "areas allowed" matrix outside of a time-stratified analysis. As I said before, if you have a single stratum, and you disallow any area in this stratum, the run will crash because one of the areas will not be occupiable and this will conflict with the geography file. And it may not be the best approach here, because the continents have not appeared/disappeared, they only connected/disconnected (which is better modeled using dispersal matrices). I would use "areas allowed" to model areas that are only available for occupation in certain periods of time (e.g., an island that has appeared in a particular time stratum). I would recommend checking in more detail how the "areas allowed" matrices work so you can decide if they are useful for testing your hypotheses.

Best,

Ivan

To view this discussion visit https://groups.google.com/d/msgid/biogeobears/CABNT8%3DFDPn6xzOAF%3DX65-%3DuOmW%3DynmpVOaAqc-TOpnK1tUYf1A%40mail.gmail.com.

[Emanuel Seculi]

Hi Ivan,

Thanks for your answer and suggestions.

After exploring different implementations in BioGeoBEARS (some papers and Nick examples), I think the best approach is indeed to use a time-stratified framework. This allows to separate two different components of the biogeographic process in my analysis: 1) the adjacency matrices + time slices add the temporal fragmentation and connectivity of continental areas through time, reflecting the tectonic history of Gondwana and Laurasia and 2) with the dispersal multiplier matrices I can evaluate different probabilities of dispersal among areas under distinct paleobiogeographic scenarios (e.g., South American, African, or European origin/dispersal routes). 

I have made some tests using both simultaneously and it seems important because the adjacency matrices constrain whether areas can be connected at all during a given interval, whereas the dispersal multipliers allow the user to test biologically distinct dispersal scenarios within those paleogeographic constraints.  

Thanks again Ivan for your help

Best regards,
Emanuel

To view this discussion visit https://groups.google.com/d/msgid/biogeobears/CAAsMhYqPDKqiTnyqw50xSnAvMZY0t7w9LwWNY7-pp1ovfbBWow%40mail.gmail.com.