20190726nigata-test-simple

set.seed(1234)

nigata data

read table & rounding to integers

dat <- read.csv("nigata.csv")
colnames(dat) <- c("city", "sosu", "tokuhyo")
data = list(N=nrow(dat), x=round(dat$tokuhyo), n=round(dat$sosu), city=as.character(dat$city))

remove jouetsu from data : THIS IS DIFFERENT FROM THE PREVIOUS ANALYSIS

In data2 (N, n, x), joetsu is removed from data. The original data is now renamed as N_new, n_new, x_new. N, n, x are used for computing posterior distribution. Then predictive distribution is computed for n_new, and posterior probability is computed at x_new.

id2 = (1:data$N)[-2] # 2=joetsu
data2 = list(N=data$N-1, x=data$x[id2], n=data$n[id2], id=id2, city=data$city[id2], N_new=data$N, x_new=data$x, n_new=data$n, city_new=data$city)

plot

plot(data2$n_new, data2$x_new, type="n", xlab="投票総数", ylab="得票数", family=family); text(data2$n_new, data2$x_new, cex=0.5,col="red") ; text(data2$n, data2$x, data2$id, cex=0.5)

Computing posterior distributions (Bayesian inference)

stan

library(rstan)

## Loading required package: StanHeaders

## Loading required package: ggplot2

## rstan (Version 2.19.2, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores()) 

fit

betabinom3 = stan_model("betabinom3.stan")
fit1 <- sampling(betabinom3, data=data2, chains=8, warmup=2000, iter=20000)

stan diagnostics

fit1

