StatsBeginner: 初学者の統計学習ノート

統計学およびR、Pythonでのプログラミングの勉強の過程をメモっていくノート。たまにMacの話題。

Rでピボットテーブル

学生の分析を手伝っていたところ、突然、Rでピボットテーブルみたいな集計をしたくなりました。
度数分布表はtable関数でつくれますが、ピボットテーブルってどうやるのかな、と。
ググるといろいろやり方が出てきますが、とりあえず、{dplyr}のgroup_byと{tidyr}のspreadを使うとできました。(もっと簡単な方法ありそうなので誰か教えて下さい。)


あと、実はその前段で、年齢のデータを「10代、20代、30代」みたいに階級別したかったので、そこには{berryFunctions}というパッケージのclassify関数を使ってみました。


まずパッケージを読み込みます。
{psych}の中に入っているTal_Orというデータセットを使います。
dplyrとtidyrは、tidyverseでもいいと思います。(両方まとめて入ります。)

> library(psych)  # データセットを使うため
> library(berryFunctions)  # classify用
> library(dplyr)  # tidyverseでまとめてもよい
> library(tidyr)  # tidyverseでまとめてもよい
> 


データセットを読み込んで中身をみてみます。

> data(Tal_Or)
> head(Tal_Or)
  cond pmi import reaction gender age
1    1 7.0      6     5.25      1  51
2    0 6.0      1     1.25      1  40
3    1 5.5      6     5.00      1  26
4    0 6.5      6     2.75      2  21
5    0 6.0      5     2.50      1  27
6    0 5.5      1     1.25      1  25
> 


これの、genderとageを軸につかって、reactionの値を集計したり、データの個数を数えたりしてみます。
その前に、年齢を、「10代、20代...」と階級に分けるために、classify関数をつかってみましょう。

> # 年齢階級ラベルを作成
> ageclass <- classify(Tal_Or$age, 
+                       method='custom', 
+                       breaks=c(10,19,29,39,49,59,69)
+                       )$index
> 


ヘルプをみると他の使い方も載っていますが、上の例では、methodには「自由に分割点を決める」という意味の'custom'を指定し、breaksに、分割点を与えています。
下記の参考例をみるとわかるように、

> # 参考:
> # 右端と左端の値は含まれる。その他の分割点は「以下」を表す。
> # なので10代、20代、30代、40代に分けたければ、両端は20と40
> # にしておいて、区切りは19、29、39を与える。
> classify(c(10,19,20,21,30,39,40), 
+          method='custom', 
+          breaks=c(10,19,29,39,49)
+          )$index
[1] 1 1 2 2 3 3 4
> 


breaksに与えた数字が両端と区切りになるので、ここでは数字を7つ与えているので6つの階級に区切られるのですが、注意点として、両端の値は「含む」であり、区切りは(未満ではなく)「以下」を表すようなので、「10代、20代...60代」とするときは、c(10,19,29,39,49,59,69)になるわけです。


以下、group_by、summarise、spreadを組み合わせて、ピボットテーブルを作ります。
簡単にいうと、軸にしたい2つの変数名でまずgroup_byする。
で、集計対象となる変数を合計したり平均とったりするのであれば、summariseの中で、変数名をつけて集計します。
最後に、spreadに、「列」にしたい軸変数の名前(そこで選ばなかったほうの軸変数が行になる)と、summariseの中で設けた集計変数名を指定します。
データの個数をカウントしたい場合(つまり2つの軸で度数のクロス集計をしたい場合)は、dplyrの中に入ってるn()という関数が使えます。

> # 年齢階級をくっつける
> Tal_Or <- data.frame(Tal_Or, ageclass)
> 
> # 性別と年齢階級でピボットテーブルを作成。
> # 値は"reaction"という変数にする。
> 
> # reactionの平均値を集計
> Tal_Or %>% 
+   dplyr::group_by(gender,ageclass) %>%
+   dplyr::summarise(mean_react = mean(reaction)) %>%
+   tidyr::spread(gender, mean_react)
# A tibble: 6 x 3
  ageclass   `1`   `2`
     <int> <dbl> <dbl>
1        1  2.5   4.92
2        2  3.53  3.44
3        3  4.25 NA   
4        4  1.25 NA   
5        5  5.25 NA   
6        6  2.38 NA   
> 
> # 行と列を入れ替え
> Tal_Or %>% 
+   dplyr::group_by(gender,ageclass) %>%
+   dplyr::summarise(mean_react = mean(reaction)) %>%
+   tidyr::spread(ageclass, mean_react)
# A tibble: 2 x 7
# Groups:   gender [2]
  gender   `1`   `2`   `3`   `4`   `5`   `6`
   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1      1  2.5   3.53  4.25  1.25  5.25  2.38
2      2  4.92  3.44 NA    NA    NA    NA   
> 
> # データの個数を集計
> Tal_Or %>% 
+   dplyr::group_by(gender,ageclass) %>%
+   dplyr::summarise(sample_size=n()) %>%
+   tidyr::spread(ageclass, sample_size)
# A tibble: 2 x 7
# Groups:   gender [2]
  gender   `1`   `2`   `3`   `4`   `5`   `6`
   <dbl> <int> <int> <int> <int> <int> <int>
1      1     1    37     1     1     1     2
2      2     3    77    NA    NA    NA    NA
> 


`1` `2` `3` `4` `5` `6`ってのは、さっき設けた年齢階級の階級番号ですね。一番下のテーブルでいうと、行が性別(gender)になっていて、1 or 2に分かれています。で、列が年齢階級になっていて、`1`が10代、`6`は60代を表しています。
つくったピボットテーブルに何か操作したいなら、下記のように、データフレームにします。

> # データフレームに
> Tal_Or %>% 
+   dplyr::group_by(gender,ageclass) %>%
+   dplyr::summarise(sample_size=n()) %>%
+   tidyr::spread(ageclass, sample_size) %>%
+   as.data.frame() -> pivot_table
> print(pivot_table)
  gender 1  2  3  4  5  6
1      1 1 37  1  1  1  2
2      2 3 77 NA NA NA NA
> 


私は、dplyrのクラスのままだとわけが分からなくなるので、as.data.frame()で標準のデータフレームに変換しています。

Rで動学的パネルデータ分析:plm、panelvarパッケージをつかったGMM推定

plmパッケージとpanelvarパッケージ

最近、パネルデータを扱うことが増えてきたのだが、パネルデータで動学的な(つまりt-1期とかのラグ項が出てくる)分析をやろうとすると最小二乗法ではなくGMM推定量を用いる必要がある。
備忘として、動学的パネルデータ分析(ダイナミックパネル分析)の基本的な考え方とRのパッケージの使いかたをここにメモしておこうと思う。といっても自分自身の理解もだいぶあやふやで、色々間違いもありそうなので、お読みになった方から指摘いただけると大変助かります…(汗)


Rの場合、結論から言うとまずは{plm}パッケージを使うのがいいと思う。いわゆる「パネルVAR」の形で分析したいなら、2018年に開発されたらしい{panelvar}パッケージを使うことができるのだが、後述するとおり{panelvar}のほうにはまだ不便なところもあって、今のところ、なるべく{plm}でできる分析をやったほうがいいと思う。まぁ、どちらも使い方自体は簡単である。なお、いずれにしても分析は線形モデルの範囲内になる。
{plm}パッケージについては日本語でもネット上に少しだけ解説があるが、たぶん{panelvar}のほうはまだ無いんじゃないかな?


なお理論的な背景については、幸い日本語でも、


千木良・早川・山本(2011)『動学的パネルデータ分析』
https://www.amazon.co.jp/dp/4862851029/


という非常に読みやすく実務的な利用にも配慮された教科書がある(神)。
また、クロスセクションと時系列の基本が頭に入ってる人なら、この奥井亮氏のスライドを読めばだいたいのストーリーは分かる(神)。
奥井亮:動学的パネルデータモデル


また、パッケージを使うときに、引数の設定や結果の解釈をするためにはヘルプ(plmpanelvar)を確認する必要があり、さらに意味が分からないところは作者の解説文書(plmpanelvar)もみることになると思われる。
しかし、とりあえずはさっきの奥井氏のスライドが大まかに理解できたら、あとは手を動かしながら並行して細かい解説を見ていったほうが早いのではないだろうか。

OLS(最小二乗法)の利用を妨げる諸問題

まず、社会調査データを扱う場合に発生する、OLSでは対処できない問題の種類、原因、対処法を大まかに知っておく必要がある。解説はネットで調べても山のように見つかるのでイメージだけまとめておく。


OLSで問題が生じる、ありがちなパターンとして、

  • 目的変数yの誤差項と説明変数xが相関している(内生性)
  • 目的変数yの誤差項どうしで相関(自己回帰の関係)がある(系列相関)
  • 目的変数yの誤差項の分散が一定じゃない(不均一分散)


といったケースがある。
内生性があるとOLSではパラメータの一致推定量が得られない*1。系列相関は、よほど大きくなければ一致推定は得られるとされるが、パラメータの標準誤差が正確に推定できず仮説検定(パラメータが有意かどうかの判断)が有効でなくなる。不均一分散の場合も、基本的に一致推定は得られるものの、仮説検定が有効でなくなる。


内生性は、「欠落変数」(交絡変数)や、「同時性」(逆の因果)や、セレクションバイアスや、測定誤差等によって引き起こされるもので、とくに、欠落変数による見せかけの相関や、同時性による逆の因果については、査読や研究発表などでも指摘されることは多いだろう。内生性の問題には、基本的には、操作変数(yの誤差項とは相関しないがxとは相関する変数)を用いた方法で対処することになる。有名なのは二段階最小二乗法(2SLS)*2で、それを一般化したものと言えるのがGMM(一般化モーメント法、一般化積率法)*3。とりあえず、「操作変数ってスゴい」ということを理屈として実感するには、まず二段階最小二乗法を理解するのがいいと思う。


系列相関は、欠落変数により引き起こされる場合もあるし、自己相関的な状態依存(惰性や反動など)によって生じる場合もある。分析後に検定して系列相関がみられる場合、可能ならモデルを変えたほうがいいとは思うが、誤差項の分散共分散構造を理論的に仮定できるのであれば、GLS(一般化最小二乗法)で対処する場合もある*4。パラメータの標準誤差が正確に推定できない問題については、ロバスト標準誤差*5を用いて対処する。


不均一分散は、たとえば、規模がかなり異なる会社や国のデータが混じっている場合や、あるいは急速に成長したりしている場合に、元の数字が大きければ分散も大きくなるといった具合で発生する。その場合、変数変換して、たとえば人口やGDPや売上規模で割るとか、あるいはべき乗分布している変数の対数を取ったりすることで、分散が均一に近づく場合がある*6。ただ、変数変換により情報が失われる面もあるから、そもそもモデルの関数型を変えるべきではないかと言われることもある。不均一の原因(たとえば人口規模の違いとか)が分かっているのであれば、その原因変数をGLSの一種である重み付き最小二乗法(WLS)のウェイトに使うことでも対処できる。標準誤差は、ロバスト標準誤差を使う必要がある。

パネルデータにおけるバイアスへの対処

ところで、パネルデータを用いると、欠落変数については、それが「個体によって異なるが時間に対して不変(iに依存するがtに依存しない)な変数」や「時間によって変化するが個体間を通じて共通(tに依存するがiに依存しない)な変数」であれば、固定効果として考慮することで消去できる。


固定効果の考慮とは、具体的には、個体ダミー・年次ダミーを入れたり、階差を取ったり、元データに対して平均除去(demeaning、縦断方向と横断方向の平均をそれぞれ引く。中心化 centeringともいう。)を行ったりした上でOLS推定すればよい。ただしパラメータの標準誤差についてはロバスト標準誤差を用いた上で、仮説検定をする必要がある。
ちなみに、時間と個体の両方に依存するなど、複雑な動きをする欠落変数の影響はこの方法では消去できないことを忘れてはいけない。たとえば、100地域×10期分のデータがあり、目的変数と説明変数が1つずつあって双方に天候が影響していると思われるのだが天候のデータが手元になく、かつ地域と期のどちらでみても天候が似るような観測単位ではない場合は、困ったことになる。*7


なお、固定効果モデルのほかにランダム効果モデルというのもある。あまりポピュラーではない気がするがこの記事で説明されていて、「内生性が考えられるが系列相関はないと言える」場合なら固定効果モデル、「内生性はないが系列相関や分散不均一性はあるかも知れない」場合ならランダム効果モデル*8を用いるとされている。


さて、このように、欠落変数バイアスにある程度対処できるというのが、パネルデータを用いる大きな利点であると言える。・・・が、しかしである。
上の奥井氏のスライドのp.37〜38あたりの説明が分かりやすいが、パネルデータモデルにラグ変数を導入して動学的な関係を組み込んだ場合、固定効果モデルと同様の方法では、固定効果を除くときに、別の内生性が新たに発生してしまう!!
たとえば一階階差を取る場合でいうと、ΔYtの誤差項に含まれているYt-1の誤差が、説明変数側のΔYt-1にも含まれているからである。


そこで、たとえばラグを1期分まで含むモデルなら、その変数からみて1期以上前の値を操作変数に用いることで、内生性に対処することになる。奥井氏のスライドのp.37でいうと、(y_it−1 − y_it−2)とε_it−1のあいだに相関があって困るという話なのだが、yの2期ラグであるy_it−2なら、(y_it−1 − yit−2)とは相関をもつけどε_it−1とは相関をもたないと考えられるから*9、これを操作変数に使うわけである。
で、具体的には回帰係数をGMM推定量として求めるのだが、細かいことは上記の千木良・早川・山本(2011)などに書いてあるからいいとして、イメージを言えば、目的変数と説明変数と操作変数の関係をいくつかのモーメント条件として記述した上で、ある最適化問題に帰着させてこれを解くと、いい感じにパラメータの一致推定量が得られるというわけである(雑過ぎ?)。


いい感じと言っても、注意点はある。

  • GMMは操作変数を大量に使う(使える)方法であり、「過剰識別制約」の検定(Sargenの検定やHansenのJ検定)を行って、要するに無駄に操作変数を取り入れ過ぎていないか(操作変数と撹乱項の直交が確保されてるか)を確認する必要がある。*10
  • 千木良・早川・山本(2011)でも注意喚起されているが、系列相関が無いことが前提になっているので、推定後に一応、系列相関がないかどうかArellano-Bond検定({plm}にはあるけど{panelvar}では提供されていない)等でチェックする。

 
 

GMMの種類

千木良・早川・山本(2011)ではGMMのバリエーションについて説明されている。まとめると・・・


一階階差GMM:最も一般的。後述のシステムGMMと違って、初期条件に仮定が必要ないところが最大の長所。短所として、撹乱項に系列相関があると一致性がなくなることが挙げられるが、系列相関は事後的に検定できるし、操作変数の工夫で回避できる(どういう工夫なのかは説明されてない)。それより、一番深刻な問題は、「弱操作変数」の問題である(弱操作変数の問題に対して頑健なのは後述のシステムGMM)。なお、1-step GMMと2-step GMMがあって、2ステップというのは1ステップ目の推定時の残差を用いて最適なウェイト行列を得てからもう1回推定するというやつなのだが、ふつうは2ステップのほうを使っておけばよい。*11


FOD-GMM:変数の変換の仕方が一階階差(first-difference)ではなくFOD(forward orthogonal deviation)になるだけで、利用可能な操作変数を全部使う場合は一階階差GMMと同じ推定値が得られる。操作変数を絞る場合は、少し結果に違いが生まれる。長所として、T≦N(クロスセクション方向の長さが時系列方向の長さを下回らない)であれば、FOD-GMMのほうが一階階差GMMよりも効率的な(分散が小さい)推定値を得られるらしい。


レベルGMM:モデルの説明変数から個別効果を取り除くのではなく、操作変数のほうから個別効果を取り除く方法。短所として、初期条件が平均定常でなければならない。しかし、レベルGMMはα(自己相関項の係数)が1に近づいても弱操作変数の問題が発生しにくいという長所がある。


システムGMM:一階階差GMMとレベルGMMをあわせて一つのシステムとみなして推定する(なお2ステップ推定が必要)。長所としては、αが1に近い場合でも弱操作変数の問題が生じないこと(αが1に近づくとレベルGMMのウェイトを増やすことで弱操作変数の問題を回避)。もうひとつは、一階階差GMMやレベルGMMよりも効率的であること。短所としては、レベルGMMのモーメント条件を含んでいることから、初期条件が平均定常でなければならない点。また、小標本ではバイアスが大きくなる可能性も指摘されている。


おそらく一階階差GMMかシステムGMMを使うのが基本なのだろうが、奥井氏の説明では、どういう場合にどっちが良いのかはいろいろ研究が進められてる段階で、議論が尽きないらしい。

{plm}パッケージの使い方

以上が前置きで、ここから具体的にRでの計算を実行していく。まずは{plm}からだ。
このパッケージの使いかたについては、以下の日本語資料で実行例が紹介されているので読むといいと思う。(ちなみにplmは「動学的」パネルデータ分析に限らず、静学的な固定効果モデルとかにも使えるもので、前半はその説明になっている。)
長倉大輔:R によるパネルデータモデルの推定


{plm}のヘルプの、pgmm()関数のところに載っている例を実行してみる。データは、{plm}に入っている練習用データのEmplUKを使う。

> data(EmplUK)
> head(EmplUK, 20)
   firm year sector    emp    wage capital   output
