의연방09) PSA 분석1(매칭)

📜 제목으로 보기

✏마지막 댓글로

• 학습목표
  • 주요용어
  • 목차
• Motivation(하는 이유)
  • Association vs. causality
  • Confounding을 다루는 법
  • Propensity score analysis
  • PSA요구하는 4가지 가정
  • Regression vs. PSA
  • Conditional vs. marginal
  • 2가지 남은 문제
• PS 추정 방법
  • PS를 추정하는 방법
  • PS model의 설명변수
• PS Matching
  • PS 추정 후 분석방법
  • PS matching
  • 매칭 방법
  • 매칭 후 balance 평가
  • Standardized mean difference (SMD)
  • 매칭 후 분석
• R 실습
  • matchit() from MatchIt
  • CreateTableOne() from tableone
  • geeglm() from geepack
  • 코드 실습
    • 이분형 결과변수 만들어주기
    • 연속형 -> 범주형변수로 한꺼번에 처리하기 by mutate_at
    • 매칭전 [치료군vs대조군(treat)별] confounder의심변수들의 balance(2군간의 차이 by p)확인 by tableone
      • 일단, 치료여부별 표1을 만들어서 confounder 유의미한 차이보기
      • smd를 이용하여 balance평가 with print( tobj, smd =T)
      • matching의 목적: 2군간의 confounder의 p차이 -> smd 차이를 0.2 미만으로 줄이는 것
    • 매칭 by matching
      • 매칭된 데이터만 뽑아내기 from mobj by match.data
      • 매칭데이터 balance(치료여부별 confounder의 p/SMD 차이) 평가 by CreateOneTable
      • random -> set.seed()
      • 1:2 매칭 by ratio=2
    • 매칭된 데이터에서 치료효과(y) 평가 by GEE
      • 0,1범주형을 -> 0,1 연속형으로 만드는 sense ( bool -> as.double )
      • 최종 회귀분석에서는 결과변수(질병여부, 연속형) ~ 설명변수(치료여부), id=매칭id변수(subclass) 2+1 개의 변수만
      • 회귀분석 결과obj로 OR를 구하기 by broom패키지의 tidy()
• 정리하기
  • 연습문제
    • 01 다음 중 Propensity score analysis에 대한 설명으로 틀린 것은?
    • 02 다음 중 PS(propensity score)를 추정한 후 자료 분석 방법으로 알맞지 않은 것은?

• R 고급 실습(ps matching with datacleaning-dplyr, MatchIt, tableone-CreateTableOne, geepack-geeglm, broom-회귀계수를 OR로 변환) 포함됨

학습목표

Propensity score analysis의 목적과 정의를 설명할 수 있다.

Propensity score의 정의를 알고, 추정하는 방법을 설명할 수 있다.

PS matching 방법을 이해하고, 직접 수행할 수 있다.

주요용어

• Propensity(프로펜서티) score analysis : confounding이 존재하는 관찰연구의 데이터를 RCT의 데이터처럼 confounding이 없도록 가공하여 분석하는 방법

• Propensity score : 치료군에 속할 확률

목차

1. Motivation
2. PS 추정 방법
3. PS Matching
   • 4가지 중 1가지 방법

Motivation(하는 이유)

Association vs. causality

• 위험요소-질병간의 상관관계 -> 인과관계를 살펴보고 있는데
  • 상관관계가 곧 인과관계를 의미하지는 않는다

1. 상관관계가 관측될 경우, 아래를 확인해야한다.
   • 이 상관관계가 진짜인가?
   • 우연이나 잘못된 연구디자인에 의한 것인가?
2. 상관관계가 진짜로 판단되는 경우
   • 인과관계가 성립하는가?
   • confounding에 의해 생겨난 상관관계인가?

Confounding을 다루는 법

• 연구 디자인 단계에서 confounding을 최대한 제거한다
  1. 랜덤화 -> RCT를 하는 것
  2. 매칭
• 데이터 분석 단계에서 confounding 효과를 보정한다
  1. Stratification
     • confounder 값별로 그룹을 다시 나눠서 분석
  2. Adjustment
     • multivariable model로 분석하는 것
  3. Propensity score analysis
     • confounding 효과를 보정하기 위해서
     • (마치 RCT에서 나온 것 같이) 데이터를 가공하여, X값을 랜덤하게 배정했을 때 관측될 법한 virtual data를 만들어 분석한다

Propensity score analysis

• Confounding이 존재하는 관찰연구의 데이터를 -> RCT의 데이터처럼 confounding이 없도록 가공하여 분석하는 것

• Propensity score (PS) : (치료군/대조군 비교 연구에서) 치료군에 속할 확률

  • 위험요소 - 질병 or 치료방법 - 질병 간의 관계를 살펴볼 때
    • 위험요소O/X or 치료O/X 2군으로 나누고, 질병 발생 빈도를 살펴보게 된다.

  • 2군으로 나누었을 때, 그 중에 1군(치료군)에 속할 확률을 PS라 한다.

    • RCT에서는, 모든 환자에 대해서 PS = 0.5(치료군 0.5 -> 대조군도 0.5)
    • 관찰연구에서는 PS(치료군에 속할 확률)는 알려지지 않은 (일정한 확률이 없는)미지의 값
    • 관찰연구에서 PS는 환자 특성(or 치료자, 환자가 속한 지역/병원)에 따라 달라짐
• PS를 사용하는 이유?
  • PS에 condition을 하고나면, baseline characteristics(환자 특성)와 treatment selection(치료군이냐 대조군이냐)이 독립이 된다.
  • 무슨 얘기냐?
    • PS가 비슷한 환자들끼리 모아놓으면, 치료군/대조군의 baseline characteristics 분포가 비슷해진다
  • 원래는 환자 특성(baseline characteristics) -> 환자가 치료군에 속할 확률(treatment selection)에 영향을 주게 되는데
    • 환자 특성을 일일히 확인해서 비슷한 사람들끼리 모아놓는게 아니라, 환자 특성에 영향 받는 PS만 알고 있다면, PS를 기준으로 비슷한 환자를 모아두기만 하면, 환자 특성 분포가 비슷해진다

PSA요구하는 4가지 가정

1. consistency(컨시스턴시)
   • 환자가 받은 치료법에 따라서, 치료군이든 대조군이든 관측될만한 outcome(질병)이 있다. 그것이 real world에서 관측한 것과 outcome(질병)이 같다
   • 환자의 potential outcome = 실제로 관측된 outcome
2. no unmeasured confounder
   • 모든 confounder를 알고 있고 관측했다는 가정
3. Positivity : 치료군이 될 확률＞0, 대조군이 될 확률＞0
   • 어떤 환자든지 한쪽에만 속하는게 아니라 2군 각각 속할 확률이 존재한다
4. no misspecification of PS model
   • PS score 추정모델이 잘못되어선 안된다.

