데이터마이너를 꿈꾸며

4.2.2 독립성 검정

예4.3) Goodman과 Kruskal은 6,800명을 대상으로 눈색과 머리색을 조사하여 얻은 자료를 가지고 아래와 같은 관찰표를 얻엇다. 눈색과 머리색에 따라 3X4 분할표를 구성할 때 눈색이 머리색에 영향을 주는가? 즉, 서로 독립저긴가?

   H0 눈색과 머리색은 독립이다

   H1 눈색과 머리색인 서로관련이 있다

1) 데이터 입력

> out = matrix(c(1768, 807, 189, 47, 946, 1387, 746, 53, 115, 438, 288, 16), nrow=3, byrow = T)

> dimnames(out) = list (eye=c('e1','e2','e3'),hair=c('h1','h2','h3','h4'))

> out

hair

eye    h1   h2  h3 h4

e1   1768  807 189 47

e2    946 1387 746 53

e3    115  438 288 16

> addmargins(out) #분할표 만들기

hair

eye      h1   h2   h3  h4  Sum

e1     1768  807  189  47 2811

e2      946 1387  746  53 3132

e3      115  438  288  16  857

Sum   2829 2632 1223 116 6800

2) 데이터 시각화

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

> dotchart(out)

> dotchart(t(out))

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

> mosaicplot(out)

결과 해석: 눈색과 머릿색이 서로 영향을 주고 있음을 알 수 있다

3) 카이제곱 검정 

> chisq.test(out)

Pearson's Chi-squared test

data:  out

X-squared = 1073.5, df = 6, p-value < 2.2e-16

결과 해석:

  귀무가설 H0 눈색과 머리색은 독립이다.

  대립가설 H1 눈색과 머리색인 서로관련이 있다.

  p-value 2.2e-16 < 0.001

  의사결정: p-value값이 0.001보다 작으므로 눈색과 머리색이 유의하게 서로 영향을 주고 있음을 알 수 있다.

4) 카이제곱 검정 결과 보기

> names(chisq.test(out))

[1] "statistic" "parameter" "p.value"   "method"    "data.name" "observed"  "expected"  "residuals" "stdres"  

> chisq.test(out)$observed #관찰도수

hair

eye    h1   h2  h3 h4

e1 1768  807 189 47

e2  946 1387 746 53

e3  115  438 288 16

> chisq.test(out)$expected #기대도수

hair

eye         h1        h2       h3       h4

e1 1169.4587 1088.0224 505.5666 47.95235

e2 1303.0041 1212.2682 563.2994 53.42824

e3  356.5372  331.7094 154.1340 14.61941

> chisq.test(out)$residuals #잔차

hair

eye          h1        h2         h3          h4

e1  17.502565 -8.519654 -14.079133 -0.13752858

e2  -9.890092  5.018483   7.697865 -0.05858643

e3 -12.791799  5.836008  10.782543  0.36107650

출처: 보건정보데이터 분석(이태림 저)

4.2 범주형 자료의 검정

#예 4.2) 비타민 c가 감기치료에 효과가 있는지 점검. 대조군(control) 그룹 140명에게는 플라시보를 주고 처리군(treat) 그룹 139명에게는 매일 비타민 c를 투여하였다. 아래 분할표를 가지고 비타민 C가 감기에 효과가 있는지 점검

H0 복용군과 비복용군의 감기 이환율 같다
H1 복용군과 비복용군의 감기 이환율 다르다

1) 자료 입력

> vitamin = matrix(c(31,109,17,122),nrow=2,byrow=T)

> dimnames(vitamin) = list(vitamin=c('ctr','trt'),flu=c('y','n'))

> vitamin

          flu

vitamin  y   n

    ctr 31 109

    trt 17 122

> round(vitamin/sum(vitamin),2)

            flu

vitamin    y    n

    ctr 0.11 0.39

    trt 0.06 0.44

> addmargins(vitamin)

          flu

vitamin  y   n Sum

    ctr 31 109 140

    trt 17 122 139

    Sum 48 231 279

2) 데이터 시각화

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

> dotchart(vitamin)

> dotchart(t(vitamin))

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

> mosaicplot(vitamin)

결과 해석: 비타민 복용군과 비복용군의 감기 이환율이 동일하지 않음을 알 수 있다.

