관측값이 질적 자료(qualitative data) 또는 어떤 속성에 따라 분류되어 범주(category)에 속하는 도수(frequency)로 주어질 경우 이를 범주형 자료(categorical data) 라고 합니다. 범주형 자료의 예로는 학력(초등졸, 중등졸, 고등졸, 대졸, 대학원졸), 연수익(극빈, 하, 중, 상, 극상) 등이 있습니다.

 

앞서의 포스팅에서는 종속변수가 연속형 자료(continuous data)인 경우에 사용하는 검정 방법으로 t-Test와 ANOVA에 대해서 소개하였습니다.

 

이번 포스팅부터는 종속변수가 범주형 자료(categorical data)인 경우에 사용하는 분석기법으로 카이제곱 검정(Chi-Squared Test)에 대해서 알아보도록 하겠습니다.

 

범주형 자료 분석은 크게 적합도 검정(goodness f fit test), 독립성 검정(test of independence), 동질성 검정(test of homogeneity)의 3가지로 분류할 수 있으며, 이번 포스팅에서는 (1) 적합도 검정에 대해서 알아보도록 하겠습니다.

 

 

 

[ 범주형 자료 분석 (categorical data test) ]

 

(1) 적합도 검정(goodness of fit test) : 관측값들이 어떤 이론적 분포를 따르고 있는지를 검정. 한 개의 요인을 대상으로 함 

 

(2) 독립성 검정(test of independence) : 서로 다른 요인들에 의해 분할되어 있는 경우 그 요인들이 관찰값에 영향을 주고 있는지 아닌지, 요인들이 서로 연관이 있는지 없는지를 검정. 두 개의 요인을 대상으로 함.

 

(3) 동질성 검정(test of homogeneity) : 관측값들이 정해진 범주 내에서 서로 비슷하게 나타나고 있는지를 검정. 속성 A, B를 가진 부모집단(subpopulation) 각각으로부터 정해진 표본의 크기만큼 자료를 추출하는 경우에 분할표에서 부모집단의 비율이 동일한가를 검정. 두 개의 요인을 대상으로 함.

 

 

 

 

 

 

적합도 검정(goodness of fit test)은 k개의 범주 (혹은 계급)을 가지는 한 개의 요인(factor)에 대해서 어떤 이론적 분포를 따르고 있는지를 검정하는 방법입니다. 

 

기본 원리는, 도수분포의 각 구간에 있는 관측도수를 O1, O2, ..., Ok 라 하고, 각 범주 ( 혹은 계급)가 일어날 확률을 p1, p2, ..., pk 라고 할 때 기대되는 관측도수 E1, E2, ..., Ek 를 계산하여 실제 관측도수와 기대 관측도수의 차이를 카이제곱 검정 통계량(Chi-squared statistics)을 활용하여 가정한 확률모형에 적합한지를 평가하게 됩니다. 만약 귀무가설 H0가 맞다면 관측도수와 기대도수가 별 차이가 없을 것이므로 검정통계량 X0^2 값이 작을 것이며, 반대로 대립가설 H1이 맞다면 관측도수와 기대도수의 차이가 클 것이므로 검정통계량 X0^2 값이 커질 것입니다.

 

 

 

[ 적합도 검정 가설 및 검정 통계량, 검정 방법 ]

(1) 가설

 

  - 귀무가설 H0 : 관측값의 도수와 가정한 이론 도수(기대 관측도수)가 동일하다
                        ( p1 = p10, p2 = p20, ..., pk = pko )

 

  - 대립가설 H1 : 적어도 하나의 범주 (혹은 계급)의 도수가 가정한 이론 도수(기대 관측도수)와 다르다

                        (적어도 하나의 pi는 가정된 값 pi0과 다르다)

 

 

(2) 검정 통계량

 

 

 

(3) 검정 방법

 

 

 

 

  • (a) chisq.test() of data from text problem

 

아래 문제에 대해서 R의 chisq.test() 함수를 사용해서 풀어보도록 하겠습니다.

 

 

 

(문제)  유전학자 멘델은 콩 교배에 대한 유전의 이론적 모형으로서 잡종비율을 A : B : C = 2 : 3 : 5 라고 주장하였다.  이 이론의 진위를 가리기 위해 두 콩 종자의 교배로 나타난 100개의 콩을 조사하였더니 A형 19개, B형 41개, C형 40개였다.  이러한 관찰값을 얻었을 때 멘델 유전학자의 이론이 맞다고 할 수 있는지를 유의수준 α = 0.05 에서 검정하여라.

 

 

 

> ##---------------------------------------------------------------------
> ## categorical data analysis - (1) goodness of fit test : chisq.test()
> ##---------------------------------------------------------------------
> 
> obs <- c(19, 41, 40)
> null.probs <- c(2/10, 3/10, 5/10)
> 
> chisq.test(obs, p=null.probs)

	Chi-squared test for given probabilities

data:  obs
X-squared = 6.0833, df = 2, p-value = 0.04776

 

 

위 분석결과를 보면 P-value가 0.04776 이므로 유의수준 α 0.05보다 작으므로 귀무가설 H0를 기각하고 대립가설 H1을 채택하여 "멘델이 주장한 콩의 잡종비율 이론적 분포는 적합하지 않다"고 판단할 수 있겠습니다.

 

 

 참고로, R로 통계분석을 하면 콘솔 창에 보여지는 내용 말고도 실제로 다양한 통계량이 계산이 되어 list 형태로 메모리상에 가지고 있으며 단지 눈에 보이지 않을 뿐인데요, indexing 기법을 활용하면 chisq.test() 함수로 카이제곱 검정 실행한 후에 다양한 통계량들을 선별해서 볼 수도 있고, 통계량을 다른 분석 혹은 애플리케이션에 input으로 넣어 재활용할 수도 있습니다.  R의 큰 장점 중의 하나이니 팁으로 알아두면 좋겠습니다. 주요 통계량 몇 개를 아래에 소개합니다.

 

 

> # To see results of chisquared test
> chisq.test_output_1 <- chisq.test(obs, p=null.probs)
> 
> chisq.test_output_1$observed # observed frequency
[1] 19 41 40
> chisq.test_output_1$expected # expected frequeycy
[1] 20 30 50
> chisq.test_output_1$residuals # residual between observed and expected frequecy
[1] -0.2236068  2.0083160 -1.4142136
> 
> chisq.test_output_1$statistic # chi-squared statistics
X-squared 
 6.083333 
> chisq.test_output_1$parameter # degrees of freedom
df 
 2 
> chisq.test_output_1$p.value # P-value
[1] 0.04775523

 

 

 

참고로 하나더, 검정통계량 X^2는 귀무가설 H0가 참이라는 가정 하에 근사적으로 자유도가 k-1 인 카이제곱분포를 따르는 것으로 알려져 있습니다.  (☞ 카이제곱 분포(Chi-squared distribution) 참고 포스팅)

 

카이제곱 검정에 사용하는 카이제곱 분포는 범주형 자료의 도수 추정에 사용되는데요, 이때 도수가 너무 작으면 "카이제곱 approximation)이 정확하지 않을 수도 있습니다" 라는 경고메시지가 뜹니다.

 

 

> # Warning message when there are not sufficient frequencies > # R will issue a warning message if any of the EFs fall below 5

> obs_2 <- c(5, 5) > null.probs_2 <- c(0.3, 0.7) > > chisq.test(obs_2, p=null.probs_2) Chi-squared test for given probabilities data: obs_2 X-squared = 1.9048, df = 1, p-value = 0.1675 Warning message: In chisq.test(obs_2, p = null.probs_2) : 카이제곱 approximation은 정확하지 않을수도 있습니다

 

 

 

 

  • (b) chisqtest() of data from a table object

위의 문제는 텍스트 문제로 도수와 확률이 주어지면 'x'와 'p'를 직접 입력하였는데요, 데이터가 Table object 형태로 주어졌을 때 카이제곱 검정으로 적합도 검정하는 방법을 소개하겠습니다.

 

R에 내장된 HairEyeColor table 데이터셋에 있는 Hair 요인 변수를 대상으로, Hair의 요인 수준(factor levels, 혹은 계급 class) 별로 생물학자가 주장하기를 확률이 Black 20%, Brown 50%, Red 10%, Blond 20% 라고 하는데요, 유의수준 α 0.05 로 검정을 해보겠습니다.

 

> ##------------
> # chisq.test() of data from a table object
> str(HairEyeColor) 
 table [1:4, 1:4, 1:2] 32 53 10 3 11 50 10 30 10 25 ...
 - attr(*, "dimnames")=List of 3
  ..$ Hair: chr [1:4] "Black" "Brown" "Red" "Blond"
  ..$ Eye : chr [1:4] "Brown" "Blue" "Hazel" "Green"
  ..$ Sex : chr [1:2] "Male" "Female"
> 
> HairEyeColor 
, , Sex = Male

       Eye
Hair    Brown Blue Hazel Green
  Black    32   11    10     3
  Brown    53   50    25    15
  Red      10   10     7     7
  Blond     3   30     5     8

, , Sex = Female

       Eye
Hair    Brown Blue Hazel Green
  Black    36    9     5     2
  Brown    66   34    29    14
  Red      16    7     7     7
  Blond     4   64     5     8

> 
> dimnames(HairEyeColor)
$Hair
[1] "Black" "Brown" "Red"   "Blond"

$Eye
[1] "Brown" "Blue"  "Hazel" "Green"

$Sex
[1] "Male"   "Female"

> 
> margin.table(HairEyeColor, 1) # Hair
Hair
Black Brown   Red Blond 
  108   286    71   127 
> margin.table(HairEyeColor, 2) # Eye
Eye
Brown  Blue Hazel Green 
  220   215    93    64 
> margin.table(HairEyeColor, 3) # Sex
Sex
  Male Female 
   279    313 
> 
> # vector of observed frequencies and probabilities
> Hair_Freq <- c(margin.table(HairEyeColor, 1))
> Hair_Freq
Black Brown   Red Blond 
  108   286    71   127 
> Hair_Prob <- c(0.2, 0.5, 0.1, 0.2)
> 
> chisq.test(x=Hair_Freq, p=Hair_Prob)

	Chi-squared test for given probabilities

data:  Hair_Freq
X-squared = 4.228, df = 3, p-value = 0.2379

 

 

 

  • (c) chisq.test() of data from Data Frame

 

데이터가 텍스트, Table object가 아니고 Data Frame일 경우에 카이제급 검정하는 방법도 소개해드리겠습니다.  MASS 패키지에 내장된 Cars93 Data Frame 을 이용하며, 자동차종(Type)의 이론상 분포 Compact 20%,  Large 10%, Midsize 20%, Small 20%, Sporty  20%, Van 10% 에 대해서 유의수준 α 0.05로 검정해보겠습니다. 

data(data frame, package="xxx"), table() 함수를 이용합니다.

 

> ##------------
> # chisq.test() of data from data frame
> data(Cars93, package="MASS")
> head(Cars93)
  Manufacturer   Model    Type Min.Price Price Max.Price MPG.city MPG.highway            AirBags
1        Acura Integra   Small      12.9  15.9      18.8       25          31               None
2        Acura  Legend Midsize      29.2  33.9      38.7       18          25 Driver & Passenger
3         Audi      90 Compact      25.9  29.1      32.3       20          26        Driver only
4         Audi     100 Midsize      30.8  37.7      44.6       19          26 Driver & Passenger
5          BMW    535i Midsize      23.7  30.0      36.2       22          30        Driver only
6        Buick Century Midsize      14.2  15.7      17.3       22          31        Driver only
  DriveTrain Cylinders EngineSize Horsepower  RPM Rev.per.mile Man.trans.avail Fuel.tank.capacity
1      Front         4        1.8        140 6300         2890             Yes               13.2
2      Front         6        3.2        200 5500         2335             Yes               18.0
3      Front         6        2.8        172 5500         2280             Yes               16.9
4      Front         6        2.8        172 5500         2535             Yes               21.1
5       Rear         4        3.5        208 5700         2545             Yes               21.1
6      Front         4        2.2        110 5200         2565              No               16.4
  Passengers Length Wheelbase Width Turn.circle Rear.seat.room Luggage.room Weight  Origin
1          5    177       102    68          37           26.5           11   2705 non-USA
2          5    195       115    71          38           30.0           15   3560 non-USA
3          5    180       102    67          37           28.0           14   3375 non-USA
4          6    193       106    70          37           31.0           17   3405 non-USA
5          4    186       109    69          39           27.0           13   3640 non-USA
6          6    189       105    69          41           28.0           16   2880     USA
           Make
