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

(1)

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

誌名

誌名 Animal behaviour and management ISSN

ISSN 18802133 著者

著者多田, 慎吾

新村, 毅巻/号

巻/号 53巻3号

掲載ページ

掲載ページ p. 105-116 発行年月

発行年月 2017年9月

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

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

(2)

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

多田慎吾1* ・新村毅

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

0 6 2 ‑ 8 5 5 5

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

1 8 3 ‑ 8 5 3 8

℃

o r r e s p o n d i n g a u t h o r . E ‑ m a i l a d d r e s s : t d s h n g @ a f f r c . g o . j p

(多田慎吾），

s h i m m u r a @ g o . t u a t . a c . j p

(新村毅）

要約

応用動物行動学および家畜管理学分野では、取得したデータが正規分布に従わないことが多々ある。

この場合には、正規分布を仮定した分散分析や回帰分析といった一般線形モデルを用いるのは適切でないため、主にデータの変数変換やノンパラメトリック手法の適用により対処されてきた。しかし、変数変換によってもデータの正規性を満たせない場合も多く、ノンパラメトリック手法は検出力が低いなどの問題も指摘されている。近年、これらの代替となる可能性のある手法として、一般線形モデルを拡張して、正規分布以外の分布も仮定できるようにした一般化線形モデルが普及してきており、応用動物行動学および家畜管理学分野でも一層の利用が期待される。本稿では、前編として、フリーのソフトウェアである R を用いて、まず一般線形モデルについて、その概念と練習用データを用いた解析例を、実際の Rのコマンドを示しつつ解説する。

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

A n i m a l B e h a v i o u r a n d M a n a g e m e n t , 5 3 ( 3 ) : 1 0 5 ‑ 1 1 6 , 2017 ( 2 0 1 6 . 2 . 24

受付；

2 0 1 7 . 5 . 24

受理）

はじめに

研究に手をつけた人が苦労してやっとのことで充分なデータを取り、さあ、取りまとめという段階になって直面する問題の一つは、統計解析ではないだろうか。そこで多くの人は、先輩から教えてもらうなり文献を参考するなりして、分散分析や回帰分析といった手法を学ぶことと思う。これらの手法は、後述するように正規分布を仮定するもので、まとめて一般線形モデル

( G e n e r a l l i n e a r m o d e l )

と呼ばれる。そして、これらの手法を習得し、使いこなして意気込んで論文を書き上げ、投稿した段階で突き当たる問題が、審査員からの次のような指摘である。「データの正規性は確認したか？」、「データを変数変換しなさい。」、

「パラメトリック検定ではなくノンパラメトリック検定を使うべきだ。」・・・。ノンパラメトリック検定を使うと、今度は「ノンパラメトリック検定だと多くの情報を削ることになる。」という指摘を受けることもある・・・。こうやって、指摘の無限ループにはまり、うんざりした人も多いのではないだろうか。しかも、パラメトリック検定

とノンパラメトリック検定という古典的な考え方では、決してこの無限ループから抜け出すことはできない。これまで、著者らも応用動物行動学および家畜管理学関連の学術誌に論文を掲載しているが、投稿した際にこういった指摘を受けてきている。というのも、応用動物行動学および家畜管理学分野で扱うデータは、正規分布を仮定した一般線形モデルをそのまま当てはめづらい性質のデータを頻繁に扱う学問だからだ。それらの指摘を解決する手法として、一般化線形モデル

( G e n e r a l i z e d l i n e a r m o d e l )

を利用できるのではないかということをモチベーションとして、著者らは

2013

年

3

月に「動物行動学における統計手法の発展：

GLM

と

GLMM

」と題した勉強会を開催し、この手法をもっと私たちの研究分野でも普及すべきではないかという感触を得た。しかし、

以上のことを初学者に分かりやすく、実例を示しながら説明した適当な日本語の文献が見当たらなかったので、勉強会の内容を踏まえ、著したのが本解説記事である。

本稿は前編として、まず世界的にも有名な無料解析ソフト「

R

」を紹介し

( RCore Team

(3)

応用動物行動学におけるGLM

2 0 1 5 )

、一般線形モデルとは何かを説明したもの

である。解説記事後編でいよいよ紹介する一般線形モデルを拡張した一般化線形モデルはとは何か、理解する導入となることを意図している。そのまま使用できるRのコマンドや練習用の実データも併せて掲載しているので、読者のみなさまも実際に手を動かしながら読み進めてもらえると理解もより進むと思われ、実際に自分のデータを解析する際にも応用できるはずである。本解説記事を通して、統計初学者が一般線形モデルや一般化線形モデルといった統計手法を身に付け、実際のデータ解析に活用し、また、使うべき統計手法を見据えて実験を実施できるようになる一助になれば幸いである。

Rについて

Rは、統計解析や作図に用いることができるソフトウェアである

( H o r n i k ,2 0 1 5 )

。ライセンス登録のような特別な手続きをすることもなく、無料でダウンロードして使用できる上に、

Windows

、

Mac

や

L i n u x

といった様々な

OS

に対応しているため、導入が容易である。また、拡張性にも優れ、

CRAN

と呼ばれるインターネットサイト

( h t t p s : / / c r a n . r ‑ p r o j e c t . o r g / )

^{から適宜必} 要なパッケージ（解析用の関数、プログラムの集まり）をインストールすることで、大抵の解析手法は

R

で実行可能である。

R

本体やこれらパッケージのソースコードは公開されており、世界中のユーザーに実際に使用されながら発展・洗練され続けている。実際、日本語でRのコマンド例

を示し、解説しているサイトも多数ある。このような背景があるため信頼性が高く、最近では、

A p p l i e d A n i m a l B e h a v i o u r S c i e n c e

など応用動物行動学の国際誌でR を用いて解析している学術論文も数多くみられる

( e . g .

一般化線形（混合）

モデル：

K i d d i e

と

C o l l i n s2 0 1 5 ; S t r a t m a n n

ら

2 0 1 5 ;

ベイズ法：

P e d e r s e n

ら

2 0 1 4 ;

ソーシャルネットワーク解析：

B o y l a n d

ら

2 0 1 6 )

。基本的に操作はコマンド入力により行なうので、

E x c e l

のように「セルをクリックして入力」、「ドラッグ＆

ドロップでデータ選択」といった直感的な操作のできない点がとっつきづらく、最初は何をどうすればよいか分からないかもしれない。しかし、本稿も含め分かりやすい解説も多く出ているので、

最初はコピー＆ペーストでもコマンドの実行例と出力結果をつき合せて進めていけば、 Rで何をやっていて何ができるのか自ずと把握できていく

と思う。

Rのインストール方法とデータの読み込み

「 TheR P r o j e c t f o r S t a t i s t i c a l C o m p u t i n g

」というインターネットサイトのダウンロードページ

(URL: h t t p : / / c r a n . r ‑ p r o j e c t . o r g / m i r r o r s . h t m l )

にはRをダウンロードできる各国のミラーサイトが一覧になっている。日本からダウンロードするのであれば、

J a p a n

の

TheI n s t i t u t e o f S t a t i s t i c a l M a t h e m a t i c s

^{(統計数理研究所、}

URL:

h t t p s : / / c r a n . i s m . a c . j p / )

を選択すれば差し支えない。あとは自身の用いる

OS

を選択し、

Windows

なら

b a s e

ファイル、

Mac

なら

pkg

ファイルをダウンロードして、ダウンロードされたファイルを実行すればインストールされる。セットアップ時の各設定も初期のままで問題ない。

インストールが成功して Rを起動すると図1 の画面となる。まず、解析の第一段階としてデー

タの読み込みを行なってみる。著者は、

E x c e l

に実験データを入力しておき、そのデータを

E x c e l

上でコピーし、それをクリップボード経由でR に読み込む方法をとっているので、ここではその方法を紹介する。他にも、コマンド入力でデータテーブルを作成する、テキストファイルや

csv

ファイルを直接読み込む、といった方法もあるので、興味のある方はそちらを選択してみてもよいかもしれない。クリップボードを用いる方法は具体的には以下の手順に従う。

1 . E x c e l

シート上に実験データを用意し、

1

行目をタイトル行とする（図2)。ここでは、

処理区

( T r e a t m e n t )

が

2

つ

( A

と

B )

存在し、

各区での行動時間（例えば

1

日あたりの摂食行動時間（時間））のデータがyの列に入力

されているものと仮定する。

2 . E x c e l

シート上のデータテーブルを選択しコピーする（図2)。これ以降、データを読み込むまでコピーは行なわないようにする。

3 . R

上で以下のコマンドを入力する（キーボードで以下の通りに文字を打つ；図3)。

