인과성 추론에서 성향점수 매칭에 대한 비판적 고찰

A Critical Review of Propensity Score Matching in Causal Inference

Article information

J Health Info Stat. 2022;47(Suppl 1):S9-S19
Publication date (electronic) : 2022 May 31
doi : https://doi.org/10.21032/jhis.2022.47.S1.S9
1 Graduate Student, Department of Public Health Sciences, Graduate School of Public Health, Seoul National University, Seoul, Korea
2 Professor, Department of Public Health Sciences, Graduate School of Public Health, Seoul National University, Seoul, Korea
유지웅1orcid_icon, 이우주,2orcid_icon
1 서울대학교 보건대학원 보건학과 석사과정생
2 서울대학교 보건대학원 보건학과 교수
Corresponding author: Woojoo Lee 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Korea E-mail: lwj221@gmail.com
No potential conflict of interest relevant to this article was reported.
Received 2022 March 26; Revised 2022 May 24; Accepted 2022 May 31.

Trans Abstract

Propensity score matching (PSM) is one of the most widely-used causal inference methods to estimate the causal estimands such as average treatment effect or average treatment effect on the treated from observational studies. To implement PSM, a researcher first selects an appropriate set of confounders, estimates the propensity score, and matches the treated group with the control group using a matching algorithm such as nearest neighborhood or optimal matching. In this paper, we highlight the importance of investigating the assumptions employed in the PSM procedure thoroughly because they strongly affect the analysis result, but are not testable using observational data. We explain how to exploit the domain knowledge to avoid the potential risks from the violation of the untestable assumptions, and show how the research purpose is linked to selecting the matching algorithm and downstream analysis after PSM. In addition, to examine the vulnerability of the causal result, we highlight the use of sensitivity analysis for the analysis after PSM. These points are demonstrated in detail using National Supported Work data.

서 론

특정 위험인자가 질병 발생에 미치는 영향이나 특정 질병에 대한 치료제의 효과처럼 처치(treatment)와 결과(outcome) 사이의 인과관계를 파악하고 그 효과를 추정하는 것이 여러 역학연구와 임상연구의 주요 목적이다. 이러한 목적을 달성하기 위한 가장 높은 수준의 증거를 갖추는 방법은 무작위 대조 시험(randomized controlled trial, RCT)으로 개개인의 특성과 관계없이 무작위로 처치를 배정하기 때문에 편향이 없는 인과 효과를 추론할 수 있다[1]. 하지만 처치를 무작위로 배정하는 것은 윤리적 측면에서 문제가 생길 수 있으며 결과를 모으는 데 시간이 걸리는 경우가 많아 RCT를 항상 구현할 수는 없다. 이러한 이유로 관찰연구(observational study)에서 인과관계를 추론하는 것이 큰 관심을 받고 있다[2].

RCT에서는 연구목적에 맞춰 설계된 연구에서 데이터를 얻지만 관찰연구에서는 이미 수집된 데이터를 연구 목적에 맞게 분석하기 때문에 인과효과가 왜곡될 수 있다[2]. 관찰연구에서는 처치의 배정에 영향을 주는 처치 전 변수로 인해 실제로 처치와 결과 사이에 인과관계가 없더라도 둘 사이에 상관관계가 존재할 수 있다. 따라서 연구자는 거짓 인과관계를 마치 실제 존재하는 인과효과로 해석하여 함정에 빠지는 것에 각별히 유의해야 한다.

관찰연구에서 올바른 인과효과를 추론하기 위해 회귀모형, 가중(weighting), 도구변수, 매칭(matching) 방법 등 다양한 방법들이 연구 되었으며 그 중에서도 매칭은 직관적이기 때문에 매우 널리 사용되고 있는 인과추론 방법이다[24]. 대표적인 매칭방법은 실험군의 처치 전 변수와 같은 값을 갖는 대조군을 찾아 짝짓는 것이다. 이는 양 쪽 그룹에 처치 전 변수가 완벽하게 똑같은 분포를 가지게 하므로 마치 무작위 배정한 것과 유사한 효과를 가지게 된다[3,5]. 이러한 매칭방법의 이론적 연구는 1970년대에 Cochran과 Rubin 등에 의해 이뤄졌으며 인과효과에 편향을 줄이기 위한 노력으로 이해되었다[2,6]. 하지만 이러한 매칭방법은 보정해야 할 교란변수가 많아질수록 실험군과 대조군에서 완전히 일치하는 개체를 찾기 어렵다는 문제점이 있어서 실제 사용에 어려움이 많다[5]. 이를 극복하기 위해 Rosenbaum and Rubin [7]은 성향점수(propensity score)를 제안하였으며 이후 성향점수를 사용한 매칭에 대해 연구가 활발히 이뤄졌다[813]. 인과추론에서 성향점수 매칭의 사용은 여러 이점을 지니고 있다. 특히 성향점수는 고차원의 교란변수를 1차원으로 축소한 것으로 교란변수의 수가 많거나 연속형 변수가 있을 때에도 짝지을 수 있는 개체를 찾기 쉬워 데이터를 좀 더 효율적으로 사용할 수 있다. 하지만 성향점수 매칭이 올바른 인과관계를 추론해내기 위해서는 여전히 데이터로 검증하기 어려운 가정이 요구되며 많은 응용연구에서 이러한 과정에 대한 검토가 적절히 이루어지지 않고 있다[4].

먼저 올바른 인과관계를 추론하기 위해서는 교란변수를 적절히 보정해야 하는 데 이 때 보정해야 할 변수의 집합을 찾는 유용한 도구로 방향성 비순환 그래프(directed acyclic graph, DAG)가 활용될 수 있다[14]. 또한 연구자는 연구 목적을 우선적으로 설정하고, 연구 대상을 고려하여 인과효과를 추론해야 한다. 예를 들어 인구집단 전체에 대한 인과효과(average treatment effect, ATE)와 처치를 받은 집단에 대한 인과효과(average treatment effect on the treated, ATT)는 서로 다른 의미를 가지므로 확실히 구분되어야 한다[15]. 어떠한 매칭 알고리즘을 선택하느냐에 따라 ATE가 구해지기도 하고 ATT가 구해지기도 하기 때문에 연구 대상에 맞는 알고리즘 선택은 중요한 문제이다. 또한 데이터 분석 결과에 대해 민감도 분석을 실시하여 데이터로 검증되지 않는 가정이 일부 위반되었을 때 현재의 결과가 얼마나 강건(robust)한지 확인하는 것이 필수적이라 할 수 있다.

본 논문의 2장과 3장에서는 각각 성향점수와 DAG에 대한 정의와 개념을 소개하고, 4장에서는 매칭 알고리즘에 대한 설명과 알고리즘 선택의 기준을 설명한다. 5장에서는 성향점수 모형 및 균형 진단, 6장에서는 짝지어진 데이터의 분석방법, 7장에서는 분석된 결과에 대한 민감도 분석 방법을 각각 소개한다. 8장과 9장에서는 성향점수 매칭 방법 자체에서 주의할 점과 최근 연구들에 대해 소개하고자 한다. 또한 성향점수 매칭을 이용한 분석을 구체적으로 설명하기 위해서 Lalonde [16]에서 사용된 국가 지원 사업 데이터를 통해 각 단계별로 분석 방법과 관련된 통계적인 이슈를 설명하고자 한다.

성향점수의 정의와 개념

