투고일：2016. 8. 5         심사일：2016. 9. 8         게재확정일：2016. 9. 19

치의학 분야에서 SPSS를 이용한
일반화 추정방정식의 단계별 안내

1)전남대학교 치의학전문대학원 치과교정학교실, 치의학 연구소, 2)전남대학교 자연과학대학 통계학과
임 회 정1), 박 수 현2)

ABSTRACT

A step-by-step guide to Generalized Estimating Equations
using SPSS in dental research

Department of Orthodontics, Chonnam National University School of Dentistry1)
Dental Science Research Institute, Chonnam National University2)
Department of Statistics, Chonnam National University3)
Hoi-Jeong Lim1, 2), Su-Hyeon Park3)

The Generalized Estimating Equations (GEE) approach is a widely used statistical method for analyzing longitudinal data and clustered data in clinical studies. In dentistry, due to multiple outcomes obtained from one patient, the outcomes produced from an individual patient are correlated with one another. This study  focused on the basic ideas of GEE and introduced the types of covariance matrix and working correlation matrix. The quasi-likelihood information criterion (QIC) and quasi-likelihood information criterion approximation (QICu) were used to select the best working correlation matrix and the best fitting model for the correlated outcomes. The purpose of this study is to show a detailed process for the GEE analysis using SPSS software along with an orthodontic miniscrew example, and to help understand how to use GEE analysis in dental research.

Key words : Generalized Estimating Equations; clustered data; dental research

Corresponding Author

Hoi-Jeong Lim, PhD
Department of Orthodontics, Chonnam National University School of Dentistry, Dental Science Research Institute
33 Yongbong-ro, Buk-gu, Gwangju 500-757, South Korea
Telephone : +82-62-530-5830, Fax : +82-62-530-5659, E-mail : hjlim@jnu.ac.kr

이 논문은 2015년도 전남대학교 학술연구비 지원 (#2014-2174) 에 의하여 연구되었음.

Ⅰ. 서론

통계 방법은 결과변수의 척도에 의해 결정되는데, 결과변수가 연속형인 경우 한 개 이상의 설명변수들을 가지고 결과변수를 예측하고자 할 때 쓰이는 통계 방법은 회귀분석이다. 그리고 사망/생존, 질병이 있다/없다 등의 범주형 자료가 결과변수가 되는 경우 이러한 형태로 이루어진 각각 독립인 일변량 결과변수를 분석하기 위해서 로지스틱 회귀분석 방법이 이용되어왔다. 그러나 결과변수가 동일한 피험자에 대하여 반복측정(longitudinal data) 되었거나 군집 표집(cluster sampling) 된 경우, 반복 측정된 결과변수들 간에, 혹은 군집된 자료로 이루어진 결과 변수들 간에 상관관계가 존재할 것이고 이를 통계적 분석에서 반드시 고려해야 한다. 치의학 분야에서는 주로 한 개체의 입안에서 여러 개의 치아를 대상으로 결과 값을 얻어내거나, 한 개체의 잇몸에 여러 개의 미니스크류의 식립 안정성을 알아보는 연구를 할 경우, 이는 군집된 자료(clustered data)로 여러 개의 미니스크류의 식립 안정성 간의 상관성이 존재한다. 따라서 이러한 경우에는 결과변수들 간의 상관관계를 모형에 포함시킬 수 없는 로지스틱 회귀분석보다는 가상관 행렬을 가정하여 상관성을 분석과정에서 고려해주는 일반화 추정방정식 분석방법(GEE, Generalized Estimating Equations) 을 적용해야 할 것이다.

이 연구의 목적은 일반화 추정방정식 분석을 수행하기 위한 절차를 자세히 알아보고 SPSS를 이용하여 치의학 연구에서 일반화 추정방정식 분석이 어떻게 사용되었는지 알아보고자 한다.

Ⅱ. GEE의 기본 개념

Liang and Zeger (1986)1), Zeger and Liang (1986)2)에 의해 개발된 종속변수들 간의 상관관계를 고려해 주는 일반화 추정방정식은 각각의 결과 값이 독립인 일반화 선형모형(GLM, Generalized Linear Model)3) 의 확장이며, 종속변수들 간의 상관관계를 나타내는 가상관 행렬(working correlation matrix)을 가정하여 반복 측정된 결과변수들 간의 상관성을 고려하여 모수를 추정하는 분석방법이다4).

