R&D研究員の小松です。本連載ではネットワーク経済学の近年の発展について、忘備録的に書かせてもらっています。

前回まではネットワークにおける相互依存関係を考慮したモデル ERGM のミクロ的基礎づけとその推定・応用について概観してきました。

今回はその話題から離れ、dyadic regression、特にその standard error の計算方法について紹介します。

Dyadic regression

Dyad とは、ネットワーク科学界隈ではノード同士のペア (リンクがあってもなくてもよい) を指します。 Dyadic data とは、そうした dyad を観察単位としたデータのことです。世界の各空港間で月間フライトが何便飛んでいるのかは dyadic data の一例です。こうした dyadic data を対象とした回帰分析、dyadic regression が社会科学ではしばしば見られます。アウトカム $y_{ij}$ が、dyad を構成する $i$ と $j$ によって決まるものであり、他の dyad との相互依存関係はここでは考慮しません。

Dyadic regression の (主にパラメトリックな) 分析手法については Graham (2020) を参照いただくとして、dyadic regression の近年の応用としては以下のようなものがあります。これらの詳細はまた別で読んでいきたいと考えています。

Harmon, Fisman, and Kamenica (2019): 欧州議会の座席ルールを用いて、立法投票における peer effect の存在を実証。

Roche, Oettl, and Catalini (2022): 企業同士間の距離が、企業間の知識スピルオーバーに与える影響について、アメリカのあるコワーキングハブにおけるスタートアップ企業群のデータを用いて実証。

本稿では、以下のような dyadic linear regression を考え、パラメーターの推定値の標準誤差をどう計算するかを考えます。

$y_{ij} = \alpha + x^{\top}_{ij} \beta + \varepsilon_{ij}$

ここでは、過去の記事で扱った ERGM のような、dyad 間の相互依存関係をモデルには組み込みません。あくまで $y_{ij}$ は $i$ と $j$ の間だけで決まるとしています。

Cluster-robust standard errors for dyadic data

回帰分析において、各観察がクラスター (e.g., 市区町村、都道府県、産業...) によってグルーピングされると考えることが妥当なケースがあります。つまり、クラスター間では各観察は独立であるが、クラスター内では各観察同士に相関があるケースです。そうしたクラスター構造を無視して標準誤差を計算すると、本来の分散よりもはるかに小さい値となってしまい、帰無仮説を over-rejection してしまうことが知られています。クラスターを考慮して計算された標準誤差が cluster-robust standard error です。

dyadic regression では観察単位が dyad、つまりある当事者同士のペアなわけですが、各 dyad はそれを構成するメンバーを媒介し、他の dyad と相互に依存します。例えば、

$y_{ij} = \alpha + x^{\top}_{ij} \beta + \varepsilon_{ij}$

で、誤差項に $i$ と $j$ の観察されない要因が含まれる、つまり $\varepsilon_{ij} = \alpha_{i} + \alpha_{j} + \epsilon_{ij}$ のようなケースを想定します。このとき、dyad $ik$ $(k \neq i, j)$ の誤差項は $\varepsilon_{ik} = \alpha_{i} + \alpha_{k} + \epsilon_{ik}$ となります。 $ij$ と $ik$ の誤差項には $\alpha_{i}$ が共通して含まれるため、 $ij$ と $ik$ の間には相関が発生します。

以上から、少なくとも以下のようなクラスター構造が存在すると考えるのが dyadic regression の場合妥当のように考えられます。*1

2つの dyad に共通して含まれるメンバー (上の例だと $i$ ) が存在するとき、それらの dyad 間には相関がある
2つの dyad に共通して含まれるメンバーが含まれないとき、それらの dyad 間には相関がない

こうした構造を Aronow, Samii, and Assenova (2015) は dyadic clustering と呼んでいます。そして彼らはそのクラスター構造を考慮した dyadic cluster-robust standard error を提案しています。

具体的な例で考えてみます。以下のような dyadic data があるとし、そして $y$ を $x$ に回帰することを考えます。

このとき、 $x$ の係数の推定値の標準誤差をどう計算するかが本稿の主題なわけですが、1つの計算方法としてナイーブに heteroskedasticity-robust standard error を計算するものがあるでしょう。

reg <- fixest::feols(y ~ x, data, se = "hetero")

しかし先の dyadic clustering の話を思い出すと、上で計算した standard error はそうしたクラスターを一切考慮していないため、帰無仮説を over-rejection してしまいそうです。

別の方法として、各 dyad を構成するメンバーの id を用いた twoway cluster-robust standard error を計算するものがあるでしょう。

reg <- fixest::feols(y ~ x, data, cluster = c("i", "j"))

この場合、同じ $i$ を含む dyad 同士、または同じ $j$ を含む dyad 同士の相関を考慮した standard error を計算します。下の画像の例でいうと、黄色の行と同じ $i$ を持つ行は1行目から3行目、黄色の行と同じ $j$ を持つ行は1行目のみです。つまり、黄色の行について相関を考慮できる行は、赤の実線で囲まれた部分になります。

