BIOINFORMATICS WITH PARK-KLEIS

위 표에서 중요한 것.

1. Event 발생 여부 (세 번째 열)

2. 언제 event가 발생했는가 (두 번째 열)

Logistic regression에서는 두 번째 열 정보를 제외한다.

따라서, A환자와 C환자 모두 동일하게 MI - YES 정보를 가져가지만, 1년 후 발생, 9년 후 발생이라는 중요한 정보를 사용하지 않는다.

Survival analysis의 핵심 개념 3가지

1. Survival Probability = S(t)

의미: t 시점까지 event 없이 살아 있을 확률

2. Hazard(t)

의미: 지금 이 순간 event가 발생할 위험 (= event 발생 속도)

3. Hazard Ratio (HR)

의미: 두 그룹의 event 발생 속도 비율

예시

Smoking HR = 2

-> 흡연자는 event가 2배 빠르게 발생. 흡연자가 비흡연자에 비해 2배 발생한다는 의미가 아니다.

Cox model 의 핵심 질문: 어떤 변수들이 hazard를 증가시키는가?

예를 들어 age, smoking, protein level의 변수들이 있을 때, cox 모델은 이 변수들이 event 발생 속도를 얼마나 바꾸는지 모델링한다.

Cox model

$$ hazard(t)=baseline hazard(t)×exp(βX) $$

이 수식을 분해해서 살펴보자.

1. baseline hazard(t)

의미: 아무 변수도 없는 기본 위험. 예를 들어 평균환자, 또는 reference group의 위험.

중요한 점은 baseline hazard가 시간에 따라 변한다는 것이다.

예를 들어 심장관련질환은 시간이 지날수록 risk가 증가한다.

2. exp( βX)

의미: 환자의 risk multiplier

X = 변수 (e.g. age, smoking, protein)

β = 각 변수의 효과

예를 들어, β_smoking = 0.69 라면, exp(0.69) = 2.

따라서 hazard ratio = 2 가 된다.

3. 예시로 살펴보는 전체 hazard

hazard(t)=baseline hazard(t)×exp(βX)

위 수식을 예시로 살펴보자.

baseline hazard = 0.01, exp(βX) = 3 이면, hazard = 0.01 × 3 = 0.03

baseline hazard는 평균적인(reference) 환자의 순간 위험이므로 reference의 위험은 0.01,

어떤 환자의 exp(βX) = 3 이면 total hazard = 0.03이므로, 이 환자의 event 발생 속도는 평균환자보다 3배 높다고 할 수 있다.

4. Hazard Ratio

HR = hazard_interest / hazard_reference

위 예시에서 Hazard Ratio는 0.03/0.01 이므로 3이 된다.

R 패키지가 다운로드 되지 않는 이유는 아주 많지만..

집에서는 잘 다운로드 되던 패키지가 회사나 학교에서 다운로드가 되지 않을 때 시도해볼 만한 방법이 있다.

가끔 보안의 문제 때문에 보안이 심한 회사나 학교에서 R패키지 다운로드가 안될 때가 있다.

특히 그 중에서도 Bioconductor 패키지 설치가 안되는 경우,

R 콘솔에 아래 명령어를 입력한 후, 다시 BioManager::install("원하는 패키지") 명령어를 이용해 재시도하면 다운로드가 된다. 

Sys.setenv(R_LIBCURL_SSL_REVOKE_BEST_EFFORT=TRUE)

Sys.setenv(R_LIBCURL_SSL_REVOKE_BEST_EFFORT=TRUE)
BiocManager::install()

만약 이렇게 해도 설치가 되지 않는다면 R버전, 혹은 패키지 버전 충돌 등 다른 방법을 찾아봐야한다.

Repeatability (반복성)