개체 i 에서 시간이나 위치가 서로 다른 j 를 가진 yij는 범주형이거나 연속형인 결과변수이고, χ=(χij1, ···,χijp)’가 p×1 벡터로 구성된 범주형이거나 연속형인 설명변수인 데이터구조를 가질 때, 계수 β를 추정하는 일반화 추정방정식은 다음과 같다1, 2).

여기서 μi=(μi1, ···,μit)’은 평균벡터이며 Yi=(yi1,yi2 ···,yit)는 위에서 설명한 결과변수이고, Vi는 Yi의 공분산 행렬 (covariance matrix)의 추정량이며 다음과 같이 쓰여 진다.

위 공식에서 Ø는 overdispersion parameter를 말하고, A는 variance functions의 대각 행렬을 말하며, R(α)는 Y의 가상관 행렬을 말한다.

1. 공분산 행렬 (Covariance matrix)

GEE에서 공분산 행렬의 종류는 모형을 근거로 하는 공분산 행렬(model-based covariance matrix) 과 로버스트 공분산 행렬(robust covariance matrix) 이 있다. 로버스트 공분산 행렬은 경험적인(empirical) 혹은 샌드위치(sandwich) 공분산 행렬이라고도 부른다. 이 두 공분산 행렬의 차이점은 모형을 근거로 하는 공분산 행렬은 가상관 행렬이 바르게 지정되었을 때에만 모형의 모수에 대한 추정치와 이 추정치의 분산에 대한 추정치가 일관되고, 로버스트 공분산 행렬은 가상관 행렬이 바르게 지정되지 않을 때에도 Y에 대한 평균 μ에 대한 모형이 올바르다면 위의 추정치들이 항상 일관된다. 연구자가 가상관 행렬을 데이터에 맞게 지정해야 하는데 실제로 적합한 가상관 행렬을 고르는데 오류가 발생할 수 있으므로 가상관 행렬이 바르게 지정되지 않았다 할지라도 항상 일관된 로버스트 공분산 행렬의 추정치를 사용하는 것이 바람직하다5). GEE의 표준오차는 로버스트 공분산 행렬의 각각의 값들을 제곱근하여 얻어진 값들이다. 그러나 이 추정량은 표본크기가 적을 때나 공변량이 치우친 분포를 가질 때 바이어스 (측정값 또는 추정량의 분포 중심(평균값)과 참값과의 편차)6)가 개입될 수 있다4).

2. 가상관 행렬 (Working correlation matrix)

가상관 행렬은 각 개체에서 다른 시간에 얻어진 종속변수들 간의 상관성을 나타내는데, 모수를 추정하기 위하여 연구자가 데이터의 성질에 맞는 상관행렬을 지정해야 한다. 하지만 참 상관행렬을 알 수 없기 때문에 아래 표에 해당하는 가짜 상관행렬을 가정하고 이 상관행렬이 진짜 상관행렬인 것처럼 간주하고 모수를 추정하게 된다. 가상관 행렬이 바르게 지정되지 않았을 때, GEE 방법에 의한 모수 추정치의 상대적 효율 (relative efficiency) 이 낮아진다. 이러한 문제를 해결하기 위해 가상관 행렬을 선택하는 새로운 기준을 여러 학자들이 제안하였다. Pan(2001)7)은 QIC (Quasi-likelihood Information Criterion) 를, Rotnitzky and Jewell (1990)8)은 RJC (Rotnitzky and Jewell Criterion)를, Hin and Wang (2009)9)은 CIC (Correlation Inform ation Criterion)를, 그리고 Gosho et al. (2011)10)은 DEW(Discrepancy between the covariance matrix Estimator and the specified Working correlation matrix)를 제안하였다. 본 연구에서는 SPSS 결과에 포함된 QIC만을 다루었다. 가장 적합한 가상관 행렬을 선택하는 기준으로 다음 섹션에서 설명하고 있는 가장 낮은 값의 QIC를 사용할 수 있다. 가상관 행렬은 아래 Table 1과 같이 (1) Independent, (2) M-dependent, (3) Exchangeable, (4) Autoregressive (AR-1), (5) Unstructured, (6) Fixed 가상관 행렬로 분류되고 있다. (1)은 각 개체에서 다른 시간에 얻어진 종속변수들 간의 서로 독립임을 가정하는 것을 의미하며, (2)는 M번째까지만 상관관계가 성립되고, 그 이후에는 서로 독립임을 가정하는 상관행렬이다. 그리고 (3)은 각각의 개체로부터 다른 시간에 얻어진 종속변수들 간의 상관정도가 α로 항상 일정하다는 것을 가정한 것이며, (4)는 시간이 지남에 따라 상관의 정도가 감소함을 의미한다. (5)는 어떤 특정한 패턴을 가지고 있지 않고, 모든 종속변수들 간의 상관정도가 각기 다르다고 가정한 것이다. 그리고 (6)은 연구자 자신이 가상관 행렬을 지정하는 것을 말한다11).
기존에는 GEE 분석에서 독립적인 가상관 행렬을 주로 사용했지만, 가상관 행렬을 어떻게 선택할지에 대한 연구는 지금까지 계속 이어져 오고 있다. 많은 경우에 독립적인 가상관 행렬을 가정한 GEE 추정치들은 높은 효율성(high efficiency) 을 가진다고 보고되어 왔고, 개체 수가 증가하면 가상관 행렬이 바르게 지정되지 않은 경우에 대한 효율성(추정량의 분산이 작으면 높은 효율성을 가진다는 것을 의미하며, 분산이 크면 낮은 효율성을 가진다는 것을 의미한다. 효율적인 추정량은 최소분산불편추정량(Minimum Variance Unbiased Estimator)인 데, 이는 불편추정량 중에 분산이 가장 작은 추정량을 말한다.)12) 손실을 줄여준다고 보고되어 왔다13, 14). 그러나 시간에 따라 변하는 공변량(time-varying covariate) 을 포함하고 있을 때, 독립적인 가상관 행렬을 사용하는 것은 비효율적인 추정치를 만들어 낼 수 있다고도 보고하고 있다15).

