Kernel density plots

이 장을 공부하기 위해 mytable 패키지를 설치해야 합니다. mytable 패키지는 github에 있으며 다음과 같이 설치할수 있읍니다.

1. 먼저 devtools패키지를 설치합니다.
2. mytable 패키지를 설치합니다.
3. mytable 패키지를 불어들입니다. acs 데이타를 불러옵니다.

install.packages("devtools")
devtools::install.github("cardiomon/mytable")

require(mytable)

## Loading required package: mytable

data(acs)

Kernel density plots

Kernel density 추정은 무작위변수에서 probability density function을 추정하는데 쓰이는 비모수적인 방법으로 hist()와 같이 그려 본 바 있다. 일반적으로 kernel desity plot은 연속변수에서 분포를 보는 효과적인 방법이다. 사용법은 다음과 같다.

plot(density(x))

여기서 x는 숫자로 된 벡터이다. NA 값이 있을 경우 error가 나므로 na.rm=TRUE로 지정해주어야 한다. 다음의 예를 보자. 위의 그림은 기본 옵션으로 그림을 그린 것이고 아래 그림은 제목과 색을 지정해서 rugplot과 같이 그린 그림이다.

par(mfrow=c(2,1))
d=density(acs$BMI,na.rm=TRUE)
plot(d)
plot(d,main="Kernel Density of Body Mass Index")
polygon(d,col="red",border="blue",lwd=2)
rug(acs$BMI,col="brown")

par(mfrow=c(1,1))

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을 그리게 하면 좋겠읍니다. 다음은 제가 만든 densityplot() 함수의 실제 사용 예입니다.

densityplot(age~Dx,data=acs)

mytable(Dx~age,data=acs)

## 
##           Descriptive Statistics by 'Dx'         
## __________________________________________________ 
##        NSTEMI       STEMI    Unstable Angina   p  
##        (N=153)     (N=304)       (N=400)    
## -------------------------------------------------- 
##  age 64.3 ± 12.3 62.1 ± 12.1   63.8 ± 11.0   0.073
## --------------------------------------------------

또다른 사용 예입니다.

densityplot(mpg~cyl,data=mtcars)

여러분 만의 densityplot 을 만들어보세요. 제가 만든 densityplot() 함수는 일주일 후에 공개할 예정입니다.