"header= T"

で、一行目

( T r e a t m e n t

と

y )

が、名前（変数名）であることを指定している。

d a t < ‑ r e a d . t a b l e ( " c l i p b o a r d " , header= T ) 4 .

キーボードの

" E n t e r "

を押す。これで

R

に

d a t

という名でデータテーブルを読み込めた。

d a t

の部分は任意の文字列に変更しても問題ない

( d a t a

、

d a t 2

など）。

5 . d a t

と入力しキーボードの

" E n t e r "

を押すことで、

d a t

の内容を確認することができる（図

4 )

。これで、

E x c e l

上でコピーしたデータが Rに読み込まれたことになる。

ちなみに、 R上でキーボードの ↑ "を押すと以前に入力したコマンドが表示される (1度押すと

1

つ前のもの、

2

度押すと

2

つ前のもの）。これを利用すれば以前に入力したコマンドをいちいちキーボードで入力する手間が省ける。そのため、

今回の

d a t < ‑ r e a d . t a b l e ( " c l i p b o a r d " , header= T )

のように、よく使うコマンドはあらかじめテキス

(4)

その他 "Hッケージウインドウヘルプへ

R ver'3ion 3. 2. 0 (2015‑04‑16) ‑‑"Full of In,;;redient,i"

、

^Platform^:^xB6̲64^‑^w64‑mio^叩³²^/x64 (64‑bit)

1 ・ ? e t l ! i i ? t ' J f i 彎彎，聾;" " . . .

