「Rゼミ/輪読GLM赤本」の編集履歴（バックアップ）一覧はこちら

「Rゼミ/輪読GLM赤本」(2009/01/01 (木) 11:00:19) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

*輪読Extending the linear model with R

このページでは、Julian J. Farawayの「Extending the linear model with R」の輪読会の様子を紹介します。発表ファイルやExerciseの回答、当日の議論の内容など。
----
Update! 7/4：9章と10章の発表ファイルを追加
----
累積訪問者数（from 08/05/2007）：&counter()
今日来てくれた人：&counter(today)

#contents

*ゼミの日時・場所
- 場所：N356
- 4/18（水）10時～：Introduction, Likelihood Theory（飯島）
- 4/25（水）10時～：Binomial Data（上野）
- 5/2（水）13時～：Count Data（高橋）
- 5/9（水）13時～：Contingency Tables上（渡辺）
- 5/16（水）13時～：Contingency Tables下（渡辺）
- 5/24（木）13時～：Generalized Linear Models（江口）
- 5/31（木）13時～：Other GLMs（森）
- 6/6（水）13時～：Random Effects上（志田）
- 6/20（水）13時～：Random Effects下（志田）
- 7/4（水）13時～：Repeated Measures and Longitudinal Data, Mixed Effect Models for Nonnormal Responses（飯島）

*Introduction, Appendix A
担当：飯島（造林）
発表ファイル：[[1andAppendix>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=1andAppendix.pdf]]

**当日議論したこと
***線形モデルの結果を出力するsummary()における、各変数の有意確率って？
各変数が相互に完全に独立であると仮定した場合（そんなことはほとんどの場合ありえないが）、各変数の推定された係数が、0と有意に異なっているかどうかを検定した結果。

***線形モデルにおけるカテゴリカル変数の扱いについて
Rでは、アルファベットで最も最初にあるカテゴリーの係数が0として扱われ、残りのカテゴリーの係数は、最初のカテゴリーに対する値として見る。

どのカテゴリーとどのカテゴリーが有意に異なる、ということは算出できるのかわからない。

***重みつき最小二乗法で、具体的には、どのように重み付けをしているのか？
残差を、重み（要は共変量的にばらつきを大きくさせる変数）で割ることで重み付けをしている。

**この章を読む上で役立つもの
- 「生物学を学ぶ人のための統計のはなし」粕谷英一著　文一総合出版：通称ぴんく本。Appendixを読むときに、6章の最尤法の説明がすばらしく役に立つ。


*Binomial Data
担当：上野（野生生物保護）
発表ファイル：[[Binomial>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=Binomial.pdf]]

**当日議論したこと
***Devianceとは何か？
より大きいモデルの対数尤度（Ll）とより小さいモデルの対数尤度（Ls）の差の2倍。具体的には、Llは完全にデータに当てはまっているモデル（従属変数そのもの）であり、Lsは自分の組んだモデルである。

RでGLMを行ったときに算出されるNull Devianceは、Lsが切片しか与えないモデルにおけるDevianceであり、Residual Devianceは、Lsが自分が組んだモデル（要は独立変数を与えたモデル）におけるDevianceである。

***ProspectiveとRetrospectiveな研究の違い
医学調査でよく問題となる。あるリスクに曝露されたときに、曝露された群とされなかった群で病気になる危険性がどれだけ高まるかということを問題にしたときに考える必要がある。
- Prospective（前向き）：あらかじめリスクに曝露している群と曝露していない群を同じ数だけ用意し、その後実際に病気が発祥したか調べる方法。
- Retrospective（後ろ向き）：実際にはProspectiveな方法をとれるシチュエーションがほとんどないので、病院などにいる患者に、リスクを曝露されたことがあるかどうか、そして現在病気を持っているかどうか訪ねる方法。容易に調査を行えるが、調査対象群にフィルターがかかっている（死んだ人に対して調査ができない、病気を持っている人の割合が多くなる）ため、リスクの曝露による病気の発症率の増加を簡単には評価できなくなる。

ただし、これらの方法の違いは、切片の推定値の違いにしか現れないとされている（p.35下段）。

***モデルの当てはまり
- Deviance：データに対する予測値の当てはまりを判断する基準である。定義は上述。
- PearsonのX2：これはDevianceと（ほぼ）同じ。
- R2値：この本で紹介されているR2値（p.41）は、いわゆる自由度調整済みR2値を、GLMの枠組みの中で使えるように一般化したものである。決めたモデルが、全体のばらつきの「何割」を説明しているかに興味があるときに使う。

***Overdispersionって？
誤差の分布として仮定した分布から予想されるばらつきよりもばらつきが大きいこと。二項分布やポアソン分布の分散は期待値の関数であり、個体差がまったく生じないものとして定義されている。そのため、実際のデータはしばしば予想されるばらつきよりも大きいばらつきとなり、仮定した分布に対するあてはまりが悪くなる。

***Matched case-control studiesって？
なにかしらの線形モデルを組むときに、独立変数のある要因1（カテゴリーデータ）内において、検討したい要因の影響を検討すること。要因1は、例えば、検討したい要因が病気の発生に与える影響を検討するときに、対象患者の人種、性別など。様々な人種、あるいは性別を対象に取られたデータでも、これらの要因の影響を考慮した形で、病気の発生しやすさへの検討したい要因の影響を検討できる。「対応」あるいは「ブロック」と呼ばれるものが要因1にあたる（と思う）。

Rの中では、library(survival)に入っているstrata()に要因1を入れて、独立変数としてモデルに投入することで行える。

**この章を読む上で役立つもの
- [[Rによる保健医療データ解析>http://phi.ypu.jp/msb/medstatbook.pdf]]：群馬大学の中澤先生の講義資料。12章のあたりが特に関連。
- 「医学統計学シリーズ2　統計モデル入門」丹後俊郎著　朝倉書店：6章の一般化線形モデルが参考になる。


*Count Regression
担当：高橋（生態管理）
発表ファイル：[[Count>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=Count.pdf]]

**当日議論したこと
***リンク関数と「変数変換」の違い
- リンク関数：構築した線形予測子にリンク関数をかませることで、期待値を算出する。つまり、「線形予測子」をリンク関数で変換しているだけである。リンク関数をかますのは、想定している誤差分布に線形予測子を合わせるためだけであり、結果は変数（つまり取ったデータそのもの）の影響として考えることができる。
- 変数変換：取ったデータそのものを本当に変換してしまう。つまり、それによって得られる結果は「変数変換した変数」（←当然自然にはこんなものは存在しない）の影響を評価していることになる。

***offsetがらみ
offset()の中に入れる変数にlogがついている（p.63の式）のは、poissonのデフォルトのリンク関数がlogだから。p.62の最後の式がわかりやすいが、式の右辺は線形予測子であるから、従属変数は当然logで囲まれたものとなる。従属変数のうち、offset()に含めるべき部分を独立変数に回すので、logがついている。p.62の式と、p.63の式で、Residual devianceの自由度が異なっている点に注意。

***negative binomialの使い方
> > library(MASS)
> > result <- glm.nb(..., df) #familyの指定がいらない
で実行できる。

***GLMの結果としてR2値があまり使われないのはなぜか？
（推測の域を出ないが）誤差構造に正規分布を仮定する回帰（重回帰）の場合、誤差の等分散性が仮定されているため、データのレンジによらず、与えられたデータのばらつきに対する説明できた分でR2を構築しても（それほど）問題ないが、誤差構造に正規分布以外を仮定するGLMの場合、誤差の等分散性はまったく仮定されていないため、R2値の意味が、データのレンジによって異なってしまうから。

また、サンプル数に対するばらつきの増え方が、線形的なのかもよくわからない。そのため、やはりR2の持つ意味は、データ数や範囲によって（かなり）異なると考えられるので、R2を出すことにそれほど意味がないから使われていない（のかもしれない）。Binomialの賞で説明した方法でむりくり算出はできるが。


*Contingency Tables
担当：渡辺（野生生物保護）
発表ファイル：[[Contingency>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=Contingency.pdf]]

**当日議論したこと
***drop1によるモデル比較
- ポアソンモデルでdrop1関数を用いてモデルを比較する際にカイ二乗検定を使用している（p.70）。前の章（p.60）で過分散があるときはカイ二乗検定よりもF検定を使用すべきと述べているが、過分散がない（あるいは考えなくてもいい）場合にはカイ二乗検定でいいということか？

***2×2表でのGLM
- 応答変数にmatrix()を使うと、ベクトルを使用した場合と比べてモデル式の作り方が若干変わる。（p.73-74）。

***4×4表でのGLM
- ポアソンモデルで分析している（p.76）。この際変数間の交互作用を予測変数に含めると自由度が不足するので、予測変数には交互作用項を含めないで分析を行い、算出されたデビアンスと自由度の数値から独立性を判断する。ただし判断に際しての数値的な基準は不明。

***三元表
- マンテル－ヘンツェル検定で帰無仮説（すべての部分表の条件付オッズ比が1）が棄却された場合には、部分表を見比べて変数間の関係をみることも重要。

***順序の情報の生かし方
- 例では、順序情報をもつ変数を1,2,3...という整数で置き換えて解析している。これらは順序変数ではなく間隔変数だが、これらを用いて解析することによってモデルの当てはまりがよくなれば、順序情報を使用しないモデルよりもよいと判断できるのであろう。


*Generalized Linear Models
担当：江口（造林）
発表ファイル：[[GLM>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=GLM.pdf]]

**当日議論したこと
***GLMのアルゴリズム
[[この辺>統計解析ログ/GLM]]に少しまとめてあります。

***誤差構造が違うDeviance（結局は対数尤度）を比較し、どちらが当てはまりがよいかを論じてよいのか？
- 結論から言うと可能。ただし、連続変数と離散変数の尤度は本質的に異なるものなので、比較できない。
- [[久保さんのページ>http://hosho.ees.hokudai.ac.jp/~kubo/index-j.html]]の[[この辺（11月8日）>http://hosho.ees.hokudai.ac.jp/~kubo/log/2003/1101.html#08]]に解説が。
- 例えば、binomialとpoissonの結果は比較できるが、gaussianとpoissonは比較できない。
- また、binomialでも、random effectsを含むglmmの結果とは比較できない（現在のRのRandom effectsとして用意されているのは正規分布であり、そのためglmmの尤度は離散変数と連続変数を含むから。

**この章を読む上で役立つもの
- [[「S‐PLUSによる統計解析」>http://www.amazon.co.jp/S‐PLUSによる統計解析-W-N-ヴェナブルズ/dp/4431708219/ref=pd_rhf_p_1/250-5207051-3307459]]　W.N. ヴェナブルズ, B.D. リプリー著　シュプリンガー・フェアラーク東京：一般化線形モデルの説明のところが関連。


*Other GLMs
担当：森（野生生物保護）
発表ファイル：[[OtherGLM>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=OtherGLM.pdf]]


*Random Effects
担当：志田（生態管理）
発表ファイル：[[Random effects>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=ramdom_effect.pdf]]

**当日議論したこと
***尤度比検定は保守的な結果を導きがち
尤度比検定は尤度が構築できれば行える便利な検定だが、p値が大きくなりやすい（有意になりにくい、保守的だ）傾向があるので、ブートストラップ法によって新たに検定統計量を構築することが多いようだ。

***特殊？な条件指定
あるベクトルで一定の条件を満たすものを探す場合、
> > d <- rnorm(100, 50, 20)
> > d[d > 70]
というように、[ ]を使うが、
> > (d > 70)
とすると、dの全データに関して、70より大きいかどうかが、TRUEとFALSEで返される。これを応用すると、
> > mean(d > 70)
とすることで、70以上の値が全体の何パーセントに当たるのか計算することができる（TRUEが1、FALSEが0）。

***Random effectの種類、構造
本文中では、
- Blocks
- Split
- Nested
- Crossed
- Multilevel
を挙げている。

***Random effect1: Blocks
- まずは基本。
- その変数が応答変数に与える影響には興味がないが、影響を与える要因としては考慮したい
- ex. 個体、プロットなど

***Random effect2: Split
- Random effectは1つ
- 各Random effect内において、2つの固定効果が存在。
- 各Random effect内において、片方の固定効果は何処理かが含まれるが、もう片方の固定効果は1種類しか含まれない。
- 各Random effectが固定効果の種類によって分割されている、と見るべきか。
- ex. 苗畑で、CO2付加と窒素処理による樹木の成長を見てるときに、例えば8個同じような処理区を設けていて（いわゆる反復）、一つの処理区で窒素処理はありなし両方作れるが、CO2処理は付加する、あるいはしないのどちらかしか設定できないような場合。
- 実験計画法的に考えると、完全無作為化処理ができない場合、として考えることができる。

***Random effect3: NestedとCrossed
- Random effectが2つ以上の場合
- Nested: 下位のRandom effectが、上位の各Random effect内で完全に完結し、上位のRandom effect間ではかぶらない場合（ex. 調査地が2つあり、各調査地では区画1から50まであり、各区画で植物の個体数を測定している。例えば区画1は両調査地にあるが、それは同じものではない）
- Crossed: 下位のRandom effectが、上位のRandom effect間にまたがっている場合（ex. 2つの地域があり、両地域を移動する（しない個体もいる）動物の体サイズを測定している。両地域で同じ個体が出現する場合がある）

***Random effect4: Multilevel
- NestedとCrossedのように、階層性のあるRandom effect（NestedとCrossedの総称？）
- NestedやCrossedよりも、より上位、下位という関係がはっきりしている場合なのかもしれない。
- 上位のRandom effectの影響を組み込んだ、新しい固定効果変数を作ることもできる。

**RにおけるRandom effectの指定方法
まずは必要な架空データの生成
> > Seedlings <- c(rnbinom(50, 1, 0.3), rnbinom(50, 5, 0.3))
> > Light <- Seedlings + rnorm(100, 20, 10)
> > Subp <- rep(c(rep(1, 25), rep(2, 25)), 2)
> > Plot <- rep(rep(1:25, 2), 2)
> > Site <- c(rep("S1", 50), rep("S2", 50))
> #これらをデータフレーム化
> > d <- data.frame(Seedlings, Light, Subp, Plot, Site)
> > d$Subp <- as.factor(d$Subp)
> > d$Plot <- as.factor(d$Plot)
- このデータフレームの意味: 調査地Site内に方形区Plotがあり、方形区は小方形区Subpに分けられている。Subpごとに明るさLightが測定されており、またある植物の個体数Seedlingsが測定されている。
- いわゆる説明変数がLightであり、応答変数がSeedlingsである。その他はRandom effectである。
- Random effectは、固定効果（Fixed effect）の「傾き」あるいは「切片」に対して設定できる。
- 上記のRandom effectの種類は、「傾き」あるいは「切片」どちらに対しても生じうる。
- われわれが扱うケースでは、切片にRandom effectを仮定するケースが多いか（ex. 調査地でそもそも個体数が違っており、調査地がRandom effectの場合）。
- 以下では、Random切片を中心に紹介。

***Blocks
単純に、調査地Site間でそもそもSeedlingsの数が違っていそうで、SiteをRandom effect（切片）とする場合。
> > library(lme4)
> > result <- lmer(Seedlings ~ Light + (1 | Site), family=poisson, d)
- ランダム切片の基本的な指定の仕方は、Random effectのラベルをREとすると、以下のようである。
> (1 | RE)

SiteごとにLightに対するSeedlingsの反応が異なると予想される場合は、SiteはLightに対するRandom傾きを与えるRandom effectである。
> > library(lme4)
> > result <- lmer(Seedlings ~ Light + (0 + Light | Site), family=poisson, d)
- ランダム傾きの基本的な指定の仕方は、Random effectのラベルをRE、Fixed effectのラベルをFEとすると、以下のようである。
> (0 + FE | RE)


ちなみに、Random effectが1個で、かつ誤差構造としてbinomialあるいはpoissonが仮定され、かつRandom切片指定の場合であれば、glmmML()がお手軽（モデル選択もstepAIC()で実行できる）。
> > library(glmmML)
> > result <- glmmML(Seedlings ~ Light, cluster=Site, family=poisson, d)

***Split


***Nested
Random effectとして、SiteかつPlotを想定するケース。
> > result <- lmer(Seedlings ~ Light + (1 | Site/Plot), family=poisson, method="Laplace", d)
- この場合、Siteが上位のRandom effect、Plotが下位のRandom effectである。
- 複数のRandom effectを指定する場合、数値積分による尤度の最適解を見つけることが数学的に不可能になるため、Laplace近似による方法を指定する必要がある（詳しくは[[この辺>http://hosho.ees.hokudai.ac.jp/~kubo/ce/LmerIntroduction.html#toc7]]を参照）
- ありがちで犯しやすいミス：Site間で、Plotに同じ番号をつけてしまう（Site AのPlot 1とSite BのPlot 1は違う！）


- これを防ぐ方法1: lmer()で指定するときに、以下の方法で「階層化した変数」を作る。
> > result <- lmer(Seedlings ~ Light + (1 | Site/Plot), family=poisson, method="Laplace", d)
> #これでも同じこと
> > result <- lmer(Seedlings ~ Light + (1 | Site) + (1 | Site:Plot), family=poisson, method="Laplace", d)
- これを防ぐ方法2: あらかじめ「階層化した変数」を作っておく
> > SitePlot <- c(100 + Plot[1:50], 200 + Plot[51:100])
> > d <- data.frame(Seedlings, Light, Subp, Plot, Site, SitePlot)
> #factor化処理は省略
> > result <- lmer(Seedlings ~ Light + (1 | Site) + (1 | SitePlot), family=poisson, method="Laplace", d)
>
> #こんな方法も有効（Plotの番号を調査地間でかぶらないようにする）
> > Plot2 <- c(rep(1:25, 2), rep(26:50, 2))
> > d <- data.frame(Seedlings, Light, Subp, Plot, Site, Plot2)
> #factor化処理は省略
> > result <- lmer(Seedlings ~ Light + (1 | Site) + (1 | Plot2), family=poisson, method="Laplace", d)

***Crossed


***Multilevel



*Repeated Measures and Longitudinal Data
担当：飯島（造林）
発表ファイル：[[Repeated>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=Repeated.pdf]]


*Mixed Effect Models for Nonnormal Responses
担当：飯島（造林）
発表ファイル：[[GLMM>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=GLMM.pdf]]

*輪読Extending the linear model with R

このページでは、Julian J. Farawayの「Extending the linear model with R」の輪読会の様子を紹介します。発表ファイルやExerciseの回答、当日の議論の内容など。
----
Update! 7/4：9章と10章の発表ファイルを追加
----
累積訪問者数（from 08/05/2007）：&counter()
今日来てくれた人：&counter(today)

#contents

*ゼミの日時・場所
- 場所：N356
- 4/18（水）10時～：Introduction, Likelihood Theory（飯島）
- 4/25（水）10時～：Binomial Data（上野）
- 5/2（水）13時～：Count Data（高橋）
- 5/9（水）13時～：Contingency Tables上（渡辺）
- 5/16（水）13時～：Contingency Tables下（渡辺）
- 5/24（木）13時～：Generalized Linear Models（江口）
- 5/31（木）13時～：Other GLMs（森）
- 6/6（水）13時～：Random Effects上（志田）
- 6/20（水）13時～：Random Effects下（志田）
- 7/4（水）13時～：Repeated Measures and Longitudinal Data, Mixed Effect Models for Nonnormal Responses（飯島）

*Introduction, Appendix A
担当：飯島（造林）
発表ファイル：[[1andAppendix>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=1andAppendix.pdf]]

**当日議論したこと
***線形モデルの結果を出力するsummary()における、各変数の有意確率って？
各変数が相互に完全に独立であると仮定した場合（そんなことはほとんどの場合ありえないが）、各変数の推定された係数が、0と有意に異なっているかどうかを検定した結果。

***線形モデルにおけるカテゴリカル変数の扱いについて
Rでは、アルファベットで最も最初にあるカテゴリーの係数が0として扱われ、残りのカテゴリーの係数は、最初のカテゴリーに対する値として見る。

どのカテゴリーとどのカテゴリーが有意に異なる、ということは算出できるのかわからない。

***重みつき最小二乗法で、具体的には、どのように重み付けをしているのか？
残差を、重み（要は共変量的にばらつきを大きくさせる変数）で割ることで重み付けをしている。

**この章を読む上で役立つもの
- 「生物学を学ぶ人のための統計のはなし」粕谷英一著　文一総合出版：通称ぴんく本。Appendixを読むときに、6章の最尤法の説明がすばらしく役に立つ。


*Binomial Data
担当：上野（野生生物保護）
発表ファイル：[[Binomial>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=Binomial.pdf]]

**当日議論したこと
***Devianceとは何か？
より大きいモデルの対数尤度（Ll）とより小さいモデルの対数尤度（Ls）の差の2倍。具体的には、Llは完全にデータに当てはまっているモデル（従属変数そのもの）であり、Lsは自分の組んだモデルである。

RでGLMを行ったときに算出されるNull Devianceは、Lsが切片しか与えないモデルにおけるDevianceであり、Residual Devianceは、Lsが自分が組んだモデル（要は独立変数を与えたモデル）におけるDevianceである。

***ProspectiveとRetrospectiveな研究の違い
医学調査でよく問題となる。あるリスクに曝露されたときに、曝露された群とされなかった群で病気になる危険性がどれだけ高まるかということを問題にしたときに考える必要がある。
- Prospective（前向き）：あらかじめリスクに曝露している群と曝露していない群を同じ数だけ用意し、その後実際に病気が発祥したか調べる方法。
- Retrospective（後ろ向き）：実際にはProspectiveな方法をとれるシチュエーションがほとんどないので、病院などにいる患者に、リスクを曝露されたことがあるかどうか、そして現在病気を持っているかどうか訪ねる方法。容易に調査を行えるが、調査対象群にフィルターがかかっている（死んだ人に対して調査ができない、病気を持っている人の割合が多くなる）ため、リスクの曝露による病気の発症率の増加を簡単には評価できなくなる。

ただし、これらの方法の違いは、切片の推定値の違いにしか現れないとされている（p.35下段）。

***モデルの当てはまり
- Deviance：データに対する予測値の当てはまりを判断する基準である。定義は上述。
- PearsonのX2：これはDevianceと（ほぼ）同じ。
- R2値：この本で紹介されているR2値（p.41）は、いわゆる自由度調整済みR2値を、GLMの枠組みの中で使えるように一般化したものである。決めたモデルが、全体のばらつきの「何割」を説明しているかに興味があるときに使う。

***Overdispersionって？
誤差の分布として仮定した分布から予想されるばらつきよりもばらつきが大きいこと。二項分布やポアソン分布の分散は期待値の関数であり、個体差がまったく生じないものとして定義されている。そのため、実際のデータはしばしば予想されるばらつきよりも大きいばらつきとなり、仮定した分布に対するあてはまりが悪くなる。

***Matched case-control studiesって？
なにかしらの線形モデルを組むときに、独立変数のある要因1（カテゴリーデータ）内において、検討したい要因の影響を検討すること。要因1は、例えば、検討したい要因が病気の発生に与える影響を検討するときに、対象患者の人種、性別など。様々な人種、あるいは性別を対象に取られたデータでも、これらの要因の影響を考慮した形で、病気の発生しやすさへの検討したい要因の影響を検討できる。「対応」あるいは「ブロック」と呼ばれるものが要因1にあたる（と思う）。

Rの中では、library(survival)に入っているstrata()に要因1を入れて、独立変数としてモデルに投入することで行える。

**この章を読む上で役立つもの
- [[Rによる保健医療データ解析>http://phi.ypu.jp/msb/medstatbook.pdf]]：群馬大学の中澤先生の講義資料。12章のあたりが特に関連。
- 「医学統計学シリーズ2　統計モデル入門」丹後俊郎著　朝倉書店：6章の一般化線形モデルが参考になる。


*Count Regression
担当：高橋（生態管理）
発表ファイル：[[Count>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=Count.pdf]]

**当日議論したこと
***リンク関数と「変数変換」の違い
- リンク関数：構築した線形予測子にリンク関数をかませることで、期待値を算出する。つまり、「線形予測子」をリンク関数で変換しているだけである。リンク関数をかますのは、想定している誤差分布に線形予測子を合わせるためだけであり、結果は変数（つまり取ったデータそのもの）の影響として考えることができる。
- 変数変換：取ったデータそのものを本当に変換してしまう。つまり、それによって得られる結果は「変数変換した変数」（←当然自然にはこんなものは存在しない）の影響を評価していることになる。

***offsetがらみ
offset()の中に入れる変数にlogがついている（p.63の式）のは、poissonのデフォルトのリンク関数がlogだから。p.62の最後の式がわかりやすいが、式の右辺は線形予測子であるから、従属変数は当然logで囲まれたものとなる。従属変数のうち、offset()に含めるべき部分を独立変数に回すので、logがついている。p.62の式と、p.63の式で、Residual devianceの自由度が異なっている点に注意。

***negative binomialの使い方
> > library(MASS)
> > result <- glm.nb(..., df) #familyの指定がいらない
で実行できる。

***GLMの結果としてR2値があまり使われないのはなぜか？
（推測の域を出ないが）誤差構造に正規分布を仮定する回帰（重回帰）の場合、誤差の等分散性が仮定されているため、データのレンジによらず、与えられたデータのばらつきに対する説明できた分でR2を構築しても（それほど）問題ないが、誤差構造に正規分布以外を仮定するGLMの場合、誤差の等分散性はまったく仮定されていないため、R2値の意味が、データのレンジによって異なってしまうから。

また、サンプル数に対するばらつきの増え方が、線形的なのかもよくわからない。そのため、やはりR2の持つ意味は、データ数や範囲によって（かなり）異なると考えられるので、R2を出すことにそれほど意味がないから使われていない（のかもしれない）。Binomialの賞で説明した方法でむりくり算出はできるが。


*Contingency Tables
担当：渡辺（野生生物保護）
発表ファイル：[[Contingency>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=Contingency.pdf]]

**当日議論したこと
***drop1によるモデル比較
- ポアソンモデルでdrop1関数を用いてモデルを比較する際にカイ二乗検定を使用している（p.70）。前の章（p.60）で過分散があるときはカイ二乗検定よりもF検定を使用すべきと述べているが、過分散がない（あるいは考えなくてもいい）場合にはカイ二乗検定でいいということか？

***2×2表でのGLM
- 応答変数にmatrix()を使うと、ベクトルを使用した場合と比べてモデル式の作り方が若干変わる。（p.73-74）。

***4×4表でのGLM
- ポアソンモデルで分析している（p.76）。この際変数間の交互作用を予測変数に含めると自由度が不足するので、予測変数には交互作用項を含めないで分析を行い、算出されたデビアンスと自由度の数値から独立性を判断する。ただし判断に際しての数値的な基準は不明。

***三元表
- マンテル－ヘンツェル検定で帰無仮説（すべての部分表の条件付オッズ比が1）が棄却された場合には、部分表を見比べて変数間の関係をみることも重要。

***順序の情報の生かし方
- 例では、順序情報をもつ変数を1,2,3...という整数で置き換えて解析している。これらは順序変数ではなく間隔変数だが、これらを用いて解析することによってモデルの当てはまりがよくなれば、順序情報を使用しないモデルよりもよいと判断できるのであろう。


*Generalized Linear Models
担当：江口（造林）
発表ファイル：[[GLM>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=GLM.pdf]]

**当日議論したこと
***GLMのアルゴリズム
[[この辺>統計解析ログ/GLM]]に少しまとめてあります。

***誤差構造が違うDeviance（結局は対数尤度）を比較し、どちらが当てはまりがよいかを論じてよいのか？
- 結論から言うと可能。ただし、連続変数と離散変数の尤度は本質的に異なるものなので、比較できない。
- [[久保さんのページ>http://hosho.ees.hokudai.ac.jp/~kubo/index-j.html]]の[[この辺（11月8日）>http://hosho.ees.hokudai.ac.jp/~kubo/log/2003/1101.html#08]]に解説が。
- 例えば、binomialとpoissonの結果は比較できるが、gaussianとpoissonは比較できない。
- また、binomialでも、random effectsを含むglmmの結果とは比較できない（現在のRのRandom effectsとして用意されているのは正規分布であり、そのためglmmの尤度は離散変数と連続変数を含むから。

**この章を読む上で役立つもの
- [[「S‐PLUSによる統計解析」>http://www.amazon.co.jp/S‐PLUSによる統計解析-W-N-ヴェナブルズ/dp/4431708219/ref=pd_rhf_p_1/250-5207051-3307459]]　W.N. ヴェナブルズ, B.D. リプリー著　シュプリンガー・フェアラーク東京：一般化線形モデルの説明のところが関連。


*Other GLMs
担当：森（野生生物保護）
発表ファイル：[[OtherGLM>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=OtherGLM.pdf]]


*Random Effects
担当：志田（生態管理）
発表ファイル：[[Random effects>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=ramdom_effect.pdf]]

**当日議論したこと
***尤度比検定は保守的な結果を導きがち
尤度比検定は尤度が構築できれば行える便利な検定だが、p値が大きくなりやすい（有意になりにくい、保守的だ）傾向があるので、ブートストラップ法によって新たに検定統計量を構築することが多いようだ。

***特殊？な条件指定
あるベクトルで一定の条件を満たすものを探す場合、
> > d <- rnorm(100, 50, 20)
> > d[d > 70]
というように、[ ]を使うが、
> > (d > 70)
とすると、dの全データに関して、70より大きいかどうかが、TRUEとFALSEで返される。これを応用すると、
> > mean(d > 70)
とすることで、70以上の値が全体の何パーセントに当たるのか計算することができる（TRUEが1、FALSEが0）。

***Random effectの種類、構造
本文中では、
- Blocks
- Split
- Nested
- Crossed
- Multilevel
を挙げている。

***Random effect1: Blocks
- まずは基本。
- その変数が応答変数に与える影響には興味がないが、影響を与える要因としては考慮したい
- ex. 個体、プロットなど

***Random effect2: Split
- Random effectは1つ
- 各Random effect内において、2つの固定効果が存在。
- 各Random effect内において、片方の固定効果は何処理かが含まれるが、もう片方の固定効果は1種類しか含まれない。
- 各Random effectが固定効果の種類によって分割されている、と見るべきか。
- ex. 苗畑で、CO2付加と窒素処理による樹木の成長を見てるときに、例えば8個同じような処理区を設けていて（いわゆる反復）、一つの処理区で窒素処理はありなし両方作れるが、CO2処理は付加する、あるいはしないのどちらかしか設定できないような場合。
- 実験計画法的に考えると、完全無作為化処理ができない場合、として考えることができる。

***Random effect3: NestedとCrossed
- Random effectが2つ以上の場合
- Nested: 下位のRandom effectが、上位の各Random effect内で完全に完結し、上位のRandom effect間ではかぶらない場合（ex. 調査地が2つあり、各調査地では区画1から50まであり、各区画で植物の個体数を測定している。例えば区画1は両調査地にあるが、それは同じものではない）
- Crossed: 下位のRandom effectが、上位のRandom effect間にまたがっている場合（ex. 2つの地域があり、両地域を移動する（しない個体もいる）動物の体サイズを測定している。両地域で同じ個体が出現する場合がある）

***Random effect4: Multilevel
- NestedとCrossedのように、階層性のあるRandom effect（NestedとCrossedの総称？）
- NestedやCrossedよりも、より上位、下位という関係がはっきりしている場合なのかもしれない。
- 上位のRandom effectの影響を組み込んだ、新しい固定効果変数を作ることもできる。

**RにおけるRandom effectの指定方法
まずは必要な架空データの生成
> > Seedlings <- c(rnbinom(50, 1, 0.3),
> + rnbinom(50, 5, 0.3))
> > Light <- Seedlings + rnorm(100, 20, 10)
> > Subp <- rep(c(rep(1, 25), rep(2, 25)), 2)
> > Plot <- rep(rep(1:25, 2), 2)
> > Site <- c(rep("S1", 50), rep("S2", 50))
> #これらをデータフレーム化
> > d <- data.frame(Seedlings, Light,
> + Subp, Plot, Site)
> > d$Subp <- as.factor(d$Subp)
> > d$Plot <- as.factor(d$Plot)
- このデータフレームの意味: 調査地Site内に方形区Plotがあり、方形区は小方形区Subpに分けられている。Subpごとに明るさLightが測定されており、またある植物の個体数Seedlingsが測定されている。
- いわゆる説明変数がLightであり、応答変数がSeedlingsである。その他はRandom effectである。
- Random effectは、固定効果（Fixed effect）の「傾き」あるいは「切片」に対して設定できる。
- 上記のRandom effectの種類は、「傾き」あるいは「切片」どちらに対しても生じうる。
- われわれが扱うケースでは、切片にRandom effectを仮定するケースが多いか（ex. 調査地でそもそも個体数が違っており、調査地がRandom effectの場合）。
- 以下では、Random切片を中心に紹介。

***Blocks
単純に、調査地Site間でそもそもSeedlingsの数が違っていそうで、SiteをRandom effect（切片）とする場合。
> > library(lme4)
> > result <- lmer(Seedlings ~ Light +
> + (1 | Site),
> + family=poisson, d)
- ランダム切片の基本的な指定の仕方は、Random effectのラベルをREとすると、以下のようである。
> (1 | RE)

SiteごとにLightに対するSeedlingsの反応が異なると予想される場合は、SiteはLightに対するRandom傾きを与えるRandom effectである。
> > library(lme4)
> > result <- lmer(Seedlings ~ Light +
> + (0 + Light | Site),
> + family=poisson, d)
- ランダム傾きの基本的な指定の仕方は、Random effectのラベルをRE、Fixed effectのラベルをFEとすると、以下のようである。
> (0 + FE | RE)


ちなみに、Random effectが1個で、かつ誤差構造としてbinomialあるいはpoissonが仮定され、かつRandom切片指定の場合であれば、glmmML()がお手軽（モデル選択もstepAIC()で実行できる）。
> > library(glmmML)
> > result <- glmmML(Seedlings ~ Light,
> + cluster=Site, family=poisson, d)

***Split


***Nested
Random effectとして、SiteかつPlotを想定するケース。
> > result <- lmer(Seedlings ~ Light +
> + (1 | Site/Plot),
> + family=poisson, method="Laplace", d)
- この場合、Siteが上位のRandom effect、Plotが下位のRandom effectである。
- 複数のRandom effectを指定する場合、数値積分による尤度の最適解を見つけることが数学的に不可能になるため、Laplace近似による方法を指定する必要がある（詳しくは[[この辺>http://hosho.ees.hokudai.ac.jp/~kubo/ce/LmerIntroduction.html#toc7]]を参照）
- ありがちで犯しやすいミス：Site間で、Plotに同じ番号をつけてしまう（Site AのPlot 1とSite BのPlot 1は違う！）


- これを防ぐ方法1: lmer()で指定するときに、以下の方法で「階層化した変数」を作る。
> > result <- lmer(Seedlings ~ Light +
> + (1 | Site/Plot),
> + family=poisson, method="Laplace", d)
> #これでも同じこと
> > result <- lmer(Seedlings ~ Light +
> + (1 | Site) + (1 | Site:Plot),
> + family=poisson, method="Laplace", d)
- これを防ぐ方法2: あらかじめ「階層化した変数」を作っておく
> > SitePlot <- c(100 + Plot[1:50], 200 + Plot[51:100])
> > d <- data.frame(Seedlings, Light, Subp,
> + Plot, Site, SitePlot)
> #factor化処理は省略
> > result <- lmer(Seedlings ~ Light +
> + (1 | Site) + (1 | SitePlot),
> + family=poisson, method="Laplace", d)
>
> #こんな方法も有効（Plotの番号を調査地間でかぶらないようにする）
> > Plot2 <- c(rep(1:25, 2), rep(26:50, 2))
> > d <- data.frame(Seedlings, Light, Subp,
> + Plot, Site, Plot2)
> #factor化処理は省略
> > result <- lmer(Seedlings ~ Light +
> + (1 | Site) + (1 | Plot2),
> + family=poisson, method="Laplace", d)

***Crossed


***Multilevel



*Repeated Measures and Longitudinal Data
担当：飯島（造林）
発表ファイル：[[Repeated>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=Repeated.pdf]]


*Mixed Effect Models for Nonnormal Responses
担当：飯島（造林）
発表ファイル：[[GLMM>http://www7.atwiki.jp/hayatoiijima?cmd=upload&act=open&pageid=33&file=GLMM.pdf]]