3. QIC vs. QICu

GEE에서 종속변수들 간의 상관관계를 설명하기 위해 어떤 가상관 행렬을 사용할지에 대한 선택과 어떤 모형이 최종 모형으로 가장 잘 적합 되었는지에 대한 선택이 필요하다.
GLM은 독립된 결과 값들의 최대 우도법 (maximum likelihood method, 모집단의 모수(θ)의 값을 추정하기 위해 우도함수[L(θ)]를 최대로 하는 모수의 값을 구하는 방법)16)에 근거하고 있고17), GEE 방법은 유사 우도법(quasi-likelihood method, 관측치에 대한 분포의 가정이 없는 경우 우도(likelihood)를 구할 수 없는데, 이때 우도 대신 쓸 수 있는 함수가 유사우도 함수(quasi-likelihood function)이다. 유사우도법은 모집단의 모수(θ)의 값을 추정하기 위해 유사우도함수를 최대로 하는 모수의 값을 구하는 방법)12, 18, 19)에 근거하고 있을 뿐 아니라 결과 값들의 분포에 대해 어떤 가정도 하고 있지 않다. 그러므로 우도(likelihood) 이론을 근거로 유도된 통계량은 GEE에 직접 적용할 수 없다. 예를 들어 AIC (Akaike's information criterion) 는 GLM에서 모형 선택 시 널리 사용되는 방법이지만, GEE에선 적합하지 않다. 그리하여 Pan7)은 AIC 방법을 적절히 수정하여 GEE 분석에 적합한 모형 선택 방법인 QIC와 QICu를 만들어 냈다. 가장 낮은 QIC는 가장 적합한 가상관 행렬을 선택하는 기준이고 가장 낮은 QICu 는 가장 잘 적합된 모형(the best-fitting model) 을 선택하는 기준이 된다는 것을 보고하였다20). 본 연구에서 사용한 예제에서도 위의 QIC와 QICu를 사용하여 가장 적합한 가상관 행렬인 비구조적인 가상관 행렬을 선택하였고 나이와 위치(상악/하악), 식립 경험 등을 포함하는 최종 모형을 구축하였다.  

Ⅲ. 치의학 논문에서 사용된 데이터에 대한 정보

교정 치료 시 과거에는 헤드기어를 이용하여 교정력을 부여해왔는데, 이 방법은 환자의 협조가 반드시 필요하였다. 이러한 헤드기어의 단점을 보완한 다른 대안으로 미니 스크류가 개발되었고 이러한 미니 스크류는 환자의 협조도 불필요, 식립의 간편성, 저렴한 비용 등의 장점과 함께 교정력을 부여하는데 주로 사용하게 되었다. 그러나 미니 스크류의 크기가 작아 식립 안정성이 떨어져 이러한 식립 안정성에 영향을 미치는 요인을 알아보는 연구들이 많아졌다. 다음 데이터는 미니스크류 식립 실패(LOOSENING)에 영향을 미치는 요인을 조사하여 얻은 가상으로 만들어진 데이터이다21~23). 미니스크류 예제 데이터는 전남대학교 치의학전문대학원>일반대학원>자료실 21번에 올려두었고, 웹사이트 주소는 아래와 같다.

