12 IP Weighting and Marginal Structural Models

이 챕터에서는 관측 데이터로부터 인과 효과를 추정하기 위해 어떻게 역확률 가중치(IP weighting, IP 가중치)를 사용하는지를 다룬다. 정확히는 추가적인 가정들 하에서 IP 가중치를 사용하여 여러 공변량과 비인산형 처치가 가지는 (변수가 많아서 나타나는) 고차원성의 문제를 완화하는 방법에 대해 살펴본다.

먼저 여기서 사용할 데이터 NHEFS를 다운받아보자:

nhefs <- read_csv("https://www.hsph.harvard.edu/miguel-hernan/wp-content/uploads/sites/1268/2019/03/nhefs.csv")
glimpse(nhefs)

Rows: 1,629
Columns: 64
$ seqn              <dbl> 233, 235, 244, 245, 252, 25…
$ qsmk              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ death             <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 1, …
$ yrdth             <dbl> NA, NA, NA, 85, NA, NA, NA,…
$ modth             <dbl> NA, NA, NA, 2, NA, NA, NA, …
$ dadth             <dbl> NA, NA, NA, 14, NA, NA, NA,…
$ sbp               <dbl> 175, 123, 115, 148, 118, 14…
$ dbp               <dbl> 96, 80, 75, 78, 77, 83, 69,…
$ sex               <dbl> 0, 0, 1, 0, 0, 1, 1, 1, 0, …
$ age               <dbl> 42, 36, 56, 68, 40, 43, 56,…
$ race              <dbl> 1, 0, 1, 1, 0, 1, 0, 0, 0, …
$ income            <dbl> 19, 18, 15, 15, 18, 11, 19,…
$ marital           <dbl> 2, 2, 3, 3, 2, 4, 2, 2, 2, …
$ school            <dbl> 7, 9, 11, 5, 11, 9, 12, 12,…
$ education         <dbl> 1, 2, 2, 1, 2, 2, 3, 3, 2, …
$ ht                <dbl> 174.2, 159.4, 168.5, 170.2,…
$ wt71              <dbl> 79.04, 58.63, 56.81, 59.42,…
$ wt82              <dbl> 68.95, 61.23, 66.22, 64.41,…
$ wt82_71           <dbl> -10.0940, 2.6050, 9.4145, 4…
$ birthplace        <dbl> 47, 42, 51, 37, 42, 34, NA,…
$ smokeintensity    <dbl> 30, 20, 20, 3, 20, 10, 20, …
$ smkintensity82_71 <dbl> -10, -10, -14, 4, 0, 10, 0,…
$ smokeyrs          <dbl> 29, 24, 26, 53, 19, 21, 39,…
$ asthma            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ bronch            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ tb                <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ hf                <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ hbp               <dbl> 1, 0, 0, 1, 0, 0, 0, 0, 0, …
$ pepticulcer       <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, …
$ colitis           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ hepatitis         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ chroniccough      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ hayfever          <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, …
$ diabetes          <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, …
$ polio             <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ tumor             <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, …
$ nervousbreak      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ alcoholpy         <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ alcoholfreq       <dbl> 1, 0, 3, 2, 2, 3, 1, 0, 1, …
$ alcoholtype       <dbl> 3, 1, 4, 3, 1, 2, 3, 2, 3, …
$ alcoholhowmuch    <dbl> 7, 4, NA, 4, 2, 1, 4, 1, 2,…
$ pica              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ headache          <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, …
$ otherpain         <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 1, …
$ weakheart         <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, …
$ allergies         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ nerves            <dbl> 0, 0, 1, 0, 0, 0, 1, 1, 0, …
$ lackpep           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ hbpmed            <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, …
$ boweltrouble      <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, …
$ wtloss            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ infection         <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, …
$ active            <dbl> 0, 0, 0, 1, 1, 1, 0, 0, 2, …
$ exercise          <dbl> 2, 0, 2, 2, 1, 1, 1, 2, 2, …
$ birthcontrol      <dbl> 2, 2, 0, 2, 2, 0, 0, 0, 2, …
$ pregnancies       <dbl> NA, NA, 2, NA, NA, 1, 1, 2,…
$ cholesterol       <dbl> 197, 301, 157, 174, 216, 21…
$ hightax82         <dbl> 0, 0, 0, 0, 0, 1, NA, NA, 0…
$ price71           <dbl> 2.184, 2.347, 1.570, 1.507,…
$ price82           <dbl> 1.740, 1.797, 1.513, 1.452,…
$ tax71             <dbl> 1.1023, 1.3650, 0.5513, 0.5…
$ tax82             <dbl> 0.4620, 0.5719, 0.2310, 0.2…
$ price71_82        <dbl> 0.44379, 0.54932, 0.05620, …
$ tax71_82          <dbl> 0.6404, 0.7930, 0.3203, 0.3…

