# Load the package (after installation, see above). library(optimx) # You need to have some version of optimx available # as it is a BioGeoBEARS dependency; however, if you # don't want to use optimx, and use optim() (from R core) # you can set: # BioGeoBEARS_run_object$use_optimx = FALSE # ...everything should work either way -- NJM 2014-01-08 library(FD) # for FD::maxent() (make sure this is up-to-date) library(snow) # (if you want to use multicore functionality; some systems/R versions prefer library(parallel), try either) library(parallel) library(BioGeoBEARS) ######################################################## # TO GET THE OPTIMX/OPTIM FIX, AND THE UPPASS FIX, # SOURCE THE REVISED FUNCTIONS WITH THESE COMMANDS ######################################################## #library(BioGeoBEARS) #source("http://phylo.wdfiles.com/local--files/biogeobears/cladoRcpp.R") # (needed now that traits model added; source FIRST!) #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_add_fossils_randomly_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_basics_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_calc_transition_matrices_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_classes_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_detection_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_DNA_cladogenesis_sim_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_extract_Qmat_COOmat_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_generics_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_models_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_on_multiple_trees_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_plots_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_readwrite_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_simulate_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_SSEsim_makePlots_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_SSEsim_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_stochastic_mapping_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_stratified_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_univ_model_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/calc_uppass_probs_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/calc_loglike_sp_v01.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/get_stratified_subbranch_top_downpass_likelihoods_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/runBSM_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/stochastic_map_given_inputs.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/summarize_BSM_tables_v1.R") #source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_traits_v1.R") # added traits model calc_loglike_sp = compiler::cmpfun(calc_loglike_sp_prebyte) # crucial to fix bug in uppass calculations calc_independent_likelihoods_on_each_branch = compiler::cmpfun(calc_independent_likelihoods_on_each_branch_prebyte) # slight speedup hopefully ####################################################### # SETUP: YOUR WORKING DIRECTORY ####################################################### # You will need to set your working directory to match your local system # Note these very handy functions! # Command "setwd(x)" sets your working directory # Command "getwd()" gets your working directory and tells you what it is. # Command "list.files()" lists the files in your working directory # To get help on any command, use "?". E.g., "?list.files" # Set your working directory for output files # default here is your home directory ("~") # Change this as you like wd = "E:/dinosaurs/WIP/Lishulong wangi" setwd(wd) # Double-check your working directory with getwd() getwd() ####################################################### # SETUP: Extension data directory ####################################################### # When R packages contain extra files, they are stored in the "extdata" directory # inside the installed package. # # BioGeoBEARS contains various example files and scripts in its extdata directory. # # Each computer operating system might install BioGeoBEARS in a different place, # depending on your OS and settings. # # However, you can find the extdata directory like this: extdata_dir = np(system.file("extdata", package="BioGeoBEARS")) extdata_dir list.files(extdata_dir) # "system.file" looks in the directory of a specified package (in this case BioGeoBEARS) # The function "np" is just a shortcut for normalizePath(), which converts the # path to the format appropriate for your system (e.g., Mac/Linux use "/", but # Windows uses "\\", if memory serves). # Even when using your own data files, you should KEEP these commands in your # script, since the plot_BioGeoBEARS_results function needs a script from the # extdata directory to calculate the positions of "corners" on the plot. This cannot # be made into a straight up BioGeoBEARS function because it uses C routines # from the package APE which do not pass R CMD check for some reason. ####################################################### # SETUP: YOUR TREE FILE AND GEOGRAPHY FILE ####################################################### # Example files are given below. To run your own data, # make the below lines point to your own files, e.g. # trfn = "/mydata/frogs/frogBGB/tree.newick" # geogfn = "/mydata/frogs/frogBGB/geog.data" ####################################################### # Phylogeny file # Notes: # 1. Must be binary/bifurcating: no polytomies # 2. No negative branchlengths (e.g. BEAST MCC consensus trees sometimes have negative branchlengths) # 3. Be careful of very short branches, as BioGeoBEARS will interpret ultrashort branches as direct ancestors # 4. You can use non-ultrametric trees, but BioGeoBEARS will interpret any tips significantly below the # top of the tree as fossils! This is only a good idea if you actually do have fossils in your tree, # as in e.g. Wood, Matzke et al. (2013), Systematic Biology. # 5. The default settings of BioGeoBEARS make sense for trees where the branchlengths are in units of # millions of years, and the tree is 1-1000 units tall. If you have a tree with a total height of # e.g. 0.00001, you will need to adjust e.g. the max values of d and e, or (simpler) multiply all # your branchlengths to get them into reasonable units. # 6. DON'T USE SPACES IN SPECIES NAMES, USE E.G. "_" ####################################################### # This is the example Newick file for Hawaiian JL # (from Ree & Smith 2008) # "trfn" = "tree file name" trfn = "time-calibrated tree.newick" # Look at the raw Newick file: moref(trfn) # Look at your phylogeny: tr = read.tree(trfn) tr ####################################################### ####################################################### library(ape) library(geoscale) library(strap) treefn = "1st MPT(four high-level taxonomic OTUs pruned).phy" subtr <- read.tree(treefn) subtr agesfn = "Taxon age_revised.txt" ag <- read.table(agesfn) ag calibratedtr <- DatePhylo(subtr, ag, rlen=5, method="equal", add.terminal=FALSE) calibratedtr ######################################################## plot(calibratedtr) title("Tree from Jia Lei") axisPhylo(root.time=calibratedtr$root.time) # plots timescale ####################################################### # Geography file # Notes: # 1. This is a PHLYIP-formatted file. This means that in the # first line, # - the 1st number equals the number of rows (species) # - the 2nd number equals the number of columns (number of areas) # 2. This is the same format used for C++ LAGRANGE geography files. # 3. All names in the geography file must match names in the phylogeny file. # 4. DON'T USE SPACES IN SPECIES NAMES, USE E.G. "_" # 5. Operational taxonomic units (OTUs) should ideally be phylogenetic lineages, # i.e. genetically isolated populations. These may or may not be identical # with species. You would NOT want to just use specimens, as each specimen # automatically can only live in 1 area, which will typically favor DEC+J # models. This is fine if the species/lineages really do live in single areas, # but you wouldn't want to assume this without thinking about it at least. # In summary, you should collapse multiple specimens into species/lineages if # data indicates they are the same genetic population. ###################################################### # This is the example geography file for Hawaiian JL # (from Ree & Smith 2008) geogfn = "Taxon geographic range_revised.txt" # Look at the raw geography text file: moref(geogfn) # Look at your geographic range data: tipranges = getranges_from_LagrangePHYLIP(lgdata_fn=geogfn) tipranges # Set the maximum number of areas any species may occupy; this cannot be larger # than the number of areas you set up, but it can be smaller. max_range_size = 8 #################################################### #################################################### # KEY HINT: The number of states (= number of different possible geographic ranges) # depends on (a) the number of areas and (b) max_range_size. # If you have more than about 500-600 states, the calculations will get REALLY slow, # since the program has to exponentiate a matrix of e.g. 600x600. Often the computer # will just sit there and crunch, and never get through the calculation of the first # likelihood. # # (this is also what is usually happening when LAGRANGE hangs: you have too many states!) # # To check the number of states for a given number of ranges, try: numstates_from_numareas(numareas=4, maxareas=4, include_null_range=TRUE) numstates_from_numareas(numareas=4, maxareas=4, include_null_range=FALSE) numstates_from_numareas(numareas=4, maxareas=3, include_null_range=TRUE) numstates_from_numareas(numareas=4, maxareas=2, include_null_range=TRUE) # Large numbers of areas have problems: numstates_from_numareas(numareas=10, maxareas=10, include_null_range=TRUE) # ...unless you limit the max_range_size: numstates_from_numareas(numareas=10, maxareas=2, include_null_range=TRUE) #################################################### #################################################### ####################################################### # DEC AND DEC+J ANALYSIS ####################################################### ####################################################### # Run DEC ####################################################### BioGeoBEARS_run_object = define_BioGeoBEARS_run() BioGeoBEARS_run_object$force_sparse=FALSE # sparse=FALSE causes pathology & isn't much faster at this scale BioGeoBEARS_run_object$speedup=TRUE # shorcuts to speed ML search; use FALSE if worried (e.g. >3 params) BioGeoBEARS_run_object$calc_ancprobs=TRUE # get ancestral states from optim run # Set up a time-stratified analysis # (un-comment to use; see example files in extdata_dir, # and BioGeoBEARS google group posts for further hints) BioGeoBEARS_run_object$timesfn = "Time periods file.txt" BioGeoBEARS_run_object$dispersal_multipliers_fn = "Dispersal multiplier matrices (harsh).txt" #BioGeoBEARS_run_object$areas_allowed_fn = "areas_allowed.txt" #BioGeoBEARS_run_object$distsfn = "distances_matrix.txt" # See notes on the distances model on PhyloWiki's BioGeoBEARS updates page. # Input the maximum range size BioGeoBEARS_run_object$max_range_size = max_range_size # Multicore processing if desired BioGeoBEARS_run_object$num_cores_to_use=4 # (use more cores to speed it up; this requires # library(parallel) and/or library(snow). Parellel, # is default on Macs in R 3.0+, but apparently still # has to be typed on Windows machines. Note: apparently # parallel works on Mac command-line R, but not R.app. # BioGeoBEARS checks for this and resets to 1 # core with R.app) # Sparse matrix exponentiation is an option for huge numbers of ranges/states (600+) # I have experimented with sparse matrix exponentiation in EXPOKIT/rexpokit, # but the results are imprecise and so I haven't explored it further. # In a Bayesian analysis, it might work OK, but the ML point estimates are # not identical. # Also, I have not implemented all functions to work with force_sparse=TRUE. # Volunteers are welcome to work on it!! BioGeoBEARS_run_object$force_sparse=FALSE # Give BioGeoBEARS the location of the geography text file BioGeoBEARS_run_object$geogfn = geogfn # Give BioGeoBEARS the location of the phylogeny Newick file BioGeoBEARS_run_object$trfn = trfn # This function loads the dispersal multiplier matrix etc. from the text files into the model object. Required for these to work! # (It also runs some checks on these inputs for certain errors.) BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object) # Divide the tree up by timeperiods/strata (uncomment this for stratified analysis) BioGeoBEARS_run_object = section_the_tree(inputs=BioGeoBEARS_run_object, make_master_table=TRUE, plot_pieces=FALSE, cut_fossils=FALSE, fossils_older_than = 0) # The stratified tree is described in this table: #BioGeoBEARS_run_object$master_table # Good default settings to get ancestral states BioGeoBEARS_run_object$return_condlikes_table = TRUE BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE BioGeoBEARS_run_object$calc_ancprobs = TRUE # Set up DEC model # (nothing to do; defaults) # Look at the BioGeoBEARS_run_object; it's just a list of settings etc. BioGeoBEARS_run_object # This contains the model object BioGeoBEARS_run_object$BioGeoBEARS_model_object # This table contains the parameters of the model BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table # Run this to check inputs. Read the error messages if you get them! check_BioGeoBEARS_run(BioGeoBEARS_run_object) # For a slow analysis, run once, then set runslow=FALSE to just # load the saved result. runslow = TRUE resfn = "JL_DEC_M3_stratified_v1.Rdata" if (runslow) { res = bears_optim_run(BioGeoBEARS_run_object) res save(res, file=resfn) resDEC = res } else { # Loads to "res" load(resfn) resDEC = res } ####################################################### # Run DEC+J ####################################################### BioGeoBEARS_run_object = define_BioGeoBEARS_run() BioGeoBEARS_run_object$trfn = trfn BioGeoBEARS_run_object$geogfn = geogfn BioGeoBEARS_run_object$max_range_size = max_range_size # Set up the stratified part BioGeoBEARS_run_object$timesfn = "Time periods file.txt" BioGeoBEARS_run_object$dispersal_multipliers_fn = "Dispersal multiplier matrices (harsh).txt" #BioGeoBEARS_run_object$areas_allowed_fn = "areas_allowed.txt" #BioGeoBEARS_run_object$distsfn = "distances_matrix.txt" # See notes on the distances model on PhyloWiki's BioGeoBEARS updates page. BioGeoBEARS_run_object$use_optimx = TRUE BioGeoBEARS_run_object$num_cores_to_use=4 BioGeoBEARS_run_object$force_sparse=FALSE # sparse=FALSE causes pathology & isn't much faster at this scale BioGeoBEARS_run_object$speedup=TRUE # shorcuts to speed ML search; use FALSE if worried (e.g. >3 params) BioGeoBEARS_run_object$calc_ancprobs=TRUE # get ancestral states from optim run # This function loads the dispersal multiplier matrix etc. from the text files into the model object. Required for these to work! # (It also runs some checks on these inputs for certain errors.) BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object) # Divide the tree up by timeperiods/strata (uncomment this for stratified analysis) BioGeoBEARS_run_object = section_the_tree(inputs=BioGeoBEARS_run_object, make_master_table=TRUE, plot_pieces=FALSE, cut_fossils=FALSE,fossils_older_than = 0) # The stratified tree is described in this table: #BioGeoBEARS_run_object$master_table # Good default settings to get ancestral states BioGeoBEARS_run_object$return_condlikes_table = TRUE BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE BioGeoBEARS_run_object$calc_ancprobs = TRUE # Set up DEC+J model # Get the ML parameter values from the 2-parameter nested model # (this will ensure that the 3-parameter model always does at least as good) dstart = resDEC$outputs@params_table["d","est"] estart = resDEC$outputs@params_table["e","est"] jstart = 0.0001 # Input starting values for d, e BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","init"] = dstart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","est"] = dstart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","init"] = estart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","est"] = estart # Add j as a free parameter BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = jstart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = jstart check_BioGeoBEARS_run(BioGeoBEARS_run_object) resfn = "JL_DEC+J_M3_stratified_v1.Rdata" runslow = TRUE if (runslow) { #sourceall("/Dropbox/_njm/__packages/BioGeoBEARS_setup/") res = bears_optim_run(BioGeoBEARS_run_object) res save(res, file=resfn) resDECj = res } else { # Loads to "res" load(resfn) resDECj = res } ####################################################### # PDF plots ####################################################### pdffn = "JL_DEC_vs_DEC+J_M3_stratified_v1.pdf" pdf(pdffn, width=12, height=24) ####################################################### # Plot ancestral states - DEC ####################################################### analysis_titletxt ="BioGeoBEARS DEC on JL M3_stratified" # Setup results_object = resDEC scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS")) res2 = plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.45, splitcex=0.45, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") ####################################################### # Plot ancestral states - DECJ ####################################################### analysis_titletxt ="BioGeoBEARS DEC+J on JL M3_stratified" # Setup results_object = resDECj scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS")) # States res1 = plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") # Pie chart plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.45, splitcex=0.45, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") dev.off() # Turn off PDF cmdstr = paste("open ", pdffn, sep="") system(cmdstr) # Plot it ####################################################### # DIVALIKE AND DIVALIKE+J ANALYSIS ####################################################### # NOTE: The BioGeoBEARS "DIVA" model is not identical with # Ronquist (1997)'s parsimony DIVA. It is a likelihood # interpretation of DIVA, constructed by modelling DIVA's # processes the way DEC does, but only allowing the # processes DIVA allows (widespread vicariance: yes; subset # sympatry: no; see Ronquist & Sanmartin 2011, Figure 4). # # I thus now call the model "DIVALIKE", and you should also. ;-) # ####################################################### # Run DIVALIKE ####################################################### BioGeoBEARS_run_object = define_BioGeoBEARS_run() BioGeoBEARS_run_object$trfn = trfn BioGeoBEARS_run_object$geogfn = geogfn BioGeoBEARS_run_object$max_range_size = max_range_size # Set up the stratified part BioGeoBEARS_run_object$timesfn = "Time periods file.txt" BioGeoBEARS_run_object$dispersal_multipliers_fn = "Dispersal multiplier matrices (harsh).txt" #BioGeoBEARS_run_object$areas_allowed_fn = "areas_allowed.txt" #BioGeoBEARS_run_object$distsfn = "distances_matrix.txt" # See notes on the distances model on PhyloWiki's BioGeoBEARS updates page. BioGeoBEARS_run_object$use_optimx = TRUE BioGeoBEARS_run_object$num_cores_to_use=4 BioGeoBEARS_run_object$force_sparse=FALSE # sparse=FALSE causes pathology & isn't much faster at this scale BioGeoBEARS_run_object$speedup=TRUE # shorcuts to speed ML search; use FALSE if worried (e.g. >3 params) BioGeoBEARS_run_object$calc_ancprobs=TRUE # get ancestral states from optim run # This function loads the dispersal multiplier matrix etc. from the text files into the model object. Required for these to work! # (It also runs some checks on these inputs for certain errors.) BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object) # Divide the tree up by timeperiods/strata (uncomment this for stratified analysis) BioGeoBEARS_run_object = section_the_tree(inputs=BioGeoBEARS_run_object, make_master_table=TRUE, plot_pieces=FALSE, cut_fossils=FALSE,fossils_older_than = 0) # The stratified tree is described in this table: #BioGeoBEARS_run_object$master_table # Good default settings to get ancestral states BioGeoBEARS_run_object$return_condlikes_table = TRUE BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE BioGeoBEARS_run_object$calc_ancprobs = TRUE # Set up DIVALIKE model # 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 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 # No jump dispersal/founder-event speciation # BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free" # BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = 0.01 # BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = 0.01 check_BioGeoBEARS_run(BioGeoBEARS_run_object) runslow = TRUE resfn = "JL_DIVALIKE_M3_stratified_v1.Rdata" if (runslow) { res = bears_optim_run(BioGeoBEARS_run_object) res save(res, file=resfn) resDIVALIKE = res } else { # Loads to "res" load(resfn) resDIVALIKE = res } ####################################################### # Run DIVALIKE+J ####################################################### BioGeoBEARS_run_object = define_BioGeoBEARS_run() BioGeoBEARS_run_object$trfn = trfn BioGeoBEARS_run_object$geogfn = geogfn BioGeoBEARS_run_object$max_range_size = max_range_size # Set up the stratified part BioGeoBEARS_run_object$timesfn = "Time periods file.txt" BioGeoBEARS_run_object$dispersal_multipliers_fn = "Dispersal multiplier matrices (harsh).txt" #BioGeoBEARS_run_object$areas_allowed_fn = "areas_allowed.txt" #BioGeoBEARS_run_object$distsfn = "distances_matrix.txt" # See notes on the distances model on PhyloWiki's BioGeoBEARS updates page. BioGeoBEARS_run_object$use_optimx = TRUE BioGeoBEARS_run_object$num_cores_to_use=4 BioGeoBEARS_run_object$force_sparse=FALSE # sparse=FALSE causes pathology & isn't much faster at this scale BioGeoBEARS_run_object$speedup=TRUE # shorcuts to speed ML search; use FALSE if worried (e.g. >3 params) BioGeoBEARS_run_object$calc_ancprobs=TRUE # get ancestral states from optim run # This function loads the dispersal multiplier matrix etc. from the text files into the model object. Required for these to work! # (It also runs some checks on these inputs for certain errors.) BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object) # Divide the tree up by timeperiods/strata (uncomment this for stratified analysis) BioGeoBEARS_run_object = section_the_tree(inputs=BioGeoBEARS_run_object, make_master_table=TRUE, plot_pieces=FALSE, cut_fossils=FALSE,fossils_older_than = 0) # The stratified tree is described in this table: #BioGeoBEARS_run_object$master_table # Good default settings to get ancestral states BioGeoBEARS_run_object$return_condlikes_table = TRUE BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE BioGeoBEARS_run_object$calc_ancprobs = TRUE # Set up DIVALIKE+J model # Get the ML parameter values from the 2-parameter nested model # (this will ensure that the 3-parameter model always does at least as good) dstart = resDIVALIKE$outputs@params_table["d","est"] estart = resDIVALIKE$outputs@params_table["e","est"] jstart = 0.0001 # Input starting values for d, e BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","init"] = dstart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","est"] = dstart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","init"] = estart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","est"] = estart # 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 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 # Add jump dispersal/founder-event speciation BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = jstart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = jstart # Under DIVALIKE+J, the max of "j" should be 2, not 3 (as is default in DEC+J) BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","min"] = 0.00001 BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","max"] = 1.99999 check_BioGeoBEARS_run(BioGeoBEARS_run_object) resfn = "JL_DIVALIKE+J_M3_stratified_v1.Rdata" runslow = TRUE if (runslow) { #sourceall("/Dropbox/_njm/__packages/BioGeoBEARS_setup/") res = bears_optim_run(BioGeoBEARS_run_object) res save(res, file=resfn) resDIVALIKEj = res } else { # Loads to "res" load(resfn) resDIVALIKEj = res } pdffn = "JL_DIVALIKE_vs_DIVALIKE+J_M3_stratified_v1.pdf" pdf(pdffn, width=12, height=24) ####################################################### # Plot ancestral states - DIVALIKE ####################################################### analysis_titletxt ="BioGeoBEARS DIVALIKE on JL M3_stratified" # Setup results_object = resDIVALIKE scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS")) # States res2 = plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") # Pie chart plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.45, splitcex=0.45, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") ####################################################### # Plot ancestral states - DIVALIKE+J ####################################################### analysis_titletxt ="BioGeoBEARS DIVALIKE+J on JL M3_stratified" # Setup results_object = resDIVALIKEj scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS")) # States res1 = plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") # Pie chart plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.45, splitcex=0.45, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") dev.off() cmdstr = paste("open ", pdffn, sep="") system(cmdstr) ####################################################### # BAYAREALIKE AND BAYAREALIKE+J ANALYSIS ####################################################### # NOTE: As with DIVA, the BioGeoBEARS BayArea-like model is # not identical with the full Bayesian model implemented # in the "BayArea" program of Landis et al. (2013). # # Instead, this is a simplified likelihood interpretation # of the model. Basically, in BayArea and BioGeoBEARS-BAYAREALIKE, # "d" and "e" work like they do in the DEC model of Lagrange # (and BioGeoBEARS), and then BayArea's cladogenesis assumption # (which is that nothing in particular happens at cladogenesis) is # replicated by BioGeoBEARS. # # This leaves out 3 important things that are in BayArea: # 1. Distance dependence (you can add this with a distances # matrix + the "x" parameter in BioGeoBEARS, however) # 2. A correction for disallowing "e" events that drive # a species extinct (a null geographic range) # 3. The neat Bayesian sampling of histories, which allows # analyses on large numbers of areas. # # The main purpose of having a "BAYAREALIKE" model is # to test the importance of the cladogenesis model on # particular datasets. Does it help or hurt the data # likelihood if there is no special cladogenesis process? # # I thus now call the model "BAYAREALIKE", and you should also. ;-) # ####################################################### # Run BAYAREALIKE ####################################################### BioGeoBEARS_run_object = define_BioGeoBEARS_run() BioGeoBEARS_run_object$trfn = trfn BioGeoBEARS_run_object$geogfn = geogfn BioGeoBEARS_run_object$max_range_size = max_range_size # Set up the stratified part BioGeoBEARS_run_object$timesfn = "Time periods file.txt" BioGeoBEARS_run_object$dispersal_multipliers_fn = "Dispersal multiplier matrices (harsh).txt" #BioGeoBEARS_run_object$areas_allowed_fn = "areas_allowed.txt" #BioGeoBEARS_run_object$distsfn = "distances_matrix.txt" # See notes on the distances model on PhyloWiki's BioGeoBEARS updates page. BioGeoBEARS_run_object$use_optimx = TRUE BioGeoBEARS_run_object$num_cores_to_use=4 BioGeoBEARS_run_object$force_sparse=FALSE # sparse=FALSE causes pathology & isn't much faster at this scale BioGeoBEARS_run_object$speedup=TRUE # shorcuts to speed ML search; use FALSE if worried (e.g. >3 params) BioGeoBEARS_run_object$calc_ancprobs=TRUE # get ancestral states from optim run # This function loads the dispersal multiplier matrix etc. from the text files into the model object. Required for these to work! # (It also runs some checks on these inputs for certain errors.) BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object) # Divide the tree up by timeperiods/strata (uncomment this for stratified analysis) BioGeoBEARS_run_object = section_the_tree(inputs=BioGeoBEARS_run_object, make_master_table=TRUE, plot_pieces=FALSE, cut_fossils=FALSE,fossils_older_than = 0) # The stratified tree is described in this table: #BioGeoBEARS_run_object$master_table # Good default settings to get ancestral states BioGeoBEARS_run_object$return_condlikes_table = TRUE BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE BioGeoBEARS_run_object$calc_ancprobs = TRUE # Set up BAYAREALIKE model # No 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 # No vicariance BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","type"] = "fixed" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","init"] = 0.0 BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","est"] = 0.0 # No jump dispersal/founder-event speciation # BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free" # BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = 0.01 # BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = 0.01 # Adjust linkage between parameters BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ysv","type"] = "1-j" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ys","type"] = "ysv*1/1" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["y","type"] = "1-j" # Only sympatric/range-copying (y) events allowed, and with # exact copying (both descendants always the same size as the ancestor) BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","type"] = "fixed" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","init"] = 0.9999 BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","est"] = 0.9999 # Check the inputs check_BioGeoBEARS_run(BioGeoBEARS_run_object) runslow = TRUE resfn = "JL_BAYAREALIKE_M3_stratified_v1.Rdata" if (runslow) { res = bears_optim_run(BioGeoBEARS_run_object) res save(res, file=resfn) resBAYAREALIKE = res } else { # Loads to "res" load(resfn) resBAYAREALIKE = res } ####################################################### # Run BAYAREALIKE+J ####################################################### BioGeoBEARS_run_object = define_BioGeoBEARS_run() BioGeoBEARS_run_object$trfn = trfn BioGeoBEARS_run_object$geogfn = geogfn BioGeoBEARS_run_object$max_range_size = max_range_size # Set up the stratified part BioGeoBEARS_run_object$timesfn = "Time periods file.txt" BioGeoBEARS_run_object$dispersal_multipliers_fn = "Dispersal multiplier matrices (harsh).txt" #BioGeoBEARS_run_object$areas_allowed_fn = "areas_allowed.txt" #BioGeoBEARS_run_object$distsfn = "distances_matrix.txt" # See notes on the distances model on PhyloWiki's BioGeoBEARS updates page. BioGeoBEARS_run_object$use_optimx = TRUE BioGeoBEARS_run_object$num_cores_to_use=4 BioGeoBEARS_run_object$force_sparse=FALSE # sparse=FALSE causes pathology & isn't much faster at this scale BioGeoBEARS_run_object$speedup=TRUE # shorcuts to speed ML search; use FALSE if worried (e.g. >3 params) BioGeoBEARS_run_object$calc_ancprobs=TRUE # get ancestral states from optim run # This function loads the dispersal multiplier matrix etc. from the text files into the model object. Required for these to work! # (It also runs some checks on these inputs for certain errors.) BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object) # Divide the tree up by timeperiods/strata (uncomment this for stratified analysis) BioGeoBEARS_run_object = section_the_tree(inputs=BioGeoBEARS_run_object, make_master_table=TRUE, plot_pieces=FALSE, cut_fossils=FALSE,fossils_older_than = 0) # The stratified tree is described in this table: #BioGeoBEARS_run_object$master_table # Good default settings to get ancestral states BioGeoBEARS_run_object$return_condlikes_table = TRUE BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE BioGeoBEARS_run_object$calc_ancprobs = TRUE # Set up BAYAREALIKE+J model # Get the ML parameter values from the 2-parameter nested model # (this will ensure that the 3-parameter model always does at least as good) dstart = resBAYAREALIKE$outputs@params_table["d","est"] estart = resBAYAREALIKE$outputs@params_table["e","est"] jstart = 0.0001 # Input starting values for d, e BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","init"] = dstart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","est"] = dstart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","init"] = estart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","est"] = estart # No 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 # No vicariance BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","type"] = "fixed" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","init"] = 0.0 BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["v","est"] = 0.0 # *DO* allow jump dispersal/founder-event speciation (set the starting value close to 0) BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","type"] = "free" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","init"] = jstart BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","est"] = jstart # Under BAYAREALIKE+J, the max of "j" should be 1, not 3 (as is default in DEC+J) or 2 (as in DIVALIKE+J) BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","max"] = 0.99999 # Adjust linkage between parameters BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ysv","type"] = "1-j" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["ys","type"] = "ysv*1/1" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["y","type"] = "1-j" # Only sympatric/range-copying (y) events allowed, and with # exact copying (both descendants always the same size as the ancestor) BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","type"] = "fixed" BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","init"] = 0.9999 BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["mx01y","est"] = 0.9999 # NOTE (NJM, 2014-04): BAYAREALIKE+J seems to crash on some computers, usually Windows # machines. I can't replicate this on my Mac machines, but it is almost certainly # just some precision under-run issue, when optim/optimx tries some parameter value # just below zero. The "min" and "max" options on each parameter are supposed to # prevent this, but apparently optim/optimx sometimes go slightly beyond # these limits. Anyway, if you get a crash, try raising "min" and lowering "max" # slightly for each parameter: BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","min"] = 0.0000001 BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["d","max"] = 4.9999999 BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","min"] = 0.0000001 BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["e","max"] = 4.9999999 BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","min"] = 0.00001 BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table["j","max"] = 0.99999 check_BioGeoBEARS_run(BioGeoBEARS_run_object) resfn = "JL_BAYAREALIKE+J_M3_stratified_v1.Rdata" runslow = TRUE if (runslow) { res = bears_optim_run(BioGeoBEARS_run_object) res save(res, file=resfn) resBAYAREALIKEj = res } else { # Loads to "res" load(resfn) resBAYAREALIKEj = res } pdffn = "JL_BAYAREALIKE_vs_BAYAREALIKE+J_M3_stratified_v1.pdf" pdf(pdffn, width=12, height=24) ####################################################### # Plot ancestral states - BAYAREALIKE ####################################################### analysis_titletxt ="BioGeoBEARS BAYAREALIKE on JL M3_stratified" # Setup results_object = resBAYAREALIKE scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS")) # States res2 = plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") # Pie chart plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.45, splitcex=0.45, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") ####################################################### # Plot ancestral states - BAYAREALIKE+J ####################################################### analysis_titletxt ="BioGeoBEARS BAYAREALIKE+J on JL M3_stratified" # Setup results_object = resBAYAREALIKEj scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS")) # States res1 = plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") # Pie chart plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.45, splitcex=0.45, titlecex=0.8, plotsplits=FALSE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) trheight = max((prt(tr, printflag=FALSE))$time_bp) vertical_lines = (trheight - c(0, res$inputs$timeperiods)) abline(v=vertical_lines, lty="dashed") dev.off() cmdstr = paste("open ", pdffn, sep="") system(cmdstr) ######################################################################### ######################################################################### ######################################################################### ######################################################################### # # CALCULATE SUMMARY STATISTICS TO COMPARE # DEC, DEC+J, DIVALIKE, DIVALIKE+J, BAYAREALIKE, BAYAREALIKE+J # ######################################################################### ######################################################################### ######################################################################### ######################################################################### ######################################################################### ######################################################################### # REQUIRED READING: # # Practical advice / notes / basic principles on statistical model # comparison in general, and in BioGeoBEARS: # http://phylo.wikidot.com/advice-on-statistical-model-comparison-in-biogeobears ######################################################################### ######################################################################### # Set up empty tables to hold the statistical results restable = NULL teststable = NULL ####################################################### # Statistics -- DEC vs. DEC+J ####################################################### # We have to extract the log-likelihood differently, depending on the # version of optim/optimx LnL_2 = get_LnL_from_BioGeoBEARS_results_object(resDEC) LnL_1 = get_LnL_from_BioGeoBEARS_results_object(resDECj) numparams1 = 3 numparams2 = 2 stats = AICstats_2models(LnL_1, LnL_2, numparams1, numparams2) stats # DEC, null model for Likelihood Ratio Test (LRT) res2 = extract_params_from_BioGeoBEARS_results_object(results_object=resDEC, returnwhat="table", addl_params=c("j"), paramsstr_digits=4) # DEC+J, alternative model for Likelihood Ratio Test (LRT) res1 = extract_params_from_BioGeoBEARS_results_object(results_object=resDECj, returnwhat="table", addl_params=c("j"), paramsstr_digits=4) # The null hypothesis for a Likelihood Ratio Test (LRT) is that two models # confer the same likelihood on the data. See: Brian O'Meara's webpage: # http://www.brianomeara.info/tutorials/aic # ...for an intro to LRT, AIC, and AICc rbind(res2, res1) tmp_tests = conditional_format_table(stats) restable = rbind(restable, res2, res1) teststable = rbind(teststable, tmp_tests) ####################################################### # Statistics -- DIVALIKE vs. DIVALIKE+J ####################################################### # We have to extract the log-likelihood differently, depending on the # version of optim/optimx LnL_2 = get_LnL_from_BioGeoBEARS_results_object(resDIVALIKE) LnL_1 = get_LnL_from_BioGeoBEARS_results_object(resDIVALIKEj) numparams1 = 3 numparams2 = 2 stats = AICstats_2models(LnL_1, LnL_2, numparams1, numparams2) stats # DIVALIKE, null model for Likelihood Ratio Test (LRT) res2 = extract_params_from_BioGeoBEARS_results_object(results_object=resDIVALIKE, returnwhat="table", addl_params=c("j"), paramsstr_digits=4) # DIVALIKE+J, alternative model for Likelihood Ratio Test (LRT) res1 = extract_params_from_BioGeoBEARS_results_object(results_object=resDIVALIKEj, returnwhat="table", addl_params=c("j"), paramsstr_digits=4) rbind(res2, res1) conditional_format_table(stats) tmp_tests = conditional_format_table(stats) restable = rbind(restable, res2, res1) teststable = rbind(teststable, tmp_tests) ####################################################### # Statistics -- BAYAREALIKE vs. BAYAREALIKE+J ####################################################### # We have to extract the log-likelihood differently, depending on the # version of optim/optimx LnL_2 = get_LnL_from_BioGeoBEARS_results_object(resBAYAREALIKE) LnL_1 = get_LnL_from_BioGeoBEARS_results_object(resBAYAREALIKEj) numparams1 = 3 numparams2 = 2 stats = AICstats_2models(LnL_1, LnL_2, numparams1, numparams2) stats # BAYAREALIKE, null model for Likelihood Ratio Test (LRT) res2 = extract_params_from_BioGeoBEARS_results_object(results_object=resBAYAREALIKE, returnwhat="table", addl_params=c("j"), paramsstr_digits=4) # BAYAREALIKE+J, alternative model for Likelihood Ratio Test (LRT) res1 = extract_params_from_BioGeoBEARS_results_object(results_object=resBAYAREALIKEj, returnwhat="table", addl_params=c("j"), paramsstr_digits=4) rbind(res2, res1) conditional_format_table(stats) tmp_tests = conditional_format_table(stats) restable = rbind(restable, res2, res1) teststable = rbind(teststable, tmp_tests) ######################################################################### # RESULTS: DEC, DEC+J, DIVALIKE, DIVALIKE+J, BAYAREALIKE, BAYAREALIKE+J ######################################################################### teststable$alt = c("DEC+J", "DIVALIKE+J", "BAYAREALIKE+J") teststable$null = c("DEC", "DIVALIKE", "BAYAREALIKE") row.names(restable) = c("DEC", "DEC+J", "DIVALIKE", "DIVALIKE+J", "BAYAREALIKE", "BAYAREALIKE+J") # Look at the results!! restable teststable ####################################################### # Save the results tables for later -- check for e.g. # convergence issues ####################################################### # Loads to "restable" save(restable, file="restable_v1.Rdata") load(file="restable_v1.Rdata") # Loads to "teststable" save(teststable, file="teststable_v1.Rdata") load(file="teststable_v1.Rdata")