# Full data model and outlier check
# For different links, also change "plogit" term
Yf=rep(1,n)
c=0
cc=rep(0,n)
for(r in 1:iterates){
for(i in 1:n){
Yf[i] = rbinom(1,N[i],plogit[i,r])}
for(k in 1:n){
if (marg(Yf[k],N[k]) <= marg(y[k],N[k]))  cc[k] = cc[k] + 1}
if (marg(Yf,N) <= marg(y,N))  c = c + 1}
Pvalue=c/iterates
Pvalue
PV=cc/iterates
PV