데이터마이너를 꿈꾸며

[예] 당뇨병에 걸린 20명의 환자에 대해 혈당을 낮추는 서로 다른 다섯가지 치료법의 효능을 비교하고자 환자 20명을 랜덤하게 5룹으로 나누어 각각의 치료법을 적용하여 한달 후의 혈당량 수치를 측정하였다.

그러나 한 달 후의 혈당량 수치가 초기 혈당량 수치에 영향을 받을 것으로 생각하여 초기 혈당량 

수치도함께 측정하였다

 - 독립변수: 치료법(trt)

 - 공변량: x(초기 혈당 수치)

 - 반응변수: y(한달 후 혈당 수치)

 - 반응변수y의 모평균에 영향을 끼칠 수 있는 또 다른 변수가 존재 -> 공변량의 영향을 고려해야 함.

1. 문제제기 - 공분산 분석의 필요성

치료법(trt)에 따른 혈당이 낮아지는 효과를 알아보려 하였으나 초기수치(x)가 낮은 경우 한달 후 수치(y)도 낮아질 것으로 예상된다면, 정확한 치료법(trt)의 효과를 알 수가 없다. 이 때 반응변수 y에 영향을 끼칠 수 있는 또 다른 변수 초기수치(x)에 대해 공변량(Covariate)라 하고, 초기수치(x)를 보정한 상태(adjusted)에서 한달 후 혈당량 수치의(y) 보정된 모평균에 차이가 있는지 보는 것을 공분산분석(Analysis of Covariance: ANCOVA)라 한다. 공분산 분석은 회귀분석과 분산분석의 결합으로 각 처리안에서 공변량을 설명변수 x로 하여 회귀분석을 실시하며, 이렇게 공변량을 고려하면 이를 고려하지 않은 분산분석보다 추정의 정도(precision)을 높일 수 있다.

2. 공분산분석을 위한 두 가지 가정

1) 각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일해야 한다. 즉 교호작용이 없어야 한다.

2) 공변량x 효과가 0이 아니다. 효과가 0이라면 분산분석을 하면 된다.

3. 공분산 모형

yij = β0 + αi + βXij + εij

 - Yij : i번째 처리에서 j번째 개체의 반응값

 - Xij : i번째 처리에서 j번째 개체의 공변량 값

 - αi : 처리의 효과

 - β : 모든 처리에 공통으로 작용하는 공변량의 효과

 - ε : 등분산을 갖는 정규분포를 따른다고 가정

4. 각 모수의 검정

 1) H01: β = 0

- 귀무가설은 처리효과를 제어한 상태에서 반응변수Y에 미치는 공변량효과가 없다는 가정을 검정.

- 만약 처리와 공변량 사이에 교호작용이 존재하면 처리간에 회귀계수가 동일하지 않다는 것을 의미하고,   귀무가설이 기각되지 않으면 분산분석 시행

 2) H02 : α1 = α2 = ... = αI

- 공변량 효과를 제어한 상태에서 처리 간 반응변수의 차이가 있는지를 검정

  -> 통계 모형을 통해 진짜 알려고 하는 내용

이상 내용에 대해 R 분석을 한 내용은 아래와 같다.

1. library 호출 및 데이터 입력

> library(HH)

> library(lsmeans)

> glucose = read.csv('혈당량자료.csv',header=T)

> head(glucose)

  관측번호 치료법 초기혈당량 치료후혈당량

1        1      A       27.2         32.6

2        2      A       22.0         36.6

3        3      A       33.0         37.7

4        4      A       26.8         31.0

5        5      B       28.6         33.8

6        6      B       26.8         31.7

> colname<-c('obs','trt','x','y')

> colnames(glucose)<-colname

> head(glucose)

  obs trt    x    y

1   1   A 27.2 32.6

2   2   A 22.0 36.6

3   3   A 33.0 37.7

4   4   A 26.8 31.0

5   5   B 28.6 33.8

6   6   B 26.8 31.7

> attach(glucose)

2. 회귀계수의 동일성 검정 (교호작용 존재 확인)

공분산분석을 하기 전에 먼저 살펴보아야 할 가정 중의 하나가 바로 처리 간 회귀계수 β의 동일성이다.즉, 처리마다 공변량 효과가 동일해야 함을 의미하는데 만약 처리와 공변량 사이에 교호작용이 존재하면 처리 간에 회귀 계수가 동일하지 않다는 것을 의미하고, 이 경우 공분산 분석을 하는 것은 바람직하지않다.

따라서 먼저 둘 사이에 교호작용이 존재하는가의 여부를 살펴본 후 공분산 분석을 해야 한다.

첫번째 귀무가설은 처리효과를 제어한 반응변수Y에 미치는 공변량 효과를 검정하기 위한 것이다.

귀무가설 H01: β = 0 

대립가설 H11: β ≠ 0

두 번째 귀무가설은 공변량 효과를 제어한 상태에서 처리 간 반응변수의 차이가 있는지를 검정.

귀무가설 H02 : α1 = α2 = … = αI

대립가설 H12 : not H02

위 두 가지 귀무가설이 기각되지 않으면 공분산분석을 하지 않고 분산분석을 해야 한다.

> model1 = aov(y~factor(trt)*x,data=glucose)