^{と入力して}^く^ださい。

Rは多くの貢猷者による共同ブロジIり卜です。詳しくは 'contributor" ()'と入力してください。

また、Rや R(J)バッケージを出版物で引用する際の形式については 'citation()'と入力してください。

'demo()'と入力すればデモをみることができます。

'help()'とすればオンラインf¥ルブが出ます。

'help.'3tart ()'で HTMLブラウザによるf¥ルブがみられます。

' • q()'と入力すれば Rを終了します。

>I

JJI,

‑'‑,,. I'"" n

図1. Rの初期画面

ここにコマンドを入力して様々な作業を行なう。この画面は Rのバージョン 3.2.0のものだがバージョンはインストール時のものに準ずる。

4、A E F G H

ぅ ―

‑‑+

― ‑ ‑

—• —• 一・一

2 基

!.l

12 13

I !

~-—----—- ‑‑ ‑ ‑ ‑

,9

4 ー

•. ‑‑・‑‑‑‑;‑‑‑‑‑‑‑‑‑c‑‑‑ ‑‑‑・‑‑‑'

15 ‑1 ‑ ‑ . ‑ ‑ ‑

--—

図2. Excelシート上の実験データ

このように読み込みたいデータを選択する。なお、ここにでは説明変数を"Treatment"、従属変数を"y"、また、処理名をAおよび Bとしたが、これらは任意の文字列としても構わない。

(5)

応用動物行動学における GLM

R ver,iion 3.2.0 (2015‑04‑16) ‑‑"Full of Ingredient:,"

Rは、自由なソフトウIアであり、「完全に燕保証」です。

一定の条件に従えば、自由に,:;ti,を再翫布することができます。

配布条件の詳細（こ閲しては、 •licen3e ()• あるいは 'licence()• と入力してください。

R(さ多くの貢猷者による共同ブロジIり卜です。

詳しくは 'contributor" ()• と入力してください。

また、Rや Rのバッケージを出版物で引用する際の形式については

•citation()' と入力してください。

'demo()• と入力すればデモをみることができます。

'help()'とすればオンラインf¥ルブが出ます。

•help. 3tart ()• でHTMLブラウザによるf¥ルブがみられます。

'q(). と入力すれば Rを終了します。

> dat← ごead.table("clュpboaごd",headeご ‑T)

I

Ill

図 3. Rへのコマンドの入力

その他バッケージウインドウヘルブ

1

_Rver,sion 3. 2. 0 (2015‑0'!‑16) -- • 和11 of Ingredienc,s•

Rは、自由なソフトウ1アであな「完全に燕保証」です。

一定の条件に従えば、自由にこれを再翫布することができます。

配布条件の詳細（こ閲しては、 •licen!le ()'あるいは 'licence()• と入力してください。

I

R (さ多くの貢猷者による共同ブロジIり卜です。詳しくは 'contributor!l ()• と入力してください。

また、Rや Rのバッケージを出版物で引用する際の形式については

・citation()'と入力してください。

'd<emo ()'と入力すればデモをみる改：ができます。 'help()'とすればオンうインヘルブが出ます．

• help. !ltart ()• でi!TMLブラウザによるf¥ルブがみられます。

'q (). と入力すれば Rを終了します。

> dat <‑read.table("clュpboaごd",headeごa T)

> dat

Treatment y

1 A 7.6

2 A 6.9

3 A 5.0

1 A 5.9

5 A 9.5

6 B 9.0

7 B 6.1

8 B 9.3

，

3 9.9

"Ia

図4.読み込んだデータの表示と確認

Excelシート上のデータテーブル（図 2)が読み込めている。

(6)

(a) (b)

赳堰

平均値平均値

図5.正規分布の形状（分散は

a<

b)

平均値は等しくても、分散が小さいほど平均値に近い値のデータが得られる頻度（確率）が高く(a)、分散が大きいほど平均値にから離れた値のデータが得られる頻度（確率）が高い (b)。トファイルや

E x c e l

ファイルなどに入力して保存

しておき、用いる際に R上にいったんコマンドを貼り付けておくと便利である。

一般線形モデルとは？

まず、統計モデルとはなんであろう？いろいろな統計モデルがあるかもしれないが、よい統計モデルの共通点として挙げられるのは、観測されたデータを精度よく説明（予測）できることであると思う。そして先にも述べたが、分散分析や回帰分析、加えて共分散分析といった手法は、一般線形モデルという枠組みに一括りにすることができ

る。というのは、これらの手法には、データの誤差が正規分布に従うことを前提としている、データを足し算の関係で表している、という共通点があるためである。一般線形モデルは、これらの枠組みのもとで、データを精度よく予測しようとす

るものである。

一つ目の「誤差が正規分布に従う」は、具体的には図

5

のように表される。データが、平均値の周囲に釣鐘状に分布している。このデータのばらつき具合は、分散と呼ばれるパラメータで表すことができる。分散が小さいと平均値の周囲にデータが集中し（図

5 a )

