4-2. Rmarkdown T-test

R Notebook

등분산성, 정규성 검정 후 t-test 종류별로 해보기

# 데이터 준비
t_data <- data.frame(
  group = c(1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2),
  score = c(100,100,80,50,40,90,20,50,50,70,30,40,30,70,30,40,30,60,30,60),
  age = c(70,70,70,60,50,60,60,60,50,50,20,30,20,20,25,10,13,10,10,10)
)

head(t_data)

# 2집단(2범주)별 따로 데이터 인덱싱 -> 그래야 등분산성 var.test 및 t.test :  
# - var.test, t.test ( 집단1의 데이터, 집단2의 데이터)
# 범주별 특정범주만 가져오기 -> 행인덱싱 자리에 넣기
t_data1 <- t_data[ t_data$group == 1, ]
t_data2 <- t_data[ t_data$group == 2, ]


# t-test 전에 확인사항 3가지 중 2가지 - 1. 등분산성 2. 정규성 => 둘다 p-value 0.05보다 커야 만족

#### 1. 등분산성 확인후 독립표본 T-검정
####    H0 : 두 집단간 평균의 차이가 없다.
# 1_1. 등분산성 by F검정( = effect(그룹간 변동의 평균) / error(그룹내 변동의 평균) )
var.test( t_data1$score , t_data2$score ) # 0.1093 -> 2그룹간 등분산 H0 기각x -> 등분산성 만족

## 
##  F test to compare two variances
## 
## data:  t_data1$score and t_data2$score
## F = 3.0787, num df = 9, denom df = 9, p-value = 0.1093
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.7647065 12.3948431
## sample estimates:
## ratio of variances 
##           3.078704

# 1_2. 등분산성 확인후 t-test는 var.equal = T 옵션을 넣어야한다.
t.test( t_data1$score , t_data2$score , var.equal = TRUE) #  0.03199 : 평균차이 유의

## 
##  Two Sample t-test
## 
## data:  t_data1$score and t_data2$score
## t = 2.3247, df = 18, p-value = 0.03199
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   2.213727 43.786273
## sample estimates:
## mean of x mean of y 
##        65        42

# 1_3. 등분산성 옵션을 안주면 결과가 달라진다. t.test default가 등분산성 = FALSE => Welch t-test
t.test( t_data1$score , t_data2$score, var.equal = F) # 0.03531

## 
##  Welch Two Sample t-test
## 
## data:  t_data1$score and t_data2$score
## t = 2.3247, df = 14.289, p-value = 0.03531
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   1.819881 44.180119
## sample estimates:
## mean of x mean of y 
##        65        42

# 1_4. t.test()함수 다른 사용법
# : t.test( 종속변수(관측치) ~ 독립변수(범주칼럼명, 요인), dataframe, 유형(등분산성, 대응표본 등))
t.test( score ~ group, data = t_data, var.equal = T)

## 
##  Two Sample t-test
## 
## data:  score by group
## t = 2.3247, df = 18, p-value = 0.03199
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   2.213727 43.786273
## sample estimates:
## mean in group 1 mean in group 2 
##              65              42

#### 2. 정규성 검정의 3가지 방법 : 0.05보다 커서 H0 기각 안되어야 만족

# 2_1. shapiro-Wilk test 
shapiro.test(t_data1$score) # 0.4604 -> 정규성 만족

## 
##  Shapiro-Wilk normality test
## 
## data:  t_data1$score
## W = 0.93126, p-value = 0.4604

shapiro.test(t_data2$score) # 0.006914 -> H0기각 -> 정규성x -> Mann whitney or Wilcoxon rank-sum?

## 
##  Shapiro-Wilk normality test
## 
## data:  t_data2$score
## W = 0.77358, p-value = 0.006914

# 2_2. nortest 패키지의 ad.test
# install.packages("nortest")
library(nortest)
ad.test(t_data1$score) # 0.4458

## 
##  Anderson-Darling normality test
## 
## data:  t_data1$score
## A = 0.32911, p-value = 0.4458

ad.test(t_data2$score) # 0.006013

## 
##  Anderson-Darling normality test
## 
## data:  t_data2$score
## A = 1.0264, p-value = 0.006013

# 2_3. Kolmogrov-Smirnov test
# ks.test (test, pnorm, mean (test), sd (test))
# pnorm : 누적 정규 분포
?ks.test

## starting httpd help server ... done

ks.test(t_data1$score, "pnorm", mean(t_data1$score), sd(t_data1$score)) # 0.7726