실험군과 대조군이 무작위로 배정되는 무작위 대조 시험과는 달리 관찰연구에서는 각 개체가 지닌 속성들에 따라 처치의 배정이 영향을 받는다[6]. Rosenbaum and Rubin [7]은 관찰연구에서 이러한 개체들의 속성으로 인해 발생되는 편향을 제거하기 위해 성향점수를 제시하였다. 처치변수 Ti는 개체 i가 처치를 배정받으면 1 아니면 0을 나타내며 성향점수는 교란변수 Ci가 주어졌을 때 처치를 배정받을 조건부 확률 e(Ci)=Pr(Ti=1 │ Ci)로 정의된다. 성향점수는 균형점수(balancing score)의 일종으로 올바르게 추정하였다면 이를 보정하였을 때 편향되지 않은 인과효과를 추론할 수 있다는 성질을 갖는다[7].

성향점수를 이용한 인과효과 추론 과정을 이해하기 위해서는 우선 잠재변수(potential outcome) [17]라는 개념을 도입해야 한다. 개체 i가 처치를 받았을 때의 결과 Yi(1)라 하고, 처치를 받지 않았을 때의 결과를 Yi(0)라고 할 때 Yi(1), Yi(0)를 묶어서 잠재변수라고 부른다. 개체 i에 대한 처치의 인과효과는 두 잠재변수의 차이(Yi(1)−Yi(0)) 또는 비율(Yi(1)/Yi(0))로 표현할 수 있지만 두 잠재결과를 동시에 관측하기가 어려워서 개인의 인과효과를 구하기는 매우 어렵다. 반면 연구의 대상이 되는 인구집단에서의 인과효과는 일반적으로 아래의 세 가지 가정을 만족하면 데이터를 통해 추론할 수 있다.

첫 번째는 stable unit treatment value assumption (SUTVA) [18] 가정으로 어느 한 개체의 잠재 결과는 다른 개체에 배정된 처치에 영향을 받지 않고, 처치에 여러가지 형태가 존재하지 않아야 한다는 것이다.

가정 1. SUTVA : 만약 Ti=Ti´ 라면, Yi(Ti)=Yi(Ti´) 이다.

두 번째는 교환가능성(exchangeability) 가정으로 처치가 개체들이 가진 속성과 무관하게 무작위로 배정되었다면 잠재변수와 독립(Ti (Yi(1), Yi(0)))이 성립하여 실험군과 대조군의 결과의 평균의 차이를 통해 인과효과를 추론할 수 있다. 하지만 관찰 연구에서는 처치가 무작위로 배정되지 않아 조건부 교환가능성 가정이 제시되는데 이는 적절히 선택된 교란변수의 집합 Ci가 주어졌을 때 잠재변수와 처치가 독립이라는 가정이다.

가정 2. 조건부 교환가능성 가정 : Ti(Yi(1), Yi(0))|Ci

세 번째는 양수성(positivity) 가정으로 교란변수가 주어졌을 때 처치를 배정받을 확률이 0과 1 사이라는 것이다.

가정 3. 양수성 가정: 0<e(Ci)=Pr(Ti=1 │ Ci)<1

위의 세 가지 가정이 모두 성립되면 교란변수 대신에 성향점수가 주어지더라도 조건부 교환가능성 가정이 성립한다[7]. 이때 성향 점수를 이용하는 방식에 따라 가중, 회귀, 매칭, 층화 등의 이름으로 인과효과 추론에 사용되고 있다. 또한 위의 세 가지 가정하에서 성향점수가 주어지면 교란변수와 처치 변수는 독립이 되어 교란변수의 분포가 처치 여부에 관계없이 동일함을 확인함으로써 성향점수를 이용한 인과추론 방법의 적절성을 평가해볼 수 있다.

본 논문에서는 성향점수 매칭 적용과정을 자세히 설명하기위해 1975년 3월부터 1977년 6월까지 미국에서 국가지원사업으로 시행된 국가 지원 데이터[16]를 이용한다. 이 데이터는 직업훈련 프로그램의 효과를 추정하기 위해 분석되었으며 이후 성향점수 매칭을 적용한 많은 논문에서 예시로 사용되고 있다[11,19]. 당시 프로그램은 자격을 갖춘 참가자들 중 일부를 선택하여 배정하였으며 프로그램 전후에 실험군과 대조군의 참가자들에게 인터뷰를 통해 데이터를 수집하였다. 논문에서는 미응답자들을 제외한 데이터[11]를 사용하여 총 445명 중 185명은 실험군, 나머지 260명은 대조군에 배정되었다. 처치변수는 실험군 여부를 의미하는 이분변수이며 결과변수는 78년도의 수입이다. 분석에서 사용된 교란변수는 연령, 교육받은 기간, 흑인, 히스패닉, 결혼 여부, 고등학교 학위 유무, 74년도 수입, 75년도 수입이 사용되었다. 분석에서 사용한 데이터는 다음 인터넷 주소에서 제공하고 있다(https://users.nber.org/~rdehejia/data/.nswdata2.html).

방향성 비순환 그래프를 이용한 변수선택

올바른 인과관계를 추론하기 위해서는 조건부 교환가능성 가정을 만족시키는 적절한 교란변수의 집합을 찾아야 한다. 모형에 포함할 설명변수를 선택하는 전통적인 통계방법으로는 p−값, Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC)와 같은 기준으로 모형을 선택하거나, LASSO나 Elastic Net과 같은 방법으로 자동으로 회귀계수를 0으로 만드는 방법들이 고려되어 왔다[20]. 하지만 이러한 방식의 변수 선택은 인과효과를 추정할 때 편향을 일으킬 수 있어 많은 문헌에서는 변수들 간의 인과관계를 표현한 그래프인 DAG를 통해 모형에 포함할 변수를 선택하는 것을 권장하고 있다[14].