## Inference for Stan model: betabinom3.
## 8 chains, each with iter=20000; warmup=2000; thin=1; 
## post-warmup draws per chain=18000, total post-warmup draws=144000.
## 
##                mean se_mean     sd      2.5%       25%       50%       75%
## p[1]           0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[2]           0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[3]           0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[4]           0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[5]           0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[6]           0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[7]           0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[8]           0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[9]           0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[10]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[11]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[12]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[13]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[14]          0.01    0.00   0.00      0.00      0.00      0.01      0.01
## p[15]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[16]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[17]          0.01    0.00   0.00      0.00      0.00      0.01      0.01
## p[18]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[19]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[20]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[21]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[22]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[23]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[24]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[25]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[26]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[27]          0.01    0.00   0.00      0.00      0.01      0.01      0.01
## p[28]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[29]          0.01    0.00   0.00      0.00      0.01      0.01      0.01
## p[30]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[31]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[32]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[33]          0.01    0.00   0.00      0.00      0.00      0.01      0.01
## p[34]          0.01    0.00   0.00      0.00      0.01      0.01      0.01
## p[35]          0.01    0.00   0.00      0.01      0.01      0.01      0.01
## p[36]          0.01    0.00   0.00      0.00      0.01      0.01      0.01
## phi            0.01    0.00   0.00      0.01      0.01      0.01      0.01
## lambda      2975.21    2.77 879.03   1573.74   2347.95   2868.29   3482.28
## alpha         23.58    0.02   7.03     12.42     18.56     22.70     27.61
## beta        2951.63    2.75 872.07   1561.07   2329.48   2845.77   3454.83
## x_rep[1]     934.65    0.54 205.69    572.00    793.00    921.00   1060.00
## x_rep[2]     684.96    0.40 151.32    418.00    581.00    674.00    777.00
## x_rep[3]     597.47    0.35 132.34    364.00    506.00    588.00    678.00
## x_rep[4]     525.36    0.31 116.62    319.00    445.00    517.00    597.00
## x_rep[5]     429.31    0.25  95.68    260.00    363.00    423.00    488.00
## x_rep[6]     348.40    0.21  78.09    211.00    294.00    343.00    396.00
## x_rep[7]     328.85    0.19  73.71    199.00    278.00    324.00    374.00
## x_rep[8]     313.24    0.18  70.44    189.00    265.00    308.00    356.00
## x_rep[9]     285.63    0.17  64.39    172.00    241.00    281.00    325.00
## x_rep[10]    258.42    0.15  58.38    155.00    218.00    254.00    294.00
## x_rep[11]    234.28    0.14  53.22    140.00    197.00    231.00    267.00
## x_rep[12]    223.64    0.13  50.89    134.00    189.00    220.00    255.00
## x_rep[13]    218.83    0.13  49.78    131.00    184.00    215.00    249.00
## x_rep[14]    218.38    0.13  49.76    130.00    184.00    215.00    249.00
## x_rep[15]    214.09    0.13  48.83    128.00    180.00    211.00    244.00
## x_rep[16]    213.71    0.13  48.73    128.00    180.00    210.00    244.00
## x_rep[17]    189.07    0.11  43.34    113.00    159.00    186.00    216.00
## x_rep[18]    180.81    0.11  41.55    107.00    152.00    178.00    206.00
## x_rep[19]    171.81    0.10  39.59    102.00    144.00    169.00    196.00
## x_rep[20]    154.86    0.09  35.95     91.00    130.00    152.00    177.00
## x_rep[21]    151.56    0.09  35.21     90.00    127.00    149.00    173.00
## x_rep[22]    145.21    0.09  33.82     86.00    122.00    143.00    166.00
## x_rep[23]    144.32    0.09  33.61     85.00    121.00    142.00    165.00
## x_rep[24]    138.00    0.08  32.26     81.00    116.00    136.00    158.00
## x_rep[25]    132.42    0.08  31.00     78.00    111.00    130.00    151.00
## x_rep[26]    111.47    0.07  26.40     65.00     93.00    110.00    128.00
## x_rep[27]     99.69    0.06  23.86     58.00     83.00     98.00    114.00
## x_rep[28]     49.76    0.03  12.90     27.00     41.00     49.00     58.00
## x_rep[29]     43.99    0.03  11.63     24.00     36.00     43.00     51.00
## x_rep[30]     43.93    0.03  11.59     24.00     36.00     43.00     51.00
## x_rep[31]     42.88    0.03  11.39     23.00     35.00     42.00     50.00
## x_rep[32]     30.76    0.02   8.67     16.00     25.00     30.00     36.00
## x_rep[33]     29.68    0.02   8.46     15.00     24.00     29.00     35.00
## x_rep[34]     25.93    0.02   7.61     13.00     21.00     25.00     31.00
## x_rep[35]     19.83    0.02   6.19      9.00     15.00     19.00     24.00
## x_rep[36]     18.43    0.02   5.87      8.00     14.00     18.00     22.00
## x_rep[37]      1.79    0.00   1.39      0.00      1.00      2.00      3.00
## lp__      -45408.98    0.02   4.53 -45418.77 -45411.83 -45408.64 -45405.76
##               97.5%  n_eff Rhat
## p[1]           0.01 201705    1
## p[2]           0.01 187530    1
## p[3]           0.01 197439    1
## p[4]           0.01 200015    1
## p[5]           0.01 205708    1
## p[6]           0.01 204975    1
## p[7]           0.01 198228    1
## p[8]           0.01 193472    1
## p[9]           0.01 197279    1
## p[10]          0.01 193216    1
## p[11]          0.01 198090    1
## p[12]          0.01 203409    1
## p[13]          0.01 187908    1
## p[14]          0.01 187676    1
## p[15]          0.01 199284    1
## p[16]          0.01 196178    1
## p[17]          0.01 174493    1
## p[18]          0.01 185510    1
## p[19]          0.01 197305    1
## p[20]          0.01 196491    1
## p[21]          0.01 201517    1
## p[22]          0.01 193835    1
## p[23]          0.01 197673    1
## p[24]          0.01 206973    1
## p[25]          0.01 204791    1
## p[26]          0.01 223610    1
## p[27]          0.01 173517    1
## p[28]          0.01 177553    1
## p[29]          0.01 185722    1
## p[30]          0.01 189775    1
## p[31]          0.01 205773    1
## p[32]          0.01 197978    1
## p[33]          0.01 162531    1
## p[34]          0.01 161178    1
## p[35]          0.01 192443    1
## p[36]          0.01 188522    1
## phi            0.01 140068    1
## lambda      5002.43 100365    1
## alpha         39.78  97174    1
## beta        4962.91 100397    1
## x_rep[1]    1379.00 145301    1
## x_rep[2]    1014.00 145101    1
## x_rep[3]     885.00 145363    1
## x_rep[4]     779.00 145845    1
## x_rep[5]     636.00 145183    1
## x_rep[6]     518.00 145001    1
## x_rep[7]     489.00 144927    1
## x_rep[8]     467.00 145737    1
## x_rep[9]     425.00 145420    1
## x_rep[10]    385.00 145620    1
## x_rep[11]    350.00 145306    1
## x_rep[12]    334.00 145574    1
## x_rep[13]    326.00 145668    1
## x_rep[14]    326.00 145406    1
## x_rep[15]    320.00 145210    1
## x_rep[16]    319.00 145012    1
## x_rep[17]    283.00 145327    1
## x_rep[18]    271.00 145111    1
## x_rep[19]    257.00 145561    1
## x_rep[20]    233.00 144994    1
## x_rep[21]    228.00 145832    1
## x_rep[22]    219.00 145311    1
## x_rep[23]    217.00 145422    1
## x_rep[24]    208.00 145234    1
## x_rep[25]    200.00 144683    1
## x_rep[26]    169.00 145394    1
## x_rep[27]    152.00 145942    1
## x_rep[28]     78.00 143969    1
## x_rep[29]     69.00 144857    1
## x_rep[30]     69.00 145239    1
## x_rep[31]     68.00 145570    1
## x_rep[32]     49.00 142616    1
## x_rep[33]     48.00 143532    1
## x_rep[34]     42.00 144524    1
## x_rep[35]     33.00 145362    1
## x_rep[36]     31.00 143997    1
## x_rep[37]      5.00 142930    1
## lp__      -45401.10  57086    1
## 
## Samples were drawn using NUTS(diag_e) at Fri Jul 26 11:29:36 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

