데이터마이너를 꿈꾸며

[예] 어떤 약물에 대한 체내 배출연구에서 얻은 자료이다. 연구자는 약의 형태에 따라 체내로부터 배출되는 약물의 양이 달라지는지를 알고자 한다. 그런데 배출되는 약물의 양은 약의 형태뿐만 아니라 배출된 약물을 측정한 시간과 각 개체의 항신진대사 점수에도 영향을 받을 것으로 생각한다. 이러한 경우에는 측정시간과 항신진대사 점수를 2개의 공변량으로 하여 이들을 제어한 약의 형태에 대한 효과를 공분산분석을 통해 알 수 있다.

공변량이 2개 이상인 경우에는 공변량이 하나인 경우의 모형을 그대로 확장해서 각 모수의 추정과 검정을 할 수 있다.

<공분산분석을 위한 두 가지 가정>

1) 각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일해야 한다. 즉, 회귀계수가 모든 약의 형태에 대해서 동일해야 하며, 교호작용이 없어야 한다.

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

2. 공분산분석에서의 검정

1) H0: β1 = 0 항신진대사의 효과가 없음

2) H0: β2 = 0 소요시간의 효과가 없음

이상 2개의 공변량 효과를 제어한 후 배출된 약물량의 모평균이 약의 형태type에 따라 차이가 있는가를 검정할 수 있는데, 이를 검정하기 위한 귀무가설은 아래와 같다.

3) H0 : α1 = α2 = ... = αI

1. library 호출 및 데이터 불러오기

> library(lsmeans)

> setwd('C:/Rwork')

> drug  = read.csv('약물배출량자료.csv')

> head(drug)

  관측번호 약형태 항신진대사점수 소요시간  약물량

1        1      1             37       61 11.3208

2        2      2             37       37 12.9151

3        3      3             45       53 18.8947

4        4      4             41       41 14.6739

5        5      5             57       41  8.6493

6        6      6             49       33  9.5238

> colname<-c('obs','type','x1','x2','y')

> colnames(drug)<-colname

> head(drug)

  obs type x1 x2       y

1   1    1 37 61 11.3208

2   2    2 37 37 12.9151

3   3    3 45 53 18.8947

4   4    4 41 41 14.6739

5   5    5 57 41  8.6493

6   6    6 49 33  9.5238

> attach(drug)

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

공분산 분석을 위한 중요한 가정으로 교호작용이 존재하는지 점검한다. 교호작용이 존재하면 공분산분석이 아닌 분산분석을 실시한다.

> model1 = aov(y~factor(type) + factor(type)*x1 + factor(type)*x2 , data = drug)

> summary(model1)

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

factor(type)     5  250.3    50.1   2.099   0.113   

x1               1  696.6   696.6  29.206 3.9e-05 ***

x2               1   54.4    54.4   2.282   0.148   

factor(type):x1  5  160.8    32.2   1.349   0.289   

factor(type):x2  5   42.3     8.5   0.355   0.872   

Residuals       18  429.3    23.9                   

---

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

-> 공변량이 2개 이므로 우변에 * 형태로 입력

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

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

p-value: x1 항신진대사점수 – 0.289

x2 소요시간 – 0.872

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

결론: 두 공변량(항신진대사점수, 소요시간)과 약의 형태 사이에 교호작용이 존재하지 않으므로 공분산분석을 할 수 있다. ( = 각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일하다.)

3. 이원공분산분석 (Two-way ANOVA)

> model2 = aov(y~factor(type) + x1 + x2, data=drug)

> summary(model2)

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

factor(type)  5  250.3    50.1   2.216   0.0808 . 

x1            1  696.6   696.6  30.837 6.13e-06 ***

x2            1   54.4    54.4   2.409   0.1319   

Residuals    28  632.5    22.6                    

---

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

-> 2개의 공변량을 보정하지 않았을 때의 F-값과 p값은 각각 2.216, 0.0808으로 유의수준 5%하에서 유의하지 않다.

귀무가설 H0 : α1 = α2 = ... = αI

대립가설 H1 : 최소한 하나 이상의 약은 효과가 있다.

p-value: 0.0808

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

결론: 약의 형태에 따라 배출된 약물량의 차이가 없다.