> summary(model1)

              Df Sum Sq Mean Sq F value   Pr(>F)   

factor(trt)    4 198.41   49.60  15.868 0.000248 ***

x              1  92.53   92.53  29.601 0.000285 ***

factor(trt):x  4  36.48    9.12   2.917 0.077290 . 

Residuals     10  31.26    3.13                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-> aov함수: 공분산분석의 가정 중 하나인 처리 간 회귀계수의 동일성을 확인하기 위해 처리와 공변량 사이의 교호작용의 유무를 검정. aov함수의 우변에 factor(trt) * 공변량 X를 입력. 치료법 trt앞에 factor를 입력하는 것은 치료법trt를 의미하는 ABCDE는 각각 명목형 변수이기 때문.

귀무가설 H01: 교호작용이 존재하지 않는다.

대립가설 H11: 교호작용이 존재한다.

p-value : 0.077290

의사결정: 유의수준 5% 하에서 귀무가설을 기각할 수 없다.

결론: 교호작용이 존재하지 않으므로 공분산분석을 할 수 있다. ( = 각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일하다.)

Ancova(HH Library)함수를 이용하면, 치료법 간 공변량의 효과가 동일하다는 가정을 할 수 있는 xyplot을 그릴 수 있다.

> ancova(y~trt*x,data=glucose)

Analysis of Variance Table

Response: y

          Df  Sum Sq Mean Sq F value    Pr(>F)   

trt        4 198.407  49.602 15.8683 0.0002484 ***

x          1  92.528  92.528 29.6012 0.0002846 ***

trt:x      4  36.476   9.119  2.9173 0.0772904 . 

Residuals 10  31.258   3.126                     

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1

-> HH라이브러리 ancova 함수는 aov함수와 동일한 기능, trellis 그림인 xyplot시각화 기능 제공.함수 출력 결과는 위 aov함수와 동일하며, F값은 2.92, p-value는 0.0773으로 유의수준 5% 하에서 유의하지 않으므로 귀무가설을 기각할 수 없다.따라서 공변량과 처리 사이에 교호작용이 존재하지 않으므로 공분산분석을 할 수 있다는 것을 알 수 있다.

3. 일원 공분산분석 (One-way ANCOVA) - 공변량 효과 제어 시 치료법의 효과 검정

처리 간 공변량의 효과가 동일하다는 가정을 확인한 후 공분산 분석을 출력

> model2 = lm(y~factor(trt)+x,data=glucose)

> summary(model2)

Call:

lm(formula = y ~ factor(trt) + x, data = glucose)

Residuals:

    Min      1Q  Median      3Q     Max

-3.1360 -1.0024 -0.2827  0.7257  6.0806

Coefficients:

             Estimate Std. Error t value Pr(>|t|)   

(Intercept)   13.9437     4.8219   2.892 0.011834 * 

factor(trt)B  -2.7685     1.5554  -1.780 0.096793 . 

factor(trt)C  -1.6660     1.6186  -1.029 0.320776   

factor(trt)D  -1.6284     1.5618  -1.043 0.314787   

factor(trt)E  -4.5903     1.9115  -2.401 0.030788 * 

x              0.7534     0.1723   4.373 0.000637 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.2 on 14 degrees of freedom

Multiple R-squared:  0.8112,    Adjusted R-squared:  0.7437

F-statistic: 12.03 on 5 and 14 DF,  p-value: 0.0001164

-> lm 함수: 공분산분석에는 lm함수를 사용. ~를 중심으로 좌변에는 반응변수y, 우변에는 치료법trt와 공변량 x를 입력하고 그 사이는 *가 아닌 +를 사용.

귀무가설 H01: β = 0 

         H02 : α1 = α2 = … = αI

대립가설 H11: β ≠ 0

 H12 : not H02

p-value : 0.0001164

의사결정: 유의수준 5% 하에서 귀무가설을 매우 강하게 기각할 수 있다.

결론: 우리가 세운 모형의 자료에 적합하다는 것을 알 수 있다. 낮아진 혈당량의 모평균이 처trt의 효과와 공변량들의 효과가 없다라고 말할 수 없다.유의수준 5%하에서 혈당량의 초기수치(x)가 즉 모든 공변량이 0이라는 귀무가설을 기각할만한 증거가 충분하다.

-> 결과 해석:

 - 치료법A를 0으로 둘 때 치료법E p-value가 0.030788으로 유의하게 차이가 난다는 것을 알 수 있다.

1종 제곱합과 회귀계수 동일성 확인 위한 xyplot 그래프는 아래와 같다.

> ancova(y~trt+x,data=glucose)

Analysis of Variance Table

Response: y

          Df  Sum Sq Mean Sq F value    Pr(>F)   

trt        4 198.407  49.602  10.252 0.0004301 ***

x          1  92.528  92.528  19.125 0.0006369 ***

Residuals 14  67.734   4.838                     

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

-> ancova(y~trt+x,data=glucose), 공분산분석이므로 *가 아닌 +를 입력.

-> 그래프 결과해석: 5개 그래프의 적합회귀식 절편은 다르지만 기울기는 동일하므로 공변량 효과가 처리마다 다르지 않음을 알 수 있다. “2. 공분산분석을 위한 두 가지 가정” 중 첫 번째 내용 “각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일해야 한다. 즉 교호작용이 없어야 한다.”을 만족시킨다.

