3.1 One-sample t-test: 알려져 있는 평균과 비교

  1. 2006년 한국인의 1인 1일 평균 알코올 섭취량이 8.1g입니다. 2008년에 알코올 섭취량이 달라졌는지 조사해봅시다.
x = c(15.5,11.21,12.67,8.87,12.15,9.88,2.06,14.50,0,4.97) # 2008년 자료x를 생성
mean(x) # 평균
## [1] 9.181
sd(x) # 표준편차
## [1] 5.234965


  1. 검정 순서는 아래와 같습니다.
  • 자료가 정규분포를 하는 지 검정: shapiro.test()
  • 정규분포하는 경우: t.test()
  • 정규분포하지 않는 경우: wilcox.test()
shapiro.test(x)
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.92341, p-value = 0.3863


  1. 정규성 검정 결과 p-value가 0.05 이상이므로 정규분포 가정을 만족합니다. 따라서 t.test()를 시행합니다. 기본적으로 양측검정을 하며 유의수준 95%에서 검정을 합니다.
t.test(x, mu=8.1)
## 
##  One Sample t-test
## 
## data:  x
## t = 0.653, df = 9, p-value = 0.5301
## alternative hypothesis: true mean is not equal to 8.1
## 95 percent confidence interval:
##   5.436132 12.925868
## sample estimates:
## mean of x 
##     9.181


  1. 만일 단측검정을 하고 싶고 유의수준을 99%로 하고 싶으면 다음과 같이 입력합니다.
t.test(x, mu=8.1, conf.level = 0.99, alternative = "greater")
## 
##  One Sample t-test
## 
## data:  x
## t = 0.653, df = 9, p-value = 0.265
## alternative hypothesis: true mean is greater than 8.1
## 99 percent confidence interval:
##  4.510276      Inf
## sample estimates:
## mean of x 
##     9.181


3.2 Paired t-test & Wilcoxon Signed-rank test

  1. 쌍을 이룬 두 변수의 차이를 보는 검정 방법입니다.
  2. 검정 순서는 아래와 같습니다.
  • 자료가 정규분포를 하는 지 검정: shapiro.test()
  • 정규분포하는 경우: t.test(…, paired=TRUE)
  • 정규분포하지 않는 경우: wilcox.test()
data(sleep) # 내장데이터 불러오기, 데이터에 대한 설명은 ?sleep 입력
head(sleep)
##   extra group ID
## 1   0.7     1  1
## 2  -1.6     1  2
## 3  -0.2     1  3
## 4  -1.2     1  4
## 5  -0.1     1  5
## 6   3.4     1  6
str(sleep)
## 'data.frame':    20 obs. of  3 variables:
##  $ extra: num  0.7 -1.6 -0.2 -1.2 -0.1 3.4 3.7 0.8 0 2 ...
##  $ group: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ID   : Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...


  1. 먼저 정규성 검정을 시행합니다. p-value 0.05 이하로 정규분포를 따르지 않으므로 wilcoxn test를 시행합니다. 그러나 비모수 방법은 자료의 순서를 이용하는 데 동일한 값이 있어 순서를 정하는 데 문제가 있어 warning이 뜹니다. warning을 없애고 싶으면 exact=FALSE 옵션을 추가합니다.
with(sleep, shapiro.test(extra[group==2]-extra[group==1]))
## 
##  Shapiro-Wilk normality test
## 
## data:  extra[group == 2] - extra[group == 1]
## W = 0.82987, p-value = 0.03334
with(sleep, wilcox.test(extra[group==2]-extra[group==1], exact = FALSE))
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  extra[group == 2] - extra[group == 1]
## V = 45, p-value = 0.009091
## alternative hypothesis: true location is not equal to 0


  1. 정규분포를 따르는 경우에는 Student’s paired t-test를 시행하면 됩니다.
with(sleep, t.test(extra[group==1],extra[group==2],paired = TRUE))
## 
##  Paired t-test
## 
## data:  extra[group == 1] and extra[group == 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. 이 분석 결과를 그림으로 나타내면 다음과 같습니다.
  • 약을 1에서 2로 바꾸었을 때 수면시간이 늘어나는가?
