ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

2値アウトカムに対する回帰モデルいろいろ(1)：glm関数とロバスト分散推定法[R]

maths/theory dichotomous tidyverse reporting R neco

アウトカムが2値変数の場合、一般的にロジスティック回帰モデルが用いられることが多い。 この場合の効果指標はオッズ比（odds ration, OR）だが、アウトカムの発生頻度が高い場合（>10%ルール）、オッズ比はリスク比（risk ratio, RR）の近似にならない。 この記事ではリスク比やリスク差（risk difference, RD）を推定する方法をまとめてみる。

標準誤差と信頼区間

点推定値だけでなく信頼区間も示すように言われると思う。95%信頼区間は標準誤差（standard error, SE）を使って、

Estimate \pm 1.96 \times SE

で計算されることが多い。 これは「最尤法で推定するときサンプル数が十分あれば正規分布に近似して信頼区間を求めれば良い」ことに基づいているが、サンプル数が少ないときはt分布を使わないといけない。

自由度（degree of freedom, df）が30以上のt分布は正規分布とほとんど同じなので、df≧30（サンプル数-パラメータ数≧30）なら正規分布近似で問題なさそう。*1

「95%信頼区間がNULL値（指数変換前なら0）ことがP<0.05と同値」となるのは、

  • 95%信頼区間が正規分布近似を用いて計算されている
  • P値がWald検定を用いて計算されている

とき。Wald検定は

Wald = \left( \frac{\beta - NULL}{SE(\beta)}\right)^2

が自由度1のカイ2乗分布 [tex: \chi ^2{df=1}] に従うことを使った検定方法だが、標準正規分布に従う変数を2乗すると [tex: \chi ^2{df=1}] 分布に従うので、結局同じ計算をしていることになるため。

パッケージとデータの準備

ここではinfertデータを使って、自然流産歴と不妊の関連をみてみる。元の研究がケースコントロール研究なので発生頻度（risk）は求められないが、脳内で適当なデザインに読み替えて...。

library(tidyverse)
d <- infert %>% 
  mutate(spont = if_else(spontaneous==0, 0, 1)) %>%
  rowid_to_column() %>%
  rename(ID=rowid) %>%
  select(ID,age,case,spont)

> head(d)
  ID age case spont
1  1  26    1     1
2  2  42    1     0
3  3  39    1     0
4  4  34    1     0
5  5  35    1     1
6  6  36    1     1
  • ID：患者ID。行番号から作成。
  • age：年齢（歳）
  • case：不妊症（0/1）
  • spont：自然流産歴。元のデータは「0回/1回/2回以上」だったが、「0回/1回以上」に変換した。

偶現表は次のとおり。

> d %>% with(table(case, spont))
    spont
case   0   1
   0 113  52
   1  28  55

不妊割合は「自然流産歴あり」で51.4%、「自然流産歴なし」で19.9％なので、crudeな結果は次のとおり。

  • cRD(あり-なし) = 0.315
  • cRR(あり/なし) = 2.59
  • cOR(あり/なし) = 4.27

発生割合が高いので、RRとORに乖離がみられている。

ORモデル

ロジスティック回帰モデルは効果指標としてORが推定されるので、ORモデルと呼ぶことがある。 ロジスティック回帰モデルを当てはめて、そこからRRを推定する方法もあるが、それは別の記事で書く予定。

ロジット-二項モデル（logit-binomial model）

一般化線形モデル（generalized linear model, GLM）ではリンク関数と分布ファミリーでモデルの型が決まるので、ロジスティック回帰モデルはロジット-二項モデル（logit-binomial model）と呼ばれる。*2

まずはglm( )でモデルを当てはめる。

fit <- glm(case ~ spont + age, data = d,
           family = binomial(link="logit"))

結果はsummary( )で取り出せばよいが、adjusted ORを求めるためには指数変換が必要。下のコードでは結果の中で必要な係数の部分をcoef( )で抽出・加工してORと95%信頼区間にしている（Std. Errorは半角スペースを含んでいるので引用符で囲う必要あり）。select( )はMASSパッケージ内の同名の関数とバッティングすることがあるので、エラーが出たら下のようにdplyr::をつける。

coef(summary(fit)) %>% as.data.frame() %>%
  mutate(OR = exp(Estimate),
         OR_LL = exp(Estimate - 1.96*`Std. Error`),
         OR_UL = exp(Estimate + 1.96*`Std. Error`)) %>%
  dplyr::select(OR, OR_LL, OR_UL)

                   OR      OR_LL     OR_UL
(Intercept) 0.1140111 0.01846273 0.7040411
spont       4.4492701 2.51130981 7.8827410
age         1.0242754 0.97001840 1.0815673