4. 공변량 효과 제어시 치료법의 효과 검정 - 모형 제곱합

1) 제1종 제곱합 Type I SS

SS(trt,X) = SS(trt) + SS(X | trt)

         처리가 기여한 부분 + 처리의 기여 후 순수 공변량 x이 기여한 부분

2) 제3종 제곱합 Type III SS

SS(trt | X) + SS(X | trt)

공변량x가 기여한 상태에서 처리가 기여한 부분 + 처리의 기여 후 순수 공변량 x이 기여한 부분

x given trt, 초기혈당수치 x가 고려된 후 치료법간 차이trt의 차이를 확인한다는 분석의 목적. 즉 x 기여 외 순수 trt 기여를 아는 것이 목표이기 때문에, 우리가 관심가지는 분야 역시 제3종 제곱함 SS(trt | X) 부분

제1종 제곱합

> ancova(y~trt+x,data=glucose)

Analysis of Variance Table

Response: y

          Df  Sum Sq Mean Sq F value    Pr(>F)   

trt        4 198.407  49.602  10.252 0.0004301 ***

x          1  92.528  92.528  19.125 0.0006369 ***

Residuals 14  67.734   4.838                     

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

제3종 제곱합

> summary(aov(y~x+factor(trt),data=glucose))

            Df Sum Sq Mean Sq F value Pr(>F)   

x            1 256.75  256.75  53.067  4e-06 ***

factor(trt)  4  34.19    8.55   1.767  0.192   

Residuals   14  67.73    4.84                  

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1

제1종 제곱합은

- F-값10.252

- p-value: 0.0004301 < 0.05

- 의사결정: ss(trt) 제1종 제곱합은 유의하게 나와 치료법이 5% 유의수준 하에서 유의하다.

제3종 제곱합

- F-값1.767

- p-value: 0.192 > 0.05

- 의사결정: ss(trt | x) x given trt, 공변량x이 기여한 상태에서 처리trt가 기여한 순수한 부분을 확인하는 제3종 제곱합은 치료법trt이 5% 유의수준 하에서 유의하지 않다. 공분산분석은 공변량x의 효과를 제어했을 때 치료법trt에 따라 혈압의 보정평균이 차이가 나는지가 관심사이기 때문에 치료법의 유의성 결과는 제3종 제곱합에 나타난 결과를 보아야 한다.

초기수치

초기수치에 대해서는 제1종 제곱합과 제3종 제곱합이 같게 나오는데 이것은 제1종 제곱합에서는 SS(X | trt)이고 고유기여분을 나타내는 제3종 제곱합에서도 SS(X | trt)이기 때문이다. 초기 수치에 대한 p-value는 0.0006으로 유의수준 5%하에서 매우 유의함을 알 수 있다.

귀무가설 H02 : α1 = α2 = … = αI

p-value: 0.0006

의사결정: 초기수치의 효과가 0이라는 귀무가설을 기각하게 되어 앞의 공분산분석을 하기 위해 만족해야 하는 첫 번째 가정이 충족되는 것을 알 수 있다.

4. 각 요인의 수준별 추정치 분석

> model2 = lm(y~factor(trt)+x,data=glucose)

> summary(model2)

Call:

lm(formula = y ~ factor(trt) + x, data = glucose)

Residuals:

    Min      1Q  Median      3Q     Max

-3.1360 -1.0024 -0.2827  0.7257  6.0806

Coefficients:

             Estimate Std. Error t value Pr(>|t|)   

(Intercept)   13.9437     4.8219   2.892 0.011834 * 

factor(trt)B  -2.7685     1.5554  -1.780 0.096793 . 

factor(trt)C  -1.6660     1.6186  -1.029 0.320776   

factor(trt)D  -1.6284     1.5618  -1.043 0.314787   

factor(trt)E  -4.5903     1.9115  -2.401 0.030788 * 

x              0.7534     0.1723   4.373 0.000637 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.2 on 14 degrees of freedom

Multiple R-squared:  0.8112,    Adjusted R-squared:  0.7437

F-statistic: 12.03 on 5 and 14 DF,  p-value: 0.0001164

위의 각 추정 계수를 바탕으로 공분산 모형식을 쓰면 아래 식과 같다.

ŷ_1j = β^₀ + α^₁ + β^x_1j = 13.9437 + 0 + 0.7534x_1j

ŷ_2j = β^₀ + α^₂ + β^x_2j = 13.9437 -2.7685 + 0.7534x_2j

ŷ_3j = β^₀ + α^₃ + β^x_3j = 13.9437 -1.6660 + 0.7534x_3j

ŷ_4j = β^₀ + α^₄ + β^x_4j = 13.9437 -1.6284 + 0.7534x_4j

ŷ_5j = β^₀ + α^₅ + β^x_5j = 13.9437 -4.5903 + 0.7534x_5j

각 처리의 회귀식을 살펴보면 공변량 초기혈당량의 회귀계수 0.7534는 모두 처리에 대해 공통이고 반응변수y의 절편에서만 값이 차이가 난다는 것을 알 수 있다. 첫 번째 치료법 A의 처리효과를 0으로놓았기 때문에 Inetercept 부분의 회귀계수 13.9437은 첫번째 치료법 A의 상수항임을 알 수 있고, 치료법B부터 치료법 E에 대한 추정치는 각각 첫 번째 

