2018-03-20

RでGLMあれこれ（ポアソン回帰）

統計学

glm()関数を使ってモデリング

データセット：http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv

0. 統計モデリングの流れ

1. データの確認とモデルの選択
　要約統計量とグラフから、どの統計モデルを使うべきか決める。
2. 統計モデルの作成：あてはめ（fitting）
　GLMであれば、「説明変数」は、いくつかの「母数」によって説明される。
　「母数」は、リンク関数（「母数」の関数）と線形予測子（「応答変数」とその「パラメータ」によ
　　る線形結合）の関係によって推定される。
　　つまり、リンク関数・線形予測子を決定することで、モデリングが可能に。
3. 統計モデルの検証
　3-1. モデルによって「説明」が可能か？（あてはまりの良さ）
　3-2. モデルによって「予測」が可能か？（偶然性の小ささ）

1. データの確認とモデルの選択

> d   <- read.csv("data3a.csv")#ファイル読み込み
> head(d)#データフレームの確認
   y     x f
1  6  8.31 C
2  6  9.44 C
3  6  9.50 C
4 12  9.07 C
5 10 10.16 C
6  4  8.32 C
> summary(d)#f列はCとTの2水準で構成
       y               x          f     
 Min.   : 2.00   Min.   : 7.190   C:50  
 1st Qu.: 6.00   1st Qu.: 9.428   T:50  
 Median : 8.00   Median :10.155         
 Mean   : 7.83   Mean   :10.089         
 3rd Qu.:10.00   3rd Qu.:10.685         
 Max.   :15.00   Max.   :12.400  
> class(d$y)#yはinteger型
[1] "integer"
> class(d$x)#xはinteger型
[1] "numeric"
> class(d$f)#fはfactor型
[1] "factor"
> plot(d$x, d$y, pch = c(21, 19)[d$f])#CとTで色分けして散布図
> legend("topleft", legend = c("C", "T"), pch = c(21, 19))#CとTの凡例を追加

f:id:yumaloop:20180320232343p:plain — 横軸にx列、縦軸にy列をとった散布図

　ポアソン分布に似ている。→ポアソン回帰でGLMが使えそう

この散布図に置いて、種子数yを説明変数、体サイズxと施肥処理fを応答変数として、統計モデル（ここではポアソン回帰）を作りたい。とりあえず、いくつかモデルを作りそれぞれの尤度をさぐる。

2. 統計モデルのあてはめ（fitting）

それぞれの仮定した統計モデルに対して、glm()関数で得られたパラメータ（ $\beta_1, \beta_2, \beta_3, ...$ ）の推定値により、母数（ $\lambda_i,$ ）を推定し、さらに得られた母数（ $\lambda_i$ ）から尤度（（ $L(\lambda_i)$ ））を求める。

※パラメータ（ $\beta_1, \beta_2, \beta_3, ...$ ）の推定値を求める過程
・ポアソン分布 $Po(\lambda)$ に従う確率変数 $x$ の確率密度関数 $f(x)$

　　 $f(x)={\large\frac{\lambda^{x}exp(-\lambda)}{x!}} {\scriptsize (x = 0, 1, 2, ...)}$

ある個体 $i$ において、種子数が $y_i$ である確率 $p(y_i | \lambda_i)$ はポアソン分布 $Po(\lambda_i)$ に従っていると仮定すると、以下の式が成り立つ。（ポアソン回帰の一般式）

　　 $p(y_i | \lambda_i)={\large\frac{\lambda_{i}^{y_i}exp(-\lambda_{i})}{y_i!}} {\scriptsize (y_i = 1, 2, 3, ...)}$

この式から、パラメータ（ $\beta_1, \beta_2, \beta_3, ...$ ）推定値は、ポアソン回帰においては以下の式によって導かれる。

　　 $\log{L(\beta_1, \beta_2, \beta_3, ...)} = \sum_{i=1}^{n} \log{\large{\frac{\lambda_{i}^{y_i}exp(-\lambda_{i})}{y_i!}}}$

すなわち、対数尤度関数 $\log{L(\beta_1, \beta_2, \beta_3, ...)}$ を最大にするようなパラメータ（ $\beta_1, \beta_2, \beta_3, ...$ ）の値を求めれば良い。

以下で異なる３種類のポアソン回帰モデルを立てて、推定結果を比較する。

（モデルA）「yはxによって説明される。」
　　　 ※各個体 $i$ について種子数 $x_i$ は、integer型で1以上の整数。
　　 $\lambda_i=exp(\beta_1+\beta_2x_i)$
　　i.e. $\log{\lambda_i}=\beta_1+\beta_2x_i$

> fit <- glm(y ~ x, data = d, family = poisson)#glm()関数でモデルの最尤推定値を得る
> fit#fitには推定結果が格納されている

Call:  glm(formula = y ~ x, family = poisson, data = d)

Coefficients:
(Intercept)            x  
    1.29172      0.07566  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    89.51 
Residual Deviance: 84.99 	AIC: 474.8

（モデルB）「yはfによって説明される。」
　　　 ※各個体 $i$ について施肥処理 $f_i$ は、factor型でC or Fをとる。
　　　 ※ $f_i$ のダミー変数 $d_i$ は、factor型で0 or 1をとる。

　　 $\lambda_i=exp(\beta_1+\beta_3d_i)$
　　i.e. $\log{\lambda_i}=\beta_1+\beta_3d_i$

> fit <- glm(y ~ x, data = d, family = poisson)
> fit

Call:  glm(formula = y ~ x, family = poisson, data = d)

Coefficients:
(Intercept)            x  
    1.29172      0.07566  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    89.51 
Residual Deviance: 84.99 	AIC: 474.8

・推定結果
β1の推定値 : 1.29172
β2の推定値 : 0.07566

・

（モデルC）「yはxとfによって説明される。」

*** 参考

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (29件) を見る

2018-03-20

【書評】『統計学入門』（東京大学出版）～統計学入門書の金字塔～

統計学書評

f:id:yumaloop:20180320170554j:plain

2017年の年末くらいから、統計学を本格的に学んでみようと思っていたので、まず読み始めたのがこの本です。Googleで「統計学本」で検索すると、絶対にこの本がヒットするくらい有名。

以下書評です。↓(´・ω・`)

1. 書籍情報

細かな情報

f:id:yumaloop:20180320010625j:plain

・書名　：統計学入門

・著者名：東京大学教養学部統計学教室

・出版社：東京大学出版会　　　　　　

・種別　：単行本, 307ページ　　　　

・発売日： 1991/7/9　　　　　　　　　

・予算　： 3024円 (※Amazonでの価格)

※Amazonのリンクはこちら

姉妹本として、おもに理系分野よりの『自然科学の統計学』、文系分野よりの『人文・社会科学の統計学』の2冊が出版されています。こちらは、より発展的な内容。

時代背景　－90年代とコンピュータの普及

　この本の初版は1991年に出ています。1991年というのは、InternetやPCが本格的に普及される以前で、R・Pythonなどの統計解析プログラムはもちろん世に出ていません。*5しかし、著者たちは、コンピュータの可能性をすでに鋭く察知しています。実際、序文の最後では、プログラミングパッケージの発展を歓迎しつつも、解析過程がブラックボックス化することを憂いた以下のような言及があります。

”しかし、方法の意味がわからずに、ただ計算をしてもその結果を正しく利用することはできないし、また場合によっては誤まった判断をくだすことになる危険がある。数理統計学の方法を形式的に説明するだけではなく、その意味を把握するようにすること、それがこの教科書の最も重要な目標である。”

著者たちが、

・「ITの発展によって大規模データの取得とその処理が可能になること」

・「統計学の適用範囲が広がり、学習への興味と重要性が増すこと」

といった時代の要請に応える形で、この本を執筆したという経緯がわかります。

2. 内容

章立てと構成

　章立てを見るとわかる通り、過不足ない網羅性に加え、各論が綺麗に分割されています。文体はとても流暢で、ストレスなく読み進められます。

・第1章　統計学の基礎
・第2章　 1次元のデータ
・第3章　 2次元のデータ
・第4章　確率
・第5章　確率変数
・第6章　確率分布
・第7章　多次元の確率分布
・第8章　大数の法則と中心極限定理
・第9章　標本分布
・第10章    正規分布からの標本
・第11章    推定
・第12章    仮説検定
・第13章    回帰分析