http://dent.jnu.ac.kr/user/indexSub.action?codyMenuSeq=6806&siteId=dent&menuUIType=top

ID = 환자 식별 번호
AGE = 환자 나이
GENDER = 환자 성별 (남자=1, 여자=0)
AREA1 = 미니 스크류 식립 부위 (상악=1, 하악=0)
AREA2 = 구체적인 미니 스크류 식립 부위
(between the mandibular 1st and 2nd molars=1,
between the maxillary 1st and 2nd molars=2,
between the mandibular 1st molar and 2nd premolar=3,
between the maxillary 1st molar and 2nd premolar=4,
between the mandibular 1st and 2nd premolars=5,
between the maxillary 1st and 2nd premolars=6,
between the mandibular 1st premolar and canine=7,
between the maxillary 1st premolar and canine=8,
between mandibular canine and lateral incisor=9,
between maxillary canine and lateral incisor=10,
between mandibular central incisors=11,
between maxillary central incisors=12)
LENGTH = 미니 스크류 길이 (6mm, 7mm, 8mm, 10mm, 12mm)
DIAMETER = 미니 스크류 직경 (1.6mm, 1.8mm)
NINSERT = 의사의 미니 스크류 식립 횟수
DOCTOR_ID = 의사 식별 번호
LOOSENING = 미니 스크류 식립 실패 여부 (실패=1, 성공=0)
AGE1 = 1 if AGE쨧20, AGE1 = 0 if AGE≤20
LENGTH1 = 0 if LENGTH쨦8mm, LEN GTH1 = 1 if LENGTH≥8mm
NINSERT1 = 1 if NINSERT쨧20, NIN SERT1 = 0 if NINSERT≤20

여기서 결과 변수는 LOOSENING이고,
설명 변수는 AGE, GENDER, AREA1, LENGTH, DIAMETER, NINSERT이다.

그 중 범주형 변수(factor)는 GENDER, AREA1, LENGTH, DIAMETER이고, 연속형 변수(covariate)는 AGE, NINSERT이다. 해석의 간편함을 위해 이러한 연속형 변수를 AGE1이나 NINSERT1의 범주형 변수로 바꾸었고 LENGTH도 8mm를 기준으로 이분형 변수 LENGTH1으로 바꾸었다.
이 데이터의 특징을 살펴보면, Table 2에서 보여주는 것처럼 환자 1명당 1개의 미니 스크류를 식립한 경우부터 2개, 3개, 4개, 5개, 그리고 6개까지의 미니 스크류를 식립한 경우가 존재했다. 환자 1명이 여러 개의 결과(미니 스크류의 안정성)를 갖는 하나의 군집(cluster)을 형성한다. 이러한 하나의 군집(여기서는 개체) 안에서 결과변수가 서로 상관관계가 높은지 알아본 결과, Table 3에서는 미니 스크류 한 개만을 식립했을 경우는 제외하고, 식립된 미니 스크류가 모두 성공한 경우(all success)는 59.6%였고, 모두 실패한 경우 (all failure)는 3.8%였다. 또, 거의 성공하거나 실패한 경우(skewed)는 18.3%였고, 성공과 실패의 비율이 같은 경우(equally distributed)도 18.3%였다. 이러한 하나의 개체 안에서 식립된 미니 스크류 하나가 성공하면 다른 미니 스크류도 높은 확률로 성공하였고, 하나가 실패하면 다른 미니 스크류도 높은 확률로 실패하였다. 이는 약 80% 이상이 개체 내에서 서로 상관관계가 높았음을 보여주는 결과이다.

Ⅳ. SPSS 절차 및 결과 설명

1. 단변수 GEE 분석을 위한 SPSS 절차