치료법 A의 추정치와의 차이가 난다.

치료법A(trt A)와 치료법 E(trt E)의 차이에 대한 p-value가 0.0308로 유의하게 나와 치료법 A가 치료법E와 유의하게 차이가 난다는 것을 알 수 있다. 

혈당량의 초기 수치(x)에 대한 모수 추정치는 0.753이고 t-값과 p-값이 각각 4.373과 0.000637으로유의수준 5%하에서 혈당량의 초기 수치(x)가 0이라는 귀무가설을 기각할 만한 증거가 충분하다.

즉 공변량 효과가 있다. 이는 공분산결과 분석을 해도 된다.

5. LSMEAN(adjusted mean) 분석

> lsmeans(model2, ~ trt)

 trt   lsmean       SE df lower.CL upper.CL

 A   32.97565 1.151991 14 30.50487 35.44642

 B   30.20716 1.148211 14 27.74449 32.66982

 C   31.30960 1.104799 14 28.94004 33.67916

 D   31.34724 1.117954 14 28.94946 33.74501

 E   28.38536 1.341631 14 25.50785 31.26287

Confidence level used: 0.95

Warning message:

In model.frame.default(trms, grid, na.action = na.pass, xlev = xlev) :

  variable 'trt' is not a factor

보정된(adjust) 평균은 각 처리마다 추정된 회귀식에서 공변량값에 각 처리 평균 대신 공변량의 전체평균을 

사용하였을 때 기대되는 반응변수의 평균이다. 즉, 공변량 효과를 보정한 상태에서의 반응변수의 평균으로, 모든 처리에서 공변량 평균이 같다고 했을 때의 평균이다. 공분산분석에서는 공변량 효과를

보정한 보정 평균 간에 차이가 있는지가 가장 큰 관심사이다. 앞에서 구한 각 처리 회귀식의 공변량에

공변량의 전체 평균값을 넣어서 처리별 보정된 평균을 다음 식과 같이 계산할 수 있다.

y bar_ad1 = β^₀ + α^₁ + β^xbar = 13.9437 + 0 + 0.7534 X 25.26 = 32.975

y bar_ad2 = β^₀ + α^₂ + β^xbar = 13.9437 -2.7685 + 0.7534 X 25.26 =30.207

y bar_ad3 = β^₀ + α^₃ + β^xbar = 13.9437 -1.6660 + 0.7534 X 25.26 =31.309

y bar_ad4 = β^₀ + α^₄ + β^xbar = 13.9437 -1.6284 + 0.7534 X 25.26 =31.347

y bar_ad5 = β^₀ + α^₅ + β^xbar = 13.9437 -4.5903 + 0.7534 X 25.26 =28.385

각 처리마다 계산된 보정된 평균값의 차이는 공변량의 값이 전체 평균으로 같기 때문에 평균값의 차이는 절편의 차이가 됨을 알 수 있다. 각 처리마다 추정된 회귀계수 값들은 첫 번째 처리(trt A)와 추정치 차이인 동시에 첫 번째 처리와의 보정된 평균 차이이기도 하다.

6. Tukey(HSD) 검정 (처리 간 다중 비교)

> model2.lsm = lsmeans(model2,pairwise ~ trt,glhargs=list())

> print(model2.lsm, omit = 2)

$lsmeans

 trt   lsmean       SE df lower.CL upper.CL

 A   32.97565 1.151991 14 30.50487 35.44642

 B   30.20716 1.148211 14 27.74449 32.66982

 C   31.30960 1.104799 14 28.94004 33.67916

 D   31.34724 1.117954 14 28.94946 33.74501

 E   28.38536 1.341631 14 25.50785 31.26287

Confidence level used: 0.95

$contrasts

 contrast    estimate       SE df t.ratio p.value

 A - B     2.76849175 1.555390 14   1.780  0.4216

 A - C     1.66604727 1.618557 14   1.029  0.8378

 A - D     1.62840923 1.561818 14   1.043  0.8316

 A - E     4.59029035 1.911531 14   2.401  0.1718

 B - C    -1.10244448 1.615029 14  -0.683  0.9570

 B - D    -1.14008252 1.560695 14  -0.730  0.9457

 B - E     1.82179860 1.904048 14   0.957  0.8696

 C - D    -0.03763804 1.585115 14  -0.024  1.0000

 C - E     2.92424307 1.690871 14   1.729  0.4485

 D - E     2.96188112 1.832554 14   1.616  0.5115

P value adjustment: tukey method for comparing a family of 5 estimates

> plot(model2.lsm[[2]])

앞서 본 모델에서는 A와 E 치료법 사이에 유의한 차이가 있었지만, LSD 검정(Tukey(HSD)) 검정에서는어느 치료법도 혈당 수치에 있어 유의한 차이가 없다. 왜일까? Tukey검정방법은 상대적으로 보수적이기 때문에 정말 큰 차이가 날 경우에만 유의한 차이를 보인다.

Dunnett vs Tukey