、分散が大きいと裾が広い分布となる（図5b)。一般線形モデルの基本的な考え方は、データにはそのデータを表す真の平均値が存在するが、実際に観測により得られたデータは様々な誤差を含むため、このように平均値の周囲にばらついているというものである。

二つ目の「データを足し算の関係で表している」

は、直感的にはわかりづらいかもしれない。ここでは分かりやすいように、真の平均値が

5

であることが分かっている測定対象データを3回観察した場合を考えてみる。その結果は以下の通りであった。

4 . 2

、5

. 3

、5

. 9

全てのデータが

5

に近いようではあるが、丁度

5

の値であることはなかった。これは先に述べた通り、観測データに誤差が含まれたためであると考える。これらのデータを足し算の関係で表すとは、

具体的には以下のようなことである。

4 . 2

=

5

+

( ‑ 0 . 8 ) 5 . 3

=

5

+

0 . 3 5 . 9

=

5

+

0 . 9

すなわち、データは真の平均値に誤差の値を加えたものとして表すことができる。記号を用いて表すと以下のように書ける。

y=μ+s

ここでμ(ミュー）が平均値、 E (イプシロン）

が正規分布に従う誤差である。この正規分布の分散を ^(J (シグマ）とすると、一般線形モデルの構築とは、観測したデータに最も当てはまるようなμおよび6の値を求めることである。簡単なデータセットを用いて Rで一般線形モデルの構築をやってみることにする（付録

E x c e l

、シート 1)。先ほどと同様にして、練習用データのシート 1のデータをコピーし、 R上で以下のコマンドを打つ。なお練習用データが記載された

E x c e l

ファイルは、本解説が掲載されている J‑STAGEの

Web

ページ

(URL:h t t p s : / / w w w . j s t a g e . j s t . g o . j p / b r o w s e / ‑ c h a r / j a / )

もしくは応用動物行動学会の

Web

ページ

(URL:h t t p : / / w w w . j s a a b . o r g / )

に掲載されているので、ダウンロードしてご利用いただきたしヽ。

1. データを

" d a t l "

として読み込む。

dat1

< ‑

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

2 . " d a t l "

のデータを用いて

y=μ+

Eとした

（平均値と誤差のみの）一般線形モデルを

" f m l "

とする。ここでは一般線形モデル用の関数Imを使って下記ように入力する。関数 Im()のカッコ内には、初めにモデルの構造（今回の場合Y

=μ+ s )

を「

y ‑ 1

」のよう

に記載し、次に使用するデータ（今回の場合

(7)

応用動物行動学における

GLM

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

> fml <‑lm(y l, datl)

> summary(fml) Call:

lm(formula = y 1, data= datl) Residuals:

Min lQ Median 3Q Max

‑2.7540 ‑1.4915 ‑0.0990 0.9635 4.4760 Coefficients:

Estimate Std. Error t value Pr(>ltl) (Intercept) 5.0540 0.4156 12.16 2.06e‑10 ***

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

図6. Rの出力結果（付録Excel、シート 1) datl)を記載する。「y‑1」の「y‑」はy=ま

でに相当し、「

1

」が直感的に分かりづらいが、

この表記がμを表すことになっている。また、

関数

I m

は一般線形モデル、すなわち、正規分布に従う誤差をもつモデルを扱う関数であり、 eについてはモデル構造として入力することなく、計算されるようになっている。なお、各関数の記載方法などの詳細は、

I m

の場合「help(lm)」と入力することで表示され

るので参照されたい。

fm1

< ‑

lm(y‑1, dat1) 3. "fml"の要約を表示する。

summary(fm1)

このようにコマンドを打つと、図

6

のような結果が表示される。図中の(Intercept)のEstimate が平均値μ、一番下段のResidualstandard error が標準偏差

(J c r )

を示している。すなわち、このデータは平均値5.0540、分散3.45 (1.858りの正規分布である母集団から得られたものとモデリ

ングできたことになる。

次にもう少し複雑なケースを考えてみる。ある処理を加えると、値が

2

増加する（すなわち、処理の効果が+2) とする。処理を加えて測定対象データを 3回観察した結果は、以下の通りであった。

対照群： 5.4、5.3、4.9 処理群： 7.2、6.8、7.5

これらについても、

5.4 = 5 + 2 X ̲Q + 0.4 5.3 = 5 + 2 X ̲Q + 0.3 4.9 =

5

+

2

X ̲Q + (‑0.1) 7.2 = 5 + 2 X

l

⁺⁰^.²

6.8 = 5 + 2 X

l

⁺⁽^‑⁰^.²⁾

7.5 = 5 + 2 X 1 + 0.5

このように、下線部は対照群なら〇、処理群なら

1

をとるとすれば、このような足し算の関係で表

すことができる。すなわち、データは真の平均値に処理の効果を加え、さらに誤差の値を加えたものとして表すことができる。記号を用いて表すと、

以下のように書ける。

y=μ+ax+s (対照群の場合x=O、処理群の場合x=l)

ここでxが、今のように試験処理のようなカテゴリー分けするような性質のもの（カテゴリー変数）

である場合は、この手法は分散分析と呼ばれ、 X

が数値で示される連続値であるならば、回帰分析と呼ばれる。さらに説明変数（ここでいう x)の種類が2つ以上あり (y=μ+ nx1 + Px2 + "fX3 + ... +

s )

、それらがいずれも連続値であるものが重回帰分析、説明変数にカテゴリー変数と連続値の両方を含むものが共分散分析である。そして、以上のいずれの手法でも共通して行なっているのはデータに最も当てはまるようなμ、 a、6を求めることである（データを最も良く説明するμ、 a、

6の値を計算する）。これらのパラメータを推定することで、精度のよい予測のできるモデルをつくることができる。実際に、 R上でこれらのパラメータ推定を行なってみる。試験処理が

2

つである正規分布に従うデータ例を示す。このデータで先頭行が xである列が"control"となっているものが対照群、 "treatment"となっているものが処理群である（付録Excel、シート 2)。これらのデータに線形モデルを当てはめると次のようになる。

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

dat2

< ‑

^rêâ^d^.^tâ^b^lê⁽^"^c^lⁱ^p^bôâ^r^d^"^,^header= T⁾

2. "dat2"のデータを用いて、説明変数としてカテゴリー変数xを含めた一般線形モデルを

"fm2"とする。ここで「yx」がモデルの構造（今回の場合y=μ+ ax + s)を表しており、説明変数として xを指定していることになる。なお、今回のように説明変数を

1

つ以上含む場合、前述の練習用データのシート

1

(8)

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

> fm2 <‑ lm(y x, dat2)

> summary(fm2) Call:

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

‑5 .1510 ‑2. 0535 0 .1055 1. 7630 4. 8305 Coefficients:

Estimate Std. Error t value Pr(>ltl) (Intercept) 5.0495 0.5875 8.595 l.92e‑10 ***

xtreatment 2.0215 0.8308 2.433 0.0198 *

Multiple R‑squared: 0.1348, Adjusted R‑squared: 0.112 F‑statistic: 5.92 on 1 and 38 DF, p‑value: 0.01978

図7.Rの出力結果（付録Excel、シート 2) での例では入力していた、μ を表す「

1

」は

省略することができる（もちろん、省略せずに「y l+x」としてもよく、その場合も出力結果は同じである）。また、説明変数内のアルファベット順の最も早いものが対照群とされ、今の場合は説明変数xのcontrolがこれにあたる。

fm2

< ‑

lm(yx, dat2) 3. "fm2"の要約を表示する。

summary(fm2)

出力結果は図7のようになり (Intercept)の Estimate (5.0495) が平均値μ、 xtreatmentの Estimate (2.0215) が処理の効果a、Residual standard error (2.627)が標準偏差

(J

cr) を示

している。 controlの平均値がμ に、 treatmentの

平均値がμに aを加えた値に一致する。すなわち、対照群のデータは平均値5.0495、分散6.90

(2.627りの正規分布である母集団から、処理群のデータは平均値7.0710、分散6.90の正規分布である母集団から、それぞれ得られたものとモデ

リングできたことになる。

次に回帰分析の一例も示す（付録Excel、シート3)。ここでは、例えばある動物について、体長