DAG는 변수들의 인과관계를 노드와 화살표로 표현한 그래프로 노드는 변수, 화살표의 방향은 변수들 간의 인과 관계를 나타내며 화살표가 시작하는 변수가 원인, 화살표를 받는 변수가 결과를 의미한다. 인과관계가 있을 것으로 예상되는 노드들끼리 단방향 화살표로 연결하며 이 때 C→ B→ A→ C와 같이 순환하는 그래프가 형성되지 않도록 주의해야 한다. 노드를 잇는 화살표의 방향을 결정하는 것은 데이터를 통해 추론하는 것이 일부 가능하나[21] 변수들 사이에 화살표를 결 하기 어렵고 일부 교란변수가 누락될 수 있어 해당 분야의 전문가와 협업하여 결정하는 것이 필요하다[22]. DAG를 통해 변수를 선택할 때는 그래프에서 어떤 변수가 처치, 결과변수인지 결정하는 것이 중요한데, 이 설정에 따라 보정해야 할 변수가 달라지기 때문이다. 먼저 Figure 1A는 T, Y, C의 인과관계를 나타낸 그래프로 T와 Y는 인과관계가 없으며 C는 T와 Y 모두에 영향을 주는 것을 나타낸다. 이 경우 C를 보정하지 않으면 T← C→ Y 경로를 통해 T와 Y 사이에 거짓 인과관계가 발생하는데 T← C→ Y 경로를 뒷문 경로(backdoor path)라 한다. 이때 C는 뒷문 경로를 막기 위해 보정해야 할 변수에 해당한다[23]. Figure 1B는 T, Y, S의 인과관계를 나타낸 그래프로 T와 Y 사이에 인과관계가 없으며 S는 T와 Y 모두에 영향을 받는 것을 나타낸다. 이 경우 S를 보정하지 않아야 올바르게 T의 인과 효과를 추론할 수 있는데, S를 보정하면 막혀 있던 T→ S← Y 경로에 뒷문 경로가 열려 거짓된 인과효과가 나타나게 된다[14]. 이처럼 보정하였을 때 뒷문 경로가 열리게 하는 변수를 충돌변수(collider)라고 한다. 이와 같은 규칙으로부터 변수들 사이의 조건부 독립성을 성립하게 하는 그래프 조건을 ‘ D-separation’ [14]이라 하며 인과 추론을 위해 보정해야 할 변수를 선택할 때 널리 활용되고 있다. 이와 관련한 자세한 설명은 Pearl [24]를 참고하기를 권장한다. 또한 DAG를 그리면 자동으로 보정해야 하는 교란변수의 집합을 찾아주는 DAGitty (http://dagitty.net)라는 소프트웨어가 있어 복잡한 변수 간의 관계에서도 손쉽게 교란변수를 선택할 수 있다[25].

Figure 1.

Three directed acyclic graphs demonstrating various relation-ships of variables. Let T, Y, C, S, U denote treatment, outcome, confounder, collider, unmeasured confounder, respectively. The first two DAGs repre-sent no causal effect from treatment to outcome with (A) an open back door path (T← C→ Y) by confounder C and (B) no open backdoor path. (C) A DAG represents the causal relationship in the National Supported Work data. DAG, directed acyclic graph.

Figure 1C는 국가 지원 사업 데이터에서 변수들의 관계를 나타낸 DAG이다. 직업훈련 프로그램(T)은 78년도 수입(Y)에 인과효과를 가질 수 있다고 보며 세 가지 교란변수의 집합의 영향을 받는데 첫번째 집합(C1)은 나이, 교육받은 기간, 흑인, 히스패닉, 결혼 여부, 고등학교 학위 유무가 포함되고, 두 번째 집합(C2)은 74년도와 75년도 수입이 포함되어 있으며, 세 번째 집합(U)은 74년도와 75년도의 수입과 78년도의 수입 모두의 원인이 되는 가상의 변수들의 집합을 의미한다. C1은 T와 Y의 공통 원인일 뿐만 아니라 개인별 소득에도 영향을 주는 변수이기 때문에 C2의 원인으로 가정하였다. C2 역시 프로그램 참가와 78년의 수입에 모두 영향을 줄 것으로 판단되어 공통 원인으로 가정하였다. U는 측정되지 않은 교란변수로 보정해야 하지만 분석에 포함되지 않은 변수를 의미한다. Figure 1C의 DAG에는 총 4개의 열려 있는 뒷문 경로(T← C1→ Y, T← C2→ Y, T← C1→ C2→ Y, T← C2← U→ Y)가 있는데 C1과 C2를 모두 보정하면 4개의 뒷문 경로를 모두 막아 올바른 인과효과를 추정할 수 있다.

성향점수 매칭 알고리즘과 선택

성향점수 매칭 방법은 성향점수를 사용하여 각각 실험군과 대조군에서 유사한 성향점수를 갖는 개체들끼리 쌍을 이루는 방법이다. 매칭을 사용할 때는 연구 목적에 맞는 인과효과를 추론하도록 매칭 방법을 고려해야 한다. 매칭 방법에는 크게 완전 매칭과 일대일 매칭, 일대다 매칭이 있다. 먼저 완전 매칭은 실험군 또는 대조군에서 하나의 개체를 추출하고 나머지 군에서도 하나 이상의 개체를 뽑아 하나의 쌍을 이루게 하는 방법으로 데이터를 버리지 않고 전체를 사용한다[26,27]. 실험군과 대조군의 모든 자료를 사용하기 때문에 완전 매칭은 집단 전체에 대한 인과효과인 ATE(=E[Yi(1)−Yi(0)])를 추론할 때 주로 사용된다. 이와는 달리 일대일 매칭과 일대다 매칭은 실험군을 기준으로 대조군을 짝짓는 방법으로 쌍을 이루지 못한 개체의 자료는 분석에 사용되지 않는다. 짝을 짓는 방법은 실험군에서 하나의 개체를 추출하고 이와 가장 근접한 성향점수를 가진 대조군의 개체를 일대일 매칭은 1개, 일대다 매칭은 둘 이상의 고정된 수만큼 뽑아 하나의 쌍을 이루게 하는 것이다. 이와 같이 실험군을 기준으로 매칭이 이루어진 방법은 처치를 받은 집단에서의 인과효과인 ATT(=E[Yi(1)-Yi(0)│ Ti=1])를 추론할 때 주로 사용된다[15,28].

국가 지원 사업 데이터에서 프로그램의 대상자 전체에 대한 프로그램의 효과를 추정하고자 한다면 완전 매칭을 하여 ATE를 추정하고, 대상자들 중 프로그램을 배정받은 참가자들에 대해 프로그램의 효과를 추정하고자 한다면 일대일 매칭 또는 일대다 매칭을 하여 ATT를 추정하는 것이 일반적이다. 실험군에 짝지어지는 대조군의 비율을 선택할 때는 인과효과 추정치의 편향과 분산 사이에 상충관계를 고려해야 한다. 대조군의 비율이 작을수록 실험군과 성향점수가 유사한 개체가 짝지어지기 때문에 편향은 작아지지만 짝지어진 데이터의 크기가 작아 분산이 커지며, 대조군의 비율이 커지면 분산은 작아지지만 성 점수가 덜 유사한 개체가 짝지어지면서 편향이 커지게 된다[29]. 따라서 일대다 매칭을 무조건 선호해야 할 이유는 없다.

매칭에서 알고리즘 선택은 짝지어진 자료에서 실험군과 대조군 간에 교란변수의 균형에 영향을 줄 수 있다[29]. 완전 매칭에서는 쌍을 이루는 개체들 간의 성향점수 차이의 평균을 최소로 하는 쌍을 찾는 최적합 방법[26,27]을 사용한다. 한 개체에 짝지어지는 다른 군의 수를 제한하지 않으면 특정 쌍에 너무 많은 개체들이 포함될 수 있어 짝을 지을 수 있는 최대 개체 수를 10개 이내로 제한하는 것을 권장하고 있다[26]. 일대일 매칭과 일대다 매칭에서는 최적합 방법이나 실험군 개체를 성향점수 관점에서 가장 근접한 대조군의 개체와 짝짓는 최근접 이웃 방법을 사용한다. 최근접 이웃 방법에서는 이미 사용한 대조군을 다시 사용하는 복원추출 방법을 고려할 수 있으며, 복원추출 하는 경우 실험군과 성향점수가 유사한 대조군의 개체가 여러 번 사용되어 인과효과 추정치의 분산은 커지지만 편향은 줄어들 수 있다[29]. 또한 캘리퍼(caliper)를 설정하면 두 개체의 성향점수 차이가 일정 수준 δ보다 큰 개체들은 쌍을 이루지 못하도록 제약을 두어 교란변수의 균형을 개선시킬 수 있다. 캘리퍼를 설정하는 것은 추정치의 편향을 줄이는데 효과적이라고 알려져 있지만[29] 본래 추정하고자 하는 인구집단의 인과효과가 아닌 성향점수가 겹치는 집단에서의 인과효과(average treatment effect on the overlap, ATO)로 그 의미가 바뀔 수 있다. 특히 ATO 는 임상적으로 처치의 효과를 판단하기 어려운 집단에서 인과효과와 밀접한 관련이 있음이 알려져 있다[30]. 이처럼 캘리퍼를 사용하면 짝지어지지 않은 데이터를 이후 분석에서 사용하지 않기 때문에 연구자는 교란변수의 균형과 인과효과의 의미 사이의 상충관계에 상당한 주의가 필요하며 짝지어진 데이터가 관심 인구집단을 적절히 대표하는가에 대해 고민해야 한다[31].

성향점수 모형 및 균형 진단

성향점수는 데이터를 통해 추정하기 때문에 적절한 모형과 변수의 형태를 고려하여 올바른 성향점수 모형을 만들어야 한다. 성향점수 모형이 참이라면 같은 성향점수 값을 가지고 있는 실험군과 대조군의 교란변수는 같은 분포를 가져야 한다[7]. 하지만 DAG를 통해 선택한 교란변수를 가지고 성향점수 모형을 만들더라도 교란변수의 비선형항 등을 적절히 고려하지 않으면 모형이 잘못 설정될 수 있다. 성향 점수 모형이 잘못 설정되면 짝지어진 데이터에서 두 군의 교란변수가 균형을 이루지 않아 편향된 인과효과를 추론할 수 있다[4]. 따라서 짝지어진 데이터에 대해 균형을 진단하는 것은 적절한 성향점수 모형을 찾기 위해 매우 중요한 의미를 갖는다[4,28,32].

성향점수를 추정할 때 기본이 되는 모형은 로지스틱 회귀모형이다[7]. 로지스틱 모형에 비선형적 관계가 있을 것으로 예상된다면 고차항이나 교호작용을 포함시키는 것이 필요한데, 이를 적절히 고려하지 못하면 짝지어진 데이터의 균형이 이뤄지지 않아 편향된 인과효과가 추론될 수 있다[33]. 이에 대해 결정 나무(decision tree) 등 여러 비모수적인 기계학습 모형들은 연구자가 직접 변수들 간의 관계를 고려하지 않고서도 자동적으로 고차 교호작용 등을 고려해주므로 전통적인 회귀 모형의 대안으로 고려되었다[34,35]. 특히 성향점수의 참인 모형에 복잡한 비선형 관계가 있는 경우에 boosted CART, random forest 등이 로지스틱 모형에 비해 두 군의 균형을 더 잘 맞추어서 인과효과에 대한 편향과 평균 제곱 오차가 더 작다고 연구된 바 있다. 또한 최근 심층 신경망(deep neural network) [36]과 생성적 대립 신경망(generative ad-versarial network) [37] 등의 딥러닝(deep learning) 방법들도 성향점수 모형을 만들기 위한 방법으로서 연구되었다.

위의 방법들을 통해 성향점수 모형을 만들었다면 짝지어진 데이터에서 두 군 간 교란변수의 균형을 진단하여 성향점수의 적절성을 판단해야 한다. 균형 진단을 위해서는 두 군의 교란변수의 분포를 비교해야 하지만 교란변수의 결합분포를 비교하는 것은 매우 어렵기 때문에 다양한 요약 통계량과 그림을 통해 각각의 주변확률 분포를 비교하는 것으로 균형을 진단하는 것이 일반적이다. 가장 간단하게는 평균, 최소/최대값, 사분위수 등의 요약 통계량을 비교하여 진단할 수 있다. 그 중 평균을 비교하기 위해서 표본 수에 영향을 받는 t 검정이나 Kol-mogorov-Smirnov 검정보다는 표준화 평균 차이가 많이 쓰이고 있으며 두 군에서 교란변수의 표준화 평균 차이의 절댓값이 0.1보다 작다면 평균에 차이가 없다고 판단한다[32,38]. 연속형 변수와 이분 변수에 대한 표준화 평균 차이는 다음과 같이 정의한다.

X¯1X¯2(s12+s22)/2(X¯1,X¯2   ,s12,s22   )p^1p^2p^1(1p^1)+p^2(1p^2))/2(p^1,p^2   )

