ケィオスのRメモランダム

Rを使いこなすためのちょっとしたメモです

Quantile regression (分位点回帰)の数値実験

数学的な背景の解説は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スクリプトを実行した例が以下です. 

Rスクリプトの実行例

なぜ,Quantile regression

 私は,最近以下の2本の論文を出版しました.

Frontiers | Age- and height-dependent bias of underweight and overweight assessment standards for children and adolescents

Allometric multi-scaling of weight-for-height relation in children and adolescents: Revisiting the theoretical basis of body mass index of thinness and obesity assessment | PLOS ONE

この論文では,子どもの痩せや肥満評価で使われている基準が,必ずしも正しくないことを示しています.例えば,日本では学校保健統計の痩身傾向,肥満傾向の基準があり,それが正しいと思ってみんな使っていますが,その判定が酷く間違った結論を導く場合があります.痩せて助けが必要な子どもが,日本にもたくさんいるのに,学校保健統計の基準がいい加減なために,見過ごされています.悲しすぎます.このことについて,今後,解説していきたいと思います.

 この論文でQuantile regressionを使ったので,まずは,参考としてRスクリプトを掲載しました.