stan_trace(fit1, pars=c("phi", "lambda"))

stan_dens(fit1, pars=c("phi", "lambda"), separate_chains = TRUE)

#stan_hist(fit1)

stan_dens(fit1, pars=c("x_rep[1]","x_rep[2]","x_rep[3]","x_rep[4]"))

stan summary

ss = summary(fit1)$summary
ss

##                    mean      se_mean           sd          2.5%
## p[1]       9.014125e-03 6.047820e-07 2.716173e-04  8.489705e-03
## p[2]       1.088730e-02 8.584954e-07 3.717687e-04  1.017161e-02
## p[3]       1.083061e-02 8.888983e-07 3.949737e-04  1.007132e-02
## p[4]       9.505385e-03 9.095543e-07 4.067806e-04  8.723831e-03
## p[5]       9.571992e-03 9.937461e-07 4.507140e-04  8.708909e-03
## p[6]       8.463164e-03 9.608331e-07 4.350092e-04  7.631932e-03
## p[7]       7.395671e-03 9.370294e-07 4.171916e-04  6.597665e-03
## p[8]       9.029976e-03 1.085402e-06 4.774193e-04  8.115497e-03
## p[9]       9.369038e-03 1.151764e-06 5.115691e-04  8.396851e-03
## p[10]      9.784958e-03 1.248195e-06 5.486611e-04  8.742600e-03
## p[11]      6.810190e-03 1.055494e-06 4.697721e-04  5.916651e-03
## p[12]      7.892500e-03 1.122259e-06 5.061483e-04  6.930823e-03
## p[13]      6.433121e-03 1.066913e-06 4.624896e-04  5.557581e-03
## p[14]      5.082705e-03 9.690159e-07 4.197932e-04  4.288485e-03
## p[15]      8.532146e-03 1.193813e-06 5.329333e-04  7.520261e-03
## p[16]      9.075030e-03 1.309257e-06 5.798957e-04  7.976521e-03
## p[17]      5.253056e-03 1.104175e-06 4.612398e-04  4.384631e-03
## p[18]      7.931471e-03 1.305217e-06 5.621681e-04  6.867704e-03
## p[19]      7.792560e-03 1.318194e-06 5.855284e-04  6.684692e-03
## p[20]      6.766097e-03 1.253092e-06 5.554622e-04  5.719078e-03
## p[21]      9.463598e-03 1.492466e-06 6.699777e-04  8.202937e-03
## p[22]      7.526555e-03 1.353468e-06 5.958866e-04  6.396627e-03
## p[23]      9.832852e-03 1.563379e-06 6.950854e-04  8.519390e-03
## p[24]      7.187088e-03 1.332732e-06 6.063169e-04  6.045204e-03
## p[25]      7.541208e-03 1.471516e-06 6.659176e-04  6.289811e-03
## p[26]      8.069500e-03 1.527482e-06 7.223072e-04  6.717121e-03
## p[27]      5.768927e-03 1.959588e-06 8.162730e-04  4.255261e-03
## p[28]      6.740283e-03 2.137767e-06 9.007918e-04  5.074705e-03
## p[29]      6.390397e-03 2.059192e-06 8.874193e-04  4.739664e-03
## p[30]      7.692500e-03 2.220376e-06 9.672674e-04  5.903528e-03
## p[31]      8.102813e-03 2.417865e-06 1.096798e-03  6.108200e-03
## p[32]      8.429501e-03 2.539376e-06 1.129889e-03  6.368611e-03
## p[33]      5.648253e-03 2.502904e-06 1.009048e-03  3.785766e-03
## p[34]      6.271659e-03 2.774713e-06 1.113964e-03  4.210043e-03
## p[35]      7.255532e-03 2.719721e-06 1.193096e-03  5.071030e-03
## p[36]      7.662140e-03 3.724939e-06 1.617336e-03  4.768977e-03
## phi        7.925731e-03 8.409554e-07 3.147329e-04  7.310684e-03
## lambda     2.975211e+03 2.774674e+00 8.790310e+02  1.573744e+03
## alpha      2.357606e+01 2.253723e-02 7.025485e+00  1.242333e+01
## beta       2.951635e+03 2.752258e+00 8.720651e+02  1.561068e+03
## x_rep[1]   9.346461e+02 5.396123e-01 2.056917e+02  5.720000e+02
## x_rep[2]   6.849638e+02 3.972577e-01 1.513241e+02  4.180000e+02
## x_rep[3]   5.974714e+02 3.471135e-01 1.323421e+02  3.640000e+02
## x_rep[4]   5.253576e+02 3.053729e-01 1.166208e+02  3.190000e+02
## x_rep[5]   4.293095e+02 2.510987e-01 9.567598e+01  2.600000e+02
## x_rep[6]   3.484045e+02 2.050854e-01 7.809446e+01  2.110000e+02
## x_rep[7]   3.288479e+02 1.936110e-01 7.370636e+01  1.990000e+02
## x_rep[8]   3.132413e+02 1.845231e-01 7.044260e+01  1.890000e+02
## x_rep[9]   2.856332e+02 1.688596e-01 6.439285e+01  1.720000e+02
## x_rep[10]  2.584196e+02 1.529761e-01 5.837598e+01  1.550000e+02
## x_rep[11]  2.342781e+02 1.396110e-01 5.321843e+01  1.400000e+02
## x_rep[12]  2.236362e+02 1.333741e-01 5.088770e+01  1.340000e+02
## x_rep[13]  2.188325e+02 1.304340e-01 4.978208e+01  1.310000e+02
## x_rep[14]  2.183837e+02 1.304977e-01 4.976149e+01  1.300000e+02
## x_rep[15]  2.140903e+02 1.281489e-01 4.883296e+01  1.280000e+02
## x_rep[16]  2.137063e+02 1.279600e-01 4.872776e+01  1.280000e+02
## x_rep[17]  1.890697e+02 1.136911e-01 4.334112e+01  1.130000e+02
## x_rep[18]  1.808051e+02 1.090828e-01 4.155332e+01  1.070000e+02
## x_rep[19]  1.718127e+02 1.037607e-01 3.958734e+01  1.020000e+02
## x_rep[20]  1.548582e+02 9.440499e-02 3.594765e+01  9.100000e+01
## x_rep[21]  1.515598e+02 9.219739e-02 3.520828e+01  9.000000e+01
## x_rep[22]  1.452092e+02 8.871714e-02 3.381866e+01  8.600000e+01
## x_rep[23]  1.443209e+02 8.812943e-02 3.360750e+01  8.500000e+01
## x_rep[24]  1.379951e+02 8.464218e-02 3.225678e+01  8.100000e+01
## x_rep[25]  1.324224e+02 8.150268e-02 3.100139e+01  7.800000e+01
## x_rep[26]  1.114741e+02 6.924808e-02 2.640471e+01  6.500000e+01
## x_rep[27]  9.969236e+01 6.245457e-02 2.385915e+01  5.800000e+01
## x_rep[28]  4.976072e+01 3.400258e-02 1.290169e+01  2.700000e+01
## x_rep[29]  4.398594e+01 3.055099e-02 1.162774e+01  2.400000e+01
## x_rep[30]  4.393278e+01 3.040001e-02 1.158553e+01  2.400000e+01
## x_rep[31]  4.287799e+01 2.984996e-02 1.138886e+01  2.300000e+01
## x_rep[32]  3.075542e+01 2.295771e-02 8.669887e+00  1.600000e+01
## x_rep[33]  2.968007e+01 2.232727e-02 8.458825e+00  1.500000e+01
## x_rep[34]  2.593451e+01 2.000492e-02 7.605142e+00  1.300000e+01
## x_rep[35]  1.982784e+01 1.623348e-02 6.189244e+00  9.000000e+00
## x_rep[36]  1.843465e+01 1.546368e-02 5.867995e+00  8.000000e+00
## x_rep[37]  1.790222e+00 3.669648e-03 1.387351e+00  0.000000e+00
## lp__      -4.540898e+04 1.897896e-02 4.534592e+00 -4.541877e+04
##                     25%           50%           75%         97.5%
## p[1]       8.830368e-03  9.010907e-03  9.196356e-03  9.552163e-03
## p[2]       1.063365e-02  1.088348e-02  1.113509e-02  1.162855e-02
## p[3]       1.056124e-02  1.082493e-02  1.109425e-02  1.162019e-02
## p[4]       9.227282e-03  9.499817e-03  9.777248e-03  1.031908e-02
## p[5]       9.265562e-03  9.564905e-03  9.871040e-03  1.047928e-02
## p[6]       8.165800e-03  8.456552e-03  8.753388e-03  9.336106e-03
## p[7]       7.110797e-03  7.387813e-03  7.673186e-03  8.233882e-03
## p[8]       8.706105e-03  9.020089e-03  9.346461e-03  9.990246e-03
## p[9]       9.017348e-03  9.357525e-03  9.709323e-03  1.039791e-02
## p[10]      9.412110e-03  9.774337e-03  1.014737e-02  1.089316e-02
## p[11]      6.488270e-03  6.799820e-03  7.121516e-03  7.757993e-03
## p[12]      7.547387e-03  7.879693e-03  8.225249e-03  8.920728e-03
## p[13]      6.117014e-03  6.419973e-03  6.739043e-03  7.373331e-03
## p[14]      4.792977e-03  5.073071e-03  5.361291e-03  5.932243e-03
## p[15]      8.167344e-03  8.521295e-03  8.883244e-03  9.611564e-03
## p[16]      8.677969e-03  9.062317e-03  9.456497e-03  1.025013e-02
## p[17]      4.935548e-03  5.240101e-03  5.555557e-03  6.195713e-03
## p[18]      7.544220e-03  7.917445e-03  8.303362e-03  9.066888e-03
## p[19]      7.390949e-03  7.778358e-03  8.178559e-03  8.988891e-03
## p[20]      6.385289e-03  6.751524e-03  7.131483e-03  7.895197e-03
## p[21]      9.004751e-03  9.446088e-03  9.906971e-03  1.081652e-02
## p[22]      7.116027e-03  7.512539e-03  7.919700e-03  8.740399e-03
## p[23]      9.356657e-03  9.814404e-03  1.029030e-02  1.124453e-02
## p[24]      6.769393e-03  7.171191e-03  7.588396e-03  8.417493e-03
## p[25]      7.083997e-03  7.523215e-03  7.979963e-03  8.901687e-03
## p[26]      7.571333e-03  8.046377e-03  8.546056e-03  9.545303e-03
## p[27]      5.201507e-03  5.740579e-03  6.302644e-03  7.447919e-03
## p[28]      6.116527e-03  6.702375e-03  7.327930e-03  8.599438e-03
## p[29]      5.777119e-03  6.358403e-03  6.968147e-03  8.222638e-03
## p[30]      7.025944e-03  7.654816e-03  8.316402e-03  9.706916e-03
## p[31]      7.337551e-03  8.051453e-03  8.809959e-03  1.039823e-02
## p[32]      7.645972e-03  8.376820e-03  9.153915e-03  1.080916e-02
## p[33]      4.946541e-03  5.605173e-03  6.304480e-03  7.739673e-03
## p[34]      5.502919e-03  6.228793e-03  6.990979e-03  8.585814e-03
## p[35]      6.434459e-03  7.201208e-03  8.018661e-03  9.764371e-03
## p[36]      6.549362e-03  7.562986e-03  8.668385e-03  1.114562e-02
## phi        7.716481e-03  7.923148e-03  8.132728e-03  8.556005e-03
## lambda     2.347950e+03  2.868286e+03  3.482281e+03  5.002432e+03
## alpha      1.855524e+01  2.269620e+01  2.761085e+01  3.978465e+01
## beta       2.329477e+03  2.845768e+03  3.454831e+03  4.962910e+03
## x_rep[1]   7.930000e+02  9.210000e+02  1.060000e+03  1.379000e+03
## x_rep[2]   5.810000e+02  6.740000e+02  7.770000e+02  1.014000e+03
## x_rep[3]   5.060000e+02  5.880000e+02  6.780000e+02  8.850000e+02
## x_rep[4]   4.450000e+02  5.170000e+02  5.970000e+02  7.790000e+02
## x_rep[5]   3.630000e+02  4.230000e+02  4.880000e+02  6.360000e+02
## x_rep[6]   2.940000e+02  3.430000e+02  3.960000e+02  5.180000e+02
## x_rep[7]   2.780000e+02  3.240000e+02  3.740000e+02  4.890000e+02
## x_rep[8]   2.650000e+02  3.080000e+02  3.560000e+02  4.670000e+02
## x_rep[9]   2.410000e+02  2.810000e+02  3.250000e+02  4.250000e+02
## x_rep[10]  2.180000e+02  2.540000e+02  2.940000e+02  3.850000e+02
## x_rep[11]  1.970000e+02  2.310000e+02  2.670000e+02  3.500000e+02
## x_rep[12]  1.890000e+02  2.200000e+02  2.550000e+02  3.340000e+02
## x_rep[13]  1.840000e+02  2.150000e+02  2.490000e+02  3.260000e+02
## x_rep[14]  1.840000e+02  2.150000e+02  2.490000e+02  3.260000e+02
## x_rep[15]  1.800000e+02  2.110000e+02  2.440000e+02  3.200000e+02
## x_rep[16]  1.800000e+02  2.100000e+02  2.440000e+02  3.190000e+02
## x_rep[17]  1.590000e+02  1.860000e+02  2.160000e+02  2.830000e+02
## x_rep[18]  1.520000e+02  1.780000e+02  2.060000e+02  2.710000e+02
## x_rep[19]  1.440000e+02  1.690000e+02  1.960000e+02  2.570000e+02
## x_rep[20]  1.300000e+02  1.520000e+02  1.770000e+02  2.330000e+02
## x_rep[21]  1.270000e+02  1.490000e+02  1.730000e+02  2.280000e+02
## x_rep[22]  1.220000e+02  1.430000e+02  1.660000e+02  2.190000e+02
## x_rep[23]  1.210000e+02  1.420000e+02  1.650000e+02  2.170000e+02
## x_rep[24]  1.160000e+02  1.360000e+02  1.580000e+02  2.080000e+02
## x_rep[25]  1.110000e+02  1.300000e+02  1.510000e+02  2.000000e+02
## x_rep[26]  9.300000e+01  1.100000e+02  1.280000e+02  1.690000e+02
## x_rep[27]  8.300000e+01  9.800000e+01  1.140000e+02  1.520000e+02
## x_rep[28]  4.100000e+01  4.900000e+01  5.800000e+01  7.800000e+01
## x_rep[29]  3.600000e+01  4.300000e+01  5.100000e+01  6.900000e+01
## x_rep[30]  3.600000e+01  4.300000e+01  5.100000e+01  6.900000e+01
## x_rep[31]  3.500000e+01  4.200000e+01  5.000000e+01  6.800000e+01
## x_rep[32]  2.500000e+01  3.000000e+01  3.600000e+01  4.900000e+01
## x_rep[33]  2.400000e+01  2.900000e+01  3.500000e+01  4.800000e+01
## x_rep[34]  2.100000e+01  2.500000e+01  3.100000e+01  4.200000e+01
## x_rep[35]  1.500000e+01  1.900000e+01  2.400000e+01  3.300000e+01
## x_rep[36]  1.400000e+01  1.800000e+01  2.200000e+01  3.100000e+01
## x_rep[37]  1.000000e+00  2.000000e+00  3.000000e+00  5.000000e+00
## lp__      -4.541183e+04 -4.540864e+04 -4.540576e+04 -4.540110e+04
##               n_eff      Rhat
## p[1]      201705.24 1.0000322
## p[2]      187529.58 0.9999753
## p[3]      197438.64 0.9999890
## p[4]      200015.29 0.9999904
## p[5]      205708.04 0.9999806
## p[6]      204975.06 0.9999863
## p[7]      198227.89 0.9999676
## p[8]      193472.16 1.0000063
## p[9]      197279.17 0.9999754
## p[10]     193216.21 0.9999665
## p[11]     198090.11 0.9999786
## p[12]     203408.73 0.9999746
## p[13]     187908.07 1.0000297
## p[14]     187676.13 0.9999876
## p[15]     199284.47 0.9999759
## p[16]     196177.70 0.9999723
## p[17]     174492.74 0.9999701
## p[18]     185509.85 0.9999955
## p[19]     197304.56 0.9999829
## p[20]     196491.34 0.9999910
## p[21]     201517.17 0.9999935
## p[22]     193834.60 0.9999848
## p[23]     197673.22 0.9999963
## p[24]     206972.92 0.9999879
## p[25]     204791.03 0.9999873
## p[26]     223610.34 0.9999674
## p[27]     173516.69 1.0000528
## p[28]     177553.05 1.0000076
## p[29]     185722.23 0.9999791
## p[30]     189775.44 0.9999938
## p[31]     205773.17 0.9999672
## p[32]     197978.44 0.9999708
## p[33]     162530.74 1.0000280
## p[34]     161178.09 1.0000414
## p[35]     192442.97 0.9999880
## p[36]     188521.91 0.9999732
## phi       140067.82 0.9999807
## lambda    100365.48 1.0001570
## alpha      97174.31 1.0001614
## beta      100396.71 1.0001569
## x_rep[1]  145301.42 1.0000178
## x_rep[2]  145101.41 1.0000242
## x_rep[3]  145362.81 1.0000213
## x_rep[4]  145844.81 1.0000336
## x_rep[5]  145183.40 1.0000112
## x_rep[6]  145001.02 1.0000209
## x_rep[7]  144927.18 1.0000182
## x_rep[8]  145736.87 1.0000068
## x_rep[9]  145419.94 1.0000025
## x_rep[10] 145620.08 1.0000270
## x_rep[11] 145306.37 1.0000187
## x_rep[12] 145573.59 1.0000377
## x_rep[13] 145668.06 1.0000075
## x_rep[14] 145405.52 1.0000157
## x_rep[15] 145209.91 1.0000127
## x_rep[16] 145012.25 1.0000111
## x_rep[17] 145327.34 1.0000062
## x_rep[18] 145110.64 1.0000004
## x_rep[19] 145561.50 1.0000224
## x_rep[20] 144994.32 1.0000215
## x_rep[21] 145831.86 1.0000054
## x_rep[22] 145310.80 1.0000202
## x_rep[23] 145422.11 1.0000079
## x_rep[24] 145234.00 1.0000049
## x_rep[25] 144683.30 1.0000044
## x_rep[26] 145394.26 1.0000081
## x_rep[27] 145942.41 1.0000028
## x_rep[28] 143969.17 1.0000290
## x_rep[29] 144857.22 1.0000082
## x_rep[30] 145239.31 1.0000183
## x_rep[31] 145570.36 1.0000044
## x_rep[32] 142616.44 1.0000001
## x_rep[33] 143532.03 1.0000104
## x_rep[34] 144524.28 0.9999892
## x_rep[35] 145362.44 1.0000229
## x_rep[36] 143997.15 1.0000319
## x_rep[37] 142930.13 0.9999999
## lp__       57086.29 1.0001649