또한 교란변수 분포를 비교하기 위해 연속형 변수의 경우 각 교란변수들의 분산을 비교하기도 하는데, 분산의 비율이 대략 0.5와 2 사이에 있거나 F 검정에서 유의하지 않다면 분산에 차이가 없다고 판단한다[10]. 또 다른 방법으로 Stuart et al. [39]은 짝지어진 데이터에서 교란변수를 설명변수로, 결과변수를 종속변수로 한 회귀모형을 만든 뒤 실험군과 대조군에서의 결과를 각각 예측하여 표준화 평균 차이를 구한 것을 예후점수로 정의하여 예후점수의 절댓값이 0.1 이하일 때 균형을 이룬다고 판단한다. 교란변수의 수가 많으면 여러가지 그림을 통해 균형을 진단하는 것이 효과적일 수 있다. 연속형 변수의 요약 통계량과 분포는 상자 그림(box plot)과 막대 그림(bar plot) 등을 통해 쉽게 비교할 수 있으며 교란변수의 주변확률 분포를 비교하기 위해서 실험군을 y축, 대조군을 x축으로 한 Q-Q 그림과 경험적 누적 분포 등을 통해 효과적으로 비교할 수 있다[40,41]. R에서는 ‘ MatchIt’ 패키지[42]를 통해 매칭과 짝지어진 데이터에 대해 균형진단을 할 수 있다. matchit() 함수는 매칭 알고리즘과 다양한 성향점수 모형을 선택할 수 있으며 캘리퍼의 설정 여부도 결정할 수 있다. matchit() 함수의 객체에는 짝지어진 데이터와 함께 기술적으로 균형을 진단할 수 있는 다양한 정보들을 제공하는데 표준화 평균 차이, 분산 비율 등을 비교할 수 있다. 또한 ‘ cobalt’ 패키지[43]는 더 다양한 균형 진단을 할 수 있는 도구를 제공하여 bal.tab(), bal.plot()과 같은 함수를 통해 교란변수의 분포를 통계량이나 그림을 통해 비교할 수 있다.

Table 1에서 캘리퍼를 사용하지 않았을 때 교육받은 기간(Educ)과 고등학교 학위 유무(Nodegree)의 표준화 평균 차이가 0.1 이상임을 확인할 수 있다. 한편 나머지 변수들의 표준화 평균 차이는 0.1 미만이며 연속형인 교란변수들의 분산비율이 0.5와 2 사이에 있고 예후점수가 0.08로 일부변수에서는 균형이 상당히 맞춰진 것으로 나타났다. 또한 Figure 2의 Adjusted sample without caliper에서 두 교란변수의 실험군과 대조군의 분포가 유사함을 확인할 수 있다. 캘리퍼를 0.1로 설정하였을 때는 교육받은 기간과 고등학교 학위 유무 두 변수도 표준화 평균 차이가 0.1 미만으로 더 균형이 잘 이루어지며 Figure 2에서 캘리퍼를 작게 설정할수록 균형이 잘 이뤄지는 모습을 확인할 수 있다. 이처럼 캘리퍼를 설정하면 교란변수의 불균형을 개선할 수 있는데, 캘리퍼가 작아질수록 짝지어진 데이터의 크기가 작아지면서 연구하고자 하는 집단을 대표하는 결과로 해석하기 어렵다는 점을 유의해야 한다.

The result of balance diagnostic of 1:1 matched samples from nearest neighborhood method without caliper and with caliper 0.1

Figure 2.

