R studio로 Unpaired two sample t-test 해 보자!

티스토리에 Unpaired t-test까지 하는 과정을 정리해보려 한다. 본인은 R studio에 관련해선 정말 쌩 초보이기 때문에, 정말 아무것도 모르시더라도 내 글만 참고하면 바로 따라할 수 있을 정도로! 

1. 먼저 구글에 R studio Unpaired t-test라고 검색했더니 아래의 사이트가 나왔다.
http://www.sthda.com/english/wiki/unpaired-two-samples-t-test-in-r

Unpaired Two-Samples T-test in R - Easy Guides - Wiki - STHDA

Q: What is unpaired two-samples t-test?
A: The unpaired two-samples t-test is used to compare the mean of two independent groups.

위 대답 속에 이미 t-test와 관련한 하나의 조건이 들어가 있는데, 그것은 비교하려고 하는 두 그룹이 서로 독립적이어야 한다는 것이다. 그래서 independent t-test라고 부르기도 하나 보다. 바로 t-test를 수행하려나 싶었는데, 두 가지 조건이 필요하단다. 그것은 정규성과 등분산성이다. 

Unpaired two-samples t-test can be used only under certain conditions: 
· when the two groups of samples (A and B), being compared, are normally distributed. This can be checked using Shapiro-Wilk test.
· and when the variances of the two groups are equal. This can be checked using F-test.

그렇다면 Shapiro-Wilk test부터 해야겠다. 

2. Shapiro-Wilk test는 아래의 사이트를 참고했다.
http://www.sthda.com/english/wiki/normality-test-in-r

Normality Test in R - Easy Guides - Wiki - STHDA

다수의 statistical tests (correlation, regression, t-test, and analysis of variance (ANOVA))는 분석 대상인 데이터가 normal (or Gaussian or Gauss or Laplace-Gauss) distribution을 했다는 가정 하에 수행된다. 즉, 이런 statistical tests의 validity는 data의 distribution에 의존적이라는 것이고, 그래서 이 tests를 parametric tests라고 부르기도 한다. 
그러므로 위와 같은 parametric tests를 수행하기 위해선, 데이터가 normal distribution을 이루고 있는지를 먼저 검정하는 preleminary tests를 수행해야 한다. 만약 데이터가 normal distribution을 이루고 있다는 가정이 틀렸을 경우 (즉, normal distribution하지 않을 경우) non-paramatric tests를 사용할 것이 권장된다. 

그러면 이제 본론으로 들어가서 구체적으로 어떻게 Shapiro-wilk test를 R에서 돌릴 수 있는지 알아보자.

① Install required R packages
먼저 필요한 패키지들을 설치한다. 위 사이트에 들어가서 그대로 코드 복붙 하시면 설치가 되시겠다. 

② Load required R packages
설치된 패키지들을 로드해 준다. 두 가지 package인 dplyr과 ggpubr이다. 여기까지는 복붙만 하면 되니 어렵지 않다.

③ Import your data into R
데이터를 R로 불러들여야 한다. 그런데 대부분의 프로그램이 그렇듯, 데이터를 받아들이는 일정한 form이 있다. 가지고 있는 데이터를 아래의 가이드라인에 맞게 만든 후, 저장하면 된다.
http://www.sthda.com/english/wiki/best-practices-in-preparing-data-files-for-importing-into-r

Best Practices in Preparing Data Files for Importing into R - Easy Guides - Wiki - STHDA

나의 경우엔 따로 엑셀 시트를 만들어서 저장하는 방법보다는, 데이터를 엑셀로 작성한 후 메모장에 복사-붙여넣기 해서 만들었다. 엑셀의 데이터를 메모장에 붙여넣을 경우 tab으로 분리되어 붙여넣어지기 때문에, 손쉽게 tab-delimited txt의 형태로 데이터를 만들어 저장할 수 있다.

이렇게 저장된 데이터를 my_data에 input 해준다. 본인은 아래의 함수를 사용하였는데, 아래의 함수를 입력하면 친숙한 window가 열려서 파일의 위치를 찾을 수 있기에, 코드에 file의 위치와 이름을 입력하지 않아도 되었다.
my_data = read.delim(file.choose()) 