parameter estimates of beta distribution (for p)

phi1 = ss["phi","mean"]  # posterior mean of phi
lambda1 = ss["lambda","mean"]  # posterior mean of lambda
alpha1 = lambda1*phi1
beta1 = lambda1*(1-phi1)
beta_mean <- function(alpha,beta) alpha/(alpha+beta)
beta_var <- function(alpha,beta) alpha*beta/((alpha+beta)^2*(alpha+beta+1))

phi1

## [1] 0.007925731

lambda1

## [1] 2975.211

alpha1

## [1] 23.58072

beta1

## [1] 2951.63

beta_mean(alpha1,beta1) # mean of p

## [1] 0.007925731

sqrt(beta_var(alpha1,beta1)) # sd of p

## [1] 0.001625399

pp = seq(0,1,length=1000)
plot(pp,dbeta(pp,alpha1,beta1),type="l") # distribution of p

extract posterior distributions (used later)

ex1 <- rstan::extract(fit1)
names(ex1)

## [1] "p"      "phi"    "lambda" "alpha"  "beta"   "x_rep"  "lp__"

results

data points and posterior predictive distribution

p0=sum(data2$x)/sum(data2$n)
plot(data2$n_new, data2$x_new, type="n", xlab="投票総数", ylab="得票数", main="予測分布の95%区間", sub=paste("平均得票率=",signif(p0,4)),family=family)
#abline(a=0,b=phi1, col="red") # parameter estimates (=mean)
abline(a=0,b=p0, col="red") # = x/n
x_conf <- ss[paste("x_rep[",1:data2$N_new,"]",sep=""), c("2.5%","97.5%","50%","mean")]
segments(data2$n_new,x_conf[,1],data2$n_new,x_conf[,2],col="green")
points(data2$n_new, x_conf[,4],col="green" )
text(data2$n_new, data2$x_new, cex=0.5,col="red")
text(data2$n, data2$x, data2$id, cex=0.5)

