# 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]); 

}

}