• 3번 제외 다른 가정들은 만족시키기 어렵다
  • 모든 통계모델은 많은 가정아래 존재한다.
  • 그런데 항상 만족시키는 것은 아니다. 현실을 완벽하게 반영은 X

Regression vs. PSA

• PSA에 대한 사람들의 의문점
  • “Confounder는 공변량으로 Multivariable regression model에 넣으면 보정 되는 것 아닌가요?”
    • 이것도 유효한 방법이나, treatment effect가 조금 다르게 나온다.

1. Regression에서 추정하는 효과 : conditional treatment effect
   • 각 환자 개개인에서의 평균 효과
2. PSA와 RCT에서 추정하는 효과 : marginal treatment effect
   • 모집단 전체에서의 모든 환자의 평균 효과

Conditional vs. marginal

• 두 개는 뭐가 다른 것일까?

1. Conditional treatment effect in Regression

   • “한 환자에게 혈압강하제를 투여하는 실험을 여러번 반복할 때, 평균적으로 수축기 혈압을 10만큼 낮춘다”
     • 여기서 -10이 Conditional treatment effect이다.
     • 이 환자에게 약을 투약하는 실험을 여러번 반복해서 평균적으로 수축기혈압을 10만큼 낮춘다.
       • 개인에서의 평균

2. Marginal treatment effect in PSA, RCT

   • “고혈압 환자 전체에서 평균적으로 수축기 혈압을 10만큼 낮춘다”
     • 여기서 -10이 ConditMarginal treatment effect이다.
     • 각 환자에게 낮춰진 수축기혈압을 모두 구해, 모든 고혈압 환자에서 평균적으로 10만큼 낮춘다.
       • 전체에서의 평균

• 한 환자당 10만큼 낮추면(conditional) -> 평균적으로 전체에서도 10만큼 낮춰질 수 밖에 없다.(marginal)

  • 이 것이 가능한 상황은,
    • 현재 연속현 변수(수축기혈압)의 평균으로 추정하거나
    • 아니면, 이분형 변수에서 치료군vs대조군 비율의 차이(Risk Differenece)로 계산할 경우에
    • conditional == marginal treamment effect 로 항상 같다
  • 하지만,
  • 효과가 연속변수의 평균 차이 또는 이분형 변수의 비율차이(RD)가 아니라 OR 또는 생존분석의 Hazard Ratio로 측정되는 경우도 있디.
    • conditional tx effect ≠ marginal tx effect로 다르다.

• 그럼 무엇을 써야할지 어떻게 결정할까?

  • 내가 추정할 치료효과가 RCT에서 나타나는 결과와 유사한가?
    • yes라면, Multivariable regression가 아닌 Propensity score analysis가 더 적합하다 판단하고 Marginal treatment effect in PSA, RCT 를 구하면 된다.

2가지 남은 문제

1. PS를 어떻게 추정할 것인가
2. PS가 비슷한 환자들은 confounder가 없는 상태라고 하는데, 비슷한 환자들은 어떻게 모아서 치료효과를 추정할 것인가?

PS 추정 방법

PS를 추정하는 방법

• 보통 로지스틱 회귀분석을 이용한다
  • Propensity score는 치료군에 속할 확률이라고 했다.
  • 치료군에 속하느냐? vs 대조군에 속하느냐?의 이분형 변수로 하는 로지스틱 회귀분석을 이용해서 추정한다.
  • 그렇다고 꼭 로지스틱 회귀분석을 사용하는 것이 아니라,
    • 이분형 결과변수를 사용하는 모든 종류의 분석방법 사용가능
      • 예) SVM, neural network, tree, random forests
      • 결과변수(종속변수, 반응변수) : 실험군(치료군)/대조군
• 현실적으로는 로지스틱 회귀분석이 가장 많이 쓰인다.

PS model의 설명변수

• propensity score를 구하기 위해 regression model에 어떤 설명변수를 넣어야할까?
  • PS 개념상으로는 이 사람이 치료군에 속할지/대조군에 속할지 결정(X, 치료여부)에 영향을 끼치는 변수들을 모두 설명변수에 변수들(아래 그림 중 가운데에서 C)을 넣는 것이 맞는 것 같지만...
  • X(어느군에 속할지 == 치료여부) -> Y(outcome)
• 추가로, 여러 시뮬레이션을 통한 여러 연구결과:
  • 치료군/대조군 결정(X)에 영향을 끼치지 않더라도 outcome(Y)에 영향을 끼치는 변수들(3번째 그림)은 넣는 것이 좋다라고 연구 결과가 나왔다.
    • 분산을 감소시킴 = p-value가 낮아짐 && 신뢰구간 좁아짐 등의 효과**
      • X(치료여부)에는 영향을 끼치지 않더라도 Y(outcome)에 영향을 끼친다면 넣자.
  • outcome(Y)에 영향을 끼치지 않으면서 치료군/대조군 결정(X)에 영향을 끼치는 변수들은 넣지 않는 것이 좋다 (분산을 증가시킴)
    • PS정의상으로는 넣으면 좋을 것 같지만
    • 연구결과 넣지 않는 것이 더 좋다고 나옴
    • 가운데 그림의 C -> 넣지 않는 것이 좋다라고 연구결과가 나온 것도 있다.

• 그럼 어떤 것을 PS의 regression model에 어떤 설명변수를 넣어야할까?

• 보통 Table 1에 들어가는 baseline characteristics는 치료군/대조군 결정(X)과 outcome(Y) 둘 다에 영향을 준다
  • → 다 포함하는 것이 좋다
  • table1의 기본특성은 거의다 confounder로 생각된다. 다 포함시키면 된다.
• PS model은 overfitting 문제가 없다
  • 원래 다변수 모델에서는
    • 설명변수가 많아지면 -> 데이터에 비해 모델이 복잡해져 overfitting이 발생하지만
  • PS model에서는 그런 overfitting 문제가 없다
    • → 설명변수 개수에 제한이 없다
  • PS model은 우리가 알고 싶어하는 것에 대한 최종 model이 아니라 confounder 보정의 중간단계이기 때문에, 넣고 싶은만큼 많이 넣어도 된다.
    • 조금이라도 confounder 같다고 의심되는 변수는 모조리 포함시키면 된다.
• 절대 포함하면 안되는 변수
  • 치료군/대조군 결정 이후 관측된 변수
  • PS의 정의= 치료군에 속할 확률
    • 치료군/대조군 결정되기 이전의 변수들로 판단되어야한다
    • 군 결정이후 발생된 변수는, 과거(군 결정)에 영향을 주면 안된다.
  • cf) 관측만 결정 이후에 되었을뿐이지, 그 변수값 자체는 치료군vs대조군 결정이전에 결정된 경우도 있다.
    • 예) 2가지 수술방법에 대한 비교 연구
      • 수술이후 조직을 biopsy한 결과 -> 치료군vs대조군 결정이후 = 수술방법 결정이후 관측되었지만
      • 수술이전에 이미 내부장기상태는 결정된 상태라고도 볼 수 있다.

