BMC 의학 연구 방법론 저널을 훑어보다 '로지스틱 회귀를 이용하여 비교 위험도를 추정하는 간단 기법'이라는 논문(원문 링크: http://www.biomedcentral.com/1471-2288/12/14)을 읽게 됐다. 흔한 결과 변수인 경우 오즈비(OR)가 과다 추정되므로 위험비(RR)를 산출해야 되는데 다변수 모형에서 값을 얻기는 생각보다 쉽지 않다. 논문의 방법보다 서론에서 언급된 대로 (1) binreg 명령어를 이용하여 binomial regression을 적용하거나(참고문헌 4), (2) (1)의 방법의 경우 데이터가 크면 수렴이 안되는 경우가 흔하므로 Cox regression을 응용하면 된다(참고문헌 6).

* Hosmer & Lemeshow의 low birth weight data 불러오기.

. webuse lbw

(Hosmer & Lemeshow data)


* oddsrisk 명령어를 이용하여 단변수 분석으로 OR과 RR을 산출해 비교하기.

. oddsrisk low smoke

---------------------------------------------------------------------

Incidence for unexposed risk group =     0.2522

---------------------------------------------------------------------

Predictor    Odds Ratio   Risk Ratio     [95% Conf. Interval]

---------------------------------------------------------------------

smoke           2.0219       1.6076       1.0591       2.2230

---------------------------------------------------------------------


* Stata에서 binomial regression을 이용하여 RR 구하기

. binreg low smoke, rr nolog

Generalized linear models                   No. of obs      =       189

Optimization     : MQL Fisher scoring       Residual df     =       187

                   (IRLS EIM)               Scale parameter =         1

Deviance         =  229.8045995             (1/df) Deviance =  1.228902

Pearson          =  188.9999911             (1/df) Pearson  =  1.010695


Variance function: V(u) = u*(1-u)           [Bernoulli]

Link function    : g(u) = ln(u)             [Log]


                                            BIC             = -750.4021


-----------------------------------------------------------------------

      |                 EIM

  low | Risk Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]

------+----------------------------------------------------------------

smoke |   1.607642   .3433245     2.22   0.026     1.057812    2.443262

 cons |   .2521739    .040495    -8.58   0.000     .1840819    .3454532

-----------------------------------------------------------------------


* Stata에서 Cox regression을 응용하여 RR 구하기

. gen time=1

. stset time, fail(low==1)

     failure event:  low == 1
obs. time interval:  (0, time]
 exit on or before:  failure

-----------------------------------------------------------------------
      189  total obs.
        0  exclusions
-----------------------------------------------------------------------
      189  obs. remaining, representing
       59  failures in single record/single failure data
      189  total analysis time at risk, at risk from t =         0
                             earliest observed entry t =         0
                                  last observed exit t =         1

. stcox smoke, vce(robust) nolog

         failure _d:  low == 1
   analysis time _t:  time

Cox regression -- Breslow method for ties

No. of subjects      =          189         Number of obs   =       189
No. of failures      =           59
Time at risk         =          189
                                            Wald chi2(1)    =      4.92
Log pseudolikelihood =   -307.61219         Prob > chi2     =    0.0266

-----------------------------------------------------------------------
      |               Robust
   _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
------+----------------------------------------------------------------
smoke |   1.607642   .3442364     2.22   0.027     1.056637     2.44598
-----------------------------------------------------------------------

다변수 모형에서 적용을 위해서는 설명 변수에 원하는 교란 변수를 포함시켜주기만 하면 된다.

Posted by cyberdoc
:

연속형 변수의 일치도 평가에 이용되는 Bland-Altman plot batplot 명령어http://goo.gl/zzZTI ) 이용해 그려보자.

. ssc install batplot
. sysuse auto, clear
. batplot mpg turn, title(Agreement between mpg and turn) ///
  info valabel(make) notrend xlab(26(4)38) moptions(mlabp(9))

 
Posted by cyberdoc
:
관련성을 평가할 때 많이 쓰는 교차비(odds ratio)는 엄밀하게는 위험비(risk ratio)와 다르다. 교차비와 위험비가 비슷하려면 (1) 관심 질환의 유병률이 낮거나(희귀 질환 가정), (2) 잘 설계된 환자-대조군 연구에서 산출됐을 경우이다. 실제 데이터에서는 이 두 가지 가정 어느 쪽도 만족하지 않는데도 교차비를 제시한 경우가 흔하다. 의학보건학 연구에서 교차비가 많이 쓰이는 이유는 범주형 결과 변수를 분석할 때 로지스틱 회귀분석을 주로 이용하고, 그 결과로 교차비를 쉽게 산출할 수 있기 때문이다. 위험비는 푸아송(Poisson)[각주:1] 회귀분석을 적용해서 계산할 수 있는데 많은 연구자들이 익숙치 않다. Joseph HilbeLogistic Regression Models 책 에서 제시한 oddsrisk 명령어로 교차비와 위험비를 둘다 산출해보고 비교해보도록 하자. (명령어는 findit oddsrisk로 찾아서 설치할 수 있고, 데이터셋은 링크한 CRC Press 책 소개 사이트에서 다운로드할 수 있다.)

