The basic unit-level model (Battese et al. 1988; Rao and Molina 2015) is given by \[ y_j = \beta' x_j + v_{i[j]} + \epsilon_j\,,\\ v_i \stackrel{\mathrm{iid}}{\sim} {\cal N}(0, \sigma_v^2) \qquad \epsilon_j \stackrel{\mathrm{iid}}{\sim} {\cal N}(0, \sigma^2) \] where \(j\) runs from 1 to \(n\), the number of unit-level observations, \(\beta\) is a vector of regression coefficients for given covariates \(x_j\), and \(v_i\) are random area intercepts.
We use the api dataset included in packages survey.
library(survey)
data(api)
apipop$cname <- as.factor(apipop$cname)
apisrs$cname <- factor(apisrs$cname, levels=levels(apipop$cname))The apipop data frame contains the complete population
whereas apisrs is a simple random sample from it. The
variable cname is the county name, and we will be
interested in estimation at the county level. Not all counties in the
population are sampled. In order to be able to make predictions for
out-of-sample areas we make sure that the levels of the sample’s
cname variable match those of its population
counterpart.
The basic unit-level model with county random area effects is fit as follows
library(mcmcsae)
mod <- api00 ~
reg(~ ell + meals + stype + hsg + col.grad + grad.sch, name="beta") +
gen(factor = ~ iid(cname), name="v")
sampler <- create_sampler(mod, data=apisrs)
sim <- MCMCsim(sampler, store.all=TRUE, verbose=FALSE)
summary(sim)## llh_ :
## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## llh_ -1090 4.41 -247 0.126 -1097 -1089 -1083 1220 1
##
## beta :
## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## (Intercept) 801.860 25.335 31.65 0.52070 760.804 802.111 843.6792 2367 1.000
## ell -1.958 0.305 -6.41 0.00620 -2.455 -1.961 -1.4450 2425 1.001
## meals -1.959 0.282 -6.96 0.00586 -2.424 -1.963 -1.4923 2310 0.999
## stypeH -106.481 13.589 -7.84 0.25457 -128.270 -106.618 -83.7890 2849 0.999
## stypeM -60.176 11.444 -5.26 0.21831 -78.938 -60.057 -41.4230 2748 0.999
## hsg -0.615 0.413 -1.49 0.00754 -1.306 -0.609 0.0622 3000 0.999
## col.grad 0.590 0.502 1.18 0.01029 -0.228 0.578 1.4104 2379 1.000
## grad.sch 2.125 0.484 4.39 0.00902 1.349 2.126 2.9364 2878 1.001
##
## v_sigma :
## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## v_sigma 22.9 6.45 3.54 0.203 13.1 22.5 34 1006 1
##
## sigma_ :
## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## sigma_ 56.2 3.04 18.5 0.0697 51.3 56.2 61.3 1904 0.999
##
## v :
## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## Alameda -25.3766 15.2 -1.66784 0.351 -51.2 -24.794 -0.931 1881 1.000
## Amador -0.8011 23.3 -0.03433 0.449 -38.8 -1.037 38.259 2699 1.000
## Butte 0.4312 24.0 0.01799 0.438 -38.0 -0.382 41.339 3000 0.999
## Calaveras 7.8602 21.4 0.36759 0.399 -26.3 7.152 42.990 2877 1.000
## Colusa -0.0848 23.4 -0.00362 0.452 -38.8 -0.118 37.954 2683 1.000
## Contra Costa -9.6494 15.1 -0.63804 0.291 -34.9 -9.476 14.215 2703 1.001
## Del Norte -0.2623 23.5 -0.01116 0.429 -39.5 -0.135 38.639 3000 1.000
## El Dorado 0.0559 23.6 0.00237 0.430 -36.5 -0.538 39.818 3000 1.001
## Fresno 0.2615 15.1 0.01737 0.279 -24.5 0.449 24.674 2902 1.000
## Glenn 0.1972 23.6 0.00836 0.439 -38.6 0.322 38.751 2879 1.002
## ... 47 elements suppressed ...We wish to estimate the area population means \[
\theta_i = \frac{1}{N_i}\sum_{j \in U_i} y_j\,,
\] where \(U_i\) is the set of
units in area \(i\) of size \(N_i\). The MCMC output in variable
sim can be used to obtain draws from the posterior
distribution for \(\theta_i\). The
\(r\)th draw can be expressed as \[
\theta_{i;r} = \frac{1}{N_i} \left(n_i \bar{y}_i + \beta_r'(t_{x;i}
- n_i \bar{x}_i) + (N_i - n_i)v_{i;r} + \sum_{j \in U_i\setminus s_i}
\epsilon_{j;r} \right)\,,
\] where \(\bar{y}_i\) is the
sample mean of \(y\) in area \(i\) and \(t_{x;i}\) is a vector of population totals
for area \(i\).
N <- table(apipop$cname)
samplesums <- tapply(apisrs$api00, apisrs$cname, sum)
samplesums[is.na(samplesums)] <- 0 # substitute 0 for out-of-sample areas
m <- match(apisrs$cds, apipop$cds) # population units in the sample
res <- predict(sim, newdata=apipop, labels=names(N),
fun=function(x, p) (samplesums + tapply(x[-m], apipop$cname[-m], sum ))/N,
show.progress=FALSE)
(summ <- summary(res))## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## Alameda 679 14.54 46.7 0.335 654 679 702 1885 1.000
## Amador 737 30.79 24.0 0.591 689 737 789 2714 0.999
## Butte 680 26.20 26.0 0.478 638 679 725 3000 1.000
## Calaveras 743 27.76 26.8 0.541 697 742 788 2634 1.001
## Colusa 552 31.03 17.8 0.605 500 551 603 2629 0.999
## Contra Costa 727 14.63 49.7 0.279 702 727 750 2747 1.001
## Del Norte 678 31.88 21.3 0.582 627 678 732 3000 0.999
## El Dorado 754 26.49 28.5 0.484 712 754 798 2989 1.002
## Fresno 595 14.89 39.9 0.272 570 595 619 3000 1.000
## Glenn 625 30.84 20.3 0.563 574 625 676 3000 1.001
## Humboldt 701 26.56 26.4 0.498 659 701 745 2840 1.000
## Imperial 557 24.75 22.5 0.507 518 557 598 2385 1.000
## Inyo 706 32.66 21.6 0.596 653 706 760 3000 1.000
## Kern 597 14.83 40.2 0.298 572 597 620 2472 1.001
## Kings 595 22.66 26.3 0.445 556 596 631 2593 0.999
## Lake 662 24.81 26.7 0.472 623 662 703 2769 1.000
## Lassen 706 25.55 27.7 0.466 664 706 748 3000 0.999
## Los Angeles 636 8.16 77.9 0.149 622 635 649 3000 1.000
## Madera 614 20.08 30.6 0.367 582 614 649 3000 1.000
## Marin 817 20.37 40.1 0.379 783 816 850 2880 1.001
## Mariposa 724 35.42 20.4 0.653 666 724 781 2945 1.000
## Mendocino 663 26.57 24.9 0.487 619 663 706 2972 1.000
## Merced 578 23.28 24.9 0.426 541 578 617 2979 1.000
## Modoc 671 30.21 22.2 0.594 623 671 721 2584 0.999
## Mono 705 42.49 16.6 0.776 634 705 775 3000 0.999
## Monterey 628 18.09 34.7 0.330 600 628 658 3000 1.000
## Napa 701 22.07 31.8 0.403 665 701 737 3000 1.000
## Nevada 798 29.35 27.2 0.536 751 799 847 3000 1.000
## Orange 692 15.52 44.6 0.312 667 692 718 2468 1.000
## Placer 772 23.34 33.1 0.426 736 772 811 3000 1.000
## Plumas 689 31.48 21.9 0.575 636 689 740 3000 1.000
## Riverside 636 14.71 43.2 0.296 611 636 660 2462 1.000
## Sacramento 668 15.23 43.9 0.278 644 669 694 3000 1.001
## San Benito 709 29.83 23.8 0.545 662 709 757 3000 0.999
## San Bernardino 647 12.83 50.4 0.234 626 647 668 3000 1.000
## San Diego 704 14.53 48.5 0.325 682 704 729 1998 1.001
## San Francisco 639 19.47 32.8 0.355 608 639 671 3000 0.999
## San Joaquin 630 17.62 35.8 0.390 600 631 658 2040 1.000
## San Luis Obispo 752 23.38 32.2 0.427 714 752 791 3000 0.999
## San Mateo 733 21.32 34.4 0.400 699 733 769 2845 1.000
## Santa Barbara 679 19.69 34.5 0.370 647 679 712 2835 1.000
## Santa Clara 733 15.98 45.9 0.292 706 733 759 3000 1.000
## Santa Cruz 679 19.42 35.0 0.374 648 680 711 2694 0.999
## Shasta 696 21.57 32.3 0.438 662 695 731 2431 1.000
## Sierra 707 41.95 16.9 0.766 639 707 775 3000 1.001
## Siskiyou 697 25.61 27.2 0.482 655 697 737 2819 1.001
## Solano 711 19.72 36.0 0.390 677 712 741 2563 0.999
## Sonoma 724 22.95 31.6 0.419 685 724 762 3000 1.000
## Stanislaus 660 18.88 35.0 0.349 630 661 691 2918 0.999
## Sutter 651 24.39 26.7 0.446 610 652 690 2994 1.001
## Tehama 660 28.59 23.1 0.530 613 660 708 2907 0.999
## Trinity 649 37.95 17.1 0.693 587 648 711 3000 1.000
## Tulare 583 21.42 27.2 0.449 547 583 616 2279 1.001
## Tuolumne 720 31.05 23.2 0.574 670 720 773 2930 0.999
## Ventura 709 17.51 40.5 0.323 682 709 739 2941 0.999
## Yolo 687 23.79 28.9 0.442 647 687 725 2902 1.000
## Yuba 626 27.84 22.5 0.508 581 626 672 3000 0.999theta <- c(tapply(apipop$api00, apipop$cname, mean)) # true population quantities
plot_coef(summ, list(est=theta), n.se=2, est.names=c("mcmcsae", "true"), maxrows=30)
A model with binomial likelihood can also be fit. We now model the
target variable sch.wide, a binary variable indicating
whether a school-wide growth target has been met. We use the same mean
model structure as above for the linear model, but now using a logistic
link function, \[
y_j \stackrel{\mathrm{iid}}{\sim} {\cal Be}(p_j)\,,\\
\mathrm{logit}(p_j) = \beta' x_j + v_{i[j]}\,,\\
v_i \stackrel{\mathrm{iid}}{\sim} {\cal N}(0, \sigma_v^2)
\]
apisrs$target.met <- as.numeric(apisrs$sch.wide == "Yes")
sampler <- create_sampler(update(mod, target.met ~ .), family="binomial", data=apisrs)
sim <- MCMCsim(sampler, store.all=TRUE, verbose=FALSE)
summary(sim)## llh_ :
## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## llh_ -83.4 2.65 -31.5 0.0704 -88.1 -83.2 -79.4 1415 1
##
## beta :
## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## (Intercept) 4.552953 1.5298 2.9762 0.054750 2.1507 4.46779 7.26596 781 1.006
## ell -0.035655 0.0150 -2.3750 0.000421 -0.0603 -0.03520 -0.01160 1271 0.999
## meals 0.000552 0.0140 0.0395 0.000400 -0.0222 0.00038 0.02386 1222 1.003
## stypeH -2.813428 0.6216 -4.5261 0.019353 -3.8435 -2.80857 -1.80443 1032 1.002
## stypeM -1.638721 0.5631 -2.9101 0.017515 -2.5650 -1.64377 -0.71009 1034 1.003
## hsg -0.025371 0.0224 -1.1312 0.000664 -0.0628 -0.02480 0.01041 1141 1.005
## col.grad -0.037490 0.0270 -1.3904 0.000849 -0.0815 -0.03746 0.00656 1009 1.004
## grad.sch 0.014183 0.0319 0.4444 0.001074 -0.0371 0.01203 0.06822 883 1.003
##
## v_sigma :
## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## v_sigma 0.369 0.275 1.34 0.0129 0.0269 0.316 0.91 453 1
##
## v :
## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## Alameda -1.83e-01 0.419 -0.437536 0.01170 -1.062 -7.46e-02 0.340 1284 1.002
## Amador 2.37e-03 0.476 0.004974 0.00924 -0.744 1.34e-03 0.747 2654 1.001
## Butte 1.19e-02 0.458 0.026060 0.00843 -0.718 3.60e-05 0.778 2944 0.999
## Calaveras 3.52e-02 0.452 0.077953 0.00859 -0.651 7.31e-03 0.817 2769 0.999
## Colusa -2.84e-04 0.472 -0.000601 0.00926 -0.755 1.64e-04 0.729 2590 1.000
## Contra Costa 1.93e-02 0.397 0.048508 0.00850 -0.635 4.06e-03 0.694 2186 1.000
## Del Norte -3.41e-03 0.458 -0.007441 0.00837 -0.739 -6.15e-04 0.731 2998 1.000
## El Dorado -5.68e-05 0.454 -0.000125 0.00936 -0.746 7.88e-05 0.741 2357 1.000
## Fresno -4.42e-02 0.392 -0.112620 0.00808 -0.735 -1.21e-02 0.561 2358 1.000
## Glenn -5.14e-03 0.464 -0.011073 0.00847 -0.765 -6.38e-05 0.751 3000 1.000
## ... 47 elements suppressed ...To predict the population fractions of schools that meet the growth target by county,
samplesums <- tapply(apisrs$target.met, apisrs$cname, sum)
samplesums[is.na(samplesums)] <- 0 # substitute 0 for out-of-sample areas
res <- predict(sim, newdata=apipop, labels=names(N),
fun=function(x, p) (samplesums + tapply(x[-m], apipop$cname[-m], sum ))/N,
show.progress=FALSE)
(summ <- summary(res))## Mean SD t-value MCSE q0.05 q0.5 q0.95 n_eff R_hat
## Alameda 0.785 0.0598 13.13 0.00155 0.670 0.792 0.867 1492 1.001
## Amador 0.842 0.1175 7.16 0.00235 0.600 0.900 1.000 2492 1.001
## Butte 0.835 0.0743 11.25 0.00141 0.708 0.833 0.938 2770 1.000
## Calaveras 0.874 0.0944 9.25 0.00191 0.700 0.900 1.000 2439 1.001
## Colusa 0.647 0.1611 4.01 0.00294 0.333 0.667 0.889 3000 0.999
## Contra Costa 0.818 0.0569 14.37 0.00129 0.721 0.821 0.905 1951 1.001
## Del Norte 0.879 0.1102 7.98 0.00205 0.625 0.875 1.000 2900 1.000
## El Dorado 0.850 0.0752 11.31 0.00152 0.725 0.850 0.950 2451 1.001
## Fresno 0.782 0.0627 12.46 0.00139 0.672 0.785 0.871 2045 1.000
## Glenn 0.797 0.1394 5.71 0.00255 0.556 0.778 1.000 3000 0.999
## Humboldt 0.852 0.0745 11.44 0.00154 0.725 0.850 0.950 2332 1.000
## Imperial 0.698 0.1016 6.87 0.00209 0.525 0.700 0.850 2369 1.001
## Inyo 0.776 0.1564 4.96 0.00301 0.571 0.857 1.000 2704 1.000
## Kern 0.803 0.0534 15.03 0.00117 0.706 0.811 0.878 2091 0.999
## Kings 0.850 0.0815 10.43 0.00176 0.720 0.840 0.960 2155 1.000
## Lake 0.821 0.0956 8.59 0.00195 0.636 0.818 0.955 2410 1.000
## Lassen 0.836 0.1086 7.70 0.00228 0.636 0.818 1.000 2265 1.000
## Los Angeles 0.798 0.0436 18.30 0.00112 0.724 0.798 0.869 1525 1.000
## Madera 0.820 0.0737 11.13 0.00145 0.677 0.839 0.935 2575 1.000
## Marin 0.837 0.0723 11.58 0.00176 0.700 0.840 0.940 1689 1.001
## Mariposa 0.861 0.1491 5.77 0.00272 0.600 0.800 1.000 3000 0.999
## Mendocino 0.815 0.0906 9.00 0.00182 0.640 0.840 0.960 2465 1.001
## Merced 0.787 0.0815 9.65 0.00171 0.651 0.794 0.905 2278 1.000
## Modoc 0.834 0.1568 5.32 0.00306 0.600 0.800 1.000 2627 1.000
## Mono 0.744 0.2449 3.04 0.00460 0.333 0.667 1.000 2838 1.001
## Monterey 0.805 0.0694 11.60 0.00160 0.687 0.807 0.916 1885 1.002
## Napa 0.843 0.0791 10.66 0.00166 0.704 0.852 0.963 2256 1.000
## Nevada 0.869 0.0967 8.99 0.00193 0.714 0.857 1.000 2517 1.000
## Orange 0.773 0.0614 12.59 0.00135 0.672 0.775 0.871 2068 0.999
## Placer 0.813 0.0779 10.43 0.00176 0.682 0.818 0.909 1949 0.999
## Plumas 0.746 0.1530 4.87 0.00284 0.444 0.778 1.000 2913 1.002
## Riverside 0.823 0.0542 15.18 0.00125 0.726 0.831 0.895 1875 1.000
## Sacramento 0.846 0.0489 17.31 0.00119 0.764 0.847 0.927 1688 1.000
## San Benito 0.819 0.1268 6.46 0.00232 0.545 0.818 1.000 2984 1.000
## San Bernardino 0.862 0.0434 19.87 0.00128 0.790 0.862 0.931 1150 1.001
## San Diego 0.848 0.0432 19.65 0.00110 0.775 0.850 0.916 1549 0.999
## San Francisco 0.762 0.0740 10.29 0.00150 0.630 0.770 0.870 2440 1.000
## San Joaquin 0.828 0.0573 14.46 0.00127 0.733 0.831 0.911 2036 1.001
## San Luis Obispo 0.855 0.0713 11.99 0.00145 0.725 0.850 0.950 2422 1.001
## San Mateo 0.798 0.0758 10.53 0.00183 0.664 0.811 0.902 1721 1.001
## Santa Barbara 0.833 0.0637 13.07 0.00141 0.716 0.840 0.926 2028 1.002
## Santa Clara 0.783 0.0709 11.04 0.00210 0.652 0.796 0.875 1141 1.003
## Santa Cruz 0.790 0.0745 10.59 0.00146 0.653 0.796 0.898 2609 1.000
## Shasta 0.815 0.0797 10.22 0.00171 0.675 0.825 0.925 2175 1.000
## Sierra 0.729 0.2311 3.15 0.00450 0.333 0.667 1.000 2635 1.001
## Siskiyou 0.844 0.1013 8.32 0.00198 0.667 0.867 1.000 2631 1.001
## Solano 0.864 0.0638 13.54 0.00144 0.750 0.867 0.967 1976 1.000
## Sonoma 0.844 0.0663 12.74 0.00143 0.731 0.852 0.944 2144 1.000
## Stanislaus 0.813 0.0725 11.21 0.00188 0.670 0.824 0.912 1493 1.000
## Sutter 0.833 0.0909 9.16 0.00173 0.700 0.850 0.950 2748 1.001
## Tehama 0.844 0.0954 8.85 0.00180 0.647 0.882 1.000 2798 1.000
## Trinity 0.753 0.2008 3.75 0.00379 0.500 0.750 1.000 2813 1.000
## Tulare 0.831 0.0638 13.03 0.00133 0.718 0.836 0.918 2302 0.999
## Tuolumne 0.879 0.0946 9.30 0.00182 0.750 0.917 1.000 2692 0.999
## Ventura 0.840 0.0510 16.45 0.00106 0.752 0.845 0.919 2336 1.001
## Yolo 0.850 0.0755 11.25 0.00152 0.714 0.857 0.971 2470 1.000
## Yuba 0.793 0.1011 7.85 0.00192 0.632 0.789 0.947 2780 1.001theta <- c(tapply(apipop$sch.wide == "Yes", apipop$cname, mean)) # true population quantities
plot_coef(summ, list(est=theta), n.se=2, est.names=c("mcmcsae", "true"), maxrows=30)