tail probability of posterior predictive distribution

ptail1 <- c()
for(i in 1:data2$N_new) ptail1[i] = mean(ex1$x_rep[,i] <= data2$x_new[i])  # lower tail
ptail1 <- 2*pmin(ptail1,1-ptail1) # posterior tail probability (two-tailed)
ptail10 <- ptail1; ptail10[ptail10==0] <- 1e-10 # set default for p=0
hist(ptail10,nclass=10,xlab="posterior probability")

plot of posterior tail probability. joetsu is actually p=0, but set to 10^-10 for log-plot.

plot(data2$n_new, ptail10, log="y", type="n", xlab="投票総数", ylab="事後確率", family=family)
text(data2$n_new, ptail10);
abline(h=0.05,col="blue") # 5%
abline(h=0.01,col="orange") # 1%
abline(h=0.001,col="red")  # 0.1%

VGAM includes beta-binomial distribution

library(VGAM)

## Loading required package: stats4

## Loading required package: splines

tail probabilities of beta-binomial using the estimated parameter value

ptail2 = pbetabinom.ab(data2$x_new,data2$n_new,alpha1,beta1)
ptail2 <- 2*pmin(ptail2,1-ptail2) # p-value (two-tailed)
hist(ptail2,nclass=10,xlab="p-value") # check if uniform on (0,1)