Comparing the distributions of confounders between unmatched and matched samples for balance diagnosis. Distributions of the confounders with |standardized mean difference| ≥0.1. The leftmost plots of the panel are from unadjusted sample before matching, the middle plots are from adjusted sample after matching without caliper and the rightmost plots are from adjusted sample after matching with caliper 0.1.

매칭 후 데이터 분석

짝지어진 데이터에서 균형이 이루어진 것을 확인하였다면 올바르게 인과관계를 추론하기 위해 매칭 방법, 복원추출 여부, 결과변수의 자료형 등을 고려하여 적절한 분석 방법을 선택해야 한다. 짝지어진 데이터를 통해 ATE 또는 ATT를 구할 수 있는 대표적인 분석 방법으로 generalized estimating equation (GEE)을 고려해 볼 수 있다[28,44,45]. GEE은 같은 집단에 속한 개체들 간의 종속성을 허락하는 준모수적 추정 방법으로 샌드위치 분산 혹은 로버스트(robust) 분산을 통해 분포에 관계없이 올바른 분산 추정치를 제공한다. 짝지어진 자료를 분석하기 위해서는 처치와 결과를 각각 설명변수와 종속변수에 넣고 쌍마다 부여되는 번호를 집단변수로 처리하였을 때 나오는 처치의 회귀계수와 표준 오차(robust standard error)가 각각 인과효과의 추정치와 표준오차를 의미한다. 만약 일대다 매칭 또는 완전 매칭을 하였거나 대조군이 복원 추출되었다면 개체들의 가중치를 고려하여 분석해야 한다[28]. 일대다 매칭과 완전 매칭의 경우에는 각 쌍에 속한 같은 군의 개체수의 역수만큼 가중치를 부여하고 복원 추출의 경우에는 개체가 복원 추출된 횟수의 역수만큼 가중치를 부여하여 분석한다.

하지만 성향점수 모형이 잘못 설정되었다면 완전히 제거되지 않은 교란변수의 효과(residual confounding) 때문에 편향된 인과효과를 추론할 수 있다. 이에 대해 교란변수로 결과를 보정하여 편향을 줄이기 위한 방법들이 연구되었다[2,19,4648]. 가장 대표적인 방법으로 매칭 방법과 회귀 방법을 함께 사용하는 것이다[2,46]. 짝지어진 데이터를 갖고 GEE 방법으로 적합시킬 때 적절한 교란변수를 설명변수로 추가하면 매칭 이후에도 남아있던 교란변수의 효과가 상당히 제거될 수 있다. 하지만 이러한 보정 방법은 결과변수가 이분변수일 때 사용하면 전혀 다른 의미의 인과효과를 추론하게 되기 때문에 주의해서 사용해야 한다[49]. 조건부 교차비(odds ratio)의 기댓값은 연속형 변수와 달리 일반적으로 주변부 교차비와 일치하지 않으므로 로지스틱 모형에서 교란변수를 보정하였을 때 추론되는 인과효과는 보정하지 않았을 때와 다른 종류의 인과효과를 의미한다. 연속형 변수와 이분변수 모두에 쓸 수 있는 다른 보정 방법으로는 실험군과 대조군으로 데이터를 나누어 결과변수와 교란변수로 회귀모형을 만든 후 각 모형에 측정되지 않은 잠재변수를 채워 인과효과의 편향을 줄이는 방법이 연구되었다[19,50,51].

국가 지원 사업 데이터에서 각각 일대일 매칭과 완전 매칭을 통해 ATT와 ATE를 추정하였다. 이 두 추정치의 차이는 샘플에 의한 변동 외에 실험군에서의 프로그램 효과와 전체 집단에서 프로그램의 효과의 차이로 이해할 수 있다. 또한 남아있는 교란변수의 효과가 의심되는 상황을 고려하여 교란변수를 설명변수로 활용했을 때의 결과를 제 공하였고 이는 Table 2에서 확인할 수 있다.

Estimating treatment effects using matched sample from National Supported Work data

짝지어진 자료 분석에 대한 민감도 평가

관찰연구에서 인과효과를 데이터로 확인하기 어려운 가정들을 바탕으로 추론하기 때문에 가정의 일부가 성립하지 않는 경우에도 분석 결과가 여전히 강건한지 살펴야 한다. 민감도 분석은 결과의 강건함을 확인하기 위한 방법 중 하나로 가정의 위반 정도에 따라 분석의 결과가 얼마나 민감하게 변하는지 확인하는 방법을 의미한다. 인과성 추론에서는 조건부 교환 가능성 가정이 위반되었을 때 크게 편향된 결과를 줄 수 있기 때문에 이에 대한 민감도 분석 방법들이 주로 연구되고 있다[28,52].

매칭에서의 대표적인 민감도 분석 방법은 조건부 교환 가능성 가정이 위반된다면 서로 다른 교란변수를 가진 개체들끼리 짝지어져 인과효과를 편향되게 추정할 수 있다는 생각에서 출발하였다. 본래 매칭을 통해 짝지어진 개체들 간에 처치 받을 확률은 동일하지만 조건부 교환 가능성 가정을 위반하면 같은 쌍에 속했더라도 측정할 수 없는 교란변수에 의해 처치에 배정될 확률이 각기 달라질 수 있다. 민감도 분석에서는 측정할 수 없는 교란변수의 영향을 정량화하기 위해 처치에 배정될 확률의 교차비가 특정 범위 (1Γ)e(Cj)(1-e(Ck))e(Ck)(1-e(Cj))Γ,Γ1, j와 k는 같은 쌍에 속한 서로 다른 개체) 내에 있을 것을 가정한다. 그리고 나서 가정한 상황에서 얻을 수 있는 인과효과 추정치와 신뢰구간의 범위를 구해 가정 위반 정도에 따라 연구 결과가 얼마나 강건한지 판단한다. 인과효과의 추정치와 신뢰구간을 구하기 위해서 Wilcoxon의 부호 순위 통계량[52]이나 M-statistic [53,54] 등을 통해 계산되며 R에서는 ‘ sensitivitymv’ 패키지[55]와 ‘ rbounds’ 패키지[56]에서 계산 방법을 제공한다. 이외에 다른 민감도 분석 방법들은 측정되지 않은 교란변수가 처치와 결과 각각에 미치는 영향의 크기에 따라 갖게 될 인과효과의 크기를 계산하여 측정되지 않은 교란변수가 주는 영향이 얼마나 커야 분석결과가 강건하지 않을지를 평가한다[13,57].

Figure 3Table 2의 일대일 매칭의 결과에 대해 ‘sensitivitymv’ 패키지로 계산한 결과로 교차비가 커짐에 따라 인과효과의 추정치와 신뢰구간의 범위가 넓어지는 것을 확인할 수 있다. 특히 측정하지 못한 교란변수가 짝지어진 두 사람의 처치가 배정될 확률의 교차비가 1.1 정도로 영향을 주면 추정치의 범위에 0을 포함하게 되어 분석 결과가 뒤바뀌기 시작한다는 것을 의미한다.

Figure 3.

Sensitivity analysis plot using Rosenbaum's method with matched sample. Sensitivity analysis after the matched data analysis using National Supported Work data. A plot shows the boundaries of the point estimate and the ranges of 95% confidence interval (CI) for the average treatment effect on the treated according to odds ratio of treatment as-signment probability.

성향점수 매칭 방법론에 대한 주의점

