#7 7일차_카이제곱검정, F분포, ANOVA

01 단일집단의 모분산 추정/가설검정 (카이제곱검정)

 

02 F분포

 -01 f분포

 -02 모분산비 추정과 관련된 분포

 

03 ANOVA

---

01 단일집단의 모분산 추정/가설검정 (카이제곱검정)

 -01 추정

  1) 점추정

    E(S²) = σ²

     => 표본분산이 모분산 추정의 불편추정량(점추정량)

 

  2) 구간추정

   V = (n-1)S² / σ² ~ Chisq(n-1)

   P(X0.025 < V < X0.975) = 95% 신뢰구간

   [ (n-1)S² / X0.975, (n-1)S² / X0.025 ]

 

 -02 가설검정

   H0 : σ² = 1

   H1 : σ² ! = 1 (양측검정)

          σ²  > 1 (오른쪽검정)

          σ²  < 1 (왼쪽검정)

 

   V = (n-1)S² / σ² ~ Chisq(n-1)

---

[ 예제 - 모분산 추정 ]
v1 <- scan('초콜릿.txt')

# 1) 점추정
var(v1)  # 24.3351

# 2) 구간추정(95%)
qchisq(0.025, df = 50-1) # 31.55492
qchisq(1 - 0.025, df = 50-1) # 70.22241

[(n-1)S^2/70.22241, (n-1)S^2/31.55492]
c(49*24.3351/70.22241, 49*24.3351/31.55492) # [16.98062 37.78872]

---

[ 예제 - 모분산 가설검정 ]
초콜릿 공장의 초콜릿 무게에 대한 모분산이 24일 것이다라는 가설에 대한 검정
유의수준 : 5%

 

H0 : σ² = 24
H1 : σ² != 24

# 1) 검정통계량
V = (n-1) * S² / σ² ~ Chisq(49)
  = 49 * 24.3351 / 24
  = 49.68416

# 채택역 : [31.55492, 70.22241] => 귀무가설 채택
# 모분산이 24일 것이다라는 가설은 유의하게 받아들일 수 있다.

# ** 시각화
vx <- seq(0, 100, 0.01)
vy <- dchisq(vx, df = 49)
ld <- 31.55492
lu <- 70.22241

plot(vx, vy, type = 'l', lwd = 1.5,
     ylim = c(0, 0.05))
polygon(c(10,seq(10, ld, 0.01), ld),
        c(0, dchisq(seq(10,ld,0.01), df=49),0),
        col = 'red')

polygon(c(lu,seq(lu, 90, 0.01), 90),
        c(0, dchisq(seq(lu, 90, 0.01), df=49),0),
        col = 'red')

 

# 2) 유의확률
qchisq(0.5, df= 49) # 48.33497
V* = 49.68416 > 이므로 
p-value = P(V > 49.68416)
        = 1 - P(V > 49.68416)
        = 1 - pchisq(49.68416, df = 49) 
        = 0.4458623 >>> 0.025 (양측검정이므로 α/2) => 귀무가설 채택!

---

[ 연습문제 - 모분산 추정 ]
df1 <- read.table('weight.txt', header = T)
df1

# 우리나라 신생아의 몸무게 분산을 추정
# 점추정, 구간추정 모두

 

s2 <- var(df1$weight)
n <- nrow(df1)

# 2) 구간추정
ld <- qchisq(0.025, df = n-1)
lu <- qchisq(1-0.025, df = n-1)

c((n-1)*s2/lu, (n-1)*s2/ld) # 190333.7 447602.0

 

[ 연습문제 - 모분산 가설검정 ]
# 위 데이터를 토대로 모분산이 280000보다 작아졌을 것이라는 가설에 대한 검정

 

# 1) 검정통계량
V* = (n-1)*s2/280000
= 42.81852

qchisq(0.05, df =n-1) # 28.96472

기각역 : [0, 28.96472]
채택역 : [28.96472, inf]
# => 귀무가설 채택

# 2) 유의확률
P(V<42.81852)
= pchisq(42.81852, df = n-1)
= 0.5208931 >>>> 0.5 => 귀무가설 채택!

# 3) varTest 
  # 단일집단 : varsTest, 단일집단 : var.test
install.packages('EnvStats')
library(EnvStats)

result <-varTest(df1$weight, alternative = 'less', sigma.squared = 280000)
result$statistic
result$p.value

> result$statistic
Chi-Squared 
   42.81852 
> result$p.value
[1] 0.5208932

 

---

# 그래프 그리기

# 설정
df <- 43
alpha <- 0.05
stat <- 42.81852

