4-2. Rmarkdown T-test
2019. 2. 20. 15:47
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
'한의대 생활 > └ 통계에 대한 나의 정리' 카테고리의 다른 글
4-4. Rmarkdown 칼럼명변경 +range() + 사용자정의함수로 scaling (0) | 2019.02.21 |
---|---|
4-3. Rmarkdown ANOVA 와 interactionplot (1) | 2019.02.21 |
4-2. 3집단 이상의 (숫자형)분석 ANOVA (0) | 2019.02.19 |
4-1. 2개 집단의 평균 비교 - t-test (4) | 2019.02.19 |
4. Rmarkdown - 상관계수와 산점도 matrix (0) | 2019.02.19 |