• 정의: 동일한 실험자가 동일한 실험 조건하에서 동일한 장비를 사용하여 짧은 시간 간격으로 여러 번 측정했을 때 나타나는 측정값들 간의 일관성
• 중요성: 이는 특정 조건 하에서 측정 방법이나 측정 장치가 얼마나 일관된 결과를 제공하는지를 보여줌. 즉, 실험 방법이나 장비의 안정성과 정밀성을 평가
• 통계적 도구: 평가를 위해 Coefficient of Variation (CV)이나 Repeatability Coefficient 같은 지표가 사용됨.

Reproducibility (재현성)

• 정의: 다른 실험자들이 다른 조건에서 다른 장비를 사용하여 같은 실험을 반복했을 때 나타나는 결과들 간의 일관성. 조건의 변화는 실험 장소, 시간, 장비, 실험자 등을 포함함.
• 중요성: 재현성은 연구 결과가 다른 연구자나 다른 조건에서도 유사한 결과를 낼 수 있는지를 평가하는 데 중요. 즉, 연구 결과의 신뢰성과 일반화 가능성을 나타냄.
• 통계적 도구: 재현성을 평가하기 위해 Intraclass Correlation Coefficient (ICC)나 ANOVA, Effect Size 등이 사용됨.

Repeatability가 필요한 경우의 예시

: 한 환자를 대상으로 같은 혈압측정 도구를 이용해 여러 번 혈압을 측정할 때, 그 혈압 측정 기기를 통해 측정한 혈압이 얼마나 일관적으로 나오는지?

Reproducibility가 필요한 경우의 예시

: 여러 환자의 MRI microbleed 개수를 여러 radiologist가 측정할 때, 각 환자의 microbleed 개수를 radiologist 간 얼마나 비슷하게 측정했는가?

이번 글에서는 통계적 검정에 있어서 많은 논란이 (?) 있어 왔던 모수적 검정과 비모수적 검정 방법에 대해 이야기 해보고자 한다.

이런건 예시를 들어서 하는게 제일 좋으므로, 바로 예시를 보면서 모수검정을 할지 비모수검정을할지, 그 의사결정 과정은 어떻게 되는지 살펴보도록 하겠다.

시중에 나와 있는 과자들의 평균 무게가 15g 이라고 한다.

이를 검정하기 위해 총 25개 종류의 과자를 모았고, 유의수준 5%에서 검정한다고 치자.

표본 = (15, 19, 14,  3, 10, 18, 22, 11, 5, 20, 14, 22, 25, 26, 27, 5, 19, 27, 25, 28, 25, 9, 29, 3, 8) 

이와 같이 25개 종류의 과자들을 샘플로 선택하였고, 평균을 구해보니 17.16g 이었다.

그럼, 이 때 귀무가설과 대립가설은 다음과 같다.

$$ H_{0} :  \overline{X}=15 $$

$$ H_{1} : not ~ H_{0} $$

이를 검정하기 위해서는 대표적인 parametric test (모수검정)인 one-sample t-test를 시행할 수 있다.

혹은 wilconxon signed-rank test와 같은 비모수 검정법을 이용 할 수 있다.

그럼 우리 예시를 검정하기 위해 어떤 test를 사용하면 좋을까?

1. 정규성 검정 후 검정법 선택

많은 통계 교과서 혹은 통계 강의에서 t-test 전에, 우선 샘플데이터들이 정규분포를 만족하는지 검정해야 한다고 서술한다.

이를 정규성 검정이라고 하는데, 정규성 검정의 대표적인 test는 shapiro-wilk's 검정이 있다.

만약 shapiro-wilk's 검정에서 p-value가 0.05보다 크면 정규성을 만족한다고 보고, 그대로 t-test를 진행하며, p-value가 0.05보다 작게 나오면 정규성을 만족하지 못하여 비모수 검정 선택을 거의 교과서적인(?) 정석으로 나타내는 경우가 많다.

뿐만 아니라 논문을 작성할 때, demographics 표에서 정규성 검정을 통과한 항목들은 평균과 표준편차로 작성, 정규성 검정을 통과하지 못한 항목들은 중앙값과 (최소값, 최대값) 혹은 중앙값과 (사분위수) 등으로 표현하기도 한다.

그러나 shapiro-wilk's 검정은 샘플사이즈의 영향을 크게 받기 때문에, 대략 히스토그램을 그린 후 해당 분포가 대략 정규분포를 따르는 것처럼 보이면 정규성을 만족한다고 보고 t-test를 시행할 수 있다. 

2. 샘플사이즈 고려 후 검정법 선택

한편 또 다른 통계 교과서, 강의들에서는 정규성을 만족하는지 살펴볼 필요 없이 샘플 크기가 대략 30개 이상이면 바로 t-test를 진행해도 무방하다고 말하기도 한다. 중심극한정리에 따라 샘플 수가 커질수록 표본평균은 정규분포를 따르게 되므로 샘플크기가 크다면 굳이 정규성 검정을 할 필요 없이 바로 power가 더 높은 모수 검정 (parametric test) 을 진행해도 좋다는 입장이다.

다만, medical data는 그 특성 상 한쪽으로 치우친 (extrememly skewed data) 가 있는데, 그런 경우는 샘플사이즈가 25를 넘어도, outlier 때문에 non-parametric test를 사용하는 것이 권장된다.

📖자 그럼 우리의 경우 어떤 검정법을 사용하는 것이 좋을까?

① 우리 데이터의 샘플 사이즈는 25이다. 30개 이상은 아니지만, 그렇다고 엄청나게 적은 수의 샘플사이즈는 아니다. 

-> 음.. 좀 애매하다. 데이터 분포를 살펴봐야겠다!

② 그럼 정규성 만족 여부를 살펴볼까? 앞서 말했듯이 shapiro-wilk's 검정보다는 히스토그램을 그려 정규성을 만족하는지 보는 것이 좋다. 우리 샘플 데이터의 히스토그램은 다음과 같다.

자, 위 히스토그램을 보고, 정규분포에 근사한다고 할 수 있을까?

-> 완전한 정규분포 모양은 아니지만.. 또 완전 치우친 분포도 아니고, outlier 값도 없는 것 같은데..

이런 경우가 제일 어렵다.

일단 최근 많은 통계논문들에서 t-test전 정규성 검정을 반드시 해야할 필요는 없다고 서술하고 있다. 

우리 예시와 같은 경우에는 t-test도 무리가 없을 것 같다. 

그 이유는

1. 데이터 분포가 한쪽으로 극단적으로 치우쳐진 것이 아니라는 점

2. 샘플사이즈가 25개 정도 된다는 점이다.

많은 데이터를 보면서 감을 익히는 것과 꾸준히 통계 최신 논문들을 읽는 것이 좋다.

Mixed effect model, 우리 말로는 혼합효과모형 (or 혼합모형) 이라는 것이 있다. 다른 말로는 Generalised Linear Mixed Model (=GLMM) 으로 표현되기도 한다.

처음 접할 때 매우 어렵다고 느껴지는 부분이기도 하고, 실제로도 어렵다. (나만 그럴수도..)

Mixed effect model 은 fixed effect (고정효과) 와 random effect (랜덤효과) 가 모두 고려된 모델인데, 고정효과는 보통 우리가 생각하는 일반적인 GLM의 설명변수들이다.

Fixed effect 는 보통 우리가 관심있어 하는 변수로, 종속변수와 어떤 연관성이 있다고 생각되는 변수들이다.

예를 들어, 알츠하이머성 치매 발병 여부를 종속변수라고 했을 때, 치매 발병과 관계가 있는 것으로 생각되는 변수들은 나이, 성별, 아밀로이드 베타 단백질 침착량 등이 고정효과가 될 것이다.

그럼 Random effect (랜덤효과) 는 무엇일까?

흔히(?) 보는 랜덤효과에 대한 설명으로는 cluster (=군집) 라고 하는데, 나에게는 그렇게 직관적인 설명으로 느껴지지 않는다.

그럼 이번엔 cluster (군집) 의 개념에 대해 살짝 짚어보고 갈 필요가 있겠다.

