Zenn
🙆‍♀️

Modified Poissonを理解する - PART2 実装編

2024/05/06に公開
R
#
統計
#
修正ポアソン
#
サンドイッチ分散
tech

ゴール

  • スクラッチ実装を通してModified poisson regression modelの理解を深める

やること

  • 通常のポアソン回帰における分散共分散行列の出力をRでスクラッチ実装する
  • {sandwich}の出力をRでスクラッチ実装する
  • 擬似ポアソン回帰や修正ポアソン回帰のパッケージとの比較　など

関連記事

Modified Poissonを理解する - PART1 理論編にGLMのおさらいからModified poisson regressionモデルについての理論部分を整理したのでご参照されたい。

vcov()の実装 - nが大きい場合

vcov出力の確認

まず、Rにおいて分散共分散行列といえばvcov()関数である。これを自前の計算と一致させることを目指す。

データを適当に生成してglmでポアソン回帰を実施し、vcov()の出力を確認する。ここではn=500としておく。このデータは、ロジスティック回帰を実行したくなるようなデータでポアソン回帰を使うというModified poissonのテーマを考えて、二項分布から生成している。

datagen <-
  function(n) {
    N <- n # サンプルサイズ

    v1 <- rnorm(N, mean = 60, 10)
    v2 <- rbinom(N, size = 1, prob = 0.5)
    v3 <- rbinom(N, size = 1, prob = 0.5)
    v4 <- rbinom(N, size = 1, prob = 0.5)
    X <- do.call(cbind, list(1,v1,v2,v3,v4))
    colnames(X) <- c('const','v1', 'v2', 'v3', 'v4')
    P <- ncol(X) # 説明変数の数（定数項含む）
    
    true_beta <- c(-2.2, 0, 0.8, 0, 0) # betaの真値
    Y <- unlist(lapply(exp(X%*%true_beta), rbinom, n = 1, size = 1))
    
    X[,2] <- scale(X[,2])　# すごく大きい正規分布の確率変数だけ標準化
    mydata <- data.frame(X[,-1],Y)
    
    return(list(X=X,mydata=mydata)) # 行列計算用にデザイン行列だけ取り出せるように
  }

# generate data
n <- 500
d <- datagen(n)
mydata <- d$mydata
X <- d$X


# Poisson regression model
myposfit <- glm(Y~v1+v2+v3+v4, data = mydata, family = poisson())

> (vcov_correct <- vcov(myposfit))
              (Intercept)            v1            v2            v3            v4
(Intercept)  0.0640516310 -0.0006341787 -0.0422844243 -0.0213589116 -0.0220302090
v1          -0.0006341787  0.0105467470  0.0007115129 -0.0009420254  0.0009834316
v2          -0.0422844243  0.0007115129  0.0578881649 -0.0011632617 -0.0014411411
v3          -0.0213589116 -0.0009420254 -0.0011632617  0.0433162535  0.0025330402
v4          -0.0220302090  0.0009834316 -0.0014411411  0.0025330402  0.0433022202

スクラッチ

ポアソン回帰のw_{ii}として、glm()の最後のイテレーションのworking weightsを用いて自前計算する。

# variance of poisson equals to its mean
mu <- myposfit$weights # weights in poisson is exp(Xb) = var = mu
W <- diag(mu)

# J_n^-1 = X'WX 
J_n_inv <- solve(crossprod(X*mu,X)) # crossprod(X*mu,X) = t(X)%*%W%*%X

colnames(J_n_inv) <- colnames(vcov_correct)
rownames(J_n_inv) <- rownames(vcov_correct)

一致した。

> all.equal(vcov_correct, J_n_inv)
[1] TRUE

working weightsを使わずにmyposfit$fitted.valuesを用いると、誤差レベルではあるが結果はずれる。

## using fitted.values give slightly different vcov
mu <- myposfit$fitted.values 

W <- diag(mu)
# J_n^-1 = X'WX 
J_n_inv <- solve(crossprod(X*mu,X))

colnames(J_n_inv) <- colnames(vcov_correct)
rownames(J_n_inv) <- rownames(vcov_correct)

> all.equal(vcov_correct, J_n_inv)
[1] "Mean relative difference: 6.463396e-08"

sandwich()の実装 - nが大きい場合

sandwich()の出力の確認

