【RでWorm plot】正規性を可視化して確認 - ケィオスの時系列解析メモランダム

データが正規分布に従うかどうかを確認するときにQQ plotっていうのを使うことがあります (QQ plotは正規分布以外でも使えますが)．最近私は，身長，体重，BMIを分析する研究をしていています．これは，痩せすぎとか，太りすぎとか，そういった人を判定することを目的にしています．この研究に関連したWHOの報告書などで正規性を可視化するのに"Worm plot"っていうのが使われていることがあります．私もWorm plotを描きたいと思い，パッケージに頼らずにWorm plotを描いてみました．

　ということで，Worm plotを描くRスクリプトはこれです．ここでは，xにデータベクトルが入っているとします．

# xにデータを格納（ここでは例として，正規乱数を生成しています）
x <- rnorm(1000)
# データの長さ
n <- length(x)
# 標準化
x.stdzd <- scale(x)
#####################
# プロット位置の分位数を計算する関数
qt.plot <- function(k,n){
 return((k-0.5)/n)
}
#####################
# normal quantileの位置
qt <- qnorm(p=qt.plot(1:n,n))
# 小さい順に並び替え
x.stdzd.o <- sort(x.stdzd)
##############
# zスコアとnから標準誤差を計算
z.SE <- function(z, n) {
  sqrt(pnorm(z)*(1-pnorm(z))/n)/dnorm(z)
}
###
x.abs.max <- ceiling(max(abs(x.stdzd),na.rm=TRUE)/0.5)*0.5
x.range <- c(-x.abs.max,x.abs.max)
z.plot <- seq(-x.abs.max,x.abs.max,length.out = 500)
sd.x <- sd(x.stdzd)
SE <- sd.x*z.SE(z.plot, n)
###
# トレンドを除く
x.detrend <- x.stdzd.o - sd.x*qt
###
par(mar=c(5,5,2,2),las=1,cex.lab=1.5,cex.axis=1.3)
plot(qt,x.detrend,pch=16,col=4,lwd=1.5,xlab="Unit normal quantile",ylab="",xlim=x.range,xaxs="i",yaxs="i",ylim=c(-0.5,0.5))
mtext("Deviation",side=2,line=3.5,cex=1.6,las=3)
lines(z.plot,1.96*SE,col="gray10",lwd=1.5,lty=2)
lines(z.plot,-1.96*SE,col="gray10",lwd=1.5,lty=2)
abline(h=0,col="gray50",lwd=1.5,lty=2)
abline(v=0,col="gray50",lwd=1.5,lty=2)

　このスクリプトを実行すると，下のような図が描かれます．