MiniProject 1. densityplot() 함수 만들기
이번 MiniProject의 목표는 Kernel density plot을 그리는 함수를 만드는 것이다.
서론
Kernel density 추정은 무작위변수에서 probability density function을 추정하는데 쓰이는 비모수적인 방법으로 hist()와 같이 그려 본 바 있다. 일반적으로 kernel desity plot은 연속변수에서 분포를 보는 효과적인 방법이다. 사용법은 다음과 같다.
plot(density(x))
여기서 x는 숫자로 된 벡터이다. NA 값이 있을 경우 error가 나므로 na.rm=TRUE로 지정해주어야 한다. 다음의 예를 보자. 위의 그림은 기본 옵션으로 그림을 그린 것이고 아래 그림은 제목과 색을 지정해서 rugplot과 같이 그린 그림이다.
require(mytable)
d=density(acs$BMI,na.rm=TRUE)
plot(d)
또한 density plot은 group 별로 density를 구별하여 보여주는 데 사용할 수 있다. 예를 들어 성별에 따른 몸무게와 BMI의 차이는 다음과 같다.
dm=density(acs$weight[acs$sex=="Male"],na.rm=T)
df=density(acs$weight[acs$sex=="Female"],na.rm=T)
plot(df,main="Density Plot of Weight by Sex",col="blue",lty=1)
lines(dm,col="red",lty=2)
legend("topright",legend=c("Male","Female"),col=c("red","blue"),lty=2:1)
다음 예로 진단명에 따른 LDCL 농도를 구별하여 보여주려면 다음과 같이 할 수 있다.
d1=density(acs$LDLC[acs$Dx=="Unstable Angina"],na.rm=TRUE)
d2=density(acs$LDLC[acs$Dx=="NSTEMI"],na.rm=TRUE)
d3=density(acs$LDLC[acs$Dx=="STEMI"],na.rm=TRUE)
plot(d1,col="green",lwd=2,lty=1,main="Density Plot of LDLC by Dx")
lines(d2,col="blue",lwd=2,lty=2)
lines(d3,col="red",lwd=2,lty=3)
legend("topright",legend=c("Unstable Angina","NSTEMI","STEMI"),
col=c("green","blue","red"),lty=1:3)
문제
R로 분석하거나 프로그래밍을 할때 같은 일을 반복하지 말라는 말이 있읍니다.
Do Not Repeat Yourself !!
densityplot이라는 함수를 하나 만들어 위와 같은 density plot을 그리게 하면 좋을 것이다. 이번 MiniProject의 목표는 Kernel density plot을 그리는 함수를 만드는 것이다.
사용법은 다음과 같다.
densityplot(age~sex,data=acs)
또다른 사용 예로 그림의 제목이나 x축의 라벨(xlab)을 정해줄 수도 있다. 또한 변수 이름을 다소 잘못 쓰더라고 고쳐서 함수를 실행한다.
densityplot(ldl~dx,data=acs,main="Distribution of LDLC in ACS patients",
xlab="Level of LDL-C (mg/dL)") ##
## ' dx ' is an invalid column name: Instead ' Dx ' is used
##
## ' ldl ' is an invalid column name: Instead ' LDLC ' is used
문제 풀이
density curve를 그리는 함수는 plot(density(x))인데 이때 x는 numeric vector 즉 숫자로 된 벡터이다. 예를 들어 density(acs$age)등과 같이 호출해주어야 한다. 하지만 우리가 만들 densityplot()함수는 densityplot(age~sex,data=acs)과 같이 formula를 써서 호출할 것이다. 따라서 이 함수의 핵심은 두가지 이다. 첫째, 함수를 호출하는 formula에서 우리가 원하는 값(“age”, “sex”, “acs”)을 추출하는 것이고 둘째, 우리가 원하는 그래프를 그리는 것이다.
먼저 formula에서 우리가 원하는 데이타프레임의 “열” 이름을 추출하기 위해 함수를 다음과 같이 정의해야 한다. 규정식을 다루기 위한 함수는 모두 이와 같은 형식을 취하는 것이 좋다.
densityplot=function(formula, data, main="",xlab="",ylab="",...){
call=paste(deparse(formula),", ","data= ",substitute(data),sep="")
f=formula
myt=terms(f,data=data)
y=as.character(f[[3]])
x=as.character(f[[2]])
...
...
} 이와 같이 함수를 만들면 densityplot(age~sex,data=acs)과 같이 호출할 경우 y에는 “sex”가 x에는 “age”가 된다. 다음으로 할 일은 y가 범주형변수(factor)일 수도 있지만 문자형 변수일 수도 있다. 따라서 범주형변수인지 체크해봐서 범주형변수가 아니면 범주형변수로 바꾸어 group이라는 변수에 넣어준다.
if(is.factor(data[[y]])) group=data[[y]]
else group=factor(data[[y]])이제 group이라는 범주형 변수의 level이 몇개인지 세서 count에 넣어주고 그 갯수에 맞추어 그림을 그릴 색깔을 rainbow()함수로 지정해준다.
count=length(levels(group))
colors=rainbow(count)그래프를 그리는 순서는 먼저 plot()함수를 써서 빈 그래프를 그린후에 lines()를 써서 density 커브를 그리면 되는데 plot()함수를 호출하려면 x축과 y축의 최대값,최소값을 입력해주어야 한다. 따라서 x축의 최대, 최소값은 우리가 그리고자 하는 변수 여기서는 acs$age의 최대값,최소값이다. acs$age는 acs[[“age”]]로 쓸수 있다.
xrange=range(data[[x]],na.rm=TRUE)density curve의 y축의 값은 지금은 알 수 없다. y축의 최대값을 저장할 변수 maxy를 0으로 초기화시킨다. 전체 데이타를 아까 만든 범주형변수인 group의 level에 따라 서브셋으로 나눈 후(예를 들어 sex가 male인 경우와 female인 경우) 그 서브셋데이타의 x열의 density를 d에 저장한다. 그 전에 count갯수만큼의 density를 저장하기 위해 density의 list를 초기화 시킨후(dl=list()) 그 dl의 i번째에 density를 저장한다.
dl=list()
maxy=0
for(i in 1:count) {
subdata=subset(data,data[[y]]==levels(group)[i])
d=density(subdata[[x]],na.rm=T)
maxy=max(maxy,max(d$y))
dl[[i]]=d
}이제 plot()함수로 빈그래프를 그리는데 여유를 주기 위해 최소값은 x의 최소치보다 약간 작게, 최대값은 x의 최대치보다 약간 크게 그린다. 그리고 아까 만든 density list에 있는 density를 따라 lines를 이용해 선을 그린다. 선의 색깔은 아까 rainbow()를 이용해 지정해준 색깔로 그리고 범례를 그려 완성한다.
plot(c(xrange[1]*0.85,xrange[2]*1.1),c(0,maxy*1.1),
main=main,xlab=xlab,ylab=ylab,type="n",...)
for(i in 1:count) lines(dl[[i]],col=colors[i],lty=i,lwd=2)
legend("topright",legend=levels(group),col=colors,lty=1:count)이상이 내가 만든 densityplot() 함수의 주요 내용이고 실제 함수는 축의 제목을 정하는 부분도 있고 함수를 호출할떄 에러를 다루는 부분, 열의 이름이 부정확할때 정확한 열이름으로 바꿔주는 부분 등이 있어 조금 길어졌다. 다음은 전체 소스이다.
require(mytable)
densityplot=function(formula,data,main="",xlab="",ylab="",...){
call=paste(deparse(formula),", ","data= ",substitute(data),sep="")
f=formula
myt=terms(f,data=data)
y=as.character(f[[3]])
y=unlist(strsplit(y,"+",fixed=TRUE))
if(length(y)>1) {
cat("\n","Only one dependent variable is permitted\n")
return(invisible())
}
y1=validColname(y,colnames(data))
if(is.na(y1)) {
cat("\n","There is no column named '",y,"' in data ",
substitute(data),"\n")
return(invisible())
}
if(!identical(y,y1)) {
cat("\n","'",y,
"' is an invalid column name: Instead '",y1,"' is used\n")
s=paste(as.character(f[[2]]),y1,sep="~")
result=densityplot(as.formula(s),data)
return(invisible())
}
x=as.character(f[[2]])
x1=validColname(x,colnames(data))
if(is.na(x1)) {
cat("\n","There is no column named '",x,"' in data ",
substitute(data),"\n")
return(invisible())
}
if(!identical(x,x1)) {
cat("\n","'",x,
"' is an invalid column name: Instead '",x1,"' is used\n")
s=paste(x1,as.character(f[[3]]),sep="~")
result=densityplot(as.formula(s),data)
return(invisible())
}
if(is.factor(data[[y1]])) group=data[[y1]]
else group=factor(data[[y1]])
count=length(levels(group))
colors=rainbow(count)
if(nchar(main)==0) main=paste("Distribution of '",x,"' by '",y1,"'")
if(nchar(xlab)==0) xlab=x
if(nchar(ylab)==0) ylab="Density"
xrange=range(data[[x]],na.rm=TRUE)
dl=list()
maxy=0
for(i in 1:count) {
subdata=subset(data,data[[y]]==levels(group)[i])
d=density(subdata[[x]],na.rm=T)
maxy=max(maxy,max(d$y))
dl[[i]]=d
}
plot(c(xrange[1]*0.85,xrange[2]*1.1),c(0,maxy*1.1),
main=main,xlab=xlab,ylab=ylab,type="n",...)
for(i in 1:count) lines(dl[[i]],col=colors[i],lty=i,lwd=2)
legend("topright",legend=levels(group),col=colors,lty=1:count)
}결론
그래프를 그리는 함수를 공개하는 이유는 저자가 만든 함수의 코드가 최선의 코드는 아니라고 생각하기 때문이다. 개선할 점이 얼마든지 있을 수 있다. 한가지 당부하고 싶은 것은 다른 사람이 만든 함수의 소스를 보기 전에 먼저 스스로 함수를 만들어보는 것이 중요하다. 함수를 만들다보면 벽에 부딪혀 잘 안되는 부분이 있을 수 있고 스스로 해결하기 위해 고민하다가 다른 사람의 소스를 보고 막힌 벽이 뚫릴 때가 있다. 마치 뉴턴이 떨어지는 사과를 보고 만유인력의 법칙을 발견했듯이. 이 함수를 만들 때에 저자는 rainbow()함수를 보고 감탄했었다. 무지개 함수를 모를 때에는 그저 색깔을 1,2,3,4,5,…순으로 정했었는데 다른 사람이 만든 다른 함수의 소스를 보다가 rainbow()함수를 보고 정말 막혔던 부분이 훤히 뚫리는 느낌이었다. 함수를 만드는 과정은 시간이 많이 걸리지만 이후에 훨씬 많은 시간이 절약되며 내가 고생해서 만든 함수가 제대로 작동하는 것을 보면 아끼는 후배/제자가 제 몫을 훌륭하게 해낼때 느끼는 보람 비슷한 것을 느낄 수가 있다. 이 소스를 보고 개선할 점이 있으면 언제든지 게시판을 통해 의견을 올려주기 바란다.