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

 

• 이전 글에 이어서 본 글에서는 성향점수 평가와 매칭을 다루겠습니다.

 

3. 성향점수 평가 및 매칭

---

• 데이터 준비 및 결측치 처리, 공변량 선택이 끝났다면 이제 어떤 성향점수 평가 및 매칭을 시작할 차례입니다.
• 성향점수를 estimation 하는 방법 역시 여러가지가 있지만 본 글에서는 로지스틱 회귀분석(LR, Logistic regression) 만을 다루도록 하겠습니다.
• Generalized boosting model, random forest, artificial neural network 등 최근 관심도가 높은 머신러닝과 연관되는 방식도 있으나 로지스틱 회귀분석이 PSM에서는 주류로 알려져 있습니다. 이유는 단순하고 해석하기도 쉬워서입니다.

 

출처: https://ourworldindata.org/how-do-the-risks-of-death-change-as-people-age

• 성향점수 평가 시 LR을 사용한다면, 주의해야할 점이 있습니다. 바로 공변량과 결과변수의 관계가 선형이 아닐 때입니다.
• 예를 들어, 위 그림과 같이 나이 (공변량) 에 따른 결과 변수 (사망률) 을 가정해봅시다. 보통 질병에 취약한 그룹은 아주 어린 사람이나 나이가 많은 사람이고, 중간 나이의 사람들은 상대적으로 건강합니다. 따라서 나이는 증가함에 따라 사망률이 증가한다 혹은 감소할수록 감소한다는 식의 선형적 관계가 성립하지 않습니다. 이런 경우 LR을 이용한 PSM estimation 에 적절하지 않습니다.
• 이 경우 변수를 Factor화해서 어린 나이, 중간 나이, 고령 과 같은 식으로 범주화한다든가하는 고민이 필요합니다.

 

data$Smoke = factor(data$Smoke, exclude = c("",NA))
data$HBP = factor(data$HBP, exclude = c("",NA))
data$HE_DM_HbA1c = factor(data$HE_DM_HbA1c, exclude = c("",NA))
data$HE_HCHOL = factor(data$HE_HCHOL, exclude = c("",NA))
data$HE_HTG = factor(data$HE_HTG, exclude = c("",NA))

• 사실 위에서 설명한 경우에 해당하지는 않지만, 성향점수 평가에 앞서 범주형 변수는 모두 Facotr 함수를 써서 Factor화 했습니다. (간혹 가다가 NA가 아닌 공란도 데이터셋에 있는 경우가 있어 " " 와 NA 모두 제외했습니다.

 

 

• 성향점수 평가 및 매칭을 수행하기 위해서는 "Matchit" library를 사용합니다. 함수 내의 중요한 argument 에 대해 간단히 설명하면 아래와 같습니다.

library(MatchIt)

matchit(formula, data = NULL, method = "nearest", distance = "logit", 
	replace = FALSE, caliper = NULL, ratio = 1)

 

• formula: A ~ X1 + X2 + X3 와 같은 형식으로 작성합니다. (A는 치료 할당 변수 (T), X1부터 X3는 공변량 (C)를 의미합니다.) 
• data: 분석하고자 하는 데이터프레임
• method: 매칭 알고리즘을 선택하는 부분입니다. 주로 아래 3가지를 많이 사용합니다.
  "nearest" - nearest neighbor matching (기본값)
  "optimal" - optimal pair matching
  "full" - optimal full matching
• distance: PSM estimation 시 사용할 model을 지정합니다. "glm", "gam", "randomforest" 등 다양하지만 위에서 언급한대로 본 글에서는 LR 을 사용할 것으로 기본값인 "logit" 을 사용합니다.
• replace: TRUE (한번 매칭되었던 것도 다시 매칭될 수 있음) or FALSE (한번 매칭된 경우 다시 매칭되지 않음)
• caliper: 특정 distance 이상 멀어진 경우 matching 하지 않도록 설정하는 일종의 제한값입니다.
• ratio: 치료군과 대조군의 매칭 비율을 의미합니다. 기본값은 1로 1:1 비율로 매칭합니다. 예를 들어 3을 입력하면 1:3으로 매칭됩니다.

https://cran.r-project.org/web/packages/MatchIt/MatchIt.pdf

• 자세한 내용이 궁금하다면 위 Matchit 패키지 문서를 확인하시면 됩니다.
• 지금부터는 매칭 알고리즘이 각각 어떤 차이가 있는지 설명드리겠습니다.

 

 

1) Nearest neighbor matching (NNM)

• 가장 기본적인 방법은 NNM 입니다.
• 치료군 환자에 대해서 하나 하나씩 가장 가까운 대조군 환자를 고릅니다.
• 예를 들어 위와 같이 1번 환자에 가장 가까운 대조군을 하나 고르고, 나머지 중에서 2번 환자에 가장 가까운 환자를 고르는 식입니다.

 