군집이란 비슷한 or 연관된 속성을 가진 개체를 그룹화한 것을 뜻한다.

예를 들면 마케팅에서 비슷한 소비 성향을 가진 소비자들을 군집으로 묶어 맞춤 마케팅 전략을 세울 수 있는데, 성별과 연령, 주거지역, 직업 등에 따라 비슷한 성향의 소비자들을 군집으로 묶는 것이다.

그래서 랜덤 효과는 군집이다. 라는 설명이 나오는데, 솔직히 이게 더 혼란스럽게 만든다. 

심지어 통계 책에 따라 랜덤효과를 제각기 설명해두기 때문에 더 혼란스럽다. 🤨

따라서 Gelman & Hill (2007) 이 발췌하여 정리한 fixed & random effect 의 개념들을 정리해보고자 한다.

📖 Gelman & Hill (2007)

1. Kreft & De Leeuw (1998)

 - Fixed effect : 모든 individual에 걸쳐 constant 하다. - Random effect : 모든 individual 에 따라 달라진다.   ⇨ 랜덤 intercept ɑᵢ 와 고정 slope β 는 각 individual i에 대해 평행한 선을 가지게 된다. 

2. Searle, Casella, and McCulloch (1992)

 - Fixed : Effects are fixed if they are interesting in themselves. - Random : Effects are random if there is interest in the underlying population.

3. Green & Tukey (1960)

 - Fixed : 샘플이 모집단을 대표할 수 있을 때, corresponding variable 은 고정. - Random : 샘플이 모집단의 small part 일 때, corresponding variable 은 랜덤.

4. LaMotte (1983)

 - Fixed : 아래 랜덤 효과의 경우가 아닐 때 - Random : 어떤 효과가 랜덤 변수로부터의 'realized value'라고 추정될 때, 그게 랜덤 효과. *참고로, 여기에서 realized value 라는 것은, 예를 들어 6면체의 주사위를 던진다고 할 때, 랜덤 변수는 경우의 수인 (1,2,3,4,5,6) 이고, 실제 주사위를 던져서 나온 값이 3이라면, realized value는 3인 것이다.

5. Robinson (1991), Snijders & Bosker (1999)

 - Fixed : least squares (최소 제곱) 또는 maximum likelihood 를 이용해서 추정되었다면 고정 효과. - Random : Shrinkage 로 추정되었다면 랜덤 효과.

자, 그럼 여기까지 위 다섯 개의 고정효과 및 랜덤효과에 대한 정의가 모두 이해 되신 분 있나요?

솔직히 위 다섯 개 정의 모두 좀 ambiguous 한데.. Gelman & Hill (2007) 은 위 정의들에 대해 다음과 같은 견해를 보였다.

1. 첫 번째 정의는 다른 네 가지 정의와는 구별되는 점이 있다.

2. 두 번째 정의의 문제점은, 데이터나 디자인이 바뀌지 않았더라도, 해석의 목적을 달리하면 고정효과가 랜덤효과가 되기도 하고, 랜덤효과가 고정효과가 되기도 한다는 점이다.

3. 세 번째 정의는 유한한 모집단을 정의했다는 점에서 다른 네 가지 정의와 다르다. (그런데, 모집단에서 small 샘플이 아니고, 모집단의 모든 특성을 갖고 있지 않은 large 샘플이라면 어떻게 정의할 것인지? 에 대한 답은 없다.)

4. 네 번째 정의는 actual population에 대한 레퍼런스가 없다는 문제가 있다.

5. 주어진 효과가 고정인지 랜덤인지로 간주해야 하는지 unclear 하다는 문제가 있지만, 수학적 정확성이 있다.

Gelman & Hill (2007) 는 "fixed" 이니 "random" 이란 용어는 피하고, 모델 그 자체를 설명하는 편을 추천하였다.

그런데 저렇게 말하면.. 뭐 어쩌라는건지.. 라는 생각밖에..

