頑健標準誤差 (robust standard error) を自力で計算する

RやStataといった統計ソフトでは、たいていの回帰モデルについて頑健な標準誤差 (robust standard error) を簡単に計算できる。
そのため、自力でrobust SEを求めなければいけないケースはほとんどない。
しかし、混合分布モデルなど複雑なモデルを用いる場合には、便利なコマンドが用意されていないこともある。
そこで、Rを用いて自分で計算できるよう備忘録として残しておこうと思う。

ちょっとだけ数式

Cameron and Trivedi (2005) の5.5.2節によると、パラメーター ${\hat{\mathbf{\theta}}}$ のrobust SEは

$\displaystyle N^{-1}\mathbf{A}^{-1}_0\mathbf{B}_0\mathbf{A}'^{-1}_0$

と表される。ここで $\mathbf{A}_0$ の推定値にはヘシアン

$\displaystyle \frac{\partial^2 Q(\mathbf{\theta})}{\partial \mathbf{\theta} \partial \mathbf{\theta}'}\Big|_{\hat{\mathbf{\theta}}}$

が、 $\mathbf{B}_0$ の推定値には各観測点の勾配の積の期待値

$\displaystyle \frac{1}{N} \sum_{i = 1}^N \frac{\partial q_i(\mathbf{\theta})}{\mathbf{\theta}}\Big|_{\hat{\mathbf{\theta}}} \frac{\partial q_i(\mathbf{\theta})}{\mathbf{\theta}'}\Big|_{\hat{\mathbf{\theta}}}$

が用いられる。 $Q(\mathbf{\theta})$ および $q_i(\mathbf{\theta})$ は対数尤度関数を表す。

Rコード

ここでは二値変数を従属変数とするロジット／プロビットモデルを例として、上の式を実装してみる。
まず、以下が二値ロジット／プロビットモデルの尤度関数である。

# log lilekihood function
llik.binary <- function(param, y, X, point = FALSE, cdf = plogis){
  tmp <- cdf(X %*% param)
  llik <- y * log(tmp) + (1 - y) * log(1 - tmp)
  if (!point){
    -1 * sum(llik)
  } else {
    llik
  }
}

ここでpoint = TRUEとすることで、pointwise log likelihoodが計算できる。
また、cdf = pnormとすればプロビットモデルの対数尤度関数となる。
次に、以下が前節の数式を実装したものである。

require(numDeriv) 
# compute the robust variance-covariance matrix
vcov.binary <- function(out, y, X, cluster = NULL){
  N <- length(y) # number of observations
  K <- length(out$par) # number of parameters
  M <- NULL # number of unique clusters
  u <- matrix(NA, nrow = N, ncol = K)
  # Caluculating u
  for (i in 1:N){
    u[i,] <- grad(func = llik.binary, x = out$par, y = y[i], X = X[i,], point = TRUE)
  }
  # for clustered SE 
  if (!is.null(cluster)){
    M <- length(unique(cluster)) 
    u.clust <- matrix(NA, nrow = M, ncol = K)
    # Then adjust u for clustering
    for (j in 1:K){
      u.clust[, j] <- tapply(u[, j], cluster, sum)
    }
  }
  # Then compute the robust vcov
  if (is.null(cluster)){
    vcov <- solve(out$hessian)%*%((N / (N - 1))*t(u)%*%(u))%*%solve(out$hessian)
  } else {
    vcov <- solve(out$hessian)%*%((M / (M - 1))*t(u.clust)%*%(u.clust))%*%solve(out$hessian)
  }
  vcov
}

勾配を計算するため、numDerivパッケージのgrad関数を用いる。
clusterにyと同じ長さのグループ変数を指定することで、clustered SEを計算することもできる。
また、Stataと近い結果を得るため、(N / (N - 1)) (もしくは(M / (M - 1))) をかけて調整している。

実データへの応用

ここではpsclパッケージ内にあるiraqVoteデータセットを用いる。
まず比較対象とするため、通常用いられる方法であるglm関数及びsandwichパッケージを利用し、回帰係数及びrobust SEを推定する。

> require(lmtest)
> require(pscl)
> require(sandwich)
> out1 <- glm(y ~ rep + gorevote, data = iraqVote, family = binomial, x = TRUE)
> coeftest(out1, vcovHC, type = "HC0")

z test of coefficients:

             Estimate Std. Error z value Pr(>|z|)   
(Intercept)  5.878586   2.714224  2.1658 0.030323 * 
repTRUE      3.018809   1.052731  2.8676 0.004136 **
gorevote    -0.113217   0.054421 -2.0804 0.037489 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

次に、前節で作った関数を用いて同じ分析を行う。

> out2 <- optim(rep(0, 3), llik.binary, y = iraqVote$y, X = out1$x, hessian = TRUE,
+               method = "BFGS")
# coefficients
> out2$par 
[1]  5.8774089  3.0188869 -0.1131937
# robust SE
> sqrt(diag(vcov.binary(out2, y = iraqVote$y, X = out1$x)))
[1] 2.72891380 1.05744123 0.05473097

glmを用いた分析とはやや数値が異なるが、これは最適化の誤差やアルゴリズムの差異に基づくものだろう。
最後に、この分析では必要ないが、州レベルでクラスター化した標準誤差も計算してみる。
まず、glmとsandwichを用いた結果が以下。

> coeftest(out1, vcovCL, cluster = ~ state.name)

z test of coefficients:

            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  5.87859    2.93595  2.0023 0.045255 * 
repTRUE      3.01881    1.06338  2.8389 0.004527 **
gorevote    -0.11322    0.06005 -1.8854 0.059381 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

次に、前節の関数を用いた結果が↓。

> sqrt(diag(vcov.binary(out2, y = iraqVote$y, X = out1$x, 
+                       cluster = iraqVote$state.name)))
[1] 2.93758429 1.06273594 0.06010004