PS Matching

• PS model을 만들어서 PS를 추정했다고 가정하자
  • 그 이후 다음 단계 분석은 어떻게 해야할까?

PS 추정 후 분석방법

- Matching
- Stratification
- Weighting
- Covariate adjustment(코베리에잇 어저스먼트)

• 예제
  • 원의 크기 = 환자별 치료군에 속할 확률(PS)
    • 대조군에는 원의 크기가 작은 환자들이 많고
    • 치료군에는 원의 크기가 큰 환자들이 많다.

PS matching

• PS가 비슷한 환자별로 치료군:대조군 1:1 쌍을 지어주는 것
  • 짝이 지어진 환자들만 분석하고, 짝이 안지어지면 분석에 사용안한다
  • 2군에 걸쳐 짝이 지어졌다 = 분포가 동일해진다. = 전체적으로 PS가 비슷해짐 = 다른 기본특성도 비슷해짐 = confounding이 제거된다.
  • 장점 : 직관적으로 이해가 쉽다
  • 단점 : 짝지어지지 않는 환자의 데이터를 잃어버린다
    • n은 클수록 좋은데, 모아진 데이터를 버려야하니 큰 단점

매칭 방법

• With replacement vs. without replacement

  • With replacement: 매칭해서 짝 지은 다음 -> 복원 추출처럼 다른 환자랑 매칭가능하게 하는 것
    • 통계적으로 더 좋은 성질을 가지지만, 여러번 들어갈 수 있어서 -> 사람들이 껄끄럽게 생각함
  • withoutreplacement: 비복원추출 -> 매칭후 분석풀에서 재고
    • 실제로 많이 사용되는 것

• Greedy vs. optimal matching vs. full matching

  • Greedy: 한명씩 순차적으로 매칭
    • 가장 많이 쓰임.
  • optimal: 전체적 2군의 PS가 비슷해지도록 한번에 매칭
  • full: 버림없이 모든 환자를 매칭

• 1:1 vs. n:1

  • 대조군이 너무 많을 경우 -> 대조군:치료군 = 4:1까지 매칭된다.

• 오스틴이 밝힌 가장 좋은 매칭 방법은?

  • nearest neighbor matching(=Greedy) with caliper
    • 치료군 환자에 대해 가장 PS가 비슷한 환자를 뽑더라도, 가장 비슷한 PS를 뽑더라도, 치료군의 PS와 차이가 많이 나는 대조군이 뽑을 수 있다.
  • PS차이의 한계를 정해놓은 것 = caliper
    1. logit of PS을 먼저 계산한다
       • log 𝑃𝑆/(1-𝑃𝑆)
    2. logit of PS의 표준편차를 계산하고, 여기에 0.2를 곱한 값을 caliper로 생각한다.
       • 0.2 * SD of logit of PS
    3. 위의 caliper를 치료군-대조군 matching시 PS차이 한계점으로 생각(넘으면 버림)하고 nearest neighbor matching을 돌리면, confounder로 인한 치료효과 추정치에 생기는 bias가 99%가 제거된다.

• 이로 인해 nearest neighbor(=Greedy) matching with caliper를 가장 많이 쓴다.

매칭 후 balance 평가

• 매칭을 한 뒤, 잘 매칭되었는지 확인하기 위해
  • PS matching -> 목적:confounding이 없어졌는지 == baseline characteristic 차이가 정말 줄어들었는지를 확인해야한다
  • 즉, baseline characteristics의 차이가 줄어들었는지 평가해야한다
• 이 때, 매칭된 dataset에 대해
  • 잘못된 방식으로서
    • 치료군과 대조군의 각 baseline characteristics의 요약통계량 계산 후
    • 2군의 비교이므로 p-value로 계산하면 될까?
      • 2군 비교시 p-value가 익숙해져 있어서 매칭데이터도 p-value로 비교하는 경우가 많다 -> 좋은 방법이 아니다.
• 매칭된 dataset을 p-value로 비교하면 왜 잘못된 방식일까?
  1. p-value는 sample size에 영향을 받지만
  2. matching의 경우, 대부분 매칭되지 않는 환자들이 제외되어
     • original sample size보다 줄어들게 되고
     • power가 줄어든다 -> 두 군간의 차이가 있었더라도 -> p-value는 높게 나올 수 있다.
     • p-value가 높게 나온 것이, 줄어든 sample size 때문인데, 마치 matching이 잘되어서 p-value가 나온 것처럼 잘못된 impression을 줄 수 있다.
  • p-value로 비교한다는 것은
    • 2군의 sample에 대해 2개의 큰 모집단이 있다는 것을 가정한 것이다.
    • 2개 모집단을 비교할 수 없어서 random sample을 뽑아, 2개 모집단을 추론하는 것이 p-value이다.
  • 그러나 우리는 matching 했다
    • 2군 뒤의 2개 모집단에 대한 관심은 잃은 상태다.
    • 2개 모집단은 애초에 baseline characteristic이 다르기 때문에, confounding으로 인한 치료효과 분석이 어렵다고 생각해서 인위적으로 matching을 하는 것
      • 인위적으로 baseline characteristic를 맞춘, 인위적은 dataset을 보는 것
      • 더이상 sample로 모집단을 추론(p-value)하는 것이 아니라 sample 자체에 관심을 가지게 되는 것이므로 p-value로 2군 비교는 적합하지 않다.
      • 개념적으로 아예 맞지 않는 것
• 그렇다면 옳은 매칭후 평가 방법은?
  • matching후 평가방법은 여러가지가 있다. 그 중에서 사용하기 편하고 이해하기 쉬운 방법인 Standardized mean difference를 통해 평가해보자

Standardized mean difference (SMD)

1. 연속형 변수의 SMD

   • 복잡해보이지만, 치료군의 평균(X_bar_tx)와 대조군의 평균(X_bar_tx)의 차이에서

   • 치료군과 대조군의 표준편차를 잘 합친 것(pooled SD)으로 나누어, 데이터의 변동을 고려했을 때의 평균 차이

     

1. 이분형 변수의 SMD
   • 2 군의 비율 차이를 적당한 표준편차 값으로 나눈 것으로 생각하자

• 두 군간 평균의 차이가 pooled SD의 몇배인가
• SMD < 0.1~0.2 이면 작은 차이로 간주
• 매칭 이후 매칭변수들(PS model에 사용된 설명변수들)에 대해 SMD를 계산하고
  • 대부분 0.1미만이면 좋다
  • 몇개는 0.1을 넘어 모두 0.2 미만이면 balance가 잘 맞는 것으로 평가 －　이정도면，매칭된　dataset에서　치료군과　대조군　사이에　baseline characteristic차이가　무시할만한　차이이다

매칭 후 분석