3-1. 공변량 효과 제어 시 치료법의 효과 검정

> model3 = lm(y~factor(type) + x1 + x2, data=drug)

> summary(model3)

Call:

lm(formula = y ~ factor(type) + x1 + x2, data = drug)

Residuals:

   Min     1Q Median     3Q    Max

-7.222 -2.634 -0.379  1.475  9.646

Coefficients:

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

(Intercept)    36.6370     8.9591   4.089 0.000331 ***

factor(type)2   0.8965     2.7823   0.322 0.749677   

factor(type)3   7.9097     2.7684   2.857 0.007973 **

factor(type)4   3.0722     2.8247   1.088 0.286025   

factor(type)5   9.5434     2.8534   3.345 0.002355 **

factor(type)6   5.8389     2.8391   2.057 0.049149 * 

x1             -0.7606     0.1375  -5.531 6.51e-06 ***

x2              0.1647     0.1061   1.552 0.131868   

---

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

Residual standard error: 4.753 on 28 degrees of freedom

Multiple R-squared:  0.6129,    Adjusted R-squared:  0.5161

F-statistic: 6.332 on 7 and 28 DF,  p-value: 0.0001641

귀무가설 H0 : β1 = β2 =  0  공변량(항신진대사, 소요시간)의 효과가 없음

        H0 : α1 = α2 = ... = αI

대립가설 not H0

p-value : 0.0001641

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

결론: 우리가 세운 모형의 자료에 적합하다는 것을 알 수 있다. 배출된 약물량의 모평균이 약 형태type의 효과와 공변량들의 효과가 없다라고 말할 수 없다.

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

1)제1종 제곱합(Type I SS): SS(type)

> model2 = aov(y~factor(type) + x1 + x2, data=drug)

> summary(model2)

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

factor(type)  5  250.3    50.1   2.216   0.0808 . 

x1            1  696.6   696.6  30.837 6.13e-06 ***

x2            1   54.4    54.4   2.409   0.1319   

Residuals    28  632.5    22.6                    

---

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

OR

> summary(aov(y~factor(type)+x1+x2,data=drug))

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

factor(type)  5  250.3    50.1   2.216   0.0808 . 

x1            1  696.6   696.6  30.837 6.13e-06 ***

x2            1   54.4    54.4   2.409   0.1319   

Residuals    28  632.5    22.6                    

---

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

-> 결과해석: 약의 형태type가 기여한 제1종 제공합의 p-value는 0.0808로 유의수준 5%하에서 유의하지 않다.

2)제3종 제곱합(Type III SS): SS(type | x1, x2)

> summary(aov(y~x1+x2+factor(type),data=drug))

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

x1            1  516.5   516.5  22.867 5.03e-05 ***

x2            1   62.8    62.8   2.779   0.1067   

factor(type)  5  422.0    84.4   3.736   0.0102 * 

Residuals    28  632.5    22.6                    

---

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

-> 결과해석: x1,2 given type, 공변량x1,x2가 기여한 상태에서 순수하게 약의 형태type가 기여한 제3종 제공합의 p-value는 0.0102로 유의수준 5%하에서 유의하다.

5. 각 수준 별 추정치를 알아보자. 

> model3 = lm(y~factor(type) + x1 + x2, data=drug)

> summary(model3)

Call:

lm(formula = y ~ factor(type) + x1 + x2, data = drug)

Residuals:

   Min     1Q Median     3Q    Max

-7.222 -2.634 -0.379  1.475  9.646

Coefficients:

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

(Intercept)    36.6370     8.9591   4.089 0.000331 ***

factor(type)2   0.8965     2.7823   0.322 0.749677   

factor(type)3   7.9097     2.7684   2.857 0.007973 **

factor(type)4   3.0722     2.8247   1.088 0.286025   

factor(type)5   9.5434     2.8534   3.345 0.002355 **

factor(type)6   5.8389     2.8391   2.057 0.049149 * 

x1             -0.7606     0.1375  -5.531 6.51e-06 ***

x2              0.1647     0.1061   1.552 0.131868   

---

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

Residual standard error: 4.753 on 28 degrees of freedom

Multiple R-squared:  0.6129,    Adjusted R-squared:  0.5161