3)카이제곱 검정 실행

#카이검정

> chisq.test(vitamin)

Pearson's Chi-squared test with Yates' continuity correction

data:  vitamin

X-squared = 4.1407, df = 1, p-value = 0.04186

결과해석
대립가설 H0 : 복용군과 비복용군의 감기 이환율 같다
귀무가설 H1 : 복용군과 비복용군의 감기 이환율 다르다

p-value : 0.04186

결정: p-value값이 0.05보다 작으므로 H0를 기각, 비타민 복용군과 비복용군 간 이환율은 다르다.

관찰도수, 기대도수, 잔차를 보는 법 

#관찰도수

> names(chisq.test(vitamin))

[1] "statistic" "parameter" "p.value"   "method"    "data.name" "observed"  "expected"

[8] "residuals" "stdres"  

> chisq.test(vitamin)$observed

flu

vitamin  y   n

ctr 31 109

trt 17 122

#기대도수

> chisq.test(vitamin)$expected

flu

vitamin        y       n

ctr 24.08602 115.914

trt 23.91398 115.086

#피어슨잔차

> chisq.test(vitamin)$residual

flu

vitamin         y          n

ctr  1.408787 -0.6421849

trt -1.413846  0.6444908

출처: 보건정보데이터 분석(이태림 저)

제4장 범주형 자료의 분석

4.1 범주형 자료와 분할표

분할표 table 만들기 연습

> medi = read.table('c:/Rwork/medication.txt',header=T)

> head(medi,3)

癤퓆o medication surv

1     1      treat    y

2     2      treat    n

3     3      treat    y

열 이름이 깨져 보임.

> colnames(medi)<-c('no','medication','surv')

> head(medi,3)

no medication surv

1  1      treat    y

2  2      treat    n

3  3      treat    y

분할표 작성, treat = 처리군 , control = 대조군

> attach(medi)

> tab <- table(medication, surv)

> colnames(tab) = c('die','survival')

> rownames(tab) = c('trt','ctr')

> tab

surv

medication die survival

trt   1        4

ctr   3        2

> addmargins(tab)

surv

medication die survival Sum

trt   1        4   5

ctr   3        2   5

Sum   4        6  10

> tab/sum(tab)

surv

medication die survival

trt 0.1      0.4

ctr 0.3      0.2

> addmargins(tab/sum(tab))

surv

medication die survival Sum

trt 0.1      0.4 0.5

ctr 0.3      0.2 0.5

Sum 0.4      0.6 1.0

출처: 보건정보데이터 분석(이태림 저)

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

제3장 연속형 자료의 분석

3.1 두 집단의 평균비교

3.1.1 독립표본의 평균비교 two sample test

예제) 흡연자 집단과 비흡연자 집단 간 폐 파괴지수를 측정하였다. 높은 수치는 폐의 손상이 크다는 것을 뜻한다. 흡연자와 비흡연자의 폐 파괴지수의 평균이 같다고 할 수 있는가? (각 그룹에서의 관측치들은 정규분포를 따르는 모집단으로부터 독립적으로 얻어진 것이며 두 그룹에서의 모분산은 같다고 가정하자. )

귀무가설 H0: 흡연자와 비흡연자의 폐 파괴지수 평균은 같다.

대립가설 H1: 흡연자와 비흡연자의 폐 파괴지수 평균은 다르다. (양측 검정)

#1. 자료 입력

> smoke=c(16.6,13.9,11.3,26.5,17.4,15.3,15.8,12.3,18.6,12,24.1,16.5,21.8,16.3,23.4,18.8)

> nonsmoke=c(18.1,6,10.8,11,7.7,17.9,8.5,13,18.9)

> sapply(list(smoke,nonsmoke),mean)

[1] 17.53750 12.43333

> sapply(list(smoke,nonsmoke),sd)

[1] 4.475247 4.849227

#2. 정규성 검정

> qqnorm(smoke,main='smoke')

> qqline(smoke,col='blue')

> shapiro.test(smoke)

Shapiro-Wilk normality test

data:  smoke

W = 0.94511, p-value = 0.4163

결과 해석: shapiro.test의 결과에 따라 p value = 0.4163 > 0.05 이므로 귀무가설 기각 못한다, 즉 정규분포를 따른다