• 분석을 하려면, 우리가 애초에 궁금한 치료효과에 대해 평가해야한다.

  • 치료여부-질병여부와의 association이 얼마나 강한가?를 평가하는 measure of association를 이용해야한다.
  • 여러가지 방법이 있었다.
    • OR말고 RD(risk diffrence) = attributable risk라던가 relative risk를 계산할 수 있지만
    • OR이 가장 계산이 쉽고, 편리하기 때문에 대표로 사용해보자.

• 매칭된 치료군과 대조군에서 결과변수의 차이

  • 이분형 결과변수의 경우: OR과 OR의 95% 신뢰구간, OR이 1인지 아닌지에 대한 p-value 정도를 사용한다

• 보통의 방법을 사용하면 되는가?

  • 매칭된 dataset은, 매칭되면서 짝끼리 서로 비슷하게 맞춰놓은 상태기 때문에, 환자가 서로 독립이라고 할 순 없다. 그래서 대부분의 통계방법은 서로 독립이라고 가정을 하는데, 일반적인 통계방법으로 분석하면 안된다. -> 매칭된 짝끼리는 독립이 아니라는 점을 고려해서 분석하는 것이 좋다는 것이 다수 의견이다.

• 매칭된 짝은 서로 독립이 아니라는 점을 고려하여 분석한다

  • 이분형 결과변수의 경우:
    • McNemar’s test : p-value를 계산해주지만, OR자체를 계산해주지 않는다.
      • OR 계산하려면 logistic regression을 해야한다.
    • conditional logistic regression
      • case-control study에서 사용했던 것. 매칭된 짝끼리 비슷한 점 = 독립이 아니라는 점을 이용할 수 있다.
    • logistic regression with correlation structure (GEE)
      • 일반 logistic regression을 하되, 매칭된 짝끼리의 correlation을 허용하면서 모델을 fit하는 것
      • GEE model이라고 하는데, 사실은 계산방법이다.
      • GEE 계산을 통한 로지스틱 회귀분석을 통해 OR을 계산하고, 신뢰구간을 구하면 된다.

• 소수의견 : RCT에서 얻어진 데이터처럼 만들었기 때문에, RCT처럼 독립성을 가정하고 분석해야 한다는 관점도 존재한다

• conditional logistic regression이나 logistic regression with correlation structure (GEE)을 수행할 때
  • 결과 변수는 질병여부
  • 설명 변수는 치료군vs대조군 여부가 되는데
  • 애초에 confounder라고 생각했던 공변량들은 설명변수로 포함할 필요가 없게 된다.
    • 왜냐면, PS matching을 통해 공변량 효과를 이미 제거된 상태가 되었기 때문이다.
    • 즉, multivarable이 아니라 univariable의 로지스틱 회분석해서 OR(오즈 레이시오)를 계산하면 된다는 것이다.
• 그럼에도 불구하고, 공변량으로 포함시켜 계산하는 방법도 존재하는데 doubly robust하다라는 이름을 얻게 되었다.
  • PS model이 잘못됬거나 아니면, 최종 모델로 생각한 회귀분석이 잘못되었어도, OR 추정치가 unbiased하게 되는 좋은 성질을 갖게 된다고 한다.
  • 하지만, 매칭 분석은 univariable model만 쓰면 된다.

R 실습

• MatchIt 패키지를 이용한 PS matching
  • matchit() 함수 이용
  • Austin이 추천한 nearest neighbor matching with caliper = 0.2*SD of logit of PS 을 구현하는 방법 실습
• tableone 패키지를 이용한 SMD 계산
  • CreateTableOne() 함수 이용
  • 매칭후 balance 계산을 위해 SMD를 계산.
• geepack 패키지를 이용한 OR 추정
  • geeglm() 함수 이용
  • 매칭된 데이터셋에서 GEE방법으로 OR를 계산해주는 패키지

matchit() from MatchIt

• PS matching이외에 여러가지 매칭이 가능한 패키지이나 필요한 것만 설정하기

  matchit(formula, data, method, replace, ratio, distance, 
  link, caliper, std.caliper, m.order)

• formula : 매칭하려는 그룹 (치료군/대조군)~ 매칭변수1+매칭변수2+...+매칭변수d
  • 치료군이냐 대조군이냐 그룹변수 ~ confounder라고 의심해서 PS model에 넣고 싶은 변수들을 더하기로 연결
• data : 데이터셋 이름
• method : 매칭방법. “nearest”라고 입력하면 nearest neighbor matching을 수행한다
• replace : 재복원여부.
  • 디폴트는 FALSE (비복원표집, 한번 매칭된 환자는 다시 매칭되지 않는다)
  • without replacement = 한번 매칭된 환자는 다시 매칭하지 않는 것이 default
• Ratio : 매칭비율 (치료군 1명당 매칭되는 대조군 환자 수). 디폴트는 1.
  • 치료군 1명당 매칭될 대조군의 수

• distance=“glm”
  • matching 중에 Propensity Score matching을 이용하는 옵션
• link=“logit”
  • caliper의 기준을 정하는데, PS의 logit를 거리측도로 사용하여 매칭을 수행한다
• caliper=0.2
  • caliper 기준에 곱해줄 값을 입력
• std.caliper = TRUE
  • caliper 의 표준편차를 기준으로 매칭
  • 위 4가지를 다 설정해야, austin의 PS login의 SD의 0.2를 caliper로 사용한 매칭을 실시하게 된다