그래서! 저자들은 다음과 같은 세 가지의 motivation을 서술해서 어떤 식으로 GLMM을 설명하면 좋을지에 대해 설명해주었다. 지금껏 꽤 ambiguous한 효과의 정의보다 저자들이 말하는 다음 세 가지를 중점적으로 살펴보는 것이 더 이해에 큰 도움이 될 것이다.

1. Accounting for individual- and group-level variation in estimating group-level regression coefficients.

2. Modeling variation among individual-level regression coefficients.

3. Estimating regression coefficients for particular groups.

자세한 설명은 다음 글에서 할 예정입니다.

지난 글에서는 여러번 비교를 할 때 발생할 수 있는 1종 오류 가능성의 증가에 대해 살펴보았다.

2023.02.01 - [STATISTICS] - 다중비교 (Multiple comparisons test) [1]

여러 번 비교를 하게 되면, 1종 오류의 가능성이 증가하게 되는데, 그럼 이를 조정할 필요가 있다.

이걸 예를 들어 다시 간략하게 설명해 보겠다.

A그룹과 B그룹, C그룹이 있고, 세 그룹의 평균을 비교한다고 하자.

그럼,

귀무가설 H0 : A = B = C

대립가설 H1 : H0가 아니다.

대립가설이 성립하려면 ① A ≠ B, ② A ≠ C, ③ B ≠ C 이 셋 중 한 경우라도 만족하면 된다. 

이게 지난 글에서 예시로 들었던 세 개의 문제가 있는 두 과목을 선택하는 것이었다.

두 교과목에서 다루는 세 개의 문제의 난이도가 같다고 할 때,

 - 한 교과목 (아랍어 과목) 은 세 개의 문제를 모두 풀어야 학점을 받을 수 있고,

 - 다른 교과목 (중국어 과목) 은 세 개의 문제 중 한 개라도 풀면 학점을 받을 수 있다고 하면,

어떤 교과목을 선택할 것인가?

웬만한 사람들은 세 개의 문제 중 한 개라도 풀면 학점을 받을 수 있는 교과목인 중국어 과목을 선택할 것이다.

다중비교도 이와 같다.

세 개의 문제를 A, B, C 라고 생각해보면,

① A ≠ B, ② A ≠ C, ③ B ≠ C

이 셋 중 한 경우라도 만족하면 귀무가설이 기각된다는 문제가 생긴다.

만약 두 교과목을 공평하게 만들고자 한다면, 아랍어 과목을 중국어 과목보다 더 쉽게 만들거나, 아니면 중국어 과목을 아랍어 과목보다 훨씬 어렵게 만들어야 한다.

제 1종의 오류 ɑ (실제로 귀무가설이 맞는데, 이를 기각하는 문제)는 중요한 문제이고, 이를 최소화 할 필요가 있다.

이를 해결하기 위한 방법은 여러가지가 있는데, 많이 사용되는 방법 두 가지를 소개하고자 한다.

1. FWER control

먼저 첫 번째 방법은

FWER (Family Wise Error Rate) 을 컨트롤 하는 Bonferroni correction이다.

Bonferroni correction은 매우 직관적이고 간단한 방법이다.

ɑ를 검정하려는 가설의 수로 나누어 주면 끝!

FWER를 ɑ(=0.05)로 유지하기 위해, 만약 가설의 수가 5개라면, ɑ(=0.05)를 5로 나누어주면 된다.

그럼 ① A ≠ B, ② A ≠ C, ③ B ≠ C 이 세 가지 가설 검정에서 p-value가 0.05/5=0.01보다 작아야지만 귀무가설을 기각할 수 있게 된다.

만약 가설의 수가 10개라면, 0.05/10=0.005보다 작아야지만 귀무가설을 기각할 수 있다.

Bonferroni correction의 단점은 매우 conservative 하다는 것인데 이는 power가 떨어진다는 말과 일맥상통한다. 만약 유전자데이터처럼 매우 많은 가설을 검정한다면, 기각을 할 수 있는 p-value가 매우 작아지게 된다.