## Warning in ks.test(t_data1$score, "pnorm", mean(t_data1$score),
## sd(t_data1$score)): ties should not be present for the Kolmogorov-Smirnov
## test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  t_data1$score
## D = 0.20947, p-value = 0.7726
## alternative hypothesis: two-sided

ks.test(t_data2$score, "pnorm", mean(t_data2$score), sd(t_data2$score)) # 0.41

## Warning in ks.test(t_data2$score, "pnorm", mean(t_data2$score),
## sd(t_data2$score)): ties should not be present for the Kolmogorov-Smirnov
## test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  t_data2$score
## D = 0.28071, p-value = 0.41
## alternative hypothesis: two-sided

ks.test(t_data1$score, pnorm)

## Warning in ks.test(t_data1$score, pnorm): ties should not be present for
## the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  t_data1$score
## D = 1, p-value = 4.122e-09
## alternative hypothesis: two-sided

ks.test(t_data2$score, pnorm)

## Warning in ks.test(t_data2$score, pnorm): ties should not be present for
## the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  t_data2$score
## D = 1, p-value = 4.122e-09
## alternative hypothesis: two-sided

ks.test(t_data2$score, pnorm)[2]$p.value

## Warning in ks.test(t_data2$score, pnorm): ties should not be present for
## the Kolmogorov-Smirnov test

## [1] 4.122307e-09

#### 3. 대응표본 T-검정 (실험전/후를 범주로 주어진 상태라 가정)
# 2_1. 대응표본은 정규성 검사가 필수이다.
# 2_1_1. 샤피로 테스트
shapiro.test(t_data1$score) # 0.4604 => H0 기각x => 정규성 만족

## 
##  Shapiro-Wilk normality test
## 
## data:  t_data1$score
## W = 0.93126, p-value = 0.4604

shapiro.test(t_data2$score) # 0.006914 => H0 기각 => 정규성 x => 비모수적 방법인 

## 
##  Shapiro-Wilk normality test
## 
## data:  t_data2$score
## W = 0.77358, p-value = 0.006914

# 2_2. 대응표본일 경우 옵션 paired = T
t.test( t_data1$score, t_data2$score, paired = TRUE) # 0.05106

## 
##  Paired t-test
## 
## data:  t_data1$score and t_data2$score
## t = 2.2493, df = 9, p-value = 0.05106
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1311024 46.1311024
## sample estimates:
## mean of the differences 
##                      23

# 2_3. t.test 다른 사용법
t.test( score ~ group, data = t_data, paired = TRUE) # 0.05106

## 
##  Paired t-test
## 
## data:  score by group
## t = 2.2493, df = 9, p-value = 0.05106
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1311024 46.1311024
## sample estimates:
## mean of the differences 
##                      23

#### 3. 보정된 t검정 : 요인1에 대한 t.test가 끝난 뒤 -> 해당 요인(범주)이외에 다른 변수까지 같이 고려한다.
#### - method로서 선형회귀 -> 2way-ANOVA를 2집단에 적용해서, 분산분석으로 -> 2집단의 평균차이의 유의성을 검증(p-value값이 t-test와 똑같이 나온다.
#### - 그러기 위해선 lm()의 선형회귀로 기울기+절편 결과값을 summary()에 넣어, ANOVA테이블 형태를 얻어야한다.

# 3_1. 1개 요인(t-test와 같은 모델)에 대해 선형회귀->ANOVA 테이블에서 독립표본 t-test와 같은 결과값 얻기
#  선형회귀한 결과값을 분산분석(2way-요인2개 가능)을 2집단에 적용해 -> 평균 차이를 알아본다
mod1 <- lm( score ~ group , data = t_data) 
anova(mod1) # 0.03199 : 아노바 테이블 상 독립표본 t-test와 같은 p 결과값이 나온다.

t.test( t_data1$score , t_data2$score , var.equal = TRUE) # 0.03199

## 
##  Two Sample t-test
## 
## data:  t_data1$score and t_data2$score
## t = 2.3247, df = 18, p-value = 0.03199
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   2.213727 43.786273
## sample estimates:
## mean of x mean of y 
##        65        42

# 3_2. t-test에서 할 수 없는, 새로운 요인(숫자형)을 추가해서 lm->anova 테이블로 p값을 보자.
#      숫자형 요인은 단독으로는 t-test를 할 수 없다(집단 못나눔). 보정시에만 들어감
#      보정변수를 + 앞에 주는 것과 뒤에 주는 것과 결과가 다르다! ㅠㅠ;
#      예제는 추가한 요인 age를 group 앞에 주었다. 그대로 해보자.
mod2 <- lm( score ~ age + group, data = t_data)
anova(mod2) # 0.58633 : age를 앞쪽 요인으로 추가해줬더니, group에 의한 차이는 유의미하지않았다 -> 아! 기존 t-test에서 유의미한 것은 group 때문이 아니라 age땜에 난 것임을 알 수 있다.