• m.order=“random”: 매칭을 수행하는 순서
  • 입력하지 않으면 거리측도 값 = PS값이 큰 개체부터 순서대로 매칭한다.
    • 그러나 PS score 큰 수대로 매칭해봐야 장점이 없다. 랜덤으로 매칭하도록 선호된다.
  • “random”으로 입력하면 랜덤한 순서로 매칭하는데, 이 경우 matchit() 함수 실행 직전에 set.seed()함수로 seed number를 정해주어야 같은 결과를 재생산할 수 있다. (그렇지 않으면 matchit() 함수 실행 때마다 다른 매칭결과가 나온다

CreateTableOne() from tableone

CreateTableOne(vars, data, strata)

• 임상과학 논문의 Table 1에 해당하는 표를 편리하게 만들 수 있는 함수
  • vars : 표에 들어갈 변수 이름의 벡터
  • data : 데이터셋 이름
  • strata : 비교할 그룹을 나타내는 변수
    • 스트레이라를 입력해주면, 그룹을 나눠 그룹별 요약통계량과 그룹별 p-value가 나온다.
• CreateTableOne()으로 만든 객체를 print() 안에 넣으면 표가 인쇄된다
  • smd : TRUE로 입력하면 SMD를 계산해준다
  • test : FALSE로 입력하면 p-value를 인쇄하지 않는다

geeglm() from geepack

geeglm(formula, id, family, corstr, std.err, data)

• GLM(generalized linear model)을 GEE(generalized estimating equation) 방법으로 correlation structure를 주면서 fit하는 함수
• formula : 결과변수~설명변수
  • fitting하고 싶은 model 입력

• id : cluster를 나타내는 변수

  • PS maching의 경우 같이 매칭된 자료의 경우 matching id를 넣어주면 된다.
  • 내가 생각하기에 correlate되어있는 클러스트를 나타내는 변수를 넣는다.

• family : GLM의 종류. 로지스틱 회귀분석의 경우 “binomial”

  • 로지스틱 회귀를 설정 옵션
• corstr : cluster 내에서의 correlation structure.
  • 어떤 correlation structure를 줄 것이냐
  • “exchangeable”: cluster내에서 모든 개체 쌍의 correlation이 동일할 경우
    • 매칭된 환자(치료군 1명에 대조군 1~4)는 같은 cluster내에 있으며, 매칭 환자들 중에 아무 2명씩 짝을 지어도 항상 correlation같다는 옵션
    • 가장 무난한 옵션
  • “independence”: cluster내에서 모든 개체가 서로 독립일 경우
    • 매칭된 cluster환자들 내(치료군+매칭된 대조군(들))이 서로 독립이라고 가정하고 싶을 때 쓰는 옵션이나
    • 매칭되었으니까 independent할텐데, 왜 쓸까?
      • 어떤 복잡한 통계적 방법에 의해 이 방법이 더 타당하다고 주장/연구결과가 있어서 옵션으로 언급만
• std.err : standard error 계산방식.
  • 디폴트 “san.se”는 robust estimator방식으로 계산하게 되는데, 일반적인 방식이 아니라 데이터가 매칭=독립이 아니다라고 가정하여, 더 보수적으로 계산한 estimator다
  • 95% 신뢰구간과 p-value 계산시 영향을 끼친다
  • robust estimator를 써야 더 보수적인 결과가 나오며 default옵션이다.
• data : 데이터셋 이름

코드 실습

# install.packages("MatchIt") -> 설치 안되서 터미널에서 아래와 같이
# conda install -c conda-forge r-matchit
# - ocnda로 검색해서 안되서, stackoverflow에서 얻은 설치주소
# - https://stackoverflow.com/questions/71041174/installing-r-package-matchit-in-anaconda-environment

library(MatchIt)

# - matchIt 다른 곳에서 설치 -> R을 4로 업데이트 -> 추가적인 것 수동 업데이트

library(dplyr)
library(tableone)

다음의 패키지를 부착합니다: 'dplyr'


The following objects are masked from 'package:stats':

    filter, lag


The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

data("lalonde")

head(lalonde)
# treat	age	educ	race	married	nodegree	re74	re75	re78
# 치료여부 / / 교육 / 결혼여부 / nodgree 고등학교졸업/ 년도별 연소득
# 의료데이터는 아니나,, 가정하고 한다.

A data.frame: 6 × 9

	treat	age	educ	race	married	nodegree	re74	re75	re78
	<int>	<int>	<int>	<fct>	<int>	<int>	<dbl>	<dbl>	<dbl>
NSW1	1	37	11	black	1	1	0	0	9930.0460
NSW2	1	22	9	hispan	0	1	0	0	3595.8940
NSW3	1	30	12	black	0	0	0	0	24909.4500
NSW4	1	27	11	black	0	1	0	0	7506.1460
NSW5	1	33	8	black	0	1	0	0	289.7899
NSW6	1	22	9	black	0	1	0	0	4056.4940

이분형 결과변수 만들어주기

• re78이 5000을 넘는지 안넘는지로

dat1 <- lalonde %>%
    # outcome 변수이름을 간단히 y로 둔다. 
    # re78이 5000달러를 넘으면 [질병여부 발생여부]로 임시로 가상으로 만든다.
    mutate(y=ifelse(re78 > 5000, 1, 0))

head(dat1)

A data.frame: 6 × 10

	treat	age	educ	race	married	nodegree	re74	re75	re78	y
	<int>	<int>	<int>	<fct>	<int>	<int>	<dbl>	<dbl>	<dbl>	<dbl>
NSW1	1	37	11	black	1	1	0	0	9930.0460	1
NSW2	1	22	9	hispan	0	1	0	0	3595.8940	0
NSW3	1	30	12	black	0	0	0	0	24909.4500	1
NSW4	1	27	11	black	0	1	0	0	7506.1460	1
NSW5	1	33	8	black	0	1	0	0	289.7899	0
NSW6	1	22	9	black	0	1	0	0	4056.4940	0

dat1 <- lalonde %>%
    mutate(y=ifelse(re78 > 5000, 1, 0))

연속형 -> 범주형변수로 한꺼번에 처리하기 by mutate_at

# - y도 0 or 1로 바꿔줬는데, 연속형변수(숫자)처럼 되어있을 것이다.
# - 연속형 변수는 min/max가 나와있다.
# - 치료여부(treat) 와 질병여부(y)가 모두 연속형 변수로 되어있다.
#   그외
summary(dat1)

     treat             age             educ           race        married      
 Min.   :0.0000   Min.   :16.00   Min.   : 0.00   black :243   Min.   :0.0000  
 1st Qu.:0.0000   1st Qu.:20.00   1st Qu.: 9.00   hispan: 72   1st Qu.:0.0000  
 Median :0.0000   Median :25.00   Median :11.00   white :299   Median :0.0000  
 Mean   :0.3013   Mean   :27.36   Mean   :10.27                Mean   :0.4153  
 3rd Qu.:1.0000   3rd Qu.:32.00   3rd Qu.:12.00                3rd Qu.:1.0000  
 Max.   :1.0000   Max.   :55.00   Max.   :18.00                Max.   :1.0000  
    nodegree           re74            re75              re78        
 Min.   :0.0000   Min.   :    0   Min.   :    0.0   Min.   :    0.0  
 1st Qu.:0.0000   1st Qu.:    0   1st Qu.:    0.0   1st Qu.:  238.3  
 Median :1.0000   Median : 1042   Median :  601.5   Median : 4759.0  
 Mean   :0.6303   Mean   : 4558   Mean   : 2184.9   Mean   : 6792.8  
 3rd Qu.:1.0000   3rd Qu.: 7888   3rd Qu.: 3249.0   3rd Qu.:10893.6  
 Max.   :1.0000   Max.   :35040   Max.   :25142.2   Max.   :60307.9  
       y         
 Min.   :0.0000  
 1st Qu.:0.0000  
 Median :0.0000  
 Mean   :0.4853  
 3rd Qu.:1.0000  
 Max.   :1.0000  

dat1 <- lalonde %>%
    mutate(y=ifelse(re78 > 5000, 1, 0)) %>%
    # vars(변수들), list(~하고싶은작업(.))   . 자리에 변수들이 들어간다.
    mutate_at(vars(treat, married, nodegree, y), 
             list(~as.factor(.)))

summary(dat1)

 treat        age             educ           race     married nodegree
 0:429   Min.   :16.00   Min.   : 0.00   black :243   0:359   0:227   
 1:185   1st Qu.:20.00   1st Qu.: 9.00   hispan: 72   1:255   1:387   
         Median :25.00   Median :11.00   white :299                   
         Mean   :27.36   Mean   :10.27                                
         3rd Qu.:32.00   3rd Qu.:12.00                                
         Max.   :55.00   Max.   :18.00                                
      re74            re75              re78         y      
 Min.   :    0   Min.   :    0.0   Min.   :    0.0   0:316  
 1st Qu.:    0   1st Qu.:    0.0   1st Qu.:  238.3   1:298  
 Median : 1042   Median :  601.5   Median : 4759.0          
 Mean   : 4558   Mean   : 2184.9   Mean   : 6792.8          
 3rd Qu.: 7888   3rd Qu.: 3249.0   3rd Qu.:10893.6          
 Max.   :35040   Max.   :25142.2   Max.   :60307.9          

매칭전 [치료군vs대조군(treat)별] confounder의심변수들의 balance(2군간의 차이 by p)확인 by tableone

• 4가지 변수를 confounder라고 가정하고 만들어준다.
• 치료군 vs 대조군 = 치료여부를 의미하는 treat별로 봐야한다.
• CreateTableOne()이란 치료군vs대조군변수를 strata로 넣어, 표1을 생성하는 메서드다.

일단, 치료여부별 표1을 만들어서 confounder 유의미한 차이보기

# install.packages("haven")
# install.packages("forcats")
# install.packages("hms")

tobj <- CreateTableOne(vars=c("age", "educ", "married", "nodegree"),
                      data = dat1)
tobj

                  
                   Overall      
  n                  614        
  age (mean (SD))  27.36 (9.88) 
  educ (mean (SD)) 10.27 (2.63) 
  married = 1 (%)    255 (41.5) 
  nodegree = 1 (%)   387 (63.0) 

# -> treat 변수별로 보도록 strata를 treat로 지정해준다.
tobj <- CreateTableOne(vars=c("age", "educ", "married", "nodegree"),
                      data = dat1, strata = "treat")
tobj

                  Stratified by treat
                   0             1             p      test
  n                  429           185                    
  age (mean (SD))  28.03 (10.79) 25.82 (7.16)   0.011     
  educ (mean (SD)) 10.24 (2.86)  10.35 (2.01)   0.633     
  married = 1 (%)    220 (51.3)     35 (18.9)  <0.001     
  nodegree = 1 (%)   256 (59.7)    131 (70.8)   0.011     

tobj <- CreateTableOne(vars=c("age", "educ", "married", "nodegree"),
                      data = dat1, strata = "treat")
print(tobj)

# treat가 0 , 즉 대조군에 대한 4가지 confounder변수의 요약통계량이 나온다.
# 치료군에 대해서도 나온다.
# 치료군 185, 대조군 429명의 정보도 나온다.

# p-value를 보고, 2군간의 confounder( baseline characteristic)의 차이를 표1을 만들어서 본다.
# -> "p가 낮으면 두 군간의 유의미한 차이"가 있다는 뜻이다.
# ->age  married  nodegree 는 2군간에 유의미하게 차이가 난다.

                  Stratified by treat
                   0             1             p      test
  n                  429           185                    
  age (mean (SD))  28.03 (10.79) 25.82 (7.16)   0.011     
  educ (mean (SD)) 10.24 (2.86)  10.35 (2.01)   0.633     
  married = 1 (%)    220 (51.3)     35 (18.9)  <0.001     
  nodegree = 1 (%)   256 (59.7)    131 (70.8)   0.011     

smd를 이용하여 balance평가 with print( tobj, smd =T)

• smd는 매칭 후에 평가한다고 했는데, 매칭전 smd를 구해보면, 차후 얼마나 줄었는지 알 수 있다.

tobj <- CreateTableOne(vars=c("age", "educ", "married", "nodegree"),
                      data = dat1, strata = "treat")
print(tobj, smd = T)

# 0.2를 넘는 것(married)은 두 군간에 상당히 다르다는 것을 의미한다.

                  Stratified by treat
                   0             1             p      test SMD   
  n                  429           185                           
  age (mean (SD))  28.03 (10.79) 25.82 (7.16)   0.011       0.242
  educ (mean (SD)) 10.24 (2.86)  10.35 (2.01)   0.633       0.045
  married = 1 (%)    220 (51.3)     35 (18.9)  <0.001       0.721
  nodegree = 1 (%)   256 (59.7)    131 (70.8)   0.011       0.235

matching의 목적: 2군간의 confounder의 p차이 -> smd 차이를 0.2 미만으로 줄이는 것

매칭 by matching

• 매칭 결과는 바로 쓰는 것이 아니라 mobj라는 객체에 저장해놓고 쓰는 것이 좋다.
• 목적은 confounder들의 차이를 smd 0.2미만으로 줄여서, 2군간의 차이를 줄여놓는 것

mobj <- matchit(treat ~ age+educ+married+nodegree , # formula : 치료여부 ~ confoiunder들
                data = dat1, # dataset
                method = "nearest", # 매칭 방법: austin의 nearest ~
                distance = "glm", # 매칭 중 PS matching
                link = "logit", 
                replace = F, # 매칭된 환자 비복원
                capliper = 0.2, std.caliper = T, # caliper 지정
                m.order = "random", #매칭순서 -> default방법 안좋음. random으로 지정
               )

# 뭐가 들었는지 일단 summary 안에 넣어본다.
summary(mobj)

# 썩 보기 좋은게 아님.
# 맨 마지막만 보고, 몇 쌍이 매칭 되어있는지 확인한다.
# Sample Sizes:
#           Control Treated
# All           429     185
# Matched       185     185
# Unmatched     244       0
# Discarded       0       0

# -> 185쌍이 매칭되었다.

Call:
matchit(formula = treat ~ age + educ + married + nodegree, data = dat1, 
    method = "nearest", distance = "glm", link = "logit", replace = F, 
    m.order = "random", std.caliper = T, capliper = 0.2)

Summary of Balance for All Data:
          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance         0.3772        0.2686          0.8386     0.7471    0.2080
age             25.8162       28.0303         -0.3094     0.4400    0.0813
educ            10.3459       10.2354          0.0550     0.4959    0.0347
married0         0.8108        0.4872          0.8263          .    0.3236
married1         0.1892        0.5128         -0.8263          .    0.3236
nodegree0        0.2919        0.4033         -0.2450          .    0.1114
nodegree1        0.7081        0.5967          0.2450          .    0.1114
          eCDF Max
distance    0.3602
age         0.1577
educ        0.1114
married0    0.3236
married1    0.3236
nodegree0   0.1114
nodegree1   0.1114


Summary of Balance for Matched Data:
          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance         0.3772        0.3745          0.0210     1.0358    0.0085
age             25.8162       25.1730          0.0899     0.5492    0.0688
educ            10.3459       10.6054         -0.1290     0.7558    0.0171
married0         0.8108        0.8108          0.0000          .    0.0000
married1         0.1892        0.1892          0.0000          .    0.0000
nodegree0        0.2919        0.3297         -0.0832          .    0.0378
nodegree1        0.7081        0.6703          0.0832          .    0.0378
          eCDF Max Std. Pair Dist.
distance    0.0865          0.0378
age         0.2378          0.8348
educ        0.0541          0.6130
married0    0.0000          0.0000
married1    0.0000          0.0000
nodegree0   0.0378          0.2735
nodegree1   0.0378          0.2735

Percent Balance Improvement:
          Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
distance             97.5       87.9      95.9     76.0
age                  70.9       27.0      15.4    -50.8
educ               -134.8       60.1      50.8     51.5
married0            100.0          .     100.0    100.0
married1            100.0          .     100.0    100.0
nodegree0            66.0          .      66.0     66.0
nodegree1            66.0          .      66.0     66.0

Sample Sizes:
          Control Treated
All           429     185
Matched       185     185
Unmatched     244       0
Discarded       0       0

매칭된 데이터만 뽑아내기 from mobj by match.data

m.data <- match.data(mobj)

summary(m.data)

# ps를 의미하는 distance가 보인다. 
# subclass -> 짝 지어진 번호. 1인 환자 -> 다른데 또 1인 매칭환자가 있다
# - 즉, matching Id를 의미한다.

 treat        age             educ           race     married nodegree
 0:185   Min.   :16.00   Min.   : 3.00   black :207   0:300   0:115   
 1:185   1st Qu.:19.00   1st Qu.: 9.00   hispan: 41   1: 70   1:255   
         Median :23.00   Median :11.00   white :122                   
         Mean   :25.49   Mean   :10.48                                
         3rd Qu.:29.00   3rd Qu.:12.00                                
         Max.   :55.00   Max.   :18.00                                
                                                                      
      re74            re75              re78         y          distance      
 Min.   :    0   Min.   :    0.0   Min.   :    0.0   0:196   Min.   :0.08913  
 1st Qu.:    0   1st Qu.:    0.0   1st Qu.:  260.9   1:174   1st Qu.:0.30744  
 Median :    0   Median :  310.6   Median : 4300.7           Median :0.41360  
 Mean   : 2979   Mean   : 1779.0   Mean   : 6340.2           Mean   :0.37586  
 3rd Qu.: 3536   3rd Qu.: 2557.9   3rd Qu.: 9588.5           3rd Qu.:0.47373  
 Max.   :35040   Max.   :25142.2   Max.   :60307.9           Max.   :0.58644  
                                                                              
    weights     subclass  
 Min.   :1   1      :  2  
 1st Qu.:1   2      :  2  
 Median :1   3      :  2  
 Mean   :1   4      :  2  
 3rd Qu.:1   5      :  2  
 Max.   :1   6      :  2  
             (Other):358  

매칭데이터 balance(치료여부별 confounder의 p/SMD 차이) 평가 by CreateOneTable

mobj <- CreateTableOne(vars=c("age", "educ", "married", "nodegree"),
                      data = m.data, strata = "treat")
print(mobj, smd = T)

                  Stratified by treat
                   0             1             p      test SMD   
  n                  185           185                           
  age (mean (SD))  25.17 (9.66)  25.82 (7.16)   0.467       0.076
  educ (mean (SD)) 10.61 (2.31)  10.35 (2.01)   0.250       0.120
  married = 1 (%)     35 (18.9)     35 (18.9)   1.000      <0.001
  nodegree = 1 (%)   124 (67.0)    131 (70.8)   0.500       0.082

# print( , test = F)
print(mobj, smd = T, test = F)

# 매칭후에는 smd가 0.2보다 다 작다

                  Stratified by treat
                   0             1             SMD   
  n                  185           185               
  age (mean (SD))  25.17 (9.66)  25.82 (7.16)   0.076
  educ (mean (SD)) 10.61 (2.31)  10.35 (2.01)   0.120
  married = 1 (%)     35 (18.9)     35 (18.9)  <0.001
  nodegree = 1 (%)   124 (67.0)    131 (70.8)   0.082

random -> set.seed()

• 원래는.. seed를 안정해주면, 매칭된 환자수가 바뀐다.
• 결과 자체는 문제가 없으나 재현성이 없어져버린다.

mobj <- matchit(treat ~ age+educ+married+nodegree , # formula : 치료여부 ~ confoiunder들
                data = dat1, # dataset
                method = "nearest", # 매칭 방법: austin의 nearest ~
                distance = "glm", # 매칭 중 PS matching
                link = "logit", 
                replace = F, # 매칭된 환자 비복원
                capliper = 0.2, std.caliper = T, # caliper 지정
                m.order= "random",#매칭순서 -> default방법 안좋음. random으로 지정
               )

m.data <- match.data(mobj)

tobj <- CreateTableOne(vars=c("age", "educ", "married", "nodegree"),
                      data = m.data, strata = "treat")

print(tobj, smd = T, test=F)

# 원래는 seed를 정해주지 않으면, 매칭 n수가 바귄다.
# -> 여기서는 m.order = random을 바꾸면 SMD가 바껴서 적용은 되는 것 같다.
# -> 혹은 seed가 자동적용되었을 가능성이 있다.

                  Stratified by treat
                   0             1             SMD   
  n                  185           185               
  age (mean (SD))  25.05 (9.51)  25.82 (7.16)   0.091
  educ (mean (SD)) 10.64 (2.30)  10.35 (2.01)   0.138
  married = 1 (%)     35 (18.9)     35 (18.9)  <0.001
  nodegree = 1 (%)   123 (66.5)    131 (70.8)   0.093

set.seed(2022)

mobj <- matchit(treat ~ age+educ+married+nodegree , # formula : 치료여부 ~ confoiunder들
                data = dat1, # dataset
                method = "nearest", # 매칭 방법: austin의 nearest ~
                distance = "glm", # 매칭 중 PS matching
                link = "logit", 
                replace = F, # 매칭된 환자 비복원
                capliper = 0.2, std.caliper = T, # caliper 지정
                m.order= "random",#매칭순서 -> default방법 안좋음. random으로 지정
               )

m.data <- match.data(mobj)

tobj <- CreateTableOne(vars=c("age", "educ", "married", "nodegree"),
                      data = m.data, strata = "treat")

print(tobj, smd = T, test=F)

# SMD가 다 0.2 이하라서 2군간의 confounder차이가 사라지도록 매칭된 데이터가 완성되었다.

                  Stratified by treat
                   0             1             SMD   
  n                  185           185               
  age (mean (SD))  25.12 (9.68)  25.82 (7.16)   0.081
  educ (mean (SD)) 10.61 (2.31)  10.35 (2.01)   0.120
  married = 1 (%)     35 (18.9)     35 (18.9)  <0.001
  nodegree = 1 (%)   124 (67.0)    131 (70.8)   0.082

1:2 매칭 by ratio=2

set.seed(2022)

mobj <- matchit(treat ~ age+educ+married+nodegree , # formula : 치료여부 ~ confoiunder들
                data = dat1, # dataset
                method = "nearest", # 매칭 방법: austin의 nearest ~
                distance = "glm", # 매칭 중 PS matching
                link = "logit", 
                replace = F, # 매칭된 환자 비복원
                capliper = 0.2, std.caliper = T, # caliper 지정
                m.order= "random",#매칭순서 -> default방법 안좋음. random으로 지정
                ratio=2, # 치료군(1) 한명당 대조군 2명을 붙인다.
               )

m.data <- match.data(mobj)

tobj <- CreateTableOne(vars=c("age", "educ", "married", "nodegree"),
                      data = m.data, strata = "treat")

print(tobj, smd = T, test=F)

# 만약, 185명의 2배가 아니라 더 적게 매칭되었다면?
# -> caliper가 있어서, PS가 비슷한 환자를 찾지 못해서 1명만 매칭된 경우가 포함되어있다.

                  Stratified by treat
                   0             1             SMD   
  n                  370           185               
  age (mean (SD))  27.88 (11.19) 25.82 (7.16)   0.220
  educ (mean (SD)) 10.35 (2.56)  10.35 (2.01)   0.002
  married = 1 (%)    161 (43.5)     35 (18.9)   0.551
  nodegree = 1 (%)   240 (64.9)    131 (70.8)   0.128

매칭된 데이터에서 치료효과(y) 평가 by GEE

• GEE를 통한 회귀분석 : geepack
• 회귀계수를 OR로 바꿔주는 패키지 : broom

library(geepack)
library(broom)

0,1범주형을 -> 0,1 연속형으로 만드는 sense ( bool -> as.double )

# -> 결과변수를 다시 연속형 변수로 매핑해줘야한다.

head(m.data$y == 1)

<ol class=list-inline>

TRUE

FALSE

TRUE

TRUE

FALSE

FALSE

</ol>

head(as.double(m.data$y == 1))

<ol class=list-inline>

1

0

1

1

0

0

</ol>

최종 회귀분석에서는 결과변수(질병여부, 연속형) ~ 설명변수(치료여부), id=매칭id변수(subclass) 2+1 개의 변수만

• match된 데이터 = 이미 2군간(treat) confounder 제거된, RCT 데이터 같은 데이터
  • 효과를 제거한 변수들 4개를 넣은 doubly robust 방법을 쓸 수 도 있지만, 여기선 simple하게 univarable model을 쓴다.

obj <- geeglm(as.double(y == 1) ~ treat, 
           data = m.data, # 
           id = subclass, # subclass칼럼에 정의된 각 매칭id를 제시하여 correlation이 있다는 것을 알려준다.
           family = "binomial", # 로지스틱 회귀분석 명시
           corstr = "exchangeable", # correlation structrue 암기...
           std.err = "san.se", # default라 안써도 되는 에러 계산시, robust estimator 사용
          )

# 회귀분석 결과는 기본적으로 summary에 넣으면 된다.
summary(obj)


#  Coefficients:
#             Estimate  Std.err  Wald Pr(>|W|)
# (Intercept) -0.01093  0.10426 0.011    0.917
# treat1      -0.17337  0.18076 0.920    0.337


# treat1에 있는 Estimate -> 회귀계수 -0.17337 과 p-value 0.337
# -> 강의랑 달라서 아래 스샷을 바탕으로...

Call:
geeglm(formula = as.double(y == 1) ~ treat, family = "binomial", 
    data = m.data, id = subclass, corstr = "exchangeable", std.err = "san.se")

 Coefficients:
            Estimate Std.err Wald Pr(>|W|)
(Intercept)  -0.0109  0.1043 0.01     0.92
treat1       -0.1734  0.1808 0.92     0.34

Correlation structure = exchangeable 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)        1  0.0047
  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha     1.01   0.104