F-statistic: 6.332 on 7 and 28 DF,  p-value: 0.0001641

-> 출력결과 유의하게 나타나는 것은 3,5,번째 약의 형태의 추정치들이 첫 번째 약의 형태와의 차이를 타나내는데 p-값이 유의수준 5% 하에서 유의하다. 항신진대사 점수(x1)는 p 값이 <0.0001로 유의수준 5%하에서 매우 유의하지만 약의 배출 소요시간(x2)은 p 값이 0.131868로 유위수준 5%하에서 유의하지 않은 것으로 나타났다. 위의 출력 결과를 바탕으로 약의 형태별 공분산모형식을 쓰면 다음과 같다.

ŷ_1j = β^₀ + α^₁ + β^1x_{11j  +} β^2x_21j  = 36.637 + 0 - 0.761 x_11j + 0.165x_21j

ŷ_2j = β^₀ + α^₂ + β^1x_{12j  +} β^2x_21j = 36.637 + 0.897- 0.761 x_12j + 0.165x_22j

ŷ_3j = β^₀ + α^₃ + β^1x_{13j  +} β^2x_21j = 36.637 + 7.908 - 0.761 x_13j + 0.165x_23j

ŷ_4j = β^₀ + α^₄ + β^1x_{14j  +} β^2x_21j = 36.637 + 3.072 - 0.761 x_14j + 0.165x_24j

ŷ_5j = β^₀ + α^₅ + β^1x_{15j  +} β^2x_21j = 36.637 +9.544 - 0.761 x_15j + 0.165x_25j

ŷ6_j = β^₀ + α^₆ + β^1x_{16j  +} β^2x_26j = 36.637 + 5.839 - 0.761 x_16j + 0.165x_26j

위의 식에 나타난 각 회귀식의 추정결과 모든 약의 형태에 대해서 항신진대사 점수(x1)와 약의 배출 소요시간(x2)의 회귀계수는 동일하다. 약의 배출 소요시간을 제어한 항신진대사 점수의 효과를 보면 항신진대사 점수가 1단위 높아짐에 따라 약의 배출량은 0.761만큼 감소하고 통계적으로 유의한 효과가 있고, 항신진대사 점수를 제어한 약의 배출 소요시간은 1단위 증가할수록 0.165만큼 증가하지만 그 효과는 유의하지 않다. 

 6. LSMEANS(Adjusted means) 계산

> lsmeans(model3,~type)

 type    lsmean       SE df  lower.CL  upper.CL

    1  5.674385 1.991087 28  1.595828  9.752942

    2  6.570912 1.946518 28  2.583650 10.558174

    3 13.584122 1.969292 28  9.550210 17.618034

    4  8.746633 1.955440 28  4.741095 12.752171

    5 15.217831 1.987611 28 11.146394 19.289268

    6 11.513284 1.980761 28  7.455879 15.570690

Confidence level used: 0.95

-> 각 약의 형태에 대한 LSMEAN(보정된 평균)값이 출력되어 있다. 보정된 평균은 항신진대사 점수x1와 소요시간x2의 효과를 제어했을 때 반응변수Y인 배출되는 약물량의 평균값이다. 앞의 식에서 각 약의 형태의 회귀식의 공변량값에 전체 평균값을 넣어서 약의 형태별 보정된 평균(adjusted mean)을 아래 식과 같이 계산할 수 있다.

y bar_ad1 = 36.387 + 0 – 0.761*51.22 + 0.165*48.56 = 5.67

y bar_ad2 = 36.387 + 0.897 – 0.761*51.22 + 0.165*48.56 = 6.57

y bar_ad3 = 36.387 + 7.908 – 0.761*51.22 + 0.165*48.56 = 13.59

y bar_ad4 = 36.387 + 3.072 – 0.761*51.22 + 0.165*48.56 = 8.75

y bar_ad5 = 36.387 + 9.544 – 0.761*51.22 + 0.165*48.56 = 15.22

y bar_ad6 = 36.387 + 5.839 – 0.761*51.22 + 0.165*48.56 = 11.51

보정된 평균값을 보면 첫 번째와 두 번째 약의 형태의 보정된 평균값이 다소 작고 세 번째, 다섯 번째, 여섯 번째 약의 형태에 대한 보정된 평균값이 크다는 것을 알 수 있다.