• NNM with replacement 입니다. 함수 내 argument 중 replace = TRUE 로 설정하면, 한번 매칭되었던 대조군도 다시 매칭될 수 있습니다.
• 다만 중복되는 대조군이 생기므로 추후 통계분석 시 weight를 고려해 판단해야합니다.  주로 깔끔하게 Pair가 생기는 것을 선호하다보니 연구에서 잘 쓰이지 않는 것 같습니다.

• NNM without replacement 에서 문제가 발생하는 경우는 위와 같은 경우입니다.
• 1이 먼저 가장 가까운 대조군 환자를 고름으로서, 2는 아주 distance가 먼 (matching 측면에서 좋지 않은) 대조군 환자를 고르게 됩니다.

 

2) optimal pair matching (OM)

• NNM without replacement와는 다르게 전체 매칭의 distance를 고려하여 치료군과 대조군이 전반적으로 최적 distance로 매칭되는 알고리즘입니다.
• NNM보다 나은 경우도 있으나, 공변량의 균형 측면에서는 NNM이 더 낫다고 알려져 있습니다.

 

3) NNM with a caliper

• NNM을 사용하지만 일정한 caliper 를 정하여 그것을 초과하는 매칭은 수행하지 않는 알고리즘입니다.
• 연구에 따르면 통상적으로 0.2 정도가 많이 쓰입니다. (Austin PC. Optimal caliper widths for propensity-score matching when estimating differences in means and differences in proportions in observational studies. Pharm Stat 2011;10:150-61.)
• caliper 를 설정한다면 균형잡힌 pair를 얻을 수 있으나, 누락되는 케이스가 생겨 매칭된 pair의 수가 적어질 수 있습니다.

 

4) 1:k matching

• 함수 내 argument 중 ratio를 조정하여 하나의 치료군이 여러 대조군을 선택하도록 하는 방법입니다.
• 대조군이 적을수록 성향점수가 유사한 짝이 생기므로 편향은 작아지지만 데이터의 크기가 작아 분산이 커집니다.
• 대조군이 많을수록 성향점수가 덜 유사한 짝이 생기므로 편향은 커지고 데이터의 크기가 커져 분산은 작아집니다.
• 이전 연구에 따르면 1:2, 1:3 매칭까지는 괜찮지만 1:5 이상으로 비율이 커지면 큰 이득이 없다고 합니다.

 

 

 

• 예제에 적용하기 앞서 본 실습에서는 가상으로 설정한 연구 목적을 다시 말씀드리면, 아래와 같습니다.

흡연 여부에 따라 고혈압 진단 유무에 차이가 있는지 평가한다.

흡연 여부에 따른 군 할당 시, 혈액 검사 변수를 공변량으로 두고 성향 점수 매칭한 "흡연군" vs "비흡연군" 을 나누어 고혈압 진단 유무 차이를 분석한다.

 

요소	설명	비고
공변량 (C)	혈액 검사 변수	혈당, 콜레스테롤, 중성지방 등
치료 할당 변수 (T)	흡연 여부	평생 5갑 이상, 없음
결과 변수 (O)	고혈압 진단 유무	있음, 없음

 

 

 

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

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

• PSM1 에서는 Optimal pair matching 을 이용했습니다.
• PSM2 에서는 NNM 을 이용했고 caliper로 0.2 를 지정해보았습니다.

 

mytable(Smoke~HE_glu + HE_HbA1c + HE_DM_HbA1c, data = data)

matched_data1 = match.data(psm1)
mytable(Smoke~HE_glu + HE_HbA1c + HE_DM_HbA1c, data = matched_data1)


## moonbook, ztable 이용한 viewer
library(ztable)
library(magrittr)

mytable(Smoke~HE_glu + HE_HbA1c + HE_DM_HbA1c, data = data) %>%
  ztable %>%
  addSigColor %>%
  print(type = "viewer")

mytable(Smoke~HE_glu + HE_HbA1c + HE_DM_HbA1c, data = matched_data1) %>%
  ztable %>%
  addSigColor %>%
  print(type = "viewer")

 

• PSM을 수행한 후 match.data 함수를 이용하면 매칭된 케이스만 남기고, distance, weights 값이 붙어 있는 데이터셋을 획득할 수 있습니다.

 

• 위은 코드를 이용해 테이블로 매칭 전후 그룹 차이를 확인해보겠습니다.
• (아래 있는 코드는 moonbook, ztable, magritrr 이 필요한 방법인데, 텍스트가 아닌 뷰어로 표를 확인하기 좋은 방법이라 덧붙여두었습니다.)

PSM 수행 전 데이터

PSM 수행 후 데이터

 

 

• HE_glue, HE_HbA1c 등의 변수에서 그룹 간 유의미한 차이가 있었는데, 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.