12.1 The causal question

금연이 체중 증가에 미치는 평균 인과효과를 추정하는 것이 목적이라고 하자. 이를 위해서 NHEFS 데이터에서 25-74세 사이의 흡연자 1,566명의 데이터를 사용한다. 아래 표에서 나타나듯, 금연자와 흡연자는 성별, 인종, 교육수준, 몸무게, 흡연 정도 등의 공변량의 분포에 있어서도 차이를 보이고 있다. 즉, 이러한 변수들이 교란변수라면 분석 이전에 이들에 대한 조정이 필요하게 된다.

# provisionally ignore subjects with missing values for weight in 1982
nhefs <- nhefs |> drop_na(wt82)

nhefs |> drop_na(qsmk) |> 
  dplyr::filter(age > 24 & age < 75) |> 
  group_by(qsmk) |>
  summarize(
    `Age, years` = mean(age, na.rm = T), 
    `Men, %` = mean(sex, na.rm = T),
    `White, %` = mean(race, na.rm = T),
    `University, %` = mean(education, na.rm = T),
    `Weight, kg` = mean(wt71, na.rm = T),
    `Cigarettes/day` = mean(smokeintensity, na.rm = T),
    `Years smoking` = mean(smokeyrs, na.rm = T),
    `Little exercise` = mean(exercise, na.rm = T),
    `Inactive life, %` = mean(active, na.rm = T),
  ) |> pivot_longer(cols = -qsmk) |> 
   pivot_wider(names_from = qsmk) |> 
  rename("Mean baseline characteristics" = name) |> 
  mutate_if(is.numeric, round, 2) |> 
  gt::gt() |> 
  gt::tab_spanner(label = "A", columns = c("0", "1"))

Mean baseline characteristics	A
Mean baseline characteristics	0	1
Age, years	42.79	46.17
Men, %	0.53	0.45
White, %	0.15	0.09
University, %	2.69	2.79
Weight, kg	70.30	72.35
Cigarettes/day	21.19	18.60
Years smoking	24.09	26.03
Little exercise	1.18	1.25
Inactive life, %	0.63	0.69

12.2 Estimating IP weights via modeling

IP 가중치는 공변량 $L$ 로부터 처치 $A$ 로의 화살표가 제거된 의사 모집단을 만드는 방법이다. 보다 정확하게는 그 의사 모집단은 다음과 같은 두 가지 특징을 가지게 된다.

$A$ 와 $L$ 은 통계적으로 독립적이다.
의사 모집단의 평균 $E_{p s} [Y | A = a]$ 은 실제 모집단에서의 표준화된 평균 $\sum_{l} E [Y | A = a, L = l] Pr [L = l]$ 과 같다.

$Y, A, L$ 에 대한 조건부 교환 가능성이 실제 모집단에서 유지된다면, 이러한 속성은 (i) $Y^{a}$ 의 평균이 두 모집단에서 동일하고, (ii) 무조건 교환 가능성(즉, 교란 없음)이 의사 모집단에서 유지됨을 의미한다. (iii) 실제 모집단에서 반사실적 평균 $E [Y^{a}]$ 는 유사 모집단의 평균 $E_{p s} [Y | A = a]$ 와 같으며, (iv) 유사 모집단에서 연관성은 인과 관계이다.

비공식적으로 의사 모집단은 각 개인이 실제로 받은 처치 수준을 받을 조건부 확률의 역으로 가중치를 부여하여 만들어진다.
$L$ 로 정의된 수백만 개의 계층 각각에서 $Pr [A = 1 | L]$ 의 모수 추정치를 구하기 위해 9개의 교란 변수를 모두 공변량으로 포함하여 금연 확률에 대한 로지스틱 회귀 모델을 적합한다.

fit <-
  glm(
    qsmk ~ sex + race + age + I(age ^ 2) +
      as.factor(education) + smokeintensity +
      I(smokeintensity ^ 2) + smokeyrs + I(smokeyrs ^ 2) +
      as.factor(exercise) + as.factor(active) + wt71 + I(wt71 ^ 2),
    family = binomial(),
    data = nhefs)
texreg::screenreg(fit)


==================================
                       Model 1    
----------------------------------
(Intercept)              -2.24    
                         (1.38)   
sex                      -0.53 ***
                         (0.15)   
race                     -0.84 ***
                         (0.21)   
age                       0.12 *  
                         (0.05)   
