Why does Bayesian p-value involve the parameters in addition to the data? On page 146 of Gelman's Bayesian Data Analysis, Gelman discusses Bayesian p-value as a way to check the fit of the model. The idea is to compare the observed data ($y$) with data that could have been generated by the model if we replicate the experiment ($y^{rep}$).
He defines Bayesian p-value as
$$
p_B = Pr(T(y^{rep}, \theta) \geq T(y, \theta) | y)
$$
I don't quite understand why it makes sense to have the test statistic be a function of the parameters, $\theta$. Indeed, if the goal is to "compare the observed data with data that could have been generated by the model", shouldn't the comparison be strictly between $y$ and $y^{rep}$?
For example, on the same page, Gelman provides an example where he checks the fit of a normal model. The test statistic is:
$$
T(y, \theta) = | y_{(61)} - \theta | - |y_{(6)} - \theta |
$$
where $\theta$ is the mean of the normal model. This test statistic is designed to ignore the model fit at the extreme tail, beyond the 6th and 61th order statistics.
Why don't we use the following test statistic instead, relying purely on the data?
$$
T(y, \theta) = | y_{(61)} - \bar y | - |y_{(6)} - \bar y |
$$
 A: The p-value is used to express the outcome in a test of a model and it's parameters (to test a hypothesis). Typically it relates to some statistic that measures a discrepancy (e.g. distance from the expected mean).
$P(T> t_{observed}|H_0)$
The probability that an observation of the statistic $T$ given the null hypothesis $H_0$ is larger than the observed value of the statistic $t_{observed}$.
The Bayesian p-value is much the same. It is used in the tests of the assumptions that are used to fit the model. And in this way is like a test for goodness of fit (for instance like a Pearson's Chi-squared test). The main difference with the Bayesian p-value/hypothesis is that the model is not with fixed parameters, but instead the parameters are variables themselves.


I don't quite understand why it makes sense to have the test statistic be a function of the parameters, $\theta$

The statistic that is used to compared the observation with the hypothetical model is often a pivotal quantity.

*

*You do not compare two observations. E.g. whether $Y_{rep} > Y$

*But instead you compare two observations in relation to the hypothetical model. E.g. whether $|Y_{rep}-\mu| > |Y-\mu|$. Whether the difference of the observation $Y$ from the mode of the model, is comparable to a likely random fluctuation $Y_{rep}$ or whether it is an unlikely value for which the we should consider the observation as a anomaly.


I hope this difference of observation and model makes intuitive sense. These statistics used for computing p-values can be somewhat arbitrary. The likelihood ratio test makes this a bit more formal.


Why don't we use the following test statistic instead, relying purely on the data?
$$
T(y, \theta) = | y_{(61)} - \bar y | - |y_{(6)} - \bar y |
$$

We do not make the statistic purely based on data because we want to test the data with relation to the model. It needs to be dependent on $\theta$ or otherwise it will not be a goodness of fit for $\theta$.
You can have in some way an expression for a statistic that tests data with "data".
$$
T(y, \theta) = | y_{(61)} - \hat y | - |y_{(6)} - \hat y |
$$
here $\hat y$ is the estimated value of $y$.
A: There is nothing stopping you from using test statistics based solely on the replicate data. 
The point in the example is that the model assumes $y_i$ is normally distributed around a mean $\theta$ not around $\overline{y}$. Thus the test statistic provided in Gelman et. al. tests something about the normality assumption whereas its not really clear what your test statistic is testing. 
