Big Bro's Studying Archive

R로 접하는 통계 : T - test 본문

Big Data展

R로 접하는 통계 : T - test

빅브로오 2020. 5. 3. 14:34

Author : Yoon Baek

들어가기 전에

이번 글을 통해 T검정을 정규분포와 유사하게 사용할 수 있음을 입증해보고,
이렇게 입증한 T검정을 통해 모표준편차를 모르는 표본의 T검정에 직접 활용해보는 연습을 해보고자 한다.
궁금한 점은 댓글을 통해 달아주시면 성심성의껏 답변해드리겠습니다.

정규분포와 t분포와의 관계

모집단의 모수 mu, sigma가 있을 경우(모집단의 분포는 상관없음),
모집단으로부터 추출된 표본평균(xbar)는 다음의 분포를 따른다.

xbar ~ N(mu,sigma^2/n)
xbar - mu / (sigma/sqrt(n)) ~ N(0,1) # 표준화된 확률변수 Z의 분포

따라서 해당 이론을 통해 모집단의 모평균(mu)를 추정, 가설검정할 수 있음

따라서 모집단의 모분산(혹은 모표준편차)를 아는 경우
표준화된 확률변수 Z의 분포를 통해 가설검정하는 경우를 Z- test라고 부름

이떄, 모분산(혹은 모표준편차)를 모르는 경우에
모표준편차(sigma ) 대신 표본표준편차(s)로 대체할 수 있고
이를 자유도가 n-1을 갖는 t분포를 따른다고 가정한다

(xbar - mu) / (s/sqrt(n)) ~ T(n-1)

1. 정규분포와 T분포와의 관계 시각화

평균이 10, 분산이 4(표준편차 2)인 정규분포를 따르는 난수를 100번 추출,
표본평균이 갖는 분포는 다음과 같은 정규분포를 따른다 알려져 있다.

xbar ~ N(10,4/n)
(xbar - mu) / (sigma/sqrt(n)) ~ N(0,1) # sigma를 알 경우
(xbar - mu) / (s/sqrt(n)) ~ T(n-1)     # sigma를 모를 경우

v_xbar <- c(); v_sd <- c()

for (i in 1:1000) {
  v_data <- rnorm(n = 100, mean = 10, sd = 2)
  v_xbar[i] = mean(v_data)
  v_sd[i] = sd(v_data)
}

반복문의 매 시행마다 평균 = 10, 표준편차는 2를 갖는 난수 100개를 추출하고
그 추출한 집단의 평균과 표준편차를 v_xbar, v_sd에 쌓아 주는 반복문이다.


위와 같은 v_data의 mean과 sd를 1000개 구해 쌓는다.

실제 분포 시각화

dev.new()
par(mfrow = c(1,2))
hist(v_xbar, prob = T)

sigma를 알고 있을 경우의 이론적 분포 시각화

v_x <- seq(min(v_xbar), max(v_xbar), 0.01)
v_y <- dnorm(v_x, mean = 10, sd = 2/sqrt(100))
lines(v_x,v_y,type = 'l', col = 'red', lwd = '2')

dnorm 함수는 벡터연산이 가능하기 때문에 v_xbar의 최솟값과 최댓값 사이의
0.01 단위의 시퀀스 값들을 촘촘히 넣어주면 v_y값에 저장이 된다.
Python과는 다르게 벡터연산이 편리하게 이뤄진다는 점은 R의 장점이다.

2) sigma를 모르고 있을 경우의 이론적 분포 시각화

(xbar - mu) / (s/sqrt(n)) ~ T(n-1)
mu <- 10
n  <- 100
x_t <- (v_xbar - mu) / (v_sd / sqrt(n))

1000개의 xbar와 sd 값을 바탕으로 표준화

2-1) 표준화된 확률변수의 실제분포(sigma 대신 s로 대체)

hist(x_t, prob = T)

총 1000개의 x_t의 값을 히스토그램으로 표현

2-2) 분포의 이론적 분포

v_x2 <- seq(min(x_t), max(x_t), 0.01)
v_y2 <- dt(v_x2, n-1)
lines(v_x2, v_y2, type = 'l', col = 'blue', lwd = '2.5')


오른쪽 그래프를 볼 때 x_t의 실제 분포(히스토그램막대)와
이론적 분포(파란색 선)이 거의 유사함을 확인할 수 있고
이는 왼쪽의 정규분포 함수 그래프의 모양과도 유사하다.

따라서 t-test를 sigma를 모를 때의 유의미한 검정 방식으로 활용할 수 있음을 알 수 있다.

[ 연습문제 : t-test ]

랜덤하게 샘플링한 초콜릿 50개의 무게의 측정 데이터가 다음과 같을 때
모평균이 200과 다르다고 할 수 있는지 유의수준 5%에서 검정하시오.

t1 <- read.table('초콜릿.txt')
v1 <- unlist(t1)


mu <- 200
xbar <- mean(v1)
s <- sd(v1)
n <- length(v1)

1. 가설 수립

H0 : mu = 200
H0 : mu != 200

2. t통계량의 기각역과 채택역

ld<- qt(0.05/2, df = n-1)     # -2.009575
lu<- qt(1-0.05/2, df = n-1)   #  2.009575

여기서 df는 자유도를 의미한다. (degree of freedom)
채택역 : [-2.009575, 2.009575]

3. t분포의 검정 통계량

(xbar - mu) / (s/sqrt(n)) ~ T(n-1)

vt <- (xbar - mu) / (s/sqrt(n))   # -0.7740378

결론 :
검정통계량 -0.7740378은 자유도가 49인 t분포의 채택역인
채택역 : [-2.009575, 2.009575] 구간 안에 포함, H0을 채택한다.

4. 유의확률

P(T > |t*|) = P(T < t*)
            = P(T < -0.7740378)
            = pt(-0.7740378, df = n-1)
            = 0.2213136 >>>> 0.025

=> 유의 확률이 유의수준(양측검정)보다 크므로 H0에 유의하다고 볼 수 있다.
위 과정을 간단하게 실행해주는 함수가 t.test 함수이다.
t.test 함수에서 p-value는 0.05 기준으로 나오게 출력해준다.
따라서 0.2213을 두배 곱한 0.4426의 pvalue를 출력해준다.

t.test(v1, mu = 200)
# data:  v1
# t = -0.77404, df = 49, p-value = 0.4426
# alternative hypothesis: true mean is not equal to 200
# 95 percent confidence interval:
#   198.058 200.862
# sample estimates:
#   mean of x 
# 199.46 

[ 연습문제 ]

여아신생아.txt파일을 읽고

data = read.table('여아신생아.txt')
data = unlist(data)
names(data) <- NULL

신생아 몸무게가 2800보다 클 것이다라는 가설에 대한 가설 검정 수행

H0 : mu <= 2800
H1 : mu > 2800

mu = 2800
xbar = mean(data)
s = sd(data)

기각역

결과값 : [1.739607, limit]
표준화한 값이 위 범위에 해당하면 H0을 기각, H1을 채택

lu = qt(0.95, df = length(data) - 1)        # 1.739607
v_t = (xbar - mu) / (s/sqrt(length(data)))  # 2.233188

유의 수준 구하기

1-pt(v_t, df = length(data-1))              # 0.01963422

0.0196의 유의수준을 가지므로 H1을 채택

t.test 함수로 위 검정이 맞는지 확인

t.test(data, mu = 2800, alternative = 'greater')

'greater'옵션을 통해 우측 검정 t test를 수행할 수 있다.
'less'는 좌측 검정, 'two-sided'는 양측 검정.

Comments