6. 처리 간 다중비교

> model3.lsm = lsmeans(model3,pairwise ~ type, glhargs = list())

> print(model3.lsm,omit=1)

$lsmeans

 type    lsmean       SE df  lower.CL  upper.CL

    1  5.674385 1.991087 28  1.595828  9.752942

    2  6.570912 1.946518 28  2.583650 10.558174

    3 13.584122 1.969292 28  9.550210 17.618034

    4  8.746633 1.955440 28  4.741095 12.752171

    5 15.217831 1.987611 28 11.146394 19.289268

    6 11.513284 1.980761 28  7.455879 15.570690

Confidence level used: 0.95

$contrasts

 contrast   estimate       SE df t.ratio p.value

 1 - 2    -0.8965269 2.782316 28  -0.322  0.9995

 1 - 3    -7.9097368 2.768397 28  -2.857  0.0771

 1 - 4    -3.0722481 2.824671 28  -1.088  0.8822

 1 - 5    -9.5434459 2.853395 28  -3.345  0.0257

 1 - 6    -5.8388992 2.839070 28  -2.057  0.3378

 2 - 3    -7.0132099 2.783087 28  -2.520  0.1525

 2 - 4    -2.1757212 2.753630 28  -0.790  0.9669

 2 - 5    -8.6469191 2.802553 28  -3.085  0.0468

 2 - 6    -4.9423723 2.758561 28  -1.792  0.4869

 3 - 4     4.8374887 2.801787 28   1.727  0.5266

 3 - 5    -1.6337091 2.782316 28  -0.587  0.9911

 3 - 6     2.0708376 2.840329 28   0.729  0.9766

 4 - 5    -6.4711978 2.783344 28  -2.325  0.2178

 4 - 6    -2.7666511 2.753890 28  -1.005  0.9125

 5 - 6     3.7045468 2.831753 28   1.308  0.7781

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

위의 결과를 보면 약의 형태 5의 경우 보정된 평균이 형태 1과 2와 유의하게 차이가 나는 것을 볼 수 있다.

> names(model3.lsm)

[1] "lsmeans"   "contrasts"

> plot(model3.lsm[[1]])

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

->Tukey(HSD) 검정 결과를 신뢰구간으로 보면, 1-5, 2-5가 유의함을 알 수 있다.

Dunnett vs Tukey

다중분석 하기 전 순서형인 type 변수를 명목형 변수로 변경

다중 분석 시작

> library(multcomp)

> tukey<-glht(model4,linfct=mcp(type='Tukey'))

> summary(tukey)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lm(formula = y ~ type + x1 + x2, data = drug)

Linear Hypotheses:

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

2 - 1 == 0   0.8965     2.7823   0.322   0.9995 

3 - 1 == 0   7.9097     2.7684   2.857   0.0771 .

4 - 1 == 0   3.0722     2.8247   1.088   0.8821 

5 - 1 == 0   9.5434     2.8534   3.345   0.0257 *

6 - 1 == 0   5.8389     2.8391   2.057   0.3374 

3 - 2 == 0   7.0132     2.7831   2.520   0.1524 

4 - 2 == 0   2.1757     2.7536   0.790   0.9669 

5 - 2 == 0   8.6469     2.8026   3.085   0.0469 *

6 - 2 == 0   4.9424     2.7586   1.792   0.4867 

4 - 3 == 0  -4.8375     2.8018  -1.727   0.5265 

5 - 3 == 0   1.6337     2.7823   0.587   0.9911 

6 - 3 == 0  -2.0708     2.8403  -0.729   0.9766 

5 - 4 == 0   6.4712     2.7833   2.325   0.2178 

6 - 4 == 0   2.7667     2.7539   1.005   0.9124 

6 - 5 == 0  -3.7045     2.8318  -1.308   0.7780 

---

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

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

> plot(tukey)

->Tukey(HSD) 검정 결과를 신뢰구간으로 보면, 1-5, 2-5가 유의함을 알 수 있다.

> dunnett <- glht(model4,linfct=mcp(type='Dunnett'))

> summary(dunnett)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts

Fit: lm(formula = y ~ type + x1 + x2, data = drug)

