統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

背景因子で調整した割合の多重比較

Toukei Consul Banner

KH Coder Consul Banner

ある疾患の標準治療不応例に対する追加療法の検討について。

三種類の追加療法を比較検討している(A療法、B療法、C療法)。

質問

エンドポイント(割合)に差があるように見えるが、 その原因としては背景因子のばらつきによると考える。

そこでエンドポイントと関連する重要な背景因子を 一つにまとめたスコアを使って補正することができないかと考えた。

三群間の比較でそのようなことはできるか?

回答

ロジスティック回帰モデルで検討することが可能。

エンドポイントを目的変数にして、三群を説明変数と考えて ロジスティック回帰分析をすれば、背景因子スコアで調整しながら、 A療法に対するB療法のオッズ比、A療法に対するC療法のオッズ比を それぞれ算出することはできる。

これでAとBおよびAとCに差があるのかどうか検討したことになる。

例題

MASSパッケージのAids2というデータを用いて解析してみる。

Aids2はオーストラリアのエイズ患者さんのデータで、 感染経路と生存・死亡との関係を見てみる。

このとき診断時の年齢が感染経路によって異なっているので、 診断時年齢を調整しながら感染経路と生存との関係を見たい。

最初にMASSパッケージを呼び出す。

library(MASS)

単変量解析

まず感染経路と生存・死亡との関係を見てみる。

Aが生存、Dが死亡だ。

行パーセンテージを示している。

感染経路(T.categ)の内訳は以下の通り。

輸血(blood)が最もリスクが高いのがわかる。

> round(with(Aids2, prop.table(table(T.categ, status),1))*100,1)
        status
T.categ     A    D
  hs     37.8 62.2
  hsid   37.5 62.5
  id     60.4 39.6
  het    58.5 41.5
  haem   37.0 63.0
  blood  19.1 80.9
  mother 57.1 42.9
  other  42.9 57.1

年齢との関係はどうか?

死亡した患者さんは平均で1歳年齢が年上だ。

> with(Aids2, t.test(age ~ status))

        Welch Two Sample t-test

data:  age by status
t = -2.6484, df = 2423, p-value = 0.008139
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.7570654 -0.2620655
sample estimates:
mean in group A mean in group D 
       36.78373        37.79330 

感染経路と年齢との関係はどうか?

平均で見てみると、 輸血(blood)によって感染した患者さんの年齢が高いことがわかる。

> round(with(Aids2, tapply(age, T.categ, mean)),1)
    hs   hsid     id    het   haem  blood mother  other 
  37.5   30.5   31.1   38.9   31.9   44.4    3.1   41.8 

となると、輸血の患者さんが感染経路によって死亡リスクが高いのか、 年齢が高いことで死亡リスクが高いのか区別がつかない。

いま検討したいのは感染経路別の死亡リスク。

感染経路別の死亡リスクを年齢で調整しながら確認したい。

共変量を調整した関連性の解析

一般化線形モデルを使う。

glm()という関数だ。

従属変数が死亡か生存という1/0データなので、 ロジスティック回帰分析を使う。

family=binomial()で指定する。

rr.95ci()は対数の結果から真数へ変換する自作の関数だ。

logistic.res <- glm(status ~ T.categ + age, data=Aids2, family=binomial())
summary(logistic.res)

rr.95ci <- function(x){
res <- exp(cbind(coef(x),confint(x)))
colnames(res) <- c("RR","LL(2.5%)","UL(97.5%)")
round(res,2)
}
rr.95ci(logistic.res)

感染経路の結果は、ホモセクシュアル(hs)と比べて、 静脈注射薬物使用者(id)、ヘテロセクシュアル(het)の患者さんは有意に死亡リスクが低く、 輸血(blood)から感染した患者さんは有意に死亡リスクが高いというものだ。

年齢を調整しても、輸血による感染患者さんはリスクが高いといえる。

> logistic.res <- glm(status ~ T.categ + age, data=Aids2, family=binomial)
> summary(logistic.res)