# 3_3. 다른 방식 : aov() -> summary()
#      기존 방식 : 기존 lm() -> anova()
summary( aov( score ~ group , data = t_data) ) # 0.032

##             Df Sum Sq Mean Sq F value Pr(>F)  
## group        1   2645  2645.0   5.404  0.032 *
## Residuals   18   8810   489.4                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary( aov( score ~ age + group , data = t_data) ) # 0.586

##             Df Sum Sq Mean Sq F value Pr(>F)  
## age          1   3405    3405   7.322  0.015 *
## group        1    143     143   0.308  0.586  
## Residuals   17   7906     465                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 3_4. my) 만약 age + group 의 자리를 바꾸면 결과값이 달라진다.
summary( aov( score ~ age + group , data = t_data) ) # 0.586

##             Df Sum Sq Mean Sq F value Pr(>F)  
## age          1   3405    3405   7.322  0.015 *
## group        1    143     143   0.308  0.586  
## Residuals   17   7906     465                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary( aov( score ~ group + age , data = t_data) ) # 0.029 *

##             Df Sum Sq Mean Sq F value Pr(>F)  
## group        1   2645  2645.0   5.687  0.029 *
## age          1    904   903.6   1.943  0.181  
## Residuals   17   7906   465.1                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

연습해보기

• COPD 데이터를 가져와서, lungC(폐암)에 대한 모든 변수(요인, id칼럼만 제외! )들 종속변수로을 자동검정해보기
• mytable( 요인 ~ 종속변수 시리즈 )
• lungC(폐암) 여부(이력, 2집단으로 이루어진 범주)이 BMI(종속변수)에 주는 영향에 대한 t-test(2집단의 평균차이 유의성)확인해보기

# 데이터 준비
getwd()

## [1] "C:/Users/is2js/R_da/myR"

data <- read.csv("week4/4주차_코드/COPD_Lung_Cancer.csv")
head(data)

# 폐암칼럼에 대한 다른칼럼들을 모두 검정해보기
# - 칼럼명이 적히는 곳에 . : 모든 칼럼을 의미 ,  .-칼럼명 : 해당칼럼 제외 모든 칼럼을 의미!!***
library(moonBook)
mytable( lungC ~ .-PERSON_ID, data = data)

## 
##           Descriptive Statistics by 'lungC'          
## ______________________________________________________ 
##                             0            1         p  
##                          (N=870)       (N=74)   
## ------------------------------------------------------ 
##  SEX                                             0.002
##    - 1                 350 (40.2%)   44 (59.5%)       
##    - 2                 520 (59.8%)   30 (40.5%)       
##  DTH                                             1.000
##    - 0                 864 (99.3%)  74 (100.0%)       
##    - 1                  6 ( 0.7%)     0 ( 0.0%)       
##  CTRB_PT_TYPE_CD        5.7 ±  3.4   5.3 ±  3.5  0.307
##  BMI                   22.7 ±  3.3  22.8 ±  2.9  0.710
##  BP_LWST               81.0 ± 12.1  82.4 ± 11.5  0.326
##  BLDS                  102.4 ± 33.2 97.9 ± 26.6  0.178
##  TOT_CHOLE             200.3 ± 39.9 193.1 ± 35.6 0.132
##  HMG                   13.1 ±  1.4  13.2 ±  1.5  0.503
##  FMLY_CANCER_PATIEN_YN                           0.805
##    - 1                 752 (95.4%)   63 (96.9%)       
##    - 2                  36 ( 4.6%)   2 ( 3.1%)        
##  SMK_STAT_TYPE                                   0.654
##    - 1                 670 (80.7%)   54 (77.1%)       
##    - 2                  69 ( 8.3%)   8 (11.4%)        
##    - 3                  91 (11.0%)   8 (11.4%)        
##  DRNK_HABIT                                      0.495
##    - 1                 673 (80.0%)   52 (74.3%)       
##    - 2                  47 ( 5.6%)   6 ( 8.6%)        
##    - 3                  44 ( 5.2%)   6 ( 8.6%)        
##    - 4                  27 ( 3.2%)   1 ( 1.4%)        
##    - 5                  50 ( 5.9%)   5 ( 7.1%)        
##  EXERCI_FREQ                                     0.051
##    - 1                 674 (80.4%)   48 (68.6%)       
##    - 2                  68 ( 8.1%)   10 (14.3%)       
##    - 3                  27 ( 3.2%)   3 ( 4.3%)        
##    - 4                  9 ( 1.1%)    3 ( 4.3%)        
##    - 5                  60 ( 7.2%)   6 ( 8.6%)        
##  copd                                            0.000
##    - 0                 534 (61.4%)   29 (39.2%)       
##    - 1                 336 (38.6%)   45 (60.8%)       
##  Asthma                                          0.006
##    - 0                 422 (48.5%)   23 (31.1%)       
##    - 1                 448 (51.5%)   51 (68.9%)       
##  DM                                              0.012
##    - 0                 335 (38.5%)   40 (54.1%)       
##    - 1                 535 (61.5%)   34 (45.9%)       
##  Tub                                             0.012
##    - 0                 807 (92.8%)   62 (83.8%)       
##    - 1                  63 ( 7.2%)   12 (16.2%)       
##  before_op_score       40.1 ±  6.0  40.1 ±  5.8  0.914
##  after_op_score        50.1 ±  6.0  50.0 ±  6.2  0.907
## ------------------------------------------------------

