複数時点の個体ごとの観測値であるパネルデータはよく見かけるものです。パネルデータには固定効果モデル（Fixed Effect Model）を用いることができ、欠落変数バイアスを抑制する手軽な方法として、定番の分析手法となっています。
1 モデル
単純な固定効果モデルは
\[ y_{it} = \alpha_i + \beta_0 + \beta_1 x_{1it} + \beta_2 x_{2it} + \cdots + \epsilon_{it} \]
と表すことができます。\(y\)が被説明変数、\(i\)は個体（もしくは個体が所属する群）を表す添字、\(t\)は時点を表す添字、\(\alpha\)は個体ごとの観測不能な効果（固定効果）、\(\beta\)は説明変数の係数、\(x\)は説明変数、\(\epsilon\)は誤差項です。
\[ y_{it} = \alpha_i + \beta_0 + \beta_1 x_{1it} + \beta_2 x_{2it} + \cdots + T_t + \epsilon_{it} \]
と、時間固定効果を入れることもできます。つまり、時点ごとの観測値のサブセットもクラスターと見なすことができます。
個体の種類で見たデータの広がりを\(N\)方向、時系列で見た広がりを\(T\)方向と呼びますが、パネルデータでは2種類の固定効果を考えることができ、また、同時に分析することができます。
式で見ると素朴ですが、実用的には欠落変数バイアスを大きく緩和できる有用な手法です。欠落変数が時間不変と見做せれば\(\alpha\)に含まれることになるので、その影響を無視できるからです。
変量効果モデルは、固定効果は説明変数と相関せず、誤差項の一部として扱えると仮定しますが、固定効果モデルは、固定効果が説明変数と相関していても使えます。現在の計量経済学では、固定効果が時間不変で説明変数と相関している場合を固定効果モデルと説明します1。
2 推定
勘のよい人はお気づきでしょうが、上のモデルはそのままは推定できません。固定効果をダミー変数として取り扱うdummy variable regressionを行うか、within変換をかけてからOLSで推定することになります。なお、OLSに持ち込めると言うことは、OLSの発展であるGLSやIVMにも応用できることになります。
2.1 データ生成プロセス
乱数から何種類かのパターンでデータセットを作成する関数を用意します。
genPanelData <- function(N = 26, T = 5, G = 3, coef = c(1, -2), sd = 5, type = "FE"){
numberToLetters <- function(ns, base = "A"){
sapply(ns, function(n){
s <- character(2)
i <- 1
repeat {
m <- n %% 26
n <- n %/% 26
s[i] <- rawToChar(as.raw(as.integer(charToRaw(base)) + m - 1*(i>1)))
i <- i + 1
if(0 >= n || i > length(s)) break
}
return(paste(s[(i-1):1], collapse = ""))
})
}
names <- numberToLetters(0:(N-1))
id <- rep(names, each = T)
id <- ordered(id, levels = names) # IDは順番をつけておく
gnames <- numberToLetters(0:(G-1), "a")
t <- ordered(rep(1:T, N))
a <- rep(runif(N, min = -5, max = 5), each = T)
x <- runif(N*T, min = 0, max = 5)
z <- runif(N*T, min = -5, max = 0)
e <- NULL
if("FE" == type){
Gamma <- matrix(0, N*T, N*T)
for(j in 1:ncol(Gamma)) Gamma[j, id == id[j]] <- 0.5
diag(Gamma) <- 1
e <- Gamma %*% rnorm(N*T, sd = sd)
} else if("FEGLS" == type) {
e <- rnorm(N*T, sd = sd)
} else if("FEWLS" == type){
s <- runif(N, min = 0.1*sd, max = 9.9*sd)
e <- rnorm(N*T, sd = rep(s, each = T))
} else if("FEIV" == type){
iv <- x
e <- rnorm(N*T, sd = sd)
x <- iv + e
y <- a + coef[1]*x + coef[2]*z + e
return(data.frame(id, t, y, x, z, iv))
} else stop("unknown type!")
y <- a + coef[1]*x + coef[2]*z + e
data.frame(id, t, y, x, z)
}2.2 データセット
単純な固定効果モデルのデータセットを作成します。
set.seed(311)
T <- 3
N <- 30
coef_DGP <- c(1, -2)
df01 <- genPanelData(N, T, coef = coef_DGP, type = "FEGLS")
summary(df01) id t y x z
A : 3 1:30 Min. :-5.130 Min. :0.006617 Min. :-4.99006
B : 3 2:30 1st Qu.: 2.692 1st Qu.:1.563161 1st Qu.:-3.33545
C : 3 3:30 Median : 8.317 Median :2.702058 Median :-2.30086
D : 3 Mean : 7.888 Mean :2.714324 Mean :-2.47702
E : 3 3rd Qu.:12.628 3rd Qu.:3.861687 3rd Qu.:-1.48571
F : 3 Max. :20.475 Max. :4.978560 Max. :-0.03297
(Other):72 [1] "ordered" "factor" [1] "ordered" "factor" idとtが要素型（かつ順序型）になっている事に注意してください。後でlm関数を使う都合でそうしています。
2.3 ダミー変数による回帰
固定効果モデルの推定ですが、個体や時点が少なければ、以下のように個体と時点のダミー変数を入れればOLSで推定できます（dummy variable regression）。
(Intercept) x z id.L id.Q id.C
1.18080338 0.88976588 -1.73281210 -4.12289392 -0.11895317 -0.88672573
id^4 id^5 id^6 id^7 id^8 id^9
0.15942619 -5.76354368 -1.14631505 2.64925743 -5.40546073 5.18340173
id^10 id^11 id^12 id^13 id^14 id^15
5.17869441 -2.94299434 -1.22683501 0.87675306 -0.38459255 1.36575971
id^16 id^17 id^18 id^19 id^20 id^21
-7.19791157 -1.41277450 0.41841239 1.64663280 -0.65201158 0.18079592
id^22 id^23 id^24 id^25 id^26 id^27
-0.35725866 0.12373388 -4.45372389 3.78269172 -4.79845193 -2.59899016
id^28 id^29 t.L t.Q
-7.35090753 2.78440688 0.05865376 -0.98844445 最近の計算機の処理速度と主記憶装置の容量を考えるとこれで間に合いそうですが、個体数が何百にもなると「直交多項式を十分正確に表現できません」とエラーになります。
2.4 within変換
個体や時点が多くても、安定して推定ができるシンプルな手法があります。within変換です。個体（もしくは時点）の各変数から、個体（〃）それぞれの各変数の平均値を引いたあとに、一般線形回帰モデル（OLS）などで推定します。
\[\begin{align} y_{it} - \bar{y_i} =& \alpha_i + \beta_0 + \beta_1 x_{1it} + \cdots + \epsilon_{it}\\ &- \bar{\alpha_i} - \bar{\beta_0} - \beta_1 \bar{x_{1it}} - \cdots\\ =& \beta_1 (x_{1it} - \bar{x_{1i}}) + \beta_2 (x_{2it} - \bar{x_{2i}}) \cdots + \epsilon_{it}\\ \end{align}\]
\(\alpha_i\)は各個体、\(\beta_0\)は全体で一定なので、この引き算で消えます。個体のwithin変換の後に、時点のwithin変換をかけることも可能です。形式的に切片項を計算したい場合は、全体の平均値を足します。
2.4.1 行列による計算
within変換をかける関数を用意しましょう。オプションIsStata = TRUEでwithin変換後に、変数それぞれの平均値をそれぞれ加算できるようにします。簡単な工夫ですが、これで切片項の推定ができるようになります。
within_transfer <- function(x, ..., IsStata = FALSE){
for(j in 1:...length()){
i <- ...elt(j)
m_i <- tapply(x, i, mean)
x <- x - m_i[i] + mean(x)*IsStata
}
c(x, use.names = FALSE)
}within変換後にOLSをかければ、固定効果モデルを推定できます。idを個体固定効果、tを時間固定効果として推定してみましょう。
IsStata <- FALSE
mf <- model.frame(y ~ x + z, df01)
X <- model.matrix(mf, mf)
y <- model.response(mf)
X <- apply(X, 2, within_transfer, df01$id, df01$t, IsStata = IsStata)
if(all(X[,1]==0)) X <- X[, -1] # IsStata = FALSEだと切片項はすべてゼロになるので除外
y <- within_transfer(y, df01$id, df01$t, IsStata = IsStata)
(beta.within <- solve(t(X) %*% X) %*% t(X) %*% y) [,1]
x 0.8897659
z -1.7328121この個体の種類ではDummy Variable
Regressionと同じ結果なので有り難みを感じませんが、N <- 200ぐらいで違いが出ると思います。IsStata = TRUEにすると切片項が計算されます。
2.4.2
plmパッケージによる計算
パッケージを使えば簡単に出ます。古くからある代表的なplmパッケージでは、以下のように推定できます。
library(plm)
r_plm <- plm(y ~ x + z, df01, effect = "twoway", model = "within", index = c("id", "t"))
summary(r_plm)Twoways effects Within Model
Call:
plm(formula = y ~ x + z, data = df01, effect = "twoway", model = "within",
index = c("id", "t"))
Balanced Panel: n = 30, T = 3, N = 90
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-9.31327 -3.31232 0.42447 3.25524 10.79302
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
x 0.88977 0.57523 1.5468 0.12754
z -1.73281 0.65003 -2.6657 0.01002 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 2229.3
Residual Sum of Squares: 1871.9
R-Squared: 0.16028
Adj. R-Squared: -0.33455
F-statistic: 5.3446 on 2 and 56 DF, p-value: 0.00751122.4.3
estimatrパッケージによる計算
比較すると最近出てきたestimatrパッケージを用いる場合は、
library(estimatr)
r_lm_robust <- lm_robust(y ~ x + z, data = df01, fixed_effects = ~ id + t, se_type = "classical")
summary(r_lm_robust)
Call:
lm_robust(formula = y ~ x + z, data = df01, fixed_effects = ~id +
t, se_type = "classical")
Standard error type: classical
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
x 0.8898 0.5752 1.547 0.12754 -0.2626 2.0421 56
z -1.7328 0.6500 -2.666 0.01002 -3.0350 -0.4306 56
Multiple R-squared: 0.4643 , Adjusted R-squared: 0.1486
Multiple R-squared (proj. model): 0.1603 , Adjusted R-squared (proj. model): -0.3345
F-statistic (proj. model): 5.345 on 2 and 56 DF, p-value: 0.007511とします。
2.5 係数の検定と信頼区間
基本的にはOLSなので、係数の検定と信頼区間も同様に出ます。また、plm::plmもestimatr::r_lm_robustもvcovやconfintをサポートし、lmと同様の使い勝手になっています。
3 モデル選択
現在では推奨されていませんが、統計的仮説検定でモデル選択を行うこともできます。
Breusch-Pagan検定で統計的に有意になればプール（非パネル）よりも変量効果モデルが推奨され、Hausman検定で統計的に有意になれば変量効果モデルよりも固定効果モデルが推奨されることになりますが、
- モデル選択をした後にパラメーターを推定した場合は係数のp値に修正を入れる必要がある
- そもそもこれらの検定によるモデル選択の正当性が微妙と言う報告もあります。
- 変量効果モデルが真のときも固定効果モデルは（有効で無くても）一致推定量が得られるので、固定効果モデルを選択しておけば大間違いはないが、固定効果モデルが真のときに変量効果モデルで推定すると一致推定量が得られない
ので、欠落変数バイアスの緩和目的であれば固定効果モデルで推定し、モデル選択に使う検定は分析の頑強性を示す補足として示すのが良さそうです。
教科書的で昔の統計解析ソフトでも出力された検定3種類をplmを使って行います。
# Random効果モデルで推定する
r_plm_re <- plm(y ~ x + z, df01, effect = "individual", model = "random", index = c("id", "t"))
# Breusch-Pagan-test (random vs pool; H₀:pool)
pcdtest(r_plm_re, test = "lm")
Breusch-Pagan LM test for cross-sectional dependence in panels
data: y ~ x + z
chisq = 641.26, df = 435, p-value = 4.11e-10
alternative hypothesis: cross-sectional dependence# 固定効果モデルを推定する
r_plm_fe <- plm(y ~ x + z, df01, effect = "individual", model = "within", index = c("id", "t"))
# プールモデルで推定する
r_plm_po <- plm(y ~ x + z, df01, effect = "individual", model = "pool", index = c("id", "t"))
# F test of stability (or Chow test) (pool vs within;H₀:pool)
pooltest(r_plm_po, r_plm_fe)
F statistic
data: y ~ x + z
F = 1.0452, df1 = 29, df2 = 58, p-value = 0.4314
alternative hypothesis: unstability
Hausman Test
data: y ~ x + z
chisq = 0.59356, df = 2, p-value = 0.7432
alternative hypothesis: one model is inconsistentBreusch-Pagan検定は誤差項とindexの関係を見ており、Chow検定は固定効果のばらつきを見ており、Hausman検定は固定効果を入れることで係数がどれぐらい変化するかを評価しています。パネルデータモデルに関する検定は他にもあり、plmパッケージでは他にWooldridge’s
Test（誤差項とクラスターの相関）とBera, Sosa-Escudero and Yoon
Locally–Robust Lagrange Multiplier Testsをサポートしています。
4 クラスター頑強標準誤差
ある個体のある時点に何かハプニングがあれば、あとの時点にも影響が見られる場合があります。ある個体にハプニングがあった場合、その個体が属するグループの他の個体にも影響が見られる場合があります。同一個体のデータや同一グループのデータは、誤差項がi.i.dとは限りません。不均一分散が生じます。
不均一分散があると、標準誤差を過小評価することになります。統計的仮説検定で偽陽性が出やすくなり、信頼区間の幅が短くなってしまいます。パネルデータ分析の場合は、個体やグループなどのクラスターを内包している事が一般的であり、不均一分散の影響が強い可能性があります。このため、近年では固定効果モデルによる推定でも、クラスター頑強標準誤差を計算することが一般的です。
この計算が意外に煩雑で、一時期は頻出の質問になっていました。現在もRのパッケージでのサポートはまちまちです。確認した限り、estimatrとdfadjustが代表的な手法であるCR₂とCR₃2をオプションで指定できましたが、plmとfixestはそうではありませんでした。
4.1 推定
個体固定効果と時間固定効果があるone-way clusteringのCR₂を例に推定してみます。
4.1.1 行列による計算
クラスター頑強標準誤差の計算をしてくれる便利なビルトイン関数はありません。
vcov_cr2 <- function(X, residuals, id = NULL, t = NULL, mode = "estimatr", ncol = NULL){
if(is.null(ncol)) ncol <- ncol(X)
nobs <- length(residuals)
B <- matrix(0, ncol(X), ncol(X))
C <- chol(crossprod(X)) # t(X)%*%Xは対称行列になるのでコレスキー分解を用いる
Pss <- X %*% backsolve(C, forwardsolve(t(C), t(X))) # = X %*% solve(t(X) %*% X) %*% t(X)
if("estimatr" == mode){
# 異なるクラスター間の共分散の期待値はゼロなので、該当箇所にゼロを入れる
for(j in 1:nobs){
if(!is.null(id)){
if(is.null(t)){
# one-way clustering
Pss[, j] <- Pss[, j] * (id[j] == id)
} else {
# two-way clustering
Pss[, j] <- Pss[, j] * (id[j] == id | t[j] == t)
}
} else if(!is.null(t)) {
# one-way clustering
Pss[, j] <- Pss[, j] * (t[j] == t)
} else stop("set either i or t or both of them!")
}
# I - Pssも対称行列になるので、逆行列が転置行列になる直交行列で対角化される
eigen <- eigen(diag(nobs) - Pss, symmetric = TRUE)
inv_sqrt_I_Pss <- eigen$vector %*% diag(sqrt(pmax(0, 1/eigen$values))) %*% t(eigen$vector) # sqrt()を無くすとCR₃になる
m_residuals <- inv_sqrt_I_Pss %*% residuals
} else if("plm" == mode) {
# 論文の説明を読むとI(n) - Pssの逆行列の1/2乗を乗じそうなのだが、対角成分だけを使う
m_residuals <- residuals / sqrt(1 - diag(Pss))
} else stop("set either \"plm\" or \"estimatr\"!")
omega <- matrix(0, nobs, nobs)
for(j in 1:nobs){
if(!is.null(id)){
if(is.null(t)){
# one-way clustering
omega[,j] <- m_residuals[j] * m_residuals * (id[j] == id)
} else {
# two-way clustering
omega[,j] <- m_residuals[j] * m_residuals * (id[j] == id | t[j] == t)
}
} else if(!is.null(t)) {
# one-way clustering
omega[,j] <- m_residuals[j] * m_residuals * (t[j] == t)
} else stop("set either i or t or both of them!")
}
B <- crossprod(X, omega) %*% X # = t(X) %*% omega %*% X
r <- backsolve(C, forwardsolve(t(C), t(backsolve(C, forwardsolve(t(C), t(B)))))) # = solve(t(X) %*% X) %*% B %*% solve(t(X) %*% X)
r <- r[1:ncol, 1:ncol]
colnames(r) <- rownames(r) <- colnames(X)[1:ncol]
r
}
eX <- cbind(X, model.matrix(~ 0 + t, df01)) # witin変換で消したダミー項を追加する
residuals <- y - X %*% beta.within
(sqrt(diag(vcov_cr2(eX, residuals, df01$id, ncol = ncol(X))))) x z
0.5236931 0.6412517 witin変換で消したダミー項を追加する処理が入っていますが、ダミー項の係数にも標準誤差はあり、それが他の係数の標準誤差に影響するからです。この処理がないと、dummy variable regressionの標準誤差と差異が出ます。
idによるダミー項を入れていないのは、データ的に推定不可能でありNAが出てエラーになるからです。idダミーの係数の標準誤差はゼロに近いと考えられるため、入れなくてもその他の係数の標準誤差に影響はないでしょう。
クラスター頑強標準誤差は厳密にはt分布に従わなくなっているので、自由度を調整する必要があります。その計算も行いましょう。
df_cr2 <- function(X, residuals, id = NULL, t = NULL, IK = FALSE, ncol = NULL){
if(is.null(ncol)) ncol <- ncol(X)
nobs <- length(residuals)
C <- chol(crossprod(X)) # t(X)%*%Xは対称行列になるのでコレスキー分解を用いる
Pss <- X %*% backsolve(C, forwardsolve(t(C), t(X))) # = X %*% solve(t(X) %*% X) %*% t(X)
I_Pss <- diag(nobs) - Pss
clusters <- NULL
if(!is.null(id)){
if(is.null(t)){
# one-way clustering
clusters <- expand.grid(levels(id), NA, stringsAsFactors = FALSE)
} else {
# two-way clustering
clusters <- expand.grid(levels(id), levels(t), stringsAsFactors = FALSE)
}
} else if(!is.null(t)) {
# one-way clustering
clusters <- expand.grid(NA, levels(t))
} else stop("set either i or t or both of them!")
df <- numeric(ncol)
names(df) <- colnames(X)[1:ncol]
Omega <- diag(nobs)
# Imbens-Kolesár (2016)'s method, employing Moulton's Omega
if(IK){
m <- 0
rho <- 0
sigma <- 0
for(i in 1:nrow(clusters)){
c_id <- clusters[i, 1]
c_t <- clusters[i, 2]
if(!is.na(c_id)){
if(!is.na(c_t)) index <- c_id == id | c_t == t
else index <- c_id == id
} else index <- c_t == t
m <- m + sum(index)^2 # TRUEは1, FALSEは0として合計される
rho <- rho + c(crossprod(residuals[index]))
}
rho <- (rho - residuals^2)/(m - nobs)
sigma <- mean(residuals^2)
for(j in 1:nobs){
if(!is.null(id)){
if(is.null(t)) Omega[, j] <- rho * (id[j] == id)
else Omega[, j] <- rho * (id[j] == id | t[j] == t)
} else if(!is.null(t)) Omega[, j] <- rho * (t[j] == t)
else stop("set either i or t or both of them!")
}
diag(Omega) <- sigma
}
for(k in 1:ncol){
e <- numeric(ncol(X))
e[k] <- 1
G <- matrix(NA, nobs, length(levels(id)))
j <- 1
for(i in 1:nrow(clusters)){
c_id <- clusters[i, 1]
c_t <- clusters[i, 2]
if(!is.na(c_id)){
if(!is.na(c_t)) index <- c_id == id | c_t == t
else index <- c_id == id
} else index <- c_t == t
X_s <- X[index, ]
I_Pss_s <- I_Pss[index, index] # = diag(nrow(X_s)) - (Pss_ss <- X_s %*% backsolve(C, forwardsolve(t(C), t(X_s))))
eigen <- eigen(I_Pss_s, symmetric = TRUE)
inv_sqrt_I_Pss_s <- eigen$vector %*% diag(sqrt(pmax(0, 1/eigen$values))) %*% t(eigen$vector) # sqrt()を無くすとCR₃用の自由度計算になる
G[, j] <- I_Pss[, index] %*% inv_sqrt_I_Pss_s %*% X_s %*% backsolve(C, forwardsolve(t(C), e))
j <- j + 1
}
GG <- crossprod(G, Omega) %*% G # = t(G) %*% Omega %*% G
# トレースの和は、固有値の和に一致する
# トレースの二乗和は、crossproductした行列の固有値の和に一致する
df[k] <- (sum(diag(GG))^2)/sum(diag(crossprod(GG))) # crossprod(GG) = t(GG) %*% GG
}
df
}
df_cr2(eX, residuals, df01$id, ncol = ncol(X)) x z
16.76919 16.99286 4.1.2
plmパッケージによる計算
本稿執筆時点ではplmパッケージでCR₂（とCR₃）の計算できません。不均一分散に頑強な標準誤差を拡張してクラスター頑強標準誤差を計算してくれるのですが、Imbens and Kolesár
(2016)で紹介されるCR₂とは異なっています。
plmでも不均一分散に頑強な標準誤差HC₀（ホワイト頑強標準誤差）、HC₁（Stata互換）、HC₂、HC₃…を使うことができるのですが、CR₂、CR₃をオプションで指定できません。plmはHC₂と同時にclustersオプションをつけると系列相関を仮定した計算にかけることができますが、残差をprojection matrixの対角成分であるlevarage値で補正してきます。惜しいところで、projection matrixの同一クラスター内の共分散を非ゼロとして残差の補正をしないとCR₂になりません。またplmは、自由度調整を行わないので小標本でバイアスが大きくなります。
Estimate Std. Error t-value Pr(>|t|)
x 0.8897659 0.5072498 1.754098 0.084884697
z -1.7328121 0.6213552 -2.788763 0.007215839クラスター頑強標準誤差が必要な場合は、plmは有力な選択肢ではなくなります。
4.1.3
estimatrパッケージによる計算
estimatrパッケージはCR₂をサポートしています。推定時にclusters = id, se_type = "CR2"をつけるだけです。
r_lm_robust <- lm_robust(y ~ x + z, data = df01, fixed_effects = ~ id + t, clusters = factor(as.character(id)), se_type = "CR2")
summary(r_lm_robust)
Call:
lm_robust(formula = y ~ x + z, data = df01, clusters = factor(as.character(id)),
fixed_effects = ~id + t, se_type = "CR2")
Standard error type: CR2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
x 0.8898 0.5237 1.699 0.10779 -0.2163 1.9958 16.77
z -1.7328 0.6413 -2.702 0.01511 -3.0858 -0.3798 16.99
Multiple R-squared: 0.4643 , Adjusted R-squared: 0.1486
Multiple R-squared (proj. model): 0.1603 , Adjusted R-squared (proj. model): -0.3345
F-statistic (proj. model): 8.654 on 2 and 29 DF, p-value: 0.001129現時点のバージョンではclustersの引数にorder型をとれないので唯のfactor型に変換しています。orderはfactorとしても振る舞えるので、エラーチェックが不適切ですね。
4.2 クラスターの選択
上の例では固定効果とクラスターを一致させていましたが、一致させる必要はありません。例えば、固定効果を市町村ごとに入れて、都道府県をクラスターとすることもできます。本稿は固定効果モデルの話ですが、固定効果モデルですらなくてもクラスター頑強標準誤差は計算できます。
5 固定効果-一般化線形回帰（FE-GLS）
クラスタ頑強標準誤差を用いずに、固定効果の入った一般化線形回帰（FE-GLS）を推定する事も出来ます。GLSなので、誤差項に対する想定が正しければ、より効率的な推定になります。誤差項に関する想定は、ランダム化効果モデルと同様です。Wooldridge (2002)の§10.5で紹介されていたのですが、その説明とちょっと違う方法3がplmパッケージで提供されているので試してみましょう。
5.1 データセット
FE-GLS向きのデータセットを生成します。
5.2
plmパッケージによる推定
FE-GLSの推定を行い、FEによる推定と各種クラスター頑強標準誤差との比較を行います。
library(plm)
frml <- y ~ x + z
r_plm <- plm(frml, data = df_fegls, effect = "individual", model = "within", index = "id")
r_pggls <- pggls(frml, data = df_fegls, effect = "individual", model = "within", index = "id")
(CMP.beta <- matrix(c(coef_DGP, coef(r_plm), coef(r_pggls)), 2, 3, dimnames = list(names(coef(r_plm)), c("DGP/True", "within", "FE-GLS")))) DGP/True within FE-GLS
x 1 1.591416 1.401775
z -2 -2.675240 -2.390549(CMP.se <- matrix(c(
sqrt(diag(vcov(r_plm))),
sqrt(diag(vcovHC(r_plm, type = 'HC0', cluster = "group"))), # cluster = "group"とtype = 'HC0'で、CR₀になる
sqrt(diag(vcovHC(r_plm, type = 'sss', cluster = "group"))),
sqrt(diag(vcovHC(r_plm, type = 'HC2', cluster = "group"))),
sqrt(diag(vcovHC(r_plm, type = 'HC3', cluster = "group"))),
sqrt(diag(vcov(r_pggls)))
), 2, 6, dimnames = list(names(coef(r_plm)), c("FE(within)", "HC0", "Stata-Like", "HC2", "HC3", "FEGLS")))) FE(within) HC0 Stata-Like HC2 HC3 FEGLS
x 0.3670927 0.4810662 0.4913590 0.4876506 0.4943582 0.3330074
z 0.3778542 0.3306870 0.3377623 0.3358966 0.3412289 0.3237355微妙に係数の推定量がDGPの係数に近づき、標準誤差が小さくなっています。ただし、set.seedをコメントアウトして何回も試してもらうとわかりますが、FEの方がよい推定量になっている事も多いです。系列相関の影響は固定効果の推定量に吸収されているわけで、改善余地はそんなに無いのでしょう。
5.3 行列計算による推定
plmパッケージの挙動を再現してみましょう。ただし、plm::pgglsでは疎行列をサポートするbdsmatrixパッケージを利用していましたが、密行列としてコレスキー分解で凌ぎます。単純に逆行列を求めないのは、特異値になるからです。
df_sorted <- df_fegls[with(df_fegls, order(id, t)), ]
# Wooldridge (2002) §10.5の説明通りであれば以下だが、pgglsの実装はそうではない
# df_sorted <- subset(df_sorted, t != T)
T_sorted <- length(levels(factor(df_sorted$t)))
mf <- model.frame(update(frml, . ~ . + 0), df_sorted)
X <- model.matrix(mf, mf)
y <- model.response(mf)
X <- apply(X, 2, within_transfer, df_sorted$id)
y <- within_transfer(y, df_sorted$id)
beta_within <- solve(t(X) %*% X) %*% t(X) %*% y
residuals <- y - X %*% beta_within
# idごとにデータを切り出し、tの分散共分散行列をつくる
M <- matrix(0.0, T_sorted, T_sorted)
for(E in tapply(residuals, df_sorted$id, \(x) x %*% t(x))) M <- M + E
M <- M / N
# Omega <- diag(N) %x% Mと描きたいが、それでは数値が不安定で逆行列の乗算ができない
# コレスキー分解を用いる
C <- chol(M) # M == t(C) %*% C
# クロネッカー積で拡大すれば、Omega = diag(N) %x% Mのコレスキー分解
Cx <- diag(N) %x% C
# beta_fegls <- solve(t(X) %*% solve(Omega) %*% X) %*% t(X) %*% solve(Omega) %*% y を行う
A <- t(X) %*% backsolve(Cx, forwardsolve(t(Cx), X)) # = t(X) %*% solve(Omega) %*% y
B <- t(X) %*% backsolve(Cx, forwardsolve(t(Cx), y)) # = t(X) %*% solve(Omega) %*% X
(beta_fegls <- as.numeric(solve(A, B)))[1] 1.401775 -2.390549[1] 0.3330074 0.3237355同じ値になりました。数値的に不安定な面があり、推定が回らないと言う事はありそうです。
6 固定効果-ウェイト付き線形回帰（FE-WLS）
系列相関ではなくクラスター（もしくは同一個体）ごとに誤差項の大きさが異なることを仮定したWLSとの組み合わせも試してみましょう。
6.1 データセット
クラスターごとに誤差項の大きさが異なるデータセットを作成します。
6.2 行列計算による推定
パッケージでサポートされていないようなので、行列で計算します。
frml <- y ~ x + z
# applyはIDで整列した順番に値を返すのでソートしておく
df_fewls_sorted <- df_fewls[with(df_fewls, order(id, t)), ]
# データフレームを行列に展開
mf <- model.frame(frml, df_fewls_sorted)
X <- model.matrix(mf, mf)
y <- model.response(mf)
# within変換をかける
X <- apply(X, 2, within_transfer, df_fewls_sorted$id, IsStata = T)
y <- within_transfer(y, df_fewls_sorted$id, IsStata = T)
# Stata互換切片項あり固定効果モデル
coef_within_stata <- solve(t(X) %*% X) %*% t(X) %*% y
residuals <- X %*% coef_within_stata - y
# 切片項なし固定効果モデル
X0 <- apply(model.matrix(mf, mf)[,-1], 2, within_transfer, df_fewls_sorted$id)
y0 <- within_transfer(y, df_fewls_sorted$id)
coef_within <- solve(t(X0) %*% X0) %*% t(X0) %*% y0
# クラスター間で誤差の分散が異なる一方、クラスター内の系列相関が無いモデル
sigma2_i <- tapply(residuals, df_fewls_sorted$id, \(x) sum(x^2)/length(x))
W <- diag(1/sqrt(sigma2_i[df_fewls_sorted$id]))
Omega <- diag(1/sigma2_i[df_fewls_sorted$id])
coef_fewls <- solve(t(X) %*% Omega %*% X) %*% t(X) %*% Omega %*% y
df <- nrow(X) - ncol(X) - N + 1 # 自由度をdummy variable estimatorにあわせる
sigma_fewls <- sqrt(sum((W %*% (y - X %*% coef_fewls))^2) / df)
vcov_fewls <- sigma_fewls^2 * solve(t(X) %*% Omega %*% X)
se_fewls <- sqrt(diag(vcov_fewls))
print(CMP <- matrix(c(NA, coef_DGP[1], coef_DGP[2], NA, coef_within, coef_within_stata, coef_fewls), 3, 4,
dimnames = list(c("Const.", "x", "z"), c("DGP/True", "within", "Stata-Like", "FE-WLS")))) DGP/True within Stata-Like FE-WLS
Const. NA NA -7.590897 -4.948839
x 1 1.740772 1.740772 1.118976
z -2 -3.653774 -3.653774 -3.184559係数が真の値に近づいていますね。set.seedを消して何回か試してみましたが、誤差項に関する仮定が正しいだけに、こちらは概ね、推定結果を改善できるようです。
クラスター（もしくは同一個体）ごとに誤差の分散が異なる不均一性があると仮定していますが、他の変数に比例するように仮定してWSLSをかけるのも、WとOmegaの計算を教科書的なFGLSと同様にすればよいので難しくはないです。
7 固定効果-操作変数法（FE-IV）
固定効果モデルはOLSの応用なので、操作変数法もほぼ同様に用いることができます。
7.1 データセット
説明変数の一部に内生性があり、操作変数があるデータセットをつくります。
set.seed(2451)
T <- 5
N <- 100
coef_DGP <- c(1, -2)
df_feiv <- genPanelData(N, T, coef = coef_DGP, sd = 3, type = "FEIV")ただのFEで推定すると、xの係数に大きなバイアスが入ります。
Call:
lm_robust(formula = y ~ x + z, data = df_feiv, fixed_effects = ~id)
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
x 1.815 0.01823 99.57 2.135e-283 1.780 1.851 398
z -2.044 0.04314 -47.37 1.119e-165 -2.129 -1.959 398
Multiple R-squared: 0.9763 , Adjusted R-squared: 0.9702
Multiple R-squared (proj. model): 0.9657 , Adjusted R-squared (proj. model): 0.957
F-statistic (proj. model): 6023 on 2 and 398 DF, p-value: < 2.2e-167.2
plmパッケージによる推定
plmパッケージはモデル式に外生変数をリストすれば、操作変数法で推定してくれます。
library(plm)
summary(r_plm_feiv <- plm(y ~ x + z | iv + z, df_feiv, effect = "individual", model = "within", index = c("id", "t")))Oneway (individual) effect Within Model
Instrumental variable estimation
Call:
plm(formula = y ~ x + z | iv + z, data = df_feiv, effect = "individual",
model = "within", index = c("id", "t"))
Balanced Panel: n = 100, T = 5, N = 500
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-7.755946 -1.618664 -0.063813 1.666907 7.256245
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
x 1.02715 0.10080 10.191 < 2.2e-16 ***
z -2.05881 0.10189 -20.207 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 19025
Residual Sum of Squares: 3438.3
R-Squared: 0.90247
Adj. R-Squared: 0.87772
Chisq: 520.154 on 2 DF, p-value: < 2.22e-16plmパッケージではCR₂やCR₃は出ませんが、Stata互換クラスター頑強標準誤差（CR₁）は出せます。
Oneway (individual) effect Within Model
Instrumental variable estimation
Note: Coefficient variance-covariance matrix supplied: vcovHC(r_plm, type = "sss")
Call:
plm(formula = y ~ x + z | iv + z, data = df_feiv, effect = "individual",
model = "within", index = c("id", "t"))
Balanced Panel: n = 100, T = 5, N = 500
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-7.755946 -1.618664 -0.063813 1.666907 7.256245
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
x 1.02715 0.49136 2.0904 0.03658 *
z -2.05881 0.33776 -6.0954 1.091e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 19025
Residual Sum of Squares: 3438.3
R-Squared: 0.90247
Adj. R-Squared: 0.87772
Chisq: 39.5684 on 2 DF, p-value: 2.5575e-097.3
estimatrパッケージによる推定
estimatrパッケージはiv_robust関数でモデル式に外生変数をリストすれば、操作変数法で推定してくれます。
library(estimatr)
summary(r_iv_robust <- iv_robust(y ~ x + z | iv + z, df_feiv, fixed_effect = ~id, se_type = "classical"))
Call:
iv_robust(formula = y ~ x + z | iv + z, data = df_feiv, fixed_effects = ~id,
se_type = "classical")
Standard error type: classical
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
x 1.027 0.1008 10.19 7.997e-22 0.829 1.225 398
z -2.059 0.1019 -20.21 5.386e-63 -2.259 -1.859 398
Multiple R-squared: 0.8748 , Adjusted R-squared: 0.843
Multiple R-squared (proj. model): 0.8193 , Adjusted R-squared (proj. model): 0.7734
F-statistic (proj. model): 260.1 on 2 and 398 DF, p-value: < 2.2e-16クラスター頑強標準誤差が必要な場合は、以下のようにします。
summary(r_iv_robust <- iv_robust(y ~ x + z | iv + z, df_feiv, fixed_effect = ~id, clusters = factor(as.character(id)), se_type = "CR2"))
Call:
iv_robust(formula = y ~ x + z | iv + z, data = df_feiv, clusters = factor(as.character(id)),
fixed_effects = ~id, se_type = "CR2")
Standard error type: CR2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
x 1.027 0.10253 10.02 2.107e-15 0.8228 1.231 73.79
z -2.059 0.09409 -21.88 2.834e-35 -2.2461 -1.872 79.01
Multiple R-squared: 0.8748 , Adjusted R-squared: 0.843
Multiple R-squared (proj. model): 0.8193 , Adjusted R-squared (proj. model): 0.7734
F-statistic (proj. model): 273.8 on 2 and 99 DF, p-value: < 2.2e-16なおse_type = "CR2"はclusters指定時のデフォルトです。
7.4 行列計算による推定
係数の推定は、within変換をかけたあとに操作変数法を行うだけです。
mf <- model.frame(y ~ x + z + 0, df_feiv)
Z <- apply(model.matrix(x ~ iv + z + 0, df_feiv), 2, within_transfer, df_feiv$id)
X <- apply(model.matrix(mf, mf), 2, within_transfer, df_feiv$id)
y <- within_transfer(model.response(mf), df_feiv$id)
(beta.feiv <- solve( t(Z) %*% X ) %*% t(Z) %*% y) [,1]
x 1.027155
z -2.058807パッケージと同じ係数になりました。
# Xの予測値を計算
pX <- Z %*% solve(t(Z) %*% Z) %*% t(Z) %*% X
# 操作変数法の場合は、予測値でクラスター頑強標準誤差とその自由度を計算する
vcov.feiv <- vcov_cr2(pX, y - X %*% beta.feiv, id = df_feiv$id)
(sqrt(diag(vcov.feiv))) x z
0.1025346 0.0940882 x z
73.78922 79.00733 estimatrと同じCR₂と補正自由度になっています。なお、時間固定効果を入れていないので、時点ダミーをpXに追加していません。
7.5 検定
標準誤差が出るので、OLSと同じ要領で検定が行えます。
plmパッケージのヴィネットに詳しいです。↩︎クラスター頑強標準誤差は、代表的なものでCR₀，CR₁，CR₂，CR₃があります。末石 (2025)「ミクロ計量分析講義資料 スライド2: 不均一分散とクラスターに頑健な統計的推測」によると、CR₀は漸近的に一致するものの、クラスターの数が少ない有限標本では十分に漸近せず、CR₁はよく使われてきた統計解析ソフトStataのデフォルトですが、統計学的に十分な正当化はないそうです。CR₂は均一分散の下で不偏性がある一方、不均一分散の下での性質は必ずしもよくないそうです。CR₃は過大推定バイアスがあり、偽陽性を出しづらいと言う意味で保守的な値になります。↩︎
新しい版のWooldridge (2010)の方は確認していないです。↩︎