[R] 기본 명령어부터 ANOVA 분석까지

한글 인코딩 깨질 때 파일 인코딩 변경 방법(MAC)

cp949 -> UTF-8 인코딩 변경
iconv -f cp949 < 원본파일 > 저장파일

iconv -f cp949 < source.txt > result.txt
출처: https://kka7.tistory.com/41 [때로는 까칠하게..:티스토리]

샘플데이터 - merge, 조건, ifelse 연습

1_EMR data_sex.txt

0.03MB

2_App data_e.txt

4.92MB

dat1<-read.table(file="C:/Users/cse/Desktop/1_EMR data_sex.txt",header=TRUE,sep="\t",fill=TRUE)
# dat1<-read.table(file="/Users/chongin12/Documents/R_Practice/practice/1_EMR data_sex.txt",header=TRUE,sep="\t",fill=TRUE)
install.packages("descr")
library(descr)

dat1<-dat1[complete.cases(dat1),] # NA 제거

freq(dat1$sex)

dat2<-read.table(file="C:/Users/cse/Desktop/2_App data_e.txt",header=TRUE,sep="\t",fill=TRUE)
# dat2<-read.table(file="/Users/chongin12/Documents/R_Practice/practice/2_App data_e.txt",header=TRUE,sep="\t",fill=TRUE)
dat22<-subset(dat2,select=c(A_myhealth_id,A_steps))
dat22$pid=dat22$A_myhealth_id
View(dat22)

dat222<-dat22[,c(2:3)]
View(dat222)

res<-merge(dat1,dat222,key=pid)
View(res)

summary(res$A_steps) # 이상함. min값이 음수가 나옴.

t.test(A_steps ~ sex, data=res, var.equal=TRUE, conf.level = 0.95) # 여자, 남자 그룹별로 평균 구하기

#조건 주기 (패키지 사용)
#install.packages("dplyr")
#library(dplyr)
#res_1 <- res %>% filter(A_steps>0 & A_steps<10000)

#조건 주기 (패키지 없이)
res_1 <- res[res$A_steps>0 & res$A_steps<10000,]

#히스토그램
hist(res_1$A_steps)
hist(res_1$A_steps, breaks = seq(0,10000,by=100)) # 좀 더 자세히

#평균값(3450)을 기준으로 그룹 짓기
res_2<-res_1
res_2$A_steps_group<-ifelse(res_2$A_steps>=3450.1,1,0)

# 0~3333, 3333~6666, 6666~10000 그룹짓기
res_3<-res_1
res_3$A_steps_3group<-ifelse(res_3$A_steps>=6666.1,2,ifelse(res_3$A_steps>=3333.1,1,0))
# with로 재사용 줄이기 res3$A_steps 대신 res3 사용한다고 선언해두고 A_steps로만 사용.
res_3$A_steps_3group<-with(res3,ifelse(A_steps>=6666.1,2,ifelse(A_steps>=3333.1,1,0)))

# aggregate 함수 사용해서 한 사람당 하나의 데이터만 남길 수 있다.(한 사람당 여러개의 데이터가 있음)
# res_4<-aggregate(res_1$sex, list(res_1$pid), mean)
# 살짝 잘못돼서 주석처리. 사용 x

 

샘플 데이터(aggregate, shapiro.test, F test, T test)

download.txt

1.38MB

# 데이터 불러오기
r2<-read.table(file="C:/Users/cse/Desktop/download.txt",header=TRUE,sep=" ",fill=TRUE)

# 결측 처리 (하지만 이미 데이터에 결측이 없어서 바뀌지 않음)
r2<-r2[complete.cases(r2),]

# 데이터를 1~9999으로 자르기
r3<-r2[r2$A_steps>0 & r2$A_steps<10000,]

# r4생성
r4<-r3

# r4에 새로운 행 생성 : A_steps_group_2 : 0~100 : 0, 101~ : 1
r4$A_steps_group_2<-ifelse(r3$A_steps>=100,1,0)

