【R stan】ベイズモデル その1：シンプルな線形回帰モデル

普段は、bayesmを使うことが良いのですが、stanを使ってみようということで、stanで色々なベイズモデルを作りました。

まずは、データの説明から。
Rに元々bayesmというパッケージについているcheeseというデータを使います。

library("bayesm")
library("rstan")
data(cheese)
> head(cheese)
              RETAILER VOLUME       DISP    PRICE
1  LOS ANGELES - LUCKY  21374 0.16200000 2.578460
2 LOS ANGELES - RALPHS   6427 0.12411257 3.727867
3   LOS ANGELES - VONS  17302 0.10200000 2.711421
4   CHICAGO - DOMINICK  13561 0.02759109 2.651206
5      CHICAGO - JEWEL  42774 0.09061273 1.986674
6       CHICAGO - OMNI   4498 0.00000000 2.386616

#### 通常の線形回帰モデル ####
model.lm <- lm(log(cheese$VOLUME)~cheese$DISP + log(cheese$PRICE), data=cheese)
summary(model.lm)
pred.reg <- predict(outreg,interval="prediction")
plot(log(cheese$VOLUME),pred.reg[,1],xlim=c(5,13),ylim=c(5,13),xlab="実測値",ylab="推定値")
abline(0,1)

RETAILERによって、回帰係数のばらつきが大きいので、通常の線形回帰モデルだと、VOLUMEの変動をうまく説明できていません。
決定係数は0.0982ととても小さいです。

次に、stanを使ったベイズモデルを作っていきます。
まずは、「ベイズモデル：個体差なしの線形回帰モデル」から。
「通常の線形回帰モデル」をベイズモデルに展開しただけなので、RETAILERの個体差を取り込めていないモデルとなります。

#### ベイズモデル：個体差なしの線形回帰モデル ####
data_stan <- list(N=nrow(cheese), DISP=cheese$DISP, PRICE=cheese$PRICE, VOLUME=log(cheese$VOLUME))
res_stan <- stan(file="./stan/cheese_model1.stan", data=data_stan, seed=123)

## heese_model1.stan ##
data{
  int N;
  real DISP[N];
  real PRICE[N];
  real VOLUME[N];
}

parameters{
  real b0;
  real b1;
  real b2;
  real sigma;
}

model{
  for(n in 1:N){
    VOLUME[n] ~ normal(b0+b1*DISP[n]+b2*PRICE[n], sigma);
  }
  // Leave out prior distribution. Uninformative prior will be used. 
}

得られた係数は、「通常の線形回帰モデル」とほぼ同じ。
まぁ、最尤法で解くかベイズで解くかの違いくらいです。

決定係数、RMSE、MAEなどの精度ももほとんど変わりません。
model 1: ベイズモデル：個体差なしの線形回帰モデル
　RMSE: 6377
　MAE: 2818
でした。

やはり、個体差（RETAILERによる違い）をモデルに取り込まないとダメみたいですね。