my_data를 console에 입력해서 데이터가 잘 들어갔는지 확인해보자. 참고로, 위처럼 쓰면 first line을 header로 처리해준다. read.delim의 syntax는 아래와 같은 듯.

read.delim(file, header=TRUE, sep="\t")
· file – A file location.
· header – Whether the first line describes the column names.
· sep – The table delimiter, often times a tab (\t) or comma.

④ Check your data
이 과정은 굳이 왜 있는지 모르겠는데, dplyr package 안에 있는 sample_n()이라는 함수를 이용해 내 data에서 10개의 rows를 무작위로 추출해서 보여주는 과정을 거친다.
set.seed(1234)
dplyr:: sample_n(my_data, 10)

⑤ Assess the normality of the data in R
먼저, sample size가 충분히 큰 경우 (n >30) data의 distribution을 점검하지 않고 그냥 parametric tests를 사용할 수 있다. 이는 중심 극한의 정리 (central limit theorem)에 근거하는데, 이는 데이터의 distribution이 sample의 수가 충분히 크다면 그 distribution은  normal인 경향이 강하다는 정리이다.

그러나 sample size가 그보다 작은 경우에는, 두 가지 접근법으로 normality를 체크하는 듯 하다.
· Visual methods: Density plot이나 Q-Q plot을 이용해 시각적으로 normality를 체크한다. 
- Density plot: distribution이 bell shaped인지 확인할 수 있게 해 준다. ggpubr이라는 library를 불러들인 후, ggdensity라는 함수를 이용해 분포를 확인한다. my_data 옆에 $을 붙이면 내 data의 변수 항목이 뜬다. main은 plot의 제목을 입력하는 것이고, xlab는 Y축의 이름을 입력하는 것이다.
ggdensity(my_data$height, main="Height plot", xlab="height")

- Q-Q plot (or quantile-quantile plot): 주어진 sample과 normal distribution 간의 correlation을 보여주며, 참고를 위해 45도의 reference line 또한 제시된다. 맞게 보는 건지 모르겠지만, 내 데이터는 대개 검은 영역 안에 들어가는 것처럼 보인다.
ggqqplot(my_data$height)

· Significance test
Visual inspection은 때로 normality에 대한 명확한 지표가 될 수 없어서 significance test를 이용해 sample distribution이 normal distribution과 비교한다. 이 때 null hypothesis는 "sample distribution is normal,"이므로 test의 결과가 0.05보다 낮을 경우 영가설이 기각되며, 이는 곧 데이터가 normal하지 않다는 의미이다. 
Kolmogorov-Smirnov (K-S) normality test와 Shapiro-Wilk's test의 두 가지 test가 있는데, 일반적으로는 Shapiro-Wilk's method가 널리 권장되며 K-S에 비해 더 높은 power를 제공한다. 

아래의 함수를 이용했다.
shapiro.test(my_data$height)

결과는 p-value가 0.3044로 영가설 채택! 이 정도면 plot을 봐도, significance test 결과를 봐도 normality를 충분히 만족한다고 이야기할 수 있을 듯 하다.

3. 다음으로 분산의 동질성을 검정하기 위해 F-test를 해 보자.
http://www.sthda.com/english/wiki/f-test-compare-two-variances-in-r

F-Test: Compare Two Variances in R - Easy Guides - Wiki - STHDA

두 그룹의 분산을 비교해보는 것은 몇 가지 상황에서 유용하다.
· two sample t-test를 하려고 할 때
· 새로운 측정 기법의 variability를 기존 측정 기법의 variability와 비교하고자 할 때. 새로운 측정 기법은 측정 치의 variability를 감소시키는가? (감소시키면 더 좋다는 뜻인 듯.)

F-test는 normal assumption으로부터의 이탈에 매우 민감하기 때문에, 이를 사용하려면 data가 normal distribution을 이루고 있어야 한다. 그러므로 일반적으로 Shapiro-Wilk test가 F-test에 선행된다. 만약 normality가 충족되지 않은 것으로 여겨질 때는 F-test 대신 Levene's test나 Fligner-Killeen test를 사용해야 한다. (이 두 tests는 상대적으로 normal assumption으로부터의 이탈에 덜 민감하다.) 