しかし、4行目と5行目の dyad には構成メンバーに 2 が存在しており、1行目の dyad の構成メンバーにも 2 がいます。したがって、4行目と5行目についても1行目と相関があるとして、standard error を計算する際にそれを考慮したいところです。しかし、上で見たように two-way cluster-robust standard error では、図の黄色の行と赤の点線で囲まれた部分について相関を考慮することができません。つまり、two-way cluster-robust standard error を dyadic data に応用したとき、 $i \lt j \lt k$ とすると、 $d_{ij}$ と $d_{ik}$ 間の相関は考慮できるが、 $d_{ij}$ と $d_{jk}$ 間の相関を考慮することができません。これはデータ構造的に $i$ と $j$ を上図のように格納しないといけないために生ずる問題で、概念的には dyadic cluster-robust standard error の計算方法は two-way clustering とほぼ同じです。しかし、上の R のコードのように単純にクラスターを指定するだけでは、dyadic cluster-robust standard error をうまく計算できません。

以上の部分を考慮できるようにしたのが、Aronow, Samii, and Assenova (2015) が提案した dyadic cluster-robust standard error です。dyad 間で同じメンバーがいるかどうかを考慮して standard error を愚直に計算するのは大変そうですが、Aronow, Samii, and Assenova は以下の式で dyadic cluster-robust standard error を計算できることを示しました。

$\begin{align} \hat{V}_{r} = \sum_{i=1}^{N} \hat{V}_{C, i} - \hat{V}_{D} - (N-2) \hat{V}_{0} \end{align}$

ただし、 $\hat{V}_{C, i}$ は、メンバー $i$ を含む全ての dyad 間には相関があり、それ以外の dyad 間には相関がないとして計算された cluster-robust variance estimator、 $\hat{V}_{D}$ は同じ dyad 内の相関を許した cluster-robust variance estimator (データがパネルデータではなければ、これは robust standard error と等しい)、 $\hat{V}_{0}$ は heteroskedasticity robust variance estimator です。 $\hat{V}_{r}$ を構成するそれぞれの推定量は既存の統計パッケージで計算できるものばかりですから、上の式に従えば dyadic cluster-robust variance estimator の計算は容易です。

この計算式は、複数の次元のクラスターを考慮して cluster-robust standard error を計算する際に、それらを各要素に分解して足し合わせ重複する部分を引くというというアイデアに基づきます (Cameron, Gelbach, and Miller, 2011)。例えばあるデータで、個人 $i$ と時間 $t$ について two-way cluster-robust variance を計算する場合、以下のように計算することができます。

$\hat{V} = \hat{V}_{i} + \hat{V}_{t} - \hat{V}_{i \times t}$

ただし、 $\hat{V}_{c}$ は $c$ に関する cluster-robust variance estimator です。 $i, t$ それぞれの次元について cluster-robust variance を計算して、その共通部分を引くイメージです。この分解式を応用して得られたものが、上で紹介した dyadic cluster-robust variance estimator の計算式です。

Monte Carlo simulation

実際にシミュレーションを行い、提案された推定量のパフォーマンスを評価してみます。そのために、dyadic cluster-robust standard error を計算するパッケージを作成しました *2

github.com

以下のようなシミュレーションを行います。

まず、Iteration $r$ $(r = 1, \ldots, R)$ において、 $y_{ij, r}$ を以下のように生成します。

$y_{ij, r} = \vert x_{i, r} - x_{j, r} \rvert + \alpha_{i, r} + \alpha_{j, r} + \varepsilon_{ij, r} \quad (i, j = 1, \ldots, N, i \lt j)$

ただし、 $x_{i, r}, x_{j, r}, \alpha_{i, r}, \alpha_{j, r}, \varepsilon_{ij, r} \sim \mathcal{N}(0, 1)$ で、互いに独立とします。

次に $\alpha_{i, r}, \alpha_{j, r}$ が unobservable だとして、 $y_{ij, r}$ を $\vert x_{i, r} - x_{j, r} \rvert$ に回帰させ、 $\vert x_{i, r} - x_{j, r} \rvert$ の係数 $\beta$ の推定値 $\hat{\beta}_{r}$ を得ます。

そして、robust standard error $SE_{RS, r}$ 、two-way cluster-robust standard error $SE_{TS, r}$ 、dyadic cluster-robust standard error $SE_{DS, r}$ を計算します。

このとき、真の standard error は

$\begin{align} S.D.(\hat{\beta}_{r}) = \sqrt{\frac{1}{R} \sum_{r=1}^{R}(\hat{\beta}_{r} - mean(\hat{\beta}_{r}))^{2}} \end{align}$

