【第765回】 An example of multilevel analysis was revisited（2014年8月21日）
- 6:00起床。今日で50歳になってしまった。豚肉を醤油と味醂とおろし生姜を混ぜた液に暫く浸けてから，くし切りにしたタマネギとともにオリーブオイルで炒め，水で溶いた片栗粉をかけ，ミニトマトをいくつか載せてから蓋をして弱火で少し炒め続けたものとご飯で朝食を済ませ（半分残ったが，しっかり火を通してフライパンに蓋をしておけば，今夜くらいまでなら十分保つだろう），ゴミ出しをしてから往路4番のバスで出勤。
- 昨日のマルチレベル分析の例だが，データの性状を見るための作図からマルチレベル分析まで，一通りの操作をするコードと作図と出力をまとめて示しておく。
- まず，介入の有無も施設の違いも無視して，前後でのコレステロール値の差があるかどうかだけを対応のあるt検定で調べてみると，
と，有意水準5%で統計学的に有意な差があるといえる。Paired t-test
data: x$cholesterol_base and x$cholesterol
t = 8.0383, df = 275, p-value = 2.713e-14
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
10.96801 18.08272
sample estimates:
mean of the differences
14.52536 
- 藤野善久・近藤尚己・竹内文乃『保健医療従事者のためのマルチレベル分析』診断と治療社，ISBN 978-4-7878-2053-2（Amazon | honto | e-hon）のp.69図9-4は，Stataによる回帰分析の結果である。ここでは，まず形式的に，介入の有無とベースラインのコレステロール値の，治療後のコレステロール値への交互作用効果を調べてみた。すると，
と，交互作用効果（cholesterol_base:as.factor(intervention)1の行）も有意水準5%で統計学的に有意であった。つまり，介入した群としなかった群の間で，ベースラインのコレステロール値と治療後のコレステロール値の関係が有意に異なっていた。もし解析の目的がベースラインコレステロール値と治療後コレステロール値の関係を調べることであれば，共分散分析によって修正平均を比較するよりも，2群別々に分析すべきということになるのだが，ここでの分析の目的は介入効果の評価なので，その方法では不適切である。Call:
lm(formula = cholesterol ~ cholesterol_base * as.factor(intervention), data = x)
Residuals:Min 1Q Median 3Q Max -52.619 -18.290 -2.475 16.547 87.592 Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 41.05775 12.48988 3.287 0.00114 ** cholesterol_base 0.78941 0.05094 15.496 < 2e-16 *** as.factor(intervention)1 111.80484 24.42985 4.577 7.19e-06 *** cholesterol_base:as.factor(intervention)1 -0.48289 0.09812 -4.922 1.49e-06 *** ---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 25.97 on 272 degrees of freedom
Multiple R-squared: 0.483, Adjusted R-squared: 0.4773
F-statistic: 84.72 on 3 and 272 DF, p-value: < 2.2e-16 - そこで，図9-4と同じ結果を得るには交互作用効果を入れずにlm()関数で線形回帰すれば良い。
のようになる。Call:
lm(formula = cholesterol ~ cholesterol_base + as.factor(intervention), data = x)
Residuals:Min 1Q Median 3Q Max -51.755 -20.040 -2.563 16.194 91.172 Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 72.49400 11.18021 6.484 4.16e-10 *** cholesterol_base 0.65923 0.04535 14.536 < 2e-16 *** as.factor(intervention)1 -7.42741 3.27763 -2.266 0.0242 * ---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 27.05 on 273 degrees of freedom
Multiple R-squared: 0.437, Adjusted R-squared: 0.4329
F-statistic: 106 on 2 and 273 DF, p-value: < 2.2e-16 
- しかし施設ごとに色を変えてプロットし直してみると，施設の効果がありそうに思えるので，これはマルチレベル分析にすべきである。lmerTestパッケージのlmer()関数により得られた結果オブジェクトをsummary()に渡せば，
という結果が得られる。Linear mixed model fit by maximum likelihood ['merModLmerTest']
Formula: cholesterol ~ cholesterol_base + intervention + (1 | clinic)
Data: xAIC BIC logLik deviance df.resid 2594.9 2613.0 -1292.4 2584.9 271 Scaled residuals:
Min 1Q Median 3Q Max -1.8741 -0.7486 -0.1038 0.5595 3.0095 Random effects:
Groups Name Variance Std.Dev. clinic (Intercept) 141.7 11.90 Residual 637.0 25.24 Number of obs: 276, groups: clinic, 10
Fixed effects:Estimate Std. Error df t value Pr(>|t|) (Intercept) 116.72600 13.27558 163.48000 8.793 2e-15 *** cholesterol_base 0.47132 0.05239 247.28000 8.997 <2e-16 *** intervention -1.74572 3.26249 275.74000 -0.535 0.593 ---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:(Intr) chlst_ cholstrl_bs -0.945 interventin 0.120 -0.241 - 係数の信頼区間は，同じ結果オブジェクトをconfint()関数に渡せば得られる。昨日書いた通り，以上の手順で，Stataで得られる結果はほぼすべてRでも得られたと思う。本にも書かれていたが，マルチレベル分析にすると，回帰分析では有意であったinterventionの効果が有意でない，というのがポイントである。この場合，介入効果があったように見えていたのは，施設間差によるartefactであったと考えられる。
Computing profile confidence intervals ...
2.5 % 97.5 % .sig01 6.4527670 21.7732086 .sigma 23.2278626 27.5583053 (Intercept) 87.9180189 145.5631716 cholesterol_base 0.3581452 0.5886708 intervention -8.4035413 4.8428531 - このREUTERSのtweetも，リンク先記事でも「最大3万人分の治療薬やワクチンが必要」と書いているが，Oliver BradyのNatureの記事と，そこからリンクされている，彼らが推定に使ったExcelの表を見ると，予防・治療用に約3万人分の薬が必要ということだ。患者数としては8月19日時点で2,240人という実数を計算の出発点にしており，患者とコンタクトして薬の投与が必要になる人がその9倍（患者1人当たり平均9人と接触という）という推定になっている。この値は"conservative"なシナリオだそうだが，Excelの表の上でスライドバーを動かしてパラメータを変えられるようなので，簡単な感度分析は試せそうだ。
- 成績処理を進め，終バスで帰宅し，朝のおかずの残りを温め直してレトルトご飯とともに晩飯とした。ドラゴンズはあと１点追いつけず阪神に３タテを食らってしまった。
△Read/Write COMMENTS
▼前【764】（今日も成績処理など（2014年8月20日） ） ▲次【766】（まだまだ成績処理（2014年8月22日） ） ●Top