Propensity scoreを勉強したのでやってみる。
JAMA. 1996 Sep 18;276(11):889-97.のRight heart catheterization datasetを使う。
Propensity scoreに必要なRパッケージについてと、スクリプトはそれぞれリンクから適当にパクる。
結論からいうと論文入手できてないし、アブストから察するにPS調整してもRHCは予後不良なのにそうならなかった。やる気があれば修正する(たぶんしない)。
データセットを読み込んで、データとして使えそうなパラメータを取ってきて、PSマッチングして、調整前と調整後の生存曲線を比較してみる。
最初の状態では、RHCをされる確率(PS)がRHC群では高く、No RHC群では低い。この状態で生存曲線を作ると、RHC群のほうが予後が悪いということになる。
これは今までの反論なら、「RHCをする人はもともと重症。だから予後が悪い」とかいちゃもんがつく。
# ネットが繋がっていれば rhc <- read.csv("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.csv") rhc1 <- rhc[, c(2, 4, 9:58, 61:62)] # NA とか変なカテゴリパラメタなどを削除したとして # 生存日数 # たぶんこれが違うから、今回の結果がおかしいのだと思う。 # sadmdte study admission date # dschdte hospital discharge date # dthdte date of death # lstctdte date of last contact # となっていて、sadmdte から換算して死んでいれば dthdte, 生きていれば lstctdte までの日にちを取ったわけだが… day <- numeric(nrow(rhc)) for(i in seq(day)){ day[i] <- c(rhc$lstctdte[i], rhc$dthdte[i])[as.integer(rhc$death[i])] - rhc$sadmdte[i] } # 生存解析 library(survival) event <- as.integer(rhc$death) group <- rhc$swang1 s1 <- survfit(Surv(time=day, event=event) ~ group) plot(s1, col=c("blue", "red"), lwd=2) legend("topright", legend=c("No RHC", "RHC"), col=c("blue", "red"), lwd=3, bty="n", cex=2) # PSの計算 # RHC, すなわち Swan-Ganz カテーテルを行うか否かを各パラメータからロジスティック回帰で定量化する。 swg <- glm(swang1 ~ . - death, family = binomial(logit), data = rhc1) D <- as.integer(rhc1$swang1)-1 # 治療する Y <- as.integer(rhc1$death)%%2 # 結果 X <- swg$fitted.values # PS dens <- mapply(function(x) density(X[D==x]), 0:1, SIMPLIFY=FALSE) xl <- c(-1, 1)*max(dens[[1]]$y, dens[[2]]$y) yl <- c(0, 1) par(mar=c(5, 4.5, 2, 2)) plot(0, type="n", xlim=xl, ylim=yl, xlab="", ylab="Propensity score (PS)", cex.lab=1.6, cex.axis=1.6) abline(v=0, lty=2) points(dens[[2]]$y, dens[[2]]$x, type="l", col="red", lwd=2) points(-dens[[1]]$y, dens[[1]]$x, type="l", col="blue", lwd=2) text(mean(c(0, par()$usr[1])), par()$usr[3], "No RHC", xpd=TRUE, adj=c(NA, 4), cex=1.5, col="blue") text(mean(c(0, par()$usr[2])), par()$usr[3], "RHC", xpd=TRUE, adj=c(NA, 4), cex=1.5, col="red")
Call: survfit(formula = Surv(time = day, event = event) ~ group) records n.max n.start events median 0.95LCL 0.95UCL group=No RHC 3551 3551 3551 2236 233 209 247 group=RHC 2184 2184 2184 1486 106 80 136
PS matching にはMatchingパッケージを使う。
PS matchingするさいのcaliperの設定だとか、プロットするときの確率密度を非調整前と対応づけるのに若干失敗しているが、こんな感じでPSがほぼ同一と考えられる人を対応づけると、もともといた集団(RHC: 2184人、No RHC: 3551人)が1332人ずつの集団に減る。
この集団で生存解析をすると、RHCの予後の悪さっぷりがちょっとだけ改善されたように見えるけれども、ログランク検定ではやはり予後不良。生存中央値で見れば有意差はなさそうになる。
?? で、有意差なくなりそうじゃん…と思いつつ、論文では30日予後や医療費についても解析しているのでもうやめた!!!
# マッチングする # Match - without replacement # method ATT ATC ATE m.obj <- Match(Y = Y, Tr = D, X = X, M = 1, replace=FALSE, estimand = "ATT", caliper=0.002) m <- m.obj$MatchLoopC D1 <- D[m[, 1:2]] Y1 <- Y[m[, 1:2]] X1 <- X[m[, 1:2]] # PS dens1 <- mapply(function(x) density(X1[D1==x]), 0:1, SIMPLIFY=FALSE) xl <- c(-1, 1)*max(dens[[1]]$y, dens[[2]]$y, dens1[[1]]$y, dens1[[2]]$y) yl <- c(0, 1) par(mar=c(5, 4.5, 2, 2)) plot(0, type="n", xlim=xl, ylim=yl, xlab="", ylab="Propensity score (PS)", cex.lab=1.6, cex.axis=1.6) abline(v=0, lty=2) x <- length(X1)/length(X) points(dens[[2]]$y, dens[[2]]$x, type="l", col="red", lwd=1, lty=3) points(-dens[[1]]$y, dens[[1]]$x, type="l", col="blue", lwd=1, lty=3) points(dens1[[2]]$y*x, dens1[[2]]$x, type="l", col="red", lwd=2) points(-dens1[[1]]$y*x, dens1[[1]]$x, type="l", col="blue", lwd=2) text(mean(c(0, par()$usr[1])), par()$usr[3], "No RHC", xpd=TRUE, adj=c(NA, 4), cex=1.5, col="blue") text(mean(c(0, par()$usr[2])), par()$usr[3], "RHC", xpd=TRUE, adj=c(NA, 4), cex=1.5, col="red") legend("topleft", legend=c("Raw", "PS match"), lty=c(3, 1), bty="n", cex=1.5, lwd=2) # PS マッチしたあとの集団での生存解析 day <- day[m[, 1:2]] event <- event[m[, 1:2]] group <- group[m[, 1:2]] s1 <- survfit(Surv(time=day, event=event) ~ group) plot(s1, col=c("blue", "red"), lwd=2) legend("topright", legend=c("No RHC", "RHC"), col=c("blue", "red"), lwd=3, bty="n", cex=2)
Call: survfit(formula = Surv(time = day, event = event) ~ group) records n.max n.start events median 0.95LCL 0.95UCL group=No RHC 1342 1342 1342 840 217 153 250 group=RHC 1342 1342 1342 904 172 128 210
ログランク検定は青木先生のサイトからパクる。
source("http://aoki2.si.gunma-u.ac.jp/R/src/logrank.R", encoding="euc-jp") logrank(group=as.integer(group), event=event, time=day)