아래는 Figure 1의 GEE 단변수 분석을 위한 SPSS 절차에 대한 설명이다.
(1) 여기서 사용한 데이터는 아래와 같다.
(2) 분석 > 일반화 선형 모형 > 일반화 추정 방정식을 클릭한다.
(3) 개체 변수에 변수 ID를, 개체 내 변수에 AREA2를 옮기고, 상관행렬 작업의 구조에서 ‘비구조적’을 선택한 후, 메뉴의 모형 유형을 클릭한다.
(4) 사용자 정의에서 분포: 이항 검정, 연결함수: 로짓 로그선형분석을 선택 후, 메뉴의 응답을 클릭한다. 혹은 이분형 응답 또는 이벤트/시도 데이터의 이분형 로지스틱을 선택 후, 메뉴의 응답을 클릭해도 같은 결과를 얻을 수 있다.
(5) 종속변수에 LOOSENING을 옮기고, 참조 범주를 클릭, 처음(가장 낮은 값)을 선택. 계속을 누른 후, 메뉴의 예측자를 클릭한다.
(6) 요인에 AREA1을 옮기고, 옵션을 클릭, 요인에 대한 범주 순서에서 내림차순을 선택, 계속을 누른다.
(7) 메뉴의 모형을 선택, 요인 및 공변량에서 AREA1을 화살표를 이용하여 모형으로 옮긴다.
(8) 메뉴의 통계량을 선택, 인쇄에서 아래와 같이 통계량들을 선택한다.

위 과정에 대해서 더 구체적으로 설명하면, Figure 1(4)에서 결과변수가 연속형(normally-distributed response)일 때, 분포는 ‘정규’를, 연결함수는 ‘동일(I)’을 선택한다. 결과변수가 0/1로 구성된 경우라면 분포는 ‘이항 검정’을, 연결함수는 ‘로짓 로그선형분석’을 선택한다. 또한 결과변수가 단위시간 당 발생 건수인 경우 (예를 들어, 치과의 단위 시간 당 도착하는 환자의 수) 분포는 ‘포아송 분포(I)’를, 연결함수는 ‘로그’를 선택한다.

Figure 1(5)는 이항분포인 종속변수 (LOOSENING=0 혹은 1)의 참조 범주를 정의하는 것인데, 참조 범주 (Reference group)란 범주끼리 비교하고자 할 때 기준이 되는 범주를 말한다. 여기서 처음 (가장 낮은 값)은 LOOSENING=0 (성공)이 참조 범주이며, 마지막 (가장 높은 값)은 LOOSENING=1 (실패)이 참조 범주가 된다.

Figure 1(6)에서 요인 아래에 있는 옵션을 클릭하면 새로운 윈도우가 뜬다. 거기서 요인에 대한 범주 순서 (AREA1=0 혹은 1) 는 요인의 참조 범주를 정의하는 것이고 순서의 뒷부분이 참조 범주가 된다. 오름차순인 경우 0과 1중에 1이 순서의 뒷부분 이므로 AREA1=1 (상악)이 참조 범주가 되며, 내림차순인 경우 1과 0중에 0이 순서의 뒷부분 이므로 AREA1=0 (하악)이 참조 범주가 된다.

Figure 2는 AREA1 변수에 대한 비구조적 가상관 행렬을 가정했을 때의 단변수 분석 결과이다. Figure 2(1)에서 종속변수의 참조 변수가 LOOSENING=0일 때, OR=Exp(B)=0.54 이며 OR에 대한 95% 신뢰구간은 (0.30, 1.00)*임을 보였다. AREA1=1 (상악)이 참조 변수이므로 하악에 비해 상악에서 실패할 위험이 0.54배 더 낮다는 것을 의미하며, Figure 2(2)는 종속변수의 참조 변수를 LOOSENING=0에서 LOOSENING=1로 바꾸는 것이기 때문에 상악에서 성공할 위험이 1.84(=1/0.54)배 하악보다 더 높다는 것으로 바뀐다.

Figure 3은 독립 가상관 행렬을 사용했을 때, Figure 1(5)의 종속변수의 참조 범주와 Figure 1(6)의 요인에 대한 참조 범주에 따른 OR값의 변화를 살펴본 것이다.

1) Figure 1(5)의 종속변수 유형 (이항분포만 해당)에서 참조 범주는 처음 (가장 낮은 값)이고, Figure 1(6)의 요인에 대한 범주에서 내림차순을 선택하면 LOOSENING=0이 종속변수의 참조 범주가 되고, AREA=0이 요인의 참조 범주가 된다. 이 경우  이다.

