time_varying_vars = c("Hb","Transfusion")
time_fix_vars = c("SubjectNo", 'AGE',"SEX", 'BMI', 'CFScore', 'Charlson_Comorbidity_index_total',
'Comorbidity_MOSIAC_SMT','CRRT_within_3days', 'input_preICU_to_ICUD1_quartile','SOFA_ICUD1','SAPS3_ICUD1', 'Septic_shock_ICUD1', 's_pulmonary', 's_abdominal', 's_urinary','TypeOfInfection')
outcome = "inhos_mortality"
time = "inhos_duration"
id.vars = c(time_fix_vars, outcome, time)소개
Marginal Structural Cox regression은 고급통계 기법으로, 생존분석에서 시간이 지남에 따라 변화하는 공변량(covariates) 또는 교란변수(confounders)를 통제/보정하기 위해 사용합니다.
예를 들어, 수혈이 중환자실에 입실한 환자들의 원내 사망(in-hospital mortality) 위험에 어떤 영향을 미치는지 파악한다고 해보겠습니다.
이 때, 중환자실에 입실해 있는 첫째날부터 셋째 날까지 반복적으로 측정된 항목의 수치, 수혈 시행 여부 등을 보정할 필요가 있습니다.
이럴 때 활용할 수 있는 것이 Marginal Structural Cox regression입니다.
원리
Marginal Structural Cox regression은 다음과 같은 특징을 갖습니다.
1. Cox regression Model
생존분석이기 때문에, 관찰기간과 임상 결과 등의 변수가 활용됩니다. (e.g. 입원 기간 내 사망 여부)
2. Time-varying 변수
실제 임상 환경에서 환자가 입원한 기간동안 환자에게 사용하는 약물이나 치료 등으로 인해 환자의 생존률이 달라질 수 있습니다.
따라서 Marginal Structural Cox regression은 시간이 지나도 변하지 않는 고정된(Time-fixed) 변수 뿐만 아니라, 기존의 전통적인 통계분석 방법으로는 다루기 어려운 시간에 따라 변화(Time-varying)하는 변수까지 보정하게 됩니다.
3. Marginal Structural Models (MSMs)
Marginal Structural Cox regression은 MSM의 방식에서 착안하여 IPTW(Inverse Probability of Treatment weighting)을 계산합니다. 이후 이 가중치를 Cox regression 분석에 적용합니다.
분석 사례
ICU에 입원한 환자가 수혈을 받은 뒤 3일정도 입원했을 때, 병원 내 사망률이 어떻게 나타나는지 확인하기
- Time-fixed 변수: AGE, SEX, BMI, Charlson Comorbidity Index, Hb ICU Day 1, SOFA ICU Day 1, SAPS3 ICU Day 1
- Time-varying 변수: ICU 입실 후 1~3일 간 Hb
- Outcome:
inhos_mortality - Time:
inhos_mortality
모델링은 다음과 같은 순서로 진행됩니다.
1. 변수 나누기
시간에 관계없는 변수와 시간에 따라 변하는 변수, 즉 Time-fixed 변수와 Time-varying 변수를 나눠줍니다.
2. 데이터 변환
Time varying 변수들이 있기 때문에, 긴 데이터의 형태로 변환합니다. melt() 함수의 measure.vars 인자에 길게 변환시킬 Time varying 변수들의 이름 또는 이름의 규칙을 넣어줍니다.
data_melt = melt.data.table(data = data,
id.vars= c('inhos_mortality','inhos_duration', time_fix_vars),
measure.vars = patterns("Transfusion_ICUD"),
value.name = 'treatment')
head(data_melt)variable이라는 열에 반복 측정된 값의 이름(Transfusion)과 일자(ICUD1~3)이 포함되어 있으므로, 두 개를 나누어 각각의 column을 구성합니다.
이 때 새롭게 만들어지는 treatment은 time varying 변수인 수혈 변수의 이름, time은 Cox regression에서 time으로 들어갈 변수입니다.
# 시간(Day 1,2,3 column 만들기)
data_melt[,c("v","time") := tstrsplit(variable,"_ICUD",2)]
data_melt[,c("v",'variable') := NULL]
setorder(data_melt, 'SubjectNo', 'time')
data_melt반복 측정되는 값이 treatment (수혈 여부) 뿐만 아니라 여러 개인 경우, 해당 변수에 대해 melt()를 진행한 후, merge()를 통해 데이터를 합쳐줍니다.
data_melt_hb = melt.data.table(data = data,
id.vars= c('inhos_mortality','inhos_duration', time_fix_vars),
measure.vars = patterns("Hb"),
value.name = 'Hb')
data_melt_hb[,c("v","time") := tstrsplit(variable,"_ICUD",2)]
data_melt_hb[,c("v",'variable') := NULL]
data_melt_final = merge(data_melt, data_melt_hb[,.(SubjectNo, time, Hb)], by=c("SubjectNo","time"))
data_melt_final3. IPTW 계산하기
treatment(수혈)에 Time-fixed, Time-varying 변수들이 미치는 영향력 확인하여 가중치를 계산합니다.
covars = c(setdiff(time_fix_vars, c('SubjectNo','hospital')), c('Hb'))
form = sprintf("treatment ~ factor(time) + %s", paste0(covars, collapse="+"))
lr = glm(as.formula(form), family=binomial, data = data_melt_final)
lr_numer = glm(treatment ~ factor(time), family=binomial, data=data_melt_final)
data_melt_final[, wt := predict(lr_numer, type="response", newdata = data_melt_final)/ predict(lr, type="response", newdata = data_melt_final)]4. Cox regression 수행
가중치(wt)를 적용하여 Cox regression을 수행합니다. robust = TRUE를 통해 가중치를 적용했을 때, 추가로 발생할 수 있는 변동성을 설명합니다.
form = sprintf("Surv(inhos_duration, inhos_mortality==1) ~ treatment + %s", paste0(covars, collapse="+"))
fit = coxph(as.formula(form),
data = data_melt_final,
weights = wt,
robust = TRUE
)
summary(fit)Call:
coxph(formula = as.formula(form), data = data_melt_final, weights = wt,
robust = TRUE)
n= 9898, number of events= 2862
(848 observations deleted due to missingness)
coef exp(coef) se(coef) robust se z
treatment -0.161344 0.850999 0.041316 0.098180 -1.643
AGE 0.024675 1.024981 0.001386 0.006205 3.976
SEX -0.111386 0.894593 0.028172 0.097745 -1.140
BMI 0.008627 1.008664 0.002505 0.008624 1.000
CFScore -0.097335 0.907252 0.007012 0.026466 -3.678
Charlson_Comorbidity_index_total -0.027526 0.972849 0.007519 0.025178 -1.093
Comorbidity_MOSIAC_SMT 0.426227 1.531469 0.038081 0.104275 4.088
CRRT_within_3days 0.778360 2.177898 0.040340 0.084549 9.206
input_preICU_to_ICUD1_quartile 0.005758 1.005775 0.012663 0.040529 0.142
SOFA_ICUD1 0.012404 1.012481 0.005710 0.017222 0.720
SAPS3_ICUD1 0.015565 1.015687 0.001388 0.004669 3.334
Septic_shock_ICUD1 0.127548 1.136040 0.031568 0.084565 1.508
s_pulmonary 0.392639 1.480884 0.038081 0.116431 3.372
s_abdominal -0.363302 0.695376 0.048356 0.108694 -3.342
s_urinary -0.769593 0.463202 0.050936 0.106356 -7.236
TypeOfInfection -0.895487 0.408408 0.044256 0.118998 -7.525
Hb -0.210082 0.810518 0.002839 0.033923 -6.193
Pr(>|z|)
treatment 0.100312
AGE 7.00e-05 ***
SEX 0.254471
BMI 0.317143
CFScore 0.000235 ***
Charlson_Comorbidity_index_total 0.274282
Comorbidity_MOSIAC_SMT 4.36e-05 ***
CRRT_within_3days < 2e-16 ***
input_preICU_to_ICUD1_quartile 0.887018
SOFA_ICUD1 0.471359
SAPS3_ICUD1 0.000857 ***
Septic_shock_ICUD1 0.131480
s_pulmonary 0.000745 ***
s_abdominal 0.000830 ***
s_urinary 4.62e-13 ***
TypeOfInfection 5.26e-14 ***
Hb 5.91e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
treatment 0.8510 1.1751 0.7020 1.0316
AGE 1.0250 0.9756 1.0126 1.0375
SEX 0.8946 1.1178 0.7386 1.0835
BMI 1.0087 0.9914 0.9918 1.0259
CFScore 0.9073 1.1022 0.8614 0.9556
Charlson_Comorbidity_index_total 0.9728 1.0279 0.9260 1.0221
Comorbidity_MOSIAC_SMT 1.5315 0.6530 1.2484 1.8787
CRRT_within_3days 2.1779 0.4592 1.8453 2.5704
input_preICU_to_ICUD1_quartile 1.0058 0.9943 0.9290 1.0889
SOFA_ICUD1 1.0125 0.9877 0.9789 1.0472
SAPS3_ICUD1 1.0157 0.9846 1.0064 1.0250
Septic_shock_ICUD1 1.1360 0.8803 0.9625 1.3408
s_pulmonary 1.4809 0.6753 1.1787 1.8605
s_abdominal 0.6954 1.4381 0.5620 0.8605
s_urinary 0.4632 2.1589 0.3760 0.5706
TypeOfInfection 0.4084 2.4485 0.3234 0.5157
Hb 0.8105 1.2338 0.7584 0.8662
Concordance= 0.972 (se = 0.025 )
Likelihood ratio test= 16858 on 17 df, p=<2e-16
Wald test = 742.7 on 17 df, p=<2e-16
Score (logrank) test = 39116 on 17 df, p=<2e-16, Robust = 556.6 p=<2e-16
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).