age^2                    -0.00    
                         (0.00)   
as.factor(education)2    -0.03    
                         (0.20)   
as.factor(education)3     0.09    
                         (0.18)   
as.factor(education)4     0.06    
                         (0.27)   
as.factor(education)5     0.48 *  
                         (0.23)   
smokeintensity           -0.08 ***
                         (0.02)   
smokeintensity^2          0.00 ***
                         (0.00)   
smokeyrs                 -0.07 ** 
                         (0.03)   
smokeyrs^2                0.00    
                         (0.00)   
as.factor(exercise)1      0.35 *  
                         (0.18)   
as.factor(exercise)2      0.40 *  
                         (0.19)   
as.factor(active)1        0.03    
                         (0.13)   
as.factor(active)2        0.18    
                         (0.21)   
wt71                     -0.02    
                         (0.03)   
wt71^2                    0.00    
                         (0.00)   
----------------------------------
AIC                    1714.90    
BIC                    1816.67    
Log Likelihood         -838.45    
Deviance               1676.90    
Num. obs.              1566       
==================================
*** p < 0.001; ** p < 0.01; * p < 0.05

로지스틱 회귀모델 적합한 이후에 추정된 IP 가중치를 이용하여 만들어진 의사 모집단에서 ${\hat{E}}_{p s} [Y | A = 1] - {\hat{E}}_{p s} [Y | A = 0]$ 차이를 컴퓨팅한다.

p.qsmk.obs <- 
  ifelse(nhefs$qsmk == 0,
         1 - predict(fit, type = "response"),
         predict(fit, type = "response"))

nhefs |> mutate(w = 1 / p.qsmk.obs) -> nhefs
summary(nhefs$w); sd(nhefs$w)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.05    1.23    1.37    2.00    1.99   16.70

[1] 1.475

흡연자일 경우에는 1에서 예측값을 빼서 금연이 아닌 흡연자에 대한 예측 확률을, 반면 금연자일 경우에는 금연자일 예측확률을 계산한 다음에 이 역수를 가중치 w로 계산한다. 그리고 이 가중치를 이용하여 금연 여부(처치)와 종속변수 간의 관계를 추정한다. 이미 이 가중치는 다른 공변량의 내용을 반영하고 있기 때문에 여기에서는 기타 공변량은 포함하지 않는다. 결과, 금연은 평균적으로 ${\hat{θ}}_{1} = 3.4 kg$ 을 증가시킨다는 것을 확인할 수 있다.

여기서 사용한 {geepack}의 geeglm() 함수는 강건분산추정량을 제공하는 일반화 추정 방정식(Generalized Estimating Equations, GEE) 모델이다.
만약 $Pr [A = 1 | L]$ 이라는 모델이 잘못 특정되었다면, $θ_{0}, θ_{1}$ 에 대한 추정치도 편향될 것이다.

# install.packages("geepack") # install package if required
library("geepack")
msm.w <- geeglm(
  wt82_71 ~ qsmk,
  data = nhefs,
  weights = w,
  id = seqn,
  corstr = "independence"
)
summary(msm.w)


Call:
geeglm(formula = wt82_71 ~ qsmk, data = nhefs, weights = w, id = seqn, 
    corstr = "independence")

 Coefficients:
            Estimate Std.err Wald           Pr(>|W|)
(Intercept)    1.780   0.225 62.7 0.0000000000000023
qsmk           3.441   0.525 42.9 0.0000000000586079
               
(Intercept) ***
qsmk        ***
---
Signif. codes:  
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation structure = independence 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)     65.1    4.22
Number of clusters:   1566  Maximum cluster size: 1

12.3 Stabilized IP weights

교란 조정의 핵심 요건은 의사 모집단에서 처치 $A$ 의 확률이 교란 변수 $L$ 에 의존하지 않는다는 것이다. 따라서 의사 모집단의 모두가 처치를 배정받을 확률을 동일하게 $p$ 로 놓는다면 이 조건을 만족시킬 수 있다. 따라서 이전 섹션에서 역수를 취해서 부여한 가중치와는 다르게 처치를 받을 확률과 받지 않을 확률이 $L$ 과 무관하게 0.5라고 할 때도 마찬가지로 IP 가중치 방법을 사용할 수 있게 된다. 이렇게 역수가 아니라 특정한 확률을 바탕으로 취한 IP 가중치 방법을 안정화 IP 가중치(stabilized IP weights)라고 한다.

만약 $W^{A}$ 와 $S W^{A}$ 를 사용한 추정치가 동일하다면, 왜 $S W^{A}$ 를 사용하는 걸까?

