LiSA

図2：propensity score のずれた分布

# 適当なPS分布の作成
s1 <- seq(0, 10, length=100)  # 0 から 10 まで長さ100 の規則的な数列を作成する
ino1 <- dnorm(s1, 6, 1.2) # 平均 6，分散 1.2 の正規分布から，s1 に対応する確率密度を計算する
ino2 <- dnorm(s1, 4, 1.4) # 平均 4，分散 1.4 の正規分布から，s1 に対応する確率密度を計算する

xl <- c(-max(ino2), max(ino1)) # x 軸のプロットする範囲を決める
yl <- c(0, 1)                  # y 軸のプロットする範囲を決める

# プロットする
# x軸みぎ側に強心剤を使用した患者の傾向スコア
# x軸ひだり側に強心剤を使用しなかった患者の傾向スコア
# y軸に確率 [0, 1]
# がくるようにプロットしたい
y <- seq(0, 1, length=length(s1)) # 0 から 1 まで長さがs1 と同じな規則的な数列を作成する
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, las=1)
abline(v=0, lty=2)
lines(ino1, y, type="l", col="red", lwd=5)    # 強心剤が投与された患者の傾向スコア
lines(-ino2, y, type="l", col="blue", lwd=5)  # 強心剤が投与されなかった患者の傾向スコア
text(mean(c(0, par()$usr[1])), par()$usr[3], "No Inotrope", xpd=TRUE, adj=c(NA, 4), cex=1.5)
text(mean(c(0, par()$usr[2])), par()$usr[3], "Inotrope", xpd=TRUE, adj=c(NA, 4), cex=1.5)

図3：propensity score が一致した集団の分布

minPS <- apply(cbind(ino1, ino2), 1, min)
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, las=1)
polygon(c(minPS, rev(-minPS)), c(y, rev(y)), border=NA, col=grey(0.85)) # 傾向スコアが似ている領域を塗りつぶす
abline(v=0, lty=2)
y <- seq(0, 1, length=length(s1))
lines(ino1, y, type="l", col="red", lwd=3, lty=3)
lines(-ino2, y, type="l", col="blue", lwd=3, lty=3)
text(mean(c(0, par()$usr[1])), par()$usr[3], "No Inotrope", xpd=TRUE, adj=c(NA, 4), cex=1.5)
text(mean(c(0, par()$usr[2])), par()$usr[3], "Inotrope", xpd=TRUE, adj=c(NA, 4), cex=1.5)
lines(minPS, y, type="l", col="red", lwd=5)
lines(-minPS, y, type="l", col="blue", lwd=5)

対応のある検定を対応のない検定としてしまったときの\(p\) 値の齟齬

一般的な帰無仮説検定で，帰無仮説\(H_0:\mu_1=\mu_2\) と対立仮説\(H_1:\mu_1\not=\mu_2\) として，帰無仮説が偽である（差があることが確実に分かっている状態）のみを原稿では提示していますが，帰無仮説が真（差がないことが確実に分かっている状態）もシミュレーションします．

帰無仮説が偽のときは，サンプルサイズが大きくなるほど，対応のないt検定が小さなp値を返す確率が大きく，結果として有意な結果として報告されやすくなります．

逆に，帰無仮説が真のときは，サンプルサイズが小さくなるほど，対応のないt検定が小さなp値を返す確率が大きく，サンプルサイズが大きくなるほど，対応のないt検定と対応のあるt検定のp値の差は半々になります．ただし，帰無仮説が真であるので，たとえ差がないとしても，偶然に有意な結果になる確率が\(p=0.05\) あります．

rhc データセットを使って傾向スコアマッチングを行います．

library(tableone) # tableone package itself
library(Matching) # PS matching
library(mice)     # 欠損値の補完

# データの読み込み
url <- "http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.csv"
rhc <- read.csv(url)              # ネットにつながっているなら
#rhc <- read.csv("rhc.csv", row=1)  # PC に保存している場合

データ総数は

nrow(rhc)

[1] 5735

RHC の割付は

table(rhc$swang1)

No RHC    RHC 
  3551   2184 

と，論文の数値と一致する．

例えば，年齢別のRHC 分布は以下の通りである．

table(cut(rhc$age, c(0, 50, 60, 70, 80, 120), right=FALSE), rhc$swang1)

          
           No RHC RHC
  [0,50)      884 540
  [50,60)     546 371
  [60,70)     812 577
  [70,80)     809 529
  [80,120)    500 167

例えば，性別によるのRHC 分布は以下の通りである．

table(rhc$sex, rhc$swang1)

        
         No RHC  RHC
  Female   1637  906
  Male     1914 1278

主要評価項目である30日死亡（dth30）は一致する．

table(rhc$dth30, rhc$swang1) # 30日死亡

     
      No RHC  RHC
  No    2463 1354
  Yes   1088  830

180日死亡（death）は一致しない．

table(rhc$death, rhc$swang1) # 180日死亡　合わない

     
      No RHC  RHC
  No    1315  698
  Yes   2236 1486

