axnuo

[R] 등분산 검정, ANOVA 검정 본문

카테고리 없음

[R] 등분산 검정, ANOVA 검정

axnuo 2023. 10. 19. 00:54

가설 검정 실습

gr2=b2[b2$F_BR==8,] ##매일 식사하는 군
mean(gr2$HT) ##키 평균
sd(gr2$HT)##표준편차

gr1=b2[(b2$F_BR<=7 & b2$F_BR>=1),]
gr1_1=b2[(b2$F_BR!=8),]
dim(gr1) ##매일 식사하지 않는 군

mean(gr1$HT)
sd(gr1$HT)

b2$always <-ifelse(b2$F_BR != 8, 1, 0)
View(b2)
b2_always <- b2[b2$always ==1, ]
b2_not_always <- b2[b2$always ==0, ]
mean(b2_always$HT)

  1. normal distribution of continuous var→ 정규분포를 따르는가? 따르지 않는가?
  2. → H0 : normal dist O, H1: X
  3. equality of variance→ H0 : 시그마1 제곱 = 시그마2 제곱, H1 : 서로 같지 않다.⇒ 두 개의 그룹의 평균이 다르다는 것을 증명하고 싶음. 만약 variance가 같다고 나오면 검정통계랑을 구하는 공식이 달라진다.
  4. 궁극적으로 하려는 것은 평균을 비교하는 것인데 왜 굳이 variance 를 비교하는가?
  5. 같은 평균을 갖고 있다고 해도 퍼진 정도가 다르다면 시그마 제곱(분산)이 다른 거!
  6. compare with true means of group(level)
    aa1 <- aa[aa$group != "trt2", ]
    shapiro.test(aa1$weight) #p-value 구할 수 있음
    hist(aa1$weight)
    qqplot(aa1$group, aa1$weight)
    qqnorm(aa1$weight)
    
    var.test(weight~group, data=aa1, conf.level=0.95)
    #x = independent variances = group, y = dep = weight
    #var.test는 등분산검정 == F검정
    
    • var.test 실행 결과⇒ p-value> 0.05 ⇒ 분산이 같다
    • F test to compare two variances data: weight by group F = 0.53974, num df = 9, denom df = 9, p-value = 0.3719 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.1340645 2.1730025 sample estimates: ratio of variances 0.5397431
  7. 두 그룹의 평균 비교하기

 데이터 분석

실습 데이터 → 데이터가 상당히 많기 때문에 중심 극한 정리가 나타남. 굳이 검정을 하지 않아도 정규분포를 따른다는 가정 하에 진행

b1$all <-ifelse(b1$F_BR ==8, 1, 0)
var.test(HT~all, data = b1, conf.level=0.95)
t.test(HT ~ all, data=b1, var.equal=FALSE, conf.level=0.95)

b1에 매일 아침밥을 먹는 군을 1로 아닌 사람을 0으로 설정해서 한 열을 만듦

그 다음에 p검정을 한당!!^_^

  • var.test 결과→ 분산이 다르다. p value가 0.05보다 작잖아 ⇒ var.equal=FALSE로 옵션을 줘야함.
  • F test to compare two variances data: HT by all F = 0.95874, num df = 36911, denom df = 13536, p-value = 0.00291 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.9323159 0.9857253 sample estimates: ratio of variances 0.958736
  • t.test 결과
  • Welch Two Sample t-test data: HT by all t = -6.0348, df = 23656, p-value = 1.615e-09 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.6748571 -0.3439559 sample estimates: mean in group 0 mean in group 1 166.1259 166.6353

결론 : 두 평균의 값은 같지 않다.

아무튼 데이터가 달라지면 분석을 할 때 적용하는 스태티스틱이 달라진다

 

 b2 = b[,c("PA_TOT","WT")]
 View(b2)
> PA_TOT_gr <-ifelse(b2$PA_TOT>=3, 1, 0)
> b2$PA_TOT_gr <-ifelse(b2$PA_TOT>=3, 1, 0)
> var.test(WT~PA_TOT_gr, data = b2, conf.level=0.95)

	F test to compare two variances

data:  WT by PA_TOT_gr
F = 0.89163, num df = 23532, denom df = 26921, p-value < 2.2e-16
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.8698581 0.9139752
sample estimates:
ratio of variances
         0.8916346

> t.test(WT ~ PA_TOT_gr, data=b2, var.equal=FALSE, conf.level=0.95)

	Welch Two Sample t-test

data:  WT by PA_TOT_gr
t = -27.29, df = 50154, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.500455 -3.031330
sample estimates:
mean in group 0 mean in group 1
       57.79467        61.06056

> hist(b2$WT)
> hist(b2$PA_TOT)
> summary(b2$PA_TOT_gr)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 0.0000  0.0000  1.0000  0.5329  1.0000  1.0000
> over3 <- b2[b2$PA_TOT_gr==1,] #
> less3<- b2[b2$PA_TOT_gr==0,]#얘네는 n 알아보려고 변수 만들음

 

aa = PlantGrowth
table(aa$group) #각 요소의 n 파악 가능
library(descr)
freq(aa$group) # 각 요소의 % 파악 가능
summary(aa$weight)#값의 분포 확인 na값 확인 가능
boxplot(aa$weight)#최댓값 최솟값 확인 가능, radian도 나오고 이상치가 발생하면 삭제
hist(aa$weight)#전체적인 값 확인
shapiro.test(aa$weight)#가설검정 확인 p-value>0.05라서 기각할 수 없으니 
#continuous variable로 사용할 수 있다. 정규분포를 따르는지 확인 == 정규성검정
qqnorm(aa$weight)#만약 0.05보다 작으면 넘어가는 게 아니라 qqnorm을 그림 일직선에 가까워야함!
qqline(aa$weight)#이걸로 일직선 확인
aggregate(weight~group, data=aa, mean)
#weight= 영향을 받는 애들, group = 영향을 주는 애들 이거로 하면 평균이 나옴
aggregate(weight~group, data=aa, sd)#sd가 나옴

n m+- sd p-value

ctrl 10 5.0 +- 0.58 0.0159
trt1 10 4.7 +- 0.79  
trt2 10 5.5 +- 0.44  

step 1. normal을 만족하는지 판단

step2. 분산이 같은지 판단

등분산 검정을 하는 이유 : statistic이 달라지기 때문에.

t검정에서는 var.test를 썼는데 여기선 bartlett.test를 씀

여기서 나온 p-value로 귀무가설을 기각할 수 없다.

aov(weigth~group, data=aa)는 p-value가 안나옴

bartlett.test(weight~group, data=aa)#등분산검정 
output=aov(weight~group, data=aa)
summary(output)

이렇게 하면… F 통계량에 의한 p-value가 나옴

step3. H0 : μ1=μ2=μ3

각각의 평균값 비교

TukeyHSD(output)

이걸로 비교를 한다. 여기서 각각의 평균값 차이를 확인할 수 있음

→ 얘를 plot으로 표현할 수 있다

plot(TukeyHSD(output))

→ 0.0 넘어가면 평균값 차이가 많이 나는 거… 얘 때문에 p-value를 기각하는 것.

⇒ 여기까지 anova였습니다^^

정규성분포를 따르지 않는 ANOVA

  • 의미가 카테고리라는 것을 알기 위해서 factor를 꼭 붙여줘 줘야된다.
out = aov(HT~factor(food), data = bb)
summary(out)

summary(aov(HT~foood, data=bb)
oneway.test(HT~factor(food), data=bb, var.equal=T)

⇒ 이 두개는 매우 다름

  • oneway.test(HT~food, data =bb, var.equal=F)
  • factor를 넣든 안넣든 자유도가 같음 one way는
  • 등분산일 때 연속형 변수일 땐 aov에서 factor를 써야 되고
  • oneway → 등분산이 아닐 때 하는 anova함수
  • 분산이 같을 땐 aov

 

  • B 데이터 → 새로운 변수에 aggregate( BW_value~pid, data = b, mean)에 따른 몸무게 평균 구하기
  • ⇒ pid 에 따른 BW_value 데이터셋이 생성됨
  • 새로운 변수에 aggregate( hab1c~pid, data = a, mean)
  • ⇒ pid에 따른 had1c 평균값 데이터셋이 생성됨
  • A하고 B subset하고 merge를 해줌dang = merge(dang, C, by=’pid’)
  • dang = merge(A, B, by=’pid’)
  • dang=na.omit(dang)
  • na인 행 삭제
  • 당뇨 판정당뇨인 그룹 1, 아닌 그룹 0
  • dang$group←ifelse(dang$hba1c≥6.5, 1, 0)
  • shapiro.test(dang$BW_value) # 정규성 검정
  • var.test(BW_value~group, data=dang, conf.level=0.95)
    • 등분산 검정
  • t.test(BW_value~group, data=dang, var.equal=TRUE, conf.level=0.95)
    • 평균이 같은지 다른지 알아볼 수 있다.
  • dang$age_group←ifelse(age<40, 0, ifelse(age<60, 1, 2))
  • shapiro.test(dang$age) #정규성검정
  • bartlett.test(BW_value~age, data=dang) #등분산검정
  • output = aov(BW_value~factor(age_group), data=dang)
  • summary(output)

n m+- sd p-value

       
그룹 N 평균 +- 표준편차 p-value
ctrl 10 5.0 +- 0.58 0.0159
trt1 10 4.7 +- 0.79  
trt2 10 5.5 +- 0.44  

 

 

 

시험 전 급하게 정리한 글이어서 난잡하지만 R 공부를 기록합니다...