毎回やる儀式

baseball <- read.csv("ball2017.csv",na.strings="*",fileEncoding = "UTF-8")

一般化線形モデル

前回は一般線形モデルというお話をしました。検定とモデリングって別物だと思ってた？でも一つの線形モデルで表現できるよね！ということに統計学者が気づいて，一般的につかえるモデルということで，一般線形モデルと呼んだのでした。

ところで，一般線形モデルは英語でなんというでしょう？ 正解は，General Liner Modelです。 それがどうしたんだ，と思われるかもしれませんが，今日のお話は一般化線形モデルです。一文字違う。この一文字の違いは，英語ではさらに微妙なGeneralized Liner Modelなのです。

これまた一般的に使えるモデルで，という意味なのですが，どこが「一般化」されたのか？それが今日の注目ポイントです。

線形モデルからみた回帰分析

前回は，2群を比較するときに，一方をベースライン\(b\)として，他方を\(b+d\)としましたね。ここで\(d\)は差分(difference)の\(d\)です。 3群以上を比較するときは，ベースライン\(b_0\)を用意して，各群の違いを\(b_0+d_j\)としたのでした。群ごとの相対的な違いを\(d_j\)で表現しました。 この形をよく見ると，\(y=ax+b\)に似ていませんか。こちらの場合だと，\(b\)がベースラインです。この式では切片(intercept)とも言います。\(ax\)の部分，この\(x\)が群の識別子で\(x=\{0,1\}\)だとすると，ある群\(x=0\)は\(y=b\)で，他方が\(y=a+b\)になっているので，\(a\)が差分を表していると考えられます。 この\(a\)は傾きと呼ばれたりしてました。

今日は，\(y=ax+b\)の回帰分析の考え方を，ベイジアン的に考えてみるとしましょう！

回帰分析

回帰分析のことを忘れた人は，第4回の講義を思い出してね！ こんなことをやったんでしたね。

library(ggplot2)
g <- ggplot(baseball,aes(x=height,y=weight)) + geom_point()
plot(g)

g <- g + geom_smooth(method="lm",se=FALSE)
plot(g)

reg <- lm(weight~height,data=baseball)
reg
## 
## Call:
## lm(formula = weight ~ height, data = baseball)
## 
## Coefficients:
## (Intercept)       height  
##    -160.906        1.377

ここでの基本的な考え方は，\(weight = \beta_0 + \beta_1 height\)，\(y=ax+b\)をちょっとカッコつけて書いてますが，要は一次関数を当てはめたということです。

このとき，誤差がありましたね。ちょっと見てみましょう。

res <- data.frame(reg$residuals)
summary(res)
##  reg.residuals     
##  Min.   :-14.0407  
##  1st Qu.: -5.9865  
##  Median : -0.4176  
##  Mean   :  0.0000  
##  3rd Qu.:  5.0788  
##  Max.   : 21.9593
g <- ggplot(res,aes(x=res))+geom_histogram(binwidth=1,alpha=0.8)+xlim(-25,25)
g

ここで見て欲しいポイントは，

  • 平均がゼロになってる
  • プロットが左右対称＝正規分布っぽい

というところです。

lm関数は最小二乗法による推定を行なっています。誤差の二乗和が最小になるように推定するのでそう呼ばれています。今回は誤差の二乗和がゼロになっていますね。すごい。

ところで，最小二乗推定は，推定法の一つだというお話をしました。統計モデルとしては，最尤法という推定方法もあります。これはモデルに確率分布の発想が入っています。誤差が平均ゼロの正規分布をするという考え方です。今回のデータはそれにあっているようです。

この考え方を使って，ベイジアンモデリングをしてみましょう。誤差が正規分布する，ということから，次のようにモデリングします。 LM

data{
  int<lower=0> N; //sample size
  real X[N];      //predictor
  real Y[N];      //predicted
}

parameters{
  real beta0;
  real beta1;
  real<lower=0> sigma;
}

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

  //prior
  beta0 ~ uniform(-10000,10000);
  beta1 ~ uniform(-10000,10000);
  sigma ~ cauchy(0,5);
}

上のコードをglm1.stanとして保存し，次のコードで推定してみましょう。

# ライブラリの読み込み
library(rstan)
## Loading required package: StanHeaders
## rstan (Version 2.16.2, packaged: 2017-07-03 09:24:58 UTC, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
# モデルのコンパイル
model <- stan_model("glm1.stan")

# データセットの準備
dat <- subset(baseball,select=c("height","weight"))
datastan <- list(N=nrow(dat),X=dat$height,Y=dat$weight)

# MCMCサンプリング(ベイズ推定)
fit <- sampling(model,datastan)
# 結果の表示
print(fit)
## Inference for Stan model: glm1.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##          mean se_mean    sd    2.5%     25%     50%     75%   97.5% n_eff
## beta0 -161.29    0.93 29.75 -221.67 -180.98 -161.21 -141.34 -103.12  1013
## beta1    1.38    0.01  0.16    1.06    1.27    1.38    1.49    1.71  1014
## sigma    8.35    0.02  0.68    7.14    7.89    8.32    8.76    9.82  1632
## lp__  -203.28    0.04  1.24 -206.53 -203.86 -202.97 -202.36 -201.87  1110
##       Rhat
## beta0 1.01
## beta1 1.01
## sigma 1.00
## lp__  1.00
## 
## Samples were drawn using NUTS(diag_e) at Thu Dec 14 16:08:42 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
# 最小二乗法の結果と比較
summary(lm(weight~height,data=dat))
## 
## Call:
## lm(formula = weight ~ height, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.0407  -5.9865  -0.4176   5.0788  21.9593 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -160.9058    30.0502  -5.355 8.81e-07 ***
## height         1.3768     0.1662   8.285 3.13e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.247 on 76 degrees of freedom
## Multiple R-squared:  0.4746, Adjusted R-squared:  0.4677 
## F-statistic: 68.65 on 1 and 76 DF,  p-value: 3.128e-12

結果の違いは，最小二乗法の結果と（ここではやってませんけど最尤法の結果も）同じような結果が得られていること，ただし，最小二乗法(や最尤法) は推定された切片や傾きの数字がただ一つ求められているのに対して，ベイズ推定の結果は確率分布している（幅を持った数字になっている）ことに注意してください。

いよいよ一般化，します。

さて，今日の話は一般化線形モデルでしたが，何を一般化しているかというと，「誤差が正規分布していないかもしれないじゃない！」という問題に対してモデルを一般化した，というのが答えです。

例えば，ホームラン数のデータを見て見ましょう。

hist(baseball$HR)

これは体重のような左右対称の分布ではなく，左に大きく偏ったデータです。これに普通の回帰分析をするのは，「正規分布を仮定する」というところが不自然になってしまいます。

実は，こういう非負のカウントデータは正規分布ではなくポアソン分布することがわかっています。ポアソン分布はパラメータが一つ，平均パラメータ\(\lambda\)だけを持つものなので，これを使ってモデリングすることを考えます。

par(mfrow = c(2, 2))
x <- 0:20
plot(dpois(x,1),type="b")
plot(dpois(x,3),type="b")
plot(dpois(x,5),type="b")
plot(dpois(x,10.5),type="b")

ここで，ラムダが体重によって予測されるとします。すなわち，\(\lambda = \beta_0 + \beta_1 weight\)です。ただし注意が必要なのは，パラメータ\(\lambda\)は正の数字でないといけないということです。左辺で足したりかけたりしてますので，これが万一，負の数が出るようなことがあれば，おかしなことになります。

そこで，数字を「必ず正の数になるように」整える必要があります。

この仕掛けは，指数関数を使うことで可能になります。

par(mfrow=c(1,1))
x <- seq(-10,10,0.1)
y <- exp(x)
plot(x,y)

モデルのコードは次のようになります。

data{
  int<lower=0> N;         //sample size
  real X[N];              //predictor
  int<lower=0> Y[N];      //predicted.正の整数に限定されていることに注意
}

parameters{
  real beta0;
  real beta1;
}

model{
  real lambda[N];
  //model
  for(n in 1:N){
    lambda[n] = exp(beta0 + beta1 * X[n]);
    Y[n] ~ poisson(lambda[n]);
  }

  //prior
  beta0 ~ uniform(-10000,10000);
  beta1 ~ uniform(-10000,10000);
}

上のコードをglm2.stanとして保存し，次のコードで推定してみましょう。

# モデルのコンパイル
model2 <- stan_model("glm2.stan")
## recompiling to avoid crashing R session
# データセットの準備
dat <- subset(baseball,select=c("HR","weight"))
datastan <- list(N=nrow(dat),X=dat$weight,Y=dat$HR)

# MCMCサンプリング(ベイズ推定)
fit2 <- sampling(model2,datastan)
# 結果の表示
print(fit2)
## Inference for Stan model: glm2.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##          mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff
## beta0   -1.12    0.01 0.23   -1.54   -1.28   -1.12   -0.96   -0.67   846
## beta1    0.04    0.00 0.00    0.04    0.04    0.04    0.04    0.04   850
## lp__  1645.08    0.03 0.95 1642.65 1644.68 1645.36 1645.78 1646.06  1000
##       Rhat
## beta0 1.01
## beta1 1.01
## lp__  1.00
## 
## Samples were drawn using NUTS(diag_e) at Thu Dec 14 16:09:43 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
# 最尤法の結果との比較
summary(glm(HR~weight,data=dat,family=poisson))
## 
## Call:
## glm(formula = HR ~ weight, family = poisson, data = dat)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -4.186  -2.279  -1.005   1.718   5.027  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.124494   0.230233  -4.884 1.04e-06 ***
## weight       0.040352   0.002427  16.624  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 685.40  on 77  degrees of freedom
## Residual deviance: 434.24  on 76  degrees of freedom
## AIC: 743
## 
## Number of Fisher Scoring iterations: 5