앞에서는 성향점수 매칭 분석을 수행하는 과정에서의 주의할 점에 대해 설명하였는데 이 장에서는 성향점수 매칭 방법론 자체에 대한 주의사항에 대해 소개하고자 한다. King and Nielsen [4]은 성향점수 매칭을 다른 매칭 방법들과 비교를 통해 성향점수 매칭이 갖는 문제점을 제기하였다. 먼저 성향점수를 사용하지 않는 매칭 방법으로 완전히 동일한 교란변수를 갖는 개체들끼리 짝짓는 정확 매칭, 각 개체의 표준화된 교란변수 사이의 거리가 가까운 것끼리 짝짓는 마할라노비스 거리 매칭(Mahalanobis distance matching) [58], 그리고 교란변수가 연속형 또는 범주형 변수일 때 적당한 구간 또는 비슷한 범주로 묶어 정확 매칭을 하는 성긴 정확 매칭(coarsened exact matching) 등이 있다. 이러한 매칭 방법은 특정한 모형을 고려하지 않고 교란변수를 그대로 활용하여 짝짓는 방법인데 반해 성향점수 매칭은 모형을 기반으로 하여 성향점수를 추정한 뒤 유사한 값을 갖는 관측치끼리 짝짓기 때문에 모에 의존하는 방법이다. 연구자는 어떤 모형이 참인지 모르기 때문에 여러 성향점수 모형을 만들고 가장 적절해 보이는 것 하나를 선택해야 하는데, 이때 필연적으로 연구자의 주관이 개입된다. 특히 King and Nielsen [4]은 연구자는 가설을 기각하기 위해 인과효과가 크게 추론되는 모형을 선택할 수 있으며 그러한 경향이 많아질수록 실제 효과와의 편향이 커지게 된다는 부분을 지적하였다.

위에서 언급된 성향점수 매칭 이외의 세 가지 방법들은 유사한 교란변수를 가진 개체들끼리 짝짓기 때문에 근본적으로 블록 무작위배정(blocked randomization)을 구현하려는 방법으로 볼 수 있다. 반면에 성향점수 매칭은 성향점수가 동일한 개체들끼리 짝지었을 때 교란변수가 서로 다른 개체들끼리도 짝지어질 수 있어 완전 무작위배정을 구현하려는 것으로 볼 수 있다[4]. 일반적으로 완전 무작위 배정은 블록 무작위배정 방법보다 교란변수의 불균형을 줄이는 데 비효율적인 방법이기에 성향점수 매칭을 다른 매칭 방법에 비해 무조건적으로 우선시 할 필요는 없게 된다. 실제로 King and Nielsen [4]은 수치 연구를 통해 매칭을 허용하는 기준을 강화시켜 나갈 때 마할라노비스 거리 매칭에서는 매칭을 허용하는 기준을 강화시켜 나갈수록, 즉 점점 버려지는 데이터가 증가하여도 인과효과 추정치의 분산과 편향이 꾸준히 감소하다 일정 수준에 도달하는 데 반해 성향점수 매칭에서는 유사한 조건하에서 인과효과 추정치 분산과 편향이 초반에 감소하다 일정수준에 도달하는 것이 아니라 다시 빠르게 증가하는 결과를 확인하였다. 이 결과가 시사하는 바는 최근접 이웃 매칭에서 매칭의 품질을 높이기 위해 캘리퍼를 어느 수준 이하로 낮추었을 때, 인과효과 추정치는 분산과 편향이 클 수 있다는 위험성을 보여주는 것으로 성향점수 기반 매칭의 품질이 우수하다고 해서 그것을 인과효과 추정치가 올바르게 구해졌다는 것으로 함부로 해석하면 안된다는 것을 보여준다.

토 론

지금까지의 성향점수와 관련된 방법들은 대부분 이분 변수의 처치에 대해 다루고 있다. 기존의 성향점수를 다중 처치 또는 연속형 처치로 확장한 일반화 성향점수 또한 연구된 바 있으며 기계학습 모델을 이용하여 일반화 성향점수를 추정하는 방법도 소프트웨어로 개발되어 있다[34,59]. 하지만 일반화 성향점수를 유한한 수의 자료에서 추정하고자 할 때 처치의 수가 많아지면 특정 교란변수의 조합에서 어떤 처치를 받을 확률이 0이 되어 양수성 가정을 위반할 가능성이 높다는 점에서 그 사용이 제한적이다[22]. 또한 다중 처치군에 대한 매칭 알고리즘 및 짝지어진 자료를 분석에 대해서는 아직 연구가 진행 중이다[6062]. 최근 연구에서 외상 치료 전문 센터의 수준과 치명률의 인과관계를 밝히기 위해 세 가지 처치군에 대한 매칭 방법을 사용하였다 [62]. 이 연구에서는 세 처치군에 대한 1:1:1 최적합 매칭 방법을 제안하였는데, 같은 짝에서 개체들 간의 거리의 합이 최소가 되도록 두 군을 먼저 짝짓고 남은 한 군을 짝 지은 이후, 한 개의 군을 돌아가며 선택하여 거리의 합이 최소화되도록 남은 두 군의 짝지어진 그룹과 다시 짝짓는 행위를 반복하는 방법을 사용하였다. 하지만 이러한 매칭 알고리즘은 처치군이 더욱 늘어나면 연산량이 기하급수적으로 증가하여 확장이 어려우며, 따라서 처치가 더 많은 경우에 대한 매칭 알고리즘, 분석 방법, 소프트웨어 개발이 필요하다.

또 다른 이슈로는 성향점수 모형이 잘못 설정되어 짝지어진 데이터의 균형이 이뤄지지 않았을 때 얼마나 잘못된 결과가 나올 수 있는지 여러 연구에서 언급된 바 있다[4]. 이에 대해 여러 비모수적인 모형들을 통해 성향점수 모형을 개선시키려는 여러 시도가 있었으나[34,35], 이러한 방법들은 애초에 최적의 균형을 찾기 위한 모형들이 아니라 높은 분류성능을 갖도록 설계되었기 때문에 균형을 맞추는 데 한계가 있을 수 있다. 이를 극복하기 위해 제안된 한 가지 방법으로 성향점수의 모형이 분포의 적률(moment)을 일치시키는 데 초점을 맞춘 공변량 균형 성향점수(covariate balancing propensity score, CBPS)가 제시되었다[63]. Imai and Ratkovic [63]은 기존의 성향점수와 성능 비교를 위해 참인 성향점수 모형을 교란변수의 고차항과 교호작용을 포함한 것으로 가정하고 교란변수의 일차항만 주어졌을 때를 가정한 수치연구로 하였다. 그 결과 성향점수를 활용하는 역 확률 가중치(inverse probability weighting)와 이중 강건(doubly robust) 방법에서 기존의 성향점수와 CBPS를 각각 사용하여 추정량을 구했을 때, 모형이 복잡할수록 성향점수는 평균 제곱 오차가 증가하는 반면 CBPS는 일정하게 유지되었다. 또한 매칭 방법에서도 표준적인 방법들에 비해 편향이 크게 줄어든 인과효과 추정치를 제공하는 것으로 보고하였다. 이처럼 교란변수의 균형에 초점을 맞춘 CBPS에 대해 최근 많이 연구되고 있으며 특히 연속형 처치로의 확장도 고려되고 있다[64]. 그러나 CBPS가 갖는 이론적인 강점이 실제 데이터 분석에서 충분히 나타나는지에 대해서는 아직 세밀한 검토가 이루어지지 않았다.

