Alternatively, we can take a different Monte Carlo approach to
Equation REF  and sample from \(p(D ^{\prime } \,|\, M_{}, D)\)  by simulating datasets.
For each posterior sample of the parameter values (conditional on all the data
under study) we can simply simulate a new dataset based on those parameter
values.
We can then compare the observed data (\(D\) ) to the sample of simulated
datasets (\(D\)  \(^{\prime }\) ) from the posterior predictive distribution.
In all but the most trivial phylogenetic datasets, it is not practical to
compare the counts of site patterns directly, because there are too many
possible patterns (e.g., four raised to the power of the number of tips for DNA
data).
Thus, we have to tolerate some loss of information by summarizing the data in
some way to reduce the dimensionality.
Once a summary statistic is chosen, perhaps the simplest way to evaluate the
fit of the model is to approximate the posterior predictive p-value by finding
the percentile of the statistic from the observed data out of the values of the
statistic calculated from the simulated datasets [1]}, [2]}.
[3]} explored this approach for phylogenetic models using
simulated data,
and found that a simple JC69 model [4]} was often rejected for data
simulated under more complex K2P [5]} and GTR [6]} models.
[7]} also used this approach to corroborate their findings
based on marginal likelihoods [8]} and cross validation that
allowing among-site variation in amino acid composition (i.e., the CAT model)
leads to a better fit.
