[R 통계] 성향 점수 매칭 (Propensity score matching) - 3

 

• 이전 글에 이어서 본 글에서는 매칭 균형 평가와 통계분석을 다루겠습니다.

 

4. 매칭 균형 평가

---

• 매칭 균형 평가는 매칭이 적절하게 잘 수행되었는지 평가하는 것으로 Balance diagnostics 라고 불립니다.
• 바로 직전 글의 말미에서 보여드린 것처럼 매칭 전에는 군 간 유의미한 차이가 있던 변수가 매칭 후 군간 유의미한 차이가 없어지는 것으로 '매칭이 작동은 했겠구나' 라고 짐작할 수는 있지만 정확한 방법은 아닙니다.
• 가장 정량적이고 직접적인 방법은 치료군과 비교군 간 공변량의 SMD (standardized mean differences) 와 VR (variance ratios) 을 비교하는 방법입니다. 쉽게 말해 군 간 공변량 차이를 보여주는 값이라고 생각하시면 됩니다.
• 문헌마다 차이가 있습니다만 SMD < 0.25 면 공변량이 균형되었다고 말하는 문헌도 있고, 그보다 좀 더 엄격한 기준인 SMD <0.1 을 주장하는 문헌도 있습니다. 
• VR은 1.0 에 가까울수록 공변량이 균형되었음을 의미하고, 0.5 미만 혹은 2.0 초과일 때 균형되지 않음을 의미합니다.

군 간 균형의 조건
SMD < 0.25 (or SMD <0.1)
0.5 < VR < 2.0

 

library(MatchIt)
library(cobalt)

psm1 = matchit(Smoke~ HE_glu + HE_HbA1c + HE_DM_HbA1c + HE_chol + HE_HDL_st2 + HE_TG + HE_LDL_drct + HE_HCHOL,
               data = data, distance = 'logit', method = 'optimal', replace = FALSE, ratio = 1)

psm2 = matchit(Smoke~ HE_glu + HE_HbA1c + HE_DM_HbA1c + HE_chol + HE_HDL_st2 + HE_TG + HE_LDL_drct + HE_HCHOL,
               data = data, distance = 'logit', method = 'nearest', replace = FALSE, ratio = 1, caliper = 0.2)


bal.tab(psm1, m.threshold = 0.1, v.threshold = 2, un = TRUE)
bal.tab(psm2, m.threshold = 0.1, v.threshold = 2, un = TRUE)

 

• 직전 글에서 만들었던 매칭인 PSM1과 PSM2 에 대하여 정량적인 매칭 균형 평가를 해보겠습니다. (일부러 균형되지 않은 부분을 만들어보려고 공변량의 개수를 조금 늘렸습니다.)
• m.threshold 는 SMD 의 제한 조건을 의미하고, v.threshold 는 VR을 의미합니다. 따라서 SMD < 0.1 , VR <2 를 확인하는 의도라고 보시면 됩니다.

 

• 먼저 PSM1 입니다. bal.tab 함수를 이용하면 이렇게 쉽게 각 공변량에 대하여 Balanced, Not Balanced 로 확인이 가능합니다.
• 전체 1478개의 쌍이 매칭되었고, HE_TG 변수에서 Not Balanced 가 나타났습니다.
• 불균형이 발생하는 경우, 매칭 알고리즘을 변경해볼 수 있고 혹은 불필요한 공변량은 없는지 확인해 공변량의 개수를 줄이는 방법도 있습니다.
• 해당 부분을 조금 교정해 보기 위해 PSM2에서는 Nearest 매칭 알고리즘을 사용하고, caliper 를 0.2로 지정했습니다.

 

 

• 직전 글에서 설명된 바와 같이 caliper 값을 지정했더니 전체 변수가 Balanced 되었지만, 그 만큼 매칭된 쌍의 수가 1394개로 줄어들었습니다.

 

 

• 정량적인 수치 말고도 Plot을 이용해서 균형 잡히기 전후의 변화를 확인할 수도 있습니다.

bal.plot(psm2, var.name = 'HE_TG', which = 'both', grid=TRUE)
bal.plot(psm2, var.name = 'HE_TG', which = 'both', grid=TRUE, type="ecdf")

 

 

 

 

• 이와 같은 bal.plot 을 이용하는 방법은 변수를 하나만 확인할 수가 있어서 논문에서는 잘 사용하지 않습니다. 그보다는 전체 covariate 균형 양상을 볼 수 있는 love.plot을 이용합니다.

 

love.plot(bal.tab(psm2, m.threshold=0.1))

 

 

• 이렇게 하면 매칭으로 조정되기 전에 각 공변량의 SMD 가 어땠는지 (빨간색 점), 이후 SMD가 0.1 안쪽으로 좁혀졌는지 시각적으로 확인이 가능합니다.
• 이건 가장 간단한 방법을 쓴 것이고, R 활용이 능숙하신 분이시라면 아시겠지만 아래와 같이 더 응용해서 논문에 실제 삽입할 Figure를 만드는 것도 가능합니다.

 

library(ggplot2)