plot(ptail1,ptail2); abline(0,1)

plot of frequntist p-value (based on the estimated parameters by stan)

plot(data2$n_new, ptail2, log="y", type="n", xlab="投票総数", ylab="p値", family=family)
text(data2$n_new, ptail2);
abline(h=0.05,col="blue") # 5%
abline(h=0.01,col="orange") # 1%
abline(h=0.001,col="red")  # 0.1%

summary table of p-values

phi1

## [1] 0.007925731

lambda1

## [1] 2975.211

alpha1 

## [1] 23.58072

beta1

## [1] 2951.63

cbind(as.character(data2$city_new), ptail1, ptail2)

##                  ptail1               ptail2                
##  [1,] "長岡市"   "0.479194444444444"  "0.472491596698264"   
##  [2,] "上越市"   "0"                  "3.84307810362392e-13"
##  [3,] "中央区"   "0.0979305555555556" "0.0830259951349011"  
##  [4,] "西区"     "0.102222222222222"  "0.0875131645362417"  
##  [5,] "東区"     "0.323833333333333"  "0.312612075219217"   
##  [6,] "三条市"   "0.303625"           "0.292350568065012"   
##  [7,] "新発田市" "0.682541666666667"  "0.682643215995566"   
##  [8,] "柏崎市"   "0.802347222222222"  "0.787952360177446"   
##  [9,] "燕市"     "0.459125"           "0.45304508396282"    
## [10,] "秋葉区"   "0.351583333333333"  "0.342320051743463"   
## [11,] "北区"     "0.248791666666667"  "0.234470201532035"   
## [12,] "村上市"   "0.513541666666667"  "0.496923224748353"   
## [13,] "南魚沼市" "0.947013888888889"  "0.953477616159549"   
## [14,] "十日町市" "0.358486111111111"  "0.339528111162337"   
## [15,] "佐渡市"   "0.0540833333333333" "0.0398253901415737"  
## [16,] "江南区"   "0.6475"             "0.645718478123923"   
## [17,] "西蒲区"   "0.436180555555556"  "0.428638473908783"   
## [18,] "糸魚川市" "0.0710416666666667" "0.0542853276023904"  
## [19,] "五泉市"   "0.927694444444444"  "0.932609616522397"   
## [20,] "阿賀野市" "0.997916666666667"  "0.995927313483905"   
## [21,] "魚沼市"   "0.487277777777778"  "0.468273520181109"   
## [22,] "見附市"   "0.313069444444444"  "0.300627661967117"   
## [23,] "小千谷市" "0.865388888888889"  "0.853733911733157"   
## [24,] "南区"     "0.222319444444445"  "0.209959191963094"   
## [25,] "妙高市"   "0.684652777777778"  "0.672860017636591"   
## [26,] "胎内市"   "0.869208333333333"  "0.860823511877214"   
## [27,] "加茂市"   "0.846222222222222"  "0.851237843641421"   
## [28,] "阿賀町"   "0.106847222222222"  "0.090162036337366"   
## [29,] "聖籠町"   "0.419888888888889"  "0.403367177244446"   
## [30,] "田上町"   "0.267972222222222"  "0.252648839998035"   
## [31,] "津南町"   "0.958888888888889"  "0.952102393525368"   
## [32,] "弥彦村"   "0.783333333333333"  "0.785525116519834"   
## [33,] "湯沢町"   "0.603736111111111"  "0.601169781524318"   
## [34,] "関川村"   "0.0455"             "0.037783651748282"   
## [35,] "刈羽村"   "0.149527777777778"  "0.137286834725477"   
## [36,] "出雲崎町" "0.655402777777778"  "0.64361905083257"    
## [37,] "粟島浦村" "0.950083333333333"  "0.946225963028565"

