t-test를 수행해보자.

코드는 곽기영, R을 이용한 통계데이터분석, 도서출판 청람을 참고했다.

필요한 라이브러리를 부른다.

library(MASS) #데이터 때문에
library(dplyr) #데이터 처리 때문에
library(gplots) #플롯 때문에
library(tidyr) #데이터 처리 때문에
library(reshape2) #데이터 처리 때문에
head(cats) #Cat 데이터
str(cats)
## 'data.frame':    144 obs. of  3 variables:
##  $ Sex: Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Bwt: num  2 2 2 2.1 2.1 2.1 2.1 2.1 2.1 2.1 ...
##  $ Hwt: num  7 7.4 9.5 7.2 7.3 7.6 8.1 8.2 8.3 8.5 ...

성별에 따라 Bwt, Hwt와 같이 무게(weight)가 다를 것 같다. 분포가 정규분포에 가깝고 모집단의 정보(평균)이 사전에 알려져 있다고 가정하자.

T-test

평균이 2.6인가?

with(
  cats,{
    t.test(x=Bwt,mu=2.6)
  }
)
## 
##  One Sample t-test
## 
## data:  Bwt
## t = 3.0565, df = 143, p-value = 0.002673
## alternative hypothesis: true mean is not equal to 2.6
## 95 percent confidence interval:
##  2.643669 2.803553
## sample estimates:
## mean of x 
##  2.723611

Bwt의 평균이 2.6이라는 가설은 기각된다(p-value < 0.01). 즉, 평균은 2.6이 아니다. 계산결과를 보면 평균이 2.72인 것을 알 수 있다.

평균이 2.7인가?

with(cats,{
  t.test(Bwt,mu=2.7)
})
## 
##  One Sample t-test
## 
## data:  Bwt
## t = 0.58382, df = 143, p-value = 0.5603
## alternative hypothesis: true mean is not equal to 2.7
## 95 percent confidence interval:
##  2.643669 2.803553
## sample estimates:
## mean of x 
##  2.723611

Bwt의 평균이 2.7이라는 가설을 기각할 수 없다(p-value >= 0.1). 샘플의 평균은 2.72로 2.7과 유사하다.

평균이 2.6보다 큰가?

with(cats,{
  t.test(Bwt,mu=2.6,alternative='greater')
})
## 
##  One Sample t-test
## 
## data:  Bwt
## t = 3.0565, df = 143, p-value = 0.001337
## alternative hypothesis: true mean is greater than 2.6
## 95 percent confidence interval:
##  2.656656      Inf
## sample estimates:
## mean of x 
##  2.723611

p-value와 신뢰구간을 따로 보면?

cats.t <- with(cats,{
  t.test(Bwt,mu=2.6)
})
data.frame(
  `p-value`=cats.t$p.value,
  `confident interval`=cats.t$conf.int
) %>% t()
##                           [,1]        [,2]
## p.value            0.002673036 0.002673036
## confident.interval 2.643669291 2.803552932

Two-Independent Samples T-Test

Welch two sample t-test

t.test(Bwt ~ Sex,data=cats)
## 
##  Welch Two Sample t-test
## 
## data:  Bwt by Sex
## t = -8.7095, df = 136.84, p-value = 8.831e-15
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6631268 -0.4177242
## sample estimates:
## mean in group F mean in group M 
##        2.359574        2.900000

성별에 따라 Bwt에 차이가 있는가? 차이가 없다는 가설이 기각된다(p-value < 0.01). 성별에 따른 Bwt가 계산되어 있다. 여성은 2.3, 남성은 2.9이다.

bars <- tapply(cats$Bwt, cats$Sex, mean)
lower <- tapply(cats$Bwt, cats$Sex, function(x) t.test(x)$conf.int[1])
upper <- tapply(cats$Bwt, cats$Sex, function(x) t.test(x)$conf.int[2])
windows(width=4.0, height=5.5)
barplot2(bars, space=0.4, ylim=c(0, 3.0),
         plot.ci=TRUE, ci.l=lower, ci.u=upper, ci.color="maroon", ci.lwd=4, 
         names.arg=c("Female", "Male"), col=c("coral", "darkkhaki"),
         xlab="Cats", ylab="Body Weight (kg)", 
         main="Body Weight by Sex\nwith Confidence Interval")

비율분석(proportion test)

smokers  <- c(83, 90, 129, 70)
patients <- c(86, 93, 136, 82)

prop.test(x=smokers, n=patients)
## 
##  4-sample test for equality of proportions without continuity
##  correction
## 
## data:  smokers out of patients
## X-squared = 12.6, df = 3, p-value = 0.005585
## alternative hypothesis: two.sided
## sample estimates:
##    prop 1    prop 2    prop 3    prop 4 
## 0.9651163 0.9677419 0.9485294 0.8536585

