#####################################################################################
# Example analysing a data set using Bayesian two-phase approach
# Simple random sampling at phase I
# Data should be laid out as in Table 1 of Ross and Wakefield (2012)
# Function "Bayes2phase.func" takes input arguments:
#   - Nvec: vector of length 2m_z of phase I sample sizes
#   - nvec: 2 x m_x x m_z matrix containing phase II data
#   - maxiter: maximum number of iterations MCMC should perform (not including iterations for tuning)
# Function returns a list containing beta parameter estimates and 95% credible intervals
# Requires "tps.q"
#####################################################################################

# Use population data from simulation study so we can compare results
# Randomly sampled from this population to obtain a phase II data set

library(mvtnorm)

source("BayesTwoPhase.R")

Ntruth<-matrix(c(3930,430,94,10,109,47,212,168), nrow=2, ncol=4, byrow=TRUE)
Nvec<-apply(Ntruth, 2, sum)

Y<-c(rep(0, Nvec[1]), rep(1, Nvec[2]), rep(0, Nvec[3]), rep(1, Nvec[4]))
X<-c(rep(0, Ntruth[1,1]), rep(1,Ntruth[2,1]), rep(0, Ntruth[1,2]), rep(1, Ntruth[2,2]), rep(0, Ntruth[1,3]), rep(1,Ntruth[2,3]), rep(0, Ntruth[1,4]), rep(1, Ntruth[2,4]))
Z<-c(rep(0, Nvec[1]+Nvec[2]), rep(1, Nvec[3]+Nvec[4]))

# true model
options(contrasts=c("contr.treatment", "contr.poly"))
fit.true<-glm(Y~X+Z, family=binomial)
fit.true.CI<-confint.default(fit.true)

nvec<-matrix(c(238,92,72,5,12,8,178,95), nrow=2, ncol=4)

results<-Bayes2phase.func(Nvec, nvec, 250000)
results

##############################################################
## To compare results to NPML method
ntot<-apply(nvec, 2, sum)

y<-c(rep(0, ntot[1]), rep(1, ntot[2]), rep(0, ntot[3]), rep(1, ntot[4]))
x<-c(rep(0, nvec[1,1]), rep(1,nvec[2,1]), rep(0, nvec[1,2]), rep(1, nvec[2,2]), rep(0, nvec[1,3]), rep(1,nvec[2,3]), rep(0, nvec[1,4]), rep(1, nvec[2,4]))
z<-c(rep(0, ntot[1]+ntot[2]), rep(1, ntot[3]+ntot[4]))

# fit NPML method to get starting values and proposal variance-covariance matrix
options(contrasts=c("contr.treatment", "contr.poly"))
gr<-z+1
fitML<-tps(y~factor(x)+factor(z), nn0=Nvec[c(1,3)], nn1=Nvec[c(2,4)], group=gr, method="ML")
fitML$coef
##############################################################