1 Acura Integra
2  Acura Legend
3       Audi 90
4      Audi 100
5      BMW 535i
6 Buick Century
> str(Cars93)
'data.frame':	93 obs. of  27 variables:
 $ Manufacturer      : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4 4 4 5 ...
 $ Model             : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1 6 24 54 74 73 35 ...
 $ Type              : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3 2 2 3 2 ...
 $ Min.Price         : num  12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ...
 $ Price             : num  15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ...
 $ Max.Price         : num  18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3 36.3 ...
 $ MPG.city          : int  25 18 20 19 22 22 19 16 19 16 ...
 $ MPG.highway       : int  31 25 26 26 30 31 28 25 27 25 ...
 $ AirBags           : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2 2 2 2 2 2 ...
 $ DriveTrain        : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3 2 2 ...
 $ Cylinders         : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4 4 4 5 ...
 $ EngineSize        : num  1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ...
 $ Horsepower        : int  140 200 172 172 208 110 170 180 170 200 ...
 $ RPM               : int  6300 5500 5500 5500 5700 5200 4800 4000 4800 4100 ...
 $ Rev.per.mile      : int  2890 2335 2280 2535 2545 2565 1570 1320 1690 1510 ...
 $ Man.trans.avail   : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1 1 ...
 $ Fuel.tank.capacity: num  13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ...
 $ Passengers        : int  5 5 5 6 4 6 6 6 5 6 ...
 $ Length            : int  177 195 180 193 186 189 200 216 198 206 ...
 $ Wheelbase         : int  102 115 102 106 109 105 111 116 108 114 ...
 $ Width             : int  68 71 67 70 69 69 74 78 73 73 ...
 $ Turn.circle       : int  37 38 37 37 39 41 42 45 41 43 ...
 $ Rear.seat.room    : num  26.5 30 28 31 27 28 30.5 30.5 26.5 35 ...
 $ Luggage.room      : int  11 15 14 17 13 16 17 21 14 18 ...
 $ Weight            : int  2705 3560 3375 3405 3640 2880 3470 4105 3495 3620 ...
 $ Origin            : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1 1 1 1 ...
 $ Make              : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ...
> 
> Car_Type <- table(Cars93$Type)
> Car_Type

Compact   Large Midsize   Small  Sporty     Van 
     16      11      22      21      14       9 
> Car_Type_Prob <- c(0.2, 0.1, 0.2, 0.2, 0.2, 0.1)
> 
> chisq.test(x=Car_Type, p=Car_Type_Prob)

	Chi-squared test for given probabilities

data:  Car_Type
X-squared = 2.7527, df = 5, p-value = 0.738

 

 

 

이상으로 다양한 형태의 데이터셋을 활용해서 적합도 검정(goodness of fit test)을 카이제곱 검정 기법을 사용해서 하는 방법을 소개하였습니다.

 

=> 카이제곱 적합도 검정으로 관측치가 포아송분포로 부터의 데이터인지 여부를 검정하는 예시는 http://rfriend.tistory.com/362 를 참고하세요.

 

다음번 포스팅에서는 (2) 독립성 검정(test of independence) , (3) 동질성 검정 (test of homogeneity)대해서 알아보도록 하겠습니다.

 

많은 도움이 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

 

728x90
반응형
Posted by Rfriend
,

2개의 모집단에 대한 평균을 비교, 분석하는 통계적 기법으로 t-Test를 활용하였다면, 비교하고자 하는 집단이 3개 이상일 경우에는 분산분석 (ANOVA : Analysis Of Variance)를 이용합니다. 

 

설명변수는 범주형 자료(categorical data)이어야 하며, 종속변수는 연속형 자료(continuous data) 일 때 3개 이상 집단 간 평균 비교분석에 분산분석(ANOVA) 을 사용하게 됩니다.

 

분산분석(ANOVA)은 기본적으로 분산의 개념을 이용하여 분석하는 방법으로서, 분산을 계산할 때처럼 편차의 각각의 제곱합을 해당 자유도로 나누어서 얻게 되는 값을 이용하여 수준평균들간의 차이가 존재하는 지를 판단하게 됩니다. 

 

지난번 포스팅에서 이원분산분석(two-way ANOVA) 중에서 (1) 관측값이 하나일 경우 (one observation in each cell)의 이원분산분석 에 대해서 알아보았다면, 이번 포스팅에서는 (2) 관측값이 두개 이상일 경우(more than 2 observations in each cell)의 이원분산분석에 대해서 알아보겠습니다.  단, 이때 각 집단 내 관측값의 개수는 동일합니다.

 

관측값이 두개 이상일 경우의 이원분산분석에서는 두 개의 요인 수준별 주 효과(main effect)와 더불어 두 요인이 서로 상호간에 영향을 주고 받으면서 나타나는 반응효과인 교호작용 효과(interaction effect)를 추가로 분석하는 것이 관측값이 하나인 경우와의 차이점입니다.

 

 

[ 관측값의 개수에 따른 이원분산분석의 분석 대상 ]

 

 

 

 

요인 A의 i번째 수준과 요인 B의 j번째 수준에서 측정된 k번째 관측값을 Yijk라고 할 때 데이터셋의 형태는 아래와 같습니다. 

 

 

[ 이원분산분석 - 데이터 형태 ]

(dataset for two-way ANOVA, factor A & B, factor levels a & b, observations k)]

 

 

 

 

요인 A의 i번째 수준과 요인 B의 j번째 수준에서 측정된 k번째 관측값을 Yijk는 요인 A와 B의 주 효과(main effect)와 함께 두 요인의 교호작용 효과(interaction effect)이 존재하게 되며, 편차 는 아래와 같이 4개의 성분합으로 나타낼 수 있습니다.

 

 

[ 편차의 4개 성분합으로의 분해 ]

 


 

 

위의 식의 양변을 제곱하여 모든 관측값에 대하여 더한 후에 정리를 하면

 

SST = SSA + SSB + SSAB + SSE 

 

의 관계식을 얻게 됩니다.

 

이를 요인과 교호작용별로 제곱합, 자유도, 평균제곱, F통계량을 알기 쉽도록 정리한 이원분산분석표는 아래와 같습니다.

 

[ 이원분산분석표 (two-way ANOVA table) ]

 

 요인

제곱합

(squared sum)

자유도

(degrees of freedom) 

평균제곱

(mean squared) 

F statistics 

 요인 A

(factor A)

 SSA

a-1

MSA

MSA/MSE

 요인 B

(factor B)

 SSB

b-1

MSB

MSB/MSE

 교호작용

(interaction effect of A, B)

 SSAB

 (a-1)(b-1)

MSAB

MSAB/MSE 

 오차

(error)

 SSE

ab(n-1)

MSE

 

 계

(total)

 SST

 nab-1

 

 

 * SST : Total Sum of Squares,  SSE : Error Sum of Squares

 

 

검정통계량으로는 F 통계량(요인A 효과 검정 = MSA/MSE, 요인B 효과 검정 = MSB/MSE, 교호작용 효과 검정 = MSAB/MSE)을 사용하며, 요인 A의 효과와 요인 B, 교호작용 효과에 대한 검정 절차 및 방법은 아래와 같습니다.

 

[ 이원분산분석 검정 방법 ]

 

 

(1) 교호작용 효과 (interaction effect)

   - 귀무가설 H0 : (αβ)ij = 0

   - 대립가설 H1 : 모든 (αβ)ij는 0이 아니다

   - 검정통계량 : F0 = MSAB/MSE

   - 판정 : P-value가 유의수준(α)보다 작으면 귀무가설 H0를 기각하고, 대립가설 H1을 채택

 

 

(2) 요인 A의 주 효과 (main effect of factor A)

   - 귀무가설 H0 : α1 = ... = αa = 0

   - 대립가설 H1 : 모든 αi 는 0이 아니다

   - 검정통계량 : F0 = MSA/MSE

   - 판정 : P-value가 유의수준(α)보다 작으면 귀무가설 H0를 기각하고, 대립가설 H1을 채택

 

 

(3) 요인 B의 주 효과 (main effect of factor B)

   - 귀무가설 H0 : α1 = ... = αb = 0

   - 대립가설 H1 : 모든 αj 는 0이 아니다

   - 검정통계량 : F0 = MSB/MSE

   - 판정 : P-value가 유의수준(α)보다 작으면 귀무가설 H0를 기각하고, 대립가설 H1을 채택

 

 

 

 

 

그럼, 아래 예제에 대해서 R의 aov() 함수를 활용해서 문제를 풀어보도록 하겠습니다.

 

 

(문제)  학급과 성별에 따른 통계학 성적이 아래와 같다고 할 때(각 cell 별 3명의 학생 성적 측정), 유의수준 α = 0.05 에서 학급과 성별 요인의 주효과와 교호작용효과에 대해 검정하시오.

 

[학급과 성별에 따른 통계학 점수표]

 

         Factor B (class)

 

Factor A (gender)

Class 1 

Class 2 

Class 3 

Mean 

남성 (M) 

 71, 77, 78

76, 77, 78 

71, 70, 69 

 75.33

 여성 (F)

 76, 76, 80

79, 78, 77

71, 71, 70 

 74.11

 Mean

76.33 

77.50 

70.33 

74.72 

 

 

 

아래에 요약통계량과 교호작용 여부를 가늠해볼 수 있는 그래프 생성 R 함수도 일부 추가하였습니다. 만약 교호작용효과 그래프(Interaction Effect Plot) 가 서로 교차를 하면 교호작용이 있다고 보면 되며, 서로 평행선을 이룬다면 교호작용이 없다고 해석하면 되겠습니다.  (교호작용도는 엑셀로 치면 "꺽은선형 차트"가 되겠습니다)

 

 

[ R aov() 함수 사용 방법 ]

 

 
> # (2) two-way ANOVA when there are more than one observation per cell (different treatment groups)
> #     (the number of observations in each cell must be equal)
> 
> gender.fac <- as.factor(c(rep("M", 9), rep("F", 9)))
> gender.fac
 [1] M M M M M M M M M F F F F F F F F F
Levels: F M
> 
> class <- c("class_1", "class_1", "class_1", "class_2", "class_2", "class_2", "class_3", "class_3", "class_3")
> class.fac <- as.factor(c(rep(class, 2)))
> class.fac
 [1] class_1 class_1 class_1 class_2 class_2 class_2 class_3 class_3 class_3 class_1 class_1
[12] class_1 class_2 class_2 class_2 class_3 class_3 class_3
Levels: class_1 class_2 class_3
> 
> score_stats <- c(71, 77, 78, 76, 77, 78, 71, 70, 69, 80, 76, 80, 79, 78, 77, 73, 71, 70)
> 
> 
> # summary statistics
> score.df <- data.frame(gender.fac, class.fac, score_stats)
> score.df
   gender.fac class.fac score_stats
1           M   class_1          71
2           M   class_1          77
3           M   class_1          78
4           M   class_2          76
5           M   class_2          77
6           M   class_2          78
7           M   class_3          71
8           M   class_3          70
9           M   class_3          69
10          F   class_1          80
11          F   class_1          76
12          F   class_1          80
13          F   class_2          79
14          F   class_2          78
15          F   class_2          77
16          F   class_3          73
17          F   class_3          71
18          F   class_3          70
> 
> install.packages("doBy")

Installing package into ‘C:/Users/user/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.2/doBy_4.5-13.zip'
Content type 'application/zip' length 3431466 bytes (3.3 MB)
downloaded 3.3 MB

package ‘doBy’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpWOKfji\downloaded_packages

 

 

> library(doBy)
> summaryBy(score_stats ~ gender.fac, data=score.df, FUN = c(mean, sd, min, max))
  gender.fac score_stats.mean score_stats.sd score_stats.min score_stats.max
1          F         76.00000       3.807887              70              80
2          M         74.11111       3.756476              69              78
> summaryBy(score_stats ~ class.fac, data=score.df, FUN = c(mean, sd, min, max))
  class.fac score_stats.mean score_stats.sd score_stats.min score_stats.max
1   class_1         77.00000       3.346640              71              80
2   class_2         77.50000       1.048809              76              79
3   class_3         70.66667       1.366260              69              73
> summary(score_stats, data=score.df)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  69.00   71.00   76.50   75.06   78.00   80.00 
> 
 
> # box plot and interaction plot
> par(mfrow = c(2, 2))
> plot(score_stats ~ gender.fac, main="box plot by gender")
> plot(score_stats ~ class.fac, main="box plot by class")
> interaction.plot(gender.fac, class.fac, score_stats, bty='l', main="interaction effect plot")
> interaction.plot(class.fac, gender.fac, score_stats, bty='l', main="interaction effect plot")
> 

 


