英国総選挙の期間中、ロンドン大政治学科の選挙ブログ(LSE:
Election Experts blog)で、
Simon Hixと
Nick Vivyanの2人の研究者が、各社の世論調査を統合して、3党の得票率を推定する記事を書いていた。
1社では飛び飛びのグラフになってしまうが、コンピュータによる統合計算で、得票率の連続的な推移と誤差と推定している。
2人は、この推計値から各選挙区の得票率を推計し、当落を判断しているのだが、最も興味深いのは、最初の「各社の世論調査を統合する」作業だ。
採用された手法は、コンピュータと豪州政治を研究しているスタンフォード大学
Simon Jackmanによる
Polling the Polls。各調査主体の癖があることを組み込んだモデルを提案している。
データ:2004年の豪州総選挙(5社の世論調査)
この件については、Jackman本の最終章の最後に書いてある。
同じ手法を日本に当てはめた例:
社会科学のベイズ統計 動的線形モデルDLM
【モデル】
- i : 世論調査のナンバー
- t : 日付(最初の調査を第1日とし、第T日目に投開票日をむかえる)
- yi : i番目の調査による政党の相対得票率(調査結果)
- αt : t日における政党の相対得票率(これが知りたい)
- μi=αt(i)+δj(i) : 世論調査の期待値はα(得票率)とδ(調査主体jによる癖、house effect)
- σ2i=yi(1-yi)/ni : サンプル数nの調査の分散期待値
- yi〜N(μi,σ2i) : 調査結果は正規分布
さらに以下のような仮定をする。
- αt〜N(αt-1,ω2) : 政党得票率は分散ω(これ自体が分散0.1の雑音)でランダムウオークする
- α1〜Uniform(0.4,0.6) : 初期値は多分40-60%で一様分布
- αT=開票結果 : 最終的にこの結果にならねばならない
このモデルをコンピュータで計算する。JAGSで使うBUGS言語では、モデルは以下のように記述される。順序はほとんど関係ない。
model{
## NPOLLS個の世論調査について
for(i in 1:NPOLLS){
## μはalphaとhouse(効果)の和である
mu[i] <- house[org[i]] + alpha[date[i]]
## 世論調査結果はμとprecで決まる正規分布
y[i] ~ dnorm(mu[i],prec[i])
}
## 推移モデル
## すべての期間について
for(i in 2:NPERIODS){ ## 一つ前のalphaを基準にランダムウオークする
mu.alpha[i] <- alpha[i-1]
alpha[i] ~ dnorm(mu.alpha[i],tau)
}
## 事前確率
sigma ~ dunif(0,.01) ## sigmaは0−1の一様乱数
tau <- 1/pow(sigma,2) ## ランダムウオークの分散
alpha[1] ~ dunif(.4,.6) ## 初期値は40-60%の一様乱数
for(i in 1:5){ ## 調査主体(この場合は5社)に応じた偏差
house[i] ~ dnorm(0,.01)
}
}
このモデルをどれだけ走らせるかを指定するのがスクリプトファイルで
model in kalman.bug #モデルファイルを指定
data in dumpdata.R #データファイルを指定
compile #コンパイル
initialize #初期化
update 1000 #最初の1000回(収束まで破棄する)
monitor alpha, thin(500) #alphaを500回ごとに保存
monitor sigma, thin(500)
monitor house, thin(500)
update 25000 #25000回繰り返す
coda * #monitorした値を書きだす
dumpdata.Rには(一部はダミーの)データをRで使う形式で準備しておく。
y <-
c(0.48, 0.48, 0.47, 0.5, 0.51, 0.54, 0.52, 0.54, 0.51, 0.54,
0.5, 0.52, 0.51, 0.52, 0.52, 0.48, 0.49, 0.49, 0.5, 0.46, 0.48,
0.5, 0.5, 0.475, 0.48, 0.505, 0.5, 0.46, 0.485, 0.46, 0.47, 0.465,
0.445, 0.44, 0.455, 0.47, 0.5, 0.485, 0.49, 0.46, 0.47, 0.48,
0.51, 0.49, 0.5, 0.525, 0.51, 0.47, 0.44, 0.45, 0.49, 0.46, 0.48,
0.48, 0.5, 0.47, 0.42, 0.46, 0.46, 0.47, 0.45, 0.455, 0.46, 0.46,
0.46, 0.46)
prec <-
c(5685.09615384615, 5665.0641025641, 5676.43516659976, 5656,
5658.26330532213, 5704.50885668277, 5596.95512820513, 8168.2769726248,
3985.59423769508, 4009.66183574879, 3984, 4006.41025641026, 4041.61664665866,
4006.41025641026, 4807.69230769231, 4627.40384615385, 4589.83593437375,
4553.82152861144, 4560, 4609.50080515298, 4587.33974358974, 6936,
6828, 6712.78195488722, 6814.90384615385, 6720.67206720672, 10000,
5088.56682769726, 4383.9455509959, 4247.18196457327, 8145.3231633882,
7673.60064315144, 7798.36015791072, 7573.05194805195, 3762.47605605404,
....
実行によって、CODAindex.txtとCODAchain1.txtが生成されるので、これを
codaライブラリのread.coda.interactive()で読み込むとmcmc.listクラスを得られる。このクラスには、データを目で点検するためのplot機能が組み込まれている。
例えば、house effectのシミュレーションの推移を見るには
res <- read.coda.interactive()
plot(res[,116:120])
によって分布が得られ、
summary(alpha[,116:120])
によって
Iterations = 1001:25501
Thinning interval = 500
Number of chains = 1
Sample size per chain = 50
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
house[1] -0.01141 0.007502 0.001061 0.001233
house[2] -0.05751 0.007933 0.001122 0.001026
house[3] -0.00205 0.009187 0.001299 0.001108
house[4] -0.02906 0.009271 0.001311 0.001570
house[5] -0.04893 0.007936 0.001122 0.001236
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
house[1] -0.02404 -0.01708 -0.011387 -0.005933 0.003332
house[2] -0.07109 -0.06355 -0.057064 -0.051590 -0.043512
house[3] -0.02038 -0.00767 -0.002666 0.004141 0.014595
house[4] -0.04683 -0.03259 -0.028943 -0.023893 -0.010237
house[5] -0.06397 -0.05327 -0.048362 -0.043595 -0.034841
が得られる。これをどう解釈するかは別の話。本人がエセックス大で行った
報告によると、上のデータから、house1(ニールセン)とhouse3(ギャラクシー)が偏りがなく、他3社は過大評価していることになる。