# migration.bug # This file contains WinBUGS model statements for the analysis of the # Thai migration data analyzed in Western and Garip (2005), "Bayesian # Analysis of Comparative Survey Data." # 5/4/05 #------------------------------------------------------------------------ # The Pooled Model #------------------------------------------------------------------------ model { for(i in 1:N) { logit(p[i]) <- alpha0 + alpha[1]*male[i] + alpha[2]*age[i] + alpha[3]*educ[i] + alpha[4]*itrip[i] + gamma[1]*vtrip[i] + gamma[2]*vgini[i]; y[i] ~ dbern(p[i]); yrep[i] ~ dbern(p[i]); } alpha0 ~ dnorm(0, 1.0E-6); alpha[1] ~ dnorm(0, 1.0E-6); alpha[2] ~ dnorm(0, 1.0E-6); alpha[3] ~ dnorm(0, 1.0E-6); alpha[4] ~ dnorm(0, 1.0E-6); gamma[1] ~ dnorm(0, 1.0E-6); gamma[2] ~ dnorm(0, 1.0E-6); # Posterior predictive disribution by village for(j in 1:J) { ppred[j] <- mean(yrep[ind[j]:(ind[j+1]-1)]); pobs[j] <- mean(y[ind[j]:(ind[j+1]-1)]); pval[j] <- step(ppred[j]-pobs[j]); } # Posterior predictive disribution by village and sex for (i in 1:N) { # males yrepm[i] <- yrep[i]*male[i]; ym[i] <- y[i]*male[i]; # females yrepf[i] <- yrep[i]*(1-male[i]); yf[i] <- y[i]*(1-male[i]); } for (j in 1:J) { # compute village sizes by sex sizem[j] <- sum(male[ind[j]:(ind[j+1]-1)]); sizef[j] <- (ind[j+1]-ind[j]) - sizem[j]; # males ppredm[j] <- sum(yrepm[ind[j]:(ind[j+1]-1)])/sizem[j] ; pobsm[j] <- sum(ym[ind[j]:(ind[j+1]-1)])/sizem[j]; pvalm[j] <- step(ppredm[j]-pobsm[j]); #females ppredf[j] <- sum(yrepf[ind[j]:(ind[j+1]-1)])/sizef[j]; pobsf[j] <- sum(yf[ind[j]:(ind[j+1]-1)])/sizef[j]; pvalf[j] <- step(ppredf[j]-pobsf[j]); } } #------------------------------------------------------------------------ # The Fixed Effects Model #------------------------------------------------------------------------ model { for(i in 1:N) { logit(p[i]) <- b0[vid[i]] + alpha[1]*male[i] + alpha[2]*age[i] + alpha[3]*educ[i] + alpha[4]*itrip[i]; y[i] ~ dbern(p[i]); yrep[i] ~ dbern(p[i]); } for (j in 1:J) { b0[j] ~ dnorm(0, 1.0E-6); } alpha[1] ~ dnorm(0, 1.0E-6); alpha[2] ~ dnorm(0, 1.0E-6); alpha[3] ~ dnorm(0, 1.0E-6); alpha[4] ~ dnorm(0, 1.0E-6); # Compute the intercept - mean of village level effects bbar <- mean(b0[1:22]); # Posterior predictive disribution for(j in 1:J) { ppred[j] <- mean(yrep[ind[j]:(ind[j+1]-1)]); pobs[j] <- mean(y[ind[j]:(ind[j+1]-1)]); pval[j] <- step(ppred[j]-pobs[j]); } # Posterior predictive disribution by village and sex for (i in 1:N) { # males yrepm[i] <- yrep[i]*male[i]; ym[i] <- y[i]*male[i]; # females yrepf[i] <- yrep[i]*(1-male[i]); yf[i] <- y[i]*(1-male[i]); } for(j in 1:J) { # compute village sizes by sex sizem[j] <- sum(male[ind[j]:(ind[j+1]-1)]); sizef[j] <- (ind[j+1]-ind[j]) - sizem[j]; # males ppredm[j] <- sum(yrepm[ind[j]:(ind[j+1]-1)])/sizem[j] ; pobsm[j] <- sum(ym[ind[j]:(ind[j+1]-1)])/sizem[j]; pvalm[j] <- step(ppredm[j]-pobsm[j]); #females ppredf[j] <- sum(yrepf[ind[j]:(ind[j+1]-1)])/sizef[j]; pobsf[j] <- sum(yf[ind[j]:(ind[j+1]-1)])/sizef[j]; pvalf[j] <- step(ppredf[j]-pobsf[j]); } } #------------------------------------------------------------------------ # The Random Effects Model #------------------------------------------------------------------------ model { for(i in 1:N) { logit(p[i]) <- b0[vid[i]] + alpha0 + alpha[1]*male[i] + alpha[2]*age[i] + alpha[3]*educ[i] + alpha[4]*itrip[i] + gamma[1]*vtrip[i] + gamma[2]*vgini[i]; y[i] ~ dbern(p[i]); yrep[i] ~ dbern(p[i]); } for(j in 1:J) { b0[j] ~ dnorm(0.0, tau0); } # Gamma priors on precision; tau0 ~ dgamma(1.0E-3, 1.0E-3); sigma0 <- 1/sqrt(tau0); alpha0 ~ dnorm(0, 1.0E-6); alpha[1] ~ dnorm(0, 1.0E-6); alpha[2] ~ dnorm(0, 1.0E-6); alpha[3] ~ dnorm(0, 1.0E-6); alpha[4] ~ dnorm(0, 1.0E-6); gamma[1] ~ dnorm(0, 1.0E-6); gamma[2] ~ dnorm(0, 1.0E-6); # Posterior predictive disribution for(j in 1:J) { ppred[j] <- mean(yrep[ind[j]:(ind[j+1]-1)]); pobs[j] <- mean(y[ind[j]:(ind[j+1]-1)]); pval[j] <- step(ppred[j]-pobs[j]); } # Posterior predictive disribution by village and sex for (i in 1:N) { # males yrepm[i] <- yrep[i]*male[i]; ym[i] <- y[i]*male[i]; # females yrepf[i] <- yrep[i]*(1-male[i]); yf[i] <- y[i]*(1-male[i]); } for(j in 1:J) { # compute village sizes by sex sizem[j] <- sum(male[ind[j]:(ind[j+1]-1)]); sizef[j] <- (ind[j+1]-ind[j]) - sizem[j]; # males ppredm[j] <- sum(yrepm[ind[j]:(ind[j+1]-1)])/sizem[j] ; pobsm[j] <- sum(ym[ind[j]:(ind[j+1]-1)])/sizem[j]; pvalm[j] <- step(ppredm[j]-pobsm[j]); #females ppredf[j] <- sum(yrepf[ind[j]:(ind[j+1]-1)])/sizef[j]; pobsf[j] <- sum(yf[ind[j]:(ind[j+1]-1)])/sizef[j]; pvalf[j] <- step(ppredf[j]-pobsf[j]); } } #------------------------------------------------------------------------ # The Random Slope and Intercept Model #------------------------------------------------------------------------ model { for(i in 1:N) { logit(p[i]) <- b0[vid[i]] + alpha0 + alpha1[vid[i]]*male[i] + alpha2[vid[i]]*age[i] + alpha3[vid[i]]*educ[i] + alpha4[vid[i]]*itrip[i] + gamma[1]*vtrip[i] + gamma[2]*vgini[i]; y[i] ~ dbern(p[i]); yrep[i] ~ dbern(p[i]); } for(j in 1:J) { b0[j] ~ dnorm(0.0, tau0); alpha1[j] ~ dnorm(mu1, tau1); alpha2[j] ~ dnorm(mu2, tau2); alpha3[j] ~ dnorm(mu3, tau3); alpha4[j] ~ dnorm(mu4, tau4); } alpha0 ~ dnorm(0, 1.0E-6); gamma[1] ~ dnorm(0, 1.0E-6); gamma[2] ~ dnorm(0, 1.0E-6); mu1 ~ dnorm(0, 1.0E-6); mu2 ~ dnorm(0, 1.0E-6); mu3 ~ dnorm(0, 1.0E-6); mu4 ~ dnorm(0, 1.0E-6); # Gamma priors on precision; tau0 ~ dgamma(1.0E-3, 1.0E-3); tau1 ~ dgamma(1.0E-3, 1.0E-3); tau2 ~ dgamma(1.0E-3, 1.0E-3); tau3 ~ dgamma(1.0E-3, 1.0E-3); tau4 ~ dgamma(1.0E-3, 1.0E-3); sigma0 <- 1/sqrt(tau0); sigma1 <- 1/sqrt(tau1); sigma2 <- 1/sqrt(tau2); sigma3 <- 1/sqrt(tau3); sigma4 <- 1/sqrt(tau4); # Posterior predictive disribution for(j in 1:J) { ppred[j] <- mean(yrep[ind[j]:(ind[j+1]-1)]); pobs[j] <- mean(y[ind[j]:(ind[j+1]-1)]); pval[j] <- step(ppred[j]-pobs[j]); } # Posterior predictive disribution by village and sex for (i in 1:N) { # males yrepm[i] <- yrep[i]*male[i]; ym[i] <- y[i]*male[i]; # females yrepf[i] <- yrep[i]*(1-male[i]); yf[i] <- y[i]*(1-male[i]); } for(j in 1:J) { # compute village sizes by sex sizem[j] <- sum(male[ind[j]:(ind[j+1]-1)]); sizef[j] <- (ind[j+1]-ind[j]) - sizem[j]; # males ppredm[j] <- sum(yrepm[ind[j]:(ind[j+1]-1)])/sizem[j] ; pobsm[j] <- sum(ym[ind[j]:(ind[j+1]-1)])/sizem[j]; pvalm[j] <- step(ppredm[j]-pobsm[j]); #females ppredf[j] <- sum(yrepf[ind[j]:(ind[j+1]-1)])/sizef[j]; pobsf[j] <- sum(yf[ind[j]:(ind[j+1]-1)])/sizef[j]; pvalf[j] <- step(ppredf[j]-pobsf[j]); } }