> 
> # two-way ANOVA : aov()
> # replicates, interaction effect
> aov_model = aov(score_stats ~ gender.fac + class.fac + gender.fac:class.fac)
> summary(aov_model)
                     Df Sum Sq Mean Sq F value   Pr(>F)    
gender.fac            1  16.06   16.06   3.853 0.073249 .  
class.fac             2 174.11   87.06  20.893 0.000123 ***
gender.fac:class.fac  2   4.78    2.39   0.573 0.578354    
Residuals            12  50.00    4.17                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

 

 

R 분석 결과 중에서 제일 아래에 있는 이원분산분석표를 가지고 해석을 해보겠습니다.  먼저 (1) 교호작용이 있는지를 살펴보면, 성별과 학급별에 따른 교호작용효과에 대한 P-value가 0.578354로서 유의수준 α 0.05보다 크므로 귀무가설 H0를 채택하여 "성별과 학급의 교호작용에 의한 효과는 없다"고 판단할 수 있습니다.

 

두번째로 (2) 성별 요인(factor A levels of "M", "F")에 대한 P-value는 0.073249로서 유의수준 α 0.05 보다 크므로 귀무가설 H0를 채택하게 되어 "성별에 따른 통계학 성적 차이는 없다"고 판단할 수 있습니다.

 

세번째로 (3) 학급 요인(factor B levels of "class_1", "class_2", "class_3")에 따른 P-value는 0.000123으로서 유의수준 α 0.05보다 작으므로 귀무가설 H0를 기각하고 대립가설 H1을 채택하여 "학급에 따른 통계학 성적의 차이가 있다"고 판단할 수 있습니다.

 

 

학급 요인에 대해서는 유의수준 α 0.05에서 차이가 있다고 나왔으므로, 어떤 학급 간에 차이가 있는지를 알아보고 싶다면, 사후 쌍을 이룬 다중 비교 분석(post-hoc pair-wise multiple comparison)을 하면 됩니다.  이에 대한 자세한 설명은 아래의 링크를 참조하시기 바랍니다.

  • 쌍을 이룬 집단 간 평균 다중비교 (multiple comparison)

Tukey's HSD(honestly significant difference) test 참조

Duncan's LSR(least significant range) test 참고

 

 

  • 대비 (contrast)

샤페 검정법 (scheffe test) 참고

 

 

많은 도움이 되었기를 바랍니다. 

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

728x90
반응형
Posted by Rfriend
,

2개의 모집단에 대한 평균을 비교, 분석하는 통계적 기법으로 t-Test를 활용하였다면, 비교하고자 하는 집단이 3개 이상일 경우에는 분산분석 (ANOVA : Analysis Of Variance)를 이용합니다. 

 

설명변수는 범주형 자료(categorical data)이어야 하며, 종속변수는 연속형 자료(continuous data) 일 때 3개 이상 집단 간 평균 비교분석에 분산분석(ANOVA) 을 사용하게 됩니다.

 

분산분석(ANOVA)은 기본적으로 분산의 개념을 이용하여 분석하는 방법으로서, 분산을 계산할 때처럼 편차의 각각의 제곱합을 해당 자유도로 나누어서 얻게 되는 값을 이용하여 수준평균들간의 차이가 존재하는 지를 판단하게 됩니다.  

 

이전 포스팅에서 '일원분산분석(one-way ANOVA)'에 대해서 알아봤는데요, 이번 포스팅에서는 '이원분산분석(two-way ANOVA)'에 대해서 소개하도록 하겠습니다. 

 

 

[ 분산분석(ANOVA)의 분류 ]

 

 

 

일원분산분석(one-way ANOVA)이 1개의 요인(factor) 내의 요인 수준(factor levels)들이 각각의 집단(group)/처리(treatment)이 되어서 이들 집단/처리 간의 평균 차이를 비교하는 것이라면, 이원분산분석(two-way ANOVA)은 2개의 요인(2 factors) 내의 요인 수준(factor levels) 간의 조합(combination)들 각 각을 개별 집단/처리(groups, treamments)로 간주하고 이들간에 평균을 비교하게 됩니다.

 

가령, 요인(factor) A '온도'가 3개의 요인 수준(factor levels, 온도 상, 중, 하)을 가지고 요인(factor) B '압력'이 2개의 요인 수준(factor levels, 압력 강, 약)을 가진다고 할 경우, 총 그룹/처리의 수는 (A.온도) 3 x (B.압력) 2 = 6 개가 됩니다.

 

이원분산분석은 (1) 관측값이 하나일 경우와 (2) 관측값이 2개 이상일 경우 (반복 실험을 할 경우)로 나누어볼 수 있습니다.  비용이나 시간 여건이 허락한다면 분석의 신뢰도를 높이기 위해서는 반복 실험 혹은 관찰을 통해 관측값을 2개 이상 확보하는 것이 좋겠습니다.

 

우선 이번 포스팅에서는 (1) 관측값이 하나일 경우의 이원분산분석(two-way ANOVA when there is one observation in each cell (different treatment groups)) 에 대해서 소개하고, 다음번 포스팅에서 (2) 관측값이 2개 이상일 경우(반복 실험을 할 경우)의 이원분산분석에 대해서 순차적으로 소개하도록 하겠습니다.

 

이원분산분석을 위한 데이터셋의 구조는 아래와 같습니다.

 

 

[ 이원분산분석을 위한 데이터셋 구조 (Dataset for two-way ANOVA) ]

 

 

 

요인 A가 i=1, 2, ..., a 개의 요인 수준(factor levels)을 가지고, 요인 B가 j=1, 2, ..., b 개의 요인 수준(factor levels)을 가진다고 했을 때, A와 B라는 2개의 요인 처리(treatment) 내의 관측값이 하나일 경우의 이원분산분석모형은 편차 (Yij - Y..bar)를 다음과 같이 3개의 성분합으로 나타낼 수 있습니다. 


 

이 식의 양변을 제곱하여 모든 관측값에 대해 더하면 다음과 같은 식을 얻게 됩니다.

 

 

 

관측값이 하나일 경우의 이원분산분석을 실시할 경우 통계패키지에서는 아래와 같은 형태로 정리된 이원분산분석표를 제시하여 줍니다.

 

 

 [ 이원분산분석표 (two-way ANOVA table) ]

 

 요인

제곱합

(squared sum)

자유도

(degrees of freedom) 

평균제곱

(mean squared) 

F statistics 

 요인 A

 SSA

a-1

MSA

MSA/MSE

 요인 B

 SSB

b-1

MSB

MSB/MSE

 오차

 SSE

(a-1)(b-1) 

MSE 

 

 계

 SST

 ab-1

 

 

 

 

검정통계량으로는 F 통계량(요인A 효과 검정 = MSA/MSE, 요인B 효과 검정 = MSB/MSE)을 사용하며, 요인 A의 효과와 요인 B의 효과에 대한 검정 방법은 아래와 같습니다.

 

 

(1) 요인 A의 효과

   - 귀무가설 H0 : α1 = ... = αa = 0

   - 대립가설 H1 : 모든 αi 는 0이 아니다

   - 검정통계량 : F0 = MSA/MSE

   - 판정 : P-value가 유의수준(α)보다 작으면 귀무가설 H0를 기각하고, 대립가설 H1을 채택

 

 

(2) 요인 B의 효과

   - 귀무가설 H0 : α1 = ... = αb = 0

   - 대립가설 H1 : 모든 αj 는 0이 아니다

   - 검정통계량 : F0 = MSB/MSE

   - 판정 : P-value가 유의수준(α)보다 작으면 귀무가설 H0를 기각하고, 대립가설 H1을 채택

 

 

 

그럼, 아래 예제에 대해서 R의 aov() 함수를 활용해서 문제를 풀어보도록 하겠습니다.

(예제는 'Excel을 이용한 실용통계분석', 배현웅 저, 교우사, 에서 인용함)

 

 

(문제) K와 M 두 보험회사의 차종 (1,000cc 이하, 1,500cc, 1,800cc)에 따른 분기별 보험료(단위: 천원)가 아래 표와 같다고 할 때 유의수준 α=0.05 로 차종(요인 A)과 회사(요인 B)의 효과에 대한 검정을 하여라.

 

[ 보험회사와 차종에 따른 보험료 ]

 

              보험회사

 차종

K 회사 

M 회사

평균 

1,000 cc 이하 

140

100

120

 1,500 cc

 210

180

 195

 1,800 cc

 220

200

 210

 평균

 190

160

 175

 

 

 

> ##--------------------------------------------------------------
> ## two-way ANOVA : aov()
> ##--------------------------------------------------------------
> 
> # (1) two-way ANOVA when there is one observation in each cell (different treatment groups)
> 
> car_type <- rep(c('1000', '1500', '1800'), 2)
> car_type <- as.factor(car_type) # transformation into factor
> car_type
[1] 1000 1500 1800 1000 1500 1800
Levels: 1000 1500 1800
> 
> insurance <- as.factor(c(rep('K', 3), rep('M', 3))) # transforamtion into factor
> insurance
[1] K K K M M M
Levels: K M
> 
> y <- c(140, 210, 220, 100, 180, 200)
> 
> 
> # two way ANOVA
> two_way_aov_model_1 <- aov(y ~ car_type + insurance) # no replicates, no interaction
> 
> # statistics of two-way ANOVA
> summary(two_way_aov_model_1)
            Df Sum Sq Mean Sq F value Pr(>F)  
car_type     2   9300    4650      93 0.0106 *
insurance    1   1350    1350      27 0.0351 *
Residuals    2    100      50                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

 

위의 분석결과를 해석해보면, 먼저 요인A 차종에 따른 분산분석 결과 P-value가 0.0106으로서 유의수준(significance level) 0.05보다 작으므로 우리는 "차종에 따라서 보험료에 차이가 있다"는 대립가설(H1)을 채택할 수 있게 되었습니다.

 

또한, 요인B 보험회사에 따른 분산분석 결과 P-value가 0.0351로서 유의수준(significance level, α) 0.05보다 역시 작으므로 "보험회사에 따라서 보험료에 차이가 있다"는 대립가설(H1)을 채택할 수 있겠습니다.

 

다음번 포스팅에서는 (2) 관측값이 2개 이상일 경우(반복 실험일 경우)의 이원분산분석에 대해서 알아보겠습니다.

 

많은 도움 되었기를 바랍니다.

 

 

일원분산분석 및 사후분석(post-hoc multiple comparison)에 대해서는 아래 링크를 참조하세요.

  • 1개 요인(factor)에 대한 3 집단 이상 집단의 평균 비교 (ANOVA)

one-way ANOVA

 

 

  • 쌍을 이룬 집단 간 평균 다중비교 (multiple comparison)

Tukey's HSD(honestly significant difference) test 참조

Duncan's LSR(least significant range) test 참고

 

  • 대비 (contrast)

샤페 검정법 (scheffe test) 참고

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

 

728x90
반응형
Posted by Rfriend
,

지난 포스팅에서 (1) 한개의 요인에 대해 집단 간 평균을 비교하는 one-way ANOVA, (2) 짝을 이룬 두 집단간 사후 다중비교(post-hoc pair-wise multiple comparison) 방법으로 튜키 검정법(Tukey HSD(honestly significant difference) Test), 던칸 검정법(Duncan LSR(least significant range) Test) 에 대해서 알아보았습니다.  복습을 해보자면, 다중비교는 설정된 유의수준의 크기를 그대로 유지한 상태에서 모든 가능한 두 수준평균간의 차이에 대한 검정을 사후에 동시에 시행하는 방법이었습니다.

 

☞ one-way ANOVA 참고

☞ Tukey's HSD(honestly significant difference) test 참고

☞ Duncan's LSR(least significant range) test 참고

 

 

이번 포스팅에서는 (a) 다중비교 시에 표본의 크기가 다를 경우에도 우수한 검정력을 보여주며 (b) 세 개 이상의 요인 수준(factor levels) 간의 평균들을 비교할 때도 사용 가능한 샤페검정법(scheffe test)을 소개하도록 하겠습니다. 

 

 

"Scheffé's method is a single-step multiple comparison procedure which applies to the set of estimates of all possible contrasts among the factor level means, not just the pairwise differences considered by the Tukey–Kramer method"

 

* source : wikipedia 

 

 

Tukey HSD test와 Duncan LSR test 와 Scheffe test 간의 차이점 비교는 아래 도식화를 참고하시기 바라며, Scheffe test는 Duncan LSR test와 함께 상당히 보수적(conservative)인 검정법으로 알려져 있습니다. 

 

 

[categorization of test method by '# of factor levels' & 'sample size']

 

 

 

 