Call:
glm(formula = status ~ T.categ + age, family = binomial, data = Aids2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9219  -1.3712   0.9492   0.9828   1.4076  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    0.199379   0.159132   1.253 0.210237    
T.categhsid    0.070244   0.248652   0.282 0.777561    
T.categid     -0.868489   0.299245  -2.902 0.003705 ** 
T.categhet    -0.853252   0.320149  -2.665 0.007695 ** 
T.categhaem    0.084695   0.309908   0.273 0.784630    
T.categblood   0.897524   0.266752   3.365 0.000766 ***
T.categmother -0.511978   0.777751  -0.658 0.510358    
T.categother  -0.242098   0.246017  -0.984 0.325082    
age            0.007920   0.004111   1.926 0.054046 .  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3777.5  on 2842  degrees of freedom
Residual deviance: 3739.4  on 2834  degrees of freedom
AIC: 3757.4

Number of Fisher Scoring iterations: 4


> rr.95ci(logistic.res)
Waiting for profiling to be done...
                RR LL(2.5%) UL(97.5%)
(Intercept)   1.22     0.89      1.67
T.categhsid   1.07     0.66      1.77
T.categid     0.42     0.23      0.75
T.categhet    0.43     0.22      0.79
T.categhaem   1.09     0.60      2.04
T.categblood  2.45     1.49      4.26
T.categmother 0.60     0.12      2.79
T.categother  0.78     0.49      1.28
age           1.01     1.00      1.02

共変量調整多重比較

では、感染経路のうちどれがもっともリスクが高いのか、 年齢を調整しながら、総当たりのTukeyの方法で検討してみよう。

multcompパッケージのglht()関数を使う。

library(multcomp)
glht.res <- glht(logistic.res, linfct = mcp(T.categ = "Tukey"))
summary(glht.res)

感染経路は8つあるので、 8つから2つを取り出す組み合わせは28通り。

組み合わせの計算は以下の通り。

 \displaystyle {}_8 \mathrm{ C }_2 = \frac{n!}{k! \times (n-k)!} = \frac{8!}{2! \times (8-2)!} = \frac{8 \times 7}{2 \times 1} = 28

統計学的有意な組み合わせは以下の4通り。

輸血から感染した患者さんが ホモセクシュアル、薬物使用者、ヘテロセクシュアル、 その他の感染経路の人からの患者さんに比べて 死亡するリスクが高いということがわかる。

blood - hs == 0 0.89752 0.26675 3.365 0.0144 *

blood - id == 0 1.76601 0.39787 4.439 <0.001 ***

blood - het == 0 1.75078 0.41220 4.247 <0.001 ***

other - blood == 0 -1.13962 0.35710 -3.191 0.0255 *

> summary(glht.res) 

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = status ~ T.categ + age, family = binomial, data = Aids2)

Linear Hypotheses:
                    Estimate Std. Error z value Pr(>|z|)    
hsid - hs == 0       0.07024    0.24865   0.282   1.0000    
id - hs == 0        -0.86849    0.29925  -2.902   0.0606 .  
het - hs == 0       -0.85325    0.32015  -2.665   0.1133    
haem - hs == 0       0.08469    0.30991   0.273   1.0000    
blood - hs == 0      0.89752    0.26675   3.365   0.0144 *  
mother - hs == 0    -0.51198    0.77775  -0.658   0.9972    
other - hs == 0     -0.24210    0.24602  -0.984   0.9699    
id - hsid == 0      -0.93873    0.38270  -2.453   0.1869    
het - hsid == 0     -0.92350    0.40159  -2.300   0.2578    
haem - hsid == 0     0.01445    0.39113   0.037   1.0000    
blood - hsid == 0    0.82728    0.36152   2.288   0.2640    
mother - hsid == 0  -0.58222    0.80947  -0.719   0.9952    
other - hsid == 0   -0.31234    0.34626  -0.902   0.9815    
het - id == 0        0.01524    0.43469   0.035   1.0000    
haem - id == 0       0.95318    0.42526   2.241   0.2889    
blood - id == 0      1.76601    0.39787   4.439   <0.001 ***
mother - id == 0     0.35651    0.82690   0.431   0.9998    
other - id == 0      0.62639    0.38407   1.631   0.6910    
haem - het == 0      0.93795    0.44209   2.122   0.3585    
blood - het == 0     1.75078    0.41220   4.247   <0.001 ***
mother - het == 0    0.34127    0.84021   0.406   0.9999    
other - het == 0     0.61115    0.39915   1.531   0.7547    
blood - haem == 0    0.81283    0.40592   2.002   0.4349    
mother - haem == 0  -0.59667    0.83091  -0.718   0.9952    
other - haem == 0   -0.32679    0.39241  -0.833   0.9884    
mother - blood == 0 -1.40950    0.82362  -1.711   0.6364    
other - blood == 0  -1.13962    0.35710  -3.191   0.0255 *  
other - mother == 0  0.26988    0.81662   0.330   1.0000    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Warning message:
In RET$pfunction("adjusted", ...) : Completion with error > abseps

