synth
パッケージ
#install.packages("Synth")
library(Synth)
データ
- 18単位からなるパネル・データフレーム:1つの処置地域(No.17;バスク地方)と16の対照地域(No.2-16,18)。
- 1つの結果変数(gdpcap):
一人当たり実質GDP(1986年、千米ドル)。
- 13の予測変数(6つの部門別生産シェア、6つの最高学歴カテゴリー、人口密度、投資率)。
- 6つの部門別生産シェア
- sec.agriculture: 農林水産業の生産高が総生産高に占める割合。
- sec.energy: 総生産に占めるエネルギー・水部門の生産比率。
- sec.industry: 総生産に占める工業部門の生産。
- sec.construction:
建設・エンジニアリング部門の生産高が総生産高に占める割合。
- sec.services.venta:
総生産量に占める市場性サービス部門の生産量の割合。
- sec.services.nonventa:
非市場サービス部門の生産高が総生産高に占める割合。
- 6つの最高学歴カテゴリー
- school.illit: 非識字者数。
- school.prim: 初等教育を受けた人の数。
- school.med: 高等教育を受けた人の数。
- school.high: 高等学校卒業者数
- school.post.high: 高等教育を受けた人の数。
- popdens: 人口密度(1平方キロメートルあたりの人口)
- invest: 総投資額がGDPに占める割合
- 地域名と地域番号は、regionno と regionname
- 42の期間(1955年~1997年)。
データ準備
dataprep()
# データの読み込み
data(basque)
# synthのためのデータの準備
dataprep.out <-
dataprep(
foo = basque # データ
,predictors= c("school.illit", # 予測変数1、非識字者数
"school.prim", # 予測変数2、初等教育
"school.med", # 予測変数3、高等教育
"school.high", # 予測変数4、高等学校卒業
"school.post.high" # 予測変数5、高等教育
,"invest" # 予測変数6、総投資額/GDP
)
,predictors.op = c("mean") # 予測変数の操作
,dependent = c("gdpcap") # 従属変数
,unit.variable = c("regionno") # 単位の変数
,time.variable = c("year") # 時間の変数
,special.predictors = list( # 特別な予測変数
list("gdpcap",1960:1969,c("mean")),
list("sec.agriculture",seq(1961,1969,2),c("mean")), # 農業部門の生産シェア
list("sec.energy",seq(1961,1969,2),c("mean")), # エネルギー部門の生産シェア
list("sec.industry",seq(1961,1969,2),c("mean")), # 工業部門の生産シェア
list("sec.construction",seq(1961,1969,2),c("mean")), # 建設部門の生産シェア
list("sec.services.venta",seq(1961,1969,2),c("mean")), # 市場性サービス部門の生産シェア
list("sec.services.nonventa",seq(1961,1969,2),c("mean")), # 非市場性サービス部門の生産シェア
list("popdens",1969,c("mean"))) # 人口密度
,treatment.identifier = 17 # 処置地域
,controls.identifier = c(2:16,18) # 対照地域
,time.predictors.prior = c(1964:1969) # プレ期間の時間範囲
,time.optimize.ssr = c(1960:1969) # 最適化の時間範囲
,unit.names.variable = c("regionname") # 単位名
,time.plot = c(1955:1997) # プロットの時間範囲
)
データの確認
(省略可)
# 制御群の予測データの行列。
dataprep.out[["X0"]]
## 2 3 4
## school.illit 863.389160 73.121226 31.488423
## school.prim 3062.424886 728.578929 670.909393
## school.med 155.565318 44.215389 46.398482
## school.high 57.266496 16.091676 14.799358
## school.post.high 27.278924 8.684416 6.424505
## invest 19.320031 21.577486 22.769643
## special.gdpcap.1960.1969 2.560747 3.699907 3.733876
## special.sec.agriculture.1961.1969 24.194000 21.726000 12.362000
## special.sec.energy.1961.1969 2.774000 6.278000 18.648000
## special.sec.industry.1961.1969 18.276000 22.780000 24.126000
## special.sec.construction.1961.1969 8.130000 7.832000 9.006000
## special.sec.services.venta.1961.1969 38.186000 34.289999 30.152000
## special.sec.services.nonventa.1961.1969 8.444000 7.096000 5.708000
## special.popdens.1969 68.510002 24.040001 98.739998
## 5 6 7
## school.illit 47.903906 128.308287 9.394911
## school.prim 300.813619 522.094955 289.571732
## school.med 20.045204 45.447489 24.106414
## school.high 5.921604 12.086849 7.304420
## school.post.high 3.680154 5.844122 2.885214
## invest 24.441712 25.954247 29.071211
## special.gdpcap.1960.1969 5.215974 3.051014 3.871173
## special.sec.agriculture.1961.1969 13.130000 19.944000 15.922000
## special.sec.energy.1961.1969 2.076000 7.818000 2.894000
## special.sec.industry.1961.1969 18.258000 9.816000 36.530000
## special.sec.construction.1961.1969 8.294000 8.670000 5.976000
## special.sec.services.venta.1961.1969 51.752000 45.278001 33.484000
## special.sec.services.nonventa.1961.1969 6.494000 8.482000 5.198000
## special.popdens.1969 104.169998 148.250000 87.389999
## 8 9 10
## school.illit 105.508144 254.449707 277.935237
## school.prim 1766.928691 1035.269368 2883.361735
## school.med 102.299929 33.369543 215.165476
## school.high 37.787317 16.539727 63.608315
## school.post.high 19.173427 6.308165 32.374611
## invest 19.652509 17.970690 22.520380
## special.gdpcap.1960.1969 2.807062 2.240688 5.147008
## special.sec.agriculture.1961.1969 29.718000 36.086001 6.936000
## special.sec.energy.1961.1969 7.818000 5.180000 2.880000
## special.sec.industry.1961.1969 16.474000 17.178000 40.056000
## special.sec.construction.1961.1969 7.144000 5.810000 6.706000
## special.sec.services.venta.1961.1969 30.844000 29.238000 39.068000
## special.sec.services.nonventa.1961.1969 8.000000 6.506000 4.356000
## special.popdens.1969 28.770000 22.379999 153.119995
## 11 12 13
## school.illit 266.967601 177.727966 234.752001
## school.prim 1695.117126 692.888550 1666.930420
## school.med 103.959094 24.048450 77.464478
## school.high 34.275772 11.271676 29.341941
## school.post.high 16.822106 4.570566 13.543060
## invest 23.702696 20.239484 20.606269
## special.gdpcap.1960.1969 3.843245 1.926478 2.516288
## special.sec.agriculture.1961.1969 18.582000 37.948000 28.862000
## special.sec.energy.1961.1969 2.698000 2.646000 5.328000
## special.sec.industry.1961.1969 28.046000 10.886000 17.882000
## special.sec.construction.1961.1969 6.170000 8.770000 7.032000
## special.sec.services.venta.1961.1969 39.064001 31.720000 33.714001
## special.sec.services.nonventa.1961.1969 5.444000 8.038000 7.188000
## special.popdens.1969 128.699997 28.969999 91.019997
## 14 15 16
## school.illit 133.158962 105.877481 13.438050
## school.prim 1856.074280 431.746806 280.002396
## school.med 269.961647 27.907434 20.450936
## school.high 62.458658 9.082540 6.550353
## school.post.high 57.704985 4.428528 4.190368
## invest 16.234543 20.706271 19.955159
## special.gdpcap.1960.1969 5.976692 2.899557 3.996087
## special.sec.agriculture.1961.1969 1.864000 19.352000 24.704000
## special.sec.energy.1961.1969 2.074000 10.228000 2.740000
## special.sec.industry.1961.1969 23.834000 19.644000 28.398000
## special.sec.construction.1961.1969 8.358000 6.290000 6.982000
## special.sec.services.venta.1961.1969 52.714000 36.177999 30.468000
## special.sec.services.nonventa.1961.1969 11.162000 8.312000 6.710000
## special.popdens.1969 442.450012 73.360001 44.169998
## 18
## school.illit 9.151832
## school.prim 152.269465
## school.med 9.760035
## school.high 3.377170
## school.post.high 1.732763
## invest 18.055932
## special.gdpcap.1960.1969 3.809205
## special.sec.agriculture.1961.1969 30.322001
## special.sec.energy.1961.1969 2.888000
## special.sec.industry.1961.1969 26.612000
## special.sec.construction.1961.1969 5.238000
## special.sec.services.venta.1961.1969 28.300000
## special.sec.services.nonventa.1961.1969 6.644000
## special.popdens.1969 46.580002
# 処置群の予測データの行列。
dataprep.out[["X1"]]
## 17
## school.illit 39.888465
## school.prim 1031.742299
## school.med 90.358668
## school.high 25.727525
## school.post.high 13.479720
## invest 24.647383
## special.gdpcap.1960.1969 5.285468
## special.sec.agriculture.1961.1969 6.844000
## special.sec.energy.1961.1969 4.106000
## special.sec.industry.1961.1969 45.082000
## special.sec.construction.1961.1969 6.150000
## special.sec.services.venta.1961.1969 33.754000
## special.sec.services.nonventa.1961.1969 4.072000
## special.popdens.1969 246.889999
# 処置前期間の処置群の結果データの行列
dataprep.out[["Z1"]]
## 17
## 1960 4.285918
## 1961 4.574336
## 1962 4.898957
## 1963 5.197015
## 1964 5.338903
## 1965 5.465153
## 1966 5.545916
## 1967 5.614896
## 1968 5.852185
## 1969 6.081405
# 処置前期間の制御群の結果データの行列
dataprep.out[["Z0"]]
## 2 3 4 5 6 7 8 9
## 1960 2.010140 2.881462 2.967295 4.058841 2.357684 3.137032 2.138817 1.667524
## 1961 2.129177 3.099543 3.143887 4.360254 2.445730 3.327621 2.239503 1.752428
## 1962 2.280348 3.359183 3.373536 4.646173 2.648243 3.555341 2.454227 1.920451
## 1963 2.431020 3.614182 3.597258 4.911525 2.844759 3.771423 2.672237 2.091902
## 1964 2.508855 3.680091 3.672594 5.050700 2.951157 3.839403 2.777778 2.182591
## 1965 2.584690 3.745287 3.743359 5.184662 3.054199 3.906098 2.882176 2.274707
## 1966 2.694444 3.883319 3.909383 5.466795 3.231791 4.032133 2.988075 2.378392
## 1967 2.802342 4.016138 4.073122 5.737646 3.403385 4.155955 3.094544 2.482362
## 1968 2.987361 4.243645 4.308626 6.161454 3.660312 4.375893 3.302271 2.709083
## 1969 3.179092 4.476221 4.549700 6.581691 3.912882 4.610826 3.520994 2.947444
## 10 11 12 13 14 15 16 18
## 1960 4.241788 3.219294 1.535847 1.983290 5.161097 2.118609 3.163525 2.969866
## 1961 4.575335 3.362468 1.596258 2.005784 5.632605 2.305484 3.335904 3.153171
## 1962 4.838046 3.569980 1.705584 2.185661 5.840831 2.521422 3.623393 3.404384
## 1963 5.081334 3.765210 1.817695 2.366395 6.024493 2.739074 3.894816 3.669238
## 1964 5.158098 3.823693 1.882819 2.458797 6.099329 2.851257 3.985147 3.803985
## 1965 5.223651 3.874179 1.948872 2.549700 6.152028 2.965938 4.072979 3.921808
## 1966 5.332477 3.978149 2.032633 2.669666 6.110469 3.099186 4.210011 4.032705
## 1967 5.429449 4.073408 2.117609 2.787846 6.057341 3.227292 4.352399 4.160311
## 1968 5.674379 4.279777 2.245501 2.978363 6.253142 3.461154 4.556984 4.373036
## 1969 5.915524 4.486290 2.381962 3.177378 6.435590 3.706155 4.765710 4.603542
synthを実行
synth()
# synthを実行
synth.out <- synth(data.prep.obj = dataprep.out)
##
## X1, X0, Z1, Z0 all come directly from dataprep object.
##
##
## ****************
## searching for synthetic control unit
##
##
## ****************
## ****************
## ****************
##
## MSPE (LOSS V): 0.008864621
##
## solution.v:
## 0.009024072 0.00872398 3e-10 2.53552e-05 9.644e-07 2.5177e-06 0.01072469 0.3076797 0.02974812 0.2928763 0.002138064 2.9514e-06 0.02188925 0.3171641
##
## solution.w:
## 3.27e-08 5.46e-08 7.45e-08 3.85e-08 3.87e-08 6.654e-07 4.3e-08 5.52e-08 0.8508223 1.271e-07 3.69e-08 6.15e-08 0.1491762 5.39e-08 1.27e-07 1.452e-07
# 結果を得る
synth.tables <- synth.tab(
dataprep.res = dataprep.out,
synth.res = synth.out
)
synthのプロット
path.plot
# results tables:
print(synth.tables)
## $tab.pred
## Treated Synthetic Sample Mean
## school.illit 39.888 256.338 170.786
## school.prim 1031.742 2730.111 1127.186
## school.med 90.359 223.339 76.260
## school.high 25.728 63.437 24.235
## school.post.high 13.480 36.153 13.478
## invest 24.647 21.583 21.424
## special.gdpcap.1960.1969 5.285 5.271 3.581
## special.sec.agriculture.1961.1969 6.844 6.179 21.353
## special.sec.energy.1961.1969 4.106 2.760 5.310
## special.sec.industry.1961.1969 45.082 37.636 22.425
## special.sec.construction.1961.1969 6.150 6.952 7.276
## special.sec.services.venta.1961.1969 33.754 41.104 36.528
## special.sec.services.nonventa.1961.1969 4.072 5.371 7.111
## special.popdens.1969 246.890 196.281 99.414
##
## $tab.v
## v.weights
## school.illit 0.009
## school.prim 0.009
## school.med 0
## school.high 0
## school.post.high 0
## invest 0
## special.gdpcap.1960.1969 0.011
## special.sec.agriculture.1961.1969 0.308
## special.sec.energy.1961.1969 0.03
## special.sec.industry.1961.1969 0.293
## special.sec.construction.1961.1969 0.002
## special.sec.services.venta.1961.1969 0
## special.sec.services.nonventa.1961.1969 0.022
## special.popdens.1969 0.317
##
## $tab.w
## w.weights unit.names unit.numbers
## 2 0.000 Andalucia 2
## 3 0.000 Aragon 3
## 4 0.000 Principado De Asturias 4
## 5 0.000 Baleares (Islas) 5
## 6 0.000 Canarias 6
## 7 0.000 Cantabria 7
## 8 0.000 Castilla Y Leon 8
## 9 0.000 Castilla-La Mancha 9
## 10 0.851 Cataluna 10
## 11 0.000 Comunidad Valenciana 11
## 12 0.000 Extremadura 12
## 13 0.000 Galicia 13
## 14 0.149 Madrid (Comunidad De) 14
## 15 0.000 Murcia (Region de) 15
## 16 0.000 Navarra (Comunidad Foral De) 16
## 18 0.000 Rioja (La) 18
##
## $tab.loss
## Loss W Loss V
## [1,] 0.3070471 0.008864621
# plot results:
# path
path.plot(synth.res = synth.out,
dataprep.res = dataprep.out,
Ylab = c("real per-capita GDP (1986 USD, thousand)"),
Xlab = c("year"),
Ylim = c(0,13),
Legend = c("Basque country","synthetic Basque country"),
)