Propensity score の項に，傾向スコアの計算に使った変数が書かれているが，primary disease category の変数と論文のTable 1 の総数が一致しないので適宜解析しやすいように編集する．

vars <- c(
  "age", "sex", "race", "edu", "income", "ninsclas", "cat1",
  "resp", "card", "neuro", "gastr", "renal", "meta", "hema", "seps", "trauma", "ortho", 
  "adld3p", "das2d3pc", "dnr1", "ca", "surv2md1", "aps1", "scoma1", "wtkilo1",
  "temp1", "meanbp1", "resp1", "hrt1", "pafi1", "paco21", "ph1", "wblc1", "hema1", "sod1", "pot1", "crea1", "bili1", "alb1", "urin1", 
  "cardiohx", "chfhx", "dementhx", "psychhx", "chrpulhx", "renalhx", "liverhx", "gibledhx", "malighx", "immunhx", "transhx", "amihx"
  )

# データの整形
rhc1 <- rhc[, c("swang1", vars)]
cat1 <- as.character(rhc1$cat1)
cat1[cat1 == "MOSF w/Malignancy"] <- "MOSF"
cat1[cat1 == "MOSF w/Sepsis"] <- "MOSF"
cat1[!is.na(match(cat1, c("COPD", "Cirrhosis", "Colon Cancer", "Coma" ,"Lung Cancer")))] <- "Other"
rhc1$cat1 <- factor(cat1)

# 傾向スコアを計算するために回帰モデルを作成する
psModel <- glm(swang1 ~ ., data=rhc1, family=binomial(link="logit"))

# urin1 やdas2d3pc に欠損値があるが論文ではまったく触れられていないので
# 適当に欠損値補完を行うことにする
imp <- complete(mice(rhc1, method="rf", printFlag=FALSE))
pRHC <- predict.glm(psModel, newdata=imp, type = "response") # 欠損値補完したうえでの，各患者の傾向スコア

# マッチングを行う
Tr <- (rhc$swang1 == "RHC")
X <- log(pRHC / (1-pRHC))
listMatch <- Match(Tr       = Tr,      # Need to be in 0,1
                   ## logit of PS,i.e., log(PS/(1-PS)) as matching scale
                   X        = X,
                   ## 1:1 matching
                   M        = 1,
                   ## caliper = 0.2 * SD(logit(PS))
                   caliper  = 0.2,
                   replace  = FALSE,
                   ties     = TRUE,
                   version  = "fast")
idx <- c(listMatch$index.treated, listMatch$index.control) # マッチングした患者ペアのindex

rhc2 <- rhc1
rhc2 <- imp
rhc2$swang1 <- c(rhc2$swang1) - 1
mb <- MatchBalance(swang1 ~ ., data=rhc2, digit=5, match.out=listMatch, nboots=100, print.level=0)

傾向スコアの分布を確認する．

d <- tapply(pRHC, rhc$swang1, density)
Y <- do.call(cbind, sapply(d, "[")["y",])
Y <- sweep(Y, 2, table(rhc$swang1)/nrow(rhc), "*")
xl <- c(-max(Y[,2]), max(Y[,1]))
xl <- c(-1, 1) * max(abs(Y))
yl <- c(0, 1)
y <- seq(0, 1, length=nrow(Y))


d1 <- tapply(pRHC[idx], rhc$swang1[idx], density)
Y1 <- do.call(cbind, sapply(d1, "[")["y",])
Y2 <- sweep(Y1, 2, table(rhc$swang1[idx])/table(rhc$swang1), "*")
Y2 <- sweep(Y1, 2, table(rhc$swang1)/nrow(rhc)*table(rhc$swang1[idx])/table(rhc$swang1), "*")

minPS <- apply(Y, 1, min)
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, las=1)
polygon(c(Y2[,2], rev(-Y2[,1])), c(y, rev(y)), border=NA, col=grey(0.85))
abline(v=0, lty=2)
points(Y[,2], y, type="l", col="red", lwd=1.5, lty=3)
points(-Y[,1], y, type="l", col="blue", lwd=1.5, lty=3)
points(Y2[,2], y, type="l", col="red", lwd=3, lty=1)
points(-Y2[,1], y, type="l", col="blue", lwd=3, lty=1)
text(mean(c(0, par()$usr[1])), par()$usr[3], levels(rhc$swang1)[1], xpd=TRUE, adj=c(NA, 4), cex=1.5)
text(mean(c(0, par()$usr[2])), par()$usr[3], levels(rhc$swang1)[2], xpd=TRUE, adj=c(NA, 4), cex=1.5)

マッチング前後の偏りを検討する．

tableone パッケージを使うと，検定によるp値もSMD も計算できる．

SMD では一部 10% 以下にならない変数も存在する．

# 2群の偏りを調べる
# マッチング前
Tab1 <- CreateTableOne(vars = vars, strata = "swang1", data = imp, test = TRUE)
# マッチング後
Tab2 <- CreateTableOne(vars = vars, strata = "swang1", data = imp[idx,], test = TRUE)