마지막으로는 민감도 분석에서 측정되지 않은 교란변수가 존재할 때 처치 받을 확률의 교차비의 범위를 어떻게 설정해야 하는가의 문제이다. Γ는 측정되지 않은 교란변수가 처치의 배정에 주는 영향과 관련이 있는 값으로 측정되지 않은 교란변수에 대한 정보가 없기 때문에 데이터로 추론하기 매우 어렵다. 다른 민감도 분석에서는 측정되지 않은 교란변수가 미칠 수 있는 영향의 기준을 설정하기 위해 측정된 교란변수를 하나씩 모형에서 제외했을 때와 포함하였을 때를 비교하여 모형이 변화하는 정도를 기준으로 삼는 방법을 제시하였다[13,65]. 그러나 이러한 방법은 고의로 제외한 변수에 의한 효과와 측정되지 않은 교란변수의 효과를 정확하게 분리해내는 방법이 아니므로 Γ 값을 설정하려는 구체적인 기준에 대해서는 추가적인 연구가 필요하다.

References

1. . Fisher R. Statistical methods for research workers 13th ed.th ed. London, UK: Oliver and Loyd, Ltd.; 1925. p. 99–101.
2. . Cochran WG, Rubin DB. Controlling bias in observational studies: a review. Sankhyā: Indian J Stat, Series A 1973;35(4):417–446.
3. . Greenwood E. Experimental sociology New York: Columbia University Press; 1945.
4. . King G, Nielsen R. Why propensity scores should not be used for matching. Political Anal 2019;27(4):435–454. 10.1017/pan.2019.11.
5. . Chapin FS. Experimental designs in sociological research New York: Harper & Row; 1947.
6. . Rubin DB. Matching to remove bias in observational studies. Biometrics 1973;29(1):159–183. 10.2307/2529684.
7. . Rosenbaum PR, Rubin DB. The central role of the propensity score in observational studies for causal effects. Biometrika 1983;70(1):41–55.
8. . Rosenbaum PR, Rubin DB. Constructing a control group using multivariate matched sampling methods that incorporate the propensity score. Am Stat 1985;39(1):33–38.
9. . Rubin DB, Thomas N. Characterizing the effect of matching using lin-ear propensity score methods with normal distributions. Biometrika 1992;79(4):797–809. 10.2307/2337235.
10. . Rubin DB, Thomas N. Matching using estimated propensity scores: relating theory to practice. Biometrics 1996;52(1):249–264. 10.2307/2533160.
11. . Dehejia RH, Wahba S. Causal effects in nonexperimental studies: re-evaluating the evaluation of training programs. J Am Stat Assoc 1999;94(448):1053–1062. 10.2307/2669919.
12. . Smith JA, Todd PE. Does matching overcome LaLonde's critique of nonexperimental estimators? J Econometrics 2005;125(1-2):305–353. 10.1016/j.jeconom.2004.04.011.
13. . Imbens GW. Sensitivity to exogeneity assumptions in program evaluation. Am Econ Rev 2003;93(2):126–132.
14. . Pearl J. Causal inference in statistics: an overview. Stat Surv 2009;3:96–146. 10.1214/09-SS057.
15. . Imbens GW. Nonparametric estimation of average treatment effects under exogeneity: a review. Rev Econ Stat 2004;86(1):4–29.
16. . LaLonde RJ. Evaluating the econometric evaluations of training programs with experimental data. Am Econ Rev 1986;76(4):604–620.
17. . Rubin DB. Estimating causal effects of treatments in randomized and nonrandomized studies. J Educ Psychol 1974;66(5):688.
18. . Rubin DB. Randomization analysis of experimental data: The Fisher randomization test comment. J Am Stat Assoc 1980;75(371):591–593. 10.2307/2287653.
19. . Abadie A, Imbens GW. Bias-corrected matching estimators for average treatment effects. J Bus Econ Stat 2011;29(1):1–11. 10.1198/jbes. 2009.07333.
20. . James G, Witten D, Hastie T, Tibshirani R. An introduction to statisti-cal learning (Vol. 112) New York: Springer; 2013.
21. . Spirtes P, Glymour CN, Scheines R, Heckerman D. Causation, prediction, and search 2nd ed.th ed. Cambridge: MIT Press; 2000.
22. . Hernán M, Robins JM. Causal inference: What If Boca Raton: CRC Press; 2020.
23. . Pearl J. Causal diagrams for empirical research. Biometrika 1995;82(4):669–688. 10.2307/2337329.
24. . Pearl J. Models, reasoning and inference Cambridge: Cambridge University Press; 2000. p. 2. 19.
25. . Textor J, Hardt J, Knüppel S. DAGitty: a graphical tool for analyzing causal diagrams. Epidemiol 2011;22(5):745. 10.1097/EDE.0b-013e318225c2be.
26. . Hansen BB. Full matching in an observational study of coaching for the SAT. J Am Stat Assoc 2004;99(467):609–618. 10.1198/01621-4504000000647.
27. . Rosenbaum PR. A characterization of optimal designs for observational studies. J R Statist Soc A 1991;53(3):597–610.
28. . Stuart EA. Matching methods for causal inference: a review and a look forward. Stat Sci 2010;25(1):1–21. 10.1214/09-STS313.
29. . Austin PC. A comparison of 12 algorithms for matching on the propensity score. Stat Med 2014;33(6):1057–1069. 10.1002/sim.6004.
30. . Greifer N, Stuart EA. Choosing the estimand when matching or weighting in observational studies arXiv:2106.10577. 2021. 10.48550/arXiv.2106.10577.
31. . Thomas L, Li F, Pencina M. Using propensity score methods to create target populations in observational clinical research. JAMA 2020;323(5):466–467. 10.1001/jama.2019.21558.
32. . Austin PC. Balance diagnostics for comparing the distribution of baseline covariates between treatment groups in propensity-score matched samples. Stat Med 2009;28(25):3083–3107. 10.1002/sim.3697.
33. . Drake C. Effects of misspecification of the propensity score on estimators of treatment effect. Biometrics 1993;49(4):1231–1236. 10.2307/2532266.
34. . McCaffrey DF, Ridgeway G, Morral AR. Propensity score estimation with boosted regression for evaluating causal effects in observational studies. Psychol Methods 2004;9(4):403–425. 10.1037/1082-989X.9.4.403.
35. . Lee BK, Lessler J, Stuart EA. Improving propensity score weighting using machine learning. Stat Med 2010;29(3):337–346. 10.1002/sim.3782.
36. . Collier ZK, Leite WL. A tutorial on artificial neural networks in propensity score analysis. J Exp Educ 2020;:1–18. 10.1080/00220973. 2020.1854158.
37. . Singh C, Balakrishnan G, Perona P. Matched sample selection with GANs for mitigating attribute confounding arXiv:2103.13455. 2021. 10.48550/arXiv.2103.13455.
38. . Normand ST, Landrum MB, Guadagnoli E, Ayanian JZ, Ryan TJ, Cleary PD, et al. Validating recommendations for coronary angiography fol-lowing acute myocardial infarction in the elderly: a matched analysis using propensity scores. J Clin Epidemiol 2001;54(4):387–398. 10.1016/s0895-4356(00)00321-8.
39. . Stuart EA, Lee BK, Leacy FP. Prognostic score–based balance measures can be a useful diagnostic for propensity score methods in comparative effectiveness research. J Clin Epidemiol 2013;66(8 Supple):S84–S90. 10.1016/j.jclinepi.2013.01.013.
40. . Ho DE, Imai K, King G, Stuart EA. Matching as nonparametric pre-processing for reducing model dependence in parametric causal inference. Political Anal 2007;15(3):199–236. 10.1093/pan/mpl013.
41. . Imai K, King G, Stuart EA. Misunderstandings between experimental-ists and observationalists about causal inference. J R Statist Soc A 2008;171(2):481–502.
42. . Ho D, Imai K, King G, Stuart E, Whitworth A. Package ‘ MatchIt'. Version 2018.
43. . Greifer N. Covariate balance tables and plots: a guide to the cobalt package. Accessed March 2020;10:2020.
44. . Austin PC, Stuart EA. Estimating the effect of treatment on binary outcomes using full matching on the propensity score. Stat Methods Med Res 2017;26(6):2505–2525. 10.1177/0962280215601134.
45. . Stuart EA, Green KM. Using full matching to estimate causal effects in nonexperimental studies: examining the relationship between adoles-cent marijuana use and adult outcomes. Dev Psychol 2008;44(2):395–406. 10.1037/0012-1649.44.2.395.
46. . Rubin DB. The use of matched sampling and regression adjustment to remove bias in observational studies. Biometrics 1973;29(1):185–203. 10.2307/2529685.
47. . Carpenter RG. Matching when covariables are normally distributed. Biometrika 1977;64(2):299–307. 10.2307/2335697.
48. . Robins JM, Rotnitzky A. Semiparametric efficiency in multivariate regression models with missing data. J Am Stat Assoc 1995;90(429):122–129. 10.2307/2291135.
49. . Austin PC. The performance of different propensity score methods for estimating marginal odds ratios. Stat Med 2007;26(16):3078–3094. 10.1002/sim.2781.
50. . Austin PC. Double propensity-score adjustment: a solution to design bias or bias due to incomplete matching. Stat Methods Med Res 2017;26(1):201–222. 10.1177/0962280214543508.
51. . Nguyen TL, Collins GS, Spence J, Daurès JP, Devereaux PJ, Landais P, et al. Double-adjustment in propensity score matching analysis: choosing a threshold for considering residual imbalance. BMC Med Res Methodol 2017;17(1):78. 10.1186/s12874-017-0338-0.
52. . Rosenbaum PR. Overt bias in observational studies. Observational studies New York: Springer; 2002. p. 71–104.
53. . Rosenbaum PR. Sensitivity analysis for m-estimates, tests, and confidence intervals in matched observational studies. Biometrics 2007;63(2):456–464. 10.1111/j.1541-0420.2006.00717.x.
54. . Rosenbaum PR. Weighted M-statistics with superior design sensitivity in matched observational studies with multiple controls. J Am Stat Assoc 2014;109(507):1145–1158.
55. . Rosenbaum PR. Two R packages for sensitivity analysis in observational studies. Observational Stud 2015;1(2):1–17.
56. . Keele LJ. Package rbounds: An R package 2014.
57. . VanderWeele TJ, Ding P. Sensitivity analysis in observational research: introducing the E-value. Ann Intern Med 2017;167(4):268–274. 10.7326/M16-2607.
58. . Rubin DB. Bias reduction using Mahalanobis metric matching. ETS Research Bulletin Series 1978;1978(2):i–10. 10.1002/j.2333-8504.1978.tb01164.x.
59. . Hirano K, Imbens GW. The propensity score with continuous treatments Hoboken: John Wiley & Sons, Ltd; 2004. p. 73–84.
60. . Wesley Routon P, Walker JK. College internships, tenure gaps, and stu-dent outcomes: a multiple-treatment matching approach. Educ Econ 2019;27(4):383–400. 10.1080/09645292.2019.1598336.
61. . Wellington M, Mortlock M, Moepeng P. Evaluating food transfers in Botswana using multiple matching methods. Development Southern Africa 2021 10.1080/0376835X.2021.1984874.
62. . Nattino G, Lu B, Shi J, Lemeshow S, Xiang H. Triplet matching for estimating causal effects with three treatment arms: a comparative study of mortality by trauma center level. J Am Stat Assoc 2021;116(533):44–53. 10.1080/01621459.2020.1737078.
63. . Imai K, Ratkovic M. Covariate balancing propensity score. J R Statist Soc B 2014;76(1):243–263.
64. . Fong C, Hazlett C, Imai K. Covariate balancing propensity score for a continuous treatment: application to the efficacy of political advertise-ments. Ann Appl Stat 2018;12(1):156–177. 10.1214/17-AOAS1101.
65. . Veitch V, Zaveri A. Sense and sensitivity analysis: simple post-hoc anal-ysis of bias due to unobserved confounding arXiv:2003.01747. 2020. 10.48550/arXiv.2003.01747.

