応用動物行動学における統計解析の進展(後編)

誌名

誌名 Animal behaviour and management ISSN

ISSN 18802133 著者

著者 多田, 慎吾

新村, 毅 巻/号

巻/号 53巻3号

掲載ページ

掲載ページ p. 117-126 発行年月

発行年月 2017年9月

農林水産省 農林水産技術会議事務局筑波産学連携支援センター

Tsukuba Business-Academia Cooperation Support Center, Agriculture, Forestry and Fisheries Research Council Secretariat

ー 解 説 一

応用動物行動学における統計解析の進展後編：一般化線形モデル

多 田 慎 吾1* . 新 村 毅

1農研機構北海道農業研究センター酪農研究領城，札幌市， 062‑8555

2東京農工大学大学院農学研究院，府中， 183‑8538

℃ orresponding author. E‑mail address: [email protected] (多田慎吾），

[email protected] (新村毅）

要 約

近年、一般線形モデルを拡張して、正規分布以外の分布も仮定できるようにした一般化線形モデルが 普及してきており、応用動物行動学および家畜管理学分野でも一層の利用が期待される。前編では、フ リーのソフトウェアである R を用いて、一般線形モデルの概念と実データを用いた解析例について解 説した。本稿では、後編として、一般線形モデルの拡張である一般化線形モデルについて、その概念と 練習用データを用いた解析例を、実際のRのコマンドを示しつつ解説する。

キーワード：分布、一般化線形モデル、一般線形モデル、行動、 R

Animal Behaviour and Management, 53 (3):  117‑126, 2017  (2016.  2.  24受付； 2017.  5.  17受理）

はじめに

本解説記事の前編では、 Rを活用しながら、正 規分布を前提とする一般線形モデルの概要を述 べた。本稿ではまず、データによってはなぜこの 手法をそのまま用いることができないのか説明す る。その後、いよいよ一般線形モデルを拡張し た一般化線形モデルはとは何か、具体例を示し ながら解説する。なお、練習用のデータおよび Rのコマンドをまとめた付録Excelファイルは、

本解説が掲載されている J‑STAGEのWebペー ジ (URL:https://www.jstage.jst.go.jp/browse/‑ char/ja/)もしくは応用動物行動学会のWebペー

ジ(URL:http:/ /www.jsaab.org/)からダウンロー ド可能である。

一般線形モデルが適用できない場合とは？

一般線形モデルの前提は「データの誤差が正規 分布に従う」であった。しかしながら、行動学分 野では誤差が正規分布に従わないデータも多く扱 う。代表的なものは、ある時間の中で観察された 行動の回数（例えば敵対行動の回数や排糞回数な ど）といったカウントデータである。カウントデー タの特徴として、そもそも 0以上の整数値しか とらないという点が挙げられる。正規分布は小数 も含めた連続量を表す分布であり、この点だけで

も正規分布はカウントデータに適合しないといえ る。また、分布の形状も正規分布のように左右対 称にならない場合が多く、この点でも正規分布を 前提としたモデルに当てはめるのはふさわしくな しヽ。

もう一つ、正規分布に従わない代表的なデー タとして0 1の値で示される割合データがあ る。例えば、観察群のうち発情行動を示した頭数 の割合などである。カウントデータとは異なり小 数値をとる連続値という点では正規分布と共通 するが、何しろこの種のデータは最小値が0、最 大値が1と定まっている。もちろん正規分布は0

より小さな値も 1より大きな値もとりうるので、

このようなデータに一般線形モデルを当てはめる と、割合データなのにモデル予測が負の値をとっ たり、 1より大きな値をとったりというような明

らかな矛盾が起きる。

これらの誤差が正規分布しないデータに対処す る方法として、これまで変数変換やノンパラメト リック法が用いられてきた。変数変換はデータの 対数をとる、正弦変換するなどの操作を行なうこ とでデータの分布を正規分布に近づける方法であ る。しかしながら、特に誤差のばらつきの大きい データで、変数変換を行なっても正規性を仮定で きない場合も多い。一方、ノンパラメトリック検 定はデータに誤差分布の仮定をしない統計手法で

0.4 

0.3 

2 n i  

梱睾

0.1 

゜ ゜

2 

二

4  6  8  観察された回数

10  12 

図13.ポアソン分布の形状

入がデータの平均値に一致する。入が小さいときは極度に右裾の広い形状を とるが、入が大きくなるにつれて左右対称の形状に近づく。

あり、どんなデータにでも適用できるのが特長で ある。しかし、データを一度順位情報に変換して 計算が行われるなど、データの情報の多くを削っ てしまうため、分布を仮定した統計方法に比べて 検出力が低い。また、複数の要因や交互作用の検 討が困難であるといった問題もある。

ここで、 一般線形モデルのデータを足し算の関 係で表すという部分はそのままにして、正規分布 以外の分布も扱えるように拡張したものが一般 化 線形モデルである。上で述べたカウントデータに はポアソン分布が（図13)、 割 合 デ ー タ に は 二 項 分 布 （ 図 14)と い っ た 分 布 に 従 う こ と が 多 い と されている。一般化線形モデルでは、これらの正 規分布以外の分布を扱うためにリンク関数といっ た概念が導入されるが、考え方は一般線形モデル

0.4 

0.3 

2 

゜

憾睾

0.1 

と変わらない。以下では、カウントデータ、割合 データそれぞれの一般化線形モデルについて実例 を示しながら説明する。

ここで、元のデータの集団がポアソン分布や二 項 分 布 、 ま た そ の 他 の 正 規 分 布 で は な い 分 布 で あったとしても、その平均値の誤差分布は中心極 限定理により正規分布に近似できる場合があると いうことを注意しておきたい。この例としてよく 挙げられるのが、 6面 サ イ コ ロ の 目 の で る 頻 度 で

ある。 普通のサイコロであれば、 1 6の目が同 程度の頻度で出るはずである。この頻度分布は、

明らかに正規分布ではない。しかし、サイコロを 10回 ふ っ て 、 そ の10回 の 平 均 値 を 算 出 す る と いった操作を何度も行なうとどうだろう？きっと 平 均 値 の 頻 度 分 布 は 、 平 均3.5の正規分布に従う

゜ ゜

²  観 察 さ れ た 回 数^{4  6} 

8  10 

図14.二項分布の形状

確率pで起こる事象を N回行なったとき、何回その事象が起きるかを表す分 布である（図は N= 10の場合）。 pX Nが期待値と一致する。 pが小さいと

き右裾の広い形状をとり、 pが0.5に近いと左右対称の形状に近づく。

多田・新村

‑< 

‑2  ‑1 

゜

2 

μ+ax  図 15.A = e<μ+ a X)とした対数リンク関数

このように入は常に正の値をとる。 eは自然対数の底を表し、値は2.718・・・である。

ものとなる。すなわち、元のデータの誤差分布が 正規分布ではなかったとしても、平均値の誤差分 布を統計解析に用いるならば、正規分布を前提と した一般線形モデルを用いるのが妥当な場合が多 い。いずれにせよ、検討するデータの誤差分布を 自分で確認することが重要である。

カウントデータの解析

ここでは、動物1頭あたりの飼養面積が異なる 場合の敵対行動の回数を表した架空データを用い て説明する（付録Excel、シート 5;通常の飼養 面積を対照群 (Tl)とし、飼養面積を狭くした 小面積群を処理群 (T2)とする；帰無仮説は「試 験処理に効果はない」である）。前述のように、

敵対行動の回数のようなカウントデータは、 0以 上の整数値しかとらないという特徴がある。この 条件を満たす分布の 1つがポアソン分布である。

正規分布は、平均値μと分散6の2つのパラメー タによって定義されるのに対し、ポアソン分布は 1つのパラメータ入（ラムダ）によって定義され る分布である。入はデータの平均値に一致し、非 負の値しかとらないパラメータで、図13に示し たようにこの入の値によって分布の形状は異な る。

一般線形モデルにおいてはデータに最も当ては まるように、

y=μ+ux+E 

の式のパラメータを推定したが、ポアソン分布を 用いる一般化線形モデルにおいてはデータに最も 当てはまるような入の値をもとめることが目的で ある。しかし、ここで一般線形モデルのときと同 様に単純に、

入=μ+llX 

のような線形式としてしまうと問題が生じる。こ

のような式の場合、μや aの値によっては、入が 負の値をとりうることになる。これに対処する方 法が、リンク関数の導入である。すなわち、線形 式を以下のようにする。

入=^e(μ+ax) 

このようにすればμ+axがどんな値であっても、

入は 0より大きな値しかとらない（図 15)。この リンク関数は、一般に対数リンク関数と呼ばれる。

以上のことを、 R を用いて解析する場合、以下 のコマンドで一般化線形モデルのパラメータ推定 をすることができる。

1.  データを "dat5"として読み込む。

datS 

< ‑

read.table("clipboard", header= T)  2.  "dat5"のデータを用いて、説明変数として

Tを含め、対数リンク関数を用いて誤差分 布にポアソン分布を指定した一般化線形モ デルを "fm5a"とする。ここでは本解説記事 前編で用いたlm関数ではなく、 glm関数を 用いている。 lm関数と同様に、カッコ内に は、初めにモデルの構造を「y T」のように 記載し、 Tを説明変数と指定している。次に 使用するデータ（今回の場合dat5) を記載 しているのも同様だが、 glm関数ではその後 ろでさらに用いる分布とリンク関数を指定 している（今回の場合family=poisson(link= 

"log"))。「family= poisson」で分布がポアソ ン分布 (Poissondistribution)であることを、

「link= "log"」でリンク関数が対数リンク関 数 (Logarithm(log)  link function) を指定し たことを示している。

fmSa 

< ‑

glm(y T,  data=datS,  family=poisson(link = "log")) 

3.  "fm5a"の要約を表示する（図16)。

> datS <‑read.table("clipboard",  header= T) 

> fmSa <‑glm(y T, data=datS, family=poisson(link = "log")) 

> summary(fmSa)  Call: 

glm(formula = y   T,  family= poisson(link = "log"),  data= datS)  Deviance Residuals: 

Min  lQ  Median  3Q  Max 

‑2.2978  ‑0.7834  0.0000  0.6579  1.7785  Coefficients: 

Estimate Std.  Error z value Pr(>lzl)  (Intercept)  0.6932  0.1000  6.931 4.16e‑12 *** 

TT2  0.2776  0.1326  2.094  0.0362 * 

Signif.  codes:  0 '***'0.001 '**'0.01 '*'0.05 '.' 0.1 ' ' 1   (Dispersion parameter for poisson family taken to be 1) 

Null deviance:  99.270  on 99  degrees of freedom  Residual deviance:  94.842  on 98  degrees of freedom  AIC: 348.01 

Number of Fisher Scoring  iterations:  5 

図16.Rの出力結果（付録Excel、シート 5その1)

summary(fm5a) 

図16の(Intercept)がμを、 TT2がa.を示してい る。各群の入は対数リンク関数にこれらの値を代 入することで求められる。すなわち、対照群では 入=e<0.6932) = 2̲0、小面積群では入=e co.6932 + 0.2116) =  2.6となる。対照群と小面積群それぞれについて、

推定したパラメータから得られたポアソン分布と データをあわせたものを図17に示した。それぞ れポアソン分布がデータに当てはまっていること が分かる。

有意差検定についても、基本的には一般線形モ デルと同様の考え方で行なう。すなわち、検討し たい要因を含んだモデルと含まないモデルでデー

対照群 0.4 

0.3 

2 

゜

憾啜

0.1 

タに対する当てはまり具合に差があるかを検討す る。ただし、データの当てはまりの指標としては、

正規分布を前提としたモデルで用いる F値では なく Devianceで示されるモデル間の対数尤度の 差を用いる。今回の例の場合、 Devianceの頻度 分布はパラメトリック・ブートストラップ法でも とめられる。また、 Devianceは

x

二乗分布に従 うことから、 x二乗分布を用いて P値を計算す ることもできる。Rでのコマンドは以下の通りで ある。

4. "dat5"のデータを用い、説明変数を含めず、

対数リンク関数を用いて誤差分布にポアソン 分布を指定した一般化線形モデルを "fm5b"

処理群

J . .  = 2.0 

},.̲= 2.6 

゜

0  1  2  3  4 5 6  7  8  9  10  0  1 2  3 4 5  6  7  8  9 10 

観察された回数 観察された回数

図17.データヘのポアソン分布の当てはめ

バーがデータの頻度、実線がポアソン分布を示す。対照群では入=eo.s932 = 2̲0、処理群では 入 =e(0.6932+0.2TT6) = 2.6と推定された。

多田・新村

> fm5b <‑glm(yl,  data=dat5,  family=poisson) 

> anova(fm5a,  fm5b,  test="Chi")  Analysis of Deviance Table  Model  1:  y   T 

Model 2:  y   1 

Resid.  Df Resid.  Dev Df Deviance Pr(>Chi)  1  98  94.842 

2  99  99.270  ‑1  ‑4.4279  0.03536 * 

Signif.  codes:  0'***'0.001'**'0.01'*'0.05'.'0.1''1 

> n <‑10000 

> DValue <‑numeric(n) 

> for  (i  in 1:n)  { 

+ dat5$y2 <‑simulate(fm5b) [,1] 

+ fm5z <‑glm(y2 T,data=dat5,family=poisson) 

+ DValue[i] <‑summary(fm5z)$null.deviance[l]‑summary(fm5z)$deviance[l] 

＋｝ 

> sum(DValue>=(summary(fm5a)$null.deviance[l]‑summary(fm5a)$deviance[l]))/n  [1]  0. 0323 

図18.Rの出力結果（付録Excel、シート 5その2)

パラメトリック・ブートストラップ法による P値の算出結果（この図では0.0323)は生成された乱 数によるのでこの図の値と一致するとは限らない。

とする。

fm5b 

< ‑

glm(y‑1, data=dat5, family=poisson)  5.  X二乗分布を用いたP値を算出する (fm5a

とfm5bのデータヘの当てはまりを x二乗検 定で比較）。

anova(fm5a, fm5b, test=℃ hi") 

6. 以下ではパラメトリック・ブートストラップ により P値を算出している。

n 

< ‑

10000 

DValue 

< ‑

numeric(n)  for (i  in  1 :n) { 

dat5$y2 

< ‑

simulate(fm5b)[, 1]  fm5z 

< ‑

glm(y2‑T,data=dat5,famil  y=poisson) 

DValue[i] 

< ‑

^{summary(fm5z)$null.} 

d e v i a n c e [ 1  ]‑

summary(fm5z)$deviance[1] 

｝ 

7.  生成データから得たx二乗値の分布におい て、実際のデータから観測された

x

二乗値よ り大きかったものの数をカウントし、繰り返 し数で除して頻度を算出する（これがP値）。 sum(DValue >= (summary(fm5a)$null. 

deviance[1 ]‑summary(fm5a)$deviance[1 ]))/  n 

x二 乗 検 定 に よ る P値 (0.03536) とパラメト リック・ブートストラップ法から算出された値 (0.0323) が ほ ぼ 等 し い こ と が 分 か る （ 図18)。 これらの結果から、今回実際のデータから得られ たDevianceの値は「試験処理に効果はない」と するならば、 3 4%程度の低頻度でしか観察さ れないものということになる。 5% (0.05) を基

準とすると、帰無仮説が棄却され、対立仮説「試 験処理に効果がある」が採択される。すなわち、

飼養面積を狭くすると敵対行動の回数が増加する という仮説を支持することができる。

（参考） glm関数で一般線形モデルの解析

一般線形モデルの解析（当解説前編）ではIm 関 数 を 用 い た が （ 付 録Excel、シート 1 4

 

)、 glm関数を用いても同様の解析を行なうことがで

きる。具体的には、 familyのところに正規分布を 表すgaussian、linkの部分にリンク関数を用いな い（そのままである）ことを表すidentityを指定

してやればよい（付録Excel、シート 6)。 1.  データを"dat6"として読み込む。

dat6 

< ‑

read.table("clipboard", header= T)  2.  "dat6"のデータを用いて、説明変数としてx

を含め、誤差分布に正規分布を指定した一般 化線形モデルの要約を表示する。

summary(glm(y‑x, dat6, family=gaussian(  link="identity"))) 

3. Im関数を用いた場合の一般線形モデルの要 約を表示する。

summary(fmくーlm(yx,dat6)) 

アウトプットがglm関数と Im関数で一致してい ることが分かる（図 19)。

割合データの解析

割合データの場合にも、基本的にはカウント データの一般化線形モデルと同様にして、一般線 形モデルの拡張を行なう。データは最小値が0、 最大値が 1であるデータを表せる分布が二項分

> dat6 <‑ read.table("clipboard",  header= T) 

> summary(glm(y x,  dat6,  family=gaussian(link="identity")))  Call: 

glm(formula = y   x,  family= gaussian(link = "identity"),  data= dat6)  Deviance Residuals: 

Min  lQ  Median  3Q  Max 

‑5.3595  ‑2.2804  0.3482  1.9505  6.2005  Coefficients: 

Estimate Std.  Error t value Pr(>ltl)  (Intercept)  9. 5495  0. 6451  14. 802  <2e‑16  *** 

xtreatment  2.0645  0.9123  2.263  0.0294  * 

Signif.  codes:  0'***'0.001'**'0.01'*'0.05'.'0.1''1  (Dispersion parameter for gaussian family taken to be 8.323794) 

Null deviance:  358.93  on 39  degrees  of freedom  Residual deviance:  316.30  on 38  degrees  of freedom  AIC:  202.23 

Number of Fisher Scoring iterations:  2 

> summary(lm(y x,  dat6))  Call: 

lm(formula = y   x,  data= dat6)  Residuals: 

Min  lQ  Median  3Q  Max 

‑5.3595  ‑2.2804  0.3483  1.9505  6.2005  Coefficients: 

Estimate Std.  Error t value Pr(>ltl)  (Intercept)  9. 5495  0. 6451  14. 802  <2e‑16  *** 

xtreatment  2.0645  0.9123  2.263  0.0294  * 

Signif.  codes:  0'***'0.001'**'0.01'*'0.05'.'0.1''1  Residual  standard error:  2.885  on 38  degrees of freedom 

Multiple R‑squared:  0.1187,  Adjusted R‑squared:  0.09556  F‑statistic:  5.12  on 1 and 38  OF,  p‑value:  0.02944 

図 19.Rの出力結果（付録 Excel、シート 6)

布である。二項分布は、正規分布がμと 8の2つ のパラメータ、ポアソン分布が入の1つのパラ メータで定義されたのと同様に、 pという 1つの パラメータで定義される。すなわち、このpの値 により、様々な分布の形状をとることができ（図 14)、二項分布を用いる一般化線形モデルにおい ては、データに最も当てはまるようなpの値をも

とめることが目的である。

しかしこのパラメータ pもどんな値でもよいと いうわけでなく、 O lの範囲の値しかとれない という制限がある。そのため、この場合も単純に、

p=μ+ax 

のような線形式としてしまうと、 pがO lの範 囲外の値をとりうることになってしまう。そのた め、二項分布を用いた一般線形モデルでは次のよ

うなロジットリンク関数を用いる。

p=l/e ^{（ー(μ+ax))} 

このようにすればpが0 1の範囲の値しかと らない（図20)。

Rでは、以下のコマンドで一般化線形モデル の パ ラ メ ー タ 推 定 を す る こ と が で き る （ 付 録 Excel、シート 7)。ここでは飼養施設が従来のも のである対照群と設備の配箇を改良した処理群と で、全頭数のうち敵対行動を示した頭数の割合に ついて検討することを考える（帰無仮説は「試験 処理に効果はない」である）。 T列がTlのもの が対照群のデータ、 T2が処理群のデータで、 A 列が敵対行動を示した頭数、 B列が示さなかった 頭数である。

1. データを "dat7"として読み込む。

dat7 

< ‑

read.table("clipboard", header= T) 

多田・新村

ー10 ‑5 

゜

5 

10  μ+ax 

図20.p = 1 / e< ‑<μ+ a x))としたロジットリンク関数 pは常に0 1の範囲の値をとる。

2.  "dat7"のデータのA/(A+B)を目的変数とし、

説明変数としてTを含め、ロジットリンク 関数を用いて誤差分布に二項分布を指定した 一般化線形モデルを "fm7a"とする。先ほど のポアソン分布や正規分布での glm関数と 同様に、モデルの構造 (cbind(A,B) T)、用 いるデータ (dat7)および分布とリンク関数 (family=binomial(link="logit")) を記載してい る。なお、前述のように二項分布では割合p で起こる事象を、この場合、 A+B回観測し た場合の分布を表すので、 AおよびBそれ ぞれの値が計算に反映されるように、今回は このようにcbindコマンドを用いた特殊な表

記をしている。また、 familyのbinomialは ポアソン分布を、 linkのlogitはロジットリ ンク関数を指定したことを示している。

fm7a<‑glm(cbind(A, B)‑T, data=dat7,  family=binomial(link="logit")) 

3.  "fm7a"の要約を表示する（図21)。 summary(fm7a) 

アウトプットの(Intercept)がμ をTT2がaを示 している。各群の敵対行動を示した頭数の割合p はロジットリンク関数にこれらの値を代入する ことで求められる。すなわち、対照群ではp= 1  I e―1.1695 

= 

0.76、設備の配置を改良した処理群で

> dat7  <‑ read.table("clipboard",  header= T) 

> fm7a<‑glm(cbind(A,  B) T,  data=dat7,  family=binomial(link="logit")) 

> summary(fm7a)  Call: 

glm(formula = cbind(A,  B)   T,  family= binomial(link = "logit"),  data= dat7) 

Deviance Residuals: 

Min  lQ  Median  3Q  Max 

‑2.27702  ‑0.54876  0.08337  0.59563  2.94172  Coefficients: 

Estimate Std.  Error z value Pr(>lzl)  (Intercept)  1.1695  0 .1054  11. 097  <2e‑16  *** 

TT2  ‑1. 2529  0 .1395  ‑8. 983  <2e‑16 *** 

Signif.  codes:  0'***'0.001'**'0.01'*'0.05'.'0.1''1  (Dispersion parameter for binomial  family taken to be 1) 

Null deviance:  185.58  on 99.  degrees of freedom  Residual deviance:  100.23  on 98  degrees  of freedom  AIC:  345.43 

Number of Fisher Scoring iterations:  4 

図21.Rの出力結果（付録Excel、シート 7その1)

対照群 処 理 群 0.4 

0.3  2 1  

n i o  

赳薯 p = 0.76 

p= 0.48 

゜

^{0  0.2  0}観察割合^{.4  0.6  0.8  1  0  0.2  0}観察割合^{.4  0.6  0.8}  図22.データヘの二項分布の当てはめ

バーがデータの頻度、実線が二項分布（仮にN = 10としたもの）を示す。対照群ではp= 1  le―1.1695=0.76、処理群ではp=1/e―(1.1695+ (‑1.2529))=0.48と推定された。

はP= 1 / e,‑,1.1695 ‑1.2s29 ))= 0.4Sとなる。対照群と 処理群それぞれについて、推定したパラメータか ら得られた二項分布とデータをあわせたものを図 22に示した。それぞれ二項分布がデータに当て はまっていることが分かる。

有意差検定についても、基本的には一般線形モ デルと同様の考え方で行なう。すなわち、検討し たい要因を含んだモデルと含まないモデルでデー タ に 対 す る 当 て は ま り 具 合 に 差 が あ る か を 検 討 する。ポアソン分布の一般 線 形 モ デ ル と 同 様 に Devianceを当てはまりの指標とする。同様にx

二乗検定により P値を計算することができる。

4.  "dat7"のデータの A/(A+B)を従属変数とし、

説明変数を含めず、ロジットリンク関数を用 いて誤差分布に二項分布を指定した一般化線 形モデルを "fm7b"とする。

fm7b 

< ‑

glm(cbind(A,  8)1, data=dat7,  family=binomial(link="logit")) 

5.  X二乗分布を用いた P値を算出する (fm7a とfm7bのデータヘの当てはまりを x二 乗 検 定で比較；図23)。

anova(fm7a, fm7b, test=℃ hi") 

x二乗検定による P値が、 2.2e‑16 (= 2.2 X 1/ 

(10の16乗））未満であり、非常に小さい値と算 出された。これは設備配置の改良に効果がないと す る な ら ば 、 今 回 デ ー タ か ら 得 ら れ た Deviance の 値 が 得 ら れ る の は0.001%ほ ど も な い 、 す な わち、ほぼありえないことを意味する。よって、

5% (0.05)を基準とすると、帰無仮説が棄却され、

対立仮説「試験処理に効果がある」が採択される。 すなわち、設備の配置を変えた処理区では敵対行 動を示す個体の頭数割合が低下するという仮説を 支持することができる。

まとめ

分散分析、回帰分析および共分散分析はいずれ も正規分布を仮定し、平均値と処理の効果と誤差 との足し算の関係で表す、一般線形モデルという 一つの枠組みに集約される。この一般線形モデル

を拡張し、リンク関数を利用して正規分布でない 分布（ポアソン分布や二項分布など）も扱えるよ うにしたものが一般化線形モデルである。そして 一般化線形モデルにおける有意差検定は仮定した 分布が何であっても、要因を含んだモデルと含ま ないモデルとでデータに対する当てはまりに差が あるかを検討するものとして捉えることができ る。多種多様に感じられる統計手法も、分布を仮

> fm7b <‑glm(cbind(A, B) l,  data=dat7, family=binomial(link="logit")) 

> anova(fm7a, fm7b, test="Chi")  Analysis of Deviance Table  Model  1:  cbind(A,  B)   T  Model 2: cbind(A,  B)   1 

Resid.  Of Resid.  Dev Df Deviance  Pr(>Chi)  1  98  100.23 

2  99  185.59 ‑1  ‑85.354 < 2.2e‑16 *** 

Signif. codes:  0 '***'0.001 '**'0.01'*'0.05'.' 0.1''1  図23.Rの出力結果（付録Excel、シート 7その2)

多田・新村

定するもの、すなわち、パラメトリックな手法は このように多くの部分が共通している。言い換え れば、碁本的な考え方さえ理解していれば、これ らの手法を適宜用いることは容易であると思われ る。本稿で解説してきたようにRを利用する場 合にはコマンドの一部を変更するだけで、モデル に含める要因、リンク関数や仮定する分布、また、

検定について適宜選択することができるので、ぜ ひ自身のデータにふさわしい統計モデルの構築を 行なっていただきたい。

おわりに

本稿ではRを用いながら、分散分析や回帰分 析を含む一般線形モデルの概要、そして、これを 拡張した一般化線形モデルを用いたカウントデー タや割合データの解析の流れを解説した。本稿に より一般線形モデルおよび一般化線形モデルの理 解が深まり、読者のみなさまの実際のデータ解析 に役立てていただければ幸いである。さらに遡れ ば、これらの手法を用いるにあたっては、実験前 に自然と統計モデルの構造を頭に置き、また、デー タの誤差分布を把握するのに充分なデータ数を考 慮することになるため、実験の組み立てがより洗 練されたものになることも期待される。

また、先に書いたことと矛盾するようである が、カウントデータであってもポアソン分布に従 わない、割合データであっても二項分布に従わな い、また、そもそもこれまで検討してきた分布の いずれにも従わないデータはいくらでも実在し、

本稿の内容だけではデータ解析の実践に不足であ ることは充分に考えられる。ポアソン分布や二項 分布を用いた一般化線形モデルで分布の当てはま りを確認する際の一つの目安としては、モデル の出力結果に表示される Residualdevianceをそ の自由度で除して算出する Dispersionparameter 

(VenablesとRipley2002)がある(1から乖離する ほど適合していない；図16の例の場合、 94.842 I 98 = 0.968)。モデル予測値および実際のデータ のプロットのチェックや、 Dispersionparameter  が大きい場合の過分散検定（ポアソン分布モデ ル で はAERパ ッ ケ ー ジ のdispersiontest関数：

KleiberとZeileis2008)により、用いた分布が 不適切と判断される場合には、ここでは述べな

かった手法である擬似ポアソン分布や擬似二項 分布の利用、その他、ガンマ分布や負の二項分 布といった分布を用いること、さらには、一般 化線形モデルをさらに拡張して個体差などの変 量効果も扱うようにした一般化線形混合モデル (Generalized linear mixed model: GLMM) の 利 用 (Rでは主にlme4パッケージ： Batesら2015) など、いくつかの対処方法が考えられる。また、

その他、本稿では特に触れなかった交互作用を含 んだモデルの検討、群間の多重比較 (Rでは主に multcompパッケージ： Hothornら2008) など、

各々が必要とする統計手法は様々であると思う。

幸運なことに、インターネット検索すると、これ らの統計手法を Rで扱うプログラムがインター ネット上に数多く公開されていることが分かる。

そのマニュアルや解説記事を読みながらプログラ ム例を実践してみることでそれらの手法に触れる ことはかなり容易であり、ぜひ試行してみていた だきたい。読者のみなさまの体験をもとに、本稿 の内容で触れなかったが加えるべき内容、また、

本稿で分かりづらかった点、さらに説明すべき項 目などについてのご意見、また、著者ら自身も統 計の専門家ではないため、本稿の内容で誤った点 のご指摘など、あればぜひご連絡をいただきたい。

今後の記事作成の参考にさせていただきたい。

参考文献

Bates D, Maechler M, Bolker B, Walker S. 2015.  Fitting Linear Mixed‑Effects Models Using  lme4. Journal of Statistical Software 67, 1‑48.  Hothorn T, Bretz F and Westfall P.  2008. 

Simultaneous Inference in  General Parametric  Models. Biometrical Journal 50, 346‑363. 

Kleiber C, Zeileis A. 2008. Applied Econometrics  with R. Springer‑Verlag, New York. 

R Core Team. 2015.  R:  A Language and  Environment for  Statistical  Computing  [homepage on the  Internet].  R Foundation  for  Statistical  Computing, Vienna, Austria.  [cited 27 January 2016]. Available from URL: 

https:/ /www.R‑project.org 

Venables WN, Ripley BD. 2002. Modern Applied  Statistics with S. 4th edn. Springer, New York. 

Progress of statistics in applied animal behaviour science (2):

Generalized linear model

Shingo Tadau, Tsuyoshi Shimmura2*

1 Dairy Production Research Division, NARO Hokkaido Agricultural Research Center, Sapporo, Hokkaido, 062-8555, Japan

2 Institute of Agriculture, Tokyo University of Agriculture and Technology, Fuchu, Tokyo, 183- 8538, Japan

*Corresponding author. E-mail address: [email protected] (S. Tada),

[email protected] (T. Shimmura)

Summary

Recently, generalized linear model have become widely used even in the research :fields of applied animal behaviour science and livestock management. The generalized linear model is the extension of the general linear model and deal with various distribution including normal distribution. In the previous paper, we introduced the concepts of general linear model and showed the examples of analysis with sample data using free software "R". In this manuscript, we introduce the concepts of generalized linear model and then show the examples of analysis with sample data and practical commands using "R".

Keywords: behaviour, distribution, general linear model, generalized linear model, R

Animal Behaviour and Management, 53 (3): 117-126, 2017 (Received 24 February 2016; Accepted for publication 17 May 2017)