結果を表示する方法その2は、broomパッケージのtidy( )を使う方法。broomパッケージはtidyverseパッケージをインストールするときに一緒にインストールされるが、読み込みのときは別に明示する必要あり。

> library(broom)
> tidy(fit, exponentiate = TRUE, conf.int = TRUE)
# A tibble: 3 x 7
  term        estimate std.error statistic     p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>       <dbl>    <dbl>     <dbl>
1 (Intercept)    0.114    0.929     -2.34  0.0194        0.0178     0.686
2 spont          4.45     0.292      5.12  0.000000313   2.54       7.98 
3 age            1.02     0.0278     0.864 0.388         0.970      1.08

微妙な数値の違いは丸めの問題だけなのか...ちょっと調べきれてない。

方法その3は、bistat3パッケージのeform( )を使う方法。関数名からするとStataと同じ結果になるのだろうか。

> library(biostat3)
> eform(fit, level=0.95, method="Delta")
            exp(beta)      2.5 %    97.5 %
(Intercept) 0.1140111 0.01846335 0.7040175
spont       4.4492701 2.51133621 7.8826582
age         1.0242754 0.97001937 1.0815662

こちらは方法その1（手計算）と同じになった。method="Profile"とすれば尤度に基づいた方法で信頼区間を計算できるみたい。

RRモデル

効果指標としてRRが推定できるモデル。対数-二項モデル（log-binomial model）と対数-ポアソンモデル（log-Poisson model）が用いられる。後者はロバスト推定法（robust estimation）*3を用いないといけないので、「修正ポアソンモデル（modified Poisson model）」と呼ばれることが多い。*4

ちなみに後者の方がいいらしい。*5

対数-二項モデル（log-binomial model）

まずはglm( )でモデルを当てはめる。

fit <- glm(case ~ spont + age, data=d,
           family=binomial(link="log"))

なぜかtidy( )でconf.int=TRUEとするとダメだったのでeform( )を使う。

> eform(fit, level=0.95, method="Delta")
             exp(beta)      2.5 %    97.5 %
(Intercept) 0.09968511 0.03419986 0.2905602
spont       2.70800603 1.85278841 3.9579785
age         1.02130048 0.99023777 1.0533376

修正ポアソンモデル（modified log-Poisson model）

方法1：glm( )+coeftest( )を使う

まずglm()でモデルを当てはめる。

fit <- glm(case ~ spont + age, data=d,
           family=poisson(link="log"))

このままだとロバスト推定法が適用されていない。Rだとこれが面倒くさい。あとで指数変換加工をするためresultという名前で結果を保持しておく。

library(lmtest)
library(sandwich)
result <- coeftest(fit, vcov = sandwich)

オプションとしては、

  • df = fit$df.residualを追加すると検定にt分布を用いることができる
  • vcov = vcovHC(fit,type="HC1")とすればStataと同じ結果になる（らしい）

先程の推定結果を指数変換するときに、tidy( )やeform( )がうまく働いてくれないので、以下のコードで手計算する。coeftest( )の結果はそのままデータフレームにできないみたい。

estimate <- result[,1]
se <- result[,2]
result_summary <- cbind(estimate, se) %>% as.data.frame() %>%
  mutate(RR = exp(estimate),
         RR_LL = exp(estimate - 1.96*se),
         RR_UL = exp(estimate + 1.96*se)) %>%
  dplyr::select(RR, RR_LL, RR_UL)

> result_summary
                   RR      RR_LL     RR_UL
(Intercept) 0.1245098 0.04032432 0.3844501
spont       2.6467396 1.80153550 3.8884777
age         1.0145279 0.98240012 1.0477065

方法2：geeglm( )を使う

glm( )では直接的にロバスト分散推定を組み込めないが、一般化推定方程式（generalized estimating equation, GEE）で用いるgeeglm( )なら組み込める。

library(geepack)
fit <- geeglm(case ~ spont + age, data=d,
              corstr="independence", id=ID, std.err="san.se",
              family=poisson(link="log"))

GEEではクラスター内の相関構造をcorstrで指定できる。ここでは"independence"を指定する。*6 所属するクラスターを示す変数はidで指定する。ここでは各患者に割り振ったIDを指定した（これなら同一クラスターにみなされる人はいない）。

この方法だとeform( )で表示可能。ちなみにtidy( )でも同じように上手く表示される。

> eform(fit, level=0.95, method="Delta")
            exp(beta)      2.5 %   97.5 %
(Intercept) 0.1245098 0.04032674 0.384427
spont       2.6467397 1.80155262 3.888441
age         1.0145279 0.98240179 1.047705

方法1と微妙に数値が違うが、丸め桁数程度か。

方法3：svyglm( )を使う

surveyパッケージに含まれている各種モデリング用関数は、基本的にロバスト推定（Horvitz-Thompson型推定法）を行うようになっているのでこれを使う。