모든 평균이 같다는 귀무가설이 기각되었다는 말은 그룹 중 최소한 하나는 0이 아니다라는 말이다. 어느 쌍의 차이로 귀무가설이 기각되었는지 조사하기 위해 다중비교를 한다. 분산분석에서 많이 쓰이는 다중비교 방법은 Dunnett와 Tukey이다. Tukey는 가능한 모든 조합의 쌍을, Dunnett는 하나의 대조군(reference)을 나머지 비교군(treatment)들과 비교한다.

6-1) Tukey

> library(multcomp)

> model3 = lm(y~trt+x,data=glucose)

> tukey = glht(model3,linfct=mcp(trt='Tukey'))

> summary(tukey)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lm(formula = y ~ trt + x, data = glucose)

Linear Hypotheses:

           Estimate Std. Error t value Pr(>|t|)

B - A == 0 -2.76849    1.55539  -1.780    0.419

C - A == 0 -1.66605    1.61856  -1.029    0.836

D - A == 0 -1.62841    1.56182  -1.043    0.830

E - A == 0 -4.59029    1.91153  -2.401    0.170

C - B == 0  1.10244    1.61503   0.683    0.956

D - B == 0  1.14008    1.56069   0.730    0.945

E - B == 0 -1.82180    1.90405  -0.957    0.868

D - C == 0  0.03764    1.58512   0.024    1.000

E - C == 0 -2.92424    1.69087  -1.729    0.446

E - D == 0 -2.96188    1.83255  -1.616    0.509

(Adjusted p values reported -- single-step method)

결과해석: 모든 치료방법 간 차이의 p-value가 0.05보다 크기 때문에 유의하지 않다.

> plot(tukey)

결과 해석: 치료 방법의 차이 신뢰구간이 0을 포함하고 있으므로 서로 유의하지 않다. 즉, 어느 치료법도 혈당수치에 있어 유의한 차이가 없다.

6-2) Dunnett

> dunnett=glht(model3,linfct=mcp(trt='Dunnett'))

> summary(dunnett)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts

Fit: lm(formula = y ~ trt + x, data = glucose)

Linear Hypotheses:

           Estimate Std. Error t value Pr(>|t|) 

B - A == 0   -2.768      1.555  -1.780   0.2700 

C - A == 0   -1.666      1.619  -1.029   0.7006 

D - A == 0   -1.628      1.562  -1.043   0.6919 

E - A == 0   -4.590      1.912  -2.401   0.0953 .

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Adjusted p values reported -- single-step method)

결과해석: 모든 치료방법 간 차이의 p-value가 0.05보다 크기 때문에 유의하지 않다.

> plot(dunnett)

결과 해석: 치료 방법의 차이 신뢰구간이 0을 포함하고 있으므로 서로 유의하지 않다. 즉, 어느 치료법도 혈당수치에 있어 유의한 차이가 없다.

7. 잔차 분석 

> par(mfrow=c(2,2))

> plot(모형2)

특이값 2번이 존재하지만, 정규분포를 따르고 등분산 가정은 큰 문제는 없다.

출처: 보건정보데이터 분석(이태림, 이재원, 김주한, 장대흥 공저), R을 이용한 누구나 하는 통계분석(안재형 저)

3.2 여러 집단의 비교
3.2.1 1개의 요인을 고려하는 경우
p.86, 자폐아, 정상아, 지진아에 대한 혈청 항원 농도에 대해 조사를 하였다. 
이들 사이에 면역 이상에 대한 차이가 있다고 할 수 있는가?
귀무가설 H0: u1 = u2 = u3, 대립가설 H1: not H0

1) 데이터 입력

> a<-c(755,365,820,900,170,300,325,385,380,215,400,343,415,345,410,460,225,440,400,360,435,450,360)

> b<-c(165,390,290,435,235,345,320,330,205,375,345,305,220,270,355,360,335,305,325,245,285,370,345,345,230,370,285,315,195,270,305,375,220)

> c<-c(380,510,315,565,715,380,390,245,155,335,295,200,105,105,245)

> boxplot(a,b,c,col='yellow',names=c('자폐아','정상아','지진아'))

> library(vioplot)

> vioplot(a,b,c,col='yellow',names=c('자폐아','정상아','지진아'))

2) 각 그룹의 평균과 분산

> sapply(list(a,b,c),mean)

[1] 419.9130 305.0000 329.3333

> sapply(list(a,b,c),var)

[1] 31693.356  4071.875 29224.524

3) 등분산 검정

> sera = c(a,b,c)

  > group = factor(rep(1:3,c(length(a),length(b),length(c))))

  > fligner.test(sera~group)   

  Fligner-Killeen test of homogeneity of variances

  data:  sera by group

  Fligner-Killeen:med chi-squared = 6.8506, df = 2, p-value = 0.03254

결과 해석: p-값이 0.03254이어서 등분산성에 조금은 문제가 있음을 알 수 있다.

4) one way Anova

> out = aov(sera~group)

> out

  Call:

    aov(formula = sera ~ group)

  Terms:

                       group Residuals

  Sum of Squares   185159.3 1236697.2

  Deg. of Freedom         2        68

  Residual standard error: 134.8582

  Estimated effects may be unbalanced

 > summary(out)

              Df   Sum Sq Mean Sq F value  Pr(>F)  

  group        2  185159   92580   5.091 0.00871 **

  Residuals   68 1236697   18187                  

  ---

    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

결과해석:

