Chapter 6 傾向スコア
傾向スコアは、コホート研究（横断研究）などで良く用いられる手法です。 コホート研究においても、暴露群と対照群を比較することがあります。 ランダム化比較試験と比べると、２群が必ずしも同じような属性を持っているとは限りません。 そのため、研究対象者の中から似たような人々をマッチングさせて、指定した属性が同じようになるようにサブセットを作ります。
傾向スコアを求めるには、ロジステイック回帰を用います。
ln OR\(_y\) = \(\alpha\) + ln OR\(_1\) \(\times\) \(x_1\) + ln OR\(_2\) \(x_2\) + \(\cdots\)
対数の底は自然対数 \(e\)。
OR \(_i\)：
- \(x_i\) が２値変数の時(ex. 男=0 or 女=1)、\(y\) の男に対する女のオッズ比。
- \(x_i\) が連続変数の時、\(y\) の \(x_i\) が1大きくなる時のオッズ比。
OR\(_y\) は、アウトカムの発生オッズ比
ln OR\(_y\) を傾向スコアとします。 このため、傾向スコアは０から１の間の数値をとります。
6.1 傾向スコアの進め方
- 通常通りに２群を比較
- 指定した説明変数で傾向スコアを計算
- 傾向スコアが近いもの同士をマッチング
- マッチングで作成された２群の2で使わなかった説明変数を比較
また、ここでできた２群を、通常通り統計解析することも可能です。
6.2 論文 Masaki S and Kawamoto T (2019)
傾向スコアを使用した論文に、以下のものがあります。
Masaki S & Kawamoto T (2019) Comparison of long-term outcomes between enteral nutrition via gastrostomy and total parenteral nutrition in older persons with dysphagia: A propensity-matched cohort study. PLoS One, 14(10), e0217120.
嚥下障害のある高齢者における人工栄養では、経皮的内視鏡下胃瘻造設術（PEG）が代表的です。しかしながら、日本では近年PEGが不要な延命と見なされていることもあり、（TPN）が不適切に選択されることもあります。そこで、PEGまたはTPNを受けた高齢嚥下障害患者の生存期間などを調査しました。その際に、アウトカム以外の評価項目を用いて1対1の傾向スコアマッチングをいました。
PEG（n＝180）またはTPN（n＝73）を受けた253名の患者を確認し、傾向スコアマッチングにより55組が作成されました。生存期間はPEG群のほうが有意に長い結果となりました（中央値、317日対195日、P = 0.017）。PEG群の方が重症肺炎 (pneumonia) の発生率が高く、逆にTPN群の方が敗血症 (sepsis) の発生率が低かった結果となりました。
この研究でもRが使われているようですが、今回は論文で使われているのとは違うパッケージを使用します。論文のサイトからデータがダウンロードできるようになっています。では、データをダウンロード、解凍した後、以下のコードでデータを読み込んでみましょう。
6.3 xlsx データの読み込み
PLoS ONE のサイトから、以下の記述を見つけてください。 ここから、データをダウンロードすることができます。
Data Availability: All relevant data are available from the Dryad Repository (DOI:10.5061/dryad.gg407h1).
ダウンロードし、zip を解凍してからファイルを読み込みます。
library(readxl)
dfPropensityScore <- read_excel("Dataset_PEG_TPN_Dryad.xlsx")ファイルパスは、ここでは書いていません。 もしうまく読み込めなかったら、フルパスを書いてみてください。
なお、Windows の場合、パスのバックスラッシュまたは円マークをスラッシュに変えてください。
Windows でのパスの例: C:C:¥Documents¥UserName
R でのパスの例: C:/Documents/UserName
dfPropensityScore$PEG <- factor(dfPropensityScore$PEG,
levels = c(1,0),
labels = c("PEG","TPN"))
dfPropensityScore$PORT <- factor(dfPropensityScore$PORT)
dfPropensityScore$NTCVC <- factor(dfPropensityScore$`NT-CVC`)
dfPropensityScore$PICC <- factor(dfPropensityScore$PICC)
dfPropensityScore$sex <- factor(dfPropensityScore$sex)
dfPropensityScore$CI <- factor(dfPropensityScore$CI)
dfPropensityScore$dement <- factor(dfPropensityScore$dement)
dfPropensityScore$NMD <- factor(dfPropensityScore$NMD)
dfPropensityScore$asp <- factor(dfPropensityScore$asp)
dfPropensityScore$IHD <- factor(dfPropensityScore$IHD)
dfPropensityScore$lung <- factor(dfPropensityScore$lung)
dfPropensityScore$liver <- factor(dfPropensityScore$liver)
dfPropensityScore$CKD <- factor(dfPropensityScore$CKD)PEG 列は、数字で 0 と 1 が割り振られていましたが、今後のために TPN, PEG としました。
NT-CVC 列は、マイナス記号があるためにそのままではエラーがでるので、バッククォート記号で囲っています。
なぜなら、- は、引き算を行うと解釈されるからです。
なお、今後これでは使いづらいので、NTCVC という列名に変更しました。
R の解説書では、tidyverse を使って型変換を行うことが多いようです。 tidyverse(dplyr) を使うと、このようになります。
library(dplyr)
dfPropensityScore <- dfPropensityScore %>%
mutate(PEG = factor(PEG,
levels = c(1,0),
labels = c("PEG","TPN")
))6.4 表１の作成
では、Table 1 を再現してみましょう。
パッケージ tableone を読み込み、関数 CreateTableOne を使います。
library(tableone)
CreateTableOne(data = dfPropensityScore,
strata="PEG",
vars=c("PORT","NTCVC","PICC", "age", "sex","CI","dement","NMD","asp","IHD","CHF","lung","liver","CKD","Alb", "TLC", "TC", "Hb", "CRP"),
factorVars=c("PORT","NTCVC","PICC", "sex","CI","dement","NMD","asp","IHD","CHF","lung","liver","CKD"))## Stratified by PEG
## PEG TPN p test
## n 180 73
## PORT = 1 (%) 0 ( 0.0) 28 (38.4) <0.001
## NTCVC = 1 (%) 0 ( 0.0) 26 (35.6) <0.001
## PICC = 1 (%) 0 ( 0.0) 19 (26.0) <0.001
## age (mean (SD)) 81.56 (9.84) 86.77 (6.72) <0.001
## sex = 1 (%) 109 (60.6) 45 (61.6) 0.985
## CI = 1 (%) 107 (59.4) 26 (35.6) 0.001
## dement = 1 (%) 57 (31.7) 45 (61.6) <0.001
## NMD = 1 (%) 10 ( 5.6) 4 ( 5.5) 1.000
## asp = 1 (%) 73 (40.6) 21 (28.8) 0.106
## IHD = 1 (%) 31 (17.2) 16 (21.9) 0.489
## CHF = 1 (%) 70 (38.9) 37 (50.7) 0.114
## lung = 1 (%) 12 ( 6.7) 7 ( 9.6) 0.592
## liver = 1 (%) 9 ( 5.0) 6 ( 8.2) 0.491
## CKD = 1 (%) 29 (16.1) 24 (32.9) 0.005
## Alb (mean (SD)) 3.24 (0.57) 2.82 (0.53) <0.001
## TLC (mean (SD)) 1383.84 (713.18) 1203.11 (682.95) 0.073
## TC (mean (SD)) 160.19 (38.06) 145.58 (43.21) 0.009
## Hb (mean (SD)) 11.29 (1.90) 10.21 (2.10) <0.001
## CRP (mean (SD)) 2.03 (2.77) 2.92 (3.04) 0.026引数
- strata: PEG と TPN の２層 (starta) に分けます。
- vars: Table 1 に掲載する列をすべて与えます。
- factorVars: 上で選択した列のうち、因子の列を指定します。
6.5 多重代入法による欠測処理
まず、欠測値があるため、論文では多重代入法で代入していました。論文では rms ライブラリを使用したとありましたが、ここでは mice を使用します。
なお、当初はデータフレーム dfPropensityScore を直接操作しようとしましたが、以下のようなエラーが出てきました。
iter imp variable
1 1 TLC eval(predvars, data, env) でエラー: オブジェクト 'NT' がありません そこで、欠測値 (NA) がある列と関連すると思われる列だけを抽出して、代入を行い、データフレーム dfPropensityScore に戻すという処理を行いました。
library(mice)
# micePS <- mice(dfPropensityScore)
# dfPropensityScore <- complete(micePS, 1)
dfTemp <- dfPropensityScore[c("ID", "TLC", "TC", "age", "sex", "Alb", "Hb", "CRP")]
names(dfTemp) <- c("ID2", "TLC2", "TC2", "age2", "sex2", "Alb2", "Hb2", "CRP2")
micePS <- mice(dfTemp)##
## iter imp variable
## 1 1 TLC2 TC2
## 1 2 TLC2 TC2
## 1 3 TLC2 TC2
## 1 4 TLC2 TC2
## 1 5 TLC2 TC2
## 2 1 TLC2 TC2
## 2 2 TLC2 TC2
## 2 3 TLC2 TC2
## 2 4 TLC2 TC2
## 2 5 TLC2 TC2
## 3 1 TLC2 TC2
## 3 2 TLC2 TC2
## 3 3 TLC2 TC2
## 3 4 TLC2 TC2
## 3 5 TLC2 TC2
## 4 1 TLC2 TC2
## 4 2 TLC2 TC2
## 4 3 TLC2 TC2
## 4 4 TLC2 TC2
## 4 5 TLC2 TC2
## 5 1 TLC2 TC2
## 5 2 TLC2 TC2
## 5 3 TLC2 TC2
## 5 4 TLC2 TC2
## 5 5 TLC2 TC2dfTemp <- complete(micePS, 1)
dfPropensityScore <- cbind(dfPropensityScore, dfTemp)次に、傾向スコアマッチングを実行します。
library("MatchIt")
library("optmatch")##
## 次のパッケージを付け加えます: 'optmatch'## 以下のオブジェクトは 'package:survival' からマスクされています:
##
## stratamiPropensityScore <- matchit(
formula = PEG ~ age + sex + CI + dement + NMD + asp + IHD + CHF + lung + liver + CKD + Alb + TLC2 + TC2 + Hb + CRP,
data = dfPropensityScore,
method = "nearest",
distance = "glm",
caliper = 0.05,
ratio = 1)
summary(miPropensityScore)##
## Call:
## matchit(formula = PEG ~ age + sex + CI + dement + NMD + asp +
## IHD + CHF + lung + liver + CKD + Alb + TLC2 + TC2 + Hb +
## CRP, data = dfPropensityScore, method = "nearest", distance = "glm",
## caliper = 0.05, ratio = 1)
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.4620 0.2182 1.2167 1.0721 0.3174
## age 86.7671 81.5556 0.7753 0.4663 0.1121
## sex0 0.3836 0.3944 -0.0224 . 0.0109
## sex1 0.6164 0.6056 0.0224 . 0.0109
## CI0 0.6438 0.4056 0.4976 . 0.2383
## CI1 0.3562 0.5944 -0.4976 . 0.2383
## dement0 0.3836 0.6833 -0.6165 . 0.2998
## dement1 0.6164 0.3167 0.6165 . 0.2998
## NMD0 0.9452 0.9444 0.0033 . 0.0008
## NMD1 0.0548 0.0556 -0.0033 . 0.0008
## asp0 0.7123 0.5944 0.2604 . 0.1179
## asp1 0.2877 0.4056 -0.2604 . 0.1179
## IHD0 0.7808 0.8278 -0.1135 . 0.0470
## IHD1 0.2192 0.1722 0.1135 . 0.0470
## CHF 0.5068 0.3889 0.2359 . 0.1180
## lung0 0.9041 0.9333 -0.0993 . 0.0292
## lung1 0.0959 0.0667 0.0993 . 0.0292
## liver0 0.9178 0.9500 -0.1172 . 0.0322
## liver1 0.0822 0.0500 0.1172 . 0.0322
## CKD0 0.6712 0.8389 -0.3569 . 0.1677
## CKD1 0.3288 0.1611 0.3569 . 0.1677
## Alb 2.8233 3.2417 -0.7854 0.8686 0.1381
## TLC2 1190.4507 1383.8361 -0.2883 0.8844 0.1035
## TC2 145.1233 159.1333 -0.3251 1.3139 0.0969
## Hb 10.2096 11.2922 -0.5161 1.2193 0.1239
## CRP 2.9194 2.0327 0.2920 1.2037 0.1416
## eCDF Max
## distance 0.5130
## age 0.2553
## sex0 0.0109
## sex1 0.0109
## CI0 0.2383
## CI1 0.2383
## dement0 0.2998
## dement1 0.2998
## NMD0 0.0008
## NMD1 0.0008
## asp0 0.1179
## asp1 0.1179
## IHD0 0.0470
## IHD1 0.0470
## CHF 0.1180
## lung0 0.0292
## lung1 0.0292
## liver0 0.0322
## liver1 0.0322
## CKD0 0.1677
## CKD1 0.1677
## Alb 0.3241
## TLC2 0.2117
## TC2 0.2016
## Hb 0.2868
## CRP 0.3149
##
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.3970 0.3950 0.0101 1.0158 0.0057
## age 86.7059 85.3333 0.2042 1.1671 0.0365
## sex0 0.3529 0.4902 -0.2823 . 0.1373
## sex1 0.6471 0.5098 0.2823 . 0.1373
## CI0 0.6667 0.6667 0.0000 . 0.0000
## CI1 0.3333 0.3333 0.0000 . 0.0000
## dement0 0.4314 0.4510 -0.0403 . 0.0196
## dement1 0.5686 0.5490 0.0403 . 0.0196
## NMD0 0.9608 0.9216 0.1723 . 0.0392
## NMD1 0.0392 0.0784 -0.1723 . 0.0392
## asp0 0.6471 0.6078 0.0866 . 0.0392
## asp1 0.3529 0.3922 -0.0866 . 0.0392
## IHD0 0.7843 0.7843 0.0000 . 0.0000
## IHD1 0.2157 0.2157 0.0000 . 0.0000
## CHF 0.4706 0.5686 -0.1961 . 0.0980
## lung0 0.9216 0.9216 0.0000 . 0.0000
## lung1 0.0784 0.0784 0.0000 . 0.0000
## liver0 0.9216 0.9216 0.0000 . 0.0000
## liver1 0.0784 0.0784 0.0000 . 0.0000
## CKD0 0.7451 0.6863 0.1252 . 0.0588
## CKD1 0.2549 0.3137 -0.1252 . 0.0588
## Alb 2.9216 2.9471 -0.0479 0.6072 0.0327
## TLC2 1193.5784 1082.3608 0.1658 1.3236 0.0826
## TC2 138.5882 145.0588 -0.1502 1.4710 0.0604
## Hb 10.5020 10.1353 0.1748 1.0127 0.0487
## CRP 3.0137 2.5931 0.1385 1.0619 0.0737
## eCDF Max Std. Pair Dist.
## distance 0.0588 0.0194
## age 0.1373 0.9801
## sex0 0.1373 0.8468
## sex1 0.1373 0.8468
## CI0 0.0000 0.2353
## CI1 0.0000 0.2353
## dement0 0.0196 0.7662
## dement1 0.0196 0.7662
## NMD0 0.0392 0.5170
## NMD1 0.0392 0.5170
## asp0 0.0392 0.8663
## asp1 0.0392 0.8663
## IHD0 0.0000 0.3529
## IHD1 0.0000 0.3529
## CHF 0.0980 0.8236
## lung0 0.0000 0.1569
## lung1 0.0000 0.1569
## liver0 0.0000 0.1176
## liver1 0.0000 0.1176
## CKD0 0.0588 0.9600
## CKD1 0.0588 0.9600
## Alb 0.1373 1.0343
## TLC2 0.1961 0.9573
## TC2 0.1765 0.9729
## Hb 0.1765 0.9637
## CRP 0.1765 1.1525
##
## Sample Sizes:
## Control Treated
## All 180 73
## Matched 51 51
## Unmatched 129 22
## Discarded 0 0関数
- matchit: ライブラリ MatchiIt の関数
引数
- formula: ロジスティック回帰のモデル。なお、‘=’ は引数をとる時に使うため、モデル内の ‘=’ は ‘~’ と書きます。
- data: データフレーム
- method: マッチングの方法。より厳密な “exact” や、サンプルサイズが大きくなる “full” （ただし、1:1にはなりません）などがあります。
- distance
- caliper: distance の数値でマッチングを行う際に許容される差の閾値。
- ratio: 両群のサンプルサイズの比率。
dfPSMatched <- match.data(miPropensityScore)CreateTableOne(data = dfPSMatched,
strata="PEG",
vars=c("PORT","NTCVC","PICC", "age", "sex","CI","dement","NMD","asp","IHD","CHF","lung","liver","CKD","Alb", "TLC", "TC", "Hb", "CRP"),
factorVars=c("PORT","NTCVC","PICC", "sex","CI","dement","NMD","asp","IHD","CHF","lung","liver","CKD"))## Stratified by PEG
## PEG TPN p test
## n 51 51
## PORT = 1 (%) 0 ( 0.0) 19 (37.3) <0.001
## NTCVC = 1 (%) 0 ( 0.0) 19 (37.3) <0.001
## PICC = 1 (%) 0 ( 0.0) 13 (25.5) <0.001
## age (mean (SD)) 85.33 (6.59) 86.71 (7.12) 0.315
## sex = 1 (%) 26 (51.0) 33 (64.7) 0.229
## CI = 1 (%) 17 (33.3) 17 (33.3) 1.000
## dement = 1 (%) 28 (54.9) 29 (56.9) 1.000
## NMD = 1 (%) 4 ( 7.8) 2 ( 3.9) 0.674
## asp = 1 (%) 20 (39.2) 18 (35.3) 0.838
## IHD = 1 (%) 11 (21.6) 11 (21.6) 1.000
## CHF = 1 (%) 29 (56.9) 24 (47.1) 0.428
## lung = 1 (%) 4 ( 7.8) 4 ( 7.8) 1.000
## liver = 1 (%) 4 ( 7.8) 4 ( 7.8) 1.000
## CKD = 1 (%) 16 (31.4) 13 (25.5) 0.661
## Alb (mean (SD)) 2.95 (0.59) 2.92 (0.46) 0.809
## TLC (mean (SD)) 1082.36 (520.14) 1205.66 (605.41) 0.281
## TC (mean (SD)) 145.11 (35.19) 139.12 (42.42) 0.453
## Hb (mean (SD)) 10.14 (1.87) 10.50 (1.88) 0.326
## CRP (mean (SD)) 2.59 (3.17) 3.01 (3.27) 0.511library(survival)
library(survminer)
objSurv <- survfit( Surv(survival,status) ~ PEG,
data = dfPSMatched)
ggsurvplot(objSurv,
data = dfPSMatched,
risk.table = TRUE,
pval = TRUE,
conf.int = TRUE,
risk.table.y.text.col = TRUE)
dfPSMatched$oral <- as.factor(dfPSMatched$oral)
dfPSMatched$home <- as.factor(dfPSMatched$home)
dfPSMatched$pneumo <- as.factor(dfPSMatched$pneumo)
dfPSMatched$sepsis <- as.factor(dfPSMatched$sepsis)
CreateTableOne(data = dfPSMatched,
strata="PEG",
vars=c("oral", "home", "pneumo", "sepsis"),
factorVars=c("oral", "home", "pneumo", "sepsis"))## Stratified by PEG
## PEG TPN p test
## n 51 51
## oral = 1 (%) 4 ( 7.8) 2 ( 3.9) 0.674
## home = 1 (%) 11 (21.6) 5 ( 9.8) 0.173
## pneumo = 1 (%) 28 (54.9) 12 (23.5) 0.002
## sepsis = 1 (%) 5 ( 9.8) 15 (30.0) 0.022