# r4에 새로운 행 생성 : A_steps_group_3 : 0~99 : 0, 100~4999 : 1, 5000~ : 2
r4$A_steps_group_3<-ifelse(r4$A_steps<100,0,ifelse(r4$A_steps<5000,1,2))

# r5 : 각 사람마다 A_steps 평균으로 압축.
r5=aggregate(A_steps ~ pid, data=r4, mean)

# r6 : 각 사람마다 성별 나타내기
r6<-aggregate(sex ~ pid, data=r4, mean)

# r7 : 각 사람마다 A_steps 평균과 성별이 나타나게 merge
r7<-merge(r5,r6,key=pid)

# shapiro.test
shapiro.test(r7$A_steps) # 결과 : W = 0.99115, p-value = 0.07031
## 특이하게 shapiro test는 H0이 정규분포가 맞다, H1이 정규분포가 아니다. 라고 세워져있음.
## p-value = 0.07031 > 0.05 이므로 정규분포라고 할 수 있음. (영가설 기각할 수 없음)
## 2000개 이하 표본에 대해 정규분포인지 테스트 해보려면 Shapiro-Wilk test(, 아니면 Kolmogorov-Smirnov test를) 하면 된다.

# F test (분포 테스트 : Variance 테스트)
var.test(A_steps ~ sex, data=r7, conf.level = 0.95)
## H0 : var1 == var2
## H1 : var1 != var2
## test 결과 : 0.05보다 큼. H0을 기각할 수 없음.
## 따라서 남자의 걸음걸이의 분포와 여자의 걸음걸이의 분포(퍼져있는 정도)는 같다.

# T test
t.test(A_steps ~ sex, data=r7, var.equal=TRUE, conf.level=0.95)
## t.test(영향을 받는 변수(종속변수) ~ 영향을 주는 변수(독립변수), dataset, ...)
## two sample T test
## p-value = 0.001105
## group 1: male, group 2: female
## group 1 mean : 3736.398
## group 2 mean : 3122.812
## H0 : male의 평균 걸음 수와 female의 평균 걸음 수는 같다
## H1 : male의 평균 걸음 수와 female의 평균 걸음 수는 다르다.
## 결과 : p-value가 0.05보다 작기 때문에 h0 기각, 평균 걸음 수는 다르다.

### 참고) p-value가 작을수록 영가설(H0)을 기각할 힘이 생긴다.
### p-value가 0.05보다 작아지면 영가설(H0)을 기각하고 귀무가설(H1)을 채택할 수 있다.
### 귀무가설(H1)이 내가 주장하고 싶은 가설이다.

t.test에서 var.equal=TRUE와 FALSE의 차이? : 등분산검정을 통해 얻은 값을 넣어주면 된다.

등분산검정? var.test(종속변수 ~ 독립변수, dataset, …) : 두 분산에 대한 F test.
H0 : 분산이(퍼져있는 정도가) 같다
H1 : 분산이 다르다.
이 결과(p-value)가 0.05보다 작으면 기각. var.equal=FALSE
0.05보다 크면 기각할 수 없음. var.equal=TRUE

#0. distribution of data
install.packages("descr")
library(descr)

#0.1. categorical variable: freq
a2<-freq(a$MH)
View(a2)

#0.2. continuous variable : graph
hist(a$HT)
head(hist(a$HT))

#0.2.1. 간격 조정(breaks)
hist(a$HT, breaks=seq(130, 200, by=1))

#0.3. 정규성 검정 (정규분포하는지?)
qqnorm(a$HT) # Q-Q plot은 직선에 가깝게 분포할수록 정규분포를 따름.
qqline(a$HT)

ANOVA : Analysis of variance

먼저 데이터가 5000개 이상인지 확인.

1. 데이터가 5000개 이상이면 중심극한정리에 의해 정규분포 한다.
2. 데이터가 5000개 이하이면 shapiro wilk 테스트를 통해 정규분포인지 검사한다.
   1. 특이하게 shapiro wilk 테스트 결과가 0.05 이상이어야 정규분포이다.