질문:이들 사이에 면역 이상에 대한 차이가 있다고 할 수 있는가?
귀무가설 H0: u1 = u2 = u3,

대립가설 H1: not H0

p – value: 0.00871

결정: 귀무가설을 기각, 세 그룹 사이 면역 이상에 대한 차이가 있다.

5) 모형 적합성 검토 = 오차검토

> par(mfrow=c(2,2))

> plot(out)

결과 해석: 정규성 분포에 약간 문제가 있지만 큰 문제는 아니다.

6) average profile plot 평균 반응 프로파일 그림 – 효과 크기를 알 수 있는 plot

> plot.design(sera~group)

결과 해석: 그룹1과 그룹2가 유의하게 서로 달랐다.

7) 다중 비교: 어느 그룹 간 차이가 있는지 보자.

> tukey = TukeyHSD(out)

> tukey

  Tukey multiple comparisons of means

  95% family-wise confidence level

  Fit: aov(formula = sera ~ group)

  $group

             diff        lwr       upr     p adj

  2-1 -114.91304 -202.68435 -27.14174 0.0070326

  3-1  -90.57971 -197.82092  16.66150 0.1142305

  3-2   24.33333  -76.28971 124.95638 0.8315432

 > plot(tukey)

결과 해석: 그룹1과 그룹2가 유의하게 서로 달랐다. 그룹1과 2의 차이만이 신뢰도구간을 0을 포함하지 안으므로 유의미하게 다르다고 결론을 내릴 수 있다.

8) LSD 최소유의차검정 test

> pairwise.t.test(sera,group)

Pairwise comparisons using t tests with pooled SD

data:  sera and group

1      2    

2 0.0076 -    

3 0.0938 0.5642

P value adjustment method: holm

결과 해석: 그룹1과 그룹2가 유의하게 서로 달랐다.

3.2.2 2개의 요인을 고려하는 경우
(1) 반복이 없을 때

예제) 장비 사용에 대한 3가지 방법을 연령별로 다르게 교육. 숙지 시간이 연령, 방법에 따라 다른가?

귀무가설 h0: u1 = u2 = u3, 대립가설 h1: not h0

1. 데이터 읽기

> setwd('c:/Rwork')

> data=read.table('device.txt',header=T)

> head(data)

ages way hour

1 under20   A    7

2   20~29   A    8

3   30~39   A    9

4   40~49   A   10

5 above50   A   11

6 under20   B    9

> tail(data)

ages way hour

10 above50   B   12

11 under20   C   10

12   20~29   C   10

13   30~39   C   12

14   40~49   C   12

15 above50   C   14

2. two way ANOVA

> out = aov(hour~ages+way,data=data)

> out

Call:

  aov(formula = hour ~ ages + way, data = data)

Terms:

                      ages        way Residuals

Sum of Squares  24.933333 18.533333  3.466667

Deg. of Freedom         4         2         8

Residual standard error: 0.6582806

Estimated effects may be unbalanced

> summary(out)

Df Sum Sq Mean Sq F value   Pr(>F)   

ages         4 24.933   6.233   14.38 0.001002 **

way          2 18.533   9.267   21.39 0.000617 ***

Residuals    8  3.467   0.433                     

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

결과해석:

귀무가설 h0: u1 = u2 = u3,

대립가설 h1: not h0

결론: p value가 0.05보다 적으므로 H0를 기각, h1를 받아들인다. 숙지 시간은 연령, 방법에 따라 서로 유의하게 다르다.

3. 모형적합성 검토 = 오차검토

> par(mfrow = c(2,2))

> plot(out)

결과 해석: 오차의 등분산성 및 정규성에 문제가 없음을 알 수 있다.

4. 다중비교, 왜 서로 유의한 차이가 났을까?

4.1) 나이 별 보기

> attach(data)

> pairwise.t.test(hour,ages)

Pairwise comparisons using t tests with pooled SD

data:  hour and ages

20~29 30~39 40~49 above50

30~39     1.00     -     -     -     

40~49     1.00  1.00     -     -     

above50   0.18  0.66   0.91   -     

under20   1.00  1.00   1.00  0.13  

P value adjustment method: holm

결과해석: 50대 이상이 다른 나이 수준보다 높았다. 결국 나이 수준이 숙지시간에 차이를 보인 것은 50대 이상이 다른 나이의 수준과 차이가 났기 때문이다.

4.2) 교육 방법 별 보기

> pairwise.t.test(hour,way)

Pairwise comparisons using t tests with pooled SD

data:  hour and way

A     B   

B 0.549 -   

C 0.061 0.125

P value adjustment method: holm

결과해석:

귀무가설 h0: u1 = u2 = u3,

대립가설 h1: not h0

결론: 방법은 A와 C간 차이가 났다.

위에서 본 두 함수 pairwise.t.test(hour,ages) ,pairwise.t.test(hour,way) 에 대하여 그래프를 그리면 아래와 같다.

>par(mfrow=c(1,2))

> plot.design(hour~ages)

> plot.design(hour~way)

5. 다중 비교

> tukey = TukeyHSD(out)

> tukey

Tukey multiple comparisons of means

95% family-wise confidence level

Fit: aov(formula = hour ~ ages + way, data = data)

$ages                  diff        lwr        upr     p adj

