Calculates Posterior Predictive P value for chi square (goodness of fit)
Appendix: p value
In order to evaluate the goodness of fit of our model to the data, we used the so-called the posterior predictive p value.
In the following, we use general conventional notations. Let \(y_{obs} \) be an observed dataset and \(f(y|\theta)\) be a model (likelihood) for future dataset \(y\). We denote a prior and a posterior distribution by \(\pi(\theta)\) and \(\pi(\theta|y) \propto f(y|\theta)\pi(\theta)\), respectively.
In our case, the data \(y\) is a pair of hits and false alarms; that is, \(y=(H_1,H_2, \dots H_C; F_1,F_2, \dots F_C)\) and \(\theta = (z_1,dz_1,dz_2,\dots,dz_{C-1},\mu, \sigma) \). We define the \(\chi^2\) discrepancy (goodness of fit statistics) to validate that our model fit the data. $$ T(y,\theta) := \sum_{c=1,.......,C} \biggr( \frac{\bigl(H_c-N_L\times p_c(\theta) \bigr)^2}{N_L\times p_c(\theta)}+ \frac{\bigl(F_c- q_{c}(\theta) \times N_{X}\bigr)^2}{ q_{c}(\theta) \times N_{X} }\biggr). $$
for a single reader and a single modality.
$$ T(y,\theta) := \sum_{r=1}^R \sum_{m=1}^M \sum_{c=1}^C \biggr( \frac{(H_{c,m,r}-N_L\times p_{c,m,r}(\theta))^2}{N_L\times p_{c,m,r}(\theta)}+ \frac{\bigl(F_c- q_{c}(\theta) \times N_{X}\bigr)^2}{ q_{c}(\theta) \times N_{X} }\biggr).$$
for multiple readers and multiple modalities.
Note that \(p_c\) and \(\lambda _{c}\) depend on \(\theta\).
In classical frequentist methods, the parameter \(\theta\) is a fixed estimate, e.g., the maximal likelihood estimator. However, in a Bayesian context, the parameter is not deterministic. In the following, we show the p value in the Bayesian sense.
Let \(y_{obs}\) be an observed dataset (in an FROC context, it is hits and false alarms). Then, the so-called posterior predictive p value is defined by
$$ p_value = \int \int \, dy\, d\theta\, I( T(y,\theta) > T(y_{obs},\theta) )f(y|\theta)\pi(\theta|y_{obs}) $$
In order to calculate the above integral, let \(\theta_1,\theta _2, ......., \theta_i,.......,\theta_I\) be samples from the posterior distribution of \( y_{obs} \), namely,
$$ \theta_1 \sim \pi(....|y_{obs} ),$$ $$ .......,$$ $$ \theta_i \sim \pi(....|y_{obs} ),$$ $$ .......,$$ $$ \theta_I \sim \pi(....|y_{obs} ).$$
we obtain a sequence of models (likelihoods), i.e., \(f(....|\theta_1),f(....|\theta_2),......., f(....|\theta_n)\). We then draw the samples \(y^1_1,....,y^i_j,.......,y^I_J \), such that each \(y^i_j\) is a sample from the distribution whose density function is \(f(....|\theta_i)\), namely,
$$ y^1_1,.......,y^1_j,.......,y^1_J \sim f(....|\theta_1),$$ $$.......,$$ $$ y^i_1,.......,y^i_j,.......,y^i_J \sim f(....|\theta_i),$$ $$.......,$$ $$ y^I_1,.......,y^I_j,.......,y^I_J \sim f(....|\theta_I).$$
Using the Monte Carlo integral twice, we calculate the integral of any function \( \phi(y,\theta)\).
$$ \int \int \, dy\, d\theta\, \phi(y,\theta)f(y|\theta)\pi(\theta|y_{obs}) $$ $$\approx \int \, \frac{1}{I}\sum_{i=1}^I \phi(y,\theta_i)f(y|\theta_i)\,dy$$ $$ \frac{1}{IJ}\sum_{i=1}^I \sum_{j=1}^J \phi(y^i_j,\theta_i)$$
In particular, substituting \(\phi(y,\theta):= I( T(y,\theta) > T(y_{obs},\theta) ) \) into the above equation, we can approximate the posterior predictive p value.
$$ p_value \approx \frac{1}{IJ}\sum_i \sum_j I( T(y^i_j,\theta_i) > T(y_{obs},\theta_i) ) $$
ppp_srsc(
StanS4class,
Colour = TRUE,
dark_theme = TRUE,
plot = TRUE,
summary = FALSE,
plot_data = TRUE,
replicate.number.from.model.for.each.MCMC.sample = 100
)An S4 object of class stanfitExtended which is an inherited class from the S4 class stanfit.
This R object is a fitted model object
as a return value of the function fit_Bayesian_FROC().
To be passed to DrawCurves(), ppp() and ... etc
Logical: TRUE of FALSE. whether Colour of curves is dark theme or not.
TRUE or FALSE
Logical, whether replicated data are drawn, in the following notation, replicated data are denoted by \(y_1,y_2,...,y_N\).
Logical: TRUE of FALSE. Whether to print the verbose summary. If TRUE then verbose summary is printed in the R console. If FALSE, the output is minimal. I regret, this variable name should be verbose.
A logical, whether data is plotted in the plot of data synthesized from the posterior predictive distribution I cannot understand what I wrote in the past. My head is crazy cuz I was MCS, head inflammation maybe let me down.
Suppose that $$\theta_1, \theta_2, \theta_3,...,\theta_N$$ are samples drawn in \(N\) times from posterior \(\pi(\theta|D)\) of given data \(D\).
So, these \(\theta_i;i=1,2,...\) are contained in a stanfit object specified as the variable StanS4class.
Let \(y_1,y_2,...,y_n\) be samples drawn as the manner
$$y_1 \sim likelihood ( . |\theta_1), $$ $$y_2 \sim likelihood ( . |\theta_2),$$ $$y_3 \sim likelihood ( .|\theta_3),$$ $$...,$$ $$y_N \sim likelihood ( . |\theta_N).$$
We repeat this in \(J\) times, namely, we draw the samples \(y_{n,j},n=1,..,N;j=1,...,J\) so that
$$y_{1,j} \sim likelihood ( . |\theta_1), $$ $$y_{2,j} \sim likelihood ( . |\theta_2),$$ $$y_{3,j} \sim likelihood ( .|\theta_3),$$ $$...,$$ $$y_{n,j} \sim likelihood ( . |\theta_n),$$ $$...,$$ $$y_{N,j} \sim likelihood ( . |\theta_N).$$
Yes, the variable replicate.number.from.model.for.each.MCMC.sample means \(J\)!
We can write it more explicitly without abbreviation as follows.
$$y_{1,1},y_{1,2},...,y_{1,j},...,y_{1,J} \sim likelihood ( . |\theta_1), $$ $$y_{2,1},y_{2,2},...,y_{2,j},...,y_{2,J} \sim likelihood ( . |\theta_2), $$ $$y_{3,1},y_{3,2},...,y_{3,j},...,y_{3,J} \sim likelihood ( . |\theta_3), $$ $$...,$$ $$y_{n,1},y_{n,2},...,y_{n,j},...,y_{n,J} \sim likelihood ( . |\theta_n), $$ $$...,$$ $$y_{N,1},y_{N,2},...,y_{N,j},...,y_{N,J} \sim likelihood ( . |\theta_N). $$
Now, my body is not so good, so, I am tired. Cuz I counld not understand what I wrote, so I reviesed in 2020 Aug 9.
You health is very bad condition, so, if the sentence is not clear, it is also for me! even if I wrote it! So, If I notice that my past brain is broken, then I will revise. Ha,,, I want be rest in peace.
A positive integer, representing \(J\) in the following notation.
A list, including p value and materials to calculate it.
Contents of the list as a return values is the following:
FPF,TPF,..etc data \(y_{n,j} \sim likelihood ( . |\theta_n),\)
chisq_at_observed_data $$\chi (D|\theta_1), \chi (D|\theta_2), \chi (D|\theta_3),...,\chi (D|\theta_n),$$
chisq_not_at_observed_data $$\chi (y_1|\theta_1), \chi (y_2|\theta_2), \chi (y_3|\theta_3),...,\chi (y_n|\theta_n), $$
Logical The i-th component is a logical vector indicating whether $$\chi (y_2|\theta_2) > \chi (D|\theta_2)$$ is satisfied or not. Oppai ga Ippai. If TRUE, then the inequality holds.
p.value From the component Logical, we calculate the so-called Posterior Predictive P value. Note that the author hate this notion!! I hate it!! Akkan Beeeee!!!
In addition, this function plots replicated datasets from model at each MCMC sample generated by HMC. Using the Hamiltonian Monte Carlo Sampling: HMC. we can draw the MCMC samples of size \(n\), say $$\theta_1, \theta_2, \theta_3,...,\theta_n $$, namely, $$\theta_1 \sim \pi(.|D), $$ $$\theta_2 \sim \pi(.|D), $$ $$\theta_3 \sim \pi(.|D),$$ $$...,$$ $$\theta_n \sim \pi(.|D).$$ where \(\pi(\theta|D)\) is the posterior for given data \(D\).
We draw samples as follows.
$$y_{1,1},y_{1,2},...,y_{1,j},...,y_{1,J} \sim likelihood ( . |\theta_1), $$ $$y_{2,1},y_{2,2},...,y_{2,j},...,y_{2,J} \sim likelihood ( . |\theta_2), $$ $$y_{3,1},y_{3,2},...,y_{3,j},...,y_{3,J} \sim likelihood ( . |\theta_3), $$ $$...,$$ $$y_{n,1},y_{n,2},...,y_{n,j},...,y_{n,J} \sim likelihood ( . |\theta_n), $$ $$...,$$ $$y_{N,1},y_{N,2},...,y_{N,j},...,y_{N,J} \sim likelihood ( . |\theta_N).$$
Then we calculates the chi-squares for each sample.
$$ \chi(y_{1,1}|\theta_1), \chi(y_{1,2}|\theta_1), \chi(y_{1,3}|\theta_1),..., \chi(y_{1,j}|\theta_1),...., \chi(y_{1,J}|\theta_1),$$ $$ \chi(y_{2,1}|\theta_2), \chi(y_{2,2}|\theta_2), \chi(y_{2,3}|\theta_2),..., \chi(y_{2,j}|\theta_2),...., \chi(y_{2,J}|\theta_2),$$ $$ \chi(y_{3,1}|\theta_3), \chi(y_{3,2}|\theta_3), \chi(y_{3,3}|\theta_3),..., \chi(y_{3,j}|\theta_3),...., \chi(y_{3,J}|\theta_3),$$ $$...,$$ $$ \chi(y_{i,1}|\theta_i), \chi(y_{i,2}|\theta_i), \chi(y_{i,3}|\theta_i),..., \chi(y_{i,j}|\theta_i),...., \chi(y_{I,J}|\theta_i),$$ $$...,$$ $$ \chi(y_{I,1}|\theta_I), \chi(y_{I,2}|\theta_I), \chi(y_{I,3}|\theta_I),..., \chi(y_{I,j}|\theta_I),...., \chi(y_{I,J}|\theta_I).$$
where \(L ( . |\theta_i)\) is a likelihood at parameter \(\theta_i\).
Let \( \chi(y|\theta) \) be a chi square goodness of fit statistics of our hierarchical Bayesian Model
$$\chi(y|\theta) := \sum_{r=1}^R \sum_{m=1}^M \sum_{c=1}^C ( \frac { ( H_{c,m,r}-N_L\times p_{c,m,r})^2}{N_L\times p_{c,m,r}} + \frac{(F_{c,m,r}-(\lambda_{c} -\lambda_{c+1} )\times N_{L})^2}{(\lambda_{c} -\lambda_{c+1} )\times N_{L} }).$$
and a chi square goodness of fit statistics of our non-hierarchical Bayesian Model
$$\chi(y|\theta) := \sum_{c=1}^C \biggr( \frac{( H_{c}-N_L\times p_{c})^2}{N_L\times p_{c}} + \frac{(F_{c}-(\lambda_{c} -\lambda_{c+1} ) )\times N_{L}]^2}{(\lambda_{c} -\lambda_{c+1} )\times N_{L} }\biggr).$$
where a dataset \(y\) denotes \( (F_{c,m,r}, H_{c,m,r}) \) in MRMC case and \( (F_{c}, H_{c}) \) in a single reader and a single modality case, and model parameter \(\theta\).
Then we can calculate the posterior predictive p value for a given dataset \(y_0\).
$$ \int \int I( \chi(y|\theta) > \chi(y_0|\theta) ) f(y|\theta) \pi(\theta|y_0) d \theta d y $$ $$ \approx \int \sum_i I( \chi(y|\theta_i) > \chi(y_0|\theta_i) ) f(y|\theta_i) d y $$ $$ \approx \sum_{j=1}^J \sum_{i=1}^I I( \chi(y_{i,j}|\theta_i) > \chi(y_0|\theta_i) ) $$
When we plot these synthesized data-sets \(y_{i,j}\), we use the jitter() which adds a small amount of noise to avoid overlapping points.
For example, jitter(c(1,1,1,1)) returns values: 1.0161940 1.0175678 0.9862400 0.9986126, which is changed
from 1,1,1,1 to be not exactly 1 by adding tiny errors to avoid overlapping. I love you. 2019 August 19
Nowadays, I cannot remove my self from some notion, such as honesty, or pain, or,.. maybe these thing is no longer with myself.
This programm is made to fix previous release calculation. Now, this programm calculates correct p value.
So... I calculate the ppp for MCMC and Graphical User Interface based on Shiny for MRMC, which should be variable such as number of readers, modalities, to generate such ID vectors automatically. Ha,... tired! Boaring, I want to die...t, diet!! Tinko, tinko unko unko. Manko manko. ha.
Leberiya, he will be die, ha... he cannot overcome, very old, old guy. I will get back to meet him. Or I cannot meet him? Liberiya,...very wisdom guy, Ary you already die? I will get back with presents for you. Ball, I have to throgh ball, and he will catch it.
The reason why the author made the plot of data drawn from Posterior Predictive likelihoods with each MCMC parameters is to understand our programm is correct, that is, each drawing is very mixed. Ha,.... when wright this,... I always think who read it. I love you, Ruikobach. Ruikobach is tiny and tiny, but,... cute. Ruikosan...Ruiko... But he has time only several years. He will die, he lives sufficiently so long, ha.
Using this function, user would get reliable posterior predictive p values, Cheers! Pretty Crowd!
We note that the calculation of posterior perdictive p value (PPP) relies on the law of large number. Thus, in order to obtain the relicable PPP, we need to enough large MCMC samples to approximate the double integral of PPP. For example, the MCMC samples is small, then R hat is far from 1 but, the low MCMC samples leads us to incorrect p value which sometimes said that the model is correct even if the R hat criteria reject the MCMC results.
# NOT RUN {
# }
# NOT RUN {
#========================================================================================
# 1) Create a fitted model object with data named "d"
#========================================================================================
fit <- fit_Bayesian_FROC( dataList = d,
ite = 222 # to restrict running time, but it is too small
)
#========================================================================================
# 2) Calculate p value and meta data
#========================================================================================
ppp <- ppp_srsc(fit)
#========================================================================================
# 3) Extract a p value
#========================================================================================
ppp$p.value
# Revised 2019 August 19
# Revised 2019 Nov 27
Close_all_graphic_devices() # 2020 August
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab