RStanで単回帰を試してみる - Gakushukun1’s diary

目的

estatsの統計データから, rstanによる単回帰を用いて各都道府県の死亡者数と出生児数の関係を調べる. 事後分布からサンプリングされたパラメータから, 情報量規準WAICと交差検証LOOを計算し, 単回帰と二次回帰どちらが真の分布に近いかを推測する.

今回はrstanを使って単回帰分析を試してみる. 今回対象としたのは, estatの中の, 平成29年10月〜平成30年9月間の都道府県別の出生児数・死亡者数のデータである.
このデータについて, 死亡者数を横軸, 出生児数を縦軸で表すと次のようなグラフとなる. それぞれの点が各都道府県のデータを表している.
f:id:gakushukun1:20190528203856p:plain

このデータから, 死亡者数から出生児数への回帰曲線を推測する. 死亡者数を $X$ , 出生児数を $Y$ とし, 次のような2つのモデルを考える.

$Y = a + bX + \epsilon \\ Y = a + bX + cX^2 + \epsilon$

それぞれのモデルについてstanでパラメータの事後分布を作り, サンプリングを行った. 元の確率モデルに対して. サンプリングで得られたパラメータを当てはめた数値を保存しておくことで, 擬似的な予測分布を作り, waicとlooの計算を行った. 実際のコードは次の通り.

線形回帰(R)

library(rstan)
library(loo)

d <- read.csv("FEH_00200524_190526235535.csv", stringsAsFactors=FALSE)
d[,2] <- d[,2] / 10000
d[,3] <- d[,3] / 10000
x <- d$死亡者数
y <- d$出生児数
data <- list(N=nrow(d), X=x, Y=y)
fit <- stan(file="kaiki.stan", data=data, seed=123)
log_lik <- extract_log_lik(fit)
ms <- extract(fit)
a <- mean(ms$a)
b <- mean(ms$b)
sigma <- mean(ms$sigma)
waic01 <- waic(log_lik)
loo01 <- loo(log_lik)
png("kaiki1.png", width=800, height=800)
par(family = "HiraKakuProN-W3")
plot(x, y, xlab="死亡者数", ylab="出生児数", xlim=c(0,13), ylim=c(0,13))
abline(a, b)
abline(a + sigma * 1.96, b)
abline(a - sigma * 1.96, b)
dev.off()

線形回帰(Stan) "kaiki.stan"

data {
    int N;
    real X[N];
    real Y[N];
}

parameters {
    real a;
    real b;
    real<lower=0> sigma;
}

model {
    for (n in 1:N) {
        Y[n] ~ normal(a + b * X[n], sigma);
    }
}

generated quantities {
    vector[N] log_lik;

    for (n in 1:N) {
        log_lik[n] = normal_lpdf(Y[n] | a + b * X[n], sigma);
    }
}

二次回帰(R)

library(rstan)
library(loo)

d <- read.csv("FEH_00200524_190526235535.csv", stringsAsFactors=FALSE)
d[,2] <- d[,2] / 10000
d[,3] <- d[,3] / 10000
x <- d$死亡者数
y <- d$出生児数
data <- list(N=nrow(d), X=x, Y=y)
fit <- stan(file="kaiki2.stan", data=data, seed=123)
log_lik <- extract_log_lik(fit)
ms <- extract(fit)
a <- mean(ms$a)
b <- mean(ms$b)
c <- mean(ms$c)
sigma <- mean(ms$sigma)
fun <- function(x) { a + b * x + c * (x ^ 2) }
fun_upr <- function(x) { a + b * x + c * (x ^ 2) + sigma * 1.96 }
fun_lwr <- function(x) { a + b * x + c * (x ^ 2) - sigma * 1.96 }
waic01 <- waic(log_lik)
loo01 <- loo(log_lik)

png("kaiki2.png", width=800, height=800)
par(family = "HiraKakuProN-W3")
plot(x, y, xlab="死亡者数", ylab="出生児数", xlim=c(0, 13), ylim=c(0, 13))
par(new=T)
plot(fun, xlab="", ylab="", xlim=c(0, 13), ylim=c(0, 13))
par(new=T)
plot(fun_upr, xlab="", ylab="", xlim=c(0, 13), ylim=c(0, 13))
par(new=T)
plot(fun_lwr, xlab="", ylab="", xlim=c(0, 13), ylim=c(0, 13))
dev.off()

二次回帰(Stan)

data {
    int N;
    real X[N];
    real Y[N];
}

parameters {
    real a;
    real b;
    real c;
    real<lower=0> sigma;
}

model {
    for (n in 1:N) {
        Y[n] ~ normal(a + b * X[n] + c * (X[n] ^ 2), sigma);
    }
}

generated quantities {
    vector[N] log_lik;

    for (n in 1:N) {
        log_lik[n] = normal_lpdf(Y[n] | a + b * X[n] + c * (X[n] ^ 2), sigma);
    }
}

元のデータの数値が大きいため, 線形, 二次回帰いずれの場合も全てのデータを10000分の1倍している. 予測分布のYについての平均を計算することで, 回帰曲線を推測し, サンプリングされたノイズの標準偏差 $\sigma$ の1.96倍の値を足し引きすることで, 与えられた $x$ に対する $y$ の95%のベイズ予測区間を示した.