1     1 1977      7  5.041 13.1516  0.5894  95.7072
2     1 1978      7  5.600 12.3018  0.6318  97.3569
3     1 1979      7  5.015 12.8395  0.6771  99.6083
4     1 1980      7  4.715 13.8039  0.6171 100.5501
5     1 1981      7  4.093 14.2897  0.5076  99.5581
6     1 1982      7  3.166 14.8681  0.4229  98.6151
7     1 1983      7  2.936 13.7784  0.3920 100.0301
8     2 1977      7 71.319 14.7909 16.9363  95.7072
9     2 1978      7 70.643 14.1036 17.2422  97.3569
10    2 1979      7 70.918 14.9534 17.5413  99.6083
11    2 1980      7 72.031 15.4910 17.6574 100.5501
12    2 1981      7 73.689 16.1969 16.7133  99.5581
13    2 1982      7 72.419 16.1314 16.2469  98.6151
14    2 1983      7 68.518 16.3051 17.3696 100.0301
15    3 1977      7 19.156 22.6920  7.0975  95.7072
16    3 1978      7 19.440 20.6938  6.9469  97.3569
17    3 1979      7 19.900 21.2048  6.8565  99.6083
18    3 1980      7 20.240 22.1970  6.6547 100.5501
19    3 1981      7 19.570 24.8714  6.2136  99.5581
20    3 1982      7 18.125 24.8447  5.7146  98.6151

これは、ArellanoとBondがGMMによる動学パネルの分析方法を提案した1990年の論文で使っているデータで、firmが企業を表しているのだが、企業番号1〜140、年次は1977年から1983年までのデータがあるというものだ。


データフレームのレイアウトについて、{plm}パッケージでは個体やグループを表す変数は数字にしておく必要がある点に注意が必要だ。で、基本的には個体変数を1列目、時間変数を2列目に書くのだが、列が別の場所にあったとしても後述のindex引数で明示的に指定すればOK。


以下では、Arellano & Bond(1990)のTABLE4の(c)列と同じ結果が得られる分析をしてみよう。元論文では変数の対数を取っていて、上の長倉氏の実行例のようにformulaの中でlog()を挟めばいいのだが、コードが見づらくなるので先に全部対数を取っておく。

# 先に対数をとっておく
for(i in 4:7){
  EmplUK[,i] <- log(EmplUK[,i])
}


そしてパネルモデルの推定を実行!

> # パネルGMMを実行
> # Arellano & Bond(1990)のTABLE4の(c)列と同じ結果を出す
> 
> z1 <- pgmm(emp ~ lag(emp, 1:2)       # 自己相関項(ラグ1〜2)
+                + lag(wage, 0:1)      # ラグ付き説明変数
+                + capital             # ラグ無し説明変数
+                + lag(output, 0:1)    # ラグ付き説明変数
+                | lag(emp, 2:99),     # 操作変数(最大限使う)
+                data = EmplUK,        # データフレームを指定
+                index=c('firm','year'), # 個体と時間変数の指定
+                effect = "twoways",   # 固定効果の設定
+                model = "twosteps",   # 2-stepGMMならtwostepとする
+                transformation = 'd'  # 'ld'ならシステムGMM、'd'なら一階階差GMM
+            )


