The Korean Journal of Appled Statstcs (2013) 26(5), 697 712 DOI: http://dx.do.org/10.5351/kjas.2013.26.5.697 Comparson of GEE Estmaton Methods for Repeated Bnary Data wth Tme-Varyng Covarates on Dfferent Mssng Mechansms Boram Park a Inkyung Jung b,1 a Bometrc Research Branch, Natonal Cancer Center b Department of Bostatstcs, Yonse Unversty College of Medcne (Receved October 25, 2012; Revsed July 31, 2013; Accepted September 6, 2013) Abstract When analyzng repeated bnary data, the generalzed estmatng equatons(gee) approach produces consstent estmates for regresson parameters even f an ncorrect workng correlaton matrx s used. However, tme-varyng covarates experence larger changes n coeffcents than tme-nvarant covarates across varous workng correlaton structures for fnte samples. In addton, the GEE approach may gve based estmates under mssng at random(mar). Weghted estmatng equatons and multple mputaton methods have been proposed to reduce bases n parameter estmates under MAR. Ths artcle studes f the two methods produce robust estmates across varous workng correlaton structures for longtudnal bnary data wth tme-varyng covarates under dfferent mssng mechansms. Through smulaton, we observe that tme-varyng covarates have greater dfferences n parameter estmates across dfferent workng correlaton structures than tme-nvarant covarates. The multple mputaton method produces more robust estmates under any workng correlaton structure and smaller bases compared to the other two methods. Keywords: Generalzed estmatng equatons, multple mputaton, weghted estmatng equatons, MCAR, MAR. 1. 서론다시점자료 (longtudnal data) 는시간에따라같은개체내에서반복측정된자료로관측값들사이에종속성이존재한다. 이러한관측값들사이의상관관계를고려하기위해일반화추정방정식 (generalzed estmatng equatons; GEE) 이많이이용되고있다. 일반화추정방정식은가상관행렬 (workng correlaton matrx) 을잘못가정하더라도모수의일치추정량 (consstent estmator) 을구할수있다 (Lang과 Zeger, 1986). 하지만, 일반화추정방정식은결측체계가완전임의결측 (MCAR) 이아닌경우에편의추정량을제공하고 (Troxel 등, 1997), 시간-종속적공변량 (tme-varyng covarate) 이포함된경우에는가상관행렬에따라회귀계수추정값이다르게나올수있다 (Lang과 Zeger, 1986). 결측체계가임의 1 Correspondng author: Assstant Professor, Department of Bostatstcs, Yonse Unversty College of Medcne, 50 Yonse-ro, Seodaemun-gu, Seoul 120-752, Korea. E-mal: jung@yuhs.ac
698 Boram Park, Inkyung Jung 결측 (MAR) 인경우에발생하는문제를해결하기위해가중방법 (Robns 등, 1995) 과다중대체 (Rubn, 1987) 방법을사용하는것이제안되었다. 본논문에서는모의실험을통하여자료의상관구조, 결측체계를변화시키면서세가지방법을적용하였을때시간-독립적공변량과시간-종속적공변량추정값이가상관행렬에따라어떤양상으로변화하는지살펴보고자한다. 또한, 가상관행렬에따른회귀계수추정값간에차이를통해각방법의로버스트성 (robustness) 을살펴보고그추정값의편차, 분산, 평균제곱오차를추정하여각방법별로정확성을비교하고자한다. 2절에서는일반화추정방정식, 결측체계, 가중방법, 그리고다중대체방법에대해소개하고, 3절에서는간질자료 (epleptc data) 에세가지방법을적용한결과를비교한다. 4절에서는자료의상관구조와결측체계를다양하게변화시켜세가지방법을적용한모의실험에대해설명하고, 5절에서는연구의결과를요약, 정리한다. 2. 이론적배경 2.1. 일반화추정방정식 시간에따라반복측정된자료는한개체에서의관측값들사이에종속성이존재한다. 이러한관측값들 사이의상관관계를고려하기위해 Lang 과 Zeger (1986) 가제안한일반화추정방정식 (GEE) 을이용하 여모형의모수를추정한다. 번째개체 ( = 1,..., K) 의 t 번째시간 (t = 1,..., n ) 에서반응변수값을 n 1 벡터인 Y = (y 1, y 2,..., y n ) T 라고하고, 설명변수값을 n p 행렬인 X = (x 1,..., x n ) T 이라하자. 지수 족분포인 y t 의주변밀도함수를식 (2.1) 과같이정의한다. f(y t ) = exp [{y t θ t a(θ t ) + b(y t )} ϕ], (2.1) 여기서 a( ), b( ) 는연결함수 (lnk functon), ϕ 는척도모수 (scale parameter), θ t = h(η t ), η t = x t β 이 고, β = (β 1, β 2,..., β p) T 은알려지지않은모수인 p 1 벡터이다. y t 의 1 차, 2 차적률을통해평균과 분산은각각다음과같다. E(y t ) = a (θ t ), var(y t ) = a (θ t ). (2.2) ϕ 반복측정된관측값의결합분포는구체화하지않고 y t 의주변분포를정의하고반응변수간의상관관계 를나타내는상관행렬을가정한후모수를추정한다. 이때, 가정하는상관행렬을가상관행렬 (workng correlaton matrx) 이라고부른다. 만약하나의개체로부터반복측정된관측값이서로독립이라면추정방정식은식 (2.3) 과같다 (Lang 과 Zeger, 1986). U I (β) = K X T S = 0, (2.3) =1 여기서 = dag(dθ t/dη t) 인 n n 행렬이고, S = Y a (θ) 인 n 1 벡터이다. D = d{a (θ)}/dβ라고하면, X T = D T V 1 이된다. 반복측정된관측값이서로독립이라고가정하였으므로 V 는단지 β의함수일뿐관측값사이의관계는고려되지않았다. 하지만반복측정된관측값사이에상관관계가존재하면그상관관계를나타내는가상관행렬을 n n 대각행렬인 R(α) 라고할때, y t 의분산 - 공분산행렬 V 는다음과같다. V = A 1 2 R(α)A 1 2 ϕ, (2.4)
Comparson of GEE Estmaton Methods 699 여기서 A = dag{a (θ t )} 인 n n 대각행렬, α 는가상관행렬인 R(α) 를설명하는 s 1 벡터이고만 약 R(α) 가 Y 의참상관행렬이라면 V 는 cov(y ) 와동일하다. 반복측정된관측값사이에존재하는상 관성을고려한분산 - 공분산행렬을이용해일반화추정방정식을표현하면식 (2.5) 와같이정의한다. K =1 D T V 1 S = 0, (2.5) 여기서 D = d{a (θ)}/dβ = A X 이고, = dag(dθ t/dη t), S = Y a (θ) 이다. 이때, 만 약 R(α) = I, 즉관측값이서로독립이라면식 (2.5) 는식 (2.3) 과동일하게된다. 그리고 U (β, α) = D T V 1 S 는유사우도함수와비슷하지만, 여기서는관측값사이에상관성을고려하므로 V 는 β의함수일뿐만아니라 α의함수이기도하다. β와 ϕ를알고있을경우, 식 (2.4) 의 α를 α의 K 1/2 -일치추정량 인 ˆα(Y, β, ϕ) 로대체함으로써식 (2.5) 를 β 만의함수로다시표현할수있다. β 를알고있을경우, ϕ 를 ϕ 의 K 1/2 - 일치추정량인 ˆϕ(Y, β) 로대체함으로써식 (2.5) 를다음과같은방정식의형태로정의할수있 다. K =1 { U [β, ˆα β, ˆϕ(β) }] = 0. (2.6) 식 (2.6) 을만족하는해를일반화추정방정식을이용한추정량 ˆβ GEE 라고하고, 이추정량은 y t 가정규 분포를따를때최대우도추정량과일치하게된다. 일반화추정방정식방법은가상관행렬을잘못가정하였을지라도모형의모수와그추정량의분산을일치 적으로 (consstently) 추정한다. 하지만, 모형에시간 - 종속적공변량이포함되면가상관행렬선택에따 라회귀계수추정값이다르게추정되므로, 가상관행렬의선택이중요해진다 (Wall 등, 2005). 2.2. 결측체계 다시점자료는같은개체내에서시간에따라반복측정된자료로결측 (mssng) 이많이발생한다. 자 료에서결측이발생하는원리에따라서결측체계 (mssng mechansms) 를크게완전임의결측 (mssng completely at random; MCAR), 임의결측 (mssng at random; MAR), 비임의결측 (not mssng at random; NMAR) 으로구분한다 (Lttle 과 Rubn, 2002). 결측체계에따라서분석방법이달라질수 있으므로결측체계를정확히파악하는것은자료를분석하는데매우중요한의미가있다. 번째개체의 t 시점까지반응변수를 Y = (Y 1, Y 2,..., Y t ) 라고할때, Y = (Y O Y O, Y M ) 는관측된자료로구성되어있다. 또한, 반응지시변수를 R = (R 1, R 2,..., R t) 라고 와관측되지않은자료 Y M 할때, 번째개체가 t시점에서관측되면 R t = 1이고, 결측이면 R t = 0으로정의한다. X 는모두관 측된설명변수라고가정하면각결측체계는아래와같은식으로표현할수있다. 완전임의결측은 R 가 Y O, Y M 적, Y M 로식 (2.9) 와같이표현할수있다. Pr (R Y, X ) = Pr (R X ), (2.7) ( ) Pr (R Y, X ) = Pr R Y O, X, (2.8) ( ) Pr (R Y, X ) = Pr R Y O, Y M, X. (2.9) 모두와독립적이므로식 (2.7) 과같고, 임의결측은 R 가 Y O 에는종속뿐만아니라 Y M 에도종속적이므 에는독립적이므로식 (2.8) 과같고, 비임의결측은 R 가 Y O
700 Boram Park, Inkyung Jung 2.3. 가중방법 일반화추정방정식방법은자료의결측체계가 MCAR 인경우모수에대한좋은추정량을제공하지만 MCAR 이아닌결측체계인경우에는편의추정량을제공한다. 이러한문제를해결하기위해자료의결 측체계가 MAR 일지라도불편추정량을제공해주는가중방법 (weghted estmatng equatons) 을사 용하는것이제안되었다 (Robns 등, 1995). 가중방법은 번째개체가 t 시점에서관측될확률의역비율값을관측값에가중을주는것이다. 번째 개체가 t 시점에서관측되었다면 R t = 1 이고, 그외에는 R t = 0 으로정의한다. 첫번째시점에서는항 상관측이되고, 한번결측이발생하면그뒤시점부터끝시점까지모두결측이발생한것으로가정한 다. 즉, 모든개체에대해서 R 1 = 1 이고, 만약 R t = 0 이면 R (t+1) = 0 임을의미한다. 시점 t 1 까지관측된자료중시점 t 에서결측이발생할확률은현재및미래관측치인 {Y t,..., Y T } 와 는상관없고, 과거관측치인 D t = {X, Y 0, Y 1,..., Y (t 1) } 에따라달라진다. 번째개체가 t 1 시 점에서관측되었다는조건하에 t 시점에서도관측될확률은 λ t = P (R t = 1 R (t 1) = 1, D t, Y ) = P (R t = 1 R (t 1) = 1, D t) 식을만족한다고가정한다. 그러나 MAR 체계에서는 λ t = P (R t = 1 R (t 1) = 1, D (T +1) ) = P (R t = 1 R (t 1) = 1, D t) 식을만족한다고가정한다 (Robns 등, 1994). 첫번째시점에서는항상관측된다고가정하므로모든 에대해서 λ 1 = 1 이다. 식 (2.10) 을통해가중 방법에서의 ˆβ 을구할수있다. S(β) = K =1 ( ) 1 µ W ϕa 1 2 β R A 1 2 (y µ ) = 0, (2.10) 여기서 W = dag{r 1 w 1,..., R T w T } 이고, w 1 t = ˆλ 1 ˆλ t 이다. 번째개체의 t시점에서가중값인 w t 는시점 t에서관측된절대적인확률의역수이고그추정량은조건부확률의누적곱의역수 로써위와같이구해진다. 이가중값 w t 와결측의유무를나타내는 R t 의곱을대각원소로갖는대각 행렬 W 가일반화추정방정식에추가된다. 시점 t 에서관측될확률이낮은관측값은큰가중이가해지 고, 관측될확률이높은관측값은낮은가중이가해진다 (Km, 2004). 이처럼 W 을통해결측처리를 함으로써가중방법은 MAR 가정하에서도유효한추정량을제공한다. 2.4. 다중대체 결측값에통계적모형을통하여어떤다른값으로채우는것을대체방법 (mputaton method) 이라고한 다. 결측값에하나의값으로채우는방법을단일대체 (sngle mputaton) 라고하는데, 이방법은관 측된값과대체된값을구분할수없어정보의양을과다추정하고, 추정량의분산을과소추정하는문 제가발생한다. 이러한문제점을해결하기위해결측값에여러개의값으로대체하고이값들간차이 의분산이추정량의분산을계산할때추가되어분산이과소추정되지않도록하는다중대체 (multple mputaton) 방법이 Rubn (1987) 에의해제안되었다. 다중대체를 m 번시행한후대체된자료각각에대하여 m 번분석하여얻어진모수의추정값들을 ˆβ 1, ˆβ 2,..., ˆβ m 라하고, 이모수들의추정된분산을 V 1, V 2,..., V m 이라하면 m 개의통합된모수의 추정값 ( ˆβ ) 과그모수의분산추정값 (V ) 을식 (2.11) 과같이정의한다. ˆβ = 1 m m =1 ( ) ˆβ, V m + 1 = W m + B m, (2.11) m 여기서통합된모수의분산추정값 V 은식 (2.12) 의 W m 과 B m 두개의분산성분을종합하여계산된
Comparson of GEE Estmaton Methods 701 다. W m = 1 m m =1 V, B m = 1 m 1 m ( ˆβ ˆβ ) 2. (2.12) 대체된자료들로부터추정된 m 개모수의분산추정값들평균인대체내분산 (wthn- mputaton varance; W m ) 과 m 개모수의추정값들사이의분산인대체간분산 (between-mputaton varance; B m ) 을 종합하여통합된모수의분산추정값을구한다 (Rubn, 1987). 한개의결측값에무한개의값으로대체한다면추정량의분산이정확하게추정될것이다. 하지만무한 개의값으로대체하는것은불가능하므로유한개의값으로대체하게된다. 이때비록대체시행횟수가 작더라도결측으로인해손실된모수에대한정보량이아주크지않다면다중대체를통해모수의분산이 거의비슷하게추정된다. 결측으로인해손실된모수에대한정보량이란결측이없는완전한자료와비 교했을때결측으로인해서발생한모수의정밀도 (precson) 의감소분을의미한다 (Song 과 An, 2009). =1 3. 실제자료를이용한분석 새로개발한항발작제 (ant-epleptc drug; AED) 의효능을살펴보기위해수집한총 89 명의간질환자 자료 (Faught 등, 1996) 를 2 절에서소개한세가지방법으로분석해보고자한다. 모든환자는약을복 용하기전필요한안정기간으로 12 주동안은약을복용하지않고매주병원을방문하여지난일주일 동안의발작횟수를측정하였고, 13 주부터 28 주까지 45 명은위약을복용하고 44 명은신약을복용하면 서매주발작횟수를측정하였다. 전체 89 명은최소 2 번에서최대 27 번씩관측이되었고, 총관측치는 1,419 로한사람당평균적으로 16 번씩관측되었다. 아래와같은모형을자료에적용하고자한다. ( X t = 1, x trt logt [P (Y t = 1)] = X tβ, = 1,..., 89, t = 1,..., 28 (3.1) ), x sex, x race, x age, x weght, x drug, β = (β0, β 1, β 2, β 3, β 4, β 5, β drug ). t Y t 는지난일주일동안의발작여부로발작이일어났으면 1, 일어나지않았으면 0 이다. X t 는시 간 - 독립적공변량 (tme-nvarant covarate) 인치료그룹 ( 위약, 신약 ), 성별, 인종, 나이, 몸무게, 그리 고시간 - 종속적공변량 (tme-varyng covarate) 인약복용여부가포함된다. 여기서, 약복용여부는 x drug t = I{t 13} 인지시변수이다. 시간 - 독립적공변량인치료그룹의회귀계수는위약군에비해신약 군이발작을일으킬위험을나타내고, 시간 - 종속적공변량인약복용여부의회귀계수는 t 시점에서약미 복용군에비해약복용군이발작을일으킬위험을의미한다. 일주일동안발작이일어났는지를반응변수로두고일반화추정방정식방법, 가중방법, 그리고다중대체 방법을적용한결과를비교하였다. 실제자료의결측은단조패턴 (monotone pattern) 으로다중대체는 식 (3.1) 로지스틱회귀모형을기반으로 10 번시행하였다. 각변수의회귀계수추정값 ( ˆβ w ), 로버스트 표준오차 ( V s ), 모형에근거한표준오차 ( V m ) 를이용하여가상관행렬에따라달라지는회귀계수추정값간의차이 (measure of dfference) 를 4 w=1 ( ˆβ w ˆβ) 2 /V s ( ˆβ w ) 에근거하여 Table 3.1 에나타내었다. 여기서, 회귀계수추정값간의차이란가상관행렬에따라도출된회귀계수추정값이얼마나다른지 정도를나타내며, 그값이클수록가상관행렬에따라회귀계수추정값의차이는크다는것을의미한다. 가상관행렬은독립적인 (ndependent) 구조, 교환가능한 (exchangeable) 구조, 자기상관 (autoregressve; AR(1)) 구조, 2- 종속적 (2-dependent; Toep(2)) 구조 4 가지를고려하였다. 회귀계수추정값간의차이 가시간 - 독립적공변량변수 5 개가비슷한패턴을보여치료그룹 (treatment) 변수의결과만표에제시 하였다.
702 Boram Park, Inkyung Jung Table 3.1. Coeffcent estmates wth robust and model-based standard errors ( V s and V m) for the epleptc data wth four dfferent workng correlaton matrces Measure of Method Parameter Ind Exch AR(1) Toep(2) Dfference ˆβ trt 0.443 0.366 0.440 0.427 Treatment Vs 0.279 0.275 0.277 0.277 0.051 Vm 0.117 0.285 0.155 0.166 GEE ˆβ drug 0.577 0.414 0.504 0.437 Drug Vs 0.128 0.117 0.124 0.130 1.041 Vm 0.123 0.106 0.151 0.155 WGEE MI Treatment Drug Treatment Drug ˆβ trt 0.497 0.429 0.485 0.465 Vs 0.261 0.294 0.267 0.270 Vm 0.026 0.056 0.032 0.034 0.034 ˆβ drug 0.702 0.476 0.632 0.587 Vs 0.136 0.124 0.136 0.138 Vm 0.028 0.026 0.034 0.035 1.619 ˆβ trt 0.438 0.446 0.438 0.437 Vs 0.195 0.192 0.195 0.195 Vm 0.128 0.203 0.140 0.147 0.001 ˆβ drug 0.563 0.563 0.545 0.524 Vs 0.159 0.159 0.159 0.159 Vm 0.139 0.135 0.149 0.154 0.042 세가지방법모두에서가상관행렬선택에따라달라지는회귀계수추정값간의차이는시간-독립적공변량 (treatment) 에비해시간-종속적공변량 (drug) 에서상대적으로더크게나타났다. 일반화추정방정식방법, 가중방법, 다중대체방법의시간-독립적공변량 (Treatment) 에서차이는각각 0.051, 0.034, 0.001인데반해, 시간-종속적공변량 (Drug) 에서차이는각각 1.041, 1.619, 0.042로더크게나타났다. 다중대체방법이가중방법보다일반화추정방정식방법의회귀계수추정값과더유사하였고, 가상관행렬에따른회귀계수추정값간의차이도더작게나타났다. 가상관행렬에따른회귀계수추정값의차이는시간-독립적공변량은일반화추정방정식방법에서, 시간-종속적공변량은가중방법에서가장크게나타났다. 간질자료의결측체계를정확하게알고있지않으므로더다양한결측체계에따라각각의방법을적용했을때가상관행렬에따른회귀계수추정값간의차이가나타나는패턴을살펴보고자 4절에서모의실험을시행하였다. 4. 모의실험 4.1. 자료에대한개요 4.1.1. 자료생성간질환자자료를근거로 100 개의완전한자료를재구성하였다. 공변량인성별, 인 종, 약복용시작시점은균등분포 (Unform dstrbuton) 를따르는난수발생을통하여생성하고, 나 이, 몸무게는정규분포 (Normal dstrbuton) 를따르는난수발생을통하여생성하였다. 전체대상자는 위약군 100명, 신약군 100명으로총 200명이각각 10주씩관측된다고가정한다. 약을복용하기시작한시점 t drug 은 4 6주로사람마다약복용시작시점을다르게설정하고, 시간-종속적공변량인약복용 여부는 x drug t = I{t t drug } 인지시변수이다. 각공변량을생성한후상관성이존재하는이항반응변수를생성하기위해서다변량이항분포 (mult-
Comparson of GEE Estmaton Methods 703 Table 4.1. Four dfferent mssng mechansms determned by α = (α 0, α 1, α 2, α 3) Mssng mechansm α = (α 0, α 1, α 2, α 3 ) MCAR α = (α 0, 0, 0, 0) MAR-weak α = (α 0, 0.2, 0, 0) MAR-strong α = (α 0, 0.5, 0, 0) MAR-2-dep α = (α 0, 0.5, 0.2, 0) varate bnary dstrbuton) 에근거한아래와같은방법을이용하였다 (Presser 등, 2002). 이방법은 n 1 평균벡터 π, n n 공분산행렬 V, 그리고 n n 상관행렬 C 에대한가정이필요하다. 먼 저, 평균벡터는간질환자실제자료를이용하여일반화추정방정식방법을적용했을때도출된로짓모형 에근거하여다음과같은로짓모형 logt[p (Y t = 1)] = 1.508 0.543x trt 0.006x age + 0.036x weght 0.577x drug t 0.295x sex + 0.199x race + 을이용해서구하였다. 공분산행렬은 V = A C A 으로여기 서 A = dag{v 1/2 t }, v t = π t (1 π t ) 이다. 상관행렬 C 은자기상관 (autoregressve; AR(1)) 행렬 (ρ = 0.4, 0.6) 과교환가능한 (exchangeable) 행렬 (ρ = 0.2, 0.4, 0.6) 로총 5 가지형태를가정하여살펴 보았다. 위와같이평균벡터, 공분산행렬, 상관행렬을가정하고 Z t = (Y 1,..., Y t 1 ), µ t = E(Z t ), G t = cov(z t ), s t = cov(z t, Y t ), b t = G 1 t s t (t = 2,..., T ) 가주어졌을때, 조건부평균 ν t 을식 (4.1) 과같이정의한다. ν t = ν t (z t ; π, V ) := P (Y t = 1 Z t = z t ) = π t + b t (z t µ t ) (4.1) t 1 = π t + b tj (y j π j ) (t = 2,..., T ). j=1 발작여부를나타내는이항반응변수 Y = (Y 1,..., Y 10 ) 에서 Y 1 은평균 π 1 을갖는베르누이분 포 (Bernoull dstrbuton) 를따르는난수발생을통하여생성하고, Y t (t = 2,..., 10) 는조건부평 균 ν t 을갖는베르누이분포 (Bernoull dstrbuton) 를따르는난수발생을통하여생성하였다. 이처럼 첫시점의반응변수는평균벡터를이용하고나머지시점에서의반응변수는조건부평균을이용하여조 건부선형성질 (condtonal lnear property) 을갖는다변량이항분포에근거하여상관성이존재하는 이항반응변수를생성하였다. 4.1.2. 결측생성본논문에서는공변량은모두관측되었다는가정하에반복측정된반응변수에서의결측만고려하여모의실험을시행하였다. 이항반응변수에결측을생성하기위하여각개체의각시점에서관측될확률 (λ t ) 을구한다. 첫번째시점에서는항상관측이되고, 한번결측이발생하면그뒤시점부터끝시점까지모두결측이발생한것으로가정한다. 즉, 모든개체에대해서 R 1 = 1이고, 만약 R t = 0이면 R,t+k = 0, k > 0임을의미한다. 아래와같은로짓모형을이용하여각시점에서관측될확률 (λ t) 을구하였다 (Presser 등, 2002). logt(λ t) = α 0 + α 1y (t 1) + α 2y(t 2)I(t > 2) + α 3yt, t = 2,..., 10, (4.2) 여기서 yt 는 yt = 2y t 1로 번째개체가 t시점에서발작이일어났으면 yt = 1, 발작이일어나지않았으면 yt = 1이다. 그리고두번째시점이후이면 I(t > 2) = 1, 나머지는 I(t > 2) = 0이며, 처음시점에서는반드시관측된다는가정을통해항상 λ 1 = 1이다. α 1, α 2, α 3 의값을 Table 4.1과같이결측체계에따라각각다르게설정하였다. MCAR의경우에는현재시점에서관측될확률은이전어느시점의관측값에도영향을받지않는다. 그에반해결측체계
704 Boram Park, Inkyung Jung Fgure 4.1. Measure of dfference of parameter estmates for tme-nvarant and tme-varyng covarates across the dfferent workng correlaton matrces (Autoregressve) 가 MAR일때는현재시점에서관측될확률이이전시점의관측값에영향을받는다. MAR-weak과 MAR-strong은한시점전의관측값에따라관측될확률이달라지며, MAR-2-dep는두시점전과한시점전두개의관측값에따라관측될확률이달라진다. 식 (4.2) 모형에서이전시점까지관측되었다는조건하에현재시점에서결측이될평균조건부확률 (average condtonal probablty) 과관련된 α 0 을결측률에따라다르게지정한다. 즉, 결측률이낮을때는 α 0 에큰값을결측률이높을때는작은값을지정하는데, 결측률이 5%, 10%, 20%, 40% 일때각각의 α 0 값을 3.0, 2.2, 1.4, 0.4로지정한다 (Presser 등, 2002). 본논문에서는모든결측체계에서결측률을 5% 라고가정하여 α 0 = 3.0으로지정하고모의실험을시행하였다. 4.2. 모의실험결과 반응변수의다양한상관구조와여러가지결측체계에서의일반화추정방정식방법 (GEE), 가중방 법 (WGEE), 다중대체방법 (MI) 의로버스트성 (robustness) 을살펴보고자한다. 각각의결측체계별 로세가지방법을적용하여가상관행렬에따른회귀계수추정값간의차이를비교하였다. 여기서, 회귀계수추정값간의차이 (measure of dfference) 는 100 번의모의실험에서 4 w=1 ( ˆβ w ˆβ) 2 /V s ( ˆβ w ) 에근거하여도출된 100 개의추정값을평균내어제시하였고, 회귀계수추정값간의차이가시간 - 독립적공변
Comparson of GEE Estmaton Methods 705 Fgure 4.2. Measure of dfference of parameter estmates for tme-nvarant and tme-varyng covarates across the dfferent workng correlaton matrces (Exchangeable) 량변수 5 개가비슷한패턴을보여치료그룹 (treatment) 변수의결과만그림에나타내었다. Fgure 4.1, Fgure 4.2 를보면결측체계, 적용방법에상관없이시간 - 종속적공변량 (drug) 은시간 - 독 립적공변량 (treatment) 에비해가상관행렬에따른회귀계수추정값의차이가항상더크게나타났다.
706 Boram Park, Inkyung Jung Fgure 4.3. Bas of parameter estmates for tme-nvarant and tme-varyng covarates (Autoregressve) 자료의상관행렬이자기상관 (AR(1)) 구조일때보다교환가능한 (exchangeable) 구조일때, 그리고동일한상관행렬내에서는 ρ값이증가할수록시간-종속적공변량과시간-독립적공변량간에차이가더뚜렷하게나타나는것을확인할수있다. 시간-독립적공변량은회귀계수추정값의차이가자료의상관행렬, 결측체계, 적용하는방법에크게의존하지않지만, 시간-종속적공변량은결측체계에따라각각의방법별로가상관행렬에따른회귀계수추정값의차이가다른경향을보였다. 결측체계가 MCAR에서 2-dep MAR로갈수록가상관행렬에따른회귀계수추정값의차이가가중방법과다중대체방법에서는크게변화가없었지만, 일반화추정방정식을적용했을때는점점증가하였다. 즉, 일반화추정방정식방법은결측발생이이전관측값과상관성이강할수록가상관행렬에따른회귀계수추정값의차이가증가하는경향이있었다. 또한, 동일한결측체계에서가중방법이일반화추정방정식방법보다시간-종속적공변량의회귀계수추정값차이를더감소시켰다. 다중대체방법은전반적으로가상관행렬에따른회귀계수추정값의차이가작아가상관행렬의형태에로버스트함을확인할수있었다. 각각의결측체계별로세가지방법을적용했을때추정의정확성을살펴보고자참값과각방법을통해구한추정값을이용해서편차 (bas), 분산 (varance), 평균제곱오차 (mean squared error; MSE) 를구해각각의경우를비교하였다. 이때, 각각의결측체계별로세가지각방법내에서가상관행렬에따른편차및평균제곱오차가거의비슷한결과가도출되어서사전에가정한상관행렬인자기상관 (AR(1)) 구
Comparson of GEE Estmaton Methods 707 Fgure 4.4. Bas of parameter estmates for tme-nvarant and tme-varyng covarates (Exchangeable) 조와교환가능한 (exchangeable) 구조만그림에제시하였다. Fgure 4.3, Fgure 4.4 를보면결측체계와적용한방법에따라편차는다른경향을보였다. 시간 - 독 립적공변량인경우에가중방법은일반화추정방정식방법이나다중대체방법보다편차가더작게나타
708 Boram Park, Inkyung Jung Fgure 4.5. Mean squared error (MSE) of parameter estmates for tme-nvarant and tme-varyng covarates (Autoregressve) 났고, 일반화추정방정식방법과다중대체방법은비슷한경향을보였다. 시간-종속적공변량인경우에 MCAR이나 Weak MAR 가정하에서는가중방법이편차가가장크게나타났지만, Strong MAR 또는 2-dep MAR 가정으로갈수록다중대체방법의편차가점점증가하는경향을보였다. 즉, 결측발생이이전관측값과상관성이강할수록다중대체방법은추정값의정확성이떨어지는것을확인할수있다. 위와같은경향은자료의상관행렬이자기상관구조일때보다교환가능한구조일때, 그리고동일한상관행렬내에서는 ρ값이증가할수록더뚜렷한형태를보이고있다. 또한, 시간-종속적공변량이시간-독립적공변량에비해편차의변화폭이더큰경향을보였다. Fgure 4.5, Fgure 4.6을보면결측체계나자료의상관구조에상관없이전반적으로일반화추정방정식방법과다중대체방법보다가중방법을통해구한추정값의평균제곱오차가더크게나타났다. 또한다중대체방법이일반화추정방정식방법에비해평균제곱오차가시간-독립적공변량에서는더작게추정되었지만, 시간-종속적공변량에서는더크게추정되었다. 이는자료의상관행렬이자기상관 (AR(1)) 구조일때보다교환가능한 (exchangeable) 구조일때, 그리고동일한상관행렬내에서는 ρ값이증가할수록더뚜렷하게나타났다. 분산은전반적으로평균제곱오차와비슷한패턴을보였다. 다른방법에비해가중방법을적용했을때분산이크게추정되었고, 이는가중방법이평균제곱오차가크게추정되는데
Comparson of GEE Estmaton Methods 709 Fgure 4.6. Mean squared error (MSE) of parameter estmates for tme-nvarant and tme-varyng covarates (Exchangeable) 영향을미쳤다. 전반적으로가상관행렬이정확하게가정되었을때잘못가정된상관행렬에비해편차가 더낮게추정되었고, 로버스트표준오차와모형에근거한표준오차간에차이도더작게나타났다.
710 Boram Park, Inkyung Jung 5. 결론다시점자료연구에서일반화추정방정식방법은가상관행렬을잘못가정하더라도모수에대한일치추정량을구할수있어서많이쓰이고있다. 하지만결측체계가완전임의결측이아닌경우에편의추정량을제공하고, 시간-종속적공변량이포함된경우에는가상관행렬에따라회귀계수추정값이다르게도출될수있다. 본논문에서는일반화추정방정식방법, 가중방법, 다중대체방법을이용하여 GEE 분석에서시간-독립적공변량과시간-종속적공변량의추정값이가상관행렬에따라어떤양상으로변화하는지연구하였다. 반응변수의다양한상관구조와여러가지결측체계에서각방법의로버스트성 (robustness) 을살펴보고, 참값과각방법을통해구한추정값을이용하여정확성 (accuracy) 을비교하였다. 정확성의척도로는편차 (bas), 분산 (varance), 평균제곱오차 (MSE) 를살펴보았고편차와분산의정보를모두이용한평균제곱오차에근거하여정확성을비교하였다. 시간-독립적공변량은반응변수의상관구조, 결측체계, 적용하는방법에크게의존하지않으며가상관행렬의형태에로버스트한반면, 시간-종속적공변량은가상관행렬에따라회귀계수추정값이다르게도출되었다. 시간-종속적공변량에일반화추정방정식방법을적용하면결측체계가이전시점의관측값과상관성이강할수록가상관행렬에따른회귀계수추정값의차이가크게나타났다. 하지만다중대체방법을적용하면가상관행렬에따른회귀계수추정값의차이가작아전반적으로가상관행렬의형태에로버스트하였다. 시간-독립적공변량은가중방법을적용하였을때편차는가장작게분산은가장크게추정되었고, 시간-종속적공변량인경우에는다중대체방법을적용하였을때편차는가장크게분산은가장작게추정되었다. 그래서편차와분산의정보를모두이용한평균제곱오차에근거하여세가지방법의정확성을비교하였다. 그결과, 전반적으로일반화추정방정식방법과다중대체방법에비해가중방법을적용하였을때평균제곱오차가가장크게추정되었다. 시간-독립적공변량은다중대체방법에서평균제곱오차가가장작게추정되었고, 시간-종속적공변량은일반화추정방정식방법에서평균제곱오차가가장작게추정되었다. 본논문에서반응변수의상관구조, 결측체계를변화시키면서세가지방법을모두적용해보았다. 일반화추정방정식방법은평균제곱오차가작아정확성은입증하였지만가상관행렬에따른회귀계수추정값의차이가크므로가상관행렬의선택이중요해진다. 가중방법은가상관행렬의형태에크게의존하지는않지만평균제곱오차가크게추정되어정확성이떨어지는경향이잇다. 다중대체방법은가상관행렬에따른회귀계수추정값의차이도가장작고, 평균제곱오차도작게추정되어다른방법에비해더좋은추정량을제공해주는것을확인할수있었다. 본논문에서는반응변수에만결측이있는경우에대해살펴보았지만, 실제자료에서는공변량에도결측이많이발생한다. 이러한자료에서는가중방법보다더정확한추정을하는다중대체방법을선호한다 (Beunckens 등, 2008). 또한시간-독립적공변량에비해시간-종속적공변량은가정된가상관행렬이모수의추정에영향을미친다. 모형에시간-종속적공변량이포함되어있을때, 비대각 (non-dagonal) 가상관행렬을사용하면 GEE와주변모형의가정에어긋나므로편의추정량을도출할수있다 (Pepe와 Anderson, 1994). 반면에, 독립적구조인가상관행렬을사용하면시간-종속적공변량의회귀계수추정에효율 (effcency) 이떨어진다 (Ftzmaurce, 1995). 가상관행렬에따라회귀계수추정값이달라지는시간-종속적공변량은가상관행렬을선택하는데있어서구체적인방법이필요하리라여겨진다. References Beunckens, C., Sotto, C. and Molenberghs, G. (2008). A smulaton study comparng weghted estmatng equatons wth multple mputaton based estmatng equatons for longtudnal bnary data, Compu-
Comparson of GEE Estmaton Methods 711 tatonal Statstcs & Data Analyss, 52, 1533 1548. Faught, E., Wlder, B. J., Ramsay, R. E., Refe, R. A., Kramer, L. D., Pledger, G. W. and Karm, R. M. (1996). Topramate placebo-controlled dose-rangng tral n refractory partal eplepsy usng 200-, 400-, and 600-mg daly dosages, Neurology, 46, 1684 1690. Ftzmaurce, G. M. (1995). A caveat concernng ndependence estmaton equatons wth multple multvarate bnary data, Bometrcs, 51, 309 317. Km, T. H. (2004). Handlng data n GEE wth mssng response, Sungkyunkwan Unversty. Lang, K. Y. and Zeger, S. L. (1986). Longtudnal data analyss usng generalzed lnear models, Bometrka, 73, 13 22. Lttle, R. J. A. and Rubn, D. B. (2002). Statstcal Analyss wth Mssng Data, John Wley & Sons. Pepe, M. S. and Anderson, G. (1994). A cautonary note on nference for margnal regresson models wth longtudnal data and general correlated response data, Communcaton n Statstcs B, 23, 939 951. Presser, J. S., Lohman, K. K. and Rathouz, P. J. (2002). Performance of weghted estmatng equatons for longtudnal bnary data wth drop-outs mssng at random, Statstcs n Medcne, 21, 3035 3054. Robns, J. M., Rotntzky, A. and Zhao, L. P. (1994). Estmaton of regresson coeffcents when some regressors are not always observed, Journal of the Amercan Statstcal Assocaton, 189, 846 866. Robns, J. M., Rotntzky, A. and Zhao, L. P. (1995). Analyss of semparametrc regresson models for repeated outcomes n the presence of mssng data, Journal of the Amercan Statstcal Assocaton, 90, 106 121. Rubn, D. B. (1987). Multple Imputaton for Nonresponse n Surveys, John Wley & Sons. Song, J. W. and An, H. (2009). Handlng and Analyss of Mssng Data, Statstcal Tranng Insttute, Seoul. Troxel, A. B., Lpstz, S. R. and Brennan, T. A. (1997). Weghted estmatng equatons wth nongnorably mssng response data, Bometrcs, 53, 857 869. Wall, M. M., Da, Y. and Eberly, L. E. (2005). GEE estmaton of a msspecfed tme-varyng covarate: An example wth the effect of alcoholsm treatment on medcal utlzaton, Statstcs n Medcne, 24, 925 939.
712 Boram Park, Inkyung Jung 시간 - 종속적공변량이포함된이분형 반복측정자료의 GEE 를이용한분석에서결측 체계에따른회귀계수추정방법비교 박보람 a 정인경 b,1 a 국립암센터바이오메트릭연구과, b 연세대학교의학통계학과 (2012 년 10 월 25 일접수, 2013 년 7 월 31 일수정, 2013 년 9 월 6 일채택 ) 요약다시점자료연구에서일반화추정방정식은가상관행렬을잘못가정하더라도모수의일치추정량을도출하므로많이이용된다. 하지만, 결측체계가완전임의결측이아닌경우에는편의추정량을제공하고, 시간 - 종속적공변량이포함된경우에는가상관행렬에따라회귀계수추정값이다르게도출될수있는문제점이있다. 결측체계가임의결측인경우에발생하는문제를해결하기위해가중방법과다중대체방법을사용하는것이제안되었다. 본논문에서는시간 - 종속적공변량이포함된이분형반복측정자료를 GEE 를이용하여분석할때다양한결측체계에서일반화추정방정식방법, 가중방법, 다중대체방법의회귀계수추정에대한로버스트성과정확성을모의실험을통하여비교해보았다. 세가지방법모두에서시간 - 종속적공변량의회귀계수가시간 - 독립적공변량의회귀계수에비해가상관행렬에따라추정값의차이가크게나타났다. 다른두방법에비해다중대체방법이가상관행렬의형태에대해더로버스트하고편의도작은추정치를도출하였다. 주요용어 : 일반화추정방정식, 다중대체, 가중방법, 완전임의결측, 임의결측. 1 교신저자 : (120-752) 서울시서대문구연세로 50, 연세대학교의학통계학과, 조교수. E-mal: jung@yuhs.ac