Help us understand the problem. What is going on with this article?

Stanで対数正規分布モデリング ーAirqualityで遊ぶ③ー

Rのサンプルデータ 「Airquality」で遊ぶ③です。
前回と前々回ではガンマ分布のGLM及びベイズモデリングを考えました。

前記事:
ガンマ分布GLMで重回帰分析 予測区間の導出とモデルチェック ーAirqualityで遊ぶ①ー
Stanでガンマ分布モデリング ーAirqualityで遊ぶ②ー
次記事:
StanでWAICとWBICの導出とモデル比較 BridgeSamplingの疑問

Ozone濃度が正の値で、線形モデルだとあまりよろしくないので、ガンマ分布を考えていました。第3部では、「正の連続値」への対処として、ガンマ分布ではなく対数正規分布を用いたモデリングを試みます。第2部と同様、RとStanコードを記します。

Rコード↓

library(rstan)
d <- airquality
d <- drop_na(d)
Ozone <- as.data.frame(d[,1])
colnames(Ozone) <- "Ozone"
X <- as.matrix(d[,2:4])
X <- cbind(1, X)

data <- list(N=111, X=X, Ozone=Ozone$Ozone)
fit <- stan(file="lognormal.stan", data=data, seed=1234)
fit

Stanコード↓(第2部のガンマ分布よりよっぽど簡単ですね)

data{
  int N;
  matrix[N,4] X;
  vector[N] Ozone;
}

parameters{
  vector[4] b;
  real<lower=0> sigma;
}

transformed parameters{
  vector[N] mu;
  mu = X*b;
}

model{
  Ozone ~ lognormal(mu, sigma);
}

generated quantities{
  vector[N] Y_pred;
  for(n in 1:N)
    Y_pred[n] = lognormal_rng(mu[n], sigma);
}

実測値-予測値プロットを行います。このコードは第2部と全く同じです。

ms <- rstan::extract(fit)
y_pred <- ms$Y_pred 
Probs <- c(0.05, 0.25, 0.5, 0.75, 0.95)
y_qua <- apply(y_pred, 2, quantile, probs=Probs) 
Y_qua <- data.frame(d, t(y_qua))
Y_qua <- rename(Y_qua, p5=X5., p25=X25., p50=X50., p75=X75., p95=X95.)

Y_qua$Month <- as.factor(Y_qua$Month)

ggplot(Y_qua, aes(x=Ozone, y=p50, ymin=p5, ymax=p95, colour=Month))+
  theme_bw()+
  geom_pointrange()+
  geom_abline(aes(slope=1, intercept=0), linetype="dashed")+
  coord_fixed(ratio=1, xlim=c(0,400), ylim=c(0,400))+
  labs(x="Observed", y="Predicted")

lognormal_bayes.png

第2部のガンマモデリングと比較してみると
comparison.png

また、対数正規分布の予測の正答率(第1部参照)は以下のようになりました。

nrow(dplyr::filter(Y_qua, Ozone > p5 & Ozone < p95))  # 予測率102/111 = 91.9% (90%理想)
nrow(dplyr::filter(Y_qua, Ozone > p25 & Ozone < p75))  # 予測率60/111 = 54.1% (50%理想)

どちらのモデルが「より精度が高い」のでしょうか?

見た目からではほとんど分からないですね。今回は、「正答率」が理想に近いのがガンマ分布だったので、ひとまずそちらを採用する、という結論にしたいと思います。

しかし、もう少し詳しく調べる必要がありそうです。

余裕があれば、WAICやWBICといった指標から、今回のモデル比較を行ってみたいと思います。

とにかく、現時点でなんとなく記事にしたかった部分は一通り書けたので、これで終わりとします。
ありがとうございました。

前記事:
ガンマ分布GLMで重回帰分析 予測区間の導出とモデルチェック ーAirqualityで遊ぶ①ー
Stanでガンマ分布モデリング ーAirqualityで遊ぶ②ー
次記事:
StanでWAICとWBICの導出とモデル比較 BridgeSamplingの疑問

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした