🙆♀️
Modified Poissonを理解する - PART2 実装編
ゴール
- スクラッチ実装を通して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
スクラッチ
ポアソン回帰の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と同じ中身であるので（スコア関数の直積を
Sandwich estimatorで新たに登場するmeat部分（理論編の表記でいうとB）を計算する。残差の二乗を重みとした
# 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
> 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
> 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