数学的な背景の解説はWikipediaに任せておき,Rで実際にQuantile regression (分位点回帰)を実行するスクリプトを載せておきます.
以下のスクリプトでは,相関係数 r
の2変数が2次元正規分布に従うことを仮定し,標本を生成します.標本数は N
で指定します.標本数を増やせば,推定結果が理論と一致しやすくなります.MASS
パッケージを使っていますので,事前にインストールしておいてください.
##################### # numerical experiment ########### # r: correlation coefficient r <- 0.6 ####################### sig.x <- 1 sig.y <- 1 ##################### # install the “MASS” package in advance. # install.packages("MASS") require(MASS) ##################### # data length N <- 1000 ### # Samples from Bivariate Normal Distribution SIG <- matrix(c(sig.x^2,r*sig.x*sig.y,r*sig.x*sig.y,sig.y^2),nrow=2,ncol=2,byrow=TRUE) eps <- mvrnorm(N,c(0,0),SIG) x <- eps[,1] y <- eps[,2] ###################### # Draw sample points plot(x,y,pch=16,col=gray(0.7),xlim=c(-4.5,4.5),ylim=c(-4.5,4.5),main=sprintf("q = %f",q)) # Estimation of contour lines hist2d <- kde2d(x,y, n = 100) contour(hist2d$x, hist2d$y, hist2d$z,levels=c(0.1,0.01),add=TRUE,col="gray10",drawlabels=FALSE) ##################### # model for quantile regression f.bmi <- function(arg) { a<- arg[1] b<- arg[2] model <- a*x+b res <- y - model res[res >= 0] <- res[res >= 0]*q res[res < 0] <- res[res < 0]*(1-q) return(sum(abs(res))) } ##################### # Quantile [0,1] ##################### q <- 0.025 ##################### # Estimated parameters tmp <- optim(c(1,2),f.bmi) a <- tmp$par[1] b <- tmp$par[2] ##################### # Confirmation sig <- sig.x*sig.y*sqrt(1-r^2) xp <- qnorm(q) # theoretical prediction: blue solid line curve(xp*sig+r*sig.y/sig.x*x,add=TRUE,col=4,lwd=2) # estimation: red dashed line curve(a*x+b, add=TRUE, col=2,lwd=2,lty=2) ### legend("bottomright",legend=c("Sample points","Theory","Quantile Regression"),pch=c(16,NA,NA),col=c(gray(0.7),4,2),lty=c(NA,1,2),lwd=c(NA,2,2))
このRスクリプトを実行した例が以下です.
なぜ,Quantile regression
私は,最近以下の2本の論文を出版しました.
この論文では,子どもの痩せや肥満評価で使われている基準が,必ずしも正しくないことを示しています.例えば,日本では学校保健統計の痩身傾向,肥満傾向の基準があり,それが正しいと思ってみんな使っていますが,その判定が酷く間違った結論を導く場合があります.痩せて助けが必要な子どもが,日本にもたくさんいるのに,学校保健統計の基準がいい加減なために,見過ごされています.悲しすぎます.このことについて,今後,解説していきたいと思います.
この論文でQuantile regressionを使ったので,まずは,参考としてRスクリプトを掲載しました.