F-test는 아래의 세 가지 질문에 대해 각 질문을 영가설로 설정한 후 테스트할 수 있다.
· whether the variance of group A is equal to the variance of group B?
· whether the variance of group A is less than the variance of group B?
· whether the variance of group A is greater than the variance of group B?
첫 번째를 가설로 삼아 테스트하는 경우를 two-tailed tests라고 하며, 두 번째와 세 번째를 가설로 삼아 테스트하는 경우를 one-tailed tests라고 한다. 

F값은 아래와 같은 수식으로 산출된다.
$ F=\frac{S_{A}^{2}}{S_{B}^{2}} $
또한 degress of freedom은 numerator의 경우 해당 그룹의 표본수-1이고, denominator의 경우도 해당 그룹의 표본수-1이다. (이것을 설명하는 이유는 뒤에 F-test의 결과를 해석함에 있어 이해를 증진하기 위함이다.)

① R function
var.test() 함수는 두 그룹의 variances를 비교하기 위해 사용되며, 아래와 같은 방식으로 사용된다.
var.test (values ~ groups, data, alternative = "two.sided")
var.test (x, y, alternative = "two.sided")
두 번째 줄의 x, y는 numeric vectors인데 어떤 상황에서 사용해야 하는지는 정확히 모르겠다. alternative에 넣을 수 있는 값은 "two.sided" (default), "greater", or "less"이다. (테스트하고자 하는 영가설에 맞게 사용한다.)

② Import and check your data into R
이 단계는 위에서 설명했으니 생략한다.

③ Preleminary test to check F-test assumptions
이 부분은 F-test 이전에 Shapiro-Wilk test를 이용해 정규성 검정을 해야 한다는 부분이다.

④ Compute F-test
res.ftest = var.test (height ~ group, data = my_data)
위와 같이 입력하여 계산하였고 결과는 아래와 같았다. 

F test to compare two variances
data: height by group
F = 0.7066, num df = 13, denom df = 5, p-value = 0.565
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval: 0.1089162 2.6615408
sample estimates: ratio of variances 0.7066024

위의 결과를 해석해보면,
· statistic: the value of the F test statistic. 즉, F 값이 처음으로 제시되어 있는데, 이는 위 공식에서 제시한 대로 두 그룹의 분산을 한 그룹은 분자, 한 그룹은 분모에 놓아 계산한 것이다.
· parameter: the degrees of the freedom of the F distribution of the test statistic. 이 역시 위에서 제시한 대로, 각 그룹의 표본수에서 1을 뺀 값이 된다.
· p.value: the p-value of the test. 이 수치는 0.05 미만일 때 영가설을 기각하고 대립가설 (alternative hypothesis)을 채택한다. 즉, 0.05 미만이면 분산의 동질성이 없다고 판단한다. 
· conf.int: a confidence interval for the ratio of the population variances. 95% CV를 제시한 것이다.
· estimate: the ratio of the sample variances. 이것은 F값과 같은 의미인 듯 한데 굳이 제시되는 이유를 모르겠다. 

어쨌든, 내 데이터의 F값은 0.7066이었고 p-value는 0.565이기 때문에 영가설을 채택하여, 두 그룹의 분산은 동일하다고 할 수 있다. 

그러면 정규성 (normality)과 등분산성 (homoscedasticity)을 모두 검정해봤으니, 드디어 t-test를 해 보자.

4. 멀리 돌아서 다시 unpaired two-samples t-test로. 
T-test도 F-test와 유사하게 아래의 세 가지 질문에 대해 각 질문을 영가설로 설정한 후 테스트할 수 있다.
· whether the variance of group A is equal to the variance of group B?
· whether the variance of group A is less than the variance of group B?
· whether the variance of group A is greater than the variance of group B?
마찬가지로 첫 번째를 가설로 삼아 테스트하는 경우를 two-tailed tests라고 하며, 두 번째와 세 번째를 가설로 삼아 테스트하는 경우를 one-tailed tests라고 한다.

