Hazy Ideas

日々の勉強の気づきを書き出しています

【R】既存データを使って相対リスク計算の試算をしてみる

疫学 R

公開されている疫学文献データを用いて、相対リスクの計算を行ってみます。主な目的は２つです。

・データに対する分布を検討する

・調整モデル、複数汚染物質モデルなどを検討する

使用するパッケージはdlnmです。このパッケージを作成した研究者が解説情報を公開しているため、それらを参考にします。２つ目の方が新しい情報が含まれています。

http://www.ag-myresearch.com/uploads/1/3/8/6/13864925/gasparrini_jss2011.pdf

https://github.com/gasparrini/2011_gasparrini_JSS_Rcode/blob/master/Rcode.R

dlnmパッケージには、chicagoNMMAPSというデータセットが含まれます。上記pdfによると「1987年から2000年のシカゴ市において、大気汚染と気温という2つの環境要因が全死因死亡率に与える影響を取り上げている」だそうです。なお原著論文は以下です。

https://www.healtheffects.org/system/files/SametI-final.pdf

データの確認

本データセットには、以下のデータが5114セット含まれています。

#date:1987年から2000年までの日付
#time：順番（トレンド調整に用いる）
#year:年
　#month:月
　#day:日
　#dow:曜日

#death: 全死亡数（非事故）
#cvd: 心血管系死亡数
#resp： 呼吸器疾患死亡数

#temp： 平均気温
#dptp：露点温度
#rhum：平均相対湿度
#pm10：PM10
#o3: オゾン

ライブラリーをインポートします。

library(dlnm)
library(dplyr)
library(broom)
library(splines)

まずはデータを眺めてみます。

#変数の項目を確認
head(chicago)

データの外観を見てみましょう。日付に関する項目はdateだけに絞り、死亡と環境に関する変数を含めて、図にしてみます。

#全体像を眺める
chi <- select(chicago, date,death,cvd,resp,temp,pm10,o3)
plot(chi)

死亡（全死亡、循環器死亡、呼吸器死亡）と大気環境（気温、PM10、オゾン）に、見たところわかりやすい相関関係はないように見えます。気温とPM10、オゾンには相関がみられるので、調整が必要に思えます。

次に数値的な情報として、記述統計表を簡単に眺めてみます。

summary(chi)

大気汚染物質については政策等によって環境が改善していることがあります。ということで、年ごとの違いがないか、年でグループ分けした情報も作成します。

#年ごとの平均値を見てみる
chi_sum <- 
  chicago %>%
  group_by(year) %>%
  summarize(
    No. = n(), #総数
    Death = mean(death),　#平均
    CvdD = mean(cvd), #平均
    RespD = mean(resp), #平均
    Temp = mean(temp), #平均
    PM10 = mean(pm10,na.rm = TRUE),#平均
    #欠損値が含まれるため、欠損を除去した平均値を求めるna.rm=T
    Ozone = mean(o3)　#平均
  )
chi_sum
#出力結果
# A tibble: 14 x 8
    year   No. Death  CvdD RespD  Temp  PM10 Ozone
   <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1  1987   365  117.  57.2  8.44 11.2   34.2  20.4
 2  1988   366  119.  58.1  8.90  9.96  35.8  22.0
 3  1989   365  118.  52.9  9.65  9.09  37.1  21.9
 4  1990   365  116.  50.8  9.23 10.9   34.4  20.8
 5  1991   365  117.  51.7  9.34 10.8   33.4  21.1
 6  1992   366  114.  49.7  8.70  9.34  34.6  18.8
 7  1993   365  119.  51.5  9.70  9.29  32.2  19.1
 8  1994   365  118.  51.1  9.14 10.0   36.1  20.0
 9  1995   365  120.  53.1  9.78 10.1   33.2  19.9
10  1996   366  114.  49.2  9.29  8.55  31.9  18.6
11  1997   365  111.  46.9  8.95  9.28  31.3  19.4
12  1998   365  110.  46.6  9.41 12.1   33.8  19.7
13  1999   365  114.  48.3  9.51 10.8   32.8  21.0
14  2000   366  110.  45.4  8.48 10.0   32.6  19.1

検定はかけてませんが、あまり違いがないように見えます。なので期間や年による階層化はしなくてもよいという前提で話を進めます。

今回はPM10を対象に、呼吸器死亡への影響がどの程度あったか、検証してみます。一般化線形モデルglmに当てはめて検討しますが、その際にデータの分布を仮定する必要があります。一般に死亡データは、全人口に対してその日の死者数の割合は少ないということで、確率が低い事象の分布として捉えられることが多いように思います。

分布を検討する前に、まずは呼吸器死亡のデータの外観を眺めます。

ggplot(chicago) +
  aes(x = resp) +
  geom_histogram(bins = 50L)

左右均等ではなく、裾野が広がっているように思います。正規分布が当てはめられるか確認するために、正規性の検定として知られるShapiro-Wilk検定を行います。帰無仮説は「正規性がある」なので、p値が0.05以下であれば、正規性がないということになります。なお、shapiro.test関数の制約でデータ5000までしか読めないので、ランダムに5000個のデータを抽出しました（論文的にはダメな気がするので、他の方法を検討したほうがよいです）。

chicago$resp %>% 
  sample(size=5000) %>%
  shapiro.test()
#結果

    Shapiro-Wilk normality test

data:  .
W = 0.96948, p-value < 2.2e-16

p<0.05であり、正規分布は仮定できないことがわかりました。こうした死亡データの場合、ポアソン分布か二項分布が用いられます。また、ばらつきが大きすぎる場合（過分散）には、疑似ポアソン分布か負の二項分布が用いられます。あてはまりがよいものを用いるこのになると思います（どの分布がよいか機械的に検定する方法はないはず？この辺りはまだ勉強不足です）。

#ポアソン分布
#PM10と呼吸器死亡
lm_pm_p <- glm(resp ~ pm10, data = chicago, family=poisson())
tidy(lm_pm_p)
#estimate -0.001320384


#疑似ポアソン分布
#PM10と呼吸器死亡
lm_pm_qp <- glm(resp ~ pm10, data = chicago, family=quasipoisson())
tidy(lm_pm_qp)
#estimate -0.001320384

どちらの分布でも回帰係数estimateの結果は同じでした。（明確な根拠がないですが）原著論文ではポアソン分布が使われているようでしたので、poisson()を用いることにします。

モデル化とRR計算

では線形モデルを仮定してglmでモデル化してみます。単汚染物質モデル、二汚染物質モデルを作成してみます。外観の確認で、気温に相関がありそうでしたのでtempを調整に入れています。

#1汚染物質モデル
single_pol <- glm(resp ~ pm10 + temp, data = chicago, family=poisson())
#tidy:推定値に関する情報をデータフレームとして返す
tidy(single_pol)
#結果
# A tibble: 3 x 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)  2.29     0.00953     241.   0        
2 pm10         0.00103  0.000263      3.94 8.16e-  5
3 temp        -0.0117   0.000460    -25.4  1.12e-142
#2汚染物質モデル
two_pol <- glm(resp ~ pm10 + o3 + temp, data = chicago, family=quasipoisson())
tidy(two_pol)
#結果
# A tibble: 4 x 5
  term         estimate std.error statistic  p.value
  <chr>           <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)  2.28      0.0138     166.    0       
2 pm10         0.00100   0.000296     3.39  6.95e- 4
3 o3           0.000596  0.000641     0.931 3.52e- 1
4 temp        -0.0120    0.000596   -20.1   3.39e-86

オゾンを調整因子に入れることで、PM10の回帰係数が0.00103から0.00100に、少しだけ下がりましたが、p値も含め、それほど結果は変わらないということがわかりました。では、dlnmパッケージの開発研究者が例示していた調整因子を踏まえて動かしてみます。


#上記に加え露点温度、相対湿度、長期トレンドで調整
#ns: 自然立方体スプライン  #ns(time)：トレンドを調整する
final <- glm(resp ~ pm10 + o3 + temp + dptp + rhum + ns(time), data = chicago, family = poisson())
tidy(final)
#結果
# A tibble: 7 x 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)  2.12     0.0360       58.8  0        
2 pm10         0.00122  0.000308      3.96 0.0000761
3 o3           0.00154  0.000692      2.23 0.0258   
4 temp        -0.00658  0.00263      -2.51 0.0122   
5 dptp        -0.00270  0.00151      -1.79 0.0740   
6 rhum         0.00258  0.000622      4.15 0.0000335
7 ns(time)     0.0448   0.0301        1.49 0.137

いよいよ相対リスクRRです。RRは、回帰係数estimateの指数を取ることで計算できます。

#相対リスクの計算
#conf.int = TRUEを指定すると信頼区間も出力
#出力結果のうち、対象の列(1,2,6,7)のみを選択

#あとで図にしたいのでデータ型をdata.frameにしておく
RR <- tidy(final, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE)[,c(1,2,6,7)] %>% as.data.frame()
RR

#結果
         term  estimate  conf.low conf.high
1 (Intercept) 8.3210961 7.7532279 8.9299969
2        pm10 1.0012190 1.0006116 1.0018199
3          o3 1.0015442 1.0001848 1.0029030
4        temp 0.9934412 0.9883428 0.9985674
5        dptp 0.9973037 0.9943543 1.0002627
6        rhum 1.0025828 1.0013613 1.0038053
7    ns(time) 1.0457732 0.9858924 1.1093132

PM10やオゾンに正の関連があることがわかりました。見やすくするために、図示します。使用するのはフォレストプロットという、各推定値・信頼区間を縦に並べて表示するグラフです。

#RRのフォレストプロットの作成
Fplot <- RR %>%
  #intercept情報を削除
  filter(term != "(Intercept)") %>%
  #グラフの構造説明
  ggplot(aes(x=estimate, y=term))+
  #軸ラベルの作成
  labs(x = "RR", y = "Factors") +
  #estimateを点表示
  geom_point()+
  #信頼区間をエラーバーとして表示
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high))
Fplot

グラフが作成できました。ですが、トレンドns(time)の幅が広くて見づらいので、取り除きたいです（関連はないとはいえ、positive方向に偏っているので無視してよいとは言えない気がしますが）。またグラフ記載の順番が、モデルに入れた順番と違うので見づらいです。これを修正します。

Fplot2 <- RR %>%
  #intercept、ns(time)情報を削除
  filter(term != "(Intercept)") %>%
  #トレンド情報をグラフから除く
  filter(term != "ns(time)") %>%
  #y軸の順番を指定する
  mutate(term = fct_relevel(
    term,
    "rhum","dptp","temp","o3","pm10")) %>%
  #グラフの構造説明
  ggplot(aes(x=estimate, y=term))+
  #軸ラベルの作成
  labs(x = "RR", y = "Factors") +
  #estimateを点表示
  geom_point()+
  #信頼区間をエラーバーとして表示
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high))+
  #RR=1を示す線を表示
  geom_vline(xintercept = 1, linetype="dashed")
Fplot2

これで見やすくなりました。呼吸器死亡の相対リスクは、PM10、オゾン、相対湿度の増加により増加し、露点温度では関連がなく、気温上昇により低下するということが見て取れました。

まとめ

今回は相対リスクを実際のデータを用いて計算してみました。Rのコードに当てはめて計算するのは簡単ですが、データの特性を考えて分布を当てはめたり、調整因子の入れ方を検討したりする必要があることが分かります。疫学研究は様々な専門分野の研究者が参加することが多いですが、検討事項が多いことがその理由だと推察できます。