4개의 실험 집단이 있고 각각 흡연자와 비흡연자 수를 기록했다. 흡연자와 비흡연자의 비율은 과연 일정할까? 일정하다는 가설이 깨진다(p-value < 0.01). 결과를 보면 4번째 그룹의 비율(85%)와 1,2번째 그룹의 비율(96.5%, 96.8%)에 큰 격차가 존재한다.

Paired-Samples t-test

t.test(extra ~ group, data=sleep,paired=TRUE)
## 
##  Paired t-test
## 
## data:  extra by group
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.4598858 -0.7001142
## sample estimates:
## mean of the differences 
##                   -1.58

같은 대상에게 두 개의 처리그룹으로 실험을 수행했다. 각 처리에 따른 차이가 존재할까? 가설은 존재하지 않는다이지만 이것이 기각된다(p-value < 0.01). 즉, group에 따라 extra 측면에서 유의미한 차이가 있다.

head(sleep)

sleep 데이터를 wide form으로 바꾸어서 처리해보자.

sleep.wide <- spread(sleep,key=group,value=extra)
head(sleep.wide)
with(sleep.wide,{
  t.test(`1`,`2`,paired=TRUE)
})
## 
##  Paired t-test
## 
## data:  1 and 2
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.4598858 -0.7001142
## sample estimates:
## mean of the differences 
##                   -1.58

가정

다음의 가정이 충족되어야 한다.

  1. 독립성(paired일 경우에는 제외)
  2. 정규성(만족되지 않으면 Mann-Whitney test를 고려한다, 제일 중요)
  3. 등분산성(paired일 때 고려, 가정이 만족되지 않으면 이분산성을 고려한 결과 보고)
with(cats,{
  shapiro.test(Bwt)
})
## 
##  Shapiro-Wilk normality test
## 
## data:  Bwt
## W = 0.95188, p-value = 6.731e-05

정규분포를 따른다는 가설이 지지되지 않는다(p-value < 0.01). 따라서 정규성 가정이 깨졌다.

with(cats,{
  bartlett.test(Bwt,Sex)
})
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Bwt and Sex
## Bartlett's K-squared = 15.075, df = 1, p-value = 0.0001033

등분산 가정을 만족한다는 가설이 지지되지 않는다(p-value < 0.01). 띠라서 등분산 가정이 깨졌다.

with(cats,{
  t.test(Bwt,var.equal=FALSE)
})
## 
##  One Sample t-test
## 
## data:  Bwt
## t = 67.346, df = 143, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  2.643669 2.803553
## sample estimates:
## mean of x 
##  2.723611

이분산 가정(var.equal=FALSE)로 처리를 했다.

정규분포가 충족되지 않을 경우에 어떻게 처리하면 좋을까? 제일 좋은 방법은 데이터를 많이 모으는 것이다. 현재 cats 데이터는 144개뿐이다. 그룹으로 나누면 잘 해야 한 그룹에 70개 정도뿐이다. 데이터가 충분하다면 Boxcox 변환을 사용하자.

with(cats,{
  hist(Bwt)
})

Bwt의 분포가 이상하다. Boxcox에서 \(\lambda=0\)이라 가정하면 \(log()\) 변환으로 처리하면 된다.

with(cats,{
  hist(log(Bwt))
})

정규분포 가정이 어떻게 될까? 아까보다는 좋다. 여전히 부족하지만.

with(cats,
     {
       x=log(Bwt)
       shapiro.test(x)
     })
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.96456, p-value = 0.0008872

MASS 패키지의 boxcox() 함수를 사용해보자.

#boxcox() - MASS
lam <- {
  r=boxcox(Bwt~Sex,data=cats)
  r$x[which.max(r$y)]
}

BwtL <- with(cats,{
  (Bwt^lam-1)/lam
})
catsNew <- cbind(cats,BwtL)

참고로 boxcox 변환을 통해 변수값을 변환해 계산하려면 다음과 같은 수식을 사용한다.

\[ \lambda=0, log(y) \] \[ \lambda \neq 0, \frac{y^{\lambda}-1}{\lambda} \]

변환 결과를 가지고 t-test를 수행해보자.

t.test(BwtL ~ Sex, data=catsNew, vars.equal=FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  BwtL by Sex
## t = -8.4079, df = 115.89, p-value = 1.228e-13
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.11382381 -0.07042139
## sample estimates:
## mean in group F mean in group M 
##       0.6182883       0.7104109