# t-test 전 준비하기(범주(이력, 여부)별로 데이터 가르기)
data_lung0 <- data[ data$lungC == 0, ]
data_lung1 <- data[ data$lungC == 1, ]

# 1. lungC 여부에 따른 BMI의 등분산성 검사
var.test(data_lung0$BMI, data_lung1$BMI) # 0.1736 -> 등분산성 만족 -> var.equal = T 넣을 것

## 
##  F test to compare two variances
## 
## data:  data_lung0$BMI and data_lung1$BMI
## F = 1.2864, num df = 869, denom df = 73, p-value = 0.1736
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.8933904 1.7632502
## sample estimates:
## ratio of variances 
##           1.286364

# 2. 각각 정규성 검사
shapiro.test( data_lung0$BMI ) # 0.297

## 
##  Shapiro-Wilk normality test
## 
## data:  data_lung0$BMI
## W = 0.99777, p-value = 0.297

shapiro.test( data_lung1$BMI ) # 0.1953 -> 정규성 만족

## 
##  Shapiro-Wilk normality test
## 
## data:  data_lung1$BMI
## W = 0.97698, p-value = 0.1953

# 3. 독립표본 t-test
t.test( data_lung0$BMI, data_lung1$BMI, var.equal = T) # 0.7096 

## 
##  Two Sample t-test
## 
## data:  data_lung0$BMI and data_lung1$BMI
## t = -0.37253, df = 942, p-value = 0.7096
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.9266823  0.6309930
## sample estimates:
## mean of x mean of y 
##  22.70189  22.84973

t.test( BMI ~ lungC , data = data, var.equal = T ) # 0.7096 ( 종속변수 ~ 요인(범주, 이력, 여부))

## 
##  Two Sample t-test
## 
## data:  BMI by lungC
## t = -0.37253, df = 942, p-value = 0.7096
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.9266823  0.6309930
## sample estimates:
## mean in group 0 mean in group 1 
##        22.70189        22.84973

# 4. 대응표본 t-test 
# t.test ( data_lung0$BMI, data_lung1$BMI, paired = T ) # 데이터가 갯수가 달라서, 실험 전/후의 데이터 가 아님



# 5. 보정된 t-test
# 5_1. 요인1개
summary( aov( BMI ~ lungC , data = data)) # 0.71 

##              Df Sum Sq Mean Sq F value Pr(>F)
## lungC         1      1   1.491   0.139   0.71
## Residuals   942  10118  10.741

anova( lm( BMI ~ lungC, data = data)) # 0.7096

# 5_2. 요인 1개 더 추가(숫자형)
summary( aov( BMI ~ TOT_CHOLE + lungC , data = data)) # TOT_CHOLE 0.00121, lung_C 0.60168

##              Df Sum Sq Mean Sq F value  Pr(>F)   
## TOT_CHOLE     1    112  111.62  10.535 0.00121 **
## lungC         1      3    2.89   0.273 0.60168   
## Residuals   931   9864   10.59                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 10 observations deleted due to missingness

summary( aov( BMI ~  lungC + TOT_CHOLE , data = data)) #lung_C 0.71773, TOT_CHOLE 0.00112

##              Df Sum Sq Mean Sq F value  Pr(>F)   
## lungC         1      1    1.39   0.131 0.71773   
## TOT_CHOLE     1    113  113.12  10.677 0.00112 **
## Residuals   931   9864   10.59                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 10 observations deleted due to missingness