4-6. Rmarkdown 카이제곱 test, Fisher's Exact test, trend test

Untitled

범주-범주의 교차표(contingency table) 만드는 2가지 방법

• table()
• descr패키지의 CrossTable()

# 데이터 준비
library(moonBook)
data(acs)


# 1. table()
table(acs$sex)

## 
## Female   Male 
##    287    570

table(acs$sex, acs$obesity)

##         
##           No Yes
##   Female 194  93
##   Male   373 197

# 2. descr패키지-CrossTable
# install.packages("descr")
library(descr)
a <- CrossTable(acs$sex, acs$obesity) # 빈도 뿐만 아니라 비율도 보여준다.
a

##    Cell Contents 
## |-------------------------|
## |                       N | 
## | Chi-square contribution | 
## |           N / Row Total | 
## |           N / Col Total | 
## |         N / Table Total | 
## |-------------------------|
## 
## ================================
##            acs$obesity
## acs$sex       No     Yes   Total
## --------------------------------
## Female       194      93     287
##            0.089   0.175        
##            0.676   0.324   0.335
##            0.342   0.321        
##            0.226   0.109        
## --------------------------------
## Male         373     197     570
##            0.045   0.088        
##            0.654   0.346   0.665
##            0.658   0.679        
##            0.435   0.230        
## --------------------------------
## Total        567     290     857
##            0.662   0.338        
## ================================

카이제곱 검정 - 관측값과 특정한 확률(기존의 기대값)와의 적합도 검정

• H0 : 관측값과 기대값이 동일하다 ( chisq.test ( 범주별 빈도, p = 기대값 ))
• 동전의 앞면/뒷면의 확률은 1/2인데, 관측결과 90/200, 110/200 이 나왔다면
• 그게 적합하게 나온 건지 확인하는 것

# 1. 관측값 
obs <- c(19, 41, 40) # 19 / 전체 수 로 확률계산이 내부에서 될 것이고, 기대값과 비교할 듯..

# 2. 특정한 확률(기대값)
null.probs <- c(2/10, 3/10, 5/10)

# 3. 적합도 검정
chisq.test(obs, p = null.probs)

## 
##  Chi-squared test for given probabilities
## 
## data:  obs
## X-squared = 6.0833, df = 2, p-value = 0.04776

#### acs데이터로 해보기
# 1. table()로 범주의 빈도를 관측값으로 가져오기
smk_type <- table(acs$smoking) # 범주1개를 table()을 이용하여, 범주별 빈도를 변수에 담음
smk_type

## 
## Ex-smoker     Never    Smoker 
##       204       332       321

# 2. 기대값(특정확률)은 직접 변수에 입력해주기 (집단의 수만큼)
smk_type_prob <- c(0.3, 0.35, 0.35)

# 3. 관측값(범주별 빈도) + 기대값으로 카이제곱 적합도 검정
chisq.test( smk_type, p = smk_type_prob)

## 
##  Chi-squared test for given probabilities
## 
## data:  smk_type
## X-squared = 15.869, df = 2, p-value = 0.0003582

카이제곱 검정(Chisq test) - 2개 범주 독립성 검정의 3가지 방법대

• chisq.test( 칼럼인덱싱, 칼럼인덱싱 )
• H0 : 두 범주 X, Y는 독립이다.
• 범주-범주의 교차표를 바탕으로 2개의 범주의 연관성(독립성)검정
• X-squared : 카이제곱 검정통계량 -> 교차표상 각 cell의 ( (x-E)^2 / E)를 구해서 다 더한 것

# 1. 교차표(table) 생성 후 검정
table( acs$sex, acs$obesity )

##         
##           No Yes
##   Female 194  93
##   Male   373 197

chisq.test(table( acs$sex, acs$obesity ))

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(acs$sex, acs$obesity)
## X-squared = 0.30627, df = 1, p-value = 0.58

# 2. 바로 검정
chisq.test(acs$sex, acs$obesity)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  acs$sex and acs$obesity
## X-squared = 0.30627, df = 1, p-value = 0.58

# 3. xtabs(교차표) & summary로 검정
xtabs( ~ sex + obesity, data = acs )

