Cox の比例ハザードモデルに多重共線性はあるか

Last modified: Oct 19, 2015

 重回帰分析では,独立変数間の相関が高すぎると多重共線性の問題が生じるということは広く知られているところである。
 では,Cox の比例ハザードモデルにおいては多重共線性はあるだろうか?という疑問に答えよう。


1 テストデータ

 以下のようなデータを使う。

test <- data.frame(time=c(4, 3, 1, 1, 2, 2, 3, 1, 2, 3, 2, 1, 1),  
              status=c(1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0),  
              sex=c(0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1))  
set.seed(123) 
test$x <- rnorm(13) 
test$y <- rnorm(13) 
test$z2 <- test$z <- test$x + test$y 
test$z2[1] <- test$z2[1] + 0.5 
test 
##    time status sex           x          y          z          z2 
## 1     4      1   0 -0.56047565  0.1106827 -0.4497929  0.05020707 
## 2     3      1   0 -0.23017749 -0.5558411 -0.7860186 -0.78601862 
## 3     1      1   0  1.55870831  1.7869131  3.3456215  3.34562145 
## 4     1      0   0  0.07050839  0.4978505  0.5683589  0.56835887 
## 5     2      1   1  0.12928774 -1.9666172 -1.8373294 -1.83732942 
## 6     2      1   1  1.71506499  0.7013559  2.4164209  2.41642089 
## 7     3      0   1  0.46091621 -0.4727914 -0.0118752 -0.01187520 
## 8     1      0   0 -1.26506123 -1.0678237 -2.3328849 -2.33288494 
## 9     2      1   0 -0.68685285 -0.2179749 -0.9048278 -0.90482777 
## 10    3      1   0 -0.44566197 -1.0260044 -1.4716664 -1.47166642 
## 11    2      0   1  1.22408180 -0.7288912  0.4951906  0.49519057 
## 12    1      1   1  0.35981383 -0.6250393 -0.2652254 -0.26522544 
## 13    1      0   1  0.40077145 -1.6866933 -1.2859219 -1.28592186 

 test$ztest$xteest$y の和であるから一次従属であり,この 3 変数を使うと間違いなく多重共線性が発生する。
 test$z2test$z の第 1 要素に 0.5 を加えたものであり,一次従属ではない。つまり,一次従属にならないように若干の補正をしたものである。この補正の量によっては多重共線性が発生したり,他の問題が生じたり,計算上は何も問題が起こらないが結果がおかしい,あるいは,何の問題もおかしなこともないというような様々な状態が生じる。


2 相関係数

 変数間の相関係数行列は以下のようになる。

round(cor(test), 3) 
##          time status    sex      x      y      z     z2 
## time    1.000  0.329 -0.161 -0.231 -0.061 -0.165 -0.113 
## status  0.329  1.000 -0.220  0.030  0.235  0.162  0.182 
## sex    -0.161 -0.220  1.000  0.548 -0.375  0.067  0.044 
## x      -0.231  0.030  0.548  1.000  0.454  0.831  0.809 
## y      -0.061  0.235 -0.375  0.454  1.000  0.873  0.886 
## z      -0.165  0.162  0.067  0.831  0.873  1.000  0.996 
## z2     -0.113  0.182  0.044  0.809  0.886  0.996  1.000 

 test$xteest$y は共に test$z および test$z2 と高い相関を持つ。test$ztest$z2 の相関は異常に高いが,分析においては同時に使用しないので何の問題もない。


3 分析例


3.1 一次従属な独立変数を含むデータの場合

 説明変数に,x, yz を使用して分析すると結果は以下のようになる。

summary(coxph(Surv(time, status) ~ x+y+z, test)) 
## Warning in coxph(Surv(time, status) ~ x + y + z, test): X matrix deemed to be singular; variable 3 
## Call: 
## coxph(formula = Surv(time, status) ~ x + y + z, data = test) 
##  
##   n= 13, number of events= 8  
##  
##     coef exp(coef) se(coef)     z Pr(>|z|) 
## x 0.3787    1.4603   0.5293 0.715    0.474 
## y 0.2468    1.2799   0.5533 0.446    0.656 
## z     NA        NA   0.0000    NA       NA 
##  
##   exp(coef) exp(-coef) lower .95 upper .95 
## x      1.46     0.6848    0.5175     4.121 
## y      1.28     0.7813    0.4328     3.786 
## z        NA         NA        NA        NA 
##  
## Concordance= 0.585  (se = 0.166 ) 
## Rsquare= 0.082   (max possible= 0.872 ) 
## Likelihood ratio test= 1.11  on 2 df,   p=0.5737 
## Wald test            = 1.2  on 2 df,   p=0.549 
## Score (logrank) test = 1.25  on 2 df,   p=0.5344 

 x, yz が共に独立変数に含められると,多重共線性の発生を R が察知し,Warning が発生する。
 Warning: X matrix deemed to be singular; variable 3
 この意味は,要するに 3 番目の変数が原因で計算された行列が特異行列である(逆行列を求めることができない—よって,回帰係数の推定値およびそれに引き続く計算もできない)ということである。
 R は問題となった変数(今の場合は z)を分析から外す。分析結果においては,z の推定値欄が全て NA になっているのはそのためである(z を除いた分析として,表示された結果は正しい)。


3.2 多重共線性のあるデータの場合

 説明変数に,x, yz2 を使用して分析すると結果は以下のようになる。多重共線性は存在するが回帰係数の推定値の計算は可能であるということで,次のような Warning が発生する。
 Warning: Loglik converged before variable 1,2,3 ; beta may be infinite.

summary(coxph(Surv(time, status) ~ x+y+z2, test)) 
## Warning in fitter(X, Y, strats, offset, init, control, weights = weights, : Loglik converged before variable 1,2,3 ; beta may be infinite. 
## Call: 
## coxph(formula = Surv(time, status) ~ x + y + z2, data = test) 
##  
##   n= 13, number of events= 8  
##  
##          coef  exp(coef)   se(coef)      z Pr(>|z|) 
## x   3.997e+01  2.286e+17  2.758e+04  0.001    0.999 
## y   4.082e+01  5.359e+17  2.758e+04  0.001    0.999 
## z2 -4.009e+01  3.875e-18  2.758e+04 -0.001    0.999 
##  
##    exp(coef) exp(-coef) lower .95 upper .95 
## x  2.286e+17  4.375e-18         0       Inf 
## y  5.359e+17  1.866e-18         0       Inf 
## z2 3.875e-18  2.580e+17         0       Inf 
##  
## Concordance= 0.756  (se = 0.166 ) 
## Rsquare= 0.277   (max possible= 0.872 ) 
## Likelihood ratio test= 4.22  on 3 df,   p=0.2386 
## Wald test            = 1.56  on 3 df,   p=0.6691 
## Score (logrank) test = 3.31  on 3 df,   p=0.3462 

 回帰係数の推定値は計算できたが,その数値自体おかしい(大きすぎる)とか標準誤差も大きく,z は小さく,p 値も 1 にきわめて近い。また,オッズ比の信頼区間も 0 から無限大ということで,全く信頼できない解である。


4 まとめ

 Cox の比例ハザードモデルにおいては,重回帰分析と同じく多重共線性は分析に影響を及ぼす。