Hybrid ModelをRでやる（特にpanelrパッケージについて）

ひとつ、ひとつだ。

Introduction
- 問題の所在
- ハイブリッドモデルに関する参考URL等
RでHybrid Modelを使うーpanelrパッケージ
Conclusion

Introduction

問題の所在

一部界隈で使われているパネルデータに対する手法としてハイブリッド・モデルというものがある。

$Y_{ij} = \beta_0 + \beta_1(X_{ij} - \bar{X}_i) + \beta_2 \bar{X}_i+ \beta_3 Z_i + \alpha_i + \epsilon_{ij}$

といった形で固体内偏差を投入する。

こうすると、Lv.1の変数（個人内効果）の推定値は固定効果モデルにおけるそれと一致するし、個体間の変動は時変変数の個人平均 $\bar{X}_i$ の係数 $\beta_2$ や、時不変変数 $Z_i$ の係数 $\beta_3$ によって説明することができるようになる。

これをRで使おうとしたときに若干色々めんどくさいので、そのめんどくさい点への対処法を備忘。

ハイブリッドモデルに関する参考URL等

RでHybrid Modelを使うーpanelrパッケージ

wbm()関数 for hybrid model

上記のpanelrのvignetteには、開発者の以下のような意気込みが宣言されている。

The original motivation to create this package was to automate the process of fitting “within-between” models, sometimes called “between-within” or “hybrid” models (see Allison, 2009; Bell & Jones, 2015).

このパッケージを作ったモチベーションはハイブリッドモデルにありました、とのことである
そしてそのためにつくられた関数が、wbm()である。

The workhorse function for within-between models is wbm(), which is built on top of lme4’s lmerMod() and glmerMod().
It is not so hard to understand how to treat your data to estimate within-between models,
but the programming can be a challenge to those who aren’t skilled with R (or whatever else they might use) and is error-prone in any case.

別にlme4の各関数でもhybrid modelは推定できるけど、色々ミスが起こりやすいからわかりやすくしたよってことらしい。
以下、その使い方等を記す。なお、テストデータとしてはpanelrパッケージに含まれているWageDataをつかう

基本：hybrid modelの推定

data(WageData)
wdf <- panel_data(  WageData, id = id, wave =t)
hyb1 <- wbm( lwage ~ occ + wks | occ + wks ,
data = wdf,
 model="w-b",
min.waves = 1)

formulaの指定の仕方が特徴的でひとつめの|の前はwithin（個体内変化）の変数、ひとつめの|のあとはbetween（個体間差異）の変数を指定する。
引数min.wavesは、一個体に関して最低いくつ以上の観測時点が必要かを指定する。ここでは後の固定効果モデルとの比較上、1にしている。
得られた結果は以下の通り。

summary(hyb1)
# MODEL FIT:
#   AIC = 2181.85, BIC = 2226.2
# Pseudo-R² (fixed effects) = 0.11
# Pseudo-R² (total) = 0.69
# Entity ICC = 0.64471
#  
# WITHIN EFFECTS:
# --------------------------------------------------------------------------
#   Est.        S.E.       t val.           d.f.           p
# --------- ------------ ----------- ------------ -------------- -----------
# occ         -0.0709817   0.0232268   -3.0560316   3568.0000042   0.0022595
# wks          0.0011319   0.0010203    1.1094052   3568.0000042   0.2673302
# --------------------------------------------------------------------------
#   
# BETWEEN EFFECTS:
# ----------------------------------------------------------------------
#   Est.        S.E.        t val.          d.f.
# ----------------- ------------ ----------- ------------- -------------
# (Intercept)          6.3582573   0.2137905    29.7406011   591.9985849
# imean(occ)          -0.3219410   0.0317292   -10.1465052   591.9999818
# imean(wks)           0.0103106   0.0045345     2.2738038   591.9985934
# ----------------------------------------------------------------------
#

plmとの比較（wbmで固定効果モデルの方法も）

参考URL①にある通り、hybrid modelでは固定効果モデルと時変変数の係数推定値が一致する。

レベル1すなわち個人内の係数推定値について、固定効果モデルでも、変量効果モデルの一種であるHybrid Modelでも、結果は一致する。
単純な変量効果モデルないしマルチレベルモデルによる分析では、しばしばレベル1変数の効果の推定が偏りをもつ危険性が指摘されてきた。
しかしHybrid Modelのように、当該レベル1変数の個人内平均をレベル2平均として同時に投入することで、偏りを避けることができ、固定効果と等しい結果を得ることができる（三輪・山本2012:72）

ということで、上と同じデータを使って、plm()で固定効果モデルを推定し、係数推定値を比較してみる。

fix1 <- plm( lwage ~ occ + wks,data = WageData  , model="within", index =c("id","t"))
summary(fix1)

# Coefficients:
#   Estimate Std. Error t-value Pr(>|t|)   
# occ -0.0709817  0.0232268 -3.0560 0.002259 **
# wks  0.0011319  0.0010203  1.1094 0.267330

確かにさきにみたHybrid ModelのLevel1効果と係数推定値は一致しているようである。
同じモデルをwbmで推定してみる。

fix2 <- wbm( lwage ~ occ + wks ,data = wdf, model="within", min.waves = 1)
summary(fix2 , digits=6)

# ------------------------------------------------------------------------------
#   Est.       S.E.       t val.          d.f.          p
# ----------------- ----------- ---------- ------------ ------------- ----------
# (Intercept)          6.676346   0.016162   413.083707    593.999993   0.000000
# occ                 -0.070982   0.023227    -3.056032   3568.000002   0.002259
# wks                  0.001132   0.001020     1.109405   3568.000002   0.267330
# ------------------------------------------------------------------------------

同様の推定値が得られているようだ。

lmerとの比較（hybrid modelの別の推定方法）

hybrid modelは要するに、時変の変数を個体内平均と個体内変動に分解して導入した変量効果モデルであるので、
むろんlmerでも推定はできる*1。

WageDataから、occ, wksの2変数に関して個体内平均、個体内変動に分割した変数をつくる

wdf_r <- as.data.frame(dplyr::ungroup(
dplyr::mutate( 
dplyr::group_by( WageData , id), 
occ_m = mean(occ, na.rm=T),
wks_m = mean(wks,na.rm=T))))
wdf_r <- dplyr::mutate( wdf_r , wks_w = wks- wks_m, occ_w = occ - occ_m)

そして、between/within両方に対応した変数を入れてlmerで推定する

library(lmerTest)
hyb2 <- lmer(lwage~ occ_w + wks_w + occ_m + wks_m + (1|id),data = wdf_r, REML=F)
summary(hyb2)$coefficients

# Estimate  Std. Error        df    t value      Pr(>|t|)
# (Intercept)  6.358257313 0.213250844  594.9996  29.815860 3.497717e-120
# occ_w       -0.070981734 0.023220260 3570.0001  -3.056888  2.253033e-03
# wks_w        0.001131949 0.001020035 3570.0001   1.109716  2.671961e-01
# occ_m       -0.321940988 0.031649160  594.9999 -10.172181  1.619229e-22
# wks_m        0.010310578 0.004523061  594.9996   2.279558  2.298709e-02