2) Figure 1(5)의 종속변수 유형 (이항분포만 해당)에서 참조 범주는 처음 (가장 낮은 값)이고, Figure 1(6)의 요인에 대한 범주에서 오름차순을 선택하면 LOOSENING=0이 종속변수의 참조 범주가 되고, AREA=1이 요인의 참조 범주가 된다. 이 경우  이다.

3) Figure 1(5)의 종속변수 유형 (이항분포만 해당)에서 참조 범주는 마지막 (가장 높은 값)이고, Figure 1(6)의 요인에 대한 범주에서 내림차순을 선택하면 LOOSENING=1이 종속변수의 참조 범주가 되고, AREA=0이 요인의 참조 범주가 된다. 이 경우  이다.

4) Figure 1(5)의 종속변수 유형 (이항분포만 해당)에서 참조 범주는 마지막 (가장 높은 값)이고, Figure 1(6)의 요인에 대한 범주에서 오름차순을 선택하면 LOOSENING=1이 종속변수의 참조 범주가 되고, AREA=1이 요인의 참조 범주가 된다. 이 경우  이다.

odds는 식립이 실패할 확률이 성공할 확률의 비 (몇 배인지)로 정의된다. odds의 개념은 우리나라에서는 많이 사용되지 않는 것이라 이해하기 어렵다. 상악에서의 odds=42/184=1/4.38=0.2283의 의미는 상악에서 식립에 실패할 확률이 성공할 확률의 0.23배이다. 즉, 4.38(=1/0.2283)번 성공할 때마다 1번의 실패를 기대할 수 있다는 의미이다. 또한 하악에서의 odds=27/52=1/1.92=0.52 의 의미는 하악에서 식립에 실패할 확률이 성공할 확률의 0.52배이다. 즉, 1.92 (=1/0.52)번 성공할 때마다 1번의 실패를 기대할 수 있다는 의미이다. odds ratio는 이 두 odds를 나누어 주는 것이므로 (OR=0.23 /0.52=0.44) 상악에서 식립에 실패할 위험이 하악보다 0.44배 더 높다는 것을 의미하며, 이는 상악에서의 성공할 위험이 하악보다 2.27배 더 높다는 것을 의미한다.

여기서 우리는 종속변수의 참조 범주와 요인의 참조 범주에 따라 OR값이 달라지고, 결과 해석을 좀 더 이해하기 쉽게 할 수 있다는 것을 알았다.
Figure 4는 독립 가상관 행렬을 가정했을 때의 AREA1 변수에 대한 단변수 분석을 참조 범주를 0에서 1로 바꾸어 SPSS를 이용하여 얻은 결과이며, Figure 3과 Figure 4의 OR값이 0.44와 2.27로 일치하는 것을 알 수 있었다.

2. 다변수 GEE 분석을 위한 SPSS 절차

아래는 Figure 5에서 GEE 다변수 분석을 위한 SPSS 절차에 대한 설명이다.

(1) 요인에 AGE1, AREA1, NINSERT1을 옮기고, 옵션을 클릭, 요인에 대한 범주 순서에서 내림차순을 선택한다.
(2) 상단 메뉴의 모형을 선택, 요인 및 공변량에서 AGE1, AREA1, NINSERT1을 모형으로 옮긴다.
(3) 확인을 누른다.

Figure 6는 SPSS를 이용하여 다변수 분석을 한 결과이며, 구체적인 결과 해석은 Table 5 아래에 있다.

3. 가장 잘 적합된 최종 모형 찾기

Table 4는 가장 적합한 가상관 행렬을 선택하기 위해
STEP1: 모든 종류의 가상관 행렬에서 각각의 공변량, AGE1, GENDER, AREA1, LENGTH1, DIAMETER, NINSERT1에 대해 QIC 값을 구한다. AREA1에서 비구조적 가상관 행렬을 가졌을 때 다른 가상관 행렬을 가졌을 때보다 QIC=323.410으로 한 개의 공변량을 포함한 모형 중에서 가장 낮은 QIC값을 가졌다.

STEP2: 두 개의 공변량을 포함하는 모형에서 각각의 가상관 행렬에 대한 QIC 값을 구해본다. 두 개의 공변량을 포함하는 모형이란 (AGE1, GENDER), (AGE1, AREA1), (AGE1, LENGTH1), … , (DIAMETER, NINSERT1)으로 구성된다. 두 개의 공변량을 포함하는 모형에서 각각의 가상관 행렬에 대한 QIC 값을 구해본 결과, AREA1과 NINSERT1을 포함하는 모형에서 비구조적 가상관 행렬을 가졌을 때 QIC=318.344로 한 개의 공변량 AREA1을 포함한 모형보다 더 낮은 QIC 값을 가졌다.

