Rで解析：クラスター誤差を考慮した操作変数法「clusterIV」パッケージの紹介

経済学や社会学などの社会科学分野では、因果関係を推定する際に「操作変数法（インストルメンタル変数法）」が使用されています。しかし、データにクラスター構造がある場合やインスツルメントが多数あると、通常の方法ではバイアスが生じることがあります。この問題を解決する手段として、「clusterIV」パッケージを紹介します。

本パッケージは、単一の内生変数を持つ単純識別設計でのバイアスを最小限に抑える、クラスター削除ジャックナイフ・インストルメンタル変数推定（CJIVE）やクラスターロブスト推定と多重インスツルメンタル変数を扱うためのコマンドが収録されています。

このパッケージの強みは、クラスター削除によるフィットが効率的に行われることで、大規模なデータにも対応できる点です。各観測値の第一段階値をそのクラスタを除いてフィットすることで、クラスター間のバイアスを解消し、さらにインスツルメントが多数ある場合でも正確な推定を行うことができるようです。

パッケージバージョンは0.1.0。Windows 11 x64 (build 26200)のR version 4.6.0で確認しています。

＜おすすめのRに関する書籍です＞

遺伝統計学の基礎 Rによる遺伝因子解析・遺伝子機能解析

遺伝因子解析・遺伝子機能解析をRで行うための教科書。ゲノム解析に特有の統計的な考え方を学べる。

www.amazon.co.jp

パッケージのインストール

下記コマンドを実行してください。

# パッケージのインストール
install.packages("clusterIV")
# パッケージの読み込み
library("clusterIV")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

cjive() の使用例

# 乱数の固定
set.seed(1)

### データの準備#####
# グループ数を設定
G  <- 40
# グループ内の観測値数を設定
ng <- 6;
# 全観測値数を計算
n <- G * ng
# クラスタリング変数を作成
cl <- rep(seq_len(G), each = ng)
# 判定IDを作成
j  <- factor(rep(rep(1:4, length.out = ng), G))   
# 階層効果を生成
u  <- rnorm(G)[cl]
# 解説変数を作成
x  <- as.numeric(j) + u + rnorm(n)
# 目的変数を作成
y  <- 1.5 * x + u + rnorm(n)
########

# クラスター削除ジャックナイフ・インストルメンタル変数推定:cjiveコマンド
# yオプション：数値ベクトルまたは、y ~ x | zという数式
# xオプション：単一の内生変数
# zオプション：インストゥメント：数値ベクトル/行列、またはジャッジデザイン用の因子/文字列グループ化ベクトル
# clusterオプション：クラスター識別子
# controlsオプション：オプションの共変量（FLMのX）：行列またはデータフレーム、または公式方法での1辺の公式
# weightsオプション：オプションの正の精度重み
# levelオプション：報告する区間の信頼水準
# interceptオプション：論理値；インターセプトを部分的に除去（デフォルトTRUE）
# methodオプション："auto", "dense", "leaveout_mean"
# dataオプション：評価対象のデータフレーム。
# digitsオプション：表示する桁数。
# parmオプション：無視される（単一の係数が推定される）。
fit <- cjive(y, x, j, cluster = cl)

# 結果を表示
print(fit)

## 数式で指定の場合
# データフレームを作成
dat <- data.frame(y = y, x = x, j = j, cl = cl)
# 実行
cjive(y ~ x | j, data = dat, cluster = ~cl)

実行結果：

# cjive(y, x, j, cluster = cl)とcjive(y ~ x | j, data = dat, cluster = ~cl)は同じ結果
Cluster-jackknife IV (CJIVE)
Call: cjive.formula(formula = y ~ x | j, data = dat, cluster = ~cl)

  coefficient = 1.536   cluster-robust SE = 0.0579
  z = 26.52   p = < 2.2e-16   95% CI = [1.422, 1.649]
  n = 240   G = 40 clusters   p = 3 instruments   path = dense
  max within-cluster leverage = 0.025

異なるIV推定器をクラスターロブустな標準誤差で比較する例です。

# 乱数の固定
set.seed(2)

### データの準備#####
# グループ数を設定
G <- 50; 
# グループ内の観測値数を設定
ng <- 5
# 全観測値数を計算
n <- G * ng
# クラスタリング変数を作成
cl <- rep(seq_len(G), each = ng)
# 工具変数行列を作成
z  <- matrix(rnorm(n * 3), n, 3)
# 階層効果を生成
u  <- rnorm(G)[cl]
# 解説変数を作成
x  <- z %*% c(1, -1, 0.5) + u + rnorm(n)
# 目的変数を作成
y  <- 2 * x + u + rnorm(n)
########

# iv_compare分析を実行:iv_compareコマンド
# yオプション：結果変量（数値ベクトル）
# xオプション：単一の内生変数（数値ベクトル）
# zオプション：インストルメント（数値行列/ベクトルまたはグループ化因子）
# clusterオプション：クラスター識別子（長さn）。1辺の公式(~ g)や列名でも可能
# controlsオプション：オプションの共変量（FLMのX）：行列またはデータフレーム、または公式方法での1辺の公式。ランク不足を許可（固定効果が可能）。インターセプトは追加されない場合がある；初期値NULL
# weightsオプション：オプションの正の精度重み；初期値NULL
# levelオプション：報告される区間の信頼水準；初期値0.95
# interceptオプション：論理値；インターセプトを部分的に除去（デフォルトTRUE）；初期値TRUE
iv_compare(y, x, z, cluster = cl)

実行結果：

estimator coefficient         se statistic       p.value conf.low conf.high
1       OLS    2.221181 0.05743338  38.67404  0.000000e+00 2.108613  2.333748
2      2SLS    2.034037 0.05567786  36.53224 3.413580e-292 1.924910  2.143164
3      JIVE    2.031592 0.05603122  36.25822 7.374238e-288 1.921773  2.141412
4     CJIVE    2.031094 0.05587636  36.34979 2.647444e-289 1.921578  2.140610

＜おすすめのRに関する書籍です＞

東京大学のデータサイエンティスト育成講座 ~Pythonで手を動かして学ぶデ―タ分析~

データ分析の一連の流れを手を動かして学ぶ定番書。ただしコードはRではなくPython中心なので、その点は注意。

www.amazon.co.jp

この記事が誰かの役に立ちますように。