一応，最尤推定法のやり方も書いておきました。同じような結果がでていること，表現の仕方が異なっていることに注意してください。

一般化のポイント

一般化線形モデルは，このように「正規分布以外の分布にも対応する」というところがポイントです。 ただし，その際に，パラメータによっては「正の数字しかとらない」というような制限を加える必要があります。そんな時，パラメータに入れたい数式(\(\beta_0 + \beta_1 X\)など)を必要に応じて変換する(\(exp(\beta_0 + \beta_1 X)\)など)ことで，基本的には説明変数=\(\beta_0 + \beta_1 X\)で統一（一般化）できる，というのが一般化線形モデルの味噌です。

変換の式(リンク関数といいます)については色々あり得ますが，本講義では知って置いてほしいのは二つ，ポアソン分布のときのリンク関数とベルヌーイ分布のときのリンク関数だけです。

ベルヌーイ分布のモデル

ベルヌーイ分布は，コイントスのような「表か裏か」「生か死か」「男か女か」「正答か誤答か」のような2値の実現値を取る分布です。パラメータはこれも一つ，\(\theta\)だけです。

ベルヌーイ分布のパラメータは0から1の範囲しかとらないため，変換するリンク関数にはロジスティック関数を使います。このため，ロジスティック回帰分析と呼ばれることもあります。

ロジスティック関数は\(f(x)=1/(1+exp(-x))\)と表現されます。

x <- seq(-10,10,0.1)
y <- 1/(1+exp(-x))
plot(x,y)

今回はこれを使って，ホームランの本数を予測変数，所属リーグ(セ・リーグか，パ・リーグか)を従属変数にしたモデルを書いてみたいと思います。 モデルコードは次のようになります。

data{
  int<lower=0> N;         //sample size
  real X[N];              //predictor
  int<lower=0,upper=1> Y[N];      //predicted.0か1に限定されていることに注意
}

parameters{
  real beta0;
  real beta1;
}

model{
  real theta[N];
  //model
  for(n in 1:N){
    theta[n] = 1/(1+exp(-(beta0 + beta1 * X[n])));
    Y[n] ~ bernoulli(theta[n]);
  }

  //prior
  beta0 ~ uniform(-10000,10000);
  beta1 ~ uniform(-10000,10000);
}

上のコードをglm3.stanとして保存し，次のコードで推定してみましょう。

# モデルのコンパイル
model3 <- stan_model("glm3.stan")
## recompiling to avoid crashing R session
# データセットの準備
dat <- subset(baseball,select=c("league","HR"))
datastan <- list(N=nrow(dat),X=dat$HR,Y=as.numeric(dat$league)-1)

# データの確認
datastan
## $N
## [1] 78
## 
## $X
##  [1] 31 31 27 35  6 35 10 32 18 15 26 18 20 19 32 30 26 15  3  5  8 21 23
## [24]  6 14  3  5  3  9  0 14  4 15 24 17 25  6  2 24  9  3  1  9  1  2  0
## [47] 19  2 14  5  8  3  3  1  5  6 27  4 16  9  1 16  4  3  2  8  3 28  3
## [70] 12 18  6  1 31  2  9 26 11
## 
## $Y
##  [1] 1 1 0 0 1 1 1 0 0 0 1 0 1 1 1 0 1 0 1 1 1 0 0 0 0 1 1 1 1 0 1 0 0 0 0
## [36] 1 1 0 1 0 1 1 1 0 0 1 1 0 0 0 1 0 1 0 1 1 1 0 1 1 1 0 0 1 1 0 1 0 1 1
## [71] 0 0 1 1 0 0 0 1
# MCMCサンプリング(ベイズ推定)
fit3 <- sampling(model3,datastan)
# 結果の表示
print(fit3)
## Inference for Stan model: glm3.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## beta0   0.26    0.01 0.35  -0.43   0.02   0.26   0.49   0.95  1437    1
## beta1  -0.01    0.00 0.02  -0.05  -0.02  -0.01   0.01   0.03  1513    1
## lp__  -54.73    0.03 0.99 -57.48 -55.08 -54.40 -54.03 -53.78  1402    1
## 
## Samples were drawn using NUTS(diag_e) at Thu Dec 14 16:10:41 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
# 最尤法の結果との比較
summary(glm(league~HR,data=dat,family=binomial))
## 
## Call:
## glm(formula = league ~ HR, family = binomial, data = dat)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.292  -1.235   1.071   1.099   1.194  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.264383   0.360331   0.734    0.463
## HR          -0.008681   0.021968  -0.395    0.693
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 107.67  on 77  degrees of freedom
## Residual deviance: 107.51  on 76  degrees of freedom
## AIC: 111.51
## 
## Number of Fisher Scoring iterations: 3

様々な分布に対応できる，それが一般化線形モデリング。Enjoy!

本日の課題