fitting the ordinary binomial distribution (THIS IS NOT A VALID ANALYSIS)

p

(x = sum(data2$x))

## [1] 7927

(n = sum(data2$n))

## [1] 917973

(p= x/n)

## [1] 0.00863533

tail probabilities of observed values (p-values) : DONT USE THIS!!!

ptail3 = pbinom(data2$x_new,data2$n_new,p)
ptail3 <- 2*pmin(ptail3,1-ptail3) # p-value (two-tailed) with Binomial
hist(ptail3,nclass=10,xlab="p-value (Binomial)") # check if uniform on (0,1)

plot(ptail2,ptail3); abline(0,1)

ptail3

##  [1]  1.292703e-01 7.634938e-213  1.318079e-11  4.471641e-10  1.648690e-02
##  [6]  1.808143e-02  7.947882e-01  5.423240e-03  3.047340e-01  8.898368e-02
## [11]  1.374743e-02  2.876663e-04  1.893969e-01  1.051758e-05  9.536274e-14
## [16]  9.806163e-01  3.151734e-01  5.054898e-11  2.774301e-01  2.030293e-01
## [21]  1.675758e-03  1.115310e-01  8.983691e-02  3.050654e-02  2.660637e-02
## [26]  1.393036e-01  5.576124e-01  4.542911e-04  4.246570e-02  1.201152e-02
## [31]  4.461637e-01  8.787387e-01  8.172285e-01  9.106200e-04  1.811491e-02
## [36]  3.002932e-01  8.362501e-01

posterior predictive distributions of all x’s

for(i in 1:data2$N_new) {
  hist(ex1$x_rep[,i], nclass=100, xlab="得票数", main=paste("事後予測分布 (",data2$city_new[i],")"), sub=paste(i,":", data2$city_new[i],", 事後確率=", signif(ptail1[i],4), ", p値=", signif(ptail2[i],4)), xlim=range(c(ex1$x_rep[,i],0,data2$x_new[i])), family= family) # histogram of posterior distribution
  abline(v = data2$x_new[i], col="red" ) # observed value
}