ギャップ確認:
gaps.plot()
## gaps
gaps.plot(synth.res = synth.out,
dataprep.res = dataprep.out,
Ylab = c("gap in real per-capita GDP (1986 USD, thousand)"),
Xlab = c("year"),
Ylim = c(-1.5,1.5),
)

補論
- CRANのマニュアルでは、以下のデータ加工をsynth実行前に行っている。
- 本質的には不要だが、計算を早くするためと考えられる。
データの加工(1)
# 1. 最高学歴と第2学歴を統合
# 最高学歴のカテゴリーを廃止し、第2学歴のカテゴリーに統合
dataprep.out$X1["school.high",] <-
dataprep.out$X1["school.high",] +
dataprep.out$X1["school.post.high",]
# which(rownames(dataprep.out$X1)=="school.post.high") は、dataprep.out$X1 の行名が "school.post.high" である行のインデックスを取得。
# マイナスがついているので、"school.post.high" でない行を抽出。
# dataprep.out$X1[...] は、この抽出した行のデータフレームを返す。
dataprep.out$X1 <-
as.matrix(dataprep.out$X1[-which(rownames(dataprep.out$X1)=="school.post.high"),])
#
dataprep.out$X0["school.high",] <-
dataprep.out$X0["school.high",] +
dataprep.out$X0["school.post.high",]
dataprep.out$X0 <-
dataprep.out$X0[-which(rownames(dataprep.out$X0)=="school.post.high"),]
データの加工(2)
# 2. 合計を作成し、各就学カテゴリーのシェアを計算。
# rownames(dataprep.out$X0) は、dataprep.out$X0 の行名(行の識別子)のベクトルを返す。
# which(rownames(dataprep.out$X0)=="school.illit") は、そのベクトル中で "school.illit" という行名に一致する要素のインデックスを取得。
lowest <- which(rownames(dataprep.out$X0)=="school.illit")
highest <- which(rownames(dataprep.out$X0)=="school.high")
# 各就学カテゴリーのシェアを計算
dataprep.out$X1[lowest:highest,] <-
(100 * dataprep.out$X1[lowest:highest,]) / sum(dataprep.out$X1[lowest:highest,])
#
dataprep.out$X0[lowest:highest,] <-
100 * scale(dataprep.out$X0[lowest:highest,],
center=FALSE,
scale=colSums(dataprep.out$X0[lowest:highest,])
)