Article information Continued

Figure 1.

Three directed acyclic graphs demonstrating various relation-ships of variables. Let T, Y, C, S, U denote treatment, outcome, confounder, collider, unmeasured confounder, respectively. The first two DAGs repre-sent no causal effect from treatment to outcome with (A) an open back door path (T← C→ Y) by confounder C and (B) no open backdoor path. (C) A DAG represents the causal relationship in the National Supported Work data. DAG, directed acyclic graph.

Table 1.

The result of balance diagnostic of 1:1 matched samples from nearest neighborhood method without caliper and with caliper 0.1

Variables NNM (prognostic score=0.08)
NNM with caliper (prognostic score=0.06)
SMD Variance ratio SMD Variance ratio
Age -0.01 0.96 0.01 1.08
Educ 0.10 1.37 0.05 1.03
Black -0.01 - -0.03 -
Hispanic -0.02 - -0.03 -
Married -0.03 - -0.02 -
Nodegree -0.14 - -0.04 -
Re74 0.09 1.18 0.01 1.02
Re75 0.05 0.93 -0.05 0.52

SMD, standardized mean difference; NNM, nearest neighborhood method. Matching results for the National Supported Work data using 1:1 nearest neighborhood method without caliper and with caliper 0.1. The balance for the 8 confounders (age [Age], year of schooling [Educ], blacks [Black], hispanics [Hispanic], marital status [Married], high school di-ploma [nodegree→ Nodegree], real earnings in 1974 [Re74], and real earnings in 1975 [Re75]) is examined.

Figure 2.

Comparing the distributions of confounders between unmatched and matched samples for balance diagnosis. Distributions of the confounders with |standardized mean difference| ≥0.1. The leftmost plots of the panel are from unadjusted sample before matching, the middle plots are from adjusted sample after matching without caliper and the rightmost plots are from adjusted sample after matching with caliper 0.1.

Table 2.

Estimating treatment effects using matched sample from National Supported Work data

Matching methods (estimand) Treatment effect (standard error) Treatment effect after confounder adjustment (standard error)
1:1 matching (ATT) 2,017.87 (675.84) 1,874.85 (658.33)
Full matching (ATE) 1,660.37 (641.06) 1,519.64 (615.13)

ATT, average treatment effect on the treated; ATE, average treatment effect.

Estimated treatment effects from matched sample after 1:1 matching (first row) and full matching (second row) methods.

Figure 3.

Sensitivity analysis plot using Rosenbaum's method with matched sample. Sensitivity analysis after the matched data analysis using National Supported Work data. A plot shows the boundaries of the point estimate and the ranges of 95% confidence interval (CI) for the average treatment effect on the treated according to odds ratio of treatment as-signment probability.