2. FDR control

두 번째 방법은

FDR (False Discovery Rate) 을 컨트롤 하는 Benjamini-Hochberg procedure 이다.

이 방법은 Bonferroni correction 보다는 덜 conservative한 방법이다.

생각보다 간단하니 다음 과정을 천천히 따라해보면 이해가 될 것이다.

검정을 하려는 m개의 가설이 있다면, m개의 p-value를 가장 작은 것부터, 큰 순서대로 나열한다.

가장 작은 p-value를 p(1)이라고 하고, 그 다음 작은 p-value는 p(2), 이런 순서대로 가장 큰 p(m)까지

p(1), p(2), p(3), ..., p(m) 로 나열한다.

FDR을 ɑ 미만으로 컨트롤 하기 위해서 p(k) < k*ɑ/m을 만족하는 가장 큰 k를 찾는 것이다.

예를 들어, ɑ를 가장 많이 사용되는 유의수준인 0.05라고 하자. (ɑ=0.05)

검정하고자 하는 가설이 다음과 같이 10개 (m=10) 이고,

$$  H_{0A}, H_{0B}, H_{0C}, H_{0D},H_{0E},H_{0F},H_{0G},H_{0H},H_{0I}, H_{0J} $$

각 가설의 p-value는 다음과 같다.

$$  P_{A}=0.034, P_{B}=0.047, P_{C}=0.0001, P_{D}=0.0004, P_{E}=0.028, $$

$$ P_{F}=0.0095, P_{G}=0.003, P_{H}=0.00005, P_{I}=0.0019, P_{J}=0.051 $$

이를 p-value가 작은 것부터 나열하면,

k=1, p = 0.00005

k=2, p = 0.0001

k=3, p = 0.0004

k=4, p = 0.0019

k=5, p = 0.003

k=6, p = 0.0095

k=7, p = 0.028

k=8, p = 0.034

k=9, p = 0.047

k=10, p=0.051 

이와 같이 정리할 수 있다.

p(k) < k*ɑ/m을 만족하는 가장 큰 k 찾는 것이라 했으므로,

p-value가 가장 크게 나온 Hr 가설에 대입해보면, cutoff = k*ɑ/m = k*0.05/10 = k*0.005이다.

k를 가설개수(=m)인 10부터 시작한다.

1. k*ɑ/m = 10*0.005=0.05 이기 때문에 cutoff가 0.05,

    k=10일 때의 귀무가설 하에서 p-value는 p(10) = 0.051이므로, 해당 가설은 기각할 수 없다.

2. k*ɑ/m = 9*0.005=0.045 이기 때문에 cutoff가 0.045,

    k=9일 때의 귀무가설 하에서 p-value는 p(9) = 0.047이므로, 해당 가설은  기각할 수 없다.

3. k*ɑ/m = 8*0.005=0.04 이기 때문에 cutoff가 0.04,

    k=8일 때의 귀무가설 하에서 p-value는 p(8) = 0.034이고, 이는 cutoff인 0.04보다 작으므로, 기각 가능하다.

따라서 가장 큰 k는 8 이고, 그 밑으로는 계산할 필요가 없이 전부 기각 가능하다고 결론 내릴 수 있다.

⇨ p-value가 0.051, 0.047인 가설 J와 가설 B을 제외한 나머지 가설들은 귀무가설을 기각할 수 있다.

위 과정을 R 로 구현하면 다음과 같다.

위 R 결과에서도 보듯이 앞에서부터 8개 가설만이 adjusted p-value가 0.05보다 작음을 알 수 있다.

그럼 다중 검정은 여기에서 끝!

보건통계 또는 의학통계를 접해본 사람들은 ANCOVA의 넓은(?) 쓰임새를 알 것이다.

ANCOVA는 정말 정말 많이 쓰이는 모델링 기법 중에 하나이다.

특히 보건통계를 할 때, age나 sex와 같은 변수들을 '보정(control)' 한다는 목적으로 많이 사용하고는 한다.

