R에서 variance ratio test (variance test)를 하는 방법

종종 두 그룹에서 데이터의 분산이 다른지 검정하고 싶을 수 있다 (궁금증이 많으면 그럴수도 있다. 이해하시라). 몇몇사람들은 t-test에서 등분산을 가정할지 말지 결정하기 위해서 t-test를 수행하기 이전에 variance ratio test를 먼저 수행해 보기도 한다.

 

R로 variance ratio test를 수행해 보기 위해, t-test에서와 마찬가지로 datasets 패키지에 내장되어있는 airquality 데이터를 이용해보자.

 

data(airquality)
str(airquality)

## 'data.frame':    153 obs. of  6 variables:
##  $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
##  $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
##  $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
##  $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

airquality$Ozone_g <- ifelse(airquality$Ozone >= mean(airquality$Ozone, na.rm=T), "high", "low")
head(airquality)

##   Ozone Solar.R Wind Temp Month Day Ozone_g
## 1    41     190  7.4   67     5   1     low
## 2    36     118  8.0   72     5   2     low
## 3    12     149 12.6   74     5   3     low
## 4    18     313 11.5   62     5   4     low
## 5    NA      NA 14.3   56     5   5    <NA>
## 6    28      NA 14.9   66     5   6     low

R에서는 아래와 같이 var.test() 함수를 이용한다.

 

(var_test <- var.test(Temp ~ Ozone_g, data=airquality))

## 
##  F test to compare two variances
## 
## data:  Temp by Ozone_g
## F = 0.45283, num df = 43, denom df = 71, p-value = 0.006092
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2682961 0.7925954
## sample estimates:
## ratio of variances 
##          0.4528301

var.test() 결과 object(아래 예시에서는 var_test) 역시 아래와 같은 요소들을 가지고 있으니, 확인해 보자.

 

var_test$statistic # F값
var_test$parameter # dfs
var_test$p.value # F값에 대한 p-value
2*pf(q=var_test$statistic, df1=var_test$parameter[1], df2=var_test$parameter[2]) # F값에 대한 p-value
var_test$conf.int # 두 그룹 variance ratio에 대한 95% CI
var_test$estimate # variance ratio
(v1 <- tapply(airquality$Temp, airquality$Ozone_g, FUN=function(x) var(x))[1])
(v2 <- tapply(airquality$Temp, airquality$Ozone_g, FUN=function(x) var(x))[2])
v1/v2 # variance ratio
var_test$null.value # 두 그룹 variance ratio가 어떤값과 다른지 평가하기 위한 null value, default는 1
var_test$alternative # "two.sided" (default), "greater" or "less"

cat("F값: ", var_test$statistic
    , "\ndfs: ", var_test$parameter
    , "\nF값에 대한 p-value: ", var_test$p.value, 2*pf(q=var_test$statistic, df1=var_test$parameter[1], df2=var_test$parameter[2])
    , "\n두 그룹 variance ratio에 대한 95% CI: ", var_test$conf.int
    , "\nvariance ratio: ", var_test$estimate
    , "\nvariance1: ", (v1 <- tapply(airquality$Temp, airquality$Ozone_g, FUN=function(x) var(x))[1])
    , "\nvariance2: ", (v2 <- tapply(airquality$Temp, airquality$Ozone_g, FUN=function(x) var(x))[2])
    , "\nvariance ratio: ", v1/v2
    , "\n두 그룹 variance ratio가 어떤값과 다른지 평가하기 위한 null value(default는 1): ", var_test$null.value
    , "\ntwo.sided (default), greater or less: ", var_test$alternative
    )

## F값:  0.4528301 
## dfs:  43 71 
## F값에 대한 p-value:  0.006092383 0.006092383 
## 두 그룹 variance ratio에 대한 95% CI:  0.2682961 0.7925954 
## variance ratio:  0.4528301 
## variance1:  25.92178 
## variance2:  57.24394 
## variance ratio:  0.4528301 
## 두 그룹 variance ratio가 어떤값과 다른지 평가하기 위한 null value(default는 1):  1 
## two.sided (default), greater or less:  two.sided

 

이상으로 R에서 variance ratio test를 하는 방법을 알아보았다.