BMhybExhaustive: Exhaustively evaluate models
In BMhyb: Comparative Methods for Phylogenetic Networks

Description Usage Arguments Details Value Examples

Fits all possible BMhyb models to your data.

1 2	BMhybExhaustive(phy.graph, traits, measurement.error = 0, ncores = max(c(1, parallel::detectCores() - 1), na.rm = TRUE), ...)

`phy.graph`	An ape::evonet object (a phylogeny stored in phylo format that also includes a reticulation matrix)
`traits`	A vector of trait values, with names equal to the names of taxa on the phylogeny
`measurement.error`	How much uncertainty there is in tip values; a single number is applied to all taxa, a vector is applied to the corresponding taxa
`ncores`	Number of cores to use. By default, uses parallel package to detect what's available and uses all but one.
`...`	All other parameters to pass to BMhyb (see ?BMhyb)

This takes an ape::evonet object. If all you have is a tree (an ape::phylo object), you can use CreateHybridlessEvonet() to convert the tree to an evonet object. You can then use the AddHybridization() function to add hybrid events to this object. Note that networks created in this way can, by chance, result in orders of nodes in the internal edge matrix that cause ape's reorder.phylo function to crash, which is called in many of the plot and write functions. You can still use the plot functions in this package, however.

This will return a list with one model result per element: you can plot these individually (see ?hybResult). By default, these results will include the information about uncertainty. We also compute a summary table so you can see the point estimates for each model and the likelihoods. It is often advisable to average across models, weighting each by its AICc weight, so this is also done automatically. We also return the single best model as an object for convenience, though for most users, we would suggest using the model average and looking at a set of fairly good models rather than look only at the single best one: there are often others that are nearly as good.

We do not expect large AIC difference between models unless you have a really large tree, and so you may get a warning if this happens. It is likely something has gone wrong with optimization. Look at all the models and examine for outliers. This issue can come up with certain combinations of networks and parameters (even, very rarely, in Brownian motion with no hybridization), where a step in the likelihood (inverting a matrix) does not yield a numerically stable result (the matrix is poorly conditioned). The 'likelihoods' in such cases are wrong, and they can look too good or too bad. Neither is ideal, but you should especially beware cases where the 'best' model has likelihoods much below some of the other models – you will often see bad parameter estimates, too. If you get this, do not believe the results – perhaps look at models with better condition.

To try to help with this, if one or more of the models has poor condition at the maximum likelihood estimate, we report this as it having an obvious problem. It is still returned in the results and the original.summary.df objects, but it is excluded from model averaging, the summary.df, and the best.model return (though note the ModelNumber column in summary.df allowing you to get the matching model in the results list). A model not having an obvious problem does *not* mean it worked well, just that it does not exhibit one particular problematic issue. Essentially we're saying, "This model does not have a lion eating its foot" – which suggest it's not unhealthy in that way, but doesn't mean there's not a crocodile eating its hand. User beware. Plotting the confidence using the plot functions can help.

Returns a list of objects of class BMhybResult (results), a summary data frame (summary.df) with parameter estimates and weights for all models where we do not see obvious problems, a summary data frame of all the models, whether or no they seemed to fail (original.summary.df), the model averaged result weighted by AICc weights of the unproblematic models (model.average), and the best unproblematic model (best.model).

## Not run: 
utils::data("cichlid")
traits.only <- cichlid$traits_and_SE$trait
names(traits.only) <- rownames(cichlid$traits_and_SE)
all.models <- BMhybExhaustive(phy.graph=cichlid$phy.graph, traits=traits.only)
print(all.models$summary.df)

## End(Not run)