Linear Hypotheses:

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

2 - 1 == 0   0.8965     2.7823   0.322   0.9974 

3 - 1 == 0   7.9097     2.7684   2.857   0.0324 *

4 - 1 == 0   3.0722     2.8247   1.088   0.7130 

5 - 1 == 0   9.5434     2.8534   3.345   0.0103 *

6 - 1 == 0   5.8389     2.8391   2.057   0.1733 

---

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

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

> plot(dunnett)

7. 모형 적합성 검토

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

> plot(model3)

-> 등분산성, 정규성 가정에는 큰 문제가 없음을 알 수 있다.

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

목표변수가 집단을 의미하는 범주형 의사결정나무 -> 분류나무모형

목표변수가 연속형 변수인 의사결정나무 -> 회귀나무모형

 목표 변수가 연속형 변수인 경우, 의사결정나무는 회귀나무라고 한다. 통계학에서 목표변수가 연속형일 때 흔히 수행하는 분석이 선형 회귀분석인 것을 연상해 본다면 명칭이 이해가 될 것이다. 일반적으로 범주형 변수들은 가변수화하여 0,1의 형태로 변형하여 회귀분석에 활용하지만, 범주형 변수의 수가 매우 많고, 각 범주형 변수의 개수도 많은 경우 해석은 어려워지고 분석도 복잡해 진다. 이러한 경우 회귀나무 모형을 사용하게 되면, 가변수를 생성할 필요 없이 범주형 입력 변수와 연속형 입력변수를 그대로 활용할 수 있게 되어 분석 및 그 해석이 용이해질 수 있다.

예제) 목표변수가 숫자인 보스턴하우징 데이터를 이용하여 cart 방법을 이용한 회귀나무 모형 구축

1) 데이터 읽기

> library(MASS)

> Boston$chas=factor(Boston$chas)

> Boston$rad=factor(Boston$rad)

> summary(Boston)

> Boston<-Boston[,-15]
> str(Boston)
'data.frame':   506 obs. of  14 variables:
 $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
 $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
 $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
 $ chas   : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
 $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
 $ rm     : num  6.58 6.42 7.18 7 7.15 ...
 $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
 $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
 $ rad    : Factor w/ 9 levels "1","2","3","4",..: 1 2 2 3 3 3 5 5 5 5 ...
 $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
 $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
 $ black  : num  397 397 393 395 397 ...
 $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
 $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ... 

2) cart 회귀나무모형 실행

> library(rpart)

> my.control<-rpart.control(xval=10,cp=0,minsplit=nrow(Boston)*0.05)

> fit.Boston=rpart(medv~.,data=Boston,method='anova',control=my.control)

> fit.Boston

n= 506

node), split, n, deviance, yval

* denotes terminal node

1) root 506 42716.30000 22.532810

   … 중간 생략…

  30) rad=2,7,8 10   140.04100 43.670000 *

     31) rad=1,3,4,5,24 17   135.72120 48.158820 *

함수 설명:

rpart.control(xval=10,cp=0,minsplit=nrow(Boston)*0.05) : 10-fold교차타당성 오분류율을 계산하고, 가지치기 방법에서 a값은 0이될때까지 순차적인 나무구조를 저장하도록 한다. 그리고 노드를 나누는 최소 자료이 수인 minsplit은 전체 데이터 수의 5% 수준에서 정의한다

3) 가지치기 수행 : ①최소 cp 값 찾기 

> printcp(fit.Boston)

Regression tree:
rpart(formula = medv ~ ., data = Boston, method = "anova", control = my.control)

Variables actually used in tree construction:
[1] age     crim    indus   lstat   nox     ptratio rad     rm    

Root node error: 42716/506 = 84.42

n= 506

           CP nsplit rel error  xerror     xstd