그러나 ANCOVA는 단순히 보정을 위해 사용하는 모델링 기법이 아니다.

그런데 지금도 아주 많은 논문에서 ANCOVA를 linear regression 기법의 한 종류로 보면서, age와 sex 보정을 위해 사용한다고 보고한다.

따라서 ANCOVA가 잘못 사용된 예시를 찾는 것은 전-혀 어려운 일이 아니다.

NEJM이나 JAMA 저널에서 ANCOVA와 control을 검색하면 잘못된 논문을 찾는 것은 일도 아니다.

아무튼 따라서 ANCOVA의 오용에 대해서 한 번 살펴 봐야겠다는 생각이 들어,

Miller & Chapman (2001)의 논문인 "Misunderstanding Analysis of Covariance" 를 리뷰해보도록 하겠다.

벌써 이 논문이 publish 된 지 20년이 넘었는데, 아직도 ANCOVA가 저렇게 오용된다니..

위 논문에서 제시한 예시를 들어보도록 하겠다.

-3학년과 4학년 사이에 농구 퍼포먼스에 대한 차이가 있는지 모델링을 해보고 싶다고 하자.

농구 퍼포먼스를 종속변수로,

age를 potential covariate으로,

학년을 grouping variable로 모형을 적합하려고 한다.

농구퍼포먼스 ~ age+ 학년

여기에서 ANCOVA를 적합하는 목적은 학년에 따른 퍼포먼스의 차이가 있음을 증명하고, 그것이 age에 의한 차이가 아님을 말하고 싶기 때문이다.

근데 생각해보면.. 나이와 학년은 선형적인 관계를 갖고 있지 않나?

나이가 많을수록 학년은 높을 것이고, 나이가 어릴수록 학년은 낮을 것이다.

바로 이게 문제이다! 

이 모델에서 Age를 공변량으로 넣으면 나이와 연관된 농구실력의 variance가 제거가 되고, age는 학년과 매우 상관관계가 높기 때문에 학년과 연관된 variance가 상당히 제거될 것이다.

따라서 이 모형의 ANCOVA 결과는 매우 의미가 없게 된다.

Covariate으로 인해 관심 있는 두 그룹 변수(독립변수)의 variance를 크게 제거하게 되면,

독립 변수에 covariate보다 구조적으로 이미 존재하는 difference가 있을 때, 문제가 발생할 수 있다.

이 때 covariate은 이러한 구조적인 differences를 그대로 두기 때문에, treatment effect 추정에 bias가 생길 수 있다. 이런 오류를 specification error라 한다.

Nonrandom Group Assignment

임상시험에서 치료군과 대조군을 정한다고 하자.

두 군에 속할 사람들을 정할 때, nonrandomly하게 정한다면 치료 전에 이미 존재했던 차이가 random error에서 오는지, 혹은 true group difference에 의해서 오는지 알 수가 없다.

이 문제는 treatment effects에 대한 해석을 어렵게 한다.

왜냐하면, 치료의 주 효과(main effect)가 ①treatment effect 및 pre-treatment의 interaction에서 오는지, 아니면 ②치료와 pre-treatment의 miningful overlap(variance shared)에 의해서 오는지 알 수 없기 때문이다.

"Pre-treatment 차이"와 "그룹 factor"의 overlap에 대한 잘못된 해석은 많은 경우 무시된다.

pre-existing group difference는 psychopathology 연구에서 매우 흔하게 발생한다고 한다.

Understanding Analysis of Covariance

ANCOVA는 ANOVA의 한 종류로 보면 된다.

ANCOVA는 독립변수 테스트의 파워를 개선하기 위해 개발된 것이지, 그 어떤 것을 "control"할 목적으로 생긴 것이 아니다.

ANOVA와 ANCOVA를 multiple regression and correlation(=MRC)의 framework 으로, 또한 General Linear Model(=GLM)로 이해하면 더 접근이 쉬울 것이다.