library(survey)
dsgn <- svydesign(~1, weights = ~1, data = d)
fit <- svyglm(case ~ spont + age, design = dsgn,
              family = poisson(link="log"))

結果はeform( )でもtidy( )でも大丈夫みたい。

> eform(fit, level=0.95, method="Delta")
            exp(beta)      2.5 %    97.5 %
(Intercept) 0.1245098 0.04023332 0.3853196
spont       2.6467396 1.80014731 3.8914763
age         1.0145279 0.98233677 1.0477740

> tidy(fit, exponentiate = TRUE, conf.int = TRUE)
# A tibble: 3 x 7
  term        estimate std.error statistic    p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>      <dbl>    <dbl>     <dbl>
1 (Intercept)    0.125    0.576     -3.61  0.000365     0.0402     0.385
2 spont          2.65     0.197      4.95  0.00000139   1.80       3.89 
3 age            1.01     0.0165     0.877 0.381        0.982      1.05

結果も概ね同じ。

RDモデル

効果指標としてRDが推定されるモデルとしては、恒等-二項モデル（identity-binomial model）と恒等-ガウスモデル（idensity-Gaussian model）がある。後者は修正ポアソンモデルと同様、ロバスト分散推定を用いる必要あり。

恒等-二項モデル（identity-binomial model）

glm()でモデルを当てはめる。

fit <- glm(case ~ spont + age, data=d,
           family=binomial(link="identity"))

今回は指数変換しなくてよい。

> tidy(fit, exponentiate = FALSE, conf.int = TRUE)
# A tibble: 3 x 7
  term        estimate std.error statistic      p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>        <dbl>    <dbl>     <dbl>
1 (Intercept)  0.106     0.172       0.616 0.538        -0.219      0.446 
2 spont        0.318     0.0592      5.38  0.0000000745  0.201      0.431 
3 age          0.00294   0.00529     0.555 0.579        -0.00735    0.0135

恒等-ガウスモデル（identity-Gaussian model）

名前は難しそうだが、何のことはない普通の線形回帰モデル。なのでlm( )でもOK。

方法1：glm()+coef()を使う

fit <- glm(case ~ spont + age, data=d,
           family = gaussian(link="identity"))
library(lmtest)
library(sandwich)
coeftest(fit, vcov = sandwich)

z test of coefficients:

            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.04841    0.17632    0.27     0.78    
spont        0.32249    0.05929    5.44  5.4e-08 ***
age          0.00467    0.00535    0.87     0.38    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

方法2：geeglm( )を使う

fit <- geeglm(case ~ spont + age, data=d,
              corstr = "independence", id=ID, std.err="san.se",
              family=gaussian(link="identity"))

> tidy(fit, exponentiate = FALSE, conf.int = TRUE)
# A tibble: 3 x 7
  term        estimate std.error statistic      p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>        <dbl>    <dbl>     <dbl>
1 (Intercept)  0.0484    0.176      0.0754 0.784        -0.297      0.394 
2 spont        0.322     0.0593    29.6    0.0000000536  0.206      0.439 
3 age          0.00467   0.00535    0.763  0.383        -0.00581    0.0152

方法3：svyglm( )を使う

fit <- geeglm(case ~ spont + age, data=d,
              corstr = "independence", id=ID, std.err="san.se",
              family=gaussian(link="identity"))

結果は割愛。

おわりに

  • family-linkに縛られ過ぎなくていいが、ロバスト分散法を使うこと
  • 長い休養を経て2シーズン目に突入。昨年からもう1匹加わりました（茶トラ・オス・10か月）

ランキング参加中
統計

参考資料

  • Ashley I. Naimi and Brian W. Whitcomb. Estimating Risk Ratios and Risk Differences Using Regression. Am J Epidemiol. 2020;189(6):508–510
  • Louise-Anne McNutt, Chuntao Wu, Xiaonan Xue, and Jean Paul Hafner. Estimating the Relative Risk in Cohort Studies and Clinical Trials of Common Outcomes. Am J Epidemiol 2003;157:940–943

*1:ちなみに統計ソフトの出力で"z-value"とか"z-distribuiton"とか出てくることがあるが、このzはz分布=標準正規分布（平均0, 標準誤差1）のこと。

*2:二項ロジットモデルの方が一般的なんだろうけど、英語の語順やこの後のモデルとの整合性を考えてこうした

*3:robust standard error estimation, robust variance estimationなど

*4:0/1しかとらない変数に、非負整数をとりうるポアソン分布を当てはめるので、分散の推定に工夫が要りますよ、という話

*5:Chen et al. Comparing performance between log- binomial and robust Poisson regression models for estimating risk ratios under model misspecification. BMC Medical Research Methodology (2018) 18:63

*6:全員のIDが違うので"exchangeable"としても同じ