www.crosshyou.info

政府統計の総合窓口のデータや、OECDやUCIやのデータを使って、Rの練習をしています。ときどき、読書記録も載せています。

都道府県別のあんま・マッサージ師、はり・きゅう師、柔道整復師数のデータの分析8 - Rのplmパッケージでパネルデータ分析をする。pooling法、first difference法、fixed effect法、random effect法の4つの方法を実行する。

f:id:cross_hyou:20220211082343j:plain

Photo by Mike Swigunski on Unsplash 

www.crosshyou.info

今回はplmパッケージを使ってパネルデータ分析をしてみたいと思います。

まずはlibrary(plm)と入力してパッケージを読み込みます。

f:id:cross_hyou:20220211082932p:plain

パネルデータ分析をするには、データフレームをパネルデータフレームに変換する必要があります。pdata.frame()関数で変換します。

f:id:cross_hyou:20220211083259p:plain

index = c("個々の識別変数名", "時系列の識別変数名")という指定をするとパネルデータフレームになります。

pdim()関数でパネルデータになっているか確認します。

f:id:cross_hyou:20220211083656p:plain

あれ?Unbalanced Panelとなっていますね。。データに欠損のある都道府県か、観測年があるのかな?ちょっと調べてみます。

f:id:cross_hyou:20220211083954p:plain

都道府県コードが4000番が観測数が25でどこか1年欠損していますね。

何年が欠損しているか調べます。

f:id:cross_hyou:20220211084325p:plain

2010年が都道府県コード4000番が欠損していました。

2010年を除いてパネルデータを作り直します。

f:id:cross_hyou:20220211084640p:plain

Balanced Panelになりました。

それでは、パネルデータを使って回帰分析をしてみます。

分析するモデルは、anmassage = b0 + b1*harikyu_pop + b2*judo_pop + u でやってみます。第5回目のクロスセクションでの回帰分析、第6回目の東京都のデータを使った時系列での回帰分析では、はり・きゅう師の数が増えるとあんま・マッサージ師の数も増えていました。パネルデータではどうなるでしょうか?

plmパッケージでは、plm()関数で回帰分析ができます。今回はパネルデータを通常のクロスセクションデータと同じように扱うpooling法、first difference法、fixed effect法、random effect法の4つの方法を試してみます。

まずは、pooling法です。

f:id:cross_hyou:20220211085633p:plain

harikyu_popの係数は0.62でプラス、judo_popの係数は-0.50でマイナス、どちらのp値も有意ですね。

次はfirst difference法です。

これは、ベースとなるモデルを

anmassage_pop(i, t) = b0 + b1*harikyu_pop(i, t) + b2*judo_pop(i, t) + a(i) + u(i, t) と想定します。i が各都道府県、tが各観測年を表して、a(i)というのが各都道府県特有の係数でこれは時間が経っても変わらずに一定だと仮定しています。そしてこのa(i)を削除するために、t=2とt=1の差分を考えます。

anmassage_pop(i, t) - anmassage_pop(i, t-1) = Δanmassage_pop(i, t) =

b1*Δharikyu_pop(i, t) + b2*Δjudo_pop(i, t) + Δu(i, t)

となります。Δa(i) = a(i) - a(1) = 0 なので差分の式だと、aが消えるということです。

plm()関数でmodel = "fd"とするとfirst difference法になります。

  

f:id:cross_hyou:20220211091539p:plain

harikyu_popの係数は0.67でプラス、judo_popの係数は-0.44でマイナス、どちらも有意です。

次は、fixed effect法です。これは差分を取るのでなくて、

avg(anmassage_pop(i) = b1*avg(harikyu_pop(i)) + b2*avg(judo_pop(i)) + a(i) + avg(u(i))とそれぞれの都道府県の平均値を計算してそれと元のモデル式の差分を取ってa(i)を消去するという方法です。

model = "within"とするとfixed effect法になります。

f:id:cross_hyou:20220211092516p:plain

harikyu_popの係数は0.61でプラス、judo_popの係数は-0.90でマイナスです。どちらの係数も有意です。

最後はrandom effect法です。model = "random" と指定します。

f:id:cross_hyou:20220211093204p:plain

harikyu_popの係数は0.61でプラス、judo_popの係数は-0.89でマイナス、これも両者有意ですね。

こうして4つの方法でharikyu_pop, judo_popの係数を推計しましたが4つともharikyu_popの係数はプラス、judo_popの係数はマイナスで、統計的に有意な値でした。

最後に4つの方法のharikyu_popの係数をグラフにして比較します。

まず、broomパッケージのtidy()関数を利用してharikyu_popの係数と信頼区間をデータフレームの形式でまとめます。

f:id:cross_hyou:20220211095110p:plain

次にggplotのgeom_errorbar()関数でグラフにします。

f:id:cross_hyou:20220211095805p:plain

f:id:cross_hyou:20220211095838p:plain

first difference法での推測が少し他の3つははずれていますが、はり・きゅう師の数が増えるとあんま・マッサージ師の数が増えるというのは間違いなさそうです。というか何か別の原因が作用してあんま・マッサージ師の数とはり・きゅう師の数に影響をしているのでしょう。そしてそれは、柔道整復師の数を減少させる作用があるのかもしれません。

今回は以上です。

初めから読むには、

 

www.crosshyou.info

です。