ANCOVA를 MRC의 context에 놓고, ANCOVA의 covariate을 regression predictor로 보는 연구들이 있다.

그런데 보통은 이러한 시각의 연구들 보다는 다른 쪽으로 이해하려는 것이 더 흔하다.

ANCOVA를 두 모델의 오차제곱합을 비교함으로써 그룹의 주효과를 테스트할 수 있다고 본다.

아래 두 모델을 살펴보자.

$$ Y_{ij}=\mu+\alpha_{j}+\beta X_{ij}+\epsilon_{ij} $$

$$ Y_{ij}=\mu+\beta X_{ij}+\epsilon_{ij} $$

$$ Y_{ij} $$

= j번째 그룹의 i번째 subject의 dependent variable

μ = grand mean

Biostat 논문을 작성하면서 한 환자에 대한 두 전문가의 의견에 대한 일치도를 통계적으로 검정해야할 일이 있다. 

즉, KAPPA 통계량은 chance에 의한 Agreement와 observed agreement를 비교한다.

Kappa 통계량을 구하는 식은 다음과 같다.

K = (Proportion of observed agreement - Proportion of agreement by chance) / (1 - Proportion of agreement by chance)

Kappa 통계량의 예시를 구하기 위해서 다음과 같은 예시를 들고 왔다.

(출처: Biostatistics and Epidemiology)

환자는 medical therapy(이하 Med) 혹은 침습 치료(이하 Sur)를 받을 수 있고,

둘 중 어떤 것이 더 각 환자에게 적절한지 두 전문가가 치료법을 선택하고자 한다.

그럼 아래와 같은 표를 얻을 수 있다.

📍총 661 명의 환자가 있었고, 두 전문가가 아래와 같은 판단을 내렸다.

그럼 두 전문가의 의견이 일치한 cell은 a 셀 및 d 셀,

두 전문가의 의견이 일치하지 않은 cell은 b 셀 c 셀 이다.

📝 a 셀과 d 셀의 expected number by chance를 구해보자.

❗위 2X2 테이블의 카테고리가 독립이라 했을 때,

환자가 expert1에서 medical일 확률 P(Expert1∩Med)=P(Expert1)*P(Med) 이므로,

expert1에서 medical일 n수는 (=Expected frequency를 구하는 식은) 다음과 같다.

 $$ N \times P(Med) \times P(Sur) = N \times [ \frac{a+c}{N} \times  \frac{a+b}{N}  ] $$

1. a 셀의 expected number by chance를 a'라고 표기한다면,

a'=(a+c)(a+b)/N

2. d 셀의 expected number by chance를 d'라고 표기한다면,

d'=(c+d)(b+d)/N

따라서 a'와 d'는 각각 다음과 같다.

$$ a' : \frac{494 \times 488}{661}=365 $$

$$ d' : \frac{167 \times 173}{661}=44 $$

이에 따라, the proportion of agreement expected by chance alone[=P(exp)]은

$$  \frac{365+44}{661}=0.619 $$

📝이번에는 Proportion of observed agreement를 구해봅시다.

이는 쉽게 구할 수 있는데, a셀과 d셀이 의견 일치 셀이므로, 

Proportion of observed agreement를 P(obs) 라고 표기한다면

$$ P(obs)= \frac{397+76}{661}=0.716 $$

📖 위에서 구한 수치를 Kappa 통계량 식에 대입해봅시다.

K = (Proportion of observed agreement - Proportion of agreement by chance) / (1 - Proportion of agreement by chance)

여기에서 Proportion of observed agreement는 0.716이고,

Proportion of agreement by chance는 0.619이므로,

$$ Kappa= \frac{0.716-0.619}{1-0.619}= \frac{0.097}{0.381}=0.25  $$

❓결론: Kappa 통계량은 0.25이다.

두 전문가의 의견이 완전 일치할 때에 Kappa는 1이 나오고, 불일치할 때에는 0이 나온다.

따라서, 0.25는 높은 일치도가 아니라고 할 수 있다.