ライブラリーとモジュールのインポート
ライブラリーとモジュールのインポート
・import ライブラリ名
・import ライブラリ名 as 省略名
・import ライブラリ名.モジュール名
・from ライブラリ名 import モジュール名
・import ライブラリ名.モジュール名 as 省略名
・from ライブラリ名 import モジュール名 as 省略名
やってみたら以外に簡単だった「時間依存性共変量」を考慮した生存時間解析
医療現場では、入院時の情報や評価データがその後の転帰(アウトカム)とどのように関連しているかを調べることは非常に多い。そのような場合、入院時の情報、観察期間とアウトカムを用いて『生存時間解析(いわゆる「Cox回帰」』を行うことが一般的であるが、入院時の情報や評価データは、自然若しくは治療介入等により時々刻々と変化する。即ち、入院時は、高血圧、頻脈、発熱ありであっても、その後(極端に言えば翌日には)、正常血圧、正常洞調律、平熱へと変化することは十分に考えらえる。そのような場合、共変量(≒交絡因子)を時間依存性として捉えることで処理が可能である。
データは、survivalパッケージ内にある「原発性胆汁性肝硬変(PBC)」のデータセットを用います。
R: Mayo Clinic Primary Biliary Cholangitis Data
＃ライブラリーのインストール library(survival)
＃データの読込み data(pbc, package="survival") head(pbc)
id time status trt age sex ascites hepato spiders edema bili chol albumin copper alk.phos ast trig platelet protime stage 1 400 2 1 58.76523 f 1 1 1 1.0 14.5 261 2.60 156 1718.0 137.95 172 190 12.2 4 2 4500 0 1 56.44627 f 0 1 1 0.0 1.1 302 4.14 54 7394.8 113.52 88 221 10.6 3 3 1012 2 1 70.07255 m 0 0 0 0.5 1.4 176 3.48 210 516.0 96.10 55 151 12.0 4 4 1925 2 1 54.74059 f 0 1 1 0.5 1.8 244 2.54 64 6121.8 60.63 92 183 10.3 4 5 1504 1 2 38.10541 f 0 1 1 0.0 3.4 279 3.53 143 671.0 113.15 72 136 10.9 3 6 2503 2 2 66.25873 f 0 1 0 0.0 0.8 248 3.98 50 944.0 93.00 63 NA 11.0 3
各変数の定義は以下の通りです。
- id: case number
- time: number of days between registration and the earlier of death, transplantion, or study analysis in July, 1986
- status: status at endpoint, 0(censored) / 1(transplant) / 2(dead)
- trt: treatment not randomised, 1(D-penicillmain) / 2(placebo) /NA(missing)
- age: in years 年齢
- sex: m/f 性別
- ascites: presence of ascites 腹水
- hepato: presence of hepatomegaly or enlarged liver 肝腫大
- spiders: blood vessel malformations in the skin 血管奇形
- edema: 0 no edema浮腫無 / 0.5 untreated or successfully treated 未治療or治療成功の浮腫 / 1 edema despite diuretic therapy 利尿治療を行ったにも関わらず浮腫有
- bili: serum bilirunbin (mg/dl) 血清ビリルビン
- chol: serum cholesterol (mg/dl) 血清コレステロール
- albumin: serum albumin (g/dl) 血清アルブミン
- copper: urine copper (ug/day) 尿中銅
- alk.phos: alkaline phosphotase (U/liter) ALP
- ast: aspartate aminotransferase(U/ml) AST
- trig: triglycerides(mg/dl) 中性脂肪
- platelet: platelet count 血小板数
- protime: standardised blood clotting time 凝固時間
- stage: histologic stage of disease (needs biopsy) 病期
#id<=312とし、項目もid~sexとstageに限定する temp <- subset(pbc, id <=312, select=c(id:sex, stage)) head(temp)
以下は、tempの一部
id time status trt age sex stage 1 1 400 2 1 58.76523 f 4 2 2 4500 0 1 56.44627 f 3 3 3 1012 2 1 70.07255 m 4 4 4 1925 2 1 54.74059 f 4 5 5 1504 1 2 38.10541 f 3 6 6 2503 2 2 66.25873 f 3
ここで、tmerge関数を用いる。
rdrr.io
tmerge関数の構造は以下の通りである。
tmerge(data1, data2, id,..., tstart, tstop, options) - data1 the primary data set, to which new variables and/or observation will be added - data2 second data set in which all the other arguments will be found - id subject identifier - ... operations that add new variables or intervals, see below - tstart optional variable to define the valid time range for each subject, only used on an initial call - tstop optional variable to define the valid time range for each subject, only used on an initial call - options a list of options. Valid ones are idname, tstartname, tstopname, delay, na.rm, and tdcstart. See the explanation below.
tempを用いて下記のように入力する。
pbc2 <- tmerge(temp, temp, id=id, death = event(time, status)) head(pbc2)
すると以下のように、tempにtstart, tstop, death が追加されている。
id time status trt age sex stage tstart tstop death 1 1 400 2 1 58.76523 f 4 0 400 2 2 2 4500 0 1 56.44627 f 3 0 4500 0 3 3 1012 2 1 70.07255 m 4 0 1012 2 4 4 1925 2 1 54.74059 f 4 0 1925 2 5 5 1504 1 2 38.10541 f 3 0 1504 1 6 6 2503 2 2 66.25873 f 3 0 2503 2 > 次に、pbcseqデータを用いて以下のように入力すると、先ほどのpbc2時系列のデータセットが完成する。 pbc3 <- tmerge(pbc2, pbcseq, id=id, ascites = tdc(day, ascites), bili = tdc(day, bili), albumin = tdc(day, albumin), protime = tdc(day, protime), alk.phos = tdc(day, alk.phos), platelet =tdc(day, platelet)) head(pbc3,n=30)
出力結果は以下の通り。
id time status trt age sex stage tstart tstop death ascites bili albumin protime alk.phos platelet 1 1 400 2 1 58.76523 f 4 0 192 0 1 14.5 2.60 12.2 1718 190 2 1 400 2 1 58.76523 f 4 192 400 2 1 21.3 2.94 11.2 1612 183 3 2 4500 0 1 56.44627 f 3 0 182 0 0 1.1 4.14 10.6 7395 221 4 2 4500 0 1 56.44627 f 3 182 365 0 0 0.8 3.60 11.0 2107 188 5 2 4500 0 1 56.44627 f 3 365 768 0 0 1.0 3.55 11.6 1711 161 6 2 4500 0 1 56.44627 f 3 768 1790 0 0 1.9 3.92 10.6 1365 122 7 2 4500 0 1 56.44627 f 3 1790 2151 0 1 2.6 3.32 11.3 1110 135 8 2 4500 0 1 56.44627 f 3 2151 2515 0 1 3.6 2.92 11.5 996 100 9 2 4500 0 1 56.44627 f 3 2515 2882 0 1 4.2 2.73 11.5 860 103 10 2 4500 0 1 56.44627 f 3 2882 3226 0 1 3.6 2.80 11.5 779 113 11 2 4500 0 1 56.44627 f 3 3226 4500 0 1 4.6 2.67 11.5 669 100 12 3 1012 2 1 70.07255 m 4 0 176 0 0 1.4 3.48 12.0 516 151 13 3 1012 2 1 70.07255 m 4 176 364 0 0 1.1 3.29 12.0 353 160 14 3 1012 2 1 70.07255 m 4 364 743 0 0 1.5 3.57 12.0 218 107 15 3 1012 2 1 70.07255 m 4 743 1012 2 0 1.8 3.25 13.3 447 109 16 4 1925 2 1 54.74059 f 4 0 188 0 0 1.8 2.54 10.3 6122 183 17 4 1925 2 1 54.74059 f 4 188 372 0 0 1.6 2.88 19.0 1175 240 18 4 1925 2 1 54.74059 f 4 372 729 0 0 1.7 2.80 11.6 1157 251 19 4 1925 2 1 54.74059 f 4 729 1254 0 0 3.2 2.92 10.8 1178 220 20 4 1925 2 1 54.74059 f 4 1254 1462 0 0 3.7 2.59 13.7 1067 338 21 4 1925 2 1 54.74059 f 4 1462 1824 0 0 4.0 2.59 12.8 1035 200 22 4 1925 2 1 54.74059 f 4 1824 1925 2 1 5.3 1.83 17.0 623 101 23 5 1504 1 2 38.10541 f 3 0 199 0 0 3.4 3.53 10.9 671 136 24 5 1504 1 2 38.10541 f 3 199 391 0 0 1.9 3.28 10.7 689 114 25 5 1504 1 2 38.10541 f 3 391 769 0 0 2.5 3.34 10.5 652 99 26 5 1504 1 2 38.10541 f 3 769 1098 0 0 5.7 3.09 11.4 554 90 27 5 1504 1 2 38.10541 f 3 1098 1455 0 0 5.2 3.02 11.3 588 82 28 5 1504 1 2 38.10541 f 3 1455 1504 1 0 19.0 2.09 13.9 377 101 29 6 2503 2 2 66.25873 f 3 0 378 0 0 0.8 3.98 11.0 944 NA 30 6 2503 2 2 66.25873 f 3 378 737 0 0 0.8 3.91 11.2 605 297
いよいよ生存時間解析(Cox回帰)を実施する。
fit1 <- coxph(Surv(time, status==2) ~ log(bili)+log(protime)+albumin+age+sex+platelet, pbc) summary(fit1)
結果は以下の通り。
Call:
coxph(formula=Surv(time, status==2) ~ log(bili) + log(protime) + albumin + age + sex + platelet, data = pbc)
n= 405, number of events= 154
(13 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
log(bili) 0.8741473 2.3968306 0.0853515 10.242 < 2e-16 ***
log(protime) 4.0245765 55.9566052 0.8760139 4.594 4.34e-06 ***
albumin -0.9667484 0.3803177 0.2034293 -4.752 2.01e-06 ***
age 0.0316387 1.0321446 0.0081862 3.865 0.000111 ***
sexf -0.1389077 0.8703083 0.2368859 -0.586 0.557613
platelet -0.0012133 0.9987875 0.0008059 -1.506 0.132181
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
log(bili) 2.3968 0.41722 2.0276 2.8333
log(protime) 55.9566 0.01787 10.0505 311.5414
albumin 0.3803 2.62938 0.2553 0.5666
age 1.0321 0.96886 1.0157 1.0488
sexf 0.8703 1.14902 0.5471 1.3846
platelet 0.9988 1.00121 0.9972 1.0004
Concordance= 0.829 (se = 0.018 )
Likelihood ratio test= 221.6 on 6 df, p=<2e-16
Wald test = 213.6 on 6 df, p=<2e-16
Score (logrank) test = 262.9 on 6 df, p=<2e-16時間依存性共変量を用いた生存時間解析(Cox回帰)は下記の通りで、time が tstart, tstop となっている。
fit2 <- coxph(tstart, tstop, death==2) ~ log(bili)+log(protime)+albumin+age+sex+platelet, pbc3) summary(fit2)
Call:
coxph(formula=Surv(tstart, tstop, death==2) ~ log(bili)+log(protime)+albumin+age+sex+platelet, data=pbc3)
n= 1803, number of events= 125
(4 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
log(bili) 1.243002 3.466003 0.127946 9.715 < 2e-16 ***
log(protime) 2.757451 15.759617 0.648227 4.254 2.10e-05 ***
albumin -1.592159 0.203486 0.212416 -7.495 6.61e-14 ***
age 0.040416 1.041243 0.010265 3.937 8.25e-05 ***
sexf -0.054067 0.947368 0.258892 -0.209 0.8346
platelet -0.002531 0.997472 0.001109 -2.283 0.0224 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
log(bili) 3.4660 0.28852 2.6973 4.4539
log(protime) 15.7596 0.06345 4.4236 56.1458
albumin 0.2035 4.91435 0.1342 0.3086
age 1.0412 0.96039 1.0205 1.0624
sexf 0.9474 1.05556 0.5704 1.5736
platelet 0.9975 1.00253 0.9953 0.9996
Concordance= 0.915 (se = 0.014 )
Likelihood ratio test= 441 on 6 df, p=<2e-16
Wald test = 281.4 on 6 df, p=<2e-16
Score (logrank) test = 554.6 on 6 df, p=<2e-16推定値は、log(bili)のexp(coef)は 2.39→3.47に、log(protime)のexp(coef)は 55.96→15.76、 albuminは0.38→0.20 と変化しています。時間共変量を考慮した方が詳細な推測ができそうです。
■
分かりやすいです！
jojoshin.hatenablog.com
Pythonでの「クロス表」の書き方
Aが縦列、Bが横列(実行結果はDataFrame形式) #合計値を出力しない pd.crosstab(X['A'], X['B']) #合計値を出力 pd.crosstab(X['A'], X['B'], margins=True)
object型の基本統計量
print(df.describe(include=["O"])) #Oは「オー」
ビニング(binning)
連続変数を、ビン(bin)と呼ばれる離散値に置き換える作業。