# x축
x <- seq(0, 80, 0.01)
y <- dchisq(x, df)

# 임계값
crit <- qchisq(alpha, df)

# 그래프
par(mar = c(5, 4, 4, 2))
plot(x, y, type = "l", lwd = 2,
     xlab = "", ylab = "",
     main = paste("Chi-square(df =", df, ")"))

# ---------------------------
# 1 기각역 (빨간색)
x_rej <- seq(0, crit, 0.01)
y_rej <- dchisq(x_rej, df)

polygon(c(0, x_rej, crit),
        c(0, y_rej, 0),
        col = rgb(1, 0, 0, 0.3), border = NA)

# ---------------------------
# 2 검정통계량 이하 영역 (파란색)
x_stat <- seq(0, stat, 0.01)
y_stat <- dchisq(x_stat, df)

polygon(c(0, x_stat, stat),
        c(0, y_stat, 0),
        col = rgb(0, 0, 1, 0.3), border = NA)

# ---------------------------
# 기준선들
abline(v = crit, col = "red", lty = 2, lwd = 2)
abline(v = stat, col = "blue", lty = 2, lwd = 2)

# ---------------------------
# x축 표시
mtext(paste0("crit = ", round(crit, 2)),
      side = 1, line = 2, at = crit, col = "red")

mtext(paste0("stat = ", round(stat, 2)),
      side = 1, line = 3, at = stat, col = "blue")

---

02 F분포

 -01 f분포

- 두 카이제곱 통계량 비가 나타내는 분포

- 양수, 오른쪽으로 꼬리가 긴 (양의 왜도) 비대칭형

- 두 자유도가 커지면 좌우 대칭형의 정규분포에 근사

 

X ~ F(k₁, k₂)

 

# ** f 분포 시각화

vx <- seq(0, 5, 0.01)
vy1 <- df(vx, 3, 5)
vy2 <- df(vx, 10, 5)
vy3 <- df(vx, 50, 5)

vx <- seq(0, 5, 0.01)
vy_1 <- df(vx, 3,1)
vy_2 <- df(vx, 3, 5)
vy_3 <- df(vx, 3, 20)

par(mfrow=c(1,2))
plot(vx,vy1, type = 'l', col = 2, 
     ylim = c(0,1))
lines(vx, vy2, type = 'l', col = 3)
lines(vx, vy3, type = 'l', col = 4)

plot(vx,vy_1, type = 'l', col = 2, 
     ylim = c(0,1))
lines(vx, vy_2, type = 'l', col = 3)
lines(vx, vy_3, type = 'l', col = 4)

---

 -02 모분산비 추정과 관련된 분포
# 각각 자유도 n, m을 갖는 카이제곱분포 통계량을 X,Y라 할 때,
X ~ Chisq(n)
Y ~ Chisq(m)

F = (X/n)/(Y/m) ~ F(n,m)

모분산비에 대한 검정통계량
V1 ~ Chisq(n-1)
V2 ~ Chisq(m-1)

F = (V1/n-1)/(V2/m-1) ~ F(n-1,m-1)
  = S1²/S2²

# ** 증명
자유도 10 카이제곱분포
자유도 15 카이제곱분포

v1 <- rchisq(1000, df = 10)
v2 <- rchisq(1000, df = 15)

vf <- (v1/10)/(v2/15)

par(mfrow = c(1,1))

# 난수로부터 얻은 실제 분포
hist(vf, prob = T, ylim = c(0, 0.9))

vx <- seq(0, 7, 0.01)
vy <- df(vx, 10, 15)
lines(vx, vy, type = 'l', col = 'red')

---

 [ 예제 - 등분산성 가설검정 ]
 weight.txt 파일을 읽고 남/여 신생아의 몸무게 비교를 위한 등분산성 가정에 대한 가설검정 진행(유의수준 5%)

 

df1 <- read.table('weight.txt', header = T)

H0 : 남아몸무게분산= 여아몸무게분산 <=> 남아몸무게분산/여아몸무게분산 = 1
H1 : 남아몸무게분산= 여아몸무게분산 <=> 남아몸무게분산/여아몸무게분산 != 1 (양측검정)

# 1) 검정통계량
v1 <- df1[df1$gender==1, 2] # 남아 몸무게
v2 <- df1[df1$gender==2, 2] # 여아 몸무게

n <- length(v1) # 18
m <- length(v2) # 26

F* = 남아표본분산/여아표본분산
   = var(v1)/var(v2)
   = 2.177104

qf(0.025, n-1, m-1) # 0.3924002
qf(1 - 0.025, n-1, m-1) # 2.359863

