데이터마이너를 꿈꾸며

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을 이용한 누구나 하는 통계분석(안재형 저)