と計算できますから、これと $SE_{RS, r}$ 、 $SE_{TS, r}$ 、 $SE_{DS, r}$ の分布を比較して、各推定量を評価します。

$N = 200$ 、 $R = 1,000$ としてシミュレーションを行った結果が以下の図です。再現用のコードは本稿の末尾に添付しています。

破線が真の standard error の値で、 $SE_{RS, r}$ 、 $SE_{TS, r}$ 、 $SE_{DS, r}$ の分布を boxplot で図示しています。

図から明らかなように、robust standard error は真の standard error と比較しかなり過小に推定されています。two-way cluster-robust standard error には一定の改善が見られ、dyadic cluster-robust standard error の中央値は真の値にかなり近づいています。two-way cluster-robust standard error では考慮できなかった dyad 間の相関を、dyadic cluster-robust standard error では考慮できていることに、両者の差は起因していると考えられます。

Concluding remarks

本稿では、dyadic regression における、dyad 間の相互依存関係を考慮した standard error の計算方法を紹介しました。次回はその応用例として、Roche, Oettl, and Catalini (2022) の内容を紹介します。

References

Aronow, P. M., Samii, C., & Assenova, V. A. (2015). Cluster–robust variance estimation for dyadic data. Political Analysis, 23(4), 564-577.
Cameron, A. C., Gelbach, J. B., & Miller, D. L. (2011). Robust inference with multiway clustering. Journal of Business & Economic Statistics, 29(2), 238-249.
Graham, B. S. (2020). Dyadic regression. The Econometric Analysis of Network Data, 23-40.
Harmon, N., Fisman, R., & Kamenica, E. (2019). Peer effects in legislative voting. American Economic Journal: Applied Economics, 11(4), 156-180.
Roche, M. P., Oettl, A., & Catalini, C. (2022). (Co-) Working in Close Proximity: Knowledge Spillovers and Social Interactions. National Bureau of Economic Research.

Code for replication

library(magrittr)
library(ggplot2)

if (!"fastDyadRobust" %in% rownames(installed.packages())) {
  devtools::install_github("komatsuna4747/fastDyadRobust")
}


simulate_data <- function(n_node) {
  node <-
    tibble::tibble(
      id = 1:n_node,
      x = rnorm(n_node),
      a = rnorm(n_node)
    )

  dyad <-
    list(src = node$id, dst = node$id) %>%
    purrr::cross_df() %>%
    dplyr::filter(src < dst)

  df <-
    dyad %>%
    dplyr::left_join(node, by = c("src" = "id")) %>%
    dplyr::rename_at(c("x", "a"), ~paste0(., "_src")) %>%
    dplyr::left_join(node, by = c("dst" = "id")) %>%
    dplyr::rename_at(c("x", "a"), ~paste0(., "_dst")) %>%
    dplyr::mutate(
      dx = abs(x_src - x_dst)
    ) %>%
    dplyr::mutate(
      de = rnorm(nrow(.)),
      dy = dx + a_src + a_dst + de,
      dyad_id = paste(src, dst, sep = "-")
    )

  df %<>%
    dplyr::select(src, dst, dyad_id, dplyr::everything())

  return(df)
}


df <-
  as.list(1:1000) %>%
  purrr::map_df(
    function(x) {
      if (x %% 50 == 0) {
        logger::log_info("Iteration: {x}")
      }
      df_sim <- withr::with_seed(x + 1000, simulate_data(200))

      regHetero <- fixest::feols(fml = dy ~ dx, data = df_sim, se = "hetero")

      hetero <-
        broom::tidy(regHetero) %>%
        dplyr::mutate(se_type = "hetero", iter = x)

      regTwoway <- fixest::feols(fml = dy ~ dx, data = df_sim, cluster = c("src", "dst"))

      twoway <-
        broom::tidy(regTwoway) %>%
        dplyr::mutate(se_type = "two-way cluster-robust", iter = x)

      regDyadRobust <- fastDyadRobust::feolsDyadRobust(fml = dy ~ dx, data = df_sim, cluster = c("src", "dst"))

      dyadRobust <-
        broom::tidy(regDyadRobust) %>%
        dplyr::mutate(se_type = "dyadic cluster-robust", iter = x)

      output <-
        dplyr::bind_rows(hetero, twoway, dyadRobust) %>%
        dplyr::filter(term != "(Intercept)")

      return(output)
    }
  )


df %>%
  dplyr::filter(term == "dx") %>%
  ggplot(aes(x = forcats::fct_relevel(se_type, c("hetero", "two-way cluster-robust", "dyadic cluster-robust")), y = std.error)) +
  geom_boxplot() +
  xlab("Type of standard error") +
  geom_hline(yintercept = sd(df %>% dplyr::filter(term == "dx") %>% dplyr::pull(estimate)), linetype = "dashed") +
  theme_classic() +
  theme(text = element_text(size = 15))