( x

列）と体重

( y

列）をそれぞれ入力したデータを考える。

dat3

< ‑

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

続値xを含めた一般線形モデルを "fm3"とす

> fm3 <‑ lm(y x, dat3)

> summary(fm3) Call:

‑4.0224 ‑1.1796 ‑0.2176 1.1602 4.9242 Coefficients:

Estimate Std. Error t value Pr(>ltl) (Intercept) 2.33924 0.91930 2.545 0.0151 *

X 1.91493 0.09832 19.477 <2e‑16 ***

Multiple R‑squared: 0.9089, Adjusted R‑squared: 0.9066 F‑statistic: 379.4 on 1 and 38 OF, p‑value: < 2.2e‑16

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

(9)

応用動物行動学における

GLM

る。

fm3

< ‑

lm(yx, dat3) 3. "fm3"の要約を表示する。

summary(fm3)

出力結果は図 8のようになり (Intercept)の Estimate (2.33924) とxのEstimate (1.91493) がそれぞれ回帰直線の切片と傾き、 Residual standard error (1.915) が標準偏差(..;

o ‑ )

を示している。すなわち、このデータでは体重は2.33924+ 体長 X 1.91493の値に分散3.67

(1.915りの正規分布に従う誤差を加えたものとしてモデリングできた。

一般線形モデルの有意差検定

本稿の趣旨からは逸れるかもしれないが、私たちの研究分野でおそらく求められるだろう有意差検定について、ここで解説する。なお、生態学分野などではモデル選択や、ベイズ推定によるパラメータの信頼区間算出といった有意差検定に代わりうる方法も提唱されているようなので (Garamszegiら2009)、興味のある方は自分の目的にふさわしい方法を適宜検討し、用いることを推奨したい。

有意差検定とは、帰無仮説を立て、それが棄却できるなら対立仮説を採用するというプロセスである。前項の試験処理の例に当てはめて考えてみる。この場合の帰無仮説は「試験処理に効果はない」であり、対立仮説は「試験処理に効果がある」

となる。この「試験処理に効果はない」と「試験処理に効果がある」を一般線形モデルで表すとそれぞれ次のようになる。

「試験処理に効果はない」：データは平均値と誤差だけで表される (y=μ+ i:)

「試験処理に効果がある」：データは平均値と処理効果と誤差で表される (y=μ+ax + s)

有意差検定は、分散分析であっても回帰分析であっても、含めた要因がそのモデルのデータヘの当てはまりの改善に寄与しているかという、一般線形モデルで統一した問題として考えることができる。今の場合、 y=μ+sのモデルにaxの項を加えた際のモデルの当てはまりの程度の差を考えることになる。一般線形モデルの当てはまりの指標の1つとして用いられるのがF値である。 F 値が大きいほどモデルの当てはまりの改善度合い

が大きいと考えて差し支えない。付録Excel、シート

4

のデータを用いて

F

値を算出してみる。

dat4

< ‑

^rêâ^d^.^tâ^b^lê⁽^"^c^lⁱ^p^bôâ^r^d^"^,^header =^T)

2. "dat4"のデータを用いて、説明変数として x を含めた一般線形モデルを "fm4a"、説明変数を含めない（平均値と誤差のみの）一般線形モデルを"fm4b"とする。

fm4a

< ‑

lm(yx, dat4) fm4b

< ‑

^l^m⁽^y^‑¹^,^d^a^t⁴⁾

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

この要約にある "F‑statistic" (5.12)がこのデータで説明変数として xを含めた場合の

F

値である。一方、平均値と誤差のみの一般線形モデルを

"fm4b"の要約は次のようになる。

4. "fm4b"の要約を表示する（図10)。 summary(fm4b)

こちらのモデルでは処理効果は存在せず、データ

> fm4a <‑ lm(y x, dat4)

> fm4b <‑ lm(y l, dat4)

> summary(fm4a) Call:

‑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 *

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

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

(10)

> summary(fm4b) Call:

lm(formula = y 1, data= dat4) Residuals:

‑6.3918 ‑2.2393 ‑0.3217 2.2957 5.2482 Coefficients:

Estimate Std. Error t value Pr(>ltl) (Intercept) 10.5818 0.4797 22.06 <2e‑16 ***

図10.Rの出力結果（付録Excel、シート4その2) は平均値 10.5818、分散9.21 (3.034りの正規分

布から得られるとしている。

帰無仮説を棄却する・しないの判断は次のように行なう。まず、「帰無仮説が正しい」として考えてみる。これはすなわち、先ほど構築した「試験処理に効果はない」としたモデルfm4bが正しく、このモデルからデータを適切に予測できるということを意味する。このモデルfm4bから予測されるデータに対して、「試験処理に効果はない」

モデルy=μ+Eと「試験処理に効果がある」モ

デルy=μ+ax+Eを当てはめ、

F

値を算出することを考えてみる。

5 .

帰無仮説が正しいとして得られるデータを乱数生成し、データ "dat4"のy2列として挿入する。

dat4$y2

< ‑

simulate(fm4b)[, 1]

6. dat4を表示する (y2列が追加されていることを確認）。

> fm4z <‑ (lm(y2 x, dat4))

> summary (fm4 z) Call:

lm(formula = y2 x, data= dat4) Residuals:

dat4

7 .

y2に関して、説明変数としてxを含めた一般線形モデルを"fm4z"とする。

fm4z

< ‑

(lm(y2x, dat4))

8. "fm4z"の要約を表示する（図11)。

summary(fm4z)

この要約にある "F‑statistic"が、帰無仮説が正しいとして今回生成されたデータでの

F

値である。

もう一度、 dat4$y2<‑simulate(fm4b)[,1]から以下の手順を繰り返すと、また異なる

F

値が得られる。

大抵はこのデータは帰無仮説が正しい「試験処理に効果はない」 fm4bモデルから得られたものなので、 y=μ+sでも y=μ+ax+sでもモデルの当てはまりに大差はないはず、すなわち、

F

値は小さいものになるはずである（図11では0.8068)。だが、たまたま処理区に割り当てられたデータに値が大きいデータ（もしくは値が小さいデータ）

が集中し、「試験処理に効果がある」モデルy=

‑5.4728 ‑1.9767 ‑0.0563 1.4466 6.1305 Coefficients:

Estimate Std. Error t value Pr(>ltl) (Intercept) 10.0879 0.5661 17.820 <2e‑16 ^＊＊＊

xtreatment 0.7191 0.8006 0.898 0.375

Multiple R‑squared: 0.02079, Adjusted R‑squared: ‑0.004978 F‑statistic: 0.8068 on 1 and 38 DF, p‑value: 0.3747

図11.Rの出力結果（付録Excel、シート 4その3)

この結果は乱数生成されたデータに対する解析結果なので実行の度に数値は異なる。

(11)

応用動物行動学におけるGLM

> n <‑10000

> FValue <‑numeric(n)