Number of clusters:   553  Maximum cluster size: 2 

회귀분석 결과obj로 OR를 구하기 by broom패키지의 tidy()

# treat1의 p-value가 0.05보다 크다  == treat1이라는 변수는 회귀변수 결과변수에 별로 안중요하다?!
# -> treat여부가 질병발생(outcome발생) 비율에 유의미한 영향을 주지 않는다?!

# 주의) 여기서의 p-value는 회귀계수 자체의 p-value이지, OR이 아니다.
# -> 우리가 해석할 것은 treat의 회귀변수의 p-value가 아니라 OR를 구해야한다.

# broom패키지의 tidy()함수로 회귀분석결과 -> OR를 바로 구할 수 있다.
tidy(obj, exponentiate = T, conf.int = T)

A tibble: 2 × 7

term	estimate	std.error	statistic	p.value	conf.low	conf.high
<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
(Intercept)	0.989	0.104	0.011	0.917	0.806	1.21
treat1	0.841	0.181	0.920	0.337	0.590	1.20

# conf.low ~ conf.high는 OR의 95%신뢰구간 하한~상한이다.

정리하기

• Propensity score analysis (PSA)가 어떻게 confounder들을 없앤 데이터를 만드는지 알아봤다.

• Propensity score analysis (PSA)는 confounding이 존재하는 관찰연구의 데이터를 RCT의 데이터처럼 confounding이 없도록 가공하여 분석하는 방법이다.
• Propensity score (PS)는 치료군에 속할 확률을 의미한다.
• PS는 보통 로지스틱 회귀분석으로 추정하며, PS 모델에는 outcome에 영향을 끼칠 것으로 생각되는 변수를 모두 넣는 것이 좋다. 단, 치료군/대조군 결정이후에 관측된 변수는 넣으면 안된다.
• PS 추정 이후 분석 방법에는 matching, stratification, weighting, covariate adjustment 등이 있다.
• PS matching 후 치료군과 대조군의 balance는 standardized mean difference (SMD)로 평가할 수 있다.
• PS matching된 데이터는 matching을 고려하여 분석한다.

연습문제

01 다음 중 Propensity score analysis에 대한 설명으로 틀린 것은?

1. 관찰연구 자료의 confounding을 없애는 방법 중 하나이다.

2. Propensity score를 추정하는 모델은 설명변수 하나당 샘플사이즈 15가 필요하다.

3. 최종결과변수(환자의 outcome)에 영향을 끼치는 변수들은 PS 추정 모델에 넣는 것이 좋다.

4. 치료군/대조군이 결정된 이후에 관측된 변수들은 PS 추정 모델에 넣지 않는다.

• 정답 : 2
• 해설 : PS 추정 모델은 overfitting의 문제에 구애받지 않고 많은 수의 설명변수를 포함할 수 있다.

02 다음 중 PS(propensity score)를 추정한 후 자료 분석 방법으로 알맞지 않은 것은?

1. Matching

2. Stratification

3. Weighting

4. Randomization

• 정답 : 4
• 해설 : PS 추정 후 자료분석 방법에는 matching, stratification, weighting, covariate adjustment가 있다.