new.names <- c(HE_glu = "공복혈당",
               HE_HbA1c = "당화혈색소",
               HE_DM_HbA1c = "당뇨병 유병 여부",
               HE_chol = "총 콜레스테롤",
               HE_HDL_st2 = "고밀도지단백 콜레스테롤",
               HE_TG = "중성 지방",
               HE_LDL_drct = "저밀도지단백 콜레스테롤",
               HE_HCHOL = "고콜레스테롤혈증 유병여부"
)


love.plot(bal.tab(psm2, m.threshold=0.1), stat = c("mean.diffs"), 
          grid=TRUE, stars="raw", binary= "raw", abs = TRUE, var.order = "unadjusted",
          var.names = new.names, drop.distance = TRUE, shapes = c("circle filled", "circle"),
          title = NULL,
          position = c(0.75, 0.25)) +
  theme(legend.box.background = element_rect(),
        legend.title = element_blank(),
        legend.box.margin = margin(1, 1, 1, 1))

 

 

• 함수 내 argument 에 대해 간단히 설명하면 아래와 같습니다.
• stats: "mean.diffs"가 기본이고 VR을 확인하고 싶다면 "variance.ratios" 를 입력하면 됩니다. stats = c("mean.diffs", "variance.ratios") 라고 입력해서 둘 다 보여주는 방법도 있습니다.
• abs: TRUE로 하면 절대값으로 표현할 수 있습니다.
• var.order: "adjusted"는 adjust 된 값을 기준으로 정렬하는 것이고, "unadjusted"는 unadjust 된 것을 기준으로 정렬한느 것입니다.
• drop.distance: TRUE를 입력해서 맨 윗 단에 있던 distane 값 삭제
• Standardized mean differences and raw mean differences are present in the same plot. Use the `stars` argument to distinguish between them and appropriately label the x-axis. 와 같은 메시지가 나오는 경우가 있는데 이 경우 stars="raw" 를 지정해주면 됩니다. (범주형 변수가 있는 경우 * 표시로 raw difference 임을 나타내줍니다.

 

5. 통계 분석

---

• 이제 매칭이 잘 되었으니 군 간 통계 분석을 수행할 단계입니다. 그 전에 먼저 ATE와 ATT의 개념에 대해 간단히 짚고 넘어가야 합니다.

처치 집단에 대한 인과 효과 (avearge treatment effect on the treated, ATT)
> 처치 집단이 치료를 받았을 때와 받지 않았을 때의 평균 치료 효과를 추정합니다.
> example: 금연 집단에 대한 금연의 효과를 추정

집단 전체에 대한 인과 효과 (average treatment effect, ATE)
> 전체 집단에 대해서 평균 치료 효과를 추정합니다
> example: 흡연 + 금연 집단에 대한 금연의 효과를 추정

 

• 아마 대부분 매칭 후 데이터 분석의 주요 목적은 ATT를 구하는 것일 겁니다. 
• 만약 ATE를 추정하고자 한다면 machit 함수에서 method = "full" 로 설정하여 optimal full matching 을 수행해, 일대일 매칭 혹은 일대다 매칭이 아니라, 모든 군이 완전 매칭되고 매칭에 대한 weight가 기록되도록 매칭을 수행해야 합니다.
• 본 글에서는 일단 ATT를 구한다고 가정하고 최대한 간단하게 해볼 수 있는 몇 가지 방법을 소개드리고자 합니다.

 

Paired 방법

• nonpaired 방법이 아니라 paired 방법으로 두 군을 비교하는 방법입니다.
• 가장 쉬운 예로 t-test 에서 indepedent t-test 가 nonpaired이고, paired t-test 가 paired 방법이라고 할 수 있습니다.
• 다만 이 경우, paired 이기 때문에 1:1 매칭에서만 사용이 가능하고, 1:k 매칭에서는 사용이 불가합니다.
• 아래와 같이, 흔히 모수적 방법으로는 paired t-test, 비모수적 방법으로는 Wilcoxon's sign rank test 가 있는데 이 부분은 다른 소개된 페이지가 많은 것 같아 넘어가도록 하겠습니다.

t.test(x, y, paired = TRUE, alternative = "two.sided")

wilcox.test(x, y, paired = TRUE, alternative = "two.sided")

 

 

 

선형 회귀적 방법

• 1:1이 아니라 1:k 에서도 쓸 수 있는 방법으로는 선형 모델에 데이터셋을 fitting 시켜 치료 효과를 추정(estimate)하는 방법입니다.
• Bootstrapping을 사용하는 방법도 있는데, 본 글에서는 간단하게 적용해볼 수 있는 아래 방법을 소개드리겠습니다.

 

library(MatchIt)
library(marginaleffects)

psm2 = matchit(Smoke~ HE_glu + HE_HbA1c + HE_DM_HbA1c + HE_chol + HE_HDL_st2 + HE_TG + HE_LDL_drct + HE_HCHOL,
               data = data, distance = 'logit', method = 'nearest', replace = FALSE, ratio = 1, caliper = 0.2)

data_m = match.data(psm2) ## 매칭된 데이터셋 생성


data_m$HBP <- as.numeric(as.character(data_m$HBP)) ## outcome인 HBP를 numeric으로 변경

fit = lm(HBP ~ Smoke * (HE_glu + HE_HbA1c + HE_DM_HbA1c + HE_chol + HE_HDL_st2 + HE_TG + HE_LDL_drct + HE_HCHOL),
          data = data_m, weights = weights) ## 선형 모델 생성

fit_c = avg_comparisons(fit, variables = "Smoke",
                         vcov = ~subclass,
                         newdata = subset(data_m, Smoke == "Smoker"),
                         wts = "weights")

fit_p = avg_predictions(fit, variables = "Smoke",
                         vcov = ~subclass,
                         newdata = subset(data_m, Smoke == "Smoker"),
                         wts = "weights")

 

• avg_comparison 은 군 간의 차이를 추정하는 방법이고, avg_prediction에서는 각 군의 추정값 보여 줍니다.
• 그러면 위와 같이 Non-smoker와 Smoker 일 때 estimated HBP 를 예측하고 비교할 수 있습니다. 흡연자의 경우 HBP의 비율이 더 높게 나타납니다.
• 다만 fit_c 에서 Z statistics 의 p 값이 0.109로 차이가 유의미하지는 않음을 확인할 수 있습니다.
• Figure로 표현한다면 아래와 같이 point estimates와 각 신뢰구간 상한값, 하한값을 이용해 표현할 수 있습니다.

 

ggplot(fit_p, aes(x = factor(Smoke), y = estimate)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  labs(x = "Smoking Group", y = "HBP", title = "Smoking Effect on HBP") +
  theme_minimal()

 

 

• 비교하고자 하는 Outcome이 연속형 변수냐 범주형 변수냐에 따라서 차이는 없습니다만, 범주형 변수인 경우 logistic regression model 에 fitting 할 것인지, 혹은 risk difference, realtive risk 를 계산할 것인지 등에 의해서 argument의 조정이 필요합니다.

 

• 이것으로 성향 점수 매칭과 R 분석에 관한 글은 마치겠습니다. 아래는 성향 점수 분석 중 흔히 생각할만한 고민에 대한 몇가지 Q&A입니다.
• 참고만 하시고, 상세한 내용을 알고 싶다면 아래 Reference 들을 보시길 추천드립니다.
• 저는 단순히 해당 Reference 에서 흔히 쓰일 법한 내용을 위주로 추린 것이라, 깊게 알고 싶으시다면 이것들을 보시는 것이 더 큰 도움이 되실 겁니다.

 

 

 

Q1. 공변량은 몇 개를 선택해야 적절할까?

>> 정해진 개수는 없지만, 관련된 분야의 연구를 참조하시길 추천합니다. 어느 정도 샘플 사이즈 (n) 수에서 몇 개의 공변량을 일치시키는지를 보면 판단이 조금 서실 것 같습니다. 공변량이 많아질 수록 아마 얻을 수 있는 균형 잡힌 쌍 (pair) 의 개수가 줄어들 것이고, 공변량이 너무 적으면 군 간 bias가 적절히 통제되지 않을 것 입니다.

가장 중요한 것은 앞선 글의 <2. 공변량 선택>에서 설명드린 원칙에 위배되지 않는 공변량을 선택하는 것입니다.

 

Q2. PSM 분석이 p-hacking을 유발하지 않을까?

>> p-hacking이란 쉽게 말해 통계 분석에서 유의미한 결과 (p <0.05 or p<0.01) 을 얻기 위해 데이터 자체나 데이터 분석 과정에서 임의의 조작을 가하는 것을 의미합니다. 혹은 정해진 계획 이상의 추가적인 통계 분석으로, 쉽게 말해 p 값이 적절히 나올 때까지 이것저것 해보는 것도 포함될 수 있습니다.

PSM은 그룹간 bias를 줄이려는 의도로 수행되기 때문에 그 자체가 p-hacking을 유발하진 않을 것 입니다. 하지만 PSM 수행 방법과 이후 통계 분석 방법을 미리 정해놓고 이것저것 반복적으로 해본다면, p-hacking이 생길 것입니다. 미리 방법을 정해놓고 분석에 들어가길 추천드립니다.

예를 들어,1) 가설에 따라 균형을 이뤄야하는 공변량은 무엇, 무엇인지, 2) PSM 후 불균형이 발생하면 매칭 알고리즘을 어떤 순서로 수정할 것인지, 3) 균형 잡힌 PSM 수행할 통계 분석 방법은 무엇인지 를 정할 수 있습니다.

 

 

[References]

1. Zhao, Qin-Yu, et al. "Propensity score matching with R: conventional methods and new features." Annals of translational medicine 9.9 (2021).
2. Randolph, Justus J., and Kristina Falbe. "A step-by-step guide to propensity score matching in R." Practical Assessment, Research & Evaluation 19 (2014).
3. Austin, Peter C. "An introduction to propensity score methods for reducing the effects of confounding in observational studies." Multivariate behavioral research 46.3 (2011): 399-424.
4. https://cran.r-project.org/web/packages/cobalt/vignettes/love.plot.html
5. https://cran.r-project.org/web/packages/MatchIt/vignettes/estimating-effects.html
6. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3110307/