check.mix: Check whether chains of a fitted evorates model adequately...

In bstaggmartin/backwards-BM-simulator: Evolving Rate Models of Continuous Trait Evolution

Check whether chains of a fitted evorates model adequately converged

Description

This function calculates and uses the Rhat diagnostic to check for adequate convergence, or mixing, of Hamiltonian Monte Carlo chains given a fitted evorates model. Rhat's should generally only be slightly above 1 for all parameters. Notably, this function takes all chains into account to calculate Rhat's–diagnostics param_block arrays instead calculate Rhat's for each individual chain.

Usage

check.mix(fit, printlen = 6)

Arguments

fit	An object of class "evorates_fit".
printlen	Sets the maximum number of parameter names to report for each category of Rhat values: "good", "okay", and "bad". For no maximum, set this to 0, a negative number, or infinity.

Details

Rhat's are calculated using rstan::Rhat(). More details on this diagnostic can be found in the associated documentation: Rhat. Empirically-calculated Rhat's can sometimes be slightly below 1 when chains mix well. Even though 1.01 or below is ideal, if Rhat's for all parameters are around 1.05 or below, the inferred posterior distribution should be pretty reliable in my experience. In any case, users should make sure all the Rhat's are okay to good before running combine.chains(fit) or trusting the estimated posterior distribution!

In the case of high model misspecification (i.e., an evorates model does not fit the data well), one could theoretically get high Rhat's and running longer chains could result in virtually no benefit. I have yet to see this situation, but it is certainly possible.

This function complements Rhat calculation via grapes-diagnostics-grapes since it calculates Rhat's for all parameters using all chains. It thus is tailored for diagnosing the "healthiness" of a fit as a whole, rather than the healthiness of chains for individual parameters/chains.

Value

Prints to console how many parameters exhibit "good" (Rhat <= 1.01), "okay" (1.01 < Rhat <= 1.05), and "bad" (Rhat > 1.05) mixing, along with some general advice depending on the worst observed Rhat. If the number of parameters in a particular category is less than or equal to print.len, the function will also print the names of the parameter in that category.

Invisibly, the function also returns a list of all calculated Rhat's, named by their corresponding parameters and split into good, okay, and bad categories.

See Also

Rhat for the function used to calculate Rhat's. Also, grapes-diagnostics-grapes for a function tailored towards diagnosing individual chains for specific parameters.

Other evorates diagnostic functions: check.ess()

Examples

#get whale/dolphin evorates fit
data("cet_fit")

#check how well the chains converged
Rhats <- check.mix(cet_fit) #yay!
Rhats #invisibly returned object

#note the differences between...
Rhats; cet_fit %diagnostics% list(".", "Rhat")

#generally you should make sure the Rhat's are all okay before running combine.chains()

#fit a new model (poorly--please never set your warmup this low in practice)
#this is just to see what happens with bad mixing
#this will take just a little bit to run...
fit <- fit.evorates(cet_fit$call$tree, cet_fit$call$trait.data, 
                    constrain.Rsig2 = TRUE,
                    warmup = 10, iter = 1010)
Rhats <- check.mix(fit) #not yay

bstaggmartin/backwards-BM-simulator documentation built on June 3, 2024, 5:51 p.m.