割合の従属変数をベイズのベータ回帰と
ベータ分布×双極関数のフィッティングで
author: Mr.Unadon (見習い飯炊き兵)
動作環境:Mac OS Sierra 10.12.1; R version３.３.2; rstan 2.10.1
はじめに
０から１の間にあるパーセントで表記される割合データに対する回帰系のモデルのお話です。
本稿では、従属変数に対してベータ分布を仮定したベータ回帰のベイズモデルと、
双極関数を平均とするベータ分布でデータを予測させるベータ分布×双極関数モデルをご紹介します。
従属変数は、「都道府県別の人口に占める高齢者福祉施設在所者数の割合」。
説明変数には、「都道府県別の可住地面積に占める人口の割合」、すなわち人口密度的なものを入れました。
都会って、他の人との物理的な距離は近いのに、心の距離は遠いと感じられますよね。
そういう直感から、何かが予測できるんじゃないかと思っての試みです。
パッケージとデータの読み込み
都道府県別の人口データはこちらから。2014年のものです。
都道府県別の可住地面積はこちらから。
高齢者福祉施設の在所者数データはこちらから。
以上のデータを読み込んで整理します。人口は1000人単位なので1000を乗算しておきます。 可住地面積の単位は100^2ｋｍです。
invisible({rm(list=ls());gc();gc()})
##########################################################
# Beta Regression and Beta Hyperbolic fitting
##########################################################
#一時保存データをクリア
library(tidyverse)
library(plotly)
library(rstan)
library(shinystan)
library(formattable)
#stan parallel
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
#read data
dat<-read.csv("dat.csv",fileEncoding = "shift-jis")
#data
d<-dat%>%dplyr::mutate(人口=人口*1000,人口密度=人口/可住地面積,高齢者福祉施設在所者率=介護老人福祉施設在所者数/人口)%>%
dplyr::select(都道府県,人口, 人口密度, 高齢者福祉施設在所者率)
#table output
formattable(d)
ローデータの可視化
都道府県別で、X軸に人口密度(本当の人口密度は面積で割ったものを指しますが)をとり、 Y軸に福祉施設在所率を取ります。
#ローデータの可視化
ggplot()+theme_set(theme_light(base_family = "Hiragino Mincho Pro W6"))
gg<-ggplot(d,aes(x=人口密度, y=高齢者福祉施設在所者率))+
geom_point(aes(colour=都道府県),size=3)+guides(colour=FALSE)+
theme(plot.subtitle = element_text(vjust = 1),
plot.caption = element_text(vjust = 1),
axis.title = element_text(size = 18),
axis.text = element_text(size = 14),
plot.title = element_text(size = 18)) +labs(x = "人口密度(人口/可住地面積)")
ggplotly(gg)
曲線のきれいなトレンドが見て取れます。
さて、これをまずはベータ回帰でモデリングしてみます。
ベータ回帰のベイジアンモデリング
ベータ回帰とは、従属変数にベータ分布(0から１の間の分布)を仮定した非線形回帰モデルのことを指します。
正規分布の様に左右対象ではないので、「線が曲がります」。
直線の線形モデルの予測値(mu)にリンク関数(下記ではinv_logit())をかませまして、ベータ分布の平均を通るようにさせます。
ベータ分布の平均はa/(a+b)です。今回はthetaで表現します。
Stanコードはこちら。
data {
int N;
vector [N] rate;
matrix [N,2]X;
}
parameters {
real <lower=0>phi;
vector [2] beta;
}
transformed parameters {
vector [N] theta;
vector [N] mu;
vector <lower=0> [N] a;
vector <lower=0> [N] b;
mu = X*beta;
for(i in 1: N){
theta[i] = inv_logit(mu[i]);
}
a = theta * phi;
b = (1-theta) * phi;
}
model {
rate ~ beta(a,b);
}
generated quantities{
vector [N]pred;
for(i in 1:N){
pred[i] = beta_rng(a[i],b[i]);
}
}
解析を実行して可視化まで持っていきます。
#data to be passed on stan
rate=c(d$高齢者福祉施設在所者率)
N=length(d$都道府県)
X=matrix(c(rep(1,N),scale(c(d$人口密度))),nrow = N)
#data to be passed on stan
datastan=list(N=N,rate=rate,X=X)
parameters<-c("theta","beta","phi","pred")
#from .stan file
fit<-stan(file="BetaRegression.stan",data=datastan,pars=parameters)
#-----------------------------------------
#stan results
theta<-data.frame(rstan::extract(fit)$pred)
#data frame
lo<-c()
up<-c()
lo50<-c()
up50<-c()
M<-c()
for(i in 1: N){
lo[i]<-quantile(theta[,i],0.975)
up[i]<-quantile(theta[,i],0.025)
lo50[i]<-quantile(theta[,i],0.75)
up50[i]<-quantile(theta[,i],0.25)
M[i]<-mean(theta[,i])
}
d_plot1<-d%>%dplyr::mutate(lo=lo,up=up,lo50=lo50,up50=up50,M=M)
ggPred<-ggplot(d_plot1,aes(x=scalePop, y=高齢者福祉施設在所者率))+
geom_ribbon(ymin=lo,ymax=up,alpha=0.2)+
geom_ribbon(ymin=lo50,ymax=up50,alpha=0.4)+
geom_point(aes(colour=都道府県),size=3)+guides(colour=FALSE)+
geom_line(aes(x=scalePop, y=M))+
theme(plot.subtitle = element_text(vjust = 1),
plot.caption = element_text(vjust = 1),
axis.title = element_text(size = 18),
axis.text = element_text(size = 14),
plot.title = element_text(size = 18)) +labs(x = "人口密度(人口/可住地面積)")
ggplotly(ggPred)
ベータ分布の平均を通る双極関数モデルのフィッティング
上記の結果でも良いのかもしれませんが、見た感じ「最初急激に下がって、後は横ばい」なトレンドを表現できているとは言いにくいところです。
そこで、急激な下降からの横ばいトレンドを持つ双極関数をこのデータにフィットさせてみることにします。
なお、双極関数は０に極限するためそのまま利用することは適当ではありません。
双極関数の底に相当する部分として、パラメータalphaを加算することでそれをクリアします。
gammaは減衰係数とでも呼べばよいでしょうか。凹み具合を決めるパラメータです。
ベータ回帰の場合の様にリンク関数は挟まず、そのまま平均を通るようにコードを描いていきます。
Stanコードは以下。
data {
int N;
vector [N] rate;
vector [N] X2;
}
parameters {
real <lower=0,upper=1>alpha;
real <lower=0,upper=1>gamma;
real<lower=0> phi;
}
transformed parameters {
vector [N] theta;
vector <lower=0> [N] a;
vector <lower=0> [N] b;
for(i in 1: N){
theta[i] =alpha+1/(1+(X2[i]*gamma));
a[i] = theta[i] * phi;
b[i] = (1-theta[i]) * phi;
}
}
model {
rate ~ beta(a,b);
}
generated quantities{
vector [N]pred;
for(i in 1:N){
pred[i] = beta_rng(a[i],b[i]);
}
}
では、推定を実行します。 可視化まで一気に参ります。
#
# Hyperbolic & beta distribution fit
#
#data to be passed on stan
N=length(d$都道府県)
X2=c(d$人口密度)
rate=c(d$高齢者福祉施設在所者率)
#data to be passed on stan
datastan=list(N=N,rate=rate,X2=X2)
parameters<-c("theta","phi","gamma","alpha","pred")
#from .stan file
fit2<-stan(file="hyperbolic.stan",data=datastan,pars=parameters)
#Gelman-Rubin test
#launch_shinystan(fit)
#-----------------------------------------
#stan results
theta2<-data.frame(rstan::extract(fit2)$pred)
#data frame
hlo<-c()
hup<-c()
hlo50<-c()
hup50<-c()
hM<-c()
for(i in 1: N){
hlo[i]<-quantile(theta2[,i],0.975)
hup[i]<-quantile(theta2[,i],0.025)
hlo50[i]<-quantile(theta2[,i],0.75)
hup50[i]<-quantile(theta2[,i],0.25)
hM[i]<-mean(theta2[,i])
}
# data frame
d_plot2<-d%>%dplyr::mutate(hlo=hlo,hup=hup,hlo50=hlo50,hup50=hup50,hM=hM)
#plot
ggPred2<-ggplot(d_plot2,aes(x=人口密度, y=高齢者福祉施設在所者率))+
geom_ribbon(ymin=hlo,ymax=hup,alpha=0.2)+
geom_ribbon(ymin=hlo50,ymax=hup50,alpha=0.4)+
geom_point(aes(colour=都道府県),size=3)+guides(colour=FALSE)+
geom_line(aes(x=人口密度, y=hM))+
theme(plot.subtitle = element_text(vjust = 1),
plot.caption = element_text(vjust = 1),
axis.title = element_text(size = 18),
axis.text = element_text(size = 14),
plot.title = element_text(size = 18)) +labs(x = "人口密度(人口/可住地面積)")
ggplotly(ggPred2)
最後に
割合データを従属変数とする非線形回帰にベータ回帰を使用してみました。
正規分布を仮定した通常の線形回帰やガンマ回帰を使うと、「割合なのに０や１を突き抜ける予測」になってしまいますので、ベータ回帰はとても理にかなっています。
ただ今回はそれだけでは不十分に見えましたので、双極関数によるフィットも行いました。 「高齢者福祉施設在所者率には、必ず”これ以上は下がらないラインがある”」といった根拠が事前にあるのであれば、
今回のようなトレンドデータに積極的に適用していくのも良いと思っています。みなさんも是非、色んなモデリングにチャレンジしてみてください。
私も勉強します。
んー、でもしかし、人口密度が上がるほど高齢者福祉施設在所者率がこのトレンドになるのはどうしてなのだろうなぁ…。
私がGWに考えてきたことのアウトプットは、ここまで。
平常業務、頑張るぞ。