メタ解析論文を読んでいるときに見かける、異質性と均質性という言葉について解説したいと思います。
異質性は英語ではHeterogeneityと言います。メタ解析では複数の文献の、複数の結果をまとめて図表として表記します。オッズ比などの比の形で表されることが多く、曝露の度合いをできるだけ合わせて並べられます。その時に、結果の正負の方向がばらばらで一貫性がない場合、曝露による効果に異質性がみられる、と言います。
均質性は英語ではHomogeneityと言います。異質性とは逆で、正負の方向が一致していたり、関連がみられなかったりと、結果が似通っている場合には、曝露による均質性がみられる、と言います。
異質性と均質性とは、メタ解析の結果に対して使う定性的な表現方法でした。実は、統計的に異質性を評価する方法として、Cochran’s Q、Higgins ＆Thompson’s I2統計量、H2統計量、異質性分散τ2などがあります。これらについては後日、取り扱いたいと思います。
今回はRで、既存のデータセットを使用して、オッズ比計算をデモンストレーションしてみます。使用するのはEpiパッケージのdietデータセットです。
このデータは、冠状動脈性心疾患イベントと食事の関係性を評価するためのデータセットです。
まずはパッケージをインポートします。
library(tidyverse)
library(Epi)
library(broom)dietデータセットをインポートします。
#オリジナルデータはClayton and Hills (1993).
data(diet)
head(diet)
#id: 被験者識別子、数値ベクトル。
#doe: フォローアップ研究への参加日、日付変数。
#dox：フォローアップ研究からの離脱日、日付変数。
#dob：生年月日、日付変数。
#y: 危険年数、数値ベクトル。(doxとdoeの差)
#fail: 終了時の状態、数値ベクトル（コード1,3,13はCHDイベントを表す）。
#job: 職業。レベル Driver Conductor Bank worker を持つ因子。
#month：食事調査の月、数値ベクトル。
#energy：総エネルギー摂取量（kCal/日/100）、数値ベクトル。
#height：（cm）、数値ベクトル。
#weight：(kg)、数値ベクトル。
#fat：脂肪摂取量（10g/日）、数値ベクトル。
#fibre：食物繊維の摂取量（10g/日）、数値ベクトル。
#energy.grp：1日のエネルギー摂取量、レベル<=2750 kCal, >2750 kCalの因子。
#chd：冠状動脈性心疾患イベント、数値ベクトル（1=CHDイベント、0=イベントなし）。平均値を見た後に、全体を可視化してみます。
summary(diet)
#全体像を眺める
di <- select(diet, y, energy, height, weight, fat, fibre, chd)
plot(di)
energy~fibreは相関がみられることは容易に想像できますね。chdについては、傾向は見づらいです。
では、モデルを作成します。変数chdは0/1なので、二項分布を仮定したモデルを組みます。モデルには職業を表す変数job（driver, conductor, bank worker）も含めることにしました。（※本データはコホート研究なので、オッズ比ではなく相対リスクが求められる気もしますが、今回はORとします）
#ロジット-二項モデル
#説明変数に脂肪、食物繊維、摂取エネルギー、体重、調査月、職業
fit <- glm(chd ~ fat + fibre + energy + weight + month + job + y,
data = diet,
family = binomial(link="logit"))
#結果
summary(fit)
Call:
glm(formula = chd ~ fat + fibre + energy + weight + month + job +
y, family = binomial(link = "logit"), data = diet)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0752 -0.4211 -0.2567 -0.1499 2.9383
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.41407 1.78929 1.349 0.17728
fat -0.25932 0.15530 -1.670 0.09496 .
fibre -0.99924 0.57815 -1.728 0.08393 .
energy 0.09327 0.09391 0.993 0.32064
weight 0.01472 0.02108 0.698 0.48509
month -0.05957 0.05790 -1.029 0.30357
jobConductor 1.03778 0.54414 1.907 0.05650 .
jobBank worker 1.55917 0.54841 2.843 0.00447 **
y -0.32176 0.04614 -6.973 3.1e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 258.88 on 328 degrees of freedom
Residual deviance: 174.04 on 320 degrees of freedom
( 8 個の観測値が欠損のため削除されました )
AIC: 192.04
Number of Fisher Scoring iterations: 6p値で有意なのはjob（Bank worker）とy（危険年数≒観察期間）であることが分かりました。偏回帰係数estimateでは読み取りづらいため、指数変換してオッズ比を計算します。
#偏回帰係数estimateを指数変換して調整オッズ比を計算する
coef(summary(fit)) %>% as.data.frame() %>%
mutate(OR = exp(Estimate),
conf.low = exp(Estimate - 1.96*`Std. Error`),
conf.high = exp(Estimate + 1.96*`Std. Error`)) %>%
dplyr::select(OR, conf.low, conf.high)
このままでもよいのですが、あとで図にしたいのでdata.frame形式として出力するためbloomパッケージを使います。
#bloomパッケージのtidyを使う
#,c(1,2,6,7)によりORと95%信頼区間のみ表示する
OR <- tidy(fit, exponentiate = T, conf.int = T)[,c(1,2,6,7)]
OR
微妙に結果が異なるのは、信頼区間の計算に用いた数値(1.96)の桁数の違いと考えられます。では、この結果をフォレストプロットとして図示します。
#フォレストプロットで表示
forest <- OR %>%
#intercept情報を削除
filter(term != "(Intercept)") %>%
#y軸の順番を指定する
mutate(term = fct_relevel(
term,
"fat","fibre","energy","weight","month","jobConductor","jobBank worker","y")) %>%
#グラフの構造説明
ggplot(aes(x=estimate, y=term))+
#estimateを点表示
geom_point()+
#信頼区間をエラーバーとして表示
geom_errorbar(aes(xmin=conf.low, xmax=conf.high))+
#RR=1を示す線を表示
geom_vline(xintercept = 1, linetype="dashed")+
#軸にラベル付け
labs(x = "OR", y = "factors")
forest
結果が出ました。冠状動脈性心疾患と脂肪、食物繊維、摂取エネルギー、体重、調査月に関連はみられないことがわかりました。また、職業がDriverの人に対して、ConductorとBank workerは冠状動脈性心疾患のリスクが高いと言えることも見て取れます。
今回は、既存データを用いてオッズ比の計算を行ってみました。参考になれば幸いです。
疫学のデータを扱うとき、もしくは文献を読むときに、データが偏りや交絡の影響を受けていないか注意して解釈する必要があります。
偏りとは、観察された結果が真実の姿からある方向にずれている状態を表します。偏りには2種類あります。
1つ目は、選択バイアスselection biasといい、調査対象とする集団の選定に偏りが生じる場合です。例えば、症例群を病院に来る患者、対照群を来院者ではない一般市民から選ぶ場合、健康への関心度や既往歴などに集団と差があることが考えられる。
2つ目は、情報バイアスinformation biasであり、アウトカム（疾病発症、検査結果、インタビューなど）の測定結果に偏りが生じることです。思い出しバイアスrecall bias（過去に関する記憶の度合いが症例群と対照群で異なる）、質問者バイアスinterviewer bias（質問をする担当者の聞き方が症例群と対照群で異なる）などがある。
こうした偏りは、データ収集後の解析段階では、制御することはできないことが特徴です。
交絡とは、アウトカムと曝露の関係に、他の要因が影響を及ぼすことを言い、その因子を交絡因子confounding factorと呼びます。交絡バイアスという呼び方をされることもあります。
交絡は偏りとは異なり、解析段階で制御することができます（これを調整と呼びます）。
偏りにはいくつか種類があり、交絡も偏りの一つとして捉えられることがあります。偏りは調整できない、交絡は調整可能と覚えるもの手かと思います。
疫学において、
のような事象が見られることがあります。実際には関連がないのに、分析結果で関連があるようにみえることを、交絡と呼びます。
調べていくと、「交絡因子」と「調整因子」で呼び方が異なることがあります。違いはあるのでしょうか。
冒頭の内容と被りますが、アウトカム（疾病発症や検査結果など）と曝露の関係に影響を及ぼすような別の影響のことを、交絡因子confounding factorと呼びます。疫学研究において、交絡因子を完全に排除することは難しいですが、可能な限りは交絡因子を検討・考慮する必要があります。ここでは詳しくは触れませんが、交絡因子の対処方法が適切である研究方法の方が、疫学研究としての妥当性が高い（より信頼性の高いデータとして認識される）と言われます。妥当性について大雑把には、ランダム化比較試験＞コホート研究＞症例対照研究＞生態学的研究というような順番です。
交絡の要因になる因子に見当がついた場合、その影響を取り除いて解析する必要があり、これを調整adjustmentと呼びます。
調整の方法にはいくつかあります。解析対象を制限する（例：禁煙者に絞る）、マッチングする（例：症例と対照で性別・年齢などの要素を合わせて対象者を選ぶ）、階層化する（例：年齢を10代ごとにわける）、重み付けする（例：症例と対照で喫煙者の割合が合うように数字を補正する）、モデルに組み込む（大変量解析モデルの因子に加えることで、因子ごとの影響をみる）などです。
交絡の因果関係を考えるうえで、DAG（Directed Acyclic Graph）と呼ばれる方法があります。これはたくさんの変数がある場合に、変数間の関係を医学的な知見から矢印で方向をつけて可視化することで、何で調整すべきか（過剰調整もよくない）検討するものです。
生物系の分野だと、関連のありそうな自然のデータを全部組み込んで、AIC（赤池情報基準）などを用いてベストモデルを決めるということをします。疫学の分野ではAICを使う文献は少なく、医学的知見から著者らが研究ごとに調整因子を決めているという印象です。
考えられる交絡因子の中から、測定可能かつ関連がありそうな因子を調整する（調整因子）という関係性であることがわかりました。文献を読むうえでは、（検討した）交絡因子と調整因子はほとんど同じものと考えてもよさそうです。
ひとつ前の記事で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関数には、この両方の情報を含みます。
まず各変数の中央値や平均値を眺めておきます。
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")
ラグ0、１日のときにRRが高いこと、全死亡に対する正の関連はラグ5日までは維持されることが見て取れます。
dlnmパッケージを使用して、ラグモデルを動かしてみました。正直、筆者にはまだ難易度が高く、チュートリアル以外のことに手を広げるには、データそのものへの理解とともに、関数のhelpをもっと熟読する必要がありそうです。
公開されている疫学文献データを用いて、相対リスクの計算を行ってみます。主な目的は２つです。
・データに対する分布を検討する
・調整モデル、複数汚染物質モデルなどを検討する
使用するパッケージは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-16p<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()を用いることにします。
では線形モデルを仮定して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.1093132PM10やオゾンに正の関連があることがわかりました。見やすくするために、図示します。使用するのはフォレストプロットという、各推定値・信頼区間を縦に並べて表示するグラフです。
#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のコードに当てはめて計算するのは簡単ですが、データの特性を考えて分布を当てはめたり、調整因子の入れ方を検討したりする必要があることが分かります。疫学研究は様々な専門分野の研究者が参加することが多いですが、検討事項が多いことがその理由だと推察できます。
今回は、適当に作成したデータを用いて、記述統計に用いる要約統計量をまとめた表と、可視化の方法として箱ひげ図を作成してみます。
方針としては、子供の肺機能と身長・体重、既往歴などを含めたデータを作成して、学年と性別ごとに並べて見比べることにします。使用するＲパッケージは、"モダンな"データ形成ができるtidyverseと、直感的な作図ができるesquisseです。
現実ではきれいなデータはあり得ませんが、今回はテストなので性比や身長、FVCの値は適当に指定していきます。
#libraryインポート
library(tidyverse)
library(esquisse)
#乱数種の指定
set.seed(1)
#条件
#中学１～３年生、100人ずつ
#性別 male/female
#身長 130～180cmくらい
#体重 40～60kgくらい
#肺機能として努力性肺活量（FVC）3.5～5.5Lくらい
#病気Aの既往歴 yes/no
#1年生 100人
id1 <- c(1:100)
grad1 <- rep(1,100)
#今回は性別はランダムに割り振ります（男女比半々）
s <- c("male", "female")
sex1 <- sample(s, 100, replace=T)
#平均身長120cm、標準偏差20とします
height1 <- rnorm(n=100, mean =120, sd=20)
#平均体重45kg、標準偏差10とします
weight1 <- rnorm(n=100, mean =45, sd=10)
#平均2.5L、標準偏差0.3とします
FVC1 <- rnorm(n=100, mean =2.5, sd=0.3)
#病気Aの既往歴ありの割合は20％とします
a <- c("yes", "no")
sickA1 <- sample(a, 100, replace=T, prob=c(0.2, 0.8))
G1 <- data.frame(id1, grad1, sex1, height1, weight1, FVC1, sickA1)sample()関数は、乱数を発生させランダムな割付を行う関数です。結果が毎回変わるため、今回はset.seed(1)を指定して、結果の再現性があるようにしました。他2学年については以下のようにしました。
#2年生 100人
id2 <- c(101:200)
grad2 <- rep(2,100)
#今回は性別はランダムに割り振ります（男女比半々）
sex2 <- sample(s, 100, replace=T)
#平均身長135cm、標準偏差20とします
height2 <- rnorm(n=100, mean =135, sd=20)
#平均体重48kg、標準偏差10とします
weight2 <- rnorm(n=100, mean =48, sd=10)
#平均3.3L、標準偏差0.4とします
FVC2 <- rnorm(n=100, mean =3.3, sd=0.4)
#病気Aの既往歴ありの割合は25％とします
sickA2 <- sample(a, 100, replace=T, prob=c(0.25, 0.8))
G2 <- data.frame(id2, grad2, sex2, height2, weight2, FVC2, sickA2)
#3年生 100人
id3 <- c(201:300)
grad3 <- rep(3,100)
#今回は性別はランダムに割り振ります（男女比半々）
sex3 <- sample(s, 100, replace=T)
#平均身長150cm、標準偏差20とします
height3 <- rnorm(n=100, mean =150, sd=20)
#平均体重52kg、標準偏差10とします
weight3 <- rnorm(n=100, mean =52, sd=10)
#平均4.0L、標準偏差0.5とします
FVC3 <- rnorm(n=100, mean =4.0, sd=0.5)
#病気Aの既往歴ありの割合は30％とします
sickA3 <- sample(a, 100, replace=T, prob=c(0.30, 0.8))
G3 <- data.frame(id3, grad3, sex3, height3, weight3, FVC3, sickA3)最後にこれらのデータをまとめて一つのデータフレームにします。その際に、列名が異なると結合できないため、同名に変更しました。
colnames(G1) <- c("ID","grade","sex","height","weight","FVC","medhis")
colnames(G2) <- c("ID","grade","sex","height","weight","FVC","medhis")
colnames(G3) <- c("ID","grade","sex","height","weight","FVC","medhis")
DF <- rbind(G1, G2, G3)これで仮データの作成が完了しました。いよいよ次に、記述統計表（要約統計量の表）を作成します。
まずはデータの全体像を眺めてみます。見やすくはないですが、基本のベースRのplot()を使います。
#全体像の確認
plot(DF)
学年（2列目）を見ると、学年が上がるごとに身長・体重・FVCが上がっていることが読み取れます（そういう風にデータ形成したので当たり前ですが）。
次にデータの数値的な様子を読み解くために、統計量を出してみます。まずはベースRのsummary()を使います。
#全体の統計値
summary(DF) ID 学年 性別
Min. : 1.00 Min. :1 Length:300
1st Qu.: 75.75 1st Qu.:1 Class :character
Median :150.50 Median :2 Mode :character
Mean :150.50 Mean :2
3rd Qu.:225.25 3rd Qu.:3
Max. :300.00 Max. :3
身長 体重 FVC
Min. : 81.71 Min. :15.87 Min. :1.800
1st Qu.:115.70 1st Qu.:41.68 1st Qu.:2.700
Median :134.16 Median :48.62 Median :3.243
Mean :134.88 Mean :48.26 Mean :3.275
3rd Qu.:151.95 3rd Qu.:55.70 3rd Qu.:3.766
Max. :193.79 Max. :74.06 Max. :5.115
既往歴
Length:300
Class :character
Mode :character これだけ見ても仕方ないですよね・・・。学年や性別などで分けてみたいと思います。パッケージtidyverseに含まれている、ggplot()を用います。パイプ演算子%>%や、使われている関数（dplyrの関数）の詳細説明は、他に譲らせていただきます。
summary <-
DF %>%
group_by(grade) %>%
summarize(
No. = n(), #総数
Male = sum(sex=="male"), #男女比として男性を数える
Height = mean(height), #平均身長
Weight =mean(weight), #平均体重
FVC = mean(FVC), #平均FVC
MedHis = sum(medhis=="yes") #病歴ありの数
)
summary
学年ごとの性別比（男性数）、身長、体重、FVCの平均値が一つの表に表示されました。
では次にグラフ化してみます。今回、学年が数字numericで、このままではグラフにできないため文字列characterへ事前に変換しておきます。
#numericからcharacterへ変換
DF2 <- DF %>% mutate(grade = as.character(grade))グラフ作成のポイントですが、どんな図を作りたいかあらかじめ考えておくことです。それを直感的操作で手助けしてくれるのが、esquisse()です。動かしてみる便利さがわかります。
#作図できるパターンを直感的に理解する
esquisser(DF2)
上記の画面が表示されるので、ｘ軸にgrade、ｙ軸にweightをドラッグ＆ドロップします。右下のCode▲を押すと、表示されているグラフのコードが表示されます。
このコードをそのまま使ってもいいですし、これをベースに書き換えていくのでもよいです。ではまずは学年別身長をグラフにしてみます。
#学年別身長
g_height <- ggplot(DF2)+
geom_boxplot(aes(x=grade, y=height))
g_height
次に学年別体重を、男女別に表示してみます。
#学年別体重、男女別
g_weight <- ggplot(DF2)+
geom_boxplot(aes(x=grade, y=weight, fill=sex))
g_weight
では最後に、学年別FVCを、男女別、かつ既往歴有無をグループ分けして表示します。
#学年別FVC,男女別,既往歴で分ける
g_FVC <- ggplot(DF2)+
geom_boxplot(aes(x=grade, y=FVC, fill=sex))+
facet_wrap(vars(medhis))
g_FVC
今回は記述統計表とグラフ作成の練習を紹介しました。今回はデータを適当に作成したため、性別や既往歴で違いは見れませんでしたが、図表作成の練習にはなるかと思いますので、ぜひ練習してみてください。