ラグを含める変数はlag(変数名,期)という書式で表現する。
期のところを0にするとラグではなく当期の値になる。つまり、0:2なら当期から2期前まで、1:2なら1期前から2期前までの変数がモデルに入るという意味である。もちろん「0」とだけ書いたらラグなし変数になる。
注意点として、{dplyr}パッケージの中にlagという関数があり、{dplyr}を読み込んでると変な誤作動が起きて「`n` must be a nonnegative integer scalar, not an integer vector of length...」みたいなメッセージが出たりする。(参考
なので、dplyrの利用を避けるか、前処理で使ってしまっている場合は(というかパネルデータの前処理でdplyrのgroup_byとかが便利なので使うことが多いと思う*12)、plmを用いる直前で下記のように封印しておかなければならない

detach(package:dplyr)


で、モデルの式を書いたあとに|を挟んで操作変数を定義する。これも書き方は同じで、上記ではlag(emp, 2:99)と書いているので、emp(雇用)の2期目から99期目までを操作変数に使うと指定していることになる。
ただここの99には意味がなく、これは要するに「使える操作変数は全部使う」みたいな意味合いである(もちろんTが100期以上あるならもっと大きくする。使えない期数のものは自動的に無視される)。ヘルプの書き方がそうなっているので、それにあわせている。操作変数の数を減らすときは、明示的に少ない数を与える。


操作変数は、右辺にある内生変数がラグ1期なら、ラグ2期以前の値を使う必要がある。内生変数がラグ1〜3期まであるなら、少なくとも操作変数としてはラグ2期〜4期(双方の個数が一致する時「丁度識別」という)のものは必ず含めておく必要がある。先程も言ったように、操作変数の数をどこまで入れるかというのは難しい問題で、分析例をみてると使える最大期数まで使っているものは多いのだが、後述する過剰識別制約の検定結果をみてアウトだったら減らす必要がある。


ところで、empの自己相関項以外の説明変数が上記では全部外生変数として考慮されてるわけだが、それらに内生性があることが懸念される場合は、それらに対する操作変数も与えればよい。具体的には、操作変数のところを、

| lag(emp, 2:99) + lag(wage, 2:99)

とかにすればよく、ヘルプにもそういう使用例が載っている。*13


indexは、何も指定しなければ、データフレームの1列目が個体変数、2列目が時間変数と解釈される。明示的に指定するときは、列名称又は列番号のベクトルで、c(個体変数,時間変数)の順に与える。


effectのところは、要するに時間効果(タイムダミー)を入れるかどうかの設定である。個体効果(ここでは企業=firmの個別効果)はもともと入っててそれが一階階差を取るときに消えるという話なんだが、ここを'twoways'にすると、時間効果もタイムダミーの形で考慮されるようになる。
modelのところは2-step GMMでやるのが普通なので、'twosteps'としておけばよいだろう。
transformationは、'd'にしておくと単なる一階階差GMM、'ld'にするとシステムGMMになる。
以下、推定結果。

> summary(z1, robust = F, time.dummies = F)
Twoways effects Two steps model

Call:
pgmm(formula = emp ~ lag(emp, 1:2) + lag(wage, 0:1) + capital + 
    lag(output, 0:1) | lag(emp, 2:99), data = EmplUK, effect = "twoways", 
    model = "twosteps", transformation = "d")

Unbalanced Panel: n = 140, T = 7-9, N = 1031

Number of Observations Used: 611

Residuals:
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-0.6190677 -0.0255683  0.0000000 -0.0001339  0.0332013  0.6410272 

Coefficients:
                   Estimate Std. Error  z-value              Pr(>|z|)    
lag(emp, 1:2)1     0.474151   0.085303   5.5584    0.0000000272221654 ***
lag(emp, 1:2)2    -0.052967   0.027284  -1.9413             0.0522200 .  
lag(wage, 0:1)0   -0.513205   0.049345 -10.4003 < 0.00000000000000022 ***
lag(wage, 0:1)1    0.224640   0.080063   2.8058             0.0050192 ** 
capital            0.292723   0.039463   7.4177    0.0000000000001191 ***
lag(output, 0:1)0  0.609775   0.108524   5.6188    0.0000000192269973 ***
lag(output, 0:1)1 -0.446373   0.124815  -3.5763             0.0003485 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Sargan test: chisq(25) = 30.11247 (p-value = 0.22011)
Autocorrelation test (1): normal = -2.427829 (p-value = 0.01519)
Autocorrelation test (2): normal = -0.3325401 (p-value = 0.73948)
Wald test for coefficients: chisq(7) = 371.9877 (p-value = < 0.000000000000000222)
Wald test for time dummies: chisq(6) = 26.9045 (p-value = 0.0001509)


coefficientsのところは分かりやすいと思うが、他に注目するところとしては、下の各種検定結果だ。
とくに、Sargen testのところは過剰識別制約の検定になっており、ここのp-valueが有意(例えば<.05)だと、操作変数入れすぎ。
またAutocorrelation testは系列相関のテストである。一階階差を取っていることから、誤差項に1次の自己相関があるのは自然なことなので、(1)のほうはp-valueが有意でも問題ない。(2)が有意(例えばp<.05)だと、系列相関があることになるので、その場合はモデルを変えたりする必要がある。


ところで、上記の論文ではロバスト標準誤差を使ってなくて、かつタイムダミーの係数の推定結果を載せてないのだが、それらを使う場合はsummary()の中で指定すればいいで。上のコードがそうなってるが、ここでタイムダミーを表示するかどうかも設定できる(表示するかどうかが違うだけで、計算自体は、effect = "twoways"としてあるならタイムダミーありで推定される)。
なおこれらは、pgmm()の引数としても設定できるのだが、summary()と矛盾する場合にsummary()側の設定が優先されるみたいなので、summary()で指定することにしとけばいいと思う。

> # ロバスト標準誤差を用い、かつタイムダミーのパラメータも表示
> summary(z1, robust = T, time.dummies = T)
Twoways effects Two steps model

Call:
pgmm(formula = emp ~ lag(emp, 1:2) + lag(wage, 0:1) + capital + 
    lag(output, 0:1) | lag(emp, 2:99), data = EmplUK, effect = "twoways", 
    model = "twosteps", transformation = "d")

Unbalanced Panel: n = 140, T = 7-9, N = 1031

Number of Observations Used: 611

Residuals:
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-0.6190677 -0.0255683  0.0000000 -0.0001339  0.0332013  0.6410272 

Coefficients:
                    Estimate Std. Error z-value    Pr(>|z|)    
lag(emp, 1:2)1     0.4741506  0.1853985  2.5575   0.0105437 *  
lag(emp, 1:2)2    -0.0529675  0.0517491 -1.0235   0.3060506    
lag(wage, 0:1)0   -0.5132048  0.1455653 -3.5256   0.0004225 ***
lag(wage, 0:1)1    0.2246398  0.1419495  1.5825   0.1135279    
capital            0.2927231  0.0626271  4.6741 0.000002953 ***
lag(output, 0:1)0  0.6097748  0.1562625  3.9022 0.000095304 ***
lag(output, 0:1)1 -0.4463726  0.2173020 -2.0542   0.0399605 *  
1979               0.0105090  0.0099019  1.0613   0.2885484    
1980               0.0246512  0.0157698  1.5632   0.1180087    
1981              -0.0158019  0.0267313 -0.5911   0.5544275    
1982              -0.0374420  0.0299934 -1.2483   0.2119056    
1983              -0.0392888  0.0346649 -1.1334   0.2570509    
1984              -0.0495094  0.0348578 -1.4203   0.1555141    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Sargan test: chisq(25) = 30.11247 (p-value = 0.22011)
Autocorrelation test (1): normal = -1.53845 (p-value = 0.12394)
Autocorrelation test (2): normal = -0.2796829 (p-value = 0.77972)
Wald test for coefficients: chisq(7) = 142.0353 (p-value = < 0.000000000000000222)
Wald test for time dummies: chisq(6) = 16.97046 (p-value = 0.0093924)
> 

 
 

{panelvar}パッケージの使い方

さて、こんどは{panelvar}パッケージを使ってみよう。
これはその名の通り、パネルデータでVARモデルを組むのに使えるもので、ふつうの時系列データのVARモデルのように、インパルス応答関数を出す機能も提供されている。
Stataのxtabond2をRに移植する目的で開発されたものらしいので、あとで暇ができたら、Stataの出力とも比較しておこうかなと。


さて、上記の{plm}と同じデータをつかって推定してみよう。変数を全部内生変数にして、4変数のVARにしてみる。

pvar.1<- pvargmm(dependent_vars = c('emp','wage','capital','output'),
                 lags = 2,                       # 内生変数のラグ次数
                 transformation = "fod",         # fdなら一階階差、fodならFOD変換
                 #predet_vars = c('',''),        # 先決内生変数
                 #exog_vars = c('',''),          # 外生変数
                 data = EmplUK,                   # 分析するデータフレーム
                 panel_identifier=c('firm','year'),  # 個体&時間変数の指定
                 steps = c("twostep"),           # 1-stepか2-stepか
                 system_instruments = T,         # TRUEだとシステムGMM
                 max_instr_dependent_vars = 3,   # 内生変数の操作変数の最大ラグ
                 max_instr_predet_vars = 3,      # 先決変数の操作変数の最大ラグ
                 min_instr_dependent_vars = 2L,  # 内生変数の操作変数の最小ラグ
                 min_instr_predet_vars = 2L      # 内生変数の操作変数の最小ラグ
                 )

summ.1 <- summary(pvar.1)


dependent_varsのところには、内生変数を与える。列名称か、列番号。ここにあたえた変数でVARモデルが構成される。lagは内生変数のラグ次数。


transformationが'fd'なら一階階差変換でのGMM、'fod'ならFOD変換をしたGMMになる。


上の例では空っぽにしてあるのだが、predet_varsのところには先決変数を指定できる。先決変数というのは、文字通り先に決まるという意味で、yの誤差項と相関はしてないと考えられるが、yの誤差項が将来のxとは相関すると考えられるようなxのこと。その場合、一階階差GMMの処理でラグを取ったときに内生性が生まれてしまうのでこれらの変数にも操作変数を与える必要があるわけである。
exog_varsはそのまんま、外生変数を与える。つまり外生変数付きVARにできるわけだ。


panel_identifierは、{plm}のindexと同じ。
stepsは'twostep'で良いだろう。(plmと違って単数形で書くみたい)
system_instrumentsをT(TRUE)にしてるとシステムGMMになる。Fにすれば単なる一階階差GMMとかになる。


maxとminがついているところは、操作変数を何期前から何期前まで使うかを指定する引数である。
minを2にしてmaxを4にしていると、t-2からt-4までの値が操作変数に使われる。
ここもさっきと同様、ヘルプの使用例では最大限使うという意味でmaxは99にされていた。
minを「2L」というふうにLをつけてるのは、明示的に整数値を与えるという意味なのだが、なんでそうする必要あるのかはわからない。ただ、ヘルプでそうなってたのでそうしている。


推定結果をみてみる。

> summary(pvar.1)
---------------------------------------------------
Dynamic Panel VAR estimation, two-step GMM 
---------------------------------------------------
Transformation: Forward orthogonal deviations 
Group variable: firm 
Time variable: year 
Number of observations = 611 
Number of groups = 140 
Obs per group: min = 4 
               avg = 4.364286 
               max = 6 
Number of instruments = 308 

================================================================
              emp          wage         capital      output     
----------------------------------------------------------------
lag1_emp       1.0963 ***   0.0690       0.1847      -0.0245    
              (0.0963)     (0.0707)     (0.1019)     (0.0326)   
lag1_wage      0.1564       0.7162 ***  -0.1312      -0.0792    
              (0.1431)     (0.0873)     (0.1271)     (0.0478)   
lag1_capital   0.1722 **   -0.0056       1.1135 ***   0.0206    
              (0.0563)     (0.0375)     (0.0508)     (0.0214)   
lag1_output    0.5364 ***  -0.6185 ***   0.8483 ***   1.1223 ***
              (0.1191)     (0.0932)     (0.1548)     (0.0550)   
lag2_emp      -0.0960      -0.0434      -0.1164       0.0437    
              (0.0721)     (0.0538)     (0.0906)     (0.0330)   
lag2_wage     -0.0597       0.0874       0.1626       0.1107 ** 
              (0.1013)     (0.0640)     (0.1013)     (0.0367)   
lag2_capital  -0.1623 ***  -0.0153      -0.1574 **   -0.0329    
              (0.0477)     (0.0369)     (0.0500)     (0.0307)   
lag2_output   -0.6652 ***   0.4143 ***  -0.8420 ***  -0.5837 ***
              (0.1564)     (0.1232)     (0.1860)     (0.0717)   
const          0.2682       1.5356 **   -0.2375       2.0043 ***
              (0.7296)     (0.5465)     (0.6384)     (0.2880)   
================================================================
*** p < 0.001, ** p < 0.01, * p < 0.05

---------------------------------------------------
Instruments for  equation
 Standard
  
 GMM-type
  Dependent vars: L(2, 3)
  Collapse =  FALSE 
---------------------------------------------------

Hansen test of overid. restrictions: chi2(272) = 132.37 Prob > chi2 = 1
(Robust, but weakened by many instruments.)
> 


一番下に、Hansenの過剰識別制約検定の結果が載っている。一番右の「= 1」というところが要するにp値で、これが有意でなければ過剰識別ではないのでOKということだが、1となってるので問題ないのだろう。
ちなみに、以下のようにそれだけを求める関数もある。

> hansen_j_test(pvar.1)
$statistic
[1] 132.3731

$p.value
[1] 1

$parameter
[1] 272

$nof_instruments
[1] 308

$method
[1] "Hansen-J-Test"


ところで、lagの次数を何期にすべきかの選びかたが難しい({vars}パッケージのように選択するための機能が提供されていない)。推定後のpvarモデルをAndrews_Lu_MMSC()関数に与えるとAIC等が得られるので、ラグを変えながらちまちまAICを確認していけば良いのかも知れない。ちなみにpvargmmの1回の推定にかなり時間がかかるのでとても大変・・・。

> # HQは予めパラメータの設定がいるが2.1にしておけばいいと思う(dn=2loglogn)
> Andrews_Lu_MMSC(pvar.1, HQ_criterion = 2.1)  
$MMSC_BIC
[1] -1785.741

$MMSC_AIC
[1] -465.6269

$MMSC_HQIC
[1] -1034.676


インパルス応答関数は以下のようにして求められる。

# インパルス応答関数(結果はfrom=impulse側の変数ごとにまとめられている)
irf.1 <- oirf(pvar.1, n.ahead=10)
plot(irf.1)


f:id:midnightseminar:20200103223845p:plain


なお、irf.1の中身は見づらいのだが、リストがimpulse側の変数ごとにまとめられていると思って見ていけば分かると思う。
インパルス応答関数の信頼区間も出すことができるのだが、これは極めて危険で、えげつない時間がかかる。

boot.1 <- bootstrap_irf(pvar.1, typeof_irf = c("OIRF"),
                        n.ahead=10, 
                        nof_Nstar_draws=100, 
                        confidence.band = 0.95)


typeof_irfのところは、OIRFにしておくと直交化インパルス応答、GIRFにすると一般化インパルス応答関数になる。
nof_Nstar_drawsはブートストラップのサンプリングの回数。これ、{vars}パッケージで求める(パネルじゃないデータの)IRFの場合、デフォルトが確か100回になってて一瞬で終わるのだが、この{panelvar}のブートストラップは1回分やるのに私のパソコンだと数分かかる場合もあるので、とてもじゃないが気軽には回せない。


なお、内生変数を1個だけにして他を外生変数や先決変数にすればVARではなくなり、{plm}でやっているのと同じような分析になるのだが、計算の速さとかからして{plm}のほうがいいと思う。また、{panelvar}では、系列相関の検定結果が提供されないので、自分で理論を確認しながら組む必要があり、そんなことするぐらいならやはり素直に{plm}を使っておけばいいのではないだろうか。
もちろん、「パネルVAR」にしたいケースもあるとは思われ(ちゃんと読んでないが、ECBがパネルVARをつかった研究のサーベイをしてた)、その場合は{panelvar}を使うことになるけど、いろいろ調べても解説がまだ少なく、私のようなこの種の分析に慣れてない人間には使いづらそうではある。

パネル単位根検定

パネルデータも、パネルじゃない時系列の場合と同様に、単位根や共和分の問題に気をつける必要がある(といっても今のところ、どういう場合にどこまで問題になるのか自分でイマイチ頭が整理できていない)。
千木良・早川・山本(2011)では、80ページぐらいを費やして様々なタイプのパネル単位根検定・定常性検定が紹介されていて、その後にパネル共和分についても解説されている。
パネル単位根検定としてはシンプルかつ古いのがLevin, Lin & Chu(LLC)の検定で、これはADF検定をパネルに拡張したもの。これに加えて個体間の差異を認める方向で拡張したのがIm, Pesaran & Shin (IPS)の検定。これをさらに画期的な方法で拡張したと言われるのがcombination検定というやつで、Maddala & Wu (1999)とChoi (2001)によって別々に提案された。
以上のものはすべてADF検定の拡張であり、帰無仮説が単位根で対立仮説が定常という「単位根検定」である。逆に、帰無仮説が定常で対立仮説が単位根である「定常性検定」というのもあって、これは時系列分析におけるKPSS検定ってやつをHardi(2000)がパネルデータに拡張したものがある。
以上は「第1世代のパネル単位根検定」と呼ばれる手法群で、{plm}ではこれらの機能が提供されている。下記のpurtest()関数で、引数testに、LLC検定は'levinlin'、IPS検定は'ips'と設定する。combination検定については、Maddalaらの検定(これはChoiの提案したいくつかの手法のうちp検定と呼ばれるもの)は'madwu'と設定する。Choiが他にいくつか提案してる方法も実装されていて、Nが大きい(クロスセクション方向が長い)場合のp検定は'Pm'、あとはChoi(2001)を読んでないのでよくわからないがinverse normal testというのが'invnormal'、logit testというのが'logit'という設定で使える。Hardiの定常性検定はそのまんま'hardi'とする。
Hardiだけは毛色が違っていて、これは「すべての個体において定常」というのが帰無仮説なので、p値が有意なら「少なくとも1つが単位根」となる。他は、p値が有意であれば、単位根が棄却される。

> # パネル単位根検定
> purtest(object= EmplUK, 
+         x     = 'emp',
+         data  = EmplUK,
+         index = c('firm','year'),  # 個体、時間の順で入れる
+         test  = 'madwu', # 'levinlin','ips','madwu','Pm','invnormal','logit'or'hadri'
+         exo   = 'trend',    # 'trend','none'or'intercept'
+         lags  = 'AIC',      # ラグはAICで適当に選んでという意味
+         pmax  = 6,          # 考慮するラグの最大次数
+         Hcons = TRUE        # test="hardi"のときだけ関係ある設定
+ )

	Maddala-Wu Unit-Root Test (ex. var.:
	Individual Intercepts and Trend)

data:  EmplUK
chisq = 597.8, df = 14, p-value < 2.2e-16
alternative hypothesis: stationarity


時系列分析のときと同じように、単位根検定にドリフト項(定数項)やトレンド項を入れるかどうかというのは大きな問題となる。上記では、exoのところを'none'にするとどっちも無し、'intercept'にするとドリフト項あり、'trend'にするとドリフト項もトレンド項もありとなる。
ところで大きな問題点として、パネル単位根検定は「すべてが単位根過程である」という帰無仮説を検定しているので、逆にいうと1つでも単位根仮定でない個体があるなら、「定常である」という対立仮説が採用されることになってしまう。これは、単位根の検出力が非常に低いことは意味していて、パネル単位根検定で単位根が棄却できるからといって、個々の個体の時系列が定常であるわけではない。しかし、千木良・早川・山本(2011)が説明してるのだが、パネルデータ分析というのはそもそも各グループ・各年に共通のパラメータを求めに行く分析なのだから、これは本質的に避けがたい問題であると考えたほうがよいようである。


以上はすべて「第一世代の単位根・定常性検定」と言われるもので、クロスセクション方向には基本的に独立性を仮定しているのだが、「第二世代」の検定手法群になるとクロスセクション方向の相関も認めるという拡張が行われる。{plm}パッケージでは第2世代のうちCross-sectionally augmented IPS testというのをcipstest()という関数で実行できる。帰無仮説は「単位根過程あり」。

> # pseriesというクラスのデータを与えないと行けないので…
> Produc.p <- pdata.frame(Produc, index=c('state', 'year'))
> 
> cipstest(Produc.p$gsp, 
+          lags = 1,
+          type = 'trend',
+          model = 'cmg', 
+          truncated = FALSE)

	Pesaran's CIPS test for unit roots

data:  Produc.p$gsp
CIPS test = -0.70295, lag order = 1, p-value
= 0.1
alternative hypothesis: Stationarity

Warning message:
In cipstest(Produc.p$gsp, lags = 1, type = "trend", model = "cmg",  :
  p-value greater than printed p-value
> 


帰無仮説が棄却されてないので、単位根ありってことだろう。

その他

上の例で書かなかったが、固定効果が含まれているかどうかの検定をHausman検定を使って行うことができるがいまいちロジックを理解してないので後ほど確認しておきます。
千木良・早川・山本(2011)の共和分検定の章をまだ読んでないのでこれもあとで確認します。
また、(一階階差GMMの時に問題になる)初期条件の平均定常性の検定なども行うことが望ましいとされます。

パッケージの比較

現時点での、{plm}と{panelvar}の比較をまとめておきます。

  • パネル単位根検定は、plmでのみ提供されている(第2世代までカバーされている)。
  • Hausman検定も、plmでのみ提供されている。
  • パネル「VAR」モデルの分析は、plmでは実行できずpanelvarが必要。
  • 系列相関の検定が、plmでは提供されているが、panelvarでは提供されていない。
  • 過剰識別制約の検定はどちらでもできる。
  • モデル全体のAICを確認する手段が、panelvarでは提供されているが、plmでは提供されていない。(Stackoverflowでそれを試みているエントリはあったが詳細確認してません。)
  • 計算速度は、plmは気になるほどの時間はかからないが、panelvarは遅い。

*1:内生性があるかどうかはまず理論的に考えるべきだが、結論を歪めるほどのものであるか否かは、HausmanやWooldridgeの方法による内生性検定で確認できる。 http://www.ier.hit-u.ac.jp/~kitamura/lecture/Hit/08Statsys5.pdf

*2:GLS(誤差の分散共分散行列に制約を仮定する「一般化最小二乗法」)を組み合わせた三段階最小二乗法(3SLS)というものもある。

*3:完全情最尤推定法というのもあるらしいがよく分かってない。

*4:コクラン=オーカット法や、プレイス=ウィンステン法。

*5:なんか種類は色々ある。

*6:対数差分を取るとそれは変化率に相当するが、経済データだと「前年比」に着目するとうまく捉えられる変数も多いので、理論的にも都合がいいとされる

*7:個別のケースごとになら、色々細かい対処法が考案されてそうな気はするが

*8:ランダム切片モデルというやつで、この場合、一種のGLSとみなすことができる。

*9:厳密に言えば、Yt-1の誤差項とYt-2が絶対に相関しないのか、つまり新しい期の値の誤差が古い期の値と相関を持たないのかというと、その変数の自己相関がうまく定式化されてなくて、古い期の値が新しい期の値に十分反映されてない場合、相関が残るとは思う。

*10:遠い過去の値はそもそも「弱操作変数」の問題を引き起こすので、千木良他(2011)では実務上は「3期まで」に限る方針がとられることが多いとしているが、どれぐらい一般的なのかは不明。

*11:ちなみに1ステップGMMはいわゆる二段階最小二乗法と同値。

*12:パネル用データフレームであるpdata.frameクラスの挙動を理解すればdplyrを使わなくても便利に計算できるのかも知れないが、確認していない。

*13:wageのほうを1:99にしたほうがいい場合もある気がするが、理論的にどういう内生性を考えてるかによるかな?ヘルプで引用されてる Blundell and Bond (1998)をあとで見て考えよう。

Rでよく忘れる、よく間違える書き方(随時追記)

よく忘れることのメモです。

  • NAかどうかの判定にはx==NAとかではなくis.na(x)を使う
  • 要素に含まれるかどうかの判定は、%in%かis.element()を使う。これはデータ全体の中から何かを抽出するときの条件を複数条件にしたい場面でも使えるときがある。(たとえば簡単なところでは、subset()でデータフレームから条件付の抽出をするとき、&で条件を並べなくてもよくなったりとか)
> is.element(3, c(1,2,3))
[1] TRUE
> 3 %in% c(1,2,3)
[1] TRUE
  • {dplyr}でデータフレームを操作した後、as.data.frame()で標準のデータフレームに戻してからじゃないと、別の関数に与えたときにエラーが出る場合がある
  • {data.table}のfread()関数でcsvとかを読み込んだ場合も、それを標準のデータフレームにしてから使わないとエラーが出ることがある
  • factor型のデータをnumericに変えると、factorの中身ではなくレベル(水準)番号が数字になる。
> x <- as.factor(c(1,10,100,1,10,100))
> as.numeric(x)
[1] 1 2 3 1 2 3
  • リストに要素をアペンドするとき、c()関数でできるのだが、付け加える方にlist()をかけとかないといけない
> x1 <- c(1,2,3)
> x2 <- c('a','b','c')
> l1 <- list(x1, x2)
> print(X)
[[1]]
[1] 1 2 3

[[2]]
[1] "a" "b" "c"

[[3]]
[1] "あ"

[[4]]
[1] "い"

[[5]]
[1] "う"

> x3 <- c('あ','い','う')
> l2 <- c(l1, x3)  # これはダメ
> print(l2)
[[1]]
[1] 1 2 3

[[2]]
[1] "a" "b" "c"

[[3]]
[1] "あ"

[[4]]
[1] "い"

[[5]]
[1] "う"

> l3 <- c(l1, list(x3))  # こうする
> print(l3)
[[1]]
[1] 1 2 3

[[2]]
[1] "a" "b" "c"

[[3]]
[1] "あ" "い" "う"
  • パネルデータをロング型にしたりワイド型にしたりするのには{tidyr}を使う。ワイドをロングにするのはgather()関数で、ロングをワイドにするのはspread()関数で。
  • ||とか&&は、or/and条件を入れ子にするときに使う。
  • ベクトルの要素が「全て◯◯という値である」という条件を書くシンプルな方法はたぶんない?最大と最小が一致するみたいな書き方をいつもしているのだが。
  • データを「中心化」(平均ゼロにする)したいときは、center()じゃなくて、scale()関数のオプションでscale=Fに設定する。変な話だが。
  • apply()の第一引数は小文字のxではなく大文字のX。まあ、書かなければいいともいえるが、書くときによく間違える。
  • logicalのベクトルでTRUEの個数を数えるには、length(x[x==TRUE])とかsum(x==TRUE)で数えられる。
  • 欠損値のない行にしぼりたいときは、d[complete.cases(d),]でOK。
  • plotで散布図の記号を文字(国の名前とか人の名前とか)にしたいなら、散布図の方は色を"white"にして消して、text(x, y, labels=)とすればよい。
  • 離散変数の値ごとのヒストグラムを得たいとき、hist(x)としてしまうと、離散値に対応するバーが目盛りの左側に来てしまう。これを回避するには、barplot(table(x))とすればよい。
  • 回帰係数の信頼区間の出し方をいつも忘れるが、confint(model, level=0.95)
  • 主成分分析の結果は符号が想定とは逆になってることがあり、必要ならマイナスをかけて戻すのを忘れないように
  • applyでFUNに第二引数以降をわたしたいときは、FUNの後ろにそのまま続けていけばよい。
d <- data.frame(
  A = c(1,2,3,4,5),
  B = c(2,4,NA,8,10),
  C = c(1,3,5,7,9)
  )
apply(d, MARGIN=1, FUN=mean, na.rm=T)
  • dplyrでパイプライン%>%を使うとき、関数定義から埋め込みたい場合は、%>% (functon(p){hoge}) %>%と、()でかこむ。
  • 意外と検索しても見つからないのだが、ggplot2で凡例が小さすぎて点線とか意味わからないのは、theme(legend.key.width=unit(2,"cm"))として解決。
  • dplyrのmutateで、複数の列に一気に同じ処理をしたいときは、mutate_at(var(hoge), funs(f(.)))とし、hogeのところにはstart_withとかが使えるし、ある変数以外を指定したいなら-をつければよい。
  • たまにRでエクセルのファイルを読みたい時は、今なら、gdataパッケージがいいのではないかと思う。
library(gdata)
df <- read.xls('ファイル名')
  • dplyrのなかで差分を取りたいときは、diffじゃなくてlagをつかって、mutate(x = x - lag(x))とする。
  • 正規表現でエスケープするとき、バックスラッシュ\は2回重ねないといけない。
  • dplyrで、少なくともselect、lead、lagは毎回"dplyr::"つけたほうがいい。
  • 空のデータフレームをつくるとき、matrixをas.data.frame()するのがよいが、これに1行目のデータを与えるときにrbindをつかってしまうと列名称が変わってしまうので、インデックスを指定して与えるようにする必要がある。
as.data.frame(matrix(c(NA,NA),nrow=1))
  • read.csv()するとき、stringsAsFactors = Fオプションをとにかく付けるように習慣づけないと、おもわぬエラーが起きる。
  • データフレームにscale()を適用すると、各列を標準化してくれるが、返り値がmatrixになってしまうので、たとえばlm関数のdata欄に直接与えることができない。すこしめんどうだが、numeric型の列を抽出し標準化してデータフレームに直す必要があるので、以下のような感じになる。

lm(y~x, data=d %>% select_if(is.numeric) %>% scale() %>% as.data.frame())

  • RStudioでggplot2を使っている時に、
Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto),  : 
  Viewport has zero dimension(s)

というようなエラーがでるときがあるが、コードが間違っているのではなく、plotが出るペインの掃除ボタン(ほうきのマーク)を押すと解消したりする。

  • dplyrでcomplete.cases()をやりたいときは、filter(complete.cases(.))とする。
  • AとBに挟まれた部分を取る正規表現
"(?<=A)(.*)(?=B)"
  • warningを消したい時、単にinvisible()をかければいいのではなくて、以下のようにcapture.output()を入れ子にする。
invisible(capture.output())
  • ggplot2で"Error in .Call.graphics(C_palette2, .Call(C_palette2, NULL)) : invalid graphics state"などよくわからないエラーがでたときは、とりあえずRStudioのplot欄を掃除する(clear all plots)
  • 特定のオブジェクトだけ残して自作のオブジェクトを全部消したいときは、
rm(list=ls()[which(ls() != 'hogehoge')])
  • forループの進捗を表示させたいとき、以下のような感じにすればよい。(1000周ごとに表示)
 if(i %% 1000 == 0) {message('\r', paste(i, ' / ', length(xxx)), appendLF=FALSE)}
  • glmでたとえばロジスティック回帰をやったとき、予測値は0,1ではなく0-1の少数の形で以下のコードで得られる。typeのところはデフォルトでは'link'になるが、これは線形予測子を得るもの。type='response'にしておくと、リンク関数を経たものが得られる。
predict(model, type='response')
predict(model, type='response')
  • lmとかglmオブジェクトから、「投入したデータフレーム」は$modelで取り出せる。stepしている場合は、最終的に選択されたモデルのデータフレームが取り出せる。で、ポイントとしては、目的変数が必ず1列目に来るようになっていることと、NAはオミットされてcomplete.casesになっているということ。predict関数を使うときに(newdataに与えるので)覚えておくと便利。
  • 重み付き最小二乗法をやるときはlm(weights=)という引数をつけるが、predictで予測値を出すときは、predictの中にweightsとか書かなくてもよい。lmオブジェクトの$model内の最後に(weights)っていう列ができててそれが勝手に考慮される。predict内にweightsって書いても無視されてるっぽい(変な値を与えても影響がない)。
  • 自作関数の引数でデータフレームの列を指定したい場合に、関数内でデータフレームをattach()して以下のようにすればよいという意見をツイッターでみた。私は、変数名はdata[,'var']と文字列型で指定するようにしているが。
get_hist <- function(data, var) {
  attach(data)
  hist(var)
  detach(data)
}
  • 文字列のベクトルを連結して1つの文字列にしたいときは、sep=ではなくcollapse = を使う。そうしないと、ベクトルを与えると要素を別々に扱ってベクトルを返そうとする。
  • 行列やベクトルの掛け算をするとき、*を使うと、線形代数的な積ではなく、要素同士の掛け算になる。挙動は以下のとおり。
> # 2×6の行列の場合
> mt1 <- matrix(c(1,3,5,6,2,3), ncol=3, byrow=F)
> vec1 <- c(10,100)
> vec2 <- c(10,100,1000)
> 
> mt1*vec1     # vec1の次元がmt1の行数と一致するので各行をvec2の要素倍する
     [,1] [,2] [,3]
[1,]   10   50   20
[2,]  300  600  300
> mt1*vec2     # vec2の次元がmt1の列数と一致するので各列をvec2の要素倍する
     [,1] [,2] [,3]
[1,]   10 5000  200
[2,]  300   60 3000
> mt1 %*% vec2 # vec2は縦ベクトルとして扱われ、積を計算する
     [,1]
[1,] 2510
[2,] 3630
> mt1 %*% vec1 # vec1の次元がmt1の列数に一致しないので計算できない
Error in mt1 %*% vec1 : non-conformable arguments
> 
> # 3×3の正方行列の場合
> mt2 <- matrix(c(1,3,5,6,2,3,7,7,8), ncol=3, byrow=F)
> mt2*vec2     # 行数とも列数とも一致するが基本的に縦扱いなので、各行をvec2の要素倍する
     [,1] [,2] [,3]
[1,]   10   60   70
[2,]  300  200  700
[3,] 5000 3000 8000
> mt2*t(vec2)  # 転置しても各列を要素倍するようになってくれるわけではなくエラーに
Error in mt2 * t(vec2) : non-conformable arrays
> t(t(mt2)*vec2)  # こうすると各列を要素倍してくれたことになる
     [,1] [,2] [,3]
[1,]   10  600 7000
[2,]   30  200 7000
[3,]   50  300 8000
> 
> # 行列同士をかける場合
> mt3 <- matrix(c(1,3,5,6), ncol=2, byrow=F)
> mt4 <- matrix(c(3,7,7,8), ncol=2, byrow=F)
> 
> mt3*mt4   # これは各要素同士の積
     [,1] [,2]
[1,]    3   35
[2,]   21   48
> mt3%*%mt4  # これは行列の積
     [,1] [,2]
[1,]   38   47
[2,]   51   69
  • 「Error: vector memory exhausted (limit reached?)」というメモリが足りないエラーは、RStudioの設定が問題になっている可能性がある(参考リンク)(参考ついったー
  • たくさんのデータフレームをlistでつないであって、それを全部一気にrbindしたいときの方法がこのページに載っていて参考になる。以下のどちらの書き方でも速い。
single.df <- bind_rows(list.of.dfs)  # これはdplyrの関数
single.df <- do.call("rbind", list.of.dfs)
  • 1列しかないデータフレームについて、行番号で抽出を行うと、単なるベクトルになってしまう。
> d <- data.frame(Term=c('春','夏','秋','冬'))
> class(d)
[1] "data.frame"
> class(d[1:2,])
[1] "character"
  • dplyrでの昇順・降順のソートはまじでクソ覚えにくい。昇順ならarrange(col)で、降順の場合はarrange(desc(col))
  • dplyrで「行の合計」に基づいてフィルターしたい場合はrowSums(.)を使えば良い。
df <- data.frame(x1=c(1,2,3),x2=c(2,3,4),x3=c(1,2,4))
df %>% filter(rowSums(.)>=7)
df %>% filter(rowSums(.)==7)
  • dplyrのgroup_byで複数のカラムを設定し、summariseすると、以下のよう謎のメッセージが出るのだが、あまり気にしなくていいようだ。summariseすると「最後に設定したカラムでのグループ化」が解除されるという仕様になっている(参考1参考2)。この警告自体は、意図どおりに2つのカラムでグループ化&サマライズできなかったという意味ではなく、グループ化が1個解除されてこうなりましたという意味かな。
`summarise()` has grouped output by 'xxxxx'. You can override using the `.groups` argument.
  • 行列の列を-インデックスで除いた場合、除いた結果が1列になる場合は行列型ではなく単なるベクトルになり、結果が0列になる(全部消える)が場合は行列のままとなる。前者に注意が必要で、1列になることもある処理があるときに、行列に対するコードを適用してるとエラーになる。
y1 <- cbind(c(1,2,3,4,5), c(6,7,8,9,10))
y2 <- y1[,-2]
class(y2)
y3 <- y1[,-c(1,2)]
class(y3)
  • リストの要素にインデックスでアクセスするときはlist1iとやるわけだが、ここにベクトルを与えると、複数の要素を取ってくるのではなく、再帰的なアクセスになる。どういうことかというと、list1c(1,2)となっている場合、list1の1つ目と2つ目の要素を取るのではなく、「1つ目の要素の中の2つ目の要素」を取りに行く。仮にlist11が行列だった場合、2個目の値ということで、その行列の1列目の2行目の値を取りに行く。list1の1つ目の要素と2つ目の要素を取りたい場合は、[1:2]と一重カッコにする(ややこしい!)。そもそも、リストの要素にアクセスするときにiとしているのは、Rのリストは中身の要素がリスト型でくるまれている(だからc()で追加するときに追加する要素をlist()する)からで、[i]だと返り値はリストになる。list1[1][1]とかにしてもムダで、list1[1]とlist[1][1][1][1][1][1]は同じ結果になる。list11とすることにより、アクセスできる。
> list1 <- list(matrix(c(1,2,3,4), nrow=2, byrow=T), 
+               matrix(c(5,6,7,8), nrow=2, byrow=T))
> list1[[1]]
     [,1] [,2]
[1,]    1    2
[2,]    3    4
> list1[[2]]
     [,1] [,2]
[1,]    5    6
[2,]    7    8
> list1[[c(1,2)]]
[1] 3
> list1[[1:2]]
[1] 3
> list1[1:2]
[[1]]
     [,1] [,2]
[1,]    1    2
[2,]    3    4

[[2]]
     [,1] [,2]
[1,]    5    6
[2,]    7    8
  • 0を0で割るとNaNになり、0以外の数字を0で割るとInfか-Infになる。これを忘れてると、たとえば「0で割る」ことになってしまうことがあり得ることを想定して、該当箇所を何かでreplaceするという処理をするときに、NaNとInfを両方探さないといけないのを忘れることになる。
>   0/0
[1] NaN
>   1/0
[1] Inf
>   -1/0
[1] -Inf
  • csvの文字コードを変換したいとき、read.csvで読み込む際はファイルに合わせてfileEncoding='CP932'とか指定し、write.csvするときにエンコーディングを何も指定せずにやると、綺麗にUTF-8に変換される。
  • Rの正規表現でURLを抽出する方法は以下のとおり。perlをTRUEにするのを忘れないように。
grep('https?://[\\w!\\?/\\+\\-_~=;\\.,\\*&@#\\$%\\(\\)\'\\[\\]]+', x, perl=TRUE)
  • 正規表現で「AAAを含み、かつBBBを含み、かつCCCを含む」(AND条件)を表現する方法は以下の通り。肯定先読みを使う。
grep('^(?=.*AAA)(?=.*BBB)(?=.*CCC).*$', text_all, perl = T)
  • 正規表現で「AAAを含み、かつBBBを含まない」は否定先読みを使う。
grep('^(?=.*AAA)(?!=.*BBB).*$', text_all, perl = T)
  • whichをつかって「◯◯であるものを除く」という絞り込みをするとき、whichでヒットするものが1つもない場合、全部の要素が残ってくれるのではなく、全部消えてしまう。なので、1つもない可能性が少しでもあるなら、「◯◯でないものを残す」ように書かないとダメ。
> x <- c('あ','い','う')
> x[which(x != '')]
[1] "あ" "い" "う"
> x[-which(x == '')]
character(0)

[メモ]Mac版のRでパッケージのインストールに失敗(X Codeのアプデが必要だった)

さっきハマったエラーを解決したので、備忘のためのメモです。
RStudioとRの両方で試したのですが、{panelvar}というパッケージをインストールしようとしたら、依存パッケージのコンパイルのエラーがでて進めなくなりました。
最初、Rの古いバージョンでやってたら{progress}というパッケージで躓いたんですが、Rを最新版にアプデしてからやると今度は{vctrs}というパッケージがインストールできてませんでした。
{vctrs}だけをインストールしようとしてもエラーが出る。

ERROR: compilation failed for package ‘vctrs’
* removing ‘/Library/Frameworks/R.framework/Versions/3.6/Resources/library/vctrs’
Warning in install.packages :
  installation of package ‘vctrs’ had non-zero exit status

The downloaded source packages are in
	‘/private/var/folders/7y/j9v534hs3hb971070xnz91l80000gp/T/Rtmp9EalBG/downloaded_packages’
Error: package or namespace load failed for ‘panelvar’:
 .onLoad failed in loadNamespace() for 'pillar', details:
  call: utils::packageVersion("vctrs")
  error: there is no package called ‘vctrs’
* installing *source* package ‘vctrs’ ...
** package ‘vctrs’ successfully unpacked and MD5 sums checked
** using staged installation
** libs
clang -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk -I/usr/local/include  -fPIC  -Wall -g -O2  -c arg-counter.c -o arg-counter.o
clang: warning: no such sysroot directory: '/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk' [-Wmissing-sysroot]
In file included from arg-counter.c:1:
In file included from ./vctrs.h:2:
/Library/Frameworks/R.framework/Resources/include/R.h:55:11: fatal error: 'stdlib.h' file not found
# include <stdlib.h> /* Not used by R itself, but widely assumed in packages */
          ^~~~~~~~~~
1 error generated.
make: *** [arg-counter.o] Error 1
ERROR: compilation failed for package ‘vctrs’
* removing ‘/Library/Frameworks/R.framework/Versions/3.6/Resources/library/vctrs’
Warning in install.packages :
  installation of package ‘vctrs’ had non-zero exit status

The downloaded source packages are in
	‘/private/var/folders/7y/j9v534hs3hb971070xnz91l80000gp/T/Rtmp9EalBG/downloaded_packages’


これは、こないだのMac OSのアプデ時に、X Code最新化(コマンドラインツールのインストール)しておくべきだったのが原因のようです。
下記のサイトで2段階の解決方法が紹介されてますが、私は1段階目の"xcode-select"で治りました。
_cgo_export.c:3:10: fatal error: &#39;stdlib.h&#39; file not found の解決方法 - Qiita

Macで画像ファイルをまとめてPDF化する(ただし名前順にソート)

Macで画像ファイルをまとめてPDFにしたいとき、従来はAutomatorを使ってワークフローを組んでいました。
ところが、いつからか知らないですが、OSに標準でそういう機能が付けられていたようです。
ついさっき、ある人から、とある論文のPDF入手を頼まれたのですが、諸事情により画像ファイルでしか閲覧できないので、これをPDF化しようとしていて気づきました。


f:id:midnightseminar:20191119120648p:plain


上のように、ファインダーで複数の画像を選択すると、右下のところに(これはサービスというやつですが)「PDFを作成」というボタンがあって、これを押すと選択した画像をまとめたPDFがそのフォルダに作成されます。


これは神!と思ったのですが、どうも、ファインダーで選択した順番にPDFに入っていくようです。
個人的には、どちらかといえば「選択したやつ全部を名前の順にソートしてPDF化」あるいは「フォルダ内の画像を名前順にまとめてPDF化」したい場合が多いので、やっぱりAutomatorで組んでおくことにしました。
作業としては2分ぐらいで終わります。


まずAutomatorを開いて新規作成で「クイックアクション」を選択。これはOS内で「サービス」と呼ばれる機能を提供するものです。


f:id:midnightseminar:20191119121220p:plain


次に、ワークフローの画面の一番上にあるところは(右の上)では、「ワークフローが受け取る現在の項目」は「ファイルまたはフォルダ」にしておき、「検索対象」は「Finder.app」を選んでおきます。
その上で、左の処理一覧から、処理したい項目を次のように選びます(右のワークフロー欄にドラッグする)。

  1. 「ファイルとフォルダ」から「選択されたFinder項目を取得」を選択
  2. 「ファイルとフォルダ」から「Finder項目にフィルタを適用」を選択し、「種類」が「イメージ」であるという条件を設定
  3. 「ファイルとフォルダ」から「Finder項目を並べ替える」を選択し、「名前」の「昇順」でソートするようにしておく
  4. 「PDF」から「イメージから新規PDFを作成」」を選択し、保存先は適当に決めておいて、出力するときの名前も決める。サイズは「ページをサイズに合わせる」にしておけば、勝手に画像が圧縮されたりしなくていいと思います。


f:id:midnightseminar:20191119121922p:plain


これだけやって保存しておけば、ファインダーで画像を選択したときに、さきほどの「PDFを作成」の隣にある「その他」の中にメニューが出てきますし、あるいは画像を選択した状態で右クリックして「サービス」から選択することもできます。


f:id:midnightseminar:20191119122405p:plain
f:id:midnightseminar:20191119122423p:plain


なお、サービスに追加したワークフローは、"~/Library/Services"という場所に入っていて、名前を変えたいときはここでファイル名を変えればいいようです。
また、フォルダを選択してフォルダ内の画像をまとめたい場合は、ワークフローのパーツとして「フォルダの内容を取得」ってのを使えばいいです。

言語学への関心は衰退している?

とある雑誌の連載記事に、

言葉というものは曖昧かつ不安定で捉えにくい対象であることもあって、とりわけ現代思想ブームが終焉し実証的社会科学が隆盛を極めているここ三十年ほどの間は、言語理論への関心は総じて低調であったと言える。「言語論的転回」は今や、むしろ逆転の渦中にあるとみたほうがよいのかも知れない。


というようなことをあまりよく考えずにさらっと書いてしまったのですが、本当に低調であると言えるのかどうか不安になったので、Google Ngram Viewer(特定のキーワードが19世紀以降の英語の書籍にどれぐらい頻繁に登場してるかを教えてくれるやつ)で確認してみました。


【言語学】
f:id:midnightseminar:20191006190050p:plain


30年というよりは40年ぐらいかけて下落していっていますが、日本における現代思想ブームは欧米から10年ぐらい遅れていたイメージがあるので、大きく間違ってはいなかったと安心することにしました。
構造主義〜ポスト構造主義(ポストモダニズム)の時代は、人文学全般が非常に盛んで、あの時代が終わると「哲学」や「思想」とかいうもの全般に対する関心が薄れていきましたよね。まぁ私は80年代生まれなのでリアルタイムでは知りませんが。
関連しそうな用語も調べてみました。(なお、言語学や心理学の盛衰に関しては、現代思想ブームと並行して起こった、50-70年代のいわゆる「認知革命」というムーブメントも関係してると思います。)


【意味論】
f:id:midnightseminar:20191006190841p:plain


【記号論】
f:id:midnightseminar:20191006190856p:plain


【哲学】
f:id:midnightseminar:20191006190911p:plain


哲学(philosophy)はとくに人文学の議論でなくても用いられる用語なので、あまり現代思想ブームとは関係なかったかもしれません。


【人文学】
f:id:midnightseminar:20191006191153p:plain


「人文学」は90年代ごろから下落が始まるようです。


【構造主義】
f:id:midnightseminar:20191006191253p:plain


「構造主義」も意外とピークは90年代だったようです。


【人類学】
f:id:midnightseminar:20191006191237p:plain


【心理学】
f:id:midnightseminar:20191006191342p:plain


「心理学」は、ずいぶん昔にも山がありますが、これは行動主義の趨勢に引っ張られてるんですかね?
30年代ごろはパブロフとワトソン、50年ごろはスキナーの時代という感じですね。
「行動主義」も検索してみました。


【行動主義】
f:id:midnightseminar:20191006191915p:plain


【社会学】
f:id:midnightseminar:20191006191937p:plain


社会学は、ルーマンやパーソンズのような大物理論化がいた時代にピークがありますね。


人文学とは言わないと思いますが(19世紀には人文学と社会科学の区別がなかったので両者は似たようなものだったわけですが)、政治学と経済学も見てみました。


【経済学】
f:id:midnightseminar:20191006192201p:plain


【政治学】
f:id:midnightseminar:20191006192319p:plain


なんとなく、文系全体が衰退してるような感じがしますねw


しかし途中で気づいたのですが、「数学」とか「物理学」とか「統計学」とかを入力しても90年代から2000年代にかけて衰退しているので、アカデミックなワード自体が割合としては減っていくような、サンプルの性質が何かあるのかもですね。(もちろん、たとえばstatisticsなんかはとくに、アカデミックな話で使われるとも限りませんが。)


【数学】
f:id:midnightseminar:20191006192953p:plain


【物理学】
f:id:midnightseminar:20191006193010p:plain


【統計学】
f:id:midnightseminar:20191006193124p:plain

メモ:Rのts型のインデックス(年・四半期・月情報)を文字列として取り出す

Rのts型(時系列型)のデータについているインデックス("1980 Q1"みたいな)を、文字列情報として取り出す方法が、ぱっとググって分からなかったのですが、とりあえず以下のようにしたらできました。

> library(vars)
> data(Canada)
> 
> t <- as.yearqtr(index(Canada))
> head(t)
[1] "1980 Q1" "1980 Q2" "1980 Q3"
[4] "1980 Q4" "1981 Q1" "1981 Q2"
> class(t)
[1] "yearqtr"
> 
> t <- as.character(t)
> head(t)
[1] "1980 Q1" "1980 Q2" "1980 Q3"
[4] "1980 Q4" "1981 Q1" "1981 Q2"
> class(t)
[1] "character"


たぶん、元データが四半期ではなく月単位のデータだったら、yearqtrの部分をyearmonにするのだと思います。