ねこすたっと

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

2値アウトカムに対する回帰モデルいろいろ(2)：デルタ法・ブートストラップ法[R]

dichotomous bootstrap maths/theory R

今回はロジスティック回帰モデル（ロジット-二項モデル）からオッズ比（odds ratio, OR）ではなくリスク比（risk ratio, RR）を推定する方法をまとめる。

ロジスティック回帰モデルとは、アウトカムの発生確率pを

P(Y|X,Z) = \frac{1}{1+e^{-(\beta_0 + \beta_1 X + \beta_2 Z)}}

のように患者背景因子で回帰するモデル。

これを使えば、

\begin{align} p_1 &= P(Y|X=1,Z) = \frac{1}{1+e^{-(\beta_0 + \beta_1 + \beta_2 Z)}} \\ p_0 &=P(Y|X=0,Z) = \frac{1}{1+e^{-(\beta_0 + \beta_2 Z)}} \end{align}

なので、

\begin{align} RR = \frac{p_1}{p_0} = \frac{1+e^{-(\beta_0 + \beta_2 Z)}}{1+e^{-(\beta_0 + \beta_1 + \beta_2 Z)}} \end{align}

としてしまえばいいんじゃないの？と考えてしまいたくなるが、推定誤差まで正しく推定できていないと適切な方法とは言えない。

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

前回と同様に自然流産歴と不妊の関連をみるinfertデータを用意する。

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

デルタ法

デルタ法とは「g(X)をXの平均の周りでTaylor展開することで、Xの関数g(X)の分散を近似的に求める方法」のこと。

cran.r-project.org

ロジスティック回帰モデルで係数βの推定誤差が求められるので、RRをβの関数と考えれば、上記の方法を使ってRRの推定誤差が求められるという話。

見やすさのため

\begin{align} e_1 &= 1+e^{-(\beta_0 + \beta_1 + \beta_2 Z)} \\ e_0 &= 1+e^{-(\beta_0 + \beta_2 Z)} \end{align}

と書くことにすると、さらにRR=g(β)と書くことにすると、

RR = g(\beta) = \frac{1+e_0}{1+e_1}

となる。

デルタ法より、

var[g(\beta)] \simeq \left\{ \frac{\partial g(\beta)}{\partial \beta} \right\}^T var(\beta) \left\{ \frac{\partial g(\beta)}{\partial \beta} \right\}

となって、var(β)はモデルを当てはめた結果から抽出できるので、g(β)を \beta_0, \beta_1, \beta_2 で偏微分すればよい。

高校数学を思い出して頑張ると、

\begin{align} \frac{\partial g(\beta)}{\partial \beta} &= \left(\frac{\partial g(\beta)}{\partial \beta_0}, \frac{\partial g(\beta)}{\partial \beta_1}, \frac{\partial g(\beta)}{\partial \beta_2} \right) \\ &= \left(\frac{e_1 - e_0}{(1+e_1)^2}, \frac{e_1(1+e_0)}{(1+e_1)^2}, \frac{z(e_1 - e_0)}{(1+e_1)^2} \right) \end{align}

となる。

これをRで実装してみる（deltamethodというパッケージはあるようだが、ここは手計算で）。まずはロジスティック回帰モデルの当てはめる。係数βの点推定値をbeta、βの共分散行列をvcov_betaとして抽出する。またRRを計算するときに、共変量ageは特定の値に決めておかないといけないので、集団の平均年齢（このデータでは31.5歳）とする。

fit <- glm(case ~ spont + age, data=d,
           family = binomial(link="logit"))
beta <- coef(summary(fit))
vcov_beta <- summary(fit)$cov.unscaled
z <- mean(d$age)

先程の数式をひたすら書く。

b1 <- beta[1]
b2 <- beta[2]
b3 <- beta[3]
e1 <- exp(-b1-b2-b3*z)
e0 <- exp(-b1-b3*z)

d_b1 <- (e1-e0)/(1+e1)^2
d_b2 <- e1*(1+e0)/(1+e1)^2
d_b3 <- z*(e1-e0)/(1+e1)^2
d_b <- as.vector(c(d_b1,d_b2,d_b3))

RR_est <- (1+e0)/(1+e1)
RR_var <- t(d_b) %*% vcov_beta %*% d_b
RR_CI <- c(RR_est - 1.96*sqrt(RR_var),RR_est + 1.96*sqrt(RR_var))

> RR_est
[1] 2.66

> RR_CI
[1] 1.63 3.69

ちなみに、修正ポアソン回帰モデルで求めたRRは2.65 (95CI: 1.80 to 3.89)だったので、それっぽい数値にはなっていると思う（しらんけど）。

cran.r-project.org

ブートストラップ法

ブートストラップ法は、元データから重複を許して同じ数だけ抽出したブートストラップサンプルにおいて目的の統計量（今はRR）を計算し、その分布をもとにして推定誤差を求める方法。何回も抽出・計算を繰り返すので、計算負荷が大きい。

logistic( )を使えば比較的簡単に実行できる。デルタ法の説明で述べたように、共変量についてはどんな値にするかを決めないといけない。デフォルトでは、カテゴリー変数の場合は0に相当するレベルに、連続変数の場合には最低値に設定されてしまうみたい。下のコードようにfixcovをlist形式で指定する（ここでも平均年齢を指定）。

カテゴリー変数の中身は数値にしておいた方がいいかも（別のデータで文字列のままだと走らなかった）。

result <- logisticRR(case ~ spont + age, data=d, 
                     fixcov = list("age"=mean(d$age)),
                     boot=TRUE, n.boot=1000)

n.bootでサンプリング回数を指定しているので、ここではRRが1000回得られていて、result$boot.rrで取り出すことができる。例えばヒストグラムを描いてみると、

hist(result$boot.rr)

サンプリングされた1000個のRRのヒストグラム

RRの平均と分散は

> mean(result$boot.rr)
[1] 2.77

> var(result$boot.rr)
[1] 0.344

なので、この分散をもとにして95%信頼区間を計算することができる（平方根をとればSEになる）が、左右非対称な分布を見てしまうと分位点をもとにする方がいいように思う。

> quantile(result$boot.rr,prob=c(0.025,0.25,0.5,0.75,0.975))
 2.5%   25%   50%   75% 97.5% 
 1.84  2.36  2.69  3.10  4.08

上記の結果から、95%信頼区間は(1.84, 4.08)となる。

おわりに

  • ブートストラップについてはこちらの記事も参考にしてください。

necostat.hatenablog.jp

ランキング参加中
統計

参考資料

  • ロジスティック回帰モデルから相対リスクを計算する方法について

cran.r-project.org

  • 日本語で書かれいてわかりやすい記事です

qiita.com