STEP3: 세 개의 공변량을 포함하는 모형에서 각각의 가상관 행렬에 대한 QIC 값을 구해본다. 세 개의 공변량을 포함하는 모형이란 (AGE1, GENDER, AREA1), (AGE1, GENDER, LENGTH1), (AGE1, GENDER, DIAMETER), … , (LENGTH1, DIAMETER, NINSERT1)으로 구성된다. 세 개의 공변량 AREA1, NINSERT1, AGE1을 가진 모형에서 한 개의 공변량을 포함한 모형보다, 그리고 두 개의 공변량을 포함하는 모형보다 더 낮은 QIC 값(=315.080)을 가졌다.

STEP4: 네 개의 공변량을 포함하는 모형, 다섯 개의 공변량을 포함하는 모형, 그리고 6개의 공변량을 포함하는 모형에서 각각의 가상관 행렬에 대한 QIC 값을 구해본다. 그 결과 세 개의 공변량 AREA1, NINSERT1, AGE1을 가진 모형에서 가장 낮은 QIC 값(=315.080)을 가져 비구조적 가상관 행렬이 가장 적합한 것으로 선택되었다.

가장 적합한 모형 선택을 위해
STEP1: 위에서 선택된 비구조적 가상관 행렬을 이용하여 한 개의 공변량을 가진 모형에서 가장 낮은 QICu 값을 찾는다. 그 결과 AREA1을 포함한 모형이 선택되었다 (QICu = 323.109).

STEP2: 비구조적 가상관 행렬을 이용하여 두 개의 공변량을 가진 모형에서 가장 낮은 QICu 값을 찾는다. 그 결과 AREA1과 NINSERT1을 포함한 모형이 선택되었다 (QICu = 318.215).

STEP3: 비구조적 가상관 행렬을 이용하여 세 개의 공변량을 가진 모형에서 가장 낮은 QICu 값을 찾는다. 그 결과 AREA1, NINSERT1과 AGE1을 포함한 모형이 선택되었다 (QICu = 315.004).

STEP4: 비구조적 가상관 행렬을 이용하여 네 개, 다섯 개, 그리고 여섯 개의 공변량을 가진 모형에서 가장 낮은 QICu 값을 찾는다. 그 결과 세 개의 공변량 AREA1, NINSERT1, AGE1을 가진 모형에서 가장 낮은 QICu(=315.004) 값을 가져 최종 모형으로 선택되었다.

4. 단변수 및 다변수 GEE 분석 결과와 해석

Table 5에서 다변수 분석 결과를 살펴보기 위해 보정된 OR값을 살펴보자. 20세를 초과하는 환자들은 20세 이하의 환자들에 비해 식립 시 0.58배 덜 실패했지만, 유의한 차이는 나타나지 않았다. 상악은 하악에 비해 0.43배 유의하게 덜 실패했고, 20회 초과 식립 경험이 있는 경우가 없는 경우에 비해 0.48배 유의하게 덜 실패했다. 즉, 하악이고, 20회 기준 식립 경험이 없을수록 미니스크류 식립에 유의하게 더 실패하는 것으로 나타났다. Age가 최종 모형에 포함된 이유는 95% 신뢰구간에 1이 가까스로 포함되어 유의성이 나타나지 않았지만, 표본크기가 약간만 더 커진다면 유의성을 보일 가능성이 많기 때문으로 생각된다. 귀무가설이 식립경험이 미니스크류 식립 실패와 관련이 없는 것이라면, 20세 전후에 따른 식립 실패 차이와 상악과 하악의 차이가 식립 실패에 미치는 영향을 배제한 후에도 유의하게 20회 초과의 식립 경험이 미니스크류 식립 성공을 이끈다는 결론을 이끌어 낼 수 있었다.

Ⅴ. 결론