등분산성의 충족 여부에 따라 두 가지 formula를 사용할 수 있다.
· Classical t-test: 분산이 동일할 때 (homoscedasticity)
$ t=\frac{m_{A}-m_{B}}{\sqrt{\frac{S^{2}}{n_{A}}+\frac{S^{2}}{n_{B}}}} $
위 공식에서 m은 각 그룹의 평균을 의미하며, n은 각 그룹의 표본 수를 의미한다. S^2는 an estimator of the pooled variance of the two groups라고 하며, 계산은 아래와 같이 한다고 한다.
$ S^{2}=\frac{\sum(x-m_{A})^{2}+\sum (x-m_{B})^{2}}{n_{A}+n_{B}-2} $
또한 degrees of freedom (df)은 두 표본수를 더한 값에서 2를 빼면 된다.

· Welch t-statistic: 분산이 동일하지 않을 때 (heteroscedasticity)
$ t=\frac{m_{A}-m_{B}}{\sqrt{\frac{S_{B}^{2}}{n_{A}}+\frac{S_{B}^{2}}{n_{B}}}} $
위 공식에서 달라진 것은 S가 각 그룹의 standard deviation을 의미한다는 점이다. 또한 이 때 degrees of freedom (df)은 아래와 같이 계산한다.
$  df=\frac{(\frac{S_{B}^{2}}{n_{A}}+\frac{S_{B}^{2}}{n_{B}})}{(\frac{S_{B}^{4}}{n_{A}^{2}(n_{A}-1)}+\frac{S_{B}^{4}}{n_{B}^{2}(n_{B}-1)})} $

① Install ggpubr R package for data visualization
사이트에서 코드를 그대로 복붙하여 패키지를 설치해주자.

② R function to compute unpaired two-samples t-test
t.test() 함수는 두 독립적인 그룹의 평균을 비교할 때 사용되며, syntax는 아래와 같다.
t.test(x, y, alternative="two.sided", var.equal=TRUE)
x, y는 numeric vectors이고, alternative에는 "two.sided" (default), "greater" or "less"가 입력될 수 있고, var.equal은 두 분산을 같게 취급할 것인지 여부인데 TRUE인 경우 pooled variance가 사용될 것이고, FALSE일 경우 Welch test가 사용될 것이다.

③ Import your data into R
위에서 했던 방법대로 data frame을 짜서 메모장에 넣어 불러들이면 된다.

④ Check your data
print(my_data)
잘 프린트되는지 확인해보자.
또한 dplyr 패키지를 이용해 group별로 summary statistics를 제시할 수 있다고 한다.
library(dplyr)
group_by(my_data, group) %>%
summarise(
count = n(),
mean = mean(height, na.rm = TRUE),
sd = sd(weight, na.rm = TRUE)
)

⑤ Visualize your data using box plots
아직 t-test를 한 건 아니고, group별로 데이터의 분포를 y축의 값에 제시할 수 있다. 색깔도 다르게 해서.
library("ggpubr")
ggboxplot(my_data, x = "group", y = "weight",
color = "group", palette = c("#00AFBB", "#E7B800"),
ylab = "Height", xlab = "Groups")

⑥ Preleminary test to check independent t-test assumptions
지금까지 했던 여러 tests는 t-test의 3가지 가정을 만족하는지 확인하기 위한 것들이었다. 
Assumption 1: Are the two samples independents?
Assumption 2: Are the data from each of the 2 groups follow a normal distribution?
Assumption 3: Do the two populations have the same variances?

⑦ Compute unpaired two-samples t-test
드디어 t-test를 한다. 데이터가 두 개의 서로 다른 numeric vectors로 저장되어 있으면 아래의 코드를 사용한다.
res = t.test(BD_height, TK_height, var.equal = TRUE)
데이터가 하나의 frame으로 저장되어 있으면 아래의 코드를 사용한다. 나의 경우에는 이 방법으로 했다.
res = t.test(height ~ group, data = my_data, var.equal = TRUE)

그 결과는 아래와 같았다.
Two Sample t-test
data: height by group
t = 0.808, df = 18, p-value = 0.4296
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval: -2.499273 5.623083
sample estimates: mean in group BD mean in group TK 174.7786 173.2167

p-value가 0.05보다 작지 않아 영가설을 채택해 두 집단의 평균이 다르지 않은 결과가 나왔다. 나의 경우에는 몸무게, 키라는 변수를 통제해야 했기에 다행스러운 결과이다!