95%信頼区間を求めることもできる。

> confint(glht.res)

         Simultaneous Confidence Intervals

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = status ~ T.categ + age, family = binomial, data = Aids2)

Quantile = 2.9706
95% family-wise confidence level
 

Linear Hypotheses:
                    Estimate lwr      upr     
hsid - hs == 0       0.07024 -0.66840  0.80888
id - hs == 0        -0.86849 -1.75742  0.02044
het - hs == 0       -0.85325 -1.80428  0.09777
haem - hs == 0       0.08469 -0.83591  1.00530
blood - hs == 0      0.89752  0.10512  1.68993
mother - hs == 0    -0.51198 -2.82234  1.79839
other - hs == 0     -0.24210 -0.97291  0.48871
id - hsid == 0      -0.93873 -2.07556  0.19810
het - hsid == 0     -0.92350 -2.11646  0.26947
haem - hsid == 0     0.01445 -1.14743  1.17634
blood - hsid == 0    0.82728 -0.24664  1.90120
mother - hsid == 0  -0.58222 -2.98681  1.82237
other - hsid == 0   -0.31234 -1.34092  0.71623
het - id == 0        0.01524 -1.27605  1.30652
haem - id == 0       0.95318 -0.31009  2.21646
blood - id == 0      1.76601  0.58412  2.94791
mother - id == 0     0.35651 -2.09986  2.81289
other - id == 0      0.62639 -0.51453  1.76731
haem - het == 0      0.93795 -0.37532  2.25121
blood - het == 0     1.75078  0.52631  2.97524
mother - het == 0    0.34127 -2.15462  2.83717
other - het == 0     0.61115 -0.57454  1.79685
blood - haem == 0    0.81283 -0.39299  2.01865
mother - haem == 0  -0.59667 -3.06496  1.87162
other - haem == 0   -0.32679 -1.49249  0.83890
mother - blood == 0 -1.40950 -3.85611  1.03710
other - blood == 0  -1.13962 -2.20040 -0.07884
other - mother == 0  0.26988 -2.15594  2.69570

結果を真数に変換すればオッズ比として表現することもできる。

スクリプトがわからずGoogle Sheetsの助けを借りた。

太字下線の文字・数値のセルが統計学的有意に異なる組み合わせ。

組み合わせ 対数 Estimate lwr upr 真数 Estimate lwr upr
hsid - hs 0.07024 -0.66816 0.80865 1.07 0.51 2.24
id - hs -0.86849 -1.75714 0.02016 0.42 0.17 1.02
het - hs -0.85325 -1.80398 0.09748 0.43 0.16 1.10
haem - hs 0.08469 -0.83562 1.00501 1.09 0.43 2.73
blood - hs 0.89752 0.10537 1.68968 2.45 1.11 5.42
mother - hs -0.51198 -2.82162 1.79766 0.60 0.06 6.04
other - hs -0.2421 -0.97268 0.48848 0.78 0.38 1.63
id - hsid -0.93873 -2.07521 0.19774 0.39 0.13 1.22
het - hsid -0.9235 -2.11609 0.2691 0.40 0.12 1.31
haem - hsid 0.01445 -1.14707 1.17597 1.01 0.32 3.24
blood - hsid 0.82728 -0.2463 1.90086 2.29 0.78 6.69
mother - hsid -0.58222 -2.98606 1.82161 0.56 0.05 6.18
other - hsid -0.31234 -1.34059 0.71591 0.73 0.26 2.05
het - id 0.01524 -1.27564 1.30612 1.02 0.28 3.69
haem - id 0.95318 -0.30969 2.21606 2.59 0.73 9.17
blood - id 1.76601 0.58449 2.94754 5.85 1.79 19.06
mother - id 0.35651 -2.09909 2.81212 1.43 0.12 16.65
other - id 0.62639 -0.51417 1.76695 1.87 0.60 5.85
haem - het 0.93795 -0.3749 2.2508 2.55 0.69 9.50
blood - het 1.75078 0.5267 2.97486 5.76 1.69 19.59
mother - het 0.34127 -2.15384 2.83639 1.41 0.12 17.05
other - het 0.61115 -0.57417 1.79648 1.84 0.56 6.03
blood - haem 0.81283 -0.39261 2.01827 2.25 0.68 7.53
mother - haem -0.59667 -3.06419 1.87084 0.55 0.05 6.49
other - haem -0.32679 -1.49212 0.83854 0.72 0.22 2.31
mother - blood -1.4095 -3.85534 1.03634 0.24 0.02 2.82
other - blood -1.13962 -2.20007 -0.07918 0.32 0.11 0.92
other - mother 0.26988 -2.15518 2.69494 1.31 0.12 14.80

