分位点ポアソン回帰モデルをRでしよう！

Rのpackageに一般化線形分位点回帰に関連するパッケージないじゃん！で困ってる人向けの記事です．
分位点回帰そのもの理論とかは書きません．いつか書きます．

#普通の分位点回帰

install.packages("quantreg")
library(quantreg)
library(MASS)
mu1 <- c(0,0)
Sigma1 <- matrix(c(1,0.99,0.99,1), nrow=2, ncol=2, byrow=T)
#mu2 <- c(100,0)
#Sigma2 <- matrix(c(0.1,0,0,0.1), nrow=2, ncol=2, byrow=T)
NUM1<-99
mgauss1 <- mvrnorm(NUM1, mu1, Sigma1)

mgauss<-rbind(mgauss1,c(30,0))
x.mgauss<-mgauss[,1]
y.mgauss<-mgauss[,2]

plot(mgauss)
abline(lm(y.mgauss ~ x.mgauss), col="red", lty = 1, lw=1)
abline(rq(y.mgauss ~ x.mgauss), col="blue", lty = 1, lw=1)
legend("bottomright", legend = c("平均値回帰","中央値回帰"), col = c(2,4), pch = c(1,1), lty = c(1,1))

こんな感じのが出来るはずです．平均値回帰はうまくfitしてないけど中央値回帰(分位点回帰)はうまく当てはまります．

#一般化線形回帰したくなるデータWBGT

暑さ指数（WBGT（湿球黒球温度）：Wet Bulb Globe Temperature）は、熱中症を予防することを目的として1954年にアメリカで提案された指標です。 単位は気温と同じ摂氏度（℃）で示されますが、その値は気温とは異なります。暑さ指数（WBGT）は人体と外気との熱のやりとり（熱収支）に着目した指標で、人体の熱収支に与える影響の大きい ①湿度、 ②日射・輻射(ふくしゃ)など周辺の熱環境、 ③気温の３つを取り入れた指標です。

厚生省熱中症予防情報サイトより引用

このWBGTという指標は熱中症の注意喚起によく使われるので扱います．
データはこちら．平成28年５月～９月の日別の救急搬送人員数と暑さ指数(WBGT)データです．

WBGTのcsv

データもとは
http://www.fdma.go.jp/neuter/topics/houdou/h28/10/281012_houdou_2.pdf
のスライド18ページ目です．

初めに一般化線形回帰します

#install.packages("lqmm")
#library(lqmm)

data<-read.csv("WBGT.csv")
data_x<-data[,1]
data_y<-data[,2]

theta_hat <- glm(data_y ~ data_x, family = poisson(link = "log"))
summary(theta_hat)
xxx<-(1:350)/10
#平均値
plot(data_x,data_y, xlab="WBGT",ylab="患者数",cex=.5, xlim=c(18, 35), ylim=c(0, 120))
par(new = T)
plot(xxx,exp(-6.6501+0.3448 * xxx),	xlim = c(18, 35), ylim = c(0, 120), xlab = "", ylab = "",type="l",col=2)

いわるるポアソン回帰です．当てはまりが見た目よいのでWBGTは良さげな指標な感じがします．
多分これで以下のグラフかけます．glmからplotの部分はもっとうまくかける気がする．

#一般化分位点回帰

Rには一般化分位点回帰モデルのパッケージがないのでoptimで頑張ります．

opt.func<-function(a){
	ans<-sum(abs(data_y-exp(a[1]+a[2]*data_x)))
	return(ans)
}


#	sum(data_y-exp(a[1]+a[2]*data_x)*(0.1-(data_y-exp(a[1]+a[2]*data_x)<0)))


opt.func.tau<-function(a){
	long<-length(data_y)
	ans<-rep(0,long)

	data.ev<-(data_y-exp(a[1]+a[2]*data_x))

	for(i in 1:long){
	ans[i]<-(abs(data.ev[i])+(2*tau-1)*(data.ev[i])/2)
	}
	anser<-sum(ans)
	return(anser)
}




optim(c(-5,-5), opt.func)
optim(c(-5,0), opt.func)
optim(c(0,0), opt.func)
optim(c(0.1,0.1), opt.func)
optim(c(0.5,0.5), opt.func)
optim(c(1,1), opt.func)



tau<-0.95
optim(c(-7,0.3), opt.func.tau)
optim(c(0.1,0.1), opt.func.tau)
#-5.9734938+0.3300356*xxx

tau<-0.75
optim(c(-7,0.3), opt.func.tau)
# -6.5207481+0.3434336*xxx

tau<-0.25
optim(c(-7,0.3), opt.func.tau)
# -9.3203104+0.4261812*xxx

tau<-0.05
optim(c(-7,0.3), opt.func.tau)
#  -13.052722+0.542819*xxx


plot(data_x,data_y, xlab="WBGT",ylab="患者数",cex=.5, xlim=c(18, 35), ylim=c(0, 120))
par(new = T)
#95
plot(xxx,exp(-5.9734938+0.3300356*xxx),	xlim = c(18, 35), ylim = c(0, 120), xlab = "", ylab = "",type="l",col=4)
par(new = T)
#75
plot(xxx,exp(-6.5207481+0.3434336*xxx),	xlim = c(18, 35), ylim = c(0, 120), xlab = "", ylab = "",type="l",col=3)
par(new = T)
#25
plot(xxx,exp(-9.3203104+0.4261812*xxx),	xlim = c(18, 35), ylim = c(0, 120), xlab = "", ylab = "",type="l",col=3)
par(new = T)
#5
plot(xxx,exp(-13.052722+0.542819*xxx),	xlim = c(18, 35), ylim = c(0, 120), xlab = "", ylab = "",type="l",col=4)
legend("topleft", legend = c("95%","50%"), col = c(4,3), pch = c(1,1), lty = c(1,1))

出来ました．

５％分位点と２５％分位点回帰もそれなりにできていそうです．