30~39-20~29      1.0000000 -0.8568723  2.8568723 0.4057524

40~49-20~29      1.3333333 -0.5235390  3.1902056 0.1877558

above50-20~29    3.3333333  1.4764610  5.1902056 0.0017351

under20-20~29   -0.3333333 -2.1902056  1.5235390 0.9676094

40~49-30~39      0.3333333 -1.5235390  2.1902056 0.9676094

above50-30~39    2.3333333  0.4764610  4.1902056 0.0154324

under20-30~39   -1.3333333 -3.1902056  0.5235390 0.1877558

above50-40~49    2.0000000  0.1431277  3.8568723 0.0348816

under20-40~49   -1.6666667 -3.5235390  0.1902056 0.0810838

under20-above50 -3.6666667 -5.5235390 -1.8097944 0.0009146

$way diff       lwr      upr      p adj

B-A  0.6 -0.5896489 1.789649 0.3666717

C-A  2.6  1.4103511 3.789649 0.0006358

C-B  2.0  0.8103511 3.189649 0.0034083

> plot(tukey) 

결과해석:

귀무가설 h0: u1 = u2 = u3,

대립가설 h1: not h0

결론: 그래프에서도 나이가 20대 미만과 50대 이상에서, 방법은 3번과 1번 그리고 3번과 2번에서 신뢰구간을 0을 포함하지 않으므로 유의한 차이가 있었음을 알 수 있다.

(2) 반복이 있을 때
예) 세 종류의 호르몬 처리와 성별에 따라 혈액 칼슘값에 차이가 있는지 알아보기 위해 남녀 각 15명씩을 선정하여 이들을 세 그룹으로 나누어 세 가지 호르몬 처리를 한 후 혈액 칼슘을 측정하였다.
성별에 따라 혈액 칼슘에 차이가 있는가? 처리와 성별에 대한 교호작용이 존재하는가?

H0: 성별간 차이가 없다. H1: 성별간 차이가 있다
H1: 처리간 차이가 없다, H1: 처리간 차이가 있다.

1. 데이터 입력

> data=read.csv('calcium.csv')

> head(data)

sex way   cal

1   M   A 16.87

2   M   A 16.18

3   M   A 17.12

4   M   A 16.83

5   M   A 17.19

6   F   A 15.86

> tail(data)

sex way   cal

25   M   C 24.46

26   F   C 30.54

27   F   C 32.41

28   F   C 28.97

29   F   C 28.46

30   F   C 29.65

2. two way anova

> out = aov(cal~sex*way,data=data)

> out

Call:

  aov(formula = cal ~ sex * way, data = data)

Terms:

  sex       way   sex:way Residuals

Sum of Squares     4.0627 1146.6420    3.8454   76.2924

Deg. of Freedom         1         2         2        24

Residual standard error: 1.782933

Estimated effects may be unbalanced