1  0.45274420      0   1.00000 1.00298 0.083114
2  0.17117244      1   0.54726 0.62762 0.058631
3  0.07165784      2   0.37608 0.43188 0.047879
4  0.03428819      3   0.30443 0.35161 0.043496
5  0.02661300      4   0.27014 0.32766 0.042850
6  0.01802372      5   0.24352 0.29534 0.041781
7  0.01348721      6   0.22550 0.27799 0.039217
8  0.01285085      7   0.21201 0.28299 0.039277
9  0.00844925      8   0.19916 0.26415 0.032335
10 0.00833821      9   0.19071 0.25667 0.031357
11 0.00726539     10   0.18238 0.25624 0.031407
12 0.00612633     11   0.17511 0.25343 0.031377
13 0.00480532     12   0.16898 0.24000 0.028645
14 0.00410785     13   0.16418 0.23897 0.027439
15 0.00394102     14   0.16007 0.23177 0.027282
16 0.00385379     15   0.15613 0.23177 0.027282
17 0.00223540     16   0.15228 0.22527 0.027268
18 0.00171691  17   0.15004 0.22441 0.027196
19 0.00153485     18   0.14832 0.22453 0.026980
20 0.00140981     19   0.14679 0.22585 0.026952
21 0.00135401     20   0.14538 0.22646 0.026946
22 0.00113725     22   0.14267 0.22674 0.026936
23 0.00098921     23   0.14153 0.22490 0.026948
24 0.00081924     24   0.14054 0.22674 0.027098
25 0.00074570     25   0.13972 0.22624 0.027097
26 0.00074096     27   0.13823 0.22671 0.027092
27 0.00040763     28   0.13749 0.22753 0.027088
28 0.00031464     29   0.13708 0.22749 0.027081
29 0.00019358     30   0.13677 0.22798 0.027098
30 0.00000000     31   0.13658 0.22816 0.027099

> names(fit.Boston)

[1] "frame"               "where"               "call"                "terms"              

[5] "cptable"             "method"              "parms"               "control"           

[9] "functions"           "numresp"             "splits"              "csplit"            

[13] "variable.importance" "y"                   "ordered"            

> which.min(fit.Boston$cp[,4])
18
18
> which.min(fit.Boston$cptable[,'xerror'])
18
18

> fit.Boston$cp[18,]
          CP       nsplit    rel error       xerror         xstd
 0.001716908 17.000000000  0.150039975  0.224411843  0.027196436
> fit.Boston$cp[18]
[1] 0.001716908

> fit.Boston$cptable[which.min(fit.Boston$cptable[,'xerror']),]
          CP       nsplit    rel error       xerror         xstd
 0.001716908 17.000000000  0.150039975  0.224411843  0.027196436
> fit.Boston$cptable[which.min(fit.Boston$cptable[,'xerror'])]
[1] 0.001716908

> 0.22441 +0.027196
[1] 0.251606

결과 해석: 가장 작은 cp값은 18번째 0.001716908이나 가장 작은 xerror의 xstd 범위 내 즉, 0.22441 + 0.027196 = 0.251606 내 작은 cp 값 적용도 가능하다.

3) 가지치기 수행 : ②가장 적은 cp값을 찾았으니 prune 함수를 이용하여 가지기치를 수행하자. 

> fit.prune.Boston<-prune(fit.Boston,cp=0.00304027)

> print(fit.prune.Boston)

> names(fit.prune.Boston)

[1] "frame"               "where"               "call"                "terms"             

[5] "cptable"             "method"              "parms"               "control"           

[9] "functions"           "numresp"             "splits"              "csplit"            

[13] "variable.importance" "y"                   "ordered"           

> fit.prune.Boston$variable.importance
              rm             lstat           indus             nox             age         ptratio             dis
23124.34862 16998.05183  5871.44959  5492.56562  5042.88012  4757.36213  4742.25575
            tax            crim              zn             rad           black            chas
 4155.77974  2258.47799  1760.50488  1390.73661   772.28149    11.40364

node), split, n, deviance, yval 는 각각 노드번호, 분할규칙, 해등 노드의 총 관칙치 수, 해당 노드의 분산, 해당 노드의 목표변수 예측치(회귀나무일 때는 평균값)의 순서로 정보를 제공하고 있고, 줄 마지막에 있는 * denotes terminal node는 * 표시가 최종노드임을 알려주고 있다.

fit.prune.Boston$variable.importance로 중요 변수도 알 수 있다. rm이 가장 중요하고 상대적으로 chas가 덜 중요한 변수이다.

4) cart 회귀나무 모형 그리기

> plot(fit.prune.Boston)