共変量調整多重比較(Bonferroni型p値調整)

Bonferroni型のp値調整、例えばHochbergの方法でもp値調整ができる。

Bonferroni型のp値調整の詳細はこちら。

toukeier.hatenablog.com

summary()で、test=adjusted("hochberg")と指定する。

Tukeyのように等分散性の仮定が不要な分、 より適用範囲が広く、適切な場面が多い。

> glht.res <- glht(logistic.res, linfct = mcp(T.categ = "Tukey"))
> summary(glht.res, test = adjusted("hochberg"))

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = status ~ T.categ + age, family = binomial(), data = Aids2)

Linear Hypotheses:
                    Estimate Std. Error z value Pr(>|z|)    
hsid - hs == 0       0.07024    0.24865   0.282 0.972038    
id - hs == 0        -0.86849    0.29925  -2.902 0.088914 .  
het - hs == 0       -0.85325    0.32015  -2.665 0.176984    
haem - hs == 0       0.08469    0.30991   0.273 0.972038    
blood - hs == 0      0.89752    0.26675   3.365 0.019927 *  
mother - hs == 0    -0.51198    0.77775  -0.658 0.972038    
other - hs == 0     -0.24210    0.24602  -0.984 0.972038    
id - hsid == 0      -0.93873    0.38270  -2.453 0.311732    
het - hsid == 0     -0.92350    0.40159  -2.300 0.442353    
haem - hsid == 0     0.01445    0.39113   0.037 0.972038    
blood - hsid == 0    0.82728    0.36152   2.288 0.442353    
mother - hsid == 0  -0.58222    0.80947  -0.719 0.972038    
other - hsid == 0   -0.31234    0.34626  -0.902 0.972038    
het - id == 0        0.01524    0.43469   0.035 0.972038    
haem - id == 0       0.95318    0.42526   2.241 0.475002    
blood - id == 0      1.76601    0.39787   4.439 0.000253 ***
mother - id == 0     0.35651    0.82690   0.431 0.972038    
other - id == 0      0.62639    0.38407   1.631 0.972038    
haem - het == 0      0.93795    0.44209   2.122 0.609664    
blood - het == 0     1.75078    0.41220   4.247 0.000584 ***
mother - het == 0    0.34127    0.84021   0.406 0.972038    
other - het == 0     0.61115    0.39915   1.531 0.972038    
blood - haem == 0    0.81283    0.40592   2.002 0.769053    
mother - haem == 0  -0.59667    0.83091  -0.718 0.972038    
other - haem == 0   -0.32679    0.39241  -0.833 0.972038    
mother - blood == 0 -1.40950    0.82362  -1.711 0.972038    
other - blood == 0  -1.13962    0.35710  -3.191 0.035401 *  
other - mother == 0  0.26988    0.81662   0.330 0.972038    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- hochberg method)

まとめ

割合の多重比較をする場合、 背景因子が異なっていて調整したいときは、 統計ソフトRなら、 glm(family=binomial())でロジスティック回帰を行い、 結果に対して、multcompパッケージのglht()関数を使う。

この二段構えで、解決できる。