//Bayesian nonparametric finite population inference when the survey outcome is binary //reference: Si, Y., Pillai, N. and Gelman, A. (2015), Bayesian nonparametric weighted sampling inference, Bayesian Analysis, 10(3), 605-625. //Author: Yajuan Si //Latest edit date: 03-25-2017 data { //input data elements real N; //population size, define it as real to be used in n/N int n; //sample size int J; //number of poststratification cells vector[J] vec_n; //vector of sample cell sizes vector[J] vec_w; //vector of cell weight values int vec_y[J]; //cell total of survey bonary response } transformed data { vector[J] x; //covariates in GP regression, logarithm of the weights vector[J] vec2_y; //change vec_y to a vector, convenience for vector operation x <- log(vec_w); for (j in 1:J){ vec2_y[j] <- vec_y[j]; } } } parameters { real beta_raw; //reparameterization of beta: beta=beta_raw * sd_beta (2.5) real eta; //partial sill in GP covariance kernel real kappa; //decay parameter in GP covariance kernel real delta_w; //weight assigned to independent prior specification vector[J] vec_N_p; //population cell size proportiong N_j/N vector[J] z; //linear combination predictor in the logistric regression with the coronical parameter } transformed parameters { real beta; //coefficients in the mean function of GP real eta_sq; //square of eta real c; //normalizing constant for (j in 1:J) beta <- 2.5 * beta_raw; eta_sq <- pow(eta, 2); c <- sum(vec_N_p) / N * n / sum(vec_N_p ./ vec_w); } model { // 1) y|w under GP prior vector[J] mu0; //mean of mu matrix[J,J] Sigma; //covariance matrix mu0 <- x * beta; for (i in 1:J) { for (j in i:J) { Sigma[i,j] <- eta_sq * (1-delta_w) * exp(-kappa * fabs(x[i] - x[j])); //squared exponential kernel Sigma[j,i] <- Sigma[i,j]; } Sigma[i,i] <- Sigma[i,i] + eta_sq * delta_w; //combined with independent prior } z ~ multi_normal(mu0, Sigma); vec_y ~ binomial_logit(vec_n, z); //logistic regression // 2) n_j|w vec_n ~ poisson(c * (vec_N_p ./ vec_w)); //this can be changes as multinomial distribution when J is small n ~ poisson(c * sum(vec_N_p ./ vec_w)); vec_n ~ multinomial(vec_N_p ./ vec_w/sum(vec_N_p ./ vec_w)); // 3) prior for hyperparameters, weakly informative eta ~ cauchy(0,5); kappa ~ gamma(0.5,0.5); beta_raw ~ normal(0,1); } generated quantities {//predictions vector[J] theta; //probability real y_est; y_est <- 0; for (j in 1:J){ theta[j] <- inv_logit(z[j]); } y_est <- y_prd - sum(theta .* vec_n- vec2_y) / N /sum(vec_N_p) ; }