1 synthパッケージ

#install.packages("Synth")

library(Synth)

2 データ

  • 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
    • regionno==1はスペイン全土の平均。
  • 42の期間(1955年~1997年)。

3 データ準備 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)  # プロットの時間範囲
    )

4 データの確認 (省略可)

# 制御群の予測データの行列。
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

5 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
                          ) 

6 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"),
          ) 

7 ギャップ確認: 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,])
                                                 )