. use anterior

. list

     +--------------------------+
     | count   death   anterior |
     |--------------------------|
  1. |   120       1          1 |
  2. |    67       1          0 |
  3. |  2005       0          1 |
  4. |  2504       0          0 |
     +--------------------------+

. oddsrisk death anterior [fw=count]

---------------------------------------------------------------------
Incidence for unexposed risk group =     0.0261
---------------------------------------------------------------------
Predictor    Odds Ratio   Risk Ratio     [95% Conf. Interval]
---------------------------------------------------------------------
anterior        2.2368       2.1670       1.6220       2.8807
---------------------------------------------------------------------

비노출군에서 발생 위험이 2.6%=67/(67+2504)로 낮으므로 교차비와 위험비가 비슷할 것이라는 예상을 해볼 수 있고, 실제 결과도 교차비가 2.24, 위험비가 2.17로 큰 차이가 없다.

poisson 명령어를 이용하여 푸아송 회귀분석을 수행하여 산출한 위험비와 비교해보자.

. poisson death anterior [fw=count], nolog irr robust

Poisson regression                                Number of obs   =       4696
                                                  Wald chi2(1)    =      26.69
                                                  Prob > chi2     =     0.0000
Log pseudolikelihood =  -776.2572                 Pseudo R2       =     0.0171

------------------------------------------------------------------------------
             |               Robust
       death |        IRR   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    anterior |   2.166953   .3243484     5.17   0.000     1.616003    2.905741
       _cons |   .0260599   .0031423   -30.25   0.000     .0205748    .0330073
------------------------------------------------------------------------------

계산된 위험비가 2.17로 같음을 알 수 있다.

  1. 지금까지 통계학 관련 서적에서는 포아송으로 번역해 왔으나 국립국어원 외국 인명 표기에 따르면 '푸아송'이 맞다. [본문으로]
Posted by cyberdoc
:
다변수 로지스틱 분석에서 범주형 설명 변수의 빈도표는 unitab 명령어를 써보자. OR과 95% CI는 logistic 명령어와 같으므로 비교해보자. (unitab 명령어는 사용자 작성 명령어므로 찾아서 설치해야 동작한다.)

. webuse lowbirth
 

. unitab low age lwt, c(ptd smoke ht ui)
------------------------------------------------------------------------------
     low    |    low=1(%)      Total(%)     OR   [95% Conf. Interval]  p-value
------------+-----------------------------------------------------------------
     age    |      56(50)       112(100)   1.000    0.917     1.090     1.000
------------+-----------------------------------------------------------------
     lwt    |      56(50)       112(100)   0.989    0.976     1.002     0.093
------------+-----------------------------------------------------------------
     ptd  0 |      38(44)        87(78)     .        .         .        0.013
          1 |      18(72)        25(22)    3.316    1.256     8.750
------------+-----------------------------------------------------------------
   smoke  0 |      26(39)        66(59)     .        .         .        0.007
          1 |      30(65)        46(41)    2.885    1.319     6.307
------------+-----------------------------------------------------------------
      ht  0 |      49(48)       102(91)     .        .         .        0.185
          1 |       7(70)        10( 9)    2.524    0.618    10.308
------------+-----------------------------------------------------------------
      ui  0 |      42(46)        92(82)     .        .         .        0.048
          1 |      14(70)        20(18)    2.778    0.981     7.864
------------------------------------------------------------------------------

unitab 명령어 바로 뒤는 결과 변수이고 연속형 변수는 결과 변수 뒤에, 범주형 변수는 옵션으로 c( )에 넣으면 끝. 연속형 변수의 low=1 (%)는 low 변수가 1로 코딩된 빈도와 백분율이다.

'통계 연습' 카테고리의 다른 글

[Stata] batplot 명령어  (0) 2011.12.26
[Stata] oddsrisk 명령어  (0) 2011.10.20
[Stata] 표준화비 신뢰구간 계산  (0) 2011.10.14
[Stata] 위험 차(risk difference) 계산  (2) 2011.10.14
[Stata] optmatch2 명령어  (0) 2011.10.14
Posted by cyberdoc
:
표준화비(standardized mortality/morbidity ratio, 이하 SMR)는 두 인구 집단의 율을 간접 표준화(indirect standardization) 기법을 이용하여 비교하는 지표이다.