##         obesity
## sex       No Yes
##   Female 194  93
##   Male   373 197

summary( xtabs( ~ sex + obesity, data = acs ) )

## Call: xtabs(formula = ~sex + obesity, data = acs)
## Number of cases in table: 857 
## Number of factors: 2 
## Test for independence of all factors:
##  Chisq = 0.3968, df = 1, p-value = 0.5288

카이제곱 검정 다중비교 - 범주-범주이나, 독립변수의 집단이 3개 이상이어서 nC2로 -> 종속변수(2집단)에 독립인지 아닌지 본다.

• fifer 패키지의 chisq.post.hoc( table() )

# 깃허브로 설치
# library(devtools)
# install_github("cran/fifer")

# 다중검정 in 3x3
library(fifer)

## Loading required package: MASS

table( acs$smoking, acs$Dx )

##            
##             NSTEMI STEMI Unstable Angina
##   Ex-smoker     42    66              96
##   Never         50    97             185
##   Smoker        61   141             119

chisq.post.hoc( table( acs$smoking, acs$Dx ) )

## Adjusted p-values used the fdr method.

##             comparison  raw.p  adj.p
## 1  Ex-smoker vs. Never 0.1112 0.1112
## 2 Ex-smoker vs. Smoker 0.0233 0.0350
## 3     Never vs. Smoker 0.0000 0.0000

chisq.post.hoc( table( acs$Dx, acs$smoking ) )

## Adjusted p-values used the fdr method.

##                   comparison  raw.p  adj.p
## 1           NSTEMI vs. STEMI 0.2964 0.2964
## 2 NSTEMI vs. Unstable Angina 0.0114 0.0172
## 3  STEMI vs. Unstable Angina 0.0000 0.0000

Fisher’s exact Test

• 범주-범주의 교차표에서 각 cell의 Expected(기대값)이 5가 안되는 것이 25%이상(1/4이상)시 사용

# 교차표
table(mtcars$carb, mtcars$cyl)

##    
##     4 6 8
##   1 5 2 0
##   2 6 0 4
##   3 0 0 3
##   4 0 4 6
##   6 0 1 0
##   8 0 0 1

# 카이검정제곱 -> 카이제곱 approximation은 정확하지 않을수도 있습니다 -> Fisher로!
chisq.test(table(mtcars$carb, mtcars$cyl))

## Warning in chisq.test(table(mtcars$carb, mtcars$cyl)): Chi-squared
## approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  table(mtcars$carb, mtcars$cyl)
## X-squared = 24.389, df = 10, p-value = 0.006632

# fisher's exact test
fisher.test(mtcars$carb, mtcars$cyl)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  mtcars$carb and mtcars$cyl
## p-value = 0.0003345
## alternative hypothesis: two.sided

Trend test

• 순서를 가진 3집단이상의 범주(독립변수)에 대해 종속변수(2집단의 범주)가
• H0 : 비율이 동일하다(일정하다)
• H1 : 증가/감소 추세가 있다

# 1. trend-test 기본 
# 독립변수가 열 / 종속변수의 2개 집단이 행으로 나열된 교차표를 생각하자.
# 독립변수의 집단이 3개라고 가정 -> 3열
# 종속변수의 집단 2개 -> 2행 
prop.trend.test( c(15, 10, 30), # 종속변수 집단1 
                 c(35, 22, 44)) # 종속변수 집단2

## 
##  Chi-squared Test for Trend in Proportions
## 
## data:  c(15, 10, 30) out of c(35, 22, 44) ,
##  using scores: 1 2 3
## X-squared = 5.2588, df = 1, p-value = 0.02184

# 2. DescTools 패키지 - CochranArmitageTest()
# install.packages("DescTools")
library("DescTools")

# 데이터(matrix) 준비  2행 3열(독립변수 집단개수 = 열)
test <- matrix(c(2372,1859,2373,                 3065,2837,3065), byrow=TRUE, nrow=2)

# Desc( test )
# trend test 
CochranArmitageTest(test)

## 
##  Cochran-Armitage test for trend
## 
## data:  test
## Z = 0.011174, dim = 3, p-value = 0.9911
## alternative hypothesis: two.sided