> for (i in 1:n) {

+ dat4$y2 <‑simulate (fm4b) [, 1] + fm4z <‑ (lm(y2 x, dat4))

+ FValue[i] <‑summary(fm4z)$fstatistic[l]

＋｝

> sum(FValue >= summary(fm4a)$fstatistic[l])/n [1] 0.0281

> hist(FValue, 50, freq=FALSE)

> abline(v = summary(fm4a)$fstatistic[l], lty = 2)

> anova(fm4a, fm4b, test="F") Analysis of Variance Table Model 1: y x

Model 2: y 1

Res.Df RSS Df Sum of Sq F Pr(>F) 1 38 316.30

2 39 358.93 ‑1 ‑42.622 5.1205 0.02944 *

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

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

μ+ax+i:: の当てはまりがよく、

F

値が大きくなることもある。実際のデータで得られた

F

値が、

このようにたまたま程度の頻度でしか発生しないものなのか、それとも、よくみられる程度の値なのかを判断する。具体的には、「帰無仮説が正しい」

モデルでデータを生成し、

F

値を算出するという過程を何度も繰り返すと、

F

値の頻度分布を得ることができ、 F値の発生確率 (P値）をもとめることができる。実際のデータで得られた

F

値より大きな値が得られる確率がとても低い場合（慣例的には多くの場合5%未満なら）、「試験処理に効果はない」とするならば、極端に珍しいことが起きたことになる。すなわち、「試験処理に効果はない」とは考えづらいので、帰無仮説は菓却され、対立仮説「試験処理に効果がある」が採択される。逆に、その

F

値が「試験処理に効果はない」としても充分な頻度で観測されるなら、帰無仮説は棄却されないということになる。この「帰無仮説が正しい」モデルから乱数生成でデータを得る過程を経た方法をパラメトリック・ブートストラップ法という。 Rでは以下の手順で実行できる。

9 .

以下の一連のコマンドでは上記のデータ乱数生成から

F

値算出までの作業を 10000回繰

り返し、

F

値の分布を得ている。

n <‑10000

FValue <‑numeric(n) for (i in 1 :n) {

dat4$y2 <‑simulate(fm4b)[, 1] fm4z <‑(lm(y2x, dat4))

FValue[i] <‑summary(fm4z)$fstatistic[1]

｝

10. 生成データから得た

F

値の分布において、

実際のデータから観測された

F

値より大きかったものの数をカウントし、繰り返し数

（今の場合 10000) で除して頻度を算出する

（これが P値：図 12中部の 0.0281)。 sum(FValue >= summary(fm4a)$fstatistic[11)/n

1 1 .

パラメトリック・ブートストラップ法によ

り得られた

F

値の分布を図示する。

hist(FValue, 50, freq=FALSE)

abline(v = summary(fm4a)$fstatistic[1 ], lty = 2) しかし、この方法はコンピュータの計算能力や時間を要するなど労力が大きいので、大抵の場合は既知である

F

分布を用いて

P

値を計算する。

F

値は、分子自由度が帰無仮説モデルと対立仮説モ

デルのパラメータ数の差、分母自由度がデータ数ー対立仮説モデルの推定パラメータ数の

F

分布に従う。

F

分布は次のコマンドで表示することがで

きる。

12. 分子自由度が 1、分母自由度が 38の

F

分布を表示する。

x <‑seq(O, 20, 0.1)

plot(x, df(x, 1, 38), type="n")

curve(df(x, 1, 38), add=T, from=O, to=20)

F

分布は、先ほどパラメトリック・ブートストラップ法でもとめた

F

値の頻度分布とほぼ一致する。

RでF検定を行なうプログラムももちろん用意されている。

(12)

1 3 . F

分布を用いた

P

値を算出する

( f m 4 a

と

fm4b

のデータヘの当てはまりを

F

検定で比較）。

a n o v a ( f m 4 a , f m 4 b , t e s t = " F " )

パラメトリック・ブートストラップ法から算出された値

( 0 . 0 2 8 1 )

と

F

検定による

P

値

( 0 . 0 2 9 4 4 )

がほぼ等しく、また低いことが分かる（図

1 2 )

。すなわち、今回実際のデータから得られた

F

値は「試験処理に効果はない」とするならば、

3%

程度の低頻度でしか観察されないものであり、多くの場合のように

5% ( 0 . 0 5 )

を基準とすると、

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

まとめ

本稿では、一般線形モデルではデータに正規分布を仮定し、平均値と処理の効果と誤差との足し算の関係で表すこと、また、この一般線形モデルにおける有意差検定の考え方について、 R を用いながら概説した。一般線形モデルの拡張である一般化線形モデルにおいても、多くの基本となる考え方は共通している。一方で、正規分布を前提とした一般線形モデルではなく、一般化線形モデルを用いるべき場面も確かに存在する。これらの共通点と相違点に注目しつつ、本解説記事後編、「応用動物行動学における統計解析の進展後編：一般化線形モデル」も読み進めていただければ幸いである。

参考文献

B o y l a n d NK, M l y n s k i DT, James R , Brent LJN, C r o f t DP. 2 0 1 6 . The s o c i a l n e t w o r k s t r u c t u r e o f a dynamic group o f dairy cows: From i n d i v i d u a l t o group l e v e l p a t t e r n s . A p p l i e d Animal B e h a v i o u r S c i e n c e 1 7 4 , 1 ‑ 1 0 .

Kiddie J , Collins L. 2015. Identifying e n v i r o n m e n t a l and management f a c t o r s t h a t may b e a s s o c i a t e d w i t h t h e q u a l i t y o f l i f e o f k e n n e l l e d d o g s ( C a n i s J a m i l i a r i s ) . A p p l i e d Animal B e h a v i o u r S c i e n c e 1 6 7 , 4 3 ‑ 5 5 .

Pedersen L J , H e r s k i n MS, Forkman B , Halekoh U, Kristensen KM, Jensen MB. 2014. How much i s e n o u g h ? The amount o f s t r a w n e c e s s a r y t o s a t i s f y pigs'need t o perform e x p l o r a t o r y b e h a v i o u r . A p p l i e d Animal B e h a v i o u r S c i e n c e 1 6 0 , 4 6 ‑ 5 5 .

R Core Team. 2015. R: A Language and Environment for S t a t i s t i c a l Computing [homepage on t h e I n t e r n e t ] . R Foundation f o r S t a t i s t i c a l Computing, V i e n n a , A u s t r i a . [ c i t e d 2 7 J a n u a r y 2 0 1 6 ] . A v a i l a b l e from URL:

h t t p s : / / w w w . R ‑ p r o j e c t . o r g

Stratmann A , F r o h l i c h EKF, Gebhardt‑Henrich

SG, Harlander‑Matauschek A, Wiirbel H ,

Toscano MJ. 2015. M o d i f i c a t i o n o f a v i a r y

d e s i g n r e d u c e s i n c i d e n c e o f f a l l s , c o l l i s i o n s a n d

k e e l b o n e damage i n l a y i n g h e n s . A p p l i e d Animal

B e h a v i o u r S c i e n c e 1 6 5 , 1 1 2 ‑ 1 2 3 .

(13)

fe.; )=IH!IJ!fWfi"

fJJ~

I: to It -Q G LM

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

General linear model

Shingo Tada1, 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

In the research fields of applied animal behaviour science and livestock management, we sometimes

meet the situation that normality of experimental data cannot be assumed. Since general linear model such as ANOVA or regression analysis assumes the normality, data transformation or non-parametric methods were conventionally applied for such data. However, the data transformation do not always ensure the normality of the transformed data, and statistical power of non-parametric methods become lower than those of parametric methods. As a possible alternative of these conventional methods, generalized linear model have become widely used. The generalized linear model is the extension of the general linear model and deal with various distribution including normal distribution.

In this manuscript,

we firstly introduce the concepts of general linear model and then show the examples of analysis with sample data and practical commands using free software "R".

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

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

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

誌名

誌名 Animal behaviour and management ISSN

ISSN 18802133 著者

著者 多田, 慎吾

新村, 毅 巻/号

巻/号 53巻3号

掲載ページ

掲載ページ p. 105-116 発行年月

発行年月 2017年9月

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

0 6 2 ‑ 8 5 5 5

1 8 3 ‑ 8 5 3 8

o r r e s p o n d i n g a u t h o r . E ‑ m a i l a d d r e s s : t d s h n g @ a f f r c . g o . j p

s h i m m u r a @ g o . t u a t . a c . j p

A n i m a l B e h a v i o u r a n d M a n a g e m e n t , 5 3 ( 3 ) : 1 0 5 ‑ 1 1 6 , 2017 ( 2 0 1 6 . 2 . 24

2 0 1 7 . 5 . 24

( G e n e r a l l i n e a r m o d e l )

( G e n e r a l i z e d l i n e a r m o d e l )

2013

3

GLM

GLMM

R

( RCore Team

2 0 1 5 )

( H o r n i k ,2 0 1 5 )

Windows

Mac

L i n u x

OS

CRAN

( h t t p s : / / c r a n . r ‑ p r o j e c t . o r g / )

R

R

A p p l i e d A n i m a l B e h a v i o u r S c i e n c e

( e . g .

K i d d i e

C o l l i n s2 0 1 5 ; S t r a t m a n n

2 0 1 5 ;

P e d e r s e n

2 0 1 4 ;

B o y l a n d

2 0 1 6 )

E x c e l

「 TheR P r o j e c t f o r S t a t i s t i c a l C o m p u t i n g

(URL: h t t p : / / c r a n . r ‑ p r o j e c t . o r g / m i r r o r s . h t m l )

J a p a n

TheI n s t i t u t e o f S t a t i s t i c a l M a t h e m a t i c s

URL:

h t t p s : / / c r a n . i s m . a c . j p / )

OS

Windows

b a s e

Mac

pkg

E x c e l

E x c e l

csv

1 . E x c e l

1

( T r e a t m e n t )

2

( A

B )

1

2 . E x c e l

3 . R

"header= T"

( T r e a t m e n t

y )

d a t < ‑ r e a d . t a b l e ( " c l i p b o a r d " , header= T ) 4 .

" E n t e r "

R

d a t

d a t

( d a t a

d a t 2

5 . d a t

著者多田, 慎吾

新村, 毅巻/号

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