더 좁은 신뢰구간을 제공하기 때문이다.
단, 이러한 통계적인 우위는 오직 모델이 포화(너무 복잡)하지 않은 경우에만 성립된다.

12.4 Marginal structural models

처치 $a$ 하에서의 결과의 평균에 대한 선형 모델이 다음과 같다고 하자:

$E [Y^{a}] = β_{0} + β_{1} a$ 이 모델의 결과변수에 대한 반사실적 결과는 일반적으로 관측 불가능하다. 따라서 모델은 실제 세계에 관한 연구의 데이터로 적합될 수 없다. 반사실적 결과의 한계평균(marginal mean)에 대한 모델들은 한계구조평균모델(marginal structural mean models)이라고 한다.

구조평균모델에서의 처치에 대한 모수는 평균 인과효과에 해당한다. 즉, $β_{1}$ 은 $E [Y^{a = 1}] - E [Y^{a = 0}]$ 이 된다.

우리는 IP 가중치를 이용해서 의사 모집단을 만들고, 모델 $E [Y | A] = θ_{0} + θ_{1} A$ 를 그 의사 모집단에 IP 가중치를 취한 최소제곱을 사용하여 적합할 수 있다. IP 가중치의 가정에 따라 의사 모집단에서의 관계는 곧 인과관계를 의미한다.

금연 $A$ 가 이산형 처치라는 점에서 한계구조모델 $E [Y^{a}] = β_{0} + β_{1} a$ 은 포화된 모델이다.

모델은 방정식의 양변에 2개의 알 수 없는 요인들을 가지고 있다: 좌변에는 $E [Y^{a = 1}], E [Y^{a = 0}]$ 를, 우변에는 $β_{0}, β_{1}$ 을 갖고 있는 것이다.
따라서 의사 모집단에서 컴퓨팅된 표본 평균은 인과효과를 추정하기에는 충분하다.

문제는 처치가 종종 다항이거나 연속형일 경우이다. 이 경우에 우리는 아무 $a$ 와 $a^{'}$ 에 대한 값에서 $E [Y^{a}] - E [Y^{a^{'}}]$ 를 추정하고자 한다. 이 경우에 처치의 값이 무수히 많을 수 있기 때문에 모든 모수를 추정한다는 것은 불가능하다. 이 경우, 처치 $A$ 가 결과 $Y$ 의 평균에 미치는 효과에 대한 무수히 많은 응답의 곡선을 특정하기 위해 비포화구조모델을 고려해볼 수 있다.

$E [Y^{a}] = β_{0} + β_{1} a β_{2} a^{2}$

우리가 흡연 정도의 증가가 가지는 인과 효과를 추정하고 싶다고 하자. $E [Y^{a = 20}] = β_{0} + 20 β_{1} + 400 β_{2}$ 가 되므로, $β_{1}, β_{2}$ 를 추정할 필요가 있다. 이를 위해서 $L$ 에 의해 교란이 나타나지 ㅇ낳는 의사 모집단을 만들기 위해 IP 가중치 $S W^{A}$ 를 추정할 필요가 있게 되는 것이다. 그리고 모델 $E [Y | A] = θ_{0} + θ_{1} A + θ_{2} A^{2}$ 를 의사 모집단에 적합한다.

12.5 Effect modification and marginal structural models

한계구조모델은 우리가 알고자 하는 모수가 모집단의 평균 인과효과일 경우에는 별도의 공변량이 필요 없다. 왜냐하면 이미 공변량의 정보량을 통해 가중치를 부여한 의사 모집단을 만든 것이니까. 하지만 효과 수정 여부를 평가하기 위해서는 공변량을 추가할 수 있다. 예를 들어, 다음과 같은 한계구조평균 모델이 있다고 하자:

$E [Y^{a} | V] = β_{0} + β_{1} a + β_{2} V a + β_{3} V$

가산효과수정은 $β_{2} \neq 0$ 일 경우에 나타난다.
연구자가 한계구조모델에 포함하기로 선택한 $L$ 의 특정 서브셋 $V$ 는 연구자의 실질적인 관심사만을 반영해야 한다.
교환가능성에 필요한 모든 변수 $L$ 을 조건으로 하는 한계구조모델을 가짜 한계 구조모델이라고 한다.

12.6 Censoring and missing data

결측치가 없는 개개인들만을 상대로 분석을 수행한다던지 하는 일종의 데이터 센서링은 선택 편향을 야기할 수 있다.

결측 등으로 인한 센서링은 선택 편향으로 이어질 수 있어서 일반적으로는 연구하고자 하는 모집단에서 누구도 센서링 되지 않은 경우의 인과효과에 관심을 가진다.