Stata에서 여러가지 방법으로 SMR과 95% 신뢰구간을 계산할 수 있는데, 로데이터가 없이 값만 얻어졌을
경우 신뢰구간을 계산하려면 사용자 작성 명령어인 eclpci 명령어를 쓰면 쉽다.

어떤 인구집단에서 관찰 사망자수가 40명이고, 기대 사망자수가 30명인 경우,

. eclpci 40 30, format(%3.2f)
SMR: 1.33  95% CI[0.95, 1.82]

SMR: 1.33이므로 연구 집단은 (기대 사망자수를 산출한 인구 집단에 비해) 사망이 33% 높지만, (95% 신뢰구간이 1을 포함하고 있으므로) 5% 유의수준에서 통계적으로 뚜렷한 차이는 없다.

. eclpci 20 50, format(%3.2f)
SMR: 0.40  95% CI[0.24, 0.62]

SMR: 0.40이므로 연구 집단은 (기대 사망자수를 산출한 인구 집단에 비해) 사망이 60% 낮고, 5% 유의수준에서 통계적으로 뚜렷한 차이가 있다.

Posted by cyberdoc
:
단면조사(crosssectional survey) 자료나 코호트 연구(cohort study) 자료에서 위험 차(risk difference, 이하 RD)를 구하는 Stata 기본 명령어는 cs, csi이다. 이 명령어를 수행하면 노출과 비노출군의 RD와 위험 비(risk ratio, 이하 RR), 그리고 각각의 95% 신뢰구간과 P-value를 계산할 수 있다.

. cs case exposed 
. cs case exposed, or
. csi a b N1 N2

RD의 95% 신뢰구간을 다양한 방법으로 계산하는 명령어는 rdci, rdcii이다. 이는 사용자 작성 명령어로 findit rdci 명령문을 내린 다음 찾아서 설치해야 한다.

. findit rdci
. rdci case exposed
. rdcii a b N1 N2

두 명령어 모두 노출 변수의 범주가 예/아니오로 되어 있는 경우만 가능하다. 만일 세 단계 이상의 노출 변수에 대한 RD를 산출하려면 csjl 명령어를 이용하면 된다. 이 또한 사용자 작성 명령어이므로 findit csjl 명령문으로 찾아서 설치해야 한다. csjl 명령어는 type(rr rd) 옵션을 주면 RR과 RD를 동시에 출력해주고, at 옵션을 주면 attributable risk도 제시해준다. 또한 by(변수) 옵션을 주면 층화분석을 수행하고 각 층에 대한 RD도 산출해준다. 각 rd 값의 통계적 검정과 P-value를 구하려면 chi나 exact 옵션을 추가한다.

. findit csjl
. csjl case exposed, rd
. csjl case exposed, type(rr rd)
. csjl case exposed, rd table chi exact
. csjl case exposed, rd by(sex)


다변수 모형을 통해 보정된(adjusted) RD와 95% 신뢰구간을 계산하려면, binreg 명령문에 rd 옵션을 사용한다.

. xi: binreg case i.exposed, rd nolog
. xi: binreg case i.exposed agegrp sex, rd nolog

Posted by cyberdoc
:
성향 점수(propensity score)를 만드는 데는 psmatch2 명령어가 가장 강력하지만, 만든 다음 짝지은 데이터셋을 분리해내려면 optmatch2 명령어가 편하다. 이 명령어는 isvar 명령어와 함께 동작하므로 미리 설치해두어야 한다.

치료군(treatm)을 성향 점수(propen)에 맞춰 1:1로 데이터셋을 추출하기 위한 명령어는 다음과 같다.

. optmatch2 treatm propen, minc(1) maxc(1) gen(mid)

----+--- 1 ---+--- 2  ---+--- 3 ---+--- 4 ---+--- 5
...........................

Matched 27 treated subjects and 27 untreated subjects
Sum of distances within matched sets = 0.0265667

. sort mid treatm

. list id treatm propen mid in 1/10

     +--------------------------------+
     |  id   treatm      propen   mid |
     |--------------------------------|
  1. | 251        0   .11095736     1 |
  2. |  55        1   .11095736     1 |
  3. |  89        0   .05489074     2 |
  4. |  86        1   .05489074     2 |
  5. | 136        0   .20299611     3 |
     |--------------------------------|
  6. |  36        1   .20299611     3 |
  7. |  58        0   .19385163     4 |
  8. | 297        1   .19615819     4 |
  9. |  52        0   .08239614     5 |
 10. | 186        1   .08174527     5 |
     +--------------------------------+

. drop if mid==.

. save propennew.dta 

치료군과 비치료군에서 성향 점수가 가장 비슷한 27명을 골랐고, 짝지은 쌍에 대해 새로운 아이디(mid)를 만들었다. 이제 치료군과 비치료군에서 원하는 결과 변수에 대해 통계적 검정을 수행하면 된다.

Posted by cyberdoc
: