【R】既存データを用いて分布ラグモデルを動かしてみる

ひとつ前の記事でdlnmパッケージのデータセットを使用しました。今回は、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

今回の狙いは、PM10が呼吸器死亡に与える影響を、ラグ0日～ラグ10日で観察するというものです。つまり、イベント10日前までのPM2.5濃度の平均値をモデルに組み込みます。

まずはパッケージを読み込みます。

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

さて、DLNMはDistributed lag non-linear modelsの略です。ラグを試算するDelayed effectsと、非線形性を記述するNon-linear effectsが含まれます。dlnmパッケージのcrossbasis関数には、この両方の情報を含みます。

PM10の影響

まず各変数の中央値や平均値を眺めておきます。

summary(chicagoNMMAPS)

PM10の中央値はおよそ30ug/m3でした。では次にPM10のラグを含むデータmatrixを作成します。

#crossbasis: dlnm用のmatrixを作る関数
#ラグ0~10日を指定
basis.pm <- crossbasis(chicagoNMMAPS$pm10, lag=10

作成した変数を含むモデルを作成します。

#目的変数を呼吸器死亡、説明変数にPM10と気温
#ポアソン分布を指定
model <- glm(resp ~ basis.pm + temp, family=poisson(), chicagoNMMAPS)

次にcrosspred関数を使用して、ラグによる効果を予測します。

#allRRfit: 全体の累積予測（ラグ0～10までの寄与の合計を試算）
#PM10濃度の中央値30ug/m3を閾値としたときの10ug/m3増加あたりのRR
pred.pm$allRRfit["30"]
#結果
      30 
1.234624 

#allRRlow／allRRhigh: RRの95%信頼区間
cbind(pred.pm$allRRfit, pred.pm$allRRlow, pred.pm$allRRhigh)["30",]
#結果
[1] 1.234624 1.035096 1.472613

中央値30ug/m3は閾値に用いた際、呼吸器死亡が有意に増加すると言えそうです。ではラグ何日のときにその有意差がみられるのか、図示して確認してみます。

#PM10のラグ毎（1日単位）のRR
#閾値の線も表示
plot(pred.pm, var=60, type="p", pch=19, cex=1.5, ci="bars", col=2,
  ylab="RR",main="Lag-specific effects")

ラグによる効果がほぼないことが分かりました。さらに、ラグごとの寄与を合計した全体効果を、横軸を閾値にした作図も行います。

#PM10の閾値毎のRR
#横軸はPM10の閾値
#例：閾値を100ug/m3→RR 1.10
plot(pred.pm, "overall", ci="lines", ylim=c(0.7,1.3), lwd=2, col=4,
  xlab="PM10", ylab="RR", main="Overall effect")

設定する閾値を下げると、RRが下がるという結果になりました。解釈や、コードの正否について、丁寧に確認する必要がありそうです。（筆者の理解の範囲を超えておりました。。。あまり参考にしないでください）

オゾンの影響（チュートリアル）

最後に、チュートリアルにあったオゾンに対する全死亡の検証を再現してみます。

#オゾンはラグ10日
#lag: ラグ、arg-: 因数、fun: 機能
#thr: 閾値を指定、40.3ug/m3
#strata：ダミー変数により自由度dfの区間を定義する
basis.o3 <- crossbasis(chicagoNMMAPS$o3, lag=10, 
  argvar=list(fun="thr",thr.value=40.3,side="h"),
  arglag=list(fun="strata",breaks=c(2,6)))

#チュートリアルの式を使用
#ns: 自然な 3 次スプラインの B スプライン基底行列
#疑似ポアソン分布を指定
model <- glm(death ~ basis.o3 + temp + ns(time,7*14) + dow,
  family=quasipoisson(), chicagoNMMAPS)

#オゾンの値を0～65の正数と、閾値(40.3)、閾値+10（50.3）を指定
#指定した閾値におけるラグ0～10日のRR等が計算される
pred.o3 <- crosspred(basis.o3, model, at=c(0:65,40.3,50.3))

#閾値＋10のときの累積予測RRを計算
cbind(pred.o3$allRRfit, pred.o3$allRRlow, pred.o3$allRRhigh)["50.3",]
#結果
[1] 1.105598 1.054706 1.158946

#オゾンのラグ毎のRR
#閾値の線も表示
plot(pred.o3, var=50.3, type="p", pch=19, cex=1.5, ci="bars", col=2,
  ylab="RR",main="Lag-specific effects")