> qqnorm(nonsmoke,main = 'nonsmoke')

> qqline(nonsmoke,col='red')

> shapiro.test(nonsmoke)

Shapiro-Wilk normality test

data:  nonsmoke

W = 0.90366, p-value = 0.274

#boxplot과 vioplot

> boxplot(smoke,nonsmoke,col='yellow',names=c('smoke','nonsmoke'))

> library(vioplot)

> vioplot(smoke,nonsmoke,col='yellow',names=c('smoke','nonsmoke'))

결과 해석:  두 집단에 차이가 있음을 알 수 있다.

#3. 두 모분산 비교 (양측검정)

#대립 가설의 형태: alternative = c('two.sided','less','greater')

> var.test(smoke,nonsmoke)

F test to compare two variances

data:  smoke and nonsmoke

F = 0.8517, num df = 15, denom df = 8, p-value = 0.7498

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

  0.2076714 2.7243799

sample estimates:

  ratio of variances

0.8517046

결과 해석: p value 가 0.7498로 분산이 같다는 귀무가설을 기각할 수 없다. 즉, 등분산 가정

#4.두 모분산 비교 (양측검정) - 등분산 가정

> t.test(smoke,nonsmoke,var.equal = T)

Two Sample t-test

data:  smoke and nonsmoke

t = 2.658, df = 23, p-value = 0.01405

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

  1.131680 9.076653

sample estimates:

  mean of x mean of y

17.53750  12.43333

결과해석:

귀무가설 H0: 흡연자와 비흡연자의 폐 파괴지수 평균은 같다.

대립가설 H1: 흡연자와 비흡연자의 폐 파괴지수 평균은 다르다. (양측 검정)

결정: p value는 p value는 0.01405로 두 모 평균이 같다는 귀무가설을 기각한다. 즉, 두 모평균이 서로 다르다.

#5.두 모분산 비교 (양측검정) - 이분산 가정

> t.test(smoke,nonsmoke)

Welch Two Sample t-test

data:  smoke and nonsmoke

t = 2.5964, df = 15.593, p-value = 0.01978

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

  0.9279143 9.2804190

sample estimates:

  mean of x mean of y

17.53750  12.43333

결과해석: 등분산 가정과 큰 차이는 없다.

3.1.2 짝지은 표본의 평균비교 paired sample test

예) 환자 15명에게 혈압강하제를 12주 투입 후 혈압을 비교하였다. 새로운 약은 효과적인가?

귀무가설 h0: u1-u2 = 0

대립가설 h1: u1 > u2

1) 데이터 입력

> before=c(90,56,49,64,65,88,62,91,74,93,55,71,54,64,54)

> after=c(72,55,56,57,62,79,55,72,73,74,58,59,58,71,61)

> diff = before - after

2) 정규성 차이: shapiro – wilk test

> qqnorm(diff)

> qqline(diff,col='red')

> shapiro.test(diff)

Shapiro-Wilk normality test

data:  diff

W = 0.90982, p-value = 0.1345

결과 해석: shapiro test 결과 p value는 0.1345로써 정규분포를 이루고 있다고 할 수 있다.

3) Paired sample test

> mean(diff) ; sd(diff)

[1] 4.533333

[1] 9.425396

> t.test(before, after, paired = T, alternative = 'greater') #μ 복용 전 > μ 복용 후

Paired t-test

data:  before and after

t = 1.8628, df = 14, p-value = 0.0418

alternative hypothesis: true difference in means is greater than 0

95 percent confidence interval:

  0.2469617       Inf

sample estimates:

  mean of the differences

4.533333

귀무가설 H0: 약 복용전과 복용 후의 혈압 수치는 같다. μ 복용 전 = μ 복용 후

대립가설 H1: 약 복용전 대비 복용 후의 혈압 수치기 다 낮다. μ 복용 전 > μ 복용 후

결론: 단측 검정에 대한 p value가 0.0418로서 유의수준 5%에서 두 그룹의 혈압 차이가 없다는 귀무가설을 기각할 만한 충분한 증거가 있으므로 새로운 약이 혈압을 내린다고 볼 수 있다.

출처: 보건 정보 데이터 분석(이태림 저자)