sandwich::sandwich()の結果を確認する。

library(sandwich)

> (sandwich_correct <- sandwich(myposfit))
              (Intercept)            v1            v2           v3            v4
(Intercept)  0.0585916806  0.0005937901 -0.0406663845 -0.018714125 -0.0182578852
v1           0.0005937901  0.0075867413  0.0008468263 -0.001654921 -0.0010496916
v2          -0.0406663845  0.0008468263  0.0494447742  0.002502528  0.0009474114
v3          -0.0187141252 -0.0016549206  0.0025025284  0.033609797  0.0013451451
v4          -0.0182578852 -0.0010496916  0.0009474114  0.001345145  0.0331144430

スクラッチ

Bread部分（理論編の表記でいうとA）については、上のvcovと同じ中身であるので（スコア関数の直積をE[(y_i - \mu_i)^2]=Var(y_i)と仮定した通常のI_nとヘシアンは同じ形）、新たに計算することはしない。

Sandwich estimatorで新たに登場するmeat部分（理論編の表記でいうとB）を計算する。残差の二乗を重みとしたX^TWXだったので、\sqrt{W}Xの直積として考えてみる。

# the meat part of sandwich 
resid <- mydata$Y - mu # Y - exp(X \beta) at the last iteration of IWLS in glm.fit

observed_score <- X*resid
V_n <- crossprod(observed_score,observed_score)

# Sandwich estimator 
mysandwich <- J_n_inv%*%V_n%*%J_n_inv

10^{-7}のオーダーでの違いということなので、概ね一致したと思って良いだろう。

> all.equal(mysandwich, sandwich_correct)
[1] "Mean relative difference: 1.461576e-07"

通常のフィッシャー逆行列の推定値を用いた場合の共分散行列とsandwich推定値との差分(%)を見てみる。全て差分が負なので、補正前の方が標準誤差が大きく出ていたようだ。また、一部-100%を超えるような要素もあり、かなりの違いが%ベースでは出ていることが分かる。

> round((sandwich_correct - vcov_correct)/vcov_correct * 100,2)
            (Intercept)      v1      v2      v3      v4
(Intercept)       -8.52 -193.63   -3.83  -12.38  -17.12
v1              -193.63  -28.07   19.02   75.68 -206.74
v2                -3.83   19.02  -14.59 -315.13 -165.74
v3               -12.38   75.68 -315.13  -22.41  -46.90
v4               -17.12 -206.74 -165.74  -46.90  -23.53

vcov()の実装 - nが小さい場合

今度はn=50という小サンプルのケースでやってみる。
大体の場合うまくいくのだが、何度かシードを変えてみたところsandwichがうまくいかないことがあったので、そのときのデータを用いた結果を記載する。

set.seed(11111111) 

# generate data
n <- 50
d <- datagen(n)
mydata <- d$mydata
X <- d$X

\eta_iは下記のような形で10-20%くらいのリスクを持っているデータではあるが、このシードでは実際のアウトカム割合が6%しかない。

> exp(X%*%true_beta)
           [,1]
 [1,] 0.2465970
 [2,] 0.2465970
 [3,] 0.2465970
 [4,] 0.2465970
 [5,] 0.2465970
 [6,] 0.1108032
 [7,] 0.1108032
 [8,] 0.2465970
 [9,] 0.1108032
[10,] 0.2465970
[11,] 0.1108032
[12,] 0.1108032
[13,] 0.1108032
[14,] 0.2465970
[15,] 0.1108032
[16,] 0.2465970
[17,] 0.2465970
[18,] 0.2465970
[19,] 0.1108032
[20,] 0.1108032
[21,] 0.2465970
[22,] 0.2465970
[23,] 0.1108032
[24,] 0.1108032
[25,] 0.1108032
[26,] 0.2465970
[27,] 0.2465970
[28,] 0.1108032
[29,] 0.2465970
[30,] 0.2465970
[31,] 0.1108032
[32,] 0.2465970
[33,] 0.1108032
[34,] 0.2465970
[35,] 0.1108032
[36,] 0.1108032
[37,] 0.1108032
[38,] 0.1108032
[39,] 0.2465970
[40,] 0.2465970
[41,] 0.1108032
[42,] 0.2465970
[43,] 0.2465970
[44,] 0.1108032
[45,] 0.1108032
[46,] 0.2465970
[47,] 0.2465970
[48,] 0.2465970
[49,] 0.2465970
[50,] 0.2465970

> mean(mydata$Y)
[1] 0.06

そのようなデータでどうなるかを見てみると：

# Poisson regression model
myposfit <- glm(Y~v1+v2+v3+v4, data = mydata, family = poisson())
summary(myposfit)
vcov_correct <- vcov(myposfit)

# variance of poisson equals to its mean
mu <- myposfit$weights # weights in poisson is exp(Xb) = var = mu


W <- diag(mu)
# J_n^-1 = X'WX 
J_n_inv <- solve(crossprod(X*mu,X) # crossprod(X*mu,X) = t(X)%*%W%*%X

colnames(J_n_inv) <- colnames(vcov_correct)
rownames(J_n_inv) <- rownames(vcov_correct)

all.equal()のデフォルト設定(差分の閾値が1.5e-8)だと一致している。

> all.equal(vcov_correct, J_n_inv)
[1] TRUE

が、試しにいろいろな逆行列の計算も含めて、差分を見てみることにしよう。SV分解、QR分解、コレスキー分解で逆行列を計算してみる。

## svd
J_svd <- svd(t(X)%*%W%*%X)
J_svd_inv <- J_svd$v%*%diag(1/J_svd$d)%*%t(J_svd$u)

## qr decomposition
myposfit_qr <- qr.X(myposfit$qr)
#colnames(myposfit_qr) <- colnames(X)
#all.equal(myposfit_qr, sqrt(W)%*%X, check.attributes = F) # true

myqr <- qr(sqrt(W)%*%X, tol = 1e-10)
R <- qr.R(myqr)
Q <- qr.Q(myqr)

inv_R <- backsolve(R, diag(ncol(R)))
J_qr_inv <- inv_R %*% t(inv_R)

## cholesky decomposition
R_c <- chol(t(X)%*%W%*%X)
J_chol_inv <- solve(R_c)%*%solve(t(R_c))

元々の方法 (下記コードではnormalという名称)も含めてvcov()との差分を見てみると、QR分解以外が一番近しかった。これはglmが実際にQR分解を行っているからと思われる。

> ## see differences - QR is very close!!
> list(normal = J_n_inv, 
+      svd = J_svd_inv, 
+      qr = J_qr_inv, 
+      chol = J_chol_inv) %>%
+   map(function(x) x-vcov_correct)
$normal
              (Intercept)            v1            v2            v3
(Intercept) -7.744908e-02 -1.888284e-09  7.744907e-02 -1.181087e-08
v1           2.545499e-10  3.330669e-16 -2.545495e-10 -8.847090e-16
v2           7.744908e-02  1.888284e-09 -7.744907e-02  1.181087e-08
v3          -1.693194e-08 -7.147061e-16  1.693194e-08  1.332268e-15
v4          -1.375296e-08 -3.608225e-16  1.375296e-08  1.332268e-15
                       v4
(Intercept) -9.155280e-09
v1          -6.106227e-16
v2           9.155279e-09
v3           1.110223e-15
v4           8.881784e-16

$svd
              (Intercept)            v1            v2            v3
(Intercept) -4.878620e-01 -7.723943e-09  4.878620e-01 -3.172430e-09
v1          -1.598029e-08  5.551115e-16  1.598029e-08 -1.925543e-16
v2           4.878619e-01  7.723944e-09 -4.878619e-01  3.172430e-09
v3           3.226237e-09 -1.039099e-15 -3.226238e-09  4.440892e-16
v4          -7.536268e-09 -1.110223e-15  7.536268e-09  0.000000e+00
                       v4
(Intercept) -1.207726e-08
v1          -1.332268e-15
v2           1.207726e-08
v3          -2.220446e-16
v4          -4.440892e-16

$qr
              (Intercept)           v1 v2            v3            v4
(Intercept)  0.000000e+00 0.000000e+00  0 -2.220446e-16 -3.330669e-16
v1           0.000000e+00 0.000000e+00  0  6.938894e-18  0.000000e+00
v2           0.000000e+00 0.000000e+00  0  0.000000e+00  0.000000e+00
v3          -2.220446e-16 6.938894e-18  0  0.000000e+00  0.000000e+00
v4          -3.330669e-16 0.000000e+00  0  0.000000e+00  0.000000e+00

$chol
              const            v1            v2            v3            v4
const  1.337590e-01 -2.415076e-09 -1.337590e-01 -2.376583e-08 -2.441613e-08
v1    -1.800929e-09  5.551115e-16  1.800929e-09 -1.033895e-15 -8.326673e-16
v2    -1.337590e-01  2.415077e-09  1.337589e-01  2.376582e-08  2.441613e-08
v3    -2.376583e-08 -9.436896e-16  2.376582e-08  2.220446e-15  1.776357e-15
v4    -2.441613e-08 -7.216450e-16  2.441613e-08  1.776357e-15  8.881784e-16

とはいえ、差分%小数点4桁までで見てみると差はないので、大きな問題ではなさそう。

> list(normal = J_n_inv, 
+      svd = J_svd_inv, 
+      qr = J_qr_inv, 
+      chol = J_chol_inv) %>%
+   map(function(x) round((x-vcov_correct)/vcov_correct*100,3))
$normal
            (Intercept) v1 v2 v3 v4
(Intercept)           0  0  0  0  0
v1                    0  0  0  0  0
v2                    0  0  0  0  0
v3                    0  0  0  0  0
v4                    0  0  0  0  0

$svd
            (Intercept) v1 v2 v3 v4
(Intercept)           0  0  0  0  0
v1                    0  0  0  0  0
v2                    0  0  0  0  0
v3                    0  0  0  0  0
v4                    0  0  0  0  0

$qr
            (Intercept) v1 v2 v3 v4
(Intercept)           0  0  0  0  0
v1                    0  0  0  0  0
v2                    0  0  0  0  0
v3                    0  0  0  0  0
v4                    0  0  0  0  0

$chol
      const v1 v2 v3 v4
const     0  0  0  0  0
v1        0  0  0  0  0
v2        0  0  0  0  0
v3        0  0  0  0  0
v4        0  0  0  0  0

sandwich()の実装 - nが小さい場合

QR分解の関数を作る。

## create a function for later use 
vcov_qr <- function(mymat) {
  
  myqr <- qr(mymat, tol = 1e-10)
  R <- qr.R(myqr)
  Q <- qr.Q(myqr)
  
  inv_R <- backsolve(R, diag(ncol(R)))
  return(inv_R %*% t(inv_R))
}

sandwich推定値を計算してみると、結構な違いが出てしまっている。

# the meat part of sandwich 
resid <- mydata$Y - mu # Y - exp(X \beta) at the last iteration of IWLS in glm.fit

# V_n
observed_score <- X*resid
V_n <- crossprod(observed_score,observed_score)

# Sandwich estimator 
J_n_inv <- vcov_qr(sqrt(W)%*%X) # use my defined function using QRD
mysandwich <- J_n_inv%*%V_n%*%J_n_inv

# the above sandwich estimate is almost identical with {sandwich}
sandwich_correct <- sandwich(myposfit)

colnames(mysandwich) <- colnames(sandwich_correct)
rownames(mysandwich) <- rownames(sandwich_correct)

> all.equal(mysandwich, sandwich_correct)
[1] "Mean relative difference: 0.026304"

中身を見てみる。

> mysandwich
            (Intercept)         v1         v2         v3         v4
(Intercept)   1.0129734  0.7616622 -0.2228760 -0.6007226  0.5686506
v1            0.7616622  0.6329959 -0.1711880 -0.3634898  0.5458729
v2           -0.3795918 -0.1711880  0.5234677  0.2605280 -0.3600211
v3           -0.6007226 -0.3634898  0.2605280  0.5555911 -0.2840999
v4            0.5686506  0.5458729 -0.3600210 -0.2840999  0.7557730
> sandwich_correct
            (Intercept)         v1         v2         v3         v4
(Intercept)   1.0129733  0.7616622 -0.2228760 -0.6007226  0.5686506
v1            0.7616622  0.6329959 -0.1711880 -0.3634898  0.5458729
v2           -0.2228760 -0.1711880  0.3667520  0.2605280 -0.3600211
v3           -0.6007226 -0.3634898  0.2605280  0.5555911 -0.2840999
v4            0.5686506  0.5458729 -0.3600211 -0.2840999  0.7557730

一部かなりズレた結果になっていた。

> round((mysandwich - sandwich_correct)/sandwich_correct * 100,2)
            (Intercept) v1    v2 v3 v4
(Intercept)        0.00  0  0.00  0  0
v1                 0.00  0  0.00  0  0
v2                70.32  0 42.73  0  0
v3                 0.00  0  0.00  0  0
v4                 0.00  0  0.00  0  0

どのくらいずれるのかの確認として、行列の対角成分の差分%が5%を超える数を出してみると、多くはなかった。また、そもそも"glm.fit: fitted rates numerically 0 occurred"というワーニングが数回出ているような状況でもあったので、一致をさせるのがそもそもナンセンスなケースも含まれてはいそう。

iter <- 500
n <- 50
out <- data.frame(matrix(rep(NA,iter), ncol = ncol(X)))

for(i in 1:iter) {
  d <- datagen(n)
  mydata <- d$mydata
  X <- d$X

  myposfit <- glm(Y~v1+v2+v3+v4, data = mydata, family = poisson())
  mu <- myposfit$weights 
  W <- diag(mu)
  J_n_inv <- vcov_qr(sqrt(W)%*%X)
  #vcov_est <- vcov(myposfit)

  resid <- mydata$Y - mu 
  
  observed_score <- X*resid
  V_n <- crossprod(observed_score,observed_score)
  mysandwich <- J_n_inv%*%V_n%*%J_n_inv
  #mysandwich <- vcov_est%*%V_n%*%vcov_est

  sandwich_correct <- sandwich(myposfit)
  
  out[i,] <- (diag(mysandwich) - diag(sandwich_correct))/diag(sandwich_correct)*100

}

> apply(out, 2, function(x){sum(x>5)})
X1 X2 X3 X4 X5 
 4  0  5  0  0

ただ気がかりなのは、sandiwh::bread()とsandwich::meat()の中身と自分の推定値が全然合っていないことである。

> sandwich::bread(myposfit)
              (Intercept)         v1            v2          v3         v4
(Intercept)  1.051702e+09 32.9764418 -1.051702e+09 -60.6497562 -44.440282
v1           3.297644e+01 29.9347626 -4.436275e+00   0.5963174   7.054558
v2          -1.051702e+09 -4.4362746  1.051702e+09 -10.3497393 -36.744725
v3          -6.064976e+01  0.5963174 -1.034974e+01 104.2656629  55.365686
v4          -4.444028e+01  7.0545579 -3.674472e+01  55.3656865 105.916293
> J_n_inv
              [,1]        [,2]          [,3]        [,4]       [,5]
[1,]  2.103403e+07  0.65952884 -2.103403e+07 -1.21299512 -0.8888056
[2,]  6.595288e-01  0.59869525 -8.872549e-02  0.01192635  0.1410912
[3,] -2.103403e+07 -0.08872549  2.103403e+07 -0.20699479 -0.7348945
[4,] -1.212995e+00  0.01192635 -2.069948e-01  2.08531326  1.1073137
[5,] -8.888056e-01  0.14109116 -7.348945e-01  1.10731373  2.1183259

> sandwich::meat(myposfit, adjust = T)
            (Intercept)           v1          v2           v3           v4
(Intercept)  0.03700569 -0.026809759  0.03700569  0.015477219  0.022310361
v1          -0.02680976  0.047013988 -0.02680976 -0.025161767 -0.002446627
v2           0.03700569 -0.026809759  0.03700569  0.015477219  0.022310361
v3           0.01547722 -0.025161767  0.01547722  0.015477219  0.001145946
v4           0.02231036 -0.002446627  0.02231036  0.001145946  0.022310361
> sandwich::meat(myposfit, adjust = F)
            (Intercept)           v1          v2           v3           v4
(Intercept)  0.03330512 -0.024128783  0.03330512  0.013929497  0.020079325
v1          -0.02412878  0.042312589 -0.02412878 -0.022645590 -0.002201965
v2           0.03330512 -0.024128783  0.03330512  0.013929497  0.020079325
v3           0.01392950 -0.022645590  0.01392950  0.013929497  0.001031351
v4           0.02007932 -0.002201965  0.02007932  0.001031351  0.020079325
> V_n
           const         v1         v2          v3          v4
const  1.6652561 -1.2064392  1.6652561  0.69647486  1.00396623
v1    -1.2064392  2.1156294 -1.2064392 -1.13227952 -0.11009823
v2     1.6652561 -1.2064392  1.6652561  0.69647486  1.00396623
v3     0.6964749 -1.1322795  0.6964749  0.69647486  0.05156757
v4     1.0039662 -0.1100982  1.0039662  0.05156757  1.00396623

{rqlm}によるModified poissonの実行

野間先生が最近公開された{rqlm}も使ってみた。sandwich()を使う場合は下記のように書くと結果は一致した。{rqlm}の方が桁が小さくなると0に落とされているのでパッと見は結果が違うように見えるが実質な同じである。

library(rqlm)

# rqlm
rqlm(Y~., data=mydata, family=poisson, eform=F)
rqlm(Y~., data=mydata, family=poisson, eform=T)

# use sandwich()
data.frame(rr = coef(myposfit),
           rr_exp = exp(coef(myposfit)),
           se_rr = sqrt(diag(sandwich::sandwich(myposfit)))
) %>%
  mutate(
    lwr = rr - qnorm(0.975)*se_rr,
    upr = rr + qnorm(0.975)*se_rr,
    lwr_exp = exp(rr - qnorm(0.975)*se_rr),
    upr_exp = exp(rr + qnorm(0.975)*se_rr)
  )

Poisson, Modified poisson, Quasi-poissonの比較

参考までに、n=50のデータで普通のpoissonとmodified poissonの結果の比較をしてみる。

> summary(myposfit)$coefficients
               Estimate  Std. Error      z value   Pr(>|z|)
(Intercept) -21.9963357 4586.287488 -0.004796109 0.99617327
v1           -1.4633147    0.773754 -1.891188526 0.05859918
v2           18.5646941 4586.287385  0.004047870 0.99677028
v3           -0.5604871    1.444061 -0.388132442 0.69791802
v4            1.1496874    1.455447  0.789920475 0.42957421

> rqlm(Y~., data=mydata, family=poisson, eform=F)
                coef     SE       CL       CU P-value
(Intercept) -21.9963 1.0065 -23.9690 -20.0237  0.0000
v1           -1.4633 0.7956  -3.0227   0.0961  0.0659
v2           18.5647 0.6056  17.3777  19.7516  0.0000
v3           -0.5605 0.7454  -2.0214   0.9004  0.4521
v4            1.1497 0.8694  -0.5542   2.8536  0.1860

同じデータでquassi-poisson回帰を実行した結果も貼ってみる。dispersion parameterが約0.3とかなりunderdispersionであったらしいが、その補正を行っても標準誤差は非常に大きく、modified poissonの凄さが分かる（怖いくらい値が違う）。

myposfit2 <- glm(Y~v1+v2+v3+v4, data = mydata, family = quasipoisson())

> (dispersion <- sum((mydata$Y - myposfit2$fitted.values)^2/myposfit2$fitted.values)/(n-ncol(X)))
[1] 0.2951821

> summary(myposfit2)$coefficients
               Estimate   Std. Error      t value    Pr(>|t|)
(Intercept) -21.9963357 2491.7604701 -0.008827628 0.992995683
v1           -1.4633147    0.4203857 -3.480886055 0.001122803
v2           18.5646941 2491.7604144  0.007450433 0.994088402
v3           -0.5604871    0.7845681 -0.714389278 0.478676084
v4            1.1496874    0.7907540  1.453912780 0.152909869

Modified poissonにおけるP値の計算

ロバスト分散が分散の正規分布を用いて片側P値×2で計算したところ{rqlm}のP値と一致した。

> data.frame(
+   var = names(myposfit$coefficients),
+   p_values = unname(myposfit$coefficients/sqrt(diag(sandwich_correct))) %>%
+     map(function(x) {ifelse(x > 0, 
+                             pnorm(x, lower.tail = F)*2, 
+                             pnorm(x, lower.tail = T)*2)}
+     ) %>%
+     unlist()
+ ) %>%
+   mutate(p_values = round(p_values, 4))
          var p_values
1 (Intercept)   0.0000
2          v1   0.0659
3          v2   0.0000
4          v3   0.4521
5          v4   0.1860

small nでの修正ポアソン回帰の統計的性質評価シミュレーション

引き続き同じn=50のデータで、バイアスあるいは真値の復元、95%カバレッジ、alphaエラーと検出力の確認を行なった。シミュレーションコードは下記。

out_meandiff <- list(NULL)
out_coverage <- list(NULL)
out_pvalues <- list(NULL)
iter <- 5000

for (i in 1:iter) {
  
  ## generate data
  n <- 50
  d <- datagen(n)
  mydata <- d$mydata
  X <- d$X
  
  ## fit poisson
  myposfit <- glm(Y~v1+v2+v3+v4, data = mydata, family = poisson())
  sandwich_correct <- sandwich::sandwich(myposfit)
  
  ## get CIs
  res <- 
  data.frame(rr = coef(myposfit),
             se_rr = sqrt(diag(sandwich_correct))
  ) %>%
    mutate(
      lwr = rr - qnorm(0.975)*se_rr,
      upr = rr + qnorm(0.975)*se_rr)
  
  # bias and coverage
  mean_diff <- res$rr - true_beta
  coverage <- res$lwr <= true_beta & res$upr >= true_beta
  
  names(mean_diff) <- names(myposfit$coefficients)
  names(coverage) <- names(myposfit$coefficients)

  ## get pvalues
  pvalues <-  unname(myposfit$coefficients/sqrt(diag(sandwich_correct))) %>%
    map(function(x) {2*(1-pnorm(abs(x)))}
    ) %>%
    unlist()
  
  names(pvalues) <- names(myposfit$coefficients)
  
  # store
  out_meandiff[[i]] <- c(iter = i, mean_diff)
  out_coverage[[i]] <- c(iter = i, coverage)
  out_pvalues[[i]] <- c(iter = i, pvalues)
}

シミレーションの結果としては、下のRコードスニペットにまとめた。
変数v2の対数リスク比の真値は0.8なので概ねバイアスなく推定できており、真値が0であるはずの変数に対しては、バイアスが~10%くらいであった。

真値が95%信頼区間に含まれた割合は90-95%で、被覆確率が少し低い変数もあった。

v1は信頼区間が狭いので棄却率が高まり10%に近かった。v3とv4は被覆確率が95%だったので概ね5%の棄却率であったが、v2については検出力が20%程度に留まっていた。nが小さいため検出力が低いのは当たり前ではある。しかしながら、普通のポアソン回帰ではSEが爆発していたため、普通のポアソンにおける検出力はほぼ0でカバレッジもほぼ100％になってしまっていたと思われる。

nが大きい場合でも普通のポアソン回帰と修正ポアソン回帰で分散に違いがあったことを考えると、nに関係なく修正回帰モデルにすることよって、今回のデータにおいてはより適切な評価につながっていると言える。

There were 50 or more warnings (use warnings() to see the first 50)
> 
> # bias or true effect detection
> out_meandiff %>% 
+   bind_rows() %>%
+   select(-iter, -`(Intercept)`) %>%
+   colMeans()
          v1           v2           v3           v4 
-0.099992599  1.014696863 -0.111147744  0.006184029 
> 
> # coverage
> out_coverage %>% 
+   bind_rows() %>%
+   select(-iter, -`(Intercept)`) %>%
+   colSums()/iter
    v1     v2     v3     v4 
0.8996 0.9060 0.9408 0.9448 
> 
> 
> # significant results at 5% level
> out_wide_pval <-
+   out_pvalues %>% 
+   bind_rows() 
> 
> # alpha error or 1 minus beta error
> colSums(out_wide_pval[,3:6] <= 0.05)/iter
    v1     v2     v3     v4 
0.1004 0.1982 0.0592 0.0552

感想

  • sandwich::sandwich()の完全な再現には至らなかったが、概ね理解を深められるレベルには到達したと思う。
  • 今回示した例だと、ロバスト分散の方が標準誤差が小さいという結果で、相関するデータで生じる「標準誤差が過小評価されているのでロバスト分散を使うと標準誤差が大きくなる」という動きとは逆だったことは気になりポイントでもう少し考察が必要。
  • sandwich::bread()とsandwich::meat()の中身は謎すぎるので今後の課題で、Object-Oriented Computation of Sandwich Estimatorsの読み込みが必要と考えている。
  • そういえば、ロバスト分散を使った場合のP値ってどうなるんだろう？{rqlm}だとP値も出ているのが良い。 上に内容として追記。

Discussion