##################################################################################### # Example analysing a data set using Bayesian two-phase approach # Case-control sampling at Phase I # Data should be laid out as in Table 1 of Ross and Wakefield (2012) # Function "Bayes2phaseCC.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 case-control example from paper using sex as confounder (collapse across age groups) # Phase II data set includes smoking status (x) library(mvtnorm) source("BayesTwoPhase_CaseControl.R") Nvec<-c(6310, 281, 5576, 528) nvec<-matrix(c(1053,34,1319,30, 91,19,645,69, 77,56,137,50, 88,125,350,280), nrow=4, ncol=4) results<-Bayes2phaseCC.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(2,nvec[3,1]), rep(3,nvec[4,1]), rep(0, nvec[1,2]), rep(1, nvec[2,2]), rep(2, nvec[3,2]), rep(3, nvec[4,2]), rep(0, nvec[1,3]), rep(1,nvec[2,3]), rep(2,nvec[3,3]), rep(3,nvec[4,3]), rep(0, nvec[1,4]), rep(1, nvec[2,4]), rep(2, nvec[3,4]), rep(3, nvec[4,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 ##############################################################