GEE 분석 연구는 한 환자 당 여러 개의 결과 관측치를 가지는 치의학 연구에서 반드시 필요한 통계 방법이라고 할 수 있다. 그러나 GEE 분석 방법은 통계를 전공하지 않은 치의학자들에게는 이해하기 쉽지 않은 고급 통계 방법에 속하기 때문에 최대한 연구자가 다가가기 쉽도록 꼭 필요한 절차와 이론만을 다루었다. 치의학 각 분야에서 SPSS를 이용한 GEE 분석을 위한 절차들에 대한 설명이 한 환자 당 여러 개의 결과를 가지는 군집된 자료(clustered data)를 분석하려는 연구자들에게 이러한 통계 방법을 알리고, 스스로 통계 분석을 수행하여 연구 결과를 얻고 해석하는데 도움이 되기를 바란다.

참고문헌
1. Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika 1986;73(1):13-22.
2. Zeger SL, Liang KY. Longitudinal data analysis for discrete and continuous outcomes. Biometrics 1986;42(1):121-130.
3. Nelder JA, Wedderburn RWM. Generalized linear models. J R Stat Soc Series A 1972;135(3):370-84.
4. 임회정. SAS를 이용한 고급의학 통계 및 예제. 1판 1쇄. 서울. 경문사. 2004.
5. Zeger SL, Liang KY, Albert PS. Models for Longitudinal Data: A Generalized Estimating Equation Approach. Biometrics 1988;44(4):1049-60.
6. 월간전자기술 편집위원회. 전자용어사전(E+). 개정판 3쇄. 성안당. 2007.
7. Pan W. Akaike's information criterion in generalized estimating equations. Biometrics 2001;57(1):120-5.
8. Rotnitzky A, Jewell NP. Hypothesis testing of regression parameters in semiparametric generalized linear models for cluster correlated data.  Biometrika 1990;77(3):485-97.
9. Hin LY, Wang YG. Working-correlation-structure identification in generalized estimating equations. Stat Med 2009;28(4):642-58.
10. Gosho M, Hamada C, Yoshimura I. Criterion for the selection of a working correlation structure in the generalized estimating equation approach for longitudinal balanced data. Commun Stat Theory Methods  2011;40(21):3839-56.
11. Hedeker D, Gibbons RD. Longitudinal Data Analysis. New Jersey. John Wiley & Sons. Inc. 2006.
12. 한국교육평가위원회. 교육평가 용어사전. 학지사. 2004.
13. Zeger SL. The analysis of discrete longitudinal data: Commentary. Statistics in Medicine 1988;7(1-2):161-8.
14. McDonald BW. Estimating logistic regression parameters for bivariate binary data. J R Stat Soc Series B 1993;55(2):391-7.
15. Fitzmaurice GM. A caveat concerning independence estimating equations with multiple multivariate binary data. Biometrics 1995;51(1):309-17.
16. 한국기상학회, 기상청. 대기과학용어사전(최신). 시그마프레스. 2015.
17. McCullagh P, Nelder JA. Generalized Linear Models. 2nd ed. London: Chapman & Hall. 1989.
18. Lindsay BG. Composite Likelihood Methods. Contemporary Mathematics. 1988;80(1):221-39.
19. Davidson R, MacKinnon J. Econometric Theory and Methods. New York, NY: Oxford University Press. 2004.
20. Cui J. QIC program and model selection in GEE analyses. Stata Journal 2007;7(2):209-20.
21. Hong SB, Kusnoto B, Kim EJ, BeGole EA, Hwang HS, Lim HJ. Prognostic factors associated with the success rates of posterior orthodontic miniscrew implants: A subgroup meta-analysis. Korean J Orthod 2016;46(2):111-26.
22. Lim HJ, Choi YJ, Evans CA, Hwang HS. Predictors of initial stability of orthodontic miniscrew implants. Eur J Orthod 2011;33(5):528-32.
23. Lim HJ, Eun CS, Cho JH, Lee KH, Hwang HS. Factors associated with initial stability of miniscrews for orthodontic treatment. Am J Orthod Dentofacial Orthop 2009;136(2):236-42.

Figure 1. Univariable analysis of Generalized Estimating Equations
Figure 1. (continued) Univariable analysis of Generalized Estimating Equations
Figure 2. SPSS outcome after univariable analysis for each reference category using unstructured working correlation matrix
Figure 3. 2x2 contingency tables by various reference categories
Figure 4. SPSS outcome after univariable analysis for each reference category using independent working correlation matrix
Figure 5. Multivariable analysis of Generalized Estimating Equations
Figure 6. SPSS outcome after multivariable analysis of Generalized Estimating Equations using unstructured working correlation matrix