みっちりGLM

Size: px

Start display at page:

Download "みっちりGLM"

あゆみさかいざわ
4 years ago
Views:

1 2015/3/27 12:00-13:00 日本草地学会若手 R 統計企画 ( 信州大学農学部 ) R と一般化線形モデル入門山梨県富士山科学研究所安田泰輔謝辞 : 日本草地学会若手の会の皆様発表の機会を頂きたいへんありがとうございます!

2 茨城大学学生時代自己紹介ベータ二項分布を用いた種の空間分布の解析所属 : 山梨県富士山科学研究所最近の研究テーマ近接リモートセンシングによる半自然草地のモニタリング手法開発富士山の長期的広域的な森林限界動態湖畔における外来植物アレチウリの分布推定と駆除活動

3 趣旨説明 R に触れてみよう! GLM をやってみよう! ざっくばらんに!( ということだと思う ) 文章は若手の会作成

4 R( 主に wikipedia より ) オープンソースフリーソフト R Development Core Team(NZ) がメンテナンス統計解析向けのプログラミング言語統計解析画像処理音声合成 GIS など各種多様なパッケージが利用できる (6000 以上! だそうです ) もちろんフリー R をダウンロードしたとき全部は入っていないので必要なパッケージ ( 様々な関数群 ) をダウンロードして使う E.g. 画像処理がしたい! 調べる聞く > install.packages( raster ) > library(raster) で使えるようになる

5 グラフィック他のソフトとの連携など統計以外のパッケージも充実ファイル操作もできます日々拡張と強化が行われているので情報交換がとても大事 ( ) ( ) ( )

6 テキストベースに慣れよう! 何度もボタンを押すテキストで記述処理内容を記録できる何度も使える一括処理できるなどの利点がある Notepad++

7 今回使用するソフト :R-Studio ここにコードを記載変数情報等 R の実行図やパッケージ情報等少しやってみます

8 R で GLM 一般化線形モデル Generalized Linear Model フィールド研究ではよく使われる方法の 1 つ良書参考 HP は非常にたくさんありますデータ解析のための統計モデリング入門 ( 久保拓弥岩波書店 ) 講義のーと ( データ解析のための統計モデリング ) 名古屋大学大学院生命環境農学研究科森林生態生理学研究分野 R で GLM をやってみようなどなど

9 他の手法との関連これらはGLMの特殊なバージョン二項分布やポアソン分布ガンマ分布も使える説明変数には連続変数カテゴリ変数両方が取れる分析名応答変数 (Y) 説明変数 (X) t 検定 (Welchでない場合) 量的変数 1つ 2 分変数 1つ一元配置分散分析量的変数 1つカテゴリ変数 1つ多元配置分散分析量的変数 1つカテゴリ変数複数回帰分析量的変数 1つ量的変数 1つ重回帰分析量的変数 1つ量的変数複数 ( カテゴリ変数はダミー変数化共分散分析量的変数 1つロジスティック回帰分析 2 分変数 1つ 2 分変数カテゴリ変数量的変数複数正準相関分析量的変数複数量的変数複数中澤 (2004) R による統計解析の基礎を改編引用

10 内容ロジスティック回帰を中心に 1. プロット 2. モデルの組み立てパラメータ推定 3. AIC によるモデル選択を行ってみる * 数式は重要ですがここでは極力使いません

11 問題設定とデータセットの作成土壌水分条件に対する植物種 S の種子生残率を知りたい水分条件が異なる場所 12か所を選ぶ各地点で種 Sの種子を25~50 個採取し生残 or 死んでいるかを調査した乾燥種子サンプリング 1 生きている :1 死んでいる :0 2 湿潤 12 x: 土壌含水率 n: 各地点で調査した総種子数 y: 生残種子数

12 1. プロット説明変数に伴って変化する確率分布の期待値 ( 赤 ) を推定この確率分布はなにか? 説明変数と期待値はどのような関係か?

13 2. モデル組み立て準備生残率は生残種子数 yi/ 総種子数 ni (i = 1,2,,12) 生残種子数 yi は 0,1,,n の値を取り得る生残率 (n 個中 y 個当たり ) は二項分布に従うと仮定二項分布 Binomial(pi, ni) の期待値 pi は説明変数 xi となんらかの関係を持っていると仮定これらの条件仮定からモデルを組み立てる

14 モデル組み立て 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 データの性質に合わせて確率分布を仮定する説明変数は連続変数 ( 今回 ) カテゴリ変数その両方が使える

15 いろんなデータ型へも適用可能 GLM を行うときは応答変数の確率分布リンク関数線形予測子を指定する確率分布 ( 指数型分布族 ) とリンク関数 ( デフォルト )>?glm family データ確率分布リンク関数 ( デフォルト ) 0,1,,n( 離散上下限あり ) 二項分布 binomial logit 0,1,, ( 離散上限なし ) ポアソン分布 poisson log 0 ~ ( 連続上限なし ) ガンマ分布 Gamma Inverse (or log) - ~ ( 連続上下限なし ) 正規分布 gaussian identity ( そのまま )

16 R で GLM 実行 ( パラメータ推定 ) > fit<-glm( cbind(y,n-y)~x, #cbindで成功数失敗数を指定 family=binomial(link= logit ), data=data_set ) 注意 cbind(y,n-y)~x の右辺は切片 + 傾き *x の意味

17 結果 logit(pi) = 線形予測子 pi=1/(1+exp(-( 線形予測子 ))) Residual deviance/df = 8.897/10=0.8897

18 結果の図示 β0=-10.7, β1=0.25 β0^=-9.53(s.e. 0.84), β1^=0.22(s.e. 0.02)

19 3. まだ報告書は終わらない AIC によるモデル選択 β0=-10.7, β1=0.25 β0^=-9.53(s.e. 0.84), β1^=0.22(s.e. 0.02) 線形予測子 β0 線形予測子 β0+β1*xi β1=0( 傾き =0) と違うのか? 有意差検定をする観方を変えるどのモデルが良いのか? AIC でモデルを選択する

20 AIC 赤池情報量規準 Akaike s Information Criterion AIC = -2lnL + 2k (lnl: 最大対数尤度 k: パラメータ数 ) 同じデータのもとで 2 つ以上のモデルの AIC を比較し AIC が最小となるモデルを採用する ( 山村光司さんの HP 相対的に AIC が低いモデルが良いモデルとされる

21 再 :R で GLM 実行 ( パラメータ推定 ) > fit<-glm( cbind(y,n-y)~1, #cbindで成功数失敗数を指定 family=binomial(link= logit ), data=data_set ) 注意 cbind(y,n-y)~1 の右辺は切片だけの意味

22 切片だけのモデル

23 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 が選ばれた

24 Odds( オッズ ) # 選択されたモデル fit に推定値を入れると logit(pi) = xi log(pi/(1-pi)) = xi pi/(1-pi) = exp( )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)= 倍, 増加する講義のーと ( 久保 2008)

25 まとめロジスティック回帰で GLM の一連の流れを示したプロットモデルの組み立てパラメータ推定 AIC によるモデル選択今回触れていないが尤度過大分散への対処 offset 項の使い方も大事研究計画段階でどのようなデータが取れそうか?( 仮説 ) その解析方法は? といったところまで計画し予習調査デザイン ( 調査区の選定枠の個数配置などなど ) はとても重要です!

26 ご清聴ありがとうございます

日心TWS

日心TWS 2017.09.22 (15:40~17:10) 日本心理学会第 81 回大会 TWS ベイジアンデータ解析入門回帰分析を例にベイジアンデータ解析を体験してみる広島大学大学院教育学研究科平川真ベイジアン分析のステップ (p.24) 1) データの特定 2) モデルの定義 ( 解釈可能な ) モデルの作成 3) パラメタの事前分布の設定 4) ベイズ推論を用いてパラメタの値に確信度を再配分ベイズ推定