sleep1 = with(sleep, extra[group==2]-extra[group==1])
summary(sleep1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    1.05    1.30    1.58    1.70    4.60


  1. 좀 복잡해보이는 구문이지만 익숙해지면 깔끔한 그림을 생성할 수 있습니다. 향후 그래프그리기는 다른 강의에서 설명하도록 하겠습니다.
stripchart(sleep1, method = 'stack', xlab = 'hours', main = 'Sleep prolongation (n = 10)')
boxplot(sleep1, horizontal = TRUE, add = TRUE, at = .6, pars = list(boxwex=0.5, staplewex=1)) # add =  그래프 결합여부, at = 그래프의 위치, boxwex, staplewex, outwex는 각각 상자의 크기, 수염의 길이, 이상치까지의 거리


3.3 두군에서의 평균비교

  1. 전체환자를 실험군과 대조군으로 나누어 실험군은 실험약물을 대조군은 위약을 투약한 후 효과를 비교합니다. 이 때, two-sample t-test를 사용합니다. 다음의 순서를 따라서 진행합니다.
  • 두 집단의 분산이 같은 지 검정: var.test()
  • 분산이 다르면 Welch의 t.test를 적용: t.test()
  • 분산이 같으면 pooled variance를 이용한 t.test()를 적용: t.test를 실시할 때 var.equal=TRUE
  1. 예를 들어 acs 데이터에서 성별에 따른 나이의 차이를 보고자 한다면 다음과 같이 그래프를 그릴 수 잇습니다.
require(moonBook)
data(acs)
# 성별에 따른 나이의 분산비교
aggregate(age~sex, data=acs, var)
##      sex      age
## 1 Female 115.1848
## 2   Male 126.0871
# 성별에 따른 나이의 박스플랏
boxplot(age~sex, data=acs, col=c("brown2","deepskyblue"))

# 원래는 정규분포가정도 확인해야합니다 (n 수가 많아서 그냥 넘어감)
# 여성의 분포는 정규분포 하지 않는 것으로 확인됩니다.
shapiro.test(acs$age[acs$sex == "Female"])
## 
##  Shapiro-Wilk normality test
## 
## data:  acs$age[acs$sex == "Female"]
## W = 0.96138, p-value = 6.34e-07
shapiro.test(acs$age[acs$sex == "Male"])
## 
##  Shapiro-Wilk normality test
## 
## data:  acs$age[acs$sex == "Male"]
## W = 0.99631, p-value = 0.2098
library(ggplot2)
ggplot(acs, aes(x=age, fill=sex)) +
  geom_histogram(binwidth=5, alpha=0.5)


  1. 분산이 같은지 검정합시다. var.test(반응변수~그룹변수, data=데이터)
var.test(age~sex, data=acs)
## 
##  F test to compare two variances
## 
## data:  age by sex
## F = 0.91353, num df = 286, denom df = 569, p-value = 0.3866
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.7495537 1.1209741
## sample estimates:
## ratio of variances 
##          0.9135342
  • 이 때 p-value는 0.3866으로 분산이 같다는 귀무가설을 기각하지 못합니다. 따라서 분산은 동일하므로 var.equal=TRUE로 주어야합니다.
  • 게다가 여성의 데이터 분포는 정규분포를 따르지 않기 때문에 원래는 비모수적 검정 방법인 Anasri-Bradley test를 통해 등분산을 검토해야합니다.
ansari.test(age~sex, acs)
## 
##  Ansari-Bradley test
## 
## data:  age by sex
## AB = 57724, p-value = 0.02265
## alternative hypothesis: true ratio of scales is not equal to 1
t.test(age~sex, data=acs, var.equal=TRUE)
## 
##  Two Sample t-test
## 
## data:  age by sex
## t = 10.071, df = 855, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  6.493487 9.637377
## sample estimates:
## mean in group Female   mean in group Male 
##             68.67596             60.61053
  • p-value가 매우 낮게 나왔습니다. 만일 분산이 다르다면 그냥 t.test()를 시행하면 됩니다.
t.test(age~sex, data=acs)
## 
##  Welch Two Sample t-test
## 
## data:  age by sex
## t = 10.222, df = 596.99, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  6.515847 9.615016
## sample estimates:
## mean in group Female   mean in group Male 
##             68.67596             60.61053
  • 시각적으로 확인해보면 더 확실히 알 수 있습니다. 2개의 평균이 서로 2 표준편차 안에 있는 지 보는 방법 입니다.
library(plyr)
ageSummary <- ddply(acs, "sex", summarize,
                    age.mean = mean(age), age.sd = sd(age),
                    Lower = age.mean - 2*age.sd /sqrt(NROW(age)),
                    Upper = age.mean + 2*age.sd /sqrt(NROW(age)))
ageSummary
##      sex age.mean   age.sd    Lower    Upper
## 1 Female 68.67596 10.73242 67.40893 69.94299
## 2   Male 60.61053 11.22885 59.66988 61.55118
ggplot(ageSummary, aes(x = age.mean, y = sex)) +
  geom_point() +
  geom_errorbarh(aes(xmin = Lower, xmax = Upper), height = 0.2)


3.4 세군 이상에서의 평균비교

  1. 세군 이상에서의 평균비교는 일원분산분석(one-way ANOVA)로 합니다. 일원분산분석은 설명변수가 연속형이 아닌 이산형인 회귀분석의 일종으로 생각하면 되는데, aov() 함수를 사용합니다. aov() 내부적으로 선형모형(lm)을 사용하여 적합시켜 분석합니다.
  2. 예를 들어 acs 데이터에서 진단명에 따른 LDL 농도를 비교하려면 다음과 같은 그래프를 그려봅니다.
aggregate(LDLC~Dx, data=acs, function(x) c(mean(x), sd(x)))
##                Dx    LDLC.1    LDLC.2
## 1          NSTEMI 126.09459  44.73373
## 2           STEMI 116.72449  39.45286
## 3 Unstable Angina 112.87724  40.38533
boxplot(LDLC~Dx, data=acs) # 박스플롯

moonBook::densityplot(LDLC~Dx, data=acs) # Densityplot


  1. aov()를 이용하여 모형적합을 실시합니다. aov(반응변수~그룹변수, data=데이터)
out=aov(LDLC~Dx, data=acs)
summary(out)
##              Df  Sum Sq Mean Sq F value  Pr(>F)   
## Dx            2   18765    9382   5.617 0.00377 **
## Residuals   830 1386305    1670                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 24 observations deleted due to missingness
  • 분산분석표에서 Dx의 F값의 p-value는 0.00377으로 세 군의 평균이 모두 같다는 귀무가설을 기각합니다. 즉, 세 군의 평균은 다릅니다(= 평균이 다른 그룹이 적어도 하나 이상 존재합니다).
  1. 다중비교 ANOVA는 세 군의 평균이 모두 같은가를 검정하지만 어느 군의 평균이 다른가를 답해주지는 않습니다. TukeyHSD() 함수는 다중비교를 통해 각각의 군에 대한 비교를 가능하게 해줍니다. 사용법은 매우 간단해서 aov() 함수의 결과를 넣어주면 됩니다.
TukeyHSD(out)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = LDLC ~ Dx, data = acs)
## 
## $Dx
##                              diff       lwr        upr     p adj
## STEMI-NSTEMI            -9.370105 -19.04130  0.3010954 0.0599378
## Unstable Angina-NSTEMI -13.217357 -22.47817 -3.9565482 0.0024227
## Unstable Angina-STEMI   -3.847252 -11.25450  3.5599980 0.4419434
  1. 결과를 보면 Unstable angina와 NSTEMI 환자 군 간에 차이가 있는 것을 알수 있습니다. 이 결과를 그래프로 나타내면 아래와 같습니다.
plot(TukeyHSD(out))

* y축 라벨이 축에 평행하게 되어 잇습니다. 축의 라벨은 par() 함수의 las 인수로 조절 합니다. 하지만 par 함수의 인수조절은 향후 모든 그래프에 영향을 주므로 반드시 원래의 default 값으로 회복시키는 것이 좋습니다.

par(las=1)
plot(TukeyHSD(out))

par(las=0) # Default
  • 축은 수평으로 바뀌었지만, y축 라벨이 길어 잘립니다. 이런 경우 mar 인수로 조절할 수 있습니다. mar 인수는 그래프의 plot region의 마진을 조절하는 함수로 c(bottom,left,top,right)의 순서로 정합니다. Default는 mar=c(5,4,4,2)로 되어있으므로 왼쪽 마진을 더 주면 축의 이름이 보입니다. mar 함수의 인수조절 또한 향후 모든 그래프에 영향을 주므로 반드시 원래의 default 값으로 회복시키는 것이 좋습니다.
par(las=1)
par(mar=c(5,12,4,2))
plot(TukeyHSD(out))

par(mar=c(5,4,4,2)) # Default
par(las=0) # Default


  1. 통계가설의 검정
  • ANOVA의 가정은 반응변수가 정규분포하며 각 군의 분산이 동일하다는 것이 가정입니다. 반응 변수의 정규분포를 확인하기 위해서 QQ plot을 그려볼 수 있습니다. car 패키지의 qqPlot() 함수를 이용해 봅시다.
