表からC-index　傾向スコア　操作変数

C-index

某論文で、傾向スコアと操作変数を使った解析をやっていた。しかし、どちらもデータベースベースドの研究にありがちなこれは調整できてへんやろ感が激しい。

傾向スコアの前提が信頼できるかどうかの指標として、c-indexは使用するべきではないとされる。これは、たとえば1:1RCTの場合、true propensity scoreが全症例で1/2になり、c-indexは0.5になることからも本質とは関係ないと想像できる。しかし、現実の治療で1:1に何となく治療が選ばれることは、極めて考えにくい。実際に、今回の論文の治療選択においても、相当に重症度など患者背景が影響するであろうことは想像できる。したがって、あまりにも低いc-indexはpropensity modelの構築がうまくいっていないことを示唆する。

ということで、c-indexは明記されていないが、図からどんなもんか計算してみたくなった。

傾向スコアの分布は図示されているので、拡大コピーの上各階級の棒グラフの長さを計測してc-indexを計算する。

##傾向スコアc-index

ct.v <- c(2,1,1,1,1,1,2,6,4,5,4,7,9,11,12,14,15,18,19,22,25,29,35,37,38,47,51,57,64,70,75,89,98,105,117,126,129,144,139,152,143,153,157,164,154,152,151,145,126,114,110,106,89,79,69,66,54,41,40,26,24,18,14,12,7,5,4,3,2,2,1)
tr.v <- c(0,1,1,1,0,1,2,1,1,3,3,3,3,5,6,5,6,10,12,15,14,22,22,26,27,34,37,43,50,60,67,72,85,96,102,116,128,137,144,154,164,170,171,185,180,175,180,173,174,153,149,142,128,116,116,111,85,72,60,53,45,35,26,21,15,12,8,6,3,2,1)
ct.n <- sum(ct.v) 
tr.n <- sum(tr.v) 
cat.n <- length(ct.v)
sn.v <- sp.m.v <- numeric(cat.n)

for (i in  1:(length(ct.v)-1)){
      sn.v[i] <- sum(tr.v[(i+1):cat.n]) / tr.n
      sp.m.v[i] <- sum(ct.v[(i+1):cat.n]) / ct.n
}


tr.yes <- ps <- numeric(0)
for (i in 1:cat.n){
      ps <- c(ps, rep(i, tr.v[i]))
      tr.yes <- c(tr.yes, rep(1, tr.v[i]))
      ps <- c(ps, rep(i, ct.v[i]))
      tr.yes <- c(tr.yes, rep(0, ct.v[i]))
}

library(pROC)

## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## 
##  以下のオブジェクトはマスクされています (from 'package:stats') : 
## 
##      cov, smooth, var

ROC1 <- roc(tr.yes, ps) 
ROC1.CI <- ci.auc(ROC1, conf.level=0.95) # AUC delong
plot(ROC1)

## 
## Call:
## roc.default(response = tr.yes, predictor = ps)
## 
## Data: ps in 4013 controls (tr.yes 0) < 4446 cases (tr.yes 1).
## Area under the curve: 0.575