smd <- do.call(cbind, lapply(list(Tab1, Tab2), ExtractSmd))
cols <- c("black", "red")
par(mar=c(5, 8, 2, 2))
plot(0, xlim=range(0, smd), ylim=c(nrow(smd)-1, 2), type="n", yaxt="n",
     xlab="Standardized mean difference", ylab="")
abline(v=0.1, lty=3)  # SMD 10% の線
for(i in seq(nrow(smd))){
  abline(h=i, lty=3)
  points(smd[i,], rep(i, 2), pch=16, col=cols)
  text(par()$usr[1], i, rownames(smd)[i], xpd=TRUE, pos=2)
}

生存解析を行う．

#Sadmdte
#Study Admission Date
#
#Dthdte
#Date of Death
#
#Lstctdte
#Date of Last Contact
#
#Dschdte
#Hospital Discharge Date
#
#Death
#Death at any time up to 180 Days

library(survival)
d0 <- rhc$lstctdte-rhc$sadmdte
d1 <- rhc$dthdte-rhc$sadmdte
interval <- apply(cbind(d0, d1), 1, max, na.rm=TRUE)
event <- (!is.na(d1)) + 0

# 30日死亡のKM
sdat <- cbind.data.frame(days=replace(d0, d0 > 30, 30), event=(d0 < 30) + 0, swang1=rhc$swang1)
sfit <- survfit(Surv(days, event) ~ swang1, data=sdat)

par(mar=c(4.5, 5, 2, 2), cex.lab=1.5, cex.axis=1.5, las=1)
cols <- c("red", "blue")
plot(sfit, ylim=c(0.6, 1), col=cols, lwd=5,
     xlab="Follow-up Time [day]", ylab="Proportion Surviving")
legend("topright", legend=levels(sdat$swang1), col=cols, lwd=5, cex=2)

マッチング後の生存曲線．

# マッチング後
idx <- c(listMatch$index.treated, listMatch$index.control)
mdat <- sdat[idx,]
mfit <- survfit(Surv(days, event) ~ swang1, data=mdat)

par(mar=c(4.5, 5, 2, 2), cex.lab=1.5, cex.axis=1.5, las=1)
cols <- c("red", "blue")
plot(mfit, ylim=c(0.6, 1), col=cols, lwd=5,
     xlab="Follow-up Time [day]", ylab="Proportion Surviving")
legend("topright", legend=levels(sdat$swang1), col=cols, lwd=5, cex=2)

ログランク検定は下記で行える．

# source("http://aoki2.si.gunma-u.ac.jp/R/src/logrank.R", encoding="euc-jp") 
# ログランク検定を行う
logrank <- function(    group,                  # 群を識別するベクトル（1, 2 のいずれか）
            event,                  # 死亡なら 1，生存なら 0 の値をとるベクトル
            time,                   # 生存期間ベクトル
            method=c("SAS", "Tominaga"))        # 一般的な SAS などの方法か，富永の方法か
{
    method <- match.arg(method)
    data.name <- sprintf("time: %s, event: %s, group: %s",
        deparse(substitute(time)), deparse(substitute(event)), deparse(substitute(group)))
    OK <- complete.cases(group, event, time)        # 欠損値を持つケースを除く
    group <- group[OK]
    event <- event[OK]
    time <- time[OK]
    len <- length(group)
    stopifnot(length(event) == len, length(time) == len)
    tg <- table(c(time, rep(NA, 4)),            # 後ろの4項目はダミー
            c(group, 1, 1, 2, 2)*10+c(event, 1, 0, 1, 0))
    k <- nrow(tg)
    nia <- table(group)[1]
    nib <- len-nia
    na <- c(nia, (rep(nia, k)-cumsum(tg[,1]+tg[,2]))[-k])
    nb <- c(nib, (rep(nib, k)-cumsum(tg[,3]+tg[,4]))[-k])
    da <- tg[,2]
    db <- tg[,4]
    dt <- da+db
    nt <- na+nb
    d <- dt/nt
    O <- c(sum(da), sum(db))
    ea <- na*d
    eb <- nb*d
    E <- c(sum(ea), sum(eb))

    result <- data.frame(da, db, dt, na, nb, nt, d, ea, eb)

    if (method == "Tominaga") {             # 富永による検定方式
        method <- "ログランク検定（富永）"
        chi <- sum((O-E)^2/E)
    }
    else {                          # SAS 等と同じ検定方式
        method <- "ログランク検定（一般的）"
        v <- sum(dt*(nt-dt)/(nt-1)*na/nt*(1-na/nt), na.rm=TRUE)
        chi <- (sum(da)-sum(na*d))^2/v
    }
    P <- pchisq(chi, 1, lower.tail=FALSE)
    return(structure(list(statistic=c("X-squared"=chi), parameter=c(df=1), p.value=P,
    method=method, data.name=data.name, result=result), class="htest"))
}