I read abstract only.
In your paper, the posterior predictive p  value (PPP) of the chi square goodness of fit are used.
And I just now try to implement it for my models.
In frequentist chi square is the following:
T(y,\theta) := \sum _{c=1,\cdots,C} \frac{\bigl\{y_c- \mathbb{E}[Y_c|\theta ] \}^2}{\mathbb{E}[Y_c|\theta ]}\\
where, y_c is a observed positive integer  of the  random variavle Y_c.
In frequentist method, to calculate the p value of chi square statistics,  it requires the notion of degree of freedom.
On the other hand, In Bayesian PPP, the degree of freedom is not required, since its definition  is  the following:
\iint \mathbb{I}[T(y,\theta) >T(y_{obs},\theta) ] p(y|\theta)p(\theta|y_{obs}) \, dy \, d\theta\\
in which there is no terms for the degree of freedom.
Is my interpretation correct? That is, in the Bayesian PPP, the notion of degree of freedom is no longer required.
Note that the  above PPP  is calculated by
\iint \mathbb{I}[T(y,\theta) >T(y_{obs},\theta) ] p(y|\theta)p(\theta|y_{obs}) \, dy \, d\theta\\
\approx \sum_i \int  \mathbb{I}[T(y,\theta_i) >T(y_{obs},\theta_i) ] p(y|\theta_i)\\
\approx \sum_i \sum_j   \mathbb{I}[T(y_i^j,\theta_i) >T(y_{obs},\theta_i) ]\\
where \theta_i is MCMC samples from posterior and  y_i^j, j=1,2,3,... is a samples drawn from likelihood at parameter \theta_i.
Moreover I think if C is larger, then the above T(y|\theta) is also larger,
and to avoid this, I think   I should use   the normalizing  multiple \frac{1}{C}  as follows:
T(y,\theta) := \frac{1}{C} \sum _{c=1,\cdots,C} \frac{\bigl\{y_c- \mathbb{E}[Y_c|\theta ] \}^2}{\mathbb{E}[Y_c|\theta ]}.
If I misunderstand about ppp, then please let me know.