require(car)
qqPlot(lm(LDLC~Dx, data=acs), main="qqPlot", simulate=TRUE, labels=TRUE)

## [1] 451 457
  • qqPlot의 결과를 보면 95% 신뢰구간이 점선으로 표시되어 있으며 정규분포의 가정을 만족시키지 못한다는 것을 알 수 있습니다.
  • shapiro.test()를 통해 잔차의 정규성을 검정할 수 있습니다. 단 표본의 개수가 5000개 이상인 경우에는 shapiro.test()는 에러가 나기때문에 Anderson-Darling test를 실시합니다. nortest 패키지의 ad.test()를 이용합니다.
shapiro.test(resid(out))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(out)
## W = 0.97137, p-value = 1.024e-11
nortest::ad.test(resid(out))
## 
##  Anderson-Darling normality test
## 
## data:  resid(out)
## A = 2.5798, p-value = 1.636e-06

  • shapiro.test()로 잔차의 정규성을 검정한 결과도 qqPlot과 마찬가지로 정규분포 가정을 만족시키지 않습니다.

  • 등분산 가정은 Bartlett’s test로 검정할 수 있습니다. bartlett.test(반응변수~그룹변수, data=데이터)

bartlett.test(LDLC~Dx, data=acs)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  LDLC by Dx
## Bartlett's K-squared = 3.3668, df = 2, p-value = 0.1857

  • Bartlett 검정결과도 p-value가 유의하지 않으므로 세 군의 분산은 유의하게 다르지 않다고 할 수 있습니다.

  • 또한 등분산 검정은 이상치(outlier)에 민감하므로 이상치가 있는지 검정해볼 필요가 있는데 이를 위해서는 car 패키지에 있는 outlierTest()를 사용할 수 있습니다.

outlierTest(out)
##     rstudent unadjusted p-value Bonferroni p
## 457 6.013561         2.7194e-09   2.2653e-06
## 451 4.998183         7.0625e-07   5.8830e-04
## 789 4.058303         5.4106e-05   4.5070e-02
  • 이상치 검정 결과, 세 개의 관측치가 이상치로 확인됩니다. 따라서 여러 검정을 종합해볼때, 이 모형은 ANOVA 모형의 가설을 만족시키지 못하므로 비모수통계인 Kruskal-Wallis Rank Sum test를 실시합니다.

3.5 비모수검정

  1. Two-sample t-test와 ANOVA는 반응변수가 정규분포하는 것을 가정으로 하고 시행합니다. 데이터가 정규분포 가정을 만족하지 못하는 경우 비모수 검정을 해야합니다. Two-sample t-test에 해당하는 비모수 검정은 Wilcoxon Rank-Sum test 이고 ANOVA에 해당하는 비모수 통계는 Kruskal-Wallis Rank Sum test 입니다.
  2. acs 데이터에서 성별에 따른 나이의 분포를 봅시다.
moonBook::densityplot(age~sex, data=acs)

* 반응변수가 정규분포하는지 검사하기 위해서 먼저 lm() 함수를 사용하여 선형모형에 적합시킨 후 잔차를 shapiro.test()에 넣어주어야 합니다.

out = lm(age~sex, data=acs)
shapiro.test(resid(out))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(out)
## W = 0.99343, p-value = 0.000808
  1. shapiro.test()의 결과 p-value가 매우 유의하게 나오므로 정규분포 가정을 만족하지 못합니다. 이런 경우 Wilcoxon Rank-Sum test를 실시합니다. 선형모형에 적합시킨 후 잔차를 shapiro.test()에 넣어주어야 합니다.
wilcox.test(age~sex, data=acs)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  age by sex
## W = 115271, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0


  1. 앞서 보았던 ANOVA에서 진단명에 따른 LDL의혈중농도 차이 모형은 ANOVA 모형의 가설을 만족시키지 못하므로 비모수 통계인 Kruskal-Wallis Rank Sum test를 실시합니다. kruskal.test(반응변수~그룹변수, data=데이터)
  • 이 때 반응변수는 연속형 변수여야하고 그룹변수는 두 개 이상의 수준을 갖는 범주주형 변수여야 합니다. 그룹변수가 두 개의 수준일 때 이 검정은 Mann-Whitney U test와 동일합니다. 이 때 주의할 점은 그룹변수가 factor가 아닌 문자열인 경우 에러가 나므로 반드시 factor로 바꿔주어야 합니다.
kruskal.test(LDLC~factor(Dx), data=acs)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  LDLC by factor(Dx)
## Kruskal-Wallis chi-squared = 10.733, df = 2, p-value = 0.004669


  1. 비모수 검정에서의 다중비교는 nparcomp 패키지에 있는 nparcomp() 함수를 이용할 수 있습니다.
require(nparcomp)
result=mctp(LDLC~Dx, data=acs)
## 
##  #----------------Nonparametric Multiple Comparisons for relative effects---------------# 
##  
##  - Alternative Hypothesis:  True differences of relative effects are not equal to 0 
##  - Estimation Method:  Global Pseudo Ranks 
##  - Type of Contrast : Tukey 
##  - Confidence Level: 95 % 
##  - Method = Fisher with 290 DF 
##  
##  #--------------------------------------------------------------------------------------# 
## 
summary(result)
## 
##  #----------------Nonparametric Multiple Comparisons for relative effects---------------# 
##  
##  - Alternative Hypothesis:  True differences of relative effects are not equal to 0 
##  - Estimation Method: Global Pseudo ranks 
##  - Type of Contrast : Tukey 
##  - Confidence Level: 95 % 
##  - Method = Fisher with 290 DF 
##  
##  #--------------------------------------------------------------------------------------# 
##  
##  #----Data Info-------------------------------------------------------------------------# 
##            Sample Size    Effect     Lower     Upper
## 1          NSTEMI  148 0.5466111 0.5189087 0.5740277
## 2           STEMI  294 0.4961078 0.4723244 0.5199088
## 3 Unstable Angina  391 0.4572811 0.4350779 0.4796554
## 
##  #----Contrast--------------------------------------------------------------------------# 
##        1  2 3
## 2 - 1 -1  1 0
## 3 - 1 -1  0 1
## 3 - 2  0 -1 1
## 
##  #----Analysis--------------------------------------------------------------------------# 
##       Estimator  Lower  Upper Statistic     p.Value
## 2 - 1    -0.051 -0.117  0.016    -1.785 0.174788691
## 3 - 1    -0.089 -0.152 -0.026    -3.310 0.002989234
## 3 - 2    -0.039 -0.092  0.014    -1.725 0.195779428
## 
##  #----Overall---------------------------------------------------------------------------# 
##   Quantile     p.Value
## 1 2.350988 0.002989234
## 
##  #--------------------------------------------------------------------------------------#
  • 다중비교 결과 ANOVA와 마찬가지로 Unstable angina와 NSTEMI 환자군 간에 차이가 있다는 것을 알 수 있습니다. mctp() 함수는 adjusted p-value를 구할 때 “Tukey”, “Dunnett”, “Sequen”, “Williams” 등의 방법을 쓸 수 있는데 특별히 지정하지 않으면 Tukey의 방법을 사용합니다.
  1. 단순하게 다중비교로 p값을 알고자 하는 경우에는 pairwise.wilcox.test를 사용할 수 있습니다. pairwise.wilcox.test(반응변수,그룹변수,p.adjusted.method)
pairwise.wilcox.test(acs$LDLC, acs$Dx, p.adjust.method = "bonferroni")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  acs$LDLC and acs$Dx 
## 
##                 NSTEMI STEMI 
## STEMI           0.2498 -     
## Unstable Angina 0.0041 0.2462
## 
## P value adjustment method: bonferroni


  1. 이상의 모든 과정은 mytable 함수를 쓰면 한번에 해결할 수 있습니다.
mytable(Dx~LDLC, data=acs, method=3)
## 
##                    Descriptive Statistics by 'Dx'                  
## ____________________________________________________________________ 
##             NSTEMI             STEMI         Unstable Angina     p  
##            (N=153)            (N=304)            (N=400)      
## -------------------------------------------------------------------- 
##  LDLC 123.0 [99.0;146.0] 115.0 [91.0;142.0] 111.0 [83.0;137.0] 0.005
## --------------------------------------------------------------------
  • method를 3으로 설정한 경우, 잔차의 정규성 검정을 통해 ANOVA 또는 K-W test를 실시합니다. 정규분포를 하는 경우 데이터는 평균과 표준편차로 제시되고 정규분포하지 않는 경우 중앙값과 사분위값으로 표시됩니다.