분산분석에서는 둘 이상의 수준(factor levels) 평균들을 포함하는 비교를 하기 위해 각 수준들의 평균 μi 를 하나의 직선식으로 표현한 것을 대비(Contrasts)라고 하며 L로 표기하면, L의 추정량(estimated contrast L^)과 추정량의 분산(estimated variance of L^), 샤폐계수 T (scheffe's T statistics)는 아래와 같습니다.

 

 

[ Scheffe's T statistics ]

 

 

 

 

Scheffe Test 의 가설과 검정 방법은 아래와 같습니다.

 

 

 

 


 

 

R에서 scheffe test를 하기 위해 agricolae package scheffe.test() 함수를 사용하였습니다. 

 

그리고 데이터는 agricolae package에 내장되어있는 sweetpotato dataframe을 사용하였습니다.  virus 라는 1개의 요인(one-way)에 대해서 cc, fc, ff, oo 의 4개의 수준(4 factor levels)에 대해 감자 생산량(yield)의 차이가 있는지를 동시에 비교 평가하기 위해 scheffe test를 적용해보았습니다.  (요인수준 4개에 대해 동시에 비교하기 위해 scheffe test 적용함.  유의수준 α = 0.05)

 

one-way ANOVA 분석 들어가기 전에 탐색적 분석 차원에서 boxplot을 그려보고, 기초통계량(doBy package의 summaryBy 함수 사용)을 뽑아보았습니다.  참고하시구요.

 

 

> ##--------------------------------------------------------------
> ## single-step multiple comparison, contrasts : scheffe.test()
> ##--------------------------------------------------------------
> 
> install.packages("agricolae")
Installing package into ‘C:/Users/user/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.2/agricolae_1.2-3.zip'
Content type 'application/zip' length 920587 bytes (899 KB)
downloaded 899 KB

package ‘agricolae’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpyASqz8\downloaded_packages
> library(agricolae)
> data(sweetpotato)
> 
> str(sweetpotato)
'data.frame':	12 obs. of  2 variables:
 $ virus: Factor w/ 4 levels "cc","fc","ff",..: 1 1 1 2 2 2 3 3 3 4 ...
 $ yield: num  28.5 21.7 23 14.9 10.6 13.1 41.8 39.2 28 38.2 ...
> 
> # boxplot of yield of sweetpatoto, dealt with different virus
> attach(sweetpotato)
> boxplot(yield ~ virus, 
+         main = "Yield of sweetpotato, Dealt with different virus", 
+         xlab = "Virus", 
+         ylab = "Yield")
> detach(sweetpotato)
>

 

 

 
> 
> # statistics summary by group("virus")
> install.packages("doBy")
Installing package into ‘C:/Users/user/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.2/doBy_4.5-13.zip'
Content type 'application/zip' length 3431380 bytes (3.3 MB)
downloaded 3.3 MB

package ‘doBy’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpyASqz8\downloaded_packages
> library(doBy)
필요한 패키지를 로딩중입니다: survival
> summaryBy(yield ~ virus, data=sweetpotato, FUN = c(mean, sd, min, max))
  virus yield.mean yield.sd yield.min yield.max
1    cc   24.40000 3.609709      21.7      28.5
2    fc   12.86667 2.159475      10.6      14.9
3    ff   36.33333 7.333030      28.0      41.8
4    oo   36.90000 4.300000      32.1      40.4
> 
> 
> # ANOVA test
> aov_sweetpotato <- aov(yield~virus, data=sweetpotato)
> summary(aov_sweetpotato)
            Df Sum Sq Mean Sq F value   Pr(>F)    
virus        3 1170.2   390.1   17.34 0.000733 ***
Residuals    8  179.9    22.5                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> 
> # scheffe test
> comparison <- scheffe.test(aov_sweetpotato, # ANOVA model 
+                            "virus", # vector treatment applied to each experimental unit
+                            alpha = 0.05, # significant level
+                            group=TRUE, 
+                            console=TRUE, # print out
+                            main="Yield of sweetpotato\nDealt with different virus")

Study: Yield of sweetpotato
Dealt with different virus

Scheffe Test for yield 

Mean Square Error  : 22.48917 

virus,  means

      yield      std r  Min  Max
cc 24.40000 3.609709 3 21.7 28.5
fc 12.86667 2.159475 3 10.6 14.9
ff 36.33333 7.333030 3 28.0 41.8
oo 36.90000 4.300000 3 32.1 40.4

alpha: 0.05 ; Df Error: 8 
Critical Value of F: 4.066181 

Minimum Significant Difference: 13.52368 

Means with the same letter are not significantly different.

Groups, Treatments and means
a 	 oo 	 36.9 
a 	 ff 	 36.33 
ab 	 cc 	 24.4 
b 	 fc 	 12.87

 

 

위의 scheffe test 수행 결과를 보니 virus (treatement) 'oo'와 'ff'는 같은 'a' group으로 묶였으며 0.05 유의수준 하에서 차이가 없는 것(귀무가설 H0 채택)으로 나왔습니다. virus (treatement) 'fc'는 'b' group으로 따로 묶였으며, 유의수준 0.05 에서 'oo', 'ff'와는 차이가 있는 것(귀무가설 H0 기각)으로 판단할 수 있습니다. virus (treatment) 'cc'는 'a'와 'b' group 사이에 애매하게 끼어서 차이가 없는 것으로 나왔네요.  위에 파란색으로 해놓은 'Minimum Significant Difference: 13.52368' 를 보면 'cc' virus 에 의한 감자 생산량 평균 '24.4'가 'fc' virus에 의한 생산량 '12.86'과의 차이가 '13.52368'보다는 작으며, 'ff' virus에 의한 생산량 '36.33'과의 차이도 역시 '13.5268'보다는 작음을 알 수 있습니다.

 

많은 도움 되셨기를 바랍니다.

 

☞ one-way ANOVA 참고

☞ Tukey's HSD(honestly significant difference) test 참고

☞ Duncan's LSR(least significant range) test 참고

 

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서는 ANOVA 검정 시 집단 간 평균에 유의미한 차이가 있다고 판단될 때 사후 분석으로 다중비교((post-hoc pair-wise multiple comparison)을 하는 방법 중에 Tukey's HSD(honestly significant difference) test에 대해서 알아보았습니다.

 

이번 포스팅에서는 다중비교(multiple comparison)을 하는 두번째 방법으로 Duncan's LSR(Least Significant Range) test에 대해서 알아보도록 하겠습니다.

 

 

 

[ 사후검정 중비교(post-hoc pair-wise multiple comparison) ]

 

 

 

 

 

 

기본 개념은 짝을 이룬(pair-wise) 수준평균간의 차이가 아래에 제시된 Duncan의 최소유의범위(Least Significant Range) 보다 크다면 "유의수준 α에서 두 수준평균간에는 차이가 있다"고 판단하겠다는 것입니다.

 

 

 

[ Duncan's LSR(least significant range) test procedure ]

 

 

Duncan's LSR test는 앞서 포스팅했던 Tukey's HSD test보다도 약간 더 보수적입니다.  그래서 Tukey's HSD test에서는 차이가 없다고 나오는 짝을이룬 수준평균간에도 Duncan's LSR test에서는 차이가 있다고 판정하는 경우도 있습니다.

 

R을 가지고 Duncan's LSR test를 하는데 있어 

 - duncan.test() function of agricolase package

 - LDuncan() function of laercio package

두가지 함수를 차례로 사용해보겠습니다.

 

온도의 조건 3가지에 따른 생산량을  각10번씩 측정한 아래의 데이터를 가지고 (1) ANOVA test를 먼저 해보고, (2) 사후 분석으로 Duncan's LSR test를 위의 2가지 패키지의 함수를 가지고 각각 실시해보겠습니다.

 

 

> ##---------------------------------------------------------- > ## multiple comparison > ## - Duncan's LSR(least significant range) test : > ##---------------------------------------------------------- > > ##--- Are there any daily outcome differences among temperature conditions? > # group 1 : temperature condition 1 > # group 2 : temperature condition 2 > # group 3 : temperature condition 3 > > # daily outcome by tmep condition (group 1/2/3) > y1 <- c(50.5, 52.1, 51.9, 52.4, 50.6, 51.4, 51.2, 52.2, 51.5, 50.8) > y2 <- c(47.5, 47.7, 46.6, 47.1, 47.2, 47.8, 45.2, 47.4, 45.0, 47.9) > y3 <- c(46.0, 47.1, 45.6, 47.1, 47.2, 46.4, 45.9, 47.1, 44.9, 46.2) > > y <- c(y1, y2, y3) > y [1] 50.5 52.1 51.9 52.4 50.6 51.4 51.2 52.2 51.5 50.8 47.5 47.7 46.6 47.1 47.2 47.8 45.2 47.4 45.0 [20] 47.9 46.0 47.1 45.6 47.1 47.2 46.4 45.9 47.1 44.9 46.2 > > n <- rep(10, 3) > n [1] 10 10 10 > > group <- rep(1:3, n) > group [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 > > # combining into data.frame > group_df <- data.frame(y, group) > group_df y group 1 50.5 1 2 52.1 1 3 51.9 1 4 52.4 1 5 50.6 1 6 51.4 1 7 51.2 1 8 52.2 1 9 51.5 1 10 50.8 1 11 47.5 2 12 47.7 2 13 46.6 2 14 47.1 2 15 47.2 2 16 47.8 2 17 45.2 2 18 47.4 2 19 45.0 2 20 47.9 2 21 46.0 3 22 47.1 3 23 45.6 3 24 47.1 3 25 47.2 3 26 46.4 3 27 45.9 3 28 47.1 3 29 44.9 3 30 46.2 3 > > sapply(group_df, class) y group "numeric" "integer" > > # transform from 'integer' to 'factor' > group_df <- transform(group_df, group = factor(group)) > sapply(group_df, class) y group "numeric" "factor" > > > # (1) ANOVA test > aov_model <- aov(y ~ group, data = group_df) > summary(aov_model) Df Sum Sq Mean Sq F value Pr(>F) group 2 156.30 78.15 108.8 1.2e-13 *** Residuals 27 19.39 0.72 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > > # (2-1) Duncan's LSR(least significant range) test

> # with duncan.test() function of agricolase package > install.packages("agricolae")

Installing package into ‘C:/Users/user/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.2/agricolae_1.2-3.zip'
Content type 'application/zip' length 921635 bytes (900 KB)
downloaded 900 KB

package ‘agricolae’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpYhI0Fo\downloaded_packages

 

> library(agricolae) > duncan.test(aov_model, "group", alpha = 0.05, console = TRUE) Study: aov_model ~ "group" Duncan's new multiple range test for y Mean Square Error: 0.7182593 group, means y std r Min Max 1 51.46 0.6834553 10 50.5 52.4 2 46.94 1.0415800 10 45.0 47.9 3 46.35 0.7763876 10 44.9 47.2 alpha: 0.05 ; Df Error: 27 Critical Range 2 3 0.7776731 0.8170522 Means with the same letter are not significantly different. Groups, Treatments and means a 1 51.46 b 2 46.94 b 3 46.35 > >

>

> # (2-2) Duncan's LSR(least significant range) test

> # with LDuncan() function of laercio package > install.packages("laercio")

Installing package into ‘C:/Users/user/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.2/laercio_1.0-1.zip'
Content type 'application/zip' length 20726 bytes (20 KB)
downloaded 20 KB

package ‘laercio’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpYhI0Fo\downloaded_packages

 

> library(laercio) > aov_model <- aov(y ~ group, data = group_df) > LDuncan(aov_model, "group") DUNCAN TEST TO COMPARE MEANS Confidence Level: 0.95 Dependent Variable: y Variation Coefficient: 1.75648 % Independent Variable: group Factors Means 1 51.46 a 2 46.94 b 3 46.35 b Restarting R session...

 

결과를 해석해보면, (1) ANOVA 검정 결과 P-value가 1.2e-13 로서 유의수준 0.05보다 매우 작으므로 "집단간 평균에 차이가 있다"는 대립가설을 채택할 수 있습니다.

 

(2) 이에 사후분석으로 "어느 집단 간에 차이가 있는가?"를 알아보기 위해 Duncan's LSR test 를 해보니,

온도의 조건에 따른 3가지 factor level 별 생산량의 평균은 아래와 같았고, factor level '2'와 '3'은 서로 차이가 없고, factor level '1'만 짝을 이룬 수준평균간에 차이가 있다고 나왔습니다.

 

 factor level

mean 

 group

significantly different?

 1

 51.46

 a ***

 different

 2

 46.94

 b

 not different

 3

 46.35

 b

 

 

차이가 있다고 판단한 근거를 보면, MSE(mean squared error) 가 0.7182593 이고, 이 값을 위에서 제시했던 Duncan's LSR(Least Significant Range)를 계산하는 공식에 넣어서 계산을 해보면, 위에 R분석결과에 Critical Range 라고 해서 LSR2 = 0.7776731, LSR3 = 0.8170522 가 나왔습니다.  이 값들보다 수준집단간 평균의 차이가 더 크면 차이가 있다고 판단하게 되는 것입니다. 

 

3개 수준평균을을 크기 순서대로 배열한 후에 각 차이를 구한 것을 표로 제시해보면 아래와 같으며, 이 값들이 LSR2 = 0.7776731, LSR3 = 0.8170522 보다 크면 *** 표시를 했습니다.

(LSR2 는 인접한 수준들간의 평균 차이를 의미, LSR3는 2단계 떨어진 수준들간의 평균 차이를 의미)

 

 

mean of factor level 3: 

46.35

mean of factor level 2:

46.94

mean of factor level 1:

51.46 

mean of factor level 3:

 46.35

 0.59

5.11 ***

mean of factor level 2: 

46.94

 

4.52 ***

mean of factor level 1:

 51.46

 

 

 

많은 도움 되었기를 바랍니다.

 

☞ one-way ANOVA 참고

☞ 다중비교(multiple comparison) Tukey's HSD(honestly significant difference) test 참고

☞ 대비(contrast) Scheffe test 참고

 

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

 

728x90
반응형
Posted by Rfriend
,

지난번 포스팅에서 일원분산분석(one-way ANOVA) 이론과 R 의 aov() 함수에 대해서 알아보았습니다.  그러면 이번 포스팅에서는 일원분산분석 후에 이어서 하게 되는 다중비교(Multiple Comparison)에 대해서 알아보겠습니다.

 

F-검정 통계량을 사용해서 한개의 요인(Factor) 내 세 개 이상의 수준(Levels)의 집단 간에 평균의 차이가 있는지 없는지를 검정하였는데, 만약 분석 결과가 "인자 수준들의 평균 μi 는 모두 같다"는 귀무가설 H0를 기각하고 "인자 수준들의 평균 μi 는 모두 같지 않다" 대립가설 H1을 채택할 경우가 생겼다고 합시다. 이럴 경우 수준(Levels)의 집단 간에 어디에서 차이가 생겼는지 알고 싶을 때 사용하는 분석법이 다중비교(multiple comparison)이 되겠습니다.  

 

수준의 수가 2일 때 일원분산분석에서 귀무가설 H0가 기각될 경우에는 두 수준 평균간의 차이에 대해서는 t-Test를 적용하면 됩니다.  하지만 수준의 수(number of Levles, groups)가 3개 이상인 경우에는 t-Test를 사용하면 안됩니다.  왜냐하면 비교해야할 수준수가 커질수록 두 수준간의 짝들 간의 비교(pair-wise comaprison)를 하게 되면 제 1종 오류(참을 거짓으로 판정할 오류)가 애초 정한 유의수준 α(significance level α)가 아니라 이보다 훨씬 커지게 되기 때문입니다.

 

t-Test를 사용할 경우 무슨 일이 벌어지는지 예를 들어서 설명해보겠습니다.  만약 수준의 수(number of Levels, groups)가 5개일 경우 각 각 쌍을 지어서 비교를 한다고 했을 때 수준평균짝들의 총 갯수는 10개(5 Combination 2 = 4 + 3 + 2 + 1 = 10) (R의 combination 함수를 사용하면, choose(5,2) 하면 10이 나옴) 이므로 t-Test를 10번 적용해야 합니다.  만약 유의수준 α (significance level α) = 0.05라고 한다면, 10번의 t-Test를 각각 독립적으로 하고 나면 제1종 오류가 생길 확률은 1 - (0.95)^10 = 0.4012631 가 됩니다.  이는 처음의 유의수준 α 0.05보다 약 8배 크며, 1종 오류가 생길 확률도 약 8배 정도 크다는 소리입니다. 

 

이러한 오류를 피하기 위해서는 유의수준을 설정된 크기로 유지하면서 모든 두 수준짝들의 평균을 동시에 비교할 수 있도록 고안된 검정방법인 다중비교(multiple comparison)으로 Tukey's HSD(honestly significant difference) testDuncan's LSR(least significant range) test 가 있는데요, 이번 포스팅에서는 Tukey's HSD test 에 대해서 소개하겠습니다.

 

Tukey's HSD(honestly significant difference) test는 studentized range distribution을 이용하여 모든 가능한 두 수준들의 평균간의 차이가 있는지를 검정(pairwise post-hoc testing using Tukey HSD test)하는 방법입니다.  Tukey's HSD test 의 통계량 및 검정 방법은 아래와 같습니다.

 

 

 

[ Tukey's HSD(honestly significant difference) test ]

 

 

 

 

그러면, (1) R의 aov() 함수를 사용해서 유의수준 0.05 하에 일원분산분석(one-way ANOVA)를 먼저 해보고 (<= 이것은 지난번 포스팅과 동일합니다), (2) Tukey's HSD test를 통해서 다중비교(multiple comparison)을 하여 쌍을 이룬 집단간의 평균 비교 중에서 어디서 차이가 발생한건지 알아보도록 하겠습니다.

 

 

> ##----------------------------------------------------------
> ## One-way ANOVA : aov(), oneway.test
> ##----------------------------------------------------------
> 
> ##--- Are there any daily outcome differences among temperature conditions?
> # group 1 : temperature condition 1 
> # group 2 : temperature condition 2
> # group 3 : temperature condition 3
> 
> # daily outcome by tmep condition (group 1/2/3)
> y1 <- c(50.5, 52.1, 51.9, 52.4, 50.6, 51.4, 51.2, 52.2, 51.5, 50.8)
> y2 <- c(47.5, 47.7, 46.6, 47.1, 47.2, 47.8, 45.2, 47.4, 45.0, 47.9)
> y3 <- c(46.0, 47.1, 45.6, 47.1, 47.2, 46.4, 45.9, 47.1, 44.9, 46.2)
> 
> y <- c(y1, y2, y3)
> y
 [1] 50.5 52.1 51.9 52.4 50.6 51.4 51.2 52.2 51.5 50.8 47.5 47.7 46.6 47.1 47.2 47.8 45.2 47.4 45.0
[20] 47.9 46.0 47.1 45.6 47.1 47.2 46.4 45.9 47.1 44.9 46.2
> 
> n <- rep(10, 3)
> n
[1] 10 10 10
> 
> group <- rep(1:3, n)
> group
 [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3
> 
> # combining into data.frame
> group_df <- data.frame(y, group)
> group_df
      y group
1  50.5     1
2  52.1     1
3  51.9     1
4  52.4     1
5  50.6     1
6  51.4     1
7  51.2     1
8  52.2     1
9  51.5     1
10 50.8     1
11 47.5     2
12 47.7     2
13 46.6     2
14 47.1     2
15 47.2     2
16 47.8     2
17 45.2     2
18 47.4     2
19 45.0     2
20 47.9     2
21 46.0     3
22 47.1     3
23 45.6     3
24 47.1     3
25 47.2     3
26 46.4     3
27 45.9     3
28 47.1     3
29 44.9     3
30 46.2     3
> 
> sapply(group_df, class)
        y     group 
"numeric" "integer" 
> 
> # transform from 'integer' to 'factor'
> group_df <- transform(group_df, group = factor(group))
> sapply(group_df, class)
        y     group 
"numeric"  "factor" 
> 
> 
> # boxplot
> attach(group_df)
The following objects are masked _by_ .GlobalEnv:

    group, y

> boxplot(y ~ group, 
+         main = "Boxplot of Daily Outcome by Temperature condition 1/2/3", 
+         xlab = "Factor Levels : Temperature condition 1/2/3", 
+         ylab = "Daily Outcome")
> 

 

 
> 
> # descriptive statistics by group
> tapply(y, group, summary)
$`1`
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  50.50   50.90   51.45   51.46   52.05   52.40 

$`2`
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  45.00   46.72   47.30   46.94   47.65   47.90 

$`3`
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  44.90   45.92   46.30   46.35   47.10   47.20 

> detach(group_df)
> 
> # one-wayANOVA
> aov(y ~ group, data = group_df)
Call:
   aov(formula = y ~ group, data = group_df)

Terms:
                  group Residuals
Sum of Squares  156.302    19.393
Deg. of Freedom       2        27

Residual standard error: 0.8475018
Estimated effects may be unbalanced
> summary(aov(y ~ group, data = group_df))
            Df Sum Sq Mean Sq F value  Pr(>F)    
group        2 156.30   78.15   108.8 1.2e-13 ***
Residuals   27  19.39    0.72                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> # Bartlett test to test the null hypothesis of equal group variances
> bartlett.test(y ~ group, data = group_df)

	Bartlett test of homogeneity of variances

data:  y by group
Bartlett's K-squared = 1.6565, df = 2, p-value = 0.4368

>  
> ##------------------------------------------------------------------
> ## multiple comparison
> ## - Tukey's HSD(honestly significant difference) test : TukeyHSD()
> ##------------------------------------------------------------------
> group_aov <- aov(y ~ group, data = group_df)
> TukeyHSD(group_aov)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = y ~ group, data = group_df)

$group
     diff       lwr        upr     p adj
2-1 -4.52 -5.459735 -3.5802652 0.0000000
3-1 -5.11 -6.049735 -4.1702652 0.0000000
3-2 -0.59 -1.529735  0.3497348 0.2813795

 

 

 

위의 일원분산분석(one-way ANOVA) 결과를 보면 P-value가 1.2e-13 으로서 대립가설 H1을 채택하게 되어 평균이 서로 다르다고 판단할 수 있으므로, 자연스럽게 그 다음 질문은 "다르다고 하는데 그러면 쌍을 이룬 평균 비교 중에서 어디에서 차이가 난다는 말인가?"가 될 것입니다.

 

이 질문에 답하기 위해 위의 Tukey's HSD test 결과를 보면, "multiple comparisons of means, 95% family-wise confidence level'이라는 설명이 나오는데요, 2개씩 쌍을 이룬 수준간 평균의 다중 비교를 95% 신뢰수준으로 상한(upr)과 하한(lwr)의 신뢰계수 구간을 구했고, P-value도 구해주었습니다.

 

위의 결과를 보면 'group 2'와 'group 1'은 평균 차이가 -4.52 이고 adj. P-value가 0.0000000이고, 'group 3'과 'group 1'의 평균 차이는 -5.11이고 adj. P-value가 0.0000000 으로서, 유의수준 0.05보다 훨씬 작으므로 이들 group (수준, factor levels) 간에는 평균의 차이가 있다고 유의수준 0.05 기준 하에 판단할 수 있겠습니다.  반면에, 'group 3'과 'group 2'는 평균 차이가 -0.59이고 adjusted P-value 가 0.2813795 로서 유의수준 0.05보다 크므로 귀무가설 H0 '두 집단 간 평균차이는 없다'를 채택하게 됩니다.  이는 저 위에 있는 Boxplot을 보면 좀더 쉽게 이해가 될 것 같습니다.

 

위 분석 결과를 보기에 좋게 표로 나타내면 아래와 같습니다.

 

 pair-wise comparison

 group 1

group 2

group 3 

 group 1

mean diff. -4.51 ***

(95% confidence interval

: -5.45 ~ -3.58)

mean diff. -5.11 ***

(95% confidence intervl

: -6.04 ~ -4.17)

 group 2

 

mean diff. -0.59 

(95% confidence interval

: -1.52 ~ 0.34)

 group 3

 

 

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

 

 

Tukey's HSD test의 이론에 대해서는 알쏭달쏭한 말들을 잔뜩 써놨는데요, R 함수는 TukeyHSD() 딱 한줄이어서 미친 듯이 간단합니다. ^^'  다음번 포스팅에서는 Duncan's LSR(least significant range) test에 대해서 알아보도록 하겠습니다.

 

많은 도움 되었기를 바랍니다.

 

☞ one-way ANOVA 참고

☞ 다중비교(multiple comparison) : Duncan's LSR(least significant range) test 참고

☞ 대비(contrast) : Scheffe test 참고

 

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

728x90
반응형
Posted by Rfriend
,

2개의 모집단에 대한 평균을 비교, 분석하는 통계적 기법으로 t-Test를 활용하였다면, 비교하고자 하는 집단이 3개 이상일 경우에는 분산분석 (ANOVA : Analysis Of Variance)를 이용합니다. 

 

설명변수는 범주형 자료(categorical data)이어야 하며, 종속변수는 연속형 자료(continuous data) 일 때 3개 이상 집단 간 평균 비교분석에 분산분석(ANOVA) 을 사용하게 됩니다.

 

분산분석(ANOVA)은 기본적으로 분산의 개념을 이용하여 분석하는 방법으로서, 분산을 계산할 때처럼 편차의 각각의 제곱합을 해당 자유도로 나누어서 얻게 되는 값을 이용하여 수준평균들간의 차이가 존재하는 지를 판단하게 됩니다.  

 

 

 

[ 일원분산분석 (one-way ANOVA)의 개념 ]

 

 

 

 

분산분석 (ANOVA)은 영국의 통계학자 피셔(R.A. Fisher)가 고안한 분석 기법으로서, 최초에는 농사실험과 관련된 연구에서 사용되었으나 요즘에는 사회과학, 공학, 의학 등 폭넓게 적용되고 있습니다.

 

 

[ 분산분석 (ANOVA)의 분류 ]

 

 

 

3개 이상의 집단에 대한 비교 시 모집단이 정규분포를 따르지 않을 경우에는 비모수 추정과 검정 기법을 이용해야 하며, 이에 대해서는 추후에 별도로 다루도록 하겠습니다.

 

 

분산분석은 측정하고자 하는 값에 영향을 미치는 요인(Factor)의 수에 따라서 구분하는데요, 가령, 작업자(Worker)에 따라서 생산량(Volume of Production)의 차이가 있는가를 비교 분석한다고 했을 때 '작업자(Worker)'가 바로 요인(Factor)이 되겠습니다. 그리고 작업자 3명 '홍길동', '김철희', '최영희'를 요인 수준 (Factor Level) 혹은 처리 (treatment) 라고 합니다.

 

측정값에 영향을 미치는 요인(Factor)이 1 개인 실험은 '일원분산분석' (one-way ANOVA) 라고 하며, 측정값에 영향을 미치는 요인(Factor)이 2 개인 실험은 '이원분산분석' (two-way ANOVA) 라고 부릅니다.

 

 

[ ANOVA Model and Description in R ]

 

n-way ANOVA

Model

Description

one-way ANOVA

y ~ x1

  y is explained by x1 only

two-way ANOVA

y ~ x1 + x2

  y is explained by x1 and x2

two-way ANOVA

y ~ x1 * x2

  y is explained by x1, x2 and the

interaction   
  between them

three-way ANOVA

y ~ x1 + x2 + x3

  y is explained by x1, x2 and x3

 

 

 

이번 포스팅에서는 요인(Factor)이 1개인 일원분산분석(one-way ANOVA)의 이론과 R의 aov() 함수의 사용법에 대해서 설명해보겠습니다.

 

요인(Factor)이 1개이면서 r개의 요인 수준(Factor Level)을 가지고 있고, 각 수준에서 박복 측정수가 n개인 일원분산분석(one-way ANOVA)는 다음과 같은 형태의 자료 형태와 모형을 가지게 됩니다.

 

 

 

[ 일원분산분석 데이터 형태 (dataset for one-way ANOVA) ]

 

 

 

 

[ 일원분산분석 모형 (one-way ANOVA model) ] 

 

one-way ANOVA model

 

분산분석(ANOVA)은 측정된 자료값들의 전체 변동을 비교하고자 하는 요인 수준(Factor Level) 간의 차이에 의해서 발생하는 변동과 그 밖의 요인에 의해서 발생하는 변동으로 나누어서 자료를 분석하는 것이 기본 원리가 되겠습니다.

 

 

[ 분산분석 기본 원리 1 ]

 

 

 

위에서 제시한 식의 양변을 제곱하여 전체 측정값들에 대해서 모두 더하여 식을 정리하면 아래와 같이 됩니다. (자세한 과정은 생략함)

 

 

[ 분산분석 기본 원리 2 ]

 

 

 

 

위의 총제곱합(SST)와 처리제곱합(SSTR), 오차제곱합(SSE) 간의 기본원리를 이용하여 분산분석에서 사용하는 통계량들을 표로 일목요연하게 정리한 것이 아래의 일원분산분석표(one-way ANOVA table) 가 되겠습니다.

 

 

[ 원분산분석표 (one-way ANOVA table) ]

 

 구분

 제곱합

(Sum of Squares)

자유도

(degrees of freedom) 

평균제곱

(Mean Square Error) 

검정통계량 F0

(F statistics)

 처리

(Treatment)

 SSTR

r-1 

MSTR 

MSTR/MSE 

 오차

(Error)

 SSE

 nT - r

 MSE

 

 전체

(Total)

 SST

 nT - 1

 

 

 

 

분산분석표의 제일 오른쪽에 있는 F0 통계량을 이용해서 전체 수준들간의 평균이 같은지 아니면 다른지를 검정합니다.  기본 개념을 설명하자면, F0 통계량은 처리평균제곱 (MSTR)의 크기에 영향을 받으므로 처리평균제곱 (MSTR)이 커지면 오차평균제곱(MSE)은 작아지게 되며, 따라서 F0 통계량은 분자가 커지고 분모가 작아지므로 당연히 커지게 됩니다. 즉, F0 통계량 값이 크다는 것은 수준 평균들간에 평균의 차이가 크다는 의미가 되겠습니다.

 

 

그러면, 요인효과에 대한 검정을 위해서 분산분석에서는 아래와 같은 귀무가설과 대립가설을 사용하며, 검정통계량으로는 F 를 사용하여 기각역 또는 P-value 에 의해서 검정을 하게 됩니다.  

 

 

  [ 가설 ]

   귀무가설 H0 : μ1 = μ2 = ... = μr

   대립가설 H1 : 모든 μi 는 같지 않다.  i = 1, 2, ..., r

 

  [ 기각역 혹은 P-value에 의한 검정 ]

   표본으로부터 계산된 검정통계량의 값 f0가 유의수준(significance level) α 에서

 

   f0 >= Fα(r - 1, nT - r) 또는 P-value = P(F >= f0) < α 이면, H0 기각 (H0 reject)

 

   f0 < Fα(r - 1, nT - r) 또는 P-value = P(F >= f0) >= α 이면, H0 채택 (H0 accept)

 

 

 

 

이론 설명이 무척 길어졌는데요, 이제 드디어, 아래 문제를 R의 aov() 함수를 사용해서 풀어보도록 하겠습니다.

 

 

 


 

 

문제 ) 정유사에서 온도(Factor, one-way)에 따라서 휘발유 생산량에 변화가 있는지 (즉, 영향이 있는지) 알아보기 위하여 온도를 3가지 조건(3 Factor Levels)으로 실험설계를 하여 10번에 걸쳐서 휘발유 생산량을 측정하였습니다. 관찰값이 아래와 같을 때 조사되었을 때 온도의 조건에 따라서 휘발유 생산량에 차이가 있는지 유의수준 α = 10% 로 검정하시오.

 

 

 

 10번 실험에 의한 측정값

 요인 수준

(Factor Level)

 1

2

 10

 온도 조건 1

(Temp condition 1)

 50.5

52.1

51.9

52.4 

50.6

51.4 

51.2 

52.2 

51.5 

50.8 

 온도 조건 2

(Temp condition 2)

 47.5

47.7 

46.6 

47.1 

47.2 

47.8 

45.2 

47.4 

45.0 

47.9 

 온도 조건 3

(Temp condition 3)

 46.0

47.1

45.6

47.1

47.2

46.4

45.9 

47.1

44.9

46.2 

 


 

R에 (1) 위의 관측값들을 벡터로 입력하고,

      (2) boxplot() 함수와 summary()를 이용하여 탐색적 분석을 해본 후에,

      (3) aov() 함수를 이용하여 one-way ANOVA 분석을 하였으며 
         : aov(y ~ group, data = dataset)

      (4) 기본 가정 중에 오차의 등분산성을 검정하기 위해 Bartlett 검정

         : bartlett.test(y ~ group, data = dataset)

을 실시하였습니다.

 

 

이때 조심해야 할 것이 있는데요, dataset을 데이터프레임으로 해서 그룹(group, factor level)에 해당하는 변수는 반드시 요인(Factor)형이어야 하므로, 만약 요인(Factor)형이 아니라면 먼저 transfrom(factor()) 함수를 이용해서 변환을 시켜주고 ANOVA 분석을 수행해야 합니다.

 

 

> ##---------------------------------------------------------- 
> ## One-way ANOVA : aov(), oneway.test 
> ##---------------------------------------------------------- 
>  
> ##--- Are there any daily outcome differences among temperature conditions? 
> # group 1 : temperature condition 1  
> # group 2 : temperature condition 2 
> # group 3 : temperature condition 3 
>  
> # daily outcome by tmep condition (group 1/2/3) 
> y1 <- c(50.5, 52.1, 51.9, 52.4, 50.6, 51.4, 51.2, 52.2, 51.5, 50.8) 
> y2 <- c(47.5, 47.7, 46.6, 47.1, 47.2, 47.8, 45.2, 47.4, 45.0, 47.9) 
> y3 <- c(46.0, 47.1, 45.6, 47.1, 47.2, 46.4, 45.9, 47.1, 44.9, 46.2) 
>  
> y <- c(y1, y2, y3) 
> y  
[1] 50.5 52.1 51.9 52.4 50.6 51.4 51.2 52.2 51.5 50.8 47.5 47.7 46.6 47.1 47.2 47.8 45.2 47.4 45.0 
[20] 47.9 46.0 47.1 45.6 47.1 47.2 46.4 45.9 47.1 44.9 46.2 
>  
> n <- rep(10, 3) > n [1] 10 10 10 
>  
> group <- rep(1:3, n) 
> group  [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 
>  
> # combining into data.frame 
> group_df <- data.frame(y, group) 
> group_df       
y group 
1  1 50.5     1 2  52.1     1 3  51.9     1 4  52.4     1 5  50.6     1 6  51.4     1 7  51.2     1 8  52.2     1 9  51.5     1 10 50.8     1 11 47.5     
2 12 47.7     2 13 46.6     2 14 47.1     2 15 47.2     2 16 47.8     2 17 45.2     2 18 47.4     2 19 45.0     2 20 47.9     2 21 46.0     
3 22 47.1     3 23 45.6     3 24 47.1     3 25 47.2     3 26 46.4     3 27 45.9     3 28 47.1     3 29 44.9     3 30 46.2     3 
>  
> sapply(group_df, class)         
y     group  
"numeric" "integer"  
>  
> # transform from 'integer' to 'factor' 
> group_df <- transform(group_df, group = factor(group)) 
> sapply(group_df, class)         
y     group  
"numeric"  "factor"  
>  
>  
> # boxplot 
> attach(group_df) 
The following objects are masked _by_ .GlobalEnv:      group, y  
> boxplot(y ~ group,  
+         main = "Boxplot of Daily Outcome by Temperature condition 1/2/3",  
+         xlab = "Factor Levels : Temperature condition 1/2/3",  
+         ylab = "Daily Outcome") 
>  

 

 

 

>  
> # descriptive statistics by group 
> tapply(y, group, summary) 
$`1`    
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    
50.50   50.90   51.45   51.46   52.05   52.40   
$`2`    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    
45.00   46.72   47.30   46.94   47.65   47.90   
$`3`    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    
44.90   45.92   46.30   46.35   47.10   47.20   
>
>
> detach(group_df) 
>  
> # one-wayANOVA 
> aov(y ~ group, data = group_df) 
Call:    aov(formula = y ~ group, data = group_df)  
Terms:                   
group Residuals Sum of Squares  156.302    19.393 
Deg. of Freedom       2        27  
Residual standard error: 0.8475018 Estimated effects may be unbalanced 

 

> summary(aov(y ~ group, data = group_df))             
Df Sum Sq Mean Sq F value  Pr(>F)     
group        2 156.30   78.15   108.8 1.2e-13 *** Residuals   27  19.39    0.72                     
--- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

일원분산분석(one-way ANOVA) 결과 검정통계량 F-value가 108.8으로 나왔으며, P-value 값은 '1.2e-13'으로서 매우 작게 나와 유의수준 10% 에서 귀무가설을 기각하고 대립가설을 채택하게 되어 "온도 조건 1/2/3에 따라서 휘발유 생산량에 차이가 있다"라고 판단할 수 있겠습니다.

 

 오차의 등분산성 가정에 대해 Bartlett 검정 결과 P-value가 0.4368로서 유의수준 10%에서 귀무가설을 채택하여 "오차의 등분산성 가정을 만족한다"고 할 수 있겠습니다.

 

많은 도움 되었기를 바랍니다.

 

질문은 댓글로 남겨주세요.

 

  • 쌍을 이룬 집단 간 평균 다중비교 (multiple comparisons)

Tukey's HSD(honestly significant difference) test 참조

Duncan's LSR(least significant range) test 참고

 

 

  • 대비 (contrast)

샤페 검정법 (scheffe test) 참고

 

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

728x90
반응형
Posted by Rfriend
,

두개의 모집단에 대한 추정과 검정 (two sample tests)에 대해서 정규분포 가정을 만족하는 경우와 그렇지 않은 경우로 나누어서 알아보겠습니다.

 

 

정규성 가정을 충족하는 경우

 

(1) 독립된 두 표본의 평균 차이에 대한 추정과 검정 : t.test()

     (indepentent two sample t-test)

 

(2) 짝을 이룬 표본에 대한 평균 차이에 대한 추정과 검정 : t.test(paired=TRUE)

     (paired sample t-test)

 

(3) 두 모집단의 모비율 차이에 대한 추정과 검정 : prop.test()

     (independent two population proportions test)

 

 

정규성 가정을 충족하지 못하는 경우, 혹은 분포형태를 모르는 경우

 

(4) 두 모집단의 중심 차이에 대한 비모수 검정 : wilcox.test()

    (non-parametric wilcoxon tests on two indepedent sample)

 

 

그동안의 포스팅에서는 "정규성 가정을 충족하는" 두 모집단의 모평균 혹은 모비율 차이에 대한 추정과 검정을 다루었는데요, 이번 포스팅에서는 정규성 가정을 만족하지 않거나 혹은 분포형태를 모르는 경우의 (4) 두 모집단의 중심 차이에 대한 비모수 검정(non-parametric wilcoxon tests on two indepedent sample), R의 wilcox.test() 함수에 대해 소개하겠습니다.

(cf. 정규성을 만족하면 t-test를 실시)

 

모수 검정과 비모수 검정법은 아래 비교표처럼 서로 짝을 지어서 보면 좀더 쉽게 이해할 수 있을 것입니다. 

 

 

 

[ 모수 검정과 비모수 검정 비교표 ]

 

 구분

 모수 검정

(Parametric Test)

비모수 검정

(Nonparametric Test) 

 When to use

 정규분포 가정 만족 시

 (normal distribution)

 정규분포 가정 불충족 시,

 혹은 모집단이 어떤 분포를 따르는지 모를 때

 (non-normal distribution, or un-know distribution,

  or very small sample size, or rankded data)

Statisctic

 평균 (mean)

 중앙값 (median) 

 1 sample

 1 sample t-test 

 1 sample Wilcoxon signed rank test

 2 samples

 2 sample t-test

 Wilcoxon rank sum test,  

 Mann-Whitney U-test

 paired 2-sample t-test

 Wilcoxon signed rank test

 more than

2 samples

 one-way ANOVA

 Kruskal-Wallis test

 

이론적인 부분에 대해서는 wikipedia (☞ https://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U_test) 를 참조하시기 바랍니다.

 

두 집단이 독립적이면 (1) Wilcoxon rank sum test (혹은 Mann-Whitney U-test), 짝을 이룬 집단이면 (2) Wilcoxon signed rank test 를 사용하게 되는되요, R에서는 wilcox.test() 함수를 동일하게 사용하고, 옵션에서 "paired = TRUE" 여부에 따라서 짝을 이룬 표본 여부를 구분하게 됩니다.  순서대로 예를 들어 설명해보겠습니다.

 

 

 

(1) independent 2 samples : Wilcoxon rank sum test (혹은 Mann-Whitney U-test)

     => wilcox.test()

 

정규분포를 따르지 않는 두 모집단에 중심에 대한 비모수 검정에 사용하는 R Wilcoxon rank sum test를 실시하면 되며, 아래에 R의 wilcox.test() 함수를 사용한 예를 소개하였습니다. 

 

 

[ R wilcox.test() 함수 사용방법 2가지 ] 

 

 

 

예제로 사용할 데이터는 MASS 패키지에 내장된 Cars93 데이터프레임의 가격(Price)과 생산국가(Origin) 입니다. 생산국이 USA vs. non-USA 2개의 group 에 대해서 차 가격(Price)의 평균이 차이가 있는지를 검정해보겠습니다.

 

> # Dataset
> library(MASS)
> str(Cars93)
'data.frame':	93 obs. of  27 variables:
 $ Manufacturer      : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4 4 4 5 ...
 $ Model             : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1 6 24 54 74 73 35 ...
 $ Type              : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3 2 2 3 2 ...
 $ Min.Price         : num  12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ...
 $ Price             : num  15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ...
 $ Max.Price         : num  18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3 36.3 ...
 $ MPG.city          : int  25 18 20 19 22 22 19 16 19 16 ...
 $ MPG.highway       : int  31 25 26 26 30 31 28 25 27 25 ...
 $ AirBags           : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2 2 2 2 2 2 ...
 $ DriveTrain        : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3 2 2 ...
 $ Cylinders         : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4 4 4 5 ...
 $ EngineSize        : num  1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ...
 $ Horsepower        : int  140 200 172 172 208 110 170 180 170 200 ...
 $ RPM               : int  6300 5500 5500 5500 5700 5200 4800 4000 4800 4100 ...
 $ Rev.per.mile      : int  2890 2335 2280 2535 2545 2565 1570 1320 1690 1510 ...
 $ Man.trans.avail   : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1 1 ...
 $ Fuel.tank.capacity: num  13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ...
 $ Passengers        : int  5 5 5 6 4 6 6 6 5 6 ...
 $ Length            : int  177 195 180 193 186 189 200 216 198 206 ...
 $ Wheelbase         : int  102 115 102 106 109 105 111 116 108 114 ...
 $ Width             : int  68 71 67 70 69 69 74 78 73 73 ...
 $ Turn.circle       : int  37 38 37 37 39 41 42 45 41 43 ...
 $ Rear.seat.room    : num  26.5 30 28 31 27 28 30.5 30.5 26.5 35 ...
 $ Luggage.room      : int  11 15 14 17 13 16 17 21 14 18 ...
 $ Weight            : int  2705 3560 3375 3405 3640 2880 3470 4105 3495 3620 ...
 $ Origin            : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1 1 1 1 ...
 $ Make              : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ... 

 

 

먼저, y(vector) ~ Factor, data=data.frame 형식을 먼저 소개합니다.

 

> 
> ##----------------------------------------------------------
> ## independent two population proportions test : prop.test()
> ##----------------------------------------------------------
> 
> ##--- way 1 : y ~ Factor
> wilcox.test(Price ~ Origin, 
+             data=Cars93, 
+             alternative = c("two.sided"), 
+             mu = 0, 
+             conf.int = FALSE, 
+             conf.level = 0.95)

	Wilcoxon rank sum test with continuity correction

data:  Price by Origin
W = 1024.5, p-value = 0.6724
alternative hypothesis: true location shift is not equal to 0

Warning message:
In wilcox.test.default(x = c(15.7, 20.8, 23.7, 26.3, 34.7, 40.1,  :
  tie가 있어 정확한 p값을 계산할 수 없습니다

 

 

P-value가 0.6724로서 유의수준 10%에서 귀무가설(H0 : Origin(USA와 non-USA) 별 두 집단 간에 자동차 가격 차이는 없다)을 채택하게 되었습니다.

 

 

두번째로, x, y numeric vector를 이용하는 방법은 아래와 같습니다.

 

> ##--- way 2. x, y numeric vectors
> 
> # x, y numeric vector indexing
> Price_USA <- Cars93[which(Cars93$Origin == c("USA")), c("Price")]
> Price_nonUSA <- Cars93[which(Cars93$Origin == c("non-USA")), c("Price")]
> wilcox.test(Price_USA, Price_nonUSA, 
+             alternative = c("two.sided"), 
+             mu = 0, 
+             conf.int = FALSE, 
+             conf.level = 0.95)

	Wilcoxon rank sum test with continuity correction

data:  Price_USA and Price_nonUSA
W = 1024.5, p-value = 0.6724
alternative hypothesis: true location shift is not equal to 0

Warning message:
In wilcox.test.default(Price_USA, Price_nonUSA, alternative = c("two.sided"),  :
  tie가 있어 정확한 p값을 계산할 수 없습니다

 

 

W통계량 및 P-value가 위의 첫번째 방법과 동일함을 알 수 있습니다.

 

 

 

 

다음으로 짝을 이룬 표본에 대한 비모수 검정 소개합니다. paired sample t-test 포스팅에서 사용했던 예제를 동일하게 사용하겠습니다.

 

(2) paired 2-sample : Wilcoxon signed rank test => wilcox.test(, paired = TRUE, ...)

 

Example 1 ) 새로운 당뇨병 치료제를 개발한 제약사의 예를 계속 들자면, 치료에 지대한 영향을 주는 외부요인을 통제하기 위해 10명의 당뇨병 환자를 선별하여 1달 동안 '위약(placebo)'을 투여한 기간의 혈압 (Xi)과 '신약(new medicine)'을 투여한 1달 기간 동안의 혈당 수치(Yi)를 측정하여 짝을 이루어 혈당 차이를 유의수준 5%에서 비교하는 방법이 짝을 이룬 표본에 대한 검정이 되겠습니다. 단, 모집단이 정규분포를 따르지 않는다고 합니다.

 

 

[ 환자 10명에 대한 당뇨병 치료제 투약 전/후 혈당 비교(before/after paired sample comparison) ]

  

(예를 들어 설명하기 위해 임의로 가공해서 만든 가짜 수치임.  혈당에 대해서는 아무것도 모름 ^^;)

 

 

귀무가설 H0 : 당뇨병 치료제는 효과가 없다 (mu1 = mu2, ie. difference = 0)

대립가설 H1 : 당뇨병 치료제는 효과가 있다 (혈당을 낮춘다, mu1 > mu2) => 우측검정(right-sided test)

 

위 문제를 R의 wilcox.test(paired=TRUE) 함수를 사용해서 풀어보면 아래와 같습니다.

 

당뇨병 치료제가 효과가 있어서 혈당을 낮추었는지를 검정하는 것이므로 'alternative = c("greater") 옵션(mu1 > mu2)을 입력하였습니다. P-value 가 0.006172 이므로 유의수준 (significance level) 5% 에서 귀무가설을 기각하고 대립가설(치료제가 효과가 있음. 혈당을 낮추었음)을 채택하게 됩니다.

(cf. paired sample t-test 에서도 대립가설 채택)

 

> ##----------------------------------------------------------
> ## paired sample : wilcox signed rank test
> ##----------------------------------------------------------
> 
> # paired 10 sample of patient's blood sugar
> x1 <- c(51.4, 52.0, 45.5, 54.5, 52.3, 50.9, 52.7, 50.3, 53.8, 53.1)
> x2 <- c(50.1, 51.5, 45.9, 53.1, 51.8, 50.3, 52.0, 49.9, 52.5, 53.0)
> # wilcox signed rank test
> wilcox.test(x1, x2, 
+        alternative = c("greater"), 
+        paired = TRUE, 
+        conf.level = 0.95)

	Wilcoxon signed rank test with continuity correction

data:  x1 and x2
V = 52.5, p-value = 0.006172
alternative hypothesis: true location shift is greater than 0

Warning message:
In wilcox.test.default(x1, x2, alternative = c("greater"), paired = TRUE,  :
  tie가 있어 정확한 p값을 계산할 수 없습니다

 

 

많은 도움이 되었기를 바랍니다. 

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

 

728x90
반응형
Posted by Rfriend
,

두개의 모집단에 대한 추정과 검정 (two sample tests)에 대해서 정규분포 가정을 만족하는 경우와 그렇지 않은 경우로 나누어서 알아보겠습니다.

 

 

정규성 가정을 충족하는 경우

 

(1) 독립된 두 표본의 평균 차이에 대한 추정과 검정 : t.test()

     (indepentent two sample t-test)

 

(2) 짝을 이룬 표본에 대한 평균 차이에 대한 추정과 검정 : t.test(paired=TRUE)

     (paired sample t-test)

 

(3) 두 모집단의 모비율 차이에 대한 추정과 검정 : prop.test()

     (independent two population proportions test)

 

 

 

정규성 가정을 충족하지 못하는 경우, 혹은 분포형태를 모르는 경우

 

(4) 두 모집단의 중심 차이에 대한 비모수 검정 : wilcox.test()

    (non-parametric wilcoxon tests on two indepedent sample)

 

 

 

지난번 포스팅에서는 독립된 두 표본, 짝을 이룬 표본에 대한 평균 차이에 대한 추정과 검정을 했는데요, 이번 포스팅에서는 정규성 가정을 만족하는 분포에서 (3) 두 모집단의 모비율 차이에 대한 추정과 검정(independent two population proportions test), R의 prop.test() 함수에 대해 소개하겠습니다.

 

두 집단의 모비율에 대한 추정과 검정은 범주형 변수(categorical variable)에 대해 두 집단의 비율의 차이를 비교해보려고 할 때 사용하게 됩니다.  가령, 두 집단 간의 특정 정당 혹은 후보에 대한 지지율의 차이, 흡연자와 비흡연자의 폐암 발병률의 차이, A학교와 B학교 야구부 선수들의 평균 타율의 차이 등과 같이 yes/no, success/failure 로 구분이 되는 범주형 변수(categorical variable)에 대해 집단 간 차이가 있는지를 통계량을 사용해서 P-value 를 가지고 (혹은 채택역, 기각역을 사용해서) 검정을 하게 됩니다.

 

 

[ 두 모집단의 모비율 차이에 대한 추정과 검정통계량 ]

 

 

 

 

R의 prop.test() 함수의 기본적인 사용법은 아래와 같습니다.

 

 

[ R prop.test() 함수 사용법 ]

 

 

 

 

아래 예를 R porp.test() 함수를 사용해서 풀어보도록 하겠습니다.

 

 

예제 )  A회사 직장인 500명과 B회사 직장인 600명을 대상으로 조사를 한 결과, A회사 직장인의 흡연율은 33%, B회사 직장인의 흡연율은 41%로 나타났다. 그러면 A회사와 B회사 직장인의 흡연율(proportion of smokers)에는 차이가 있다고 할 수 있는지 유의수준 (significance level) 5% 에서 검정하시오.

 

 

 

귀무가설 Ho : A회사와 B회사의 흡연율은 차이가 없다 (p1 - p2 = 0)

대립가설 H1 : A회사와 B회사의 흡연율은 차이가 있다 (p1 - p2 != 0)

 

 

prop.test(x, ...)의 x에는 'number of events'를 집어넣어야 하므로 아래 문제에서는 비율(proportion) 과 관측값수(n)을 곱해서 x <- prop*n 으로 계산을 해서 새로운 벡터 x를 생성해서 prop.test(x,...)에 집어넣었습니다.

 

> ##----------------------------------------------------------
> ## independent two population proportions test : prop.test()
> ##----------------------------------------------------------
> 
> prop <- c(0.33, 0.41) # proportion of events
> n <- c(500, 600) # number of trials
> x <- prop*n # number of events
> x 
[1] 165 246
> 
> 
> prop.test(x = x, # number of events
+           n = n, # number of trials
+           alternative = c("two.sided"), # alternative = c("two.sided", "less", "greater")
+           conf.level = 0.95) # confidence level (= 1- significance level 'alpha')

	2-sample test for equality of proportions with continuity correction

data:  x out of n
X-squared = 7.1203, df = 1, p-value = 0.007622
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.13881966 -0.02118034
sample estimates:
prop 1 prop 2 
  0.33   0.41

 

 

P-value가 0.007622 이므로 유의수준(significance level) 5%에서 귀무가설을 기각하고 대립가설을 채택하여 A회사와 B회사의 흡연율에 차이가 있다고 할 수 있겠습니다. A회사보다 B회사가 스트레스를 더 받는 회사인 모양입니다. ^^; 

 

많은 도움 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

728x90
반응형
Posted by Rfriend
,

두개의 모집단에 대한 추정과 검정 (two sample tests)에 대해서 정규분포 가정을 만족하는 경우와 그렇지 않은 경우로 나누어서 알아보겠습니다.

 

 

정규성 가정을 충족하는 경우

 

(1) 독립된 두 표본의 평균 차이에 대한 추정과 검정 : t.test()

     (indepentent two sample t-test)

 

(2) 짝을 이룬 표본에 대한 평균 차이에 대한 추정과 검정 : t.test(paired=TRUE)

     (paired sample t-test)

 

(3) 두 모집단의 모비율 차이에 대한 추정과 검정 : prop.test()

     (independent two population proportions test)

 

 

 

정규성 가정을 충족하지 못하는 경우, 혹은 분포형태를 모르는 경우

 

(4) 두 모집단의 중심 차이에 대한 비모수 검정 : wilcox.test()

    (non-parametric wilcoxon tests on two indepedent sample)

 

 

 

지난번 포스팅에서는 독립된 두 표본에 대한 평균 차이에 대한 추정과 검정을 했는데요, 이번 포스팅에서는 정규성 가정을 만족하는 분포에서 (2) 짝을 이룬 두 표본의 평균 차이에 대한 추정과 검정(paired sample t-test)를 R의 t.test(paired=TRUE) 함수에 대해 소개하겠습니다.

 

 언제 짝을 이룬 표본에 대한 평균 차이에 대한 추정과 검정이 필요한지가 궁금할 텐데요, 두 표본이 상관관계가 있을 때, 전/후 비교 (paired before-after comparison)을 할 때 사용합니다. 

 

 

 

 

가령, 새로운 당뇨병 의약품을 만든 제약사가 신약의 효과를 측정하고 싶을 텐데요, 치료의 경과와 수준은 환자의 건강상태, 성별, 나이, 식생활 환경, 체질, 스트레스 수준,... 통제해야 할 수준이 무척 많으므로 치료의 차도가 있다고 하더라도 이게 신약 때문에 효과가 있는 것인지, 아니면 앞서 말씀드린 신약 외의 다양한 외부변수에 의한 치료 개선효과인지 분간하기가 힘들게 됩니다.  이런 문제를 해결하기 위해서, 즉 외부 요인을 최대한 통제할 수 있는 방안으로 검증하고자 하는 요인 외의 요인은 모두 동일하게 만들 수 있는 방법이 '짝을 이룬 두 표본(paired sample)'에 대해서 검증하고자 하는 요인만을 달리해서 그 차이(paired comparison)를 가지고 검정을 하는 방법이 있습니다.  이해를 돕기 위해서 아래에 2개의 예를 들어보았습니다.

 

 

 

 

 

Example 1 ) 새로운 당뇨병 치료제를 개발한 제약사의 예를 계속 들자면, 치료에 지대한 영향을 주는 외부요인을 통제하기 위해 10명의 당뇨병 환자를 선별하여 1달 동안 '위약(placebo)'을 투여한 기간의 혈압 (Xi)과 '신약(new medicine)'을 투여한 1달 기간 동안의 혈당 수치(Yi)를 측정하여 짝을 이루어 혈당 차이를 유의수준 5%에서 비교하는 방법이 짝을 이룬 표본에 대한 검정이 되겠습니다.

 

 

[ 환자 10명에 대한 당뇨병 치료제 투약 전/후 혈당 비교(before/after paired sample comparison) ]

 

(예를 들어 설명하기 위해 임의로 가공해서 만든 가짜 수치임.  혈당에 대해서는 아무것도 모름 ^^;)

 

 

귀무가설 H0 : 당뇨병 치료제는 효과가 없다 (mu1 = mu2, ie. difference = 0)

대립가설 H1 : 당뇨병 치료제는 효과가 있다 (혈당을 낮춘다, mu1 > mu2) => 우측검정(right-sided test)

 

위 문제를 R의 t.test(paired=TRUE) 함수를 사용해서 풀어보면 아래와 같습니다.

 

위의 t-검정통계량 공식에 따라서 계산한 아래 값과, t.test(paired=TRUE) 함수로 계산한 t-검정통계량 값이 3.5507로 서로 같음을 알 수 있습니다.  당뇨병 치료제가 효과가 있어서 혈당을 낮추었는지를 검정하는 것이므로 'alternative = c("greater") 옵션(mu1 > mu2)을 입력하였습니다. P-value 가 0.003105 이므로 유의수준 (significance level) 5% 에서 귀무가설을 기각하고 대립가설(치료제가 효과가 있음. 혈당을 낮추었음)을 채택하게 됩니다.

 

> 
> ##----------------------------------------------------------
> ## paired sample t-test
> ##----------------------------------------------------------
> 
> # paired 10 sample of patient's blood sugar
> x1 <- c(51.4, 52.0, 45.5, 54.5, 52.3, 50.9, 52.7, 50.3, 53.8, 53.1)
> x2 <- c(50.1, 51.5, 45.9, 53.1, 51.8, 50.3, 52.0, 49.9, 52.5, 53.0)
> 
> 
> # difference of x1, x2
> diff_x <- x1 - x2 
> diff_x
 [1]  1.3  0.5 -0.4  1.4  0.5  0.6  0.7  0.4  1.3  0.1
> 
> # mean of difference
> mean_diff_x <- mean(diff_x) 
> mean_diff_x
[1] 0.64
> 
> # standard deviation of difference
> sd_diff_x <- sd(diff_x) 
> sd_diff_x
[1] 0.5699903
> 
> # t statistics
> t_x <- mean_diff_x/(sd_diff_x / sqrt(length(diff_x))) 
> t_x
[1] 3.550688
> 
> # paired sample t-test : greater
> t.test(x1, x2, 
+        alternative = c("greater"), 
+        paired = TRUE, 
+        conf.level = 0.95)

	Paired t-test

data:  x1 and x2
t = 3.5507, df = 9, p-value = 0.003105
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 0.3095874       Inf
sample estimates:
mean of the differences 
                   0.64

 

 

 

 

 

Example 2) 또 하나의 예를 들자면 두 종류의 신발 밑창의 원재료가 닳는 정도가 차이가 있는지를 검정하기 위해서 10명의 소년에게 한쪽은 A라는 원재료로 만든 신발을 신기고, 다른 한쪽은 B라는 원재료로 만든 신발을 신긴 후에, 일정 기간이 지난후에 신발을 수거하여 10명의 각 소년의 왼쪽 신발 밑창의 닳은 정도와 오른쪽 신발 밑창의 닳은 정도의 차이를 비교하여 두 종류 원재료의 재질이 다른지를 검정하는 방법이 짝을 이룬 표본에 대한 검정에 해당되겠습니다.

 

 

[ 신발 밑창 원재료 A, B별 닳은 정도 비교 (shoes material paired sample comparison) ]

 

 

귀무가설 H0 : 두 신발 원재료 A, B는 닳는 정도가 같다 (mu1 = mu2, ie. difference = 0)

대립가설 H1 : 두 신발 원재료 A, B는 닳는 정도가 다르다(mu1 != mu2, ie. difference != 0)
                    => 양측검정(two-sided test)

 

유의수준(significance level) 5%에서 위 가설을 검정하시오.

 

 

이번 예제는 MASS 패키지에 들어있는 shoes list 데이터셋을 이용해서 분석을 해보겠습니다.  위의 t-검정통계량 공식에 넣어 계산한 값과 아래의 R을 통해 계산한 t-검정통계량 값이 -3.3489 로서 서로 같음을 알 수 있습니다.  두 짝을 이룬 표본의 차이의 평균은 -0.41 이며, t-test 결과의 P-value가 0.008539 로서 유의수준 5% 하에서 귀무가설을 기각하고 대립가설(A, B 두 원재료는 차이가 있다)을 채택하게 되었습니다.

 

> ##-- Example 2
> 
> library(MASS)
> 
> str(shoes)
List of 2
 $ A: num [1:10] 13.2 8.2 10.9 14.3 10.7 6.6 9.5 10.8 8.8 13.3
 $ B: num [1:10] 14 8.8 11.2 14.2 11.8 6.4 9.8 11.3 9.3 13.6
> 
> # difference of shoes A and shoes B
> diff_shoes <- shoes$A - shoes$B
> diff_shoes
 [1] -0.8 -0.6 -0.3  0.1 -1.1  0.2 -0.3 -0.5 -0.5 -0.3
> 
> 
> # mean of difference
> mean_diff_shoes <- mean(diff_shoes)
> mean_diff_shoes
[1] -0.41
> 
> 
> # standard deviation of difference
> sd_diff_shoes <- sd(diff_shoes)
> sd_diff_shoes
[1] 0.3871549
> 
> 
> # t statistics
> t_shoes <- mean_diff_shoes/ (sd_diff_shoes/sqrt(length(diff_shoes)))
> t_shoes
[1] -3.348877
> 
> 
> # paired sample t-test : two.sided
> t.test(shoes$A, shoes$B, 
+        alternative = c("two.sided"), 
+        paired = TRUE, 
+        conf.level = 0.95)

	Paired t-test

data:  shoes$A and shoes$B
t = -3.3489, df = 9, p-value = 0.008539
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.6869539 -0.1330461
sample estimates:
mean of the differences 
                  -0.41 

 

 

많은 도움 되었기를 바랍니다.

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^

 

 

728x90
반응형
Posted by Rfriend
,