Use Bayesian hierarchical model to predict new data points I have a data set $(n_i,y_i),i=0,...,10$. I modeled it as a Bayesian hierarchical beta-binomial model.
$y_i∼Binomial(n_i,p_i)$ and $p_i∼Beta(\alpha,\beta)$.
I have used MCMC to estimate $\alpha$ and $\beta$ (use the median as the estimated value).
Given a new $(n_j, y_j)$, I want to judge if it comes from the same model as the training set $(i=0,...,10)$, namely $p_j$ is significantly greater than $\theta$ (shrinkage $p_i$ from all group).
How can I use the estimated $\alpha$ and $\beta$ to compute it?
 A: In your notation, if you have a sample of, say, 1000 points from the posterior distribution of $\alpha$ and $\beta$ (which you probably have, since you've MCMCed the model) stored in vectors alpha and beta, you can get a sample of 1000 points from the posterior predictive distribution $y_{11}\mid y_1,\dots,y_{10}$ doing rbinom(1000, n_11, rbeta(1000, alpha, beta)). With this sample, you can compute any desired summaries, such as the posterior predictive mean, build histograms, etc. Take a careful look at these notes (especially page 7).
A: In Bayesian statistics, such predictions as usually done using the posterior predictive distribution, not by estimating the parameters directly. See the Wikipedia page on the posterior predictive distribution.
If the prior distribution is beta, and the likelihood is binomial, then the posterior predictive distribution is a beta-binomial distribution. See
An introduction to the Beta-Binomial model, which has details and formulas. The previous sentence beginning "If the prior distribution" is a direct quote from pg 8 of this document. Also see the Wikipedia page on the Beta-binomial distribution.
You probably want the value of the posterior predictive distribution for your new data point to start with.
A: I would start by doing a simulation to get a sample of $y$ values when $n=n_{new}$ from the training distribution. Then subtract of $y_{new}$ from the $y$ sample you got, and check if it's credibly non-zero. Here's some code to do that:
modelDescription = "
model {
    for( i in 1 : ny ) {
    y[i] ~ dbin(p[i],n[i])
    p[i] ~ dbeta(alpha,beta)
    }
    alpha ~ dunif(0,1e3)
    beta ~ dunif(0,1e3) #some priors...
}
"


sampleSize=500
y=rbinom(sampleSize,10,0.5)
n=round(y*10)-rbinom(sampleSize,2,0.5)+rbinom(sampleSize,2,0.5) #set p=1/10 with someRandomFluctuactions
n<-ifelse(n<1,1,n) #ensure all n>0

dataList = list(
    ny = length(y) ,
    y = y ,
    n = n
)

parameters = c("alpha","beta")

#this is just a wrapper I have for rjags::coda.samples and dclone::codaSamples
test1 = bjmCreate("model1",parameters,modelDescription,dataList,60000)

require(coda)
test1p <- window(test1,start(test1)+1000) #Burn-in
chain<-as.matrix(test1p)
alpha<-chain[,1]
beta<-chain[,2]
chainSize<-nrow(chain) #note that nrow(chain)==length(alpha)==length(beta)

In the code above, we got an alpha and beta sample from MCMC, now we are going to generate a sample for $y$ when $n=50$
We model $y$ as $y \sim Binomial(n,p)$. So to generate sample for $y$, we need to do y<-rbinom(chainSize,newN,p).
We know the value of $n_{new}$, but we don't have $p$, but we can get a sample of it by remebering $p \sim Beta(\alpha,\beta)$, so we get a $p$ sample from p<-rbeta(chainSize,alpha,beta).
So, by pluging one thing into the other:
newN = 50 #setting here the value of N for which we'll predict y
p<-rbeta(chainSize,alpha,beta)
y<-rbinom(chainSize,newN,p)

Now we have a sample for $y$ when $n=50$, and our work would be done if we just wanted to predict it. But instead, you want to check if it an specific value of $y$ is from that distribution, so we do
newY = 5
yDiff <- y-newY

And check if the difference between both is credibly zero.
#Analisys of HPD
print(HPDinterval(mcmc(yDiff)))
densityplot(mcmc(yDiff))

I got the following result from densityplot(mcmc(yDiff)):
     lower upper
var1    -4     4
attr(,"Probability")
[1] 0.95

That is, zero is within the 95% HDP, and thus, we can not reject the hypothesis that the proposed $(y_{new},n_{new})$ "comes from the same model as the training set".
This, in fact, makes sense since the our "true" value p=1/10 is exactly the values we used ($y_{new}=5$ and $n_{new}=50$).