> summary(out

Df Sum Sq Mean Sq F value   Pr(>F)   

sex          1    4.1     4.1   1.278    0.269   

way          2 1146.6   573.3 180.355 3.47e-15 ***

sex:way      2    3.8     1.9   0.605    0.554   

Residuals   24   76.3     3.2                    

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

결과해석: 처리에 대한 p-value가 0.0001보다 적게 나와 처리수준 간 모평균 차이가 없다라는 귀무가설을 기각. 성별과 처리도 p value가 0.05를 넘어 교호작용은 없다.

3. 모형적합성 검토 = 오차검토

> par(mfrow=c(2,2))

> plot(out)

결과 해석: 모형적합성 검토, 잔차도를 그려본 결과 오차의 등분산성에 약간의 문제는 있으나 큰 문제는 없음

4. 교호작용 검토 

> par(mfrow=c(1,1))

> with(data,interaction.plot(sex,way,cal))

결과 해석: with(data,interaction.plot(sex,way,cal)) #두 개의 선이 비슷한 거리를 유지하면서 평행에 가까우므로 interaction  교호작용이 없음을 알 수 있다. interaction.plot은 두 그룹변수의 조합으로 y의 평균을 그래프에 넣어 두 그룹 변수가 서로 y의 평균에 영향을 주는지 보는 방법

5. 다중비교

> attach(data)

The following object is masked _by_ .GlobalEnv:

  sex

The following objects are masked from data (pos = 3):

  cal, sex, way

The following objects are masked from data (pos = 4):

  cal, sex, way

> pairwise.t.test(cal,sex)

Error in tapply(x, g, mean, na.rm = TRUE) :

  arguments must have same length

왜 오류나는지 모르겠다. 더 공부 필요.

 > pairwise.t.test(cal,way)

Pairwise comparisons using t tests with pooled SD

data:  cal and way

A       B     

B 0.052   -     

C 8.4e-16 1.2e-14

P value adjustment method: holm

결과 해석:C와 A, 그리고 C와 B간 방법이 유의하게 차이가 났다.

> tukey = TukeyHSD(out)

> tukey

Tukey multiple comparisons of means

95% family-wise confidence level

Fit: aov(formula = cal ~ sex * way, data = data)

$sex

diff        lwr     upr    p adj

M-F 0.736 -0.6076702 2.07967 0.269434

$way

diff        lwr       upr     p adj

B-A  1.609 -0.3822165  3.600217 0.1295236

C-A 13.845 11.8537835 15.836217 0.0000000

C-B 12.236 10.2447835 14.227217 0.0000000

$`sex:way`

diff        lwr       upr     p adj

M:A-F:A  1.548 -1.9385413  5.034541 0.7421633

F:B-F:A  1.956 -1.5305413  5.442541 0.5236718

M:B-F:A  2.810 -0.6765413  6.296541 0.1661169

F:C-F:A 14.716 11.2294587 18.202541 0.0000000

M:C-F:A 14.522 11.0354587 18.008541 0.0000000

F:B-M:A  0.408 -3.0785413  3.894541 0.9990770

M:B-M:A  1.262 -2.2245413  4.748541 0.8686490

F:C-M:A 13.168  9.6814587 16.654541 0.0000000

M:C-M:A 12.974  9.4874587 16.460541 0.0000000

M:B-F:B  0.854 -2.6325413  4.340541 0.9720701

F:C-F:B 12.760  9.2734587 16.246541 0.0000000

M:C-F:B 12.566  9.0794587 16.052541 0.0000000

F:C-M:B 11.906  8.4194587 15.392541 0.0000000

M:C-M:B 11.712  8.2254587 15.198541 0.0000000

M:C-F:C -0.194 -3.6805413  3.292541 0.9999760

결과 해석: 

귀무가설: H0: 성별간 차이가 없다. H1: 성별간 차이가 있다
대립가설: H1: 처리간 차이가 없다, H1: 처리간 차이가 있다.

결론: 성별간에는 유의한 차이는 없지만 방법에는 유의한 차이가 났다. C와 A, 그리고 C와 B간 방법이 유의하게 차이가 났다.

> par(mfrow=c(2,2))

> plot(tukey)

결과 해석: 

귀무가설: H0: 성별간 차이가 없다. H1: 성별간 차이가 있다
대립가설: H1: 처리간 차이가 없다, H1: 처리간 차이가 있다.

결론: 처리 C와 A, C와 B간 유의하게 서로 달랐다.

<또 다른 방법> 

위의 R 코드를 다른 방법으로 해보면 아래와 같다.(출처: R을 이용한 통계 분석, 안재형 지음)

> boxplot(cal~way+sex,col='red',data=data)

교호작용이 있는지 본 후

> with(data,interaction.plot(sex,way,cal))

결과 해석: 두개의 선이 서로 만나지 않으므로 교호작용이 존재하지 않는다는 것을 알 수 있다. (교호작용: 두 그룹 변수가 서로 y의 평균에 영향을 주는지 보는 방법)

분산분석표를 구한다. 교호작용이 존재하면 곱하기(sex*way), 존재하지 않으면 더하기(sex+way)

> out2=lm(cal~sex+way,data=data)

> anova(out2)

Analysis of Variance Table

Response: cal

Df  Sum Sq Mean Sq  F value    Pr(>F)   

sex        1    4.06    4.06   1.3181    0.2614   

way        2 1146.64  573.32 186.0089 3.944e-16 ***

Residuals 26   80.14    3.08                      

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> summary(out2)

Call:

  lm(formula = cal ~ sex + way, data = data)

Residuals:

  Min      1Q  Median      3Q     Max

-5.8170 -0.5815 -0.0335  0.6623  4.3730

Coefficients:

 Estimate Std. Error t value Pr(>|t|)   

(Intercept)  15.6960     0.6411  24.484  < 2e-16 ***

sexM          0.7360     0.6411   1.148   0.2614   

wayB          1.6090     0.7851   2.049   0.0506 . 

wayC         13.8450     0.7851  17.634 5.53e-16 ***

  ---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.756 on 26 degrees of freedom

Multiple R-squared:  0.9349,   Adjusted R-squared:  0.9274

F-statistic: 124.4 on 3 and 26 DF,  p-value: 1.532e-15

결과 해석:

-   sexM(M-F, p-value는 0.2614)의 추정치가 0.7360으로 유의하지는 않다.

-   wayC(C-A, p-value는 5.533-16)의 추청치가 13.8450으로 유의하고 평균은 A보다 높다.

다중비교

> library(multcomp)

> tukey2=glht(out2,linfct=mcp(way='Tukey'))

> tukey2

General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Linear Hypotheses:

  Estimate

B - A == 0    1.609

C - A == 0   13.845

C - B == 0   12.236

> summary(tukey2)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lm(formula = cal ~ sex + way, data = data)

Linear Hypotheses:

  Estimate Std. Error t value Pr(>|t|)   

B - A == 0   1.6090     0.7851   2.049    0.121   

C - A == 0  13.8450     0.7851  17.634   <0.001 ***

C - B == 0  12.2360     0.7851  15.584   <0.001 ***

  ---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Adjusted p values reported -- single-step method)

> plot(tukey2)

결과 해석: 

귀무가설: H0: 성별간 차이가 없다. H1: 성별간 차이가 있다
대립가설: H1: 처리간 차이가 없다, H1: 처리간 차이가 있다.

결론: 방법 C와 A, C와 B는 신뢰구간을 0을 포함하지 않으므로 유의한 차이가 있다는 결론을 내린다 (p-value < 0.001)

출처: 보건정보데이터 분석(이태림 저), R을 이용한 누구나 하는 통계분석(안재형 저)