BayesianFROC (version 0.5.0)

ppp_srsc: Calculates PPP for Models of a single reader and a single modality (Calculation is correct! :'-D)

Description

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,d{z}_1,d{z}_2,\dots ,d{z}_{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}(\frac{(H_c-N_L\times p_c(\theta ){)}^2}{N_L\times p_c(\theta )}+\frac{(F_c-q_c(\theta )\times N_X{)}^2}{q_c(\theta )\times N_X}).$$

for a single reader and a single modality.

$$T(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}(\theta ){)}^2}{N_L\times p_{c,m,r}(\theta )}+\frac{(F_c-q_c(\theta )\times N_X{)}^2}{q_c(\theta )\times N_X}).$$

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_{v}alue=\int \int dyd\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_j^i,.......,y_J^I$, such that each $y_j^i$ is a sample from the distribution whose density function is $f(....|{\theta }_i)$, namely,

$$y_1^1,.......,y_j^1,.......,y_J^1\sim f(....|{\theta }_1),$$ $$.......,$$ $$y_1^i,.......,y_j^i,.......,y_J^i\sim f(....|{\theta }_i),$$ $$.......,$$ $$y_1^I,.......,y_j^I,.......,y_J^I\sim f(....|{\theta }_I).$$

Using the Monte Carlo integral twice, we calculate the integral of any function $\varphi (y,\theta )$.

$$\int \int dyd\theta \varphi (y,\theta )f(y|\theta )\pi (\theta |y_{obs})$$ $$\approx \int \frac{1}{I}\sum _{i=1}^I\varphi (y,{\theta }_i)f(y|{\theta }_i)dy$$ $$\frac{1}{IJ}\sum _{i=1}^I\sum _{j=1}^J\varphi (y_j^i,{\theta }_i)$$

In particular, substituting $\varphi (y,\theta ):=I(T(y,\theta )>T(y_{obs},\theta ))$ into the above equation, we can approximate the posterior predictive p value.

$$p_{v}alue\approx \frac{1}{IJ}\sum _i\sum _{j}I(T(y_j^i,{\theta }_i)>T(y_{obs},{\theta }_i))$$

Usage

ppp_srsc(
  StanS4class,
  Colour = TRUE,
  dark_theme = TRUE,
  plot = TRUE,
  summary = FALSE,
  plot_data = TRUE,
  replicate.number.from.model.for.each.MCMC.sample = 100
)

Arguments

StanS4class

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

Colour

Logical: TRUE of FALSE. whether Colour of curves is dark theme or not.

dark_theme

TRUE or FALSE

plot

Logical, whether replicated data are drawn, in the following notation, replicated data are denoted by $y_1,y_2,...,y_N$.

summary

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.

plot_data

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.

replicate.number.from.model.for.each.MCMC.sample

A positive integer, representing $J$ in the following notation.

Value

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!!!

Details

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(\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}).$$

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 dy$$ $$\approx \int \sum _{i}I(\chi (y|{\theta }_i)>\chi (y_0|{\theta }_i))f(y|{\theta }_i)dy$$ $$\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.

Examples

# 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 {


# }