rm(list = ls())
# While only doing the computation for logistic
# models, we include some code for probic and
# complementary log-log
logit <- function(p)
{
lg = log(p/(1-p))
return(lg)}

logistic <- function(w)
{
lg = exp(w)/(1+exp(w))
return(lg)}

Gum <- function(w)
{
lg = 1-(exp(-exp(w)))
return(lg)}

cll <- function(p)
{
lg = log(-log(1-p))
return(lg)}

### Fletch used BIDA not LOGLIN3 Priors
#   so I have changed them
#a <- c( 1.6, 1 )
#b <- c( 1, 1.6 )
# LOGLIN3 priors
a <- c( 1, .577 )
b <- c( .577, 1 )
iterates=2000    # increase

ptilde1 <- rbeta(iterates, a[1], b[1])
ptilde2 <- rbeta(iterates, a[2], b[2])
# either
# ptildea=t(cbind(ptilde1,ptilde2))
# or
ptilde=c(ptilde1,ptilde2)
ptilde=t(matrix(ptilde,ncol=2))

Xtilde=matrix(c(1,1,55,75),ncol=2)
Xtildeinv=solve(Xtilde)

# enter data
n=23
y=c(1,1,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0,0,0,0,0)
N=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
tau=c(53,57,58,63,66,67,67,67,68,69,70,70,70,70,72,
73,75,75,76,76,78,79,81)
X=matrix(c(rep(1,23),tau),ncol=2)
X

pprobit=pnorm(X %*% Xtildeinv %*% qnorm(ptilde))

plogit=logistic(X %*% Xtildeinv %*% logit(ptilde))

pcll=Gum(X %*% Xtildeinv %*% cll(ptilde))


# define likelihood function
den <- function(p,y,N)
{
fp = prod(c(p^y, (1-p)^(N-y)))
# For nonbinary data
#fp = prod(c(p^y, (1-p)^(N-y),choose(N,y)))
return(fp)}
#choose(N, y) or factorial(N)/(factorial(y)*factorial(N-y))

# Marginal density function
# In Link Selection, computed but not a function
marg <- function(y,N)
{

Mprobitd=rep(1,iterates)
Mlogitd=rep(1,iterates)
Mclld=rep(1,iterates)

for(r in 1:iterates){
Mprobitd[r]=den(pprobit[,r],y,N)
Mlogitd[r]=den(plogit[,r],y,N)
Mclld[r]=den(pcll[,r],y,N)}

Mprobit=mean(Mprobitd)
Mlogit=mean(Mlogitd)
Mcll=mean(Mclld)
# change "return" statement for testing
# alternative links
return(Mlogit)
}
# Everything up to this point has been similar to
# the link selection program

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