데이터가 정규분포한다면, t 검정 또는 AVOVA 검정을 할 수 있다. 즉 통계적으로 유의미한 차이가 있는지 검사한다.

그 전에 먼저 var test를 해야 한다.

1. var test의 결과(p-value)가 0.05보다 작으면 기각. (false)
2. var test의 결과가 0.05보다 크면 true

이제 통계적으로 유의미한 차이가 있는지 검사한다.

1. 그룹이 2개이면 (ex. 매일 식사하는 군 vs 매일 식사하지 않는 군) T 검정을 한다.
2. 그룹이 3개 이상이면 (ex. 일주일에 0번 식사 vs 일주일에 1~6번 식사 vs 일주일에 7번 식사) ANOVA 검정을 한다.

# ANOVA 분석
out = aov(HT ~ factor(F_BR_3), data=a) # group 변수가 숫자로 되어 있으므로 factor로 해주어야 함.
summary(out)

# 사후 검정
#install.packages('agricolae')
#library(agricolae)
TukeyHSD(out)

# 그래프 그려보기
plot(TukeyHSD(out))

ANOVA 분석 후 나오는 결과에서 Df는 자유도를 의미한다. 결과는 F값을 보면 되는데 0.05보다 작으면 기각, 즉 서로 다른 그룹 간의 평균 차이가 유의미하다는 것을 뜻한다.

사후 검정으로는 Tukey를 사용한다. 각 그룹간 p adj를 보면 되는데 각각 0.05보다 작으면 유의미한 차이가 있다고 생각하면 되고, 그래프로 그리면 구간이 나오는데 모두 0을 포함하지 않으면 모두 유의미한 차이가 있다는 것이다.

---

내가 직접 해보기

dat2의 A_goal_steps에 따른 A_steps 유의미한 차이가 있는가?

A_goal_caloire는 0, 7500, 10000임.

1. 먼저, dat1에 대해 A_myhealth_id, A_steps, A_goal_steps만 남기기.
2. 결측 모두 제거
3. steps를 1~9999까지 자르기
4. A_myhealth_id와 A_steps를 기준으로 aggregate
5. A_goal_steps를 numeric으로 변환
6. A_myhealth_id와 A_goal_steps를 기준으로 aggregate
7. a4와 a6 merge
8. A_goal_steps_group을 생성하고, 0=0 ~ 7499, 1=7500 ~ 9999, 2=10000 ~ 나누기
9. A_goal_steps를 없애기
10. A_steps가 정규분포인지 검사 → 정규분포가 아님 따라서 ANOVA 말고 다른걸로 해야 하지만 그냥 연습이니 ANOVA까지.

#1
a1<-subset(dat2,select=c(A_myhealth_id,A_steps,A_goal_steps))

#2
a2<-a1[complete.cases(a1),]

#3
a3<-a2[a2$A_steps>0 & a2$A_steps<10000,]

#4
a4<-aggregate(A_steps ~ A_myhealth_id, data=a3, mean)

#5
a5<-a3
a5$A_goal_steps<-as.numeric(a5$A_goal_steps)

#6
a6<-aggregate(A_goal_steps ~ A_myhealth_id, data=a5, mean)

#7
a7<-merge(a4,a6,key=A_myhealth_id)

#8
a8<-a7
a8$A_goal_steps_group<-with(a8,ifelse(A_goal_steps<7500,0,ifelse(A_goal_steps<10000,1,2)))

#9
a9<-subset(a8,select=c(A_myhealth_id,A_steps,A_goal_steps_group))

#10
qqnorm(a9$A_steps)
qqline(a9$A_steps)

shapiro.test(a9$A_steps) # 불가. 원래 다른걸로 해야하지만 연습이니 그냥 ANOVA 고고

aout = aov(A_steps ~ factor(A_goal_steps_group), data=a9)
summary(aout)

TukeyHSD(aout)
plot(TukeyHSD(aout))