# 채택역 : [0.3924002, 2.359863]
# => 귀무가설 채택!

# 2) 유의확률
# 양측검정에서는 p-value 정의가
# F* < 0.5(누적확률이 0.5인 임계값) => P(F < F*)
# F* > 0.5 => P(F > F*)

qf(0.5, n-1, n-1)

p-value = P(F > F*)
        = P(F > 2.177104)
        = 1 - P(F > 2.177104)
        = 1 - pf(2.177104, n-1, m-1) # 0.03763133
        = 0.03763133 > 0.025 (양측검정이므로) => 귀무가설 채택

# 3) var.test()
var.test(x,         # 첫번째 집단 데이터 
         y,         # 두번째 집단 데이터터
         ratio = 1, # 모분산비의 영가설 가정값
         alternative = 'two.sided', 'greater', 'less', # 대립가설
         conf.level = 0.95 ...) # 신뢰수준

var.test(v1, v2)

> var.test(v1, v2)

	F test to compare two variances

data:  v1 and v2
F = 2.1771, num df = 17, denom df = 25, p-value =
0.07526
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.9225552 5.5481739
sample estimates:
ratio of variances 
          2.177104 

---

03 ANOVA

 

 

[ 예제 - 세 집단의 모평균 비교 가설검정 ]
df1 <- read.csv('age.csv')

> head(df1)
  scale sex score age
1     1   2     8  56
2     1   2     5  33
3     1   2     7  49
4     1   1     4  53
5     1   1     5  74
6     1   2     3  42

---

[ 예제 - 세 집단의 모평균 비교 가설검정 ]
df1 <- read.csv('age.csv')

# 지역별 평균 나이의  차이가 있는지 가설검정
# age.csv 파일을 읽고 지역별(scale) 평균 연령의 차이가 발생하는지를 유의수준 5% 가설검정


y1 <- df1[df1$scale == 1 , 'age']
y2 <- df1[df1$scale == 2 , 'age']
y3 <- df1[df1$scale == 3 , 'age']

# 1) 가설
# H0 : 지역별 평균연령이 같다
# H1 : 지역별 평균연령이 같지 않다(차이가 있다) (오른쪽 검정)

# 2) 검정통계량
SST(총편차제곱의합) = SSE(오차제곱합) + SSt(처리제곱합)
자유도 n- 1                 n-k              k-1

F = MSt / MSE ~ F(k-1, n-k)

## (1) 오차제곱합 (각 집단의 편차제곱합의 총합)
sse1 <- sum((y1 - mean(y1))^2)
sse2 <- sum((y2 - mean(y2))^2)
sse3 <- sum((y3 - mean(y3))^2)
n <- nrow(df1)
k <- 3


sse <- sse1 + sse2 + sse3
mse <- sse/(n-k)

## (2) 처리제곱합 (집단간 편차제곱합의 총합)
n1 <- length(y1)
n2 <- length(y2)
n3 <- length(y3)

sst1 <- n1 * (mean(y1) - mean(df1$age))^2
sst2 <- n2 * (mean(y2) - mean(df1$age))^2
sst3 <- n3 * (mean(y3) - mean(df1$age))^2

sst <- sst1 + sst2 + sst3
mst <- sst / (k-1)

## (3) F통계량
F = MSt / MSE
  = mst/mse 
  = 0.3660281
 
qf(1-0.05, k-1, n-k) # 3.057621

# 기각역 : [3.057621, inf]
# 채택역 : [0, 3.057621]
# 귀무가설 채택

## (4) 유의확률
p-value = P(F > 0.3660281)
        = 1 - P(F < 0.3660281)
        = 1-pf(0.3660281, k-1, n-k)
        = 0.6941136 >>> 0.05 => 귀무가설 채택!
  
# 3) ANOVA

# * R에서는 요인을 반드시 factor 선언해야함!  
df1$scale <- factor(df1$scale)  

lm(formula = , # Y ~ 요인(X)
   data = )

m1 <- lm(age ~ scale, data = df1)
summary(aov(m1))

> summary(aov(m1))

> summary(aov(m1))
             Df Sum Sq Mean Sq F value Pr(>F)
scale         2    150   75.05   0.366  0.694
Residuals   147  30139  205.03   

  

	Df	Sum Sq	Mean Sq	F value	Pr(>F) (오른쪽방향이므로 큰 쪽 검정)
scale(요인)	2(k-1)	150(SSt)	75.05(MSt)	0.366(MEt/MES)	0.694
Residuals(오차)	147(n-k)	30139(SSE)	205.03(MES)