1. No category

2015/3/27 12:00-13:00
日本草地学会 若手R統計企画 （信州大学農学部）
Rと一般化線形モデル入門
山梨県富士山科学研究所
安田 泰輔
謝辞：日本草地学会若手の会の皆様、発表の機会を頂き、
たいへんありがとうございます！
自己紹介
• 茨城大学・学生時代
– ベータ二項分布を用いた種の空間分布の解析
• 所属：山梨県富士山科学研究所
• 最近の研究テーマ
– 近接リモートセンシングによる半自然草地のモニタリング手法開発
– 富士山の長期的、広域的な森林限界動態
– 湖畔における外来植物アレチウリの分布推定と駆除活動
趣旨説明
• Rに触れてみよう！、GLMをやってみよう！ざっくばらんに！（ということだ
と思う）
文章は若手の会作成
R（主にwikipediaより）
•
•
•
•
オープンソース・フリーソフト
R Development Core Team（NZ）がメンテナンス
統計解析向けのプログラミング言語
統計解析、画像処理、音声合成、GISなど各種多様なパッケージが利用
できる（6000以上！だそうです）。もちろんフリー。
• Rをダウンロードしたとき、全部は入っていない
• ので、必要なパッケージ（様々な関数群）をダウンロードして使う。
E.g. 画像処理がしたい！→調べる、聞く
> install.packages(“raster”)
> library(raster)
で使えるようになる。
R
• グラフィック、他のソフトとの連携など、統計以外のパッケージも充実
• ファイル操作もできます。
• 日々、拡張と強化が行われているので、情報交換がとても大事
Viewer（google map, google earth）
Analysis
Editor（R-Studio, notepad++）
GIS（SAGA, GLASS GIS, QGIS）
テキストベースに慣れよう！
• “何度もボタンを押す” → “テキストで記述”
– 処理内容を記録できる、何度も使える、一括処理できる、などの利点がある。
Notepad++
今回使用するソフト：R-Studio
ここにコードを記載
Rの実行
変数情報等
図やパッケージ情報等
少しやってみます。
RでGLM
• 一般化線形モデル Generalized Linear Model
• フィールド研究ではよく使われる方法の１つ
• 良書、参考HPは非常にたくさんあります。
– データ解析のための統計モデリング入門（久保拓弥、岩波書店）
– 講義のーと（データ解析のための統計モデリング）
http://eprints.lib.hokudai.ac.jp/dspace/handle/2115/49477
– 名古屋大学大学院生命環境農学研究科 森林生態生理学研究分野
RでGLMをやってみよう
http://www.agr.nagoya-u.ac.jp/~seitai/document/R2009/t090614glmIntro.pdf
などなど
他の手法との関連
• これらはGLMの特殊なバージョン
• 二項分布やポアソン分布、ガンマ分布も使える。
• 説明変数には、連続変数、カテゴリ変数、両方が取れる。
分析名
応答変数（Y）
説明変数（X）
t検定（Welchでない場合）
量的変数１つ
2分変数１つ
一元配置分散分析
量的変数１つ
カテゴリ変数１つ
多元配置分散分析
量的変数１つ
カテゴリ変数複数
回帰分析
量的変数１つ
量的変数１つ
重回帰分析
量的変数１つ
量的変数複数（カテゴリ変数はダミー変数化
共分散分析
量的変数１つ
・・・
ロジスティック回帰分析
2分変数１つ
2分変数、カテゴリ変数、量的変数複数
正準相関分析
量的変数複数
量的変数複数
中澤（2004） “Rによる統計解析の基礎” を改編・引用
内容
ロジスティック回帰を中心に、
1. プロット
2. モデルの組み立て、パラメータ推定
3. AICによるモデル選択
を行ってみる。 ＊数式は重要ですが、ここでは極力使いません。
問題設定とデータセットの作成
• 土壌水分条件に対する植物種Sの種子生残率を知りたい
– 水分条件が異なる場所12か所を選ぶ
– 各地点で種Sの種子を25～50個採取し、生残or死んでいるかを調査した。
種子サンプリング
乾燥
1
生きている：１
死んでいる：0
2
・・・
12
湿潤
x:土壌含水率
n:各地点で調査した総種子数
y：生残種子数
1.プロット
• 説明変数に伴って変化する確率分布の期待値（赤）を推定。
“この確率分布はなにか？”
“説明変数と期待値はどのような関係か？”
2.モデル組み立て準備
• 生残率は 生残種子数yi/総種子数ni （i = 1,2,…,12）。
• 生残種子数yiは0,1,・・・,nの値を取り得る。
• 生残率（n個中ｙ個当たり）は、二項分布に従うと仮定。
• 二項分布Binomial(pi, ni)の期待値piは説明変数xiとなんらかの
関係を持っていると仮定。
• これらの条件・仮定から、モデルを組み立てる。
モデル組み立て
1. 確率分布： yi ~ Binomial(pi, ni)
i=1,2,…,12
2. リンク関数： logit(pi) = log(pi/(1-pi))=線形予測子
pi=1/(1+exp(-(線形予測子)))、 e.g. exp(3) = e^3
3. 線形予測子：β0 + β1*xi
• データの性質に合わせて、確率分布を仮定する。
• 説明変数は、連続変数（今回）、カテゴリ変数、その両方が使える。
いろんなデータ型へも適用可能
• GLMを行うときは
応答変数の確率分布、リンク関数、線形予測子
を指定する。
• 確率分布（指数型分布族 ）とリンク関数（デフォルト）>？glm → family
データ
確率分布
リンク関数
（デフォルト）
0,1,…,n（離散、上下限あり）
二項分布 binomial
logit
0,1,…,∞（離散、上限なし）
ポアソン分布 poisson
log
0 ~ ∞（連続、上限なし）
ガンマ分布 Gamma
Inverse (or log)
-∞ ~ ∞（連続、上下限なし）
正規分布 gaussian
identity （そのまま）
RでGLM実行（パラメータ推定）
>
fit<-glm(
cbind(y,n-y)~x, #cbindで成功数・失敗数を指定
family=binomial(link=“logit”),
data=data_set
)
注意
cbind(y,n-y)~x の右辺は “切片＋傾き*x” の意味
結果
logit(pi) =線形予測子
pi=1/(1+exp(-(線形予測子)))
Residual deviance/df = 8.897/10=0.8897
結果の図示
β0=-10.7, β1=0.25
β0^=-9.53(s.e. 0.84), β1^=0.22(s.e. 0.02)
3.まだ報告書は終わらない・・・AICによるモデル選択
β0=-10.7, β1=0.25
β0^=-9.53(s.e. 0.84), β1^=0.22(s.e. 0.02)
線形予測子
β0+β1*xi
線形予測子
β0
β1=0(傾き=0)と違うのか？
有意差検定をする
↓観方を変える
どのモデルが“良い”のか？
AICでモデルを選択する。
AIC
• 赤池情報量規準 Akaike’s Information Criterion
AIC = -2lnL + 2k
（lnL: 最大対数尤度、k：パラメータ数）
• 同じデータのもとで、2つ以上のモデルのAICを比較し、AICが最小と
なるモデルを採用する。(山村光司さんのHP http://cse.niaes.affrc.go.jp/yamamura/)
• 相対的にAICが低いモデルが“良い”モデルとされる。
再：RでGLM実行（パラメータ推定）
>
fit<-glm(
cbind(y,n-y)~1, #cbindで成功数・失敗数を指定
family=binomial(link=“logit”),
data=data_set
)
注意
cbind(y,n-y)~1 の右辺は “切片だけ” の意味
“切片だけ”のモデル
AICの比較
β0=-10.7, β1=0.25
fit β0^=-9.53(s.e. 0.84), β1^=0.22(s.e. 0.02)
fit0 β0^=-0.29(s.e. 0.09)
AICによるモデル選択の結果、
モデル”fit”が選ばれた。
Odds（オッズ）
＃選択されたモデルfitに推定値を入れると…
logit(pi) = -9.5294+0.2234xi
log(pi/(1-pi)) = -9.5294+0.2234xi
pi/(1-pi) = exp(-9.5294)exp(0.2234xi)
pi/(1-pi) ∝ exp(0.2234xi) 左辺は“オッズ”
#土壌含水率xiが“1単位”増加したとき、種子生残率のオッズは
pi/(1-pi) ∝ exp(0.2234(xi+1))
∝ exp(0.2234xi)exp(0.2234)
exp(0.2234)=1.2503倍,増加する。
講義のーと（久保 2008）http://eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/49477/5/kubostat2008d.pdf
まとめ
• ロジスティック回帰で、GLMの一連の流れを示した。
– プロット
– モデルの組み立て、パラメータ推定
– AICによるモデル選択
• 今回触れていないが、尤度、過大分散への対処、offset項の使い
方も大事。
• 研究計画段階で、どのようなデータが取れそうか？（仮説）、その
解析方法は？ といったところまで計画し、予習。
• 調査デザイン（調査区の選定、枠の個数・配置などなど）はとても
重要です！
ご清聴、ありがとうございます。