R을 쓰면서 중요한 것 중 하나는 ‘~’ 에 대한 이해다.
A ~ B 라는 것은, (내가 이해하기로) B라는 요인에 대한 A 값 이라고 해석할 수 있는데,(사실 맞는지는 모름)
One-way ANOVA라는 것이 결국 명목형 자료에 대한 연속형 자료의 통계라는 것,
아무튼 햇갈리면 다음과 같이 그래프를 만들어서 비교해보자
boxplot(A~B, data=myDataFrame); //myDataFrame은 read.csv, read.xlsx등으로 불러온 자료(데이터프레임)이름을 붙이면된다.
(boxplot(myDataFrame$A~myDataFrame$B)와 같은 말이다.)
boxplot(B~A, data=myDataFrame);
이해가 될 것이다.
1. DataFrame factor 만들기
사실 R 자체에서는 통계를 구함에 있어 이 변수가 명목형인지 연속형인지 신경써주지 않는다.(내가 ‘의도하는 바’는 위의 0,1,2,3,4가 실제 숫자아 아니라 index값이라는 것, 그러나 프로그램은 그냥 숫자 0,1,2,3,4로 이해한다)
따라서 내가 원하는 것은 앞의 그래프이고, x축을 다르게(명목변수를 제대로된 label을 붙여서!) 표현하고 싶으니,
> data$ss<-factor(data$ss,levels=c(‘0′,’1′,’2′,’3′,’4’),labels=c(‘worst’,’worse’,’same’,’better’,’best’))
해서 보면 각각 0,1,2,3,4를 factor로 만들고, 이에 대한 label을 붙이게 된다. 아마도(!) 내부적으로는 0,1,2,3,4라는 값은 살아있는 듯 하다.
변수를 확인할 수 있는, str함수를 쓰면
값은 살아있는 듯 하다.
아무튼 이제 이 상태에서 boxplot을 내보면
ss값을 먼저 쓴 boxplot은 에러가 난다.
2. One-way ANOVA
이제 위에서 만든 factor들에 대한 cc의 값들이 통계적인 유의성을 보이는지 확인해보자.
라고 해봤자, ANOVA는 SPSS로 돌렸을테니까…및 R로 ANOVA를 확인할 사람은 이미 더 실력이 있다는 가정하에 생략
통계는 잘 모르지만, 예전에 통계강좌에서 post huc으로 Scheffe를 이용하라고 했어서, R로 확인한 한번 해본다.
참고로 아래 함수를 쓰려면 libary(asbio) 로딩이 필요하다.
> result = scheffeCI(data$cc,data$ss,conf.level=0.95,MSE=NULL,df.err=NULL)
역시 (내가) SPSS에서 봤던 것 처럼 same과 better는 유의한 차이를 보이지 않는다.
참고로 이 결과의 표에서 “P-value를 추출하고 싶다!”면 (추후에 그래프에 포함시키기 위해서라든지)
찾아가는 방법은 다음과 같다.
> str(result);
보면 우리가 원하는 항목은 ‘result 변수’의 항목 중 $summary라는 곳이다.
그럼 $summary를 보자
그럼 여기에 Adj. P-value가 있다.
> result$summary$Diff 라든지 result$summary$Lower 는 그냥 치면 결과가 잘 나온다.
그럼 원하는 위치의 값을 result$summary$Diff[1] , Diff[2]등으로 써서 읽어올수가 있다.
그러나!! Adj. P-value는 안됨.
이때 등장하는 것이 dataframe에 대한 access의 특징
띄어쓰기가 있는 항목은 다음과 같이 접근하면 된다.
> result$summary[,”Adj. P-value”]
[1] 0 0 0 0 0 0.999994 0 0 0.000518 0.000388
Levels: 0 0.000388 0.000518 0.999994
> result$summary[6,”Adj. P-value”]
[1] 0.999994
Levels: 0 0.000388 0.000518 0.999994
> result$summary[,”Adj. P-value”][6]
[1] 0.999994
Levels: 0 0.000388 0.000518 0.999994
>
참고로 [,”Adj. P-value”]앞에 ‘,’ 앞에 아무것도 안쓰면 모든 테이블 값을 가져오는 것이고,
result$summary[6,”Adj. P-value”] 와 result$summary[,”Adj. P-value”][6] 의 값이 왜 같은지는 쉽게 알 수 있으리라 생각한다.
그럼 다시 돌아가서 one-way anova그래프를 완성해보자
3. One-way ANOVA 그래프 그리기
위에서 잠깐 사용했던 Boxplot을, 조금더 적절하게 그리기 위해 ggplot을 사용한다.
library(ggplot2)를 로드해야한다.
일단 사용할 데이터를 정의한다.
ggplot(myDataFrame, aes(x=ss,y=cc)) 와 같은 방식이다.
현재 표현하는 방식에 대해 정의하지 않았으므로 나오는 것이 없다.
위의 두줄은 동일한 결과를 나타낸다.이후
>fig 라고 치면
(참고로 보이는 검은색 띠는 아직 출판되지 않은 자료라 가려놓은 것이다. 원래는 숫자가 나온다.)
Box-plot이 그려진다.
이제 여기에 몇가지 설정을 할 수 있는데
fig <- fig + …. 이런 방식으로 설정해나가면 된다.
예를 들자면 아까 geom_boxplot() 처럼 + geom_text 함수를 이용해서 값을 표시할 수 있고, (자세한 것은 구글을 참조할 것)
위에 그래프에서 보이는 오차범위 밖의 검은색 점들을 지우려면 다음과 같이 쓰면 된다.
> fig <- ggplot(data,aes(x=ss,y=cc)) + geom_boxplot(outlier.colour = NA) //혹은 outlier.size=0 으로 해버리는 방법도 있다.
밖에 있던 점들을 outlier라고 부르는 것 같다.
이제 배경을 하얗게 바꾸고 가로실선을 회색으로 바꾸고, 세로 실선은 없어지게 한다.
> fig <- ggplot(data,aes(x=ss,y=cc)) + geom_boxplot(outlier.colour = NA)
> fig <- fig + theme_bw()+
+ theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_line(colour=”grey60″,linetype=”dashed”))
> fig
여기에 추가로 Y 축 표시를 figure 창 사이즈와 무관하게 (windows R에서는 어떻게 표시는지 잘 모르겠지만, 맥에서 R을 돌리면 figure 창 사이트 조절에 따라 y축 표시 숫자가 바뀐다)
fig <- fig + scale_y_continuous(breaks=seq(0,13,by=2)) (y축 표시를 0부터 13까지 2씩 증가시켜 표시함)
fig <- fig + scale_y_continuous(breaks=c(1:10)) (y축 표시를 1부터 10까지 1씩 증가시킴)
등등으로 바꿔서 사용
4. 통계적 의미를 위한 처리
“박스를 그려놓고 각 박스별로 통계적으로 의미있는 경우 선으로 보여줘서 처리하고 싶다.” 라고 하면
http://stackoverflow.com/questions/19658388/error-with-adding-geom-path-to-boxplot-using-ggplot-function
예제를 참조.
다시는 못찾고 있는데, 이전에 다른 예제에서는 각 박스별로 y좌표의 최고점을 찾아 편하게 처리했던 예제도 있었다.
만약에 Y-maximum값이 Y-axis의 값에 거의 근접해서 그래프 위쪽에 여유가 별로 없다면?
다음과 같은 꼼수를 이용한다.
fig2<- fig + coord_cartesian(ylim=c(-1,20)) + scale_y_continuous(breaks=seq(0,12,by=2))
coord_cartesian을 이용해서 y-axis의 범위를 정하고 뒤의 scale_y를 이용해서 표시해주는 영역은 0에서 12까지만으로 제한
그럼 위쪽 공간을 활용하기가 쉬워진다.
————–
참고 예제
http://docs.ggplot2.org/0.9.3.1/geom_boxplot.html
http://stats.stackexchange.com/questions/8137/how-to-add-horizontal-lines-to-ggplot2-boxplot
