untitled

動物生態学

生態学データの収集と解析

（Rを使った実習・入門編）

清田雅史

[email protected]

水産総合研究センター

国際水産資源研究所

1

ねらいと実習内容

• ヒトの手の大きさの雌雄差を題材とした実習：

▶

作業仮説を立てる

▶

サンプリング（データ収集）の立案

• 定義と測定要領を明確にする

▶

データの入力，チェック

• ここはExcelが便利

▶

Ｒを使って

• 基本のデータ処理 • 作図 • 統計検定 • モデル選択

１．研究立案∼データ収集

• 作業仮説を立てる

• どういう解析をするか予想しておく

作業仮説

• ヒトの手の大きさに雌雄差があるか？

必要なデータ

• 対象とする年齢範囲，地域，１個体の計測部位

解析方法

• データの定義，測定方法を明記することは大切

（測定誤差を小さくする）

この３つは （本来は） セットで考え るべきもの

データの種類

• 連続データ（小数や分数で表される実数）

1, 2.6, 2/3, 105.2, ‐0.0043

• カウントデータ（指折り数えられる自然数）

1, 3, 0, 245

• カテゴリカルデータ（名前）

名義： A, B, C... ;  オス，メス 順位あり：小，中，大

• 比率データ（0∼1の間の小数）

0.2, 3/4, 0.99

データの種類によって取り得る値の範囲や

分布型が異なり，適切なデータ解析法も異なる

扱いにくいので，なるべく 使わない方が良い

データの入力，チェック

• MS‐Excelを使うのが一般的

• 注意事項

• データ入力段階で表を美しくまとめるのは逆効果（後々使 いにくくなる） • セルの結合は使わない • 見出しはなるべく１行，名前の重複を避ける » 手，長さ，幅 » 足，長さ，幅，ではなく手長，手幅，足長，足幅 • Rで処理するなら見出しは英語の方がよい（特に Macintoshユーザー） • 半角英数と全角英数を混ぜない，半角カナは使わない

MS‐Excelの功罪

• Excelの良い点

– データ表を見たまま，思うままに加工できる – オートフィルターはデータチェックに便利 – ピボットテーブルも役に立つ

• Excelの良くない点

– データの入力，点検，修正，作図，解析を繰り返すと， 色々な状態のファイルが沢山できてしまう – セルに入力されているのが値か式か確認しにくい（他人 が何をやっているか読み取れない） – どの状態のファイル（どのデータ）でどの解析をやったか， 後から確認できない（自分が過去にやったことを思い出せない）

MS‐Excel

（もしくはMS‐Access） ・データ入力 ・エラーチェック ・データ修正

私のデータ取り扱い術

Excel形式の原データは１ 個のファイルに固定する （微妙に内容が違う複数の ファイルを作らない） 修正履歴を記録する データ修正

R

・データ読み込み ・データ集計，部分抽出，加工 ・作図 ・データ解析 ・一連のデータ処理をスクリプト （プログラム）として保存すること で作業の透明性，再現性を確保

について

• 解析作業毎にフォルダを作り，そこにデータ，Rス

クリプト，作業履歴などを保存すると良い

• 最初はRを起動して，ディレクトリの変更を行い，

新しいスクリプトを作製して保存する。

• 次からはR.Dataをダブルクリックすればそこから

スタートできる（win）。Macならｒのスクリプトをダ

ブルクリックすれば良い

• コンソール，エディタ，グラフィックの３ウィンドウ

• エディタからctrl + R（Win）, command + リターン

（Mac）で実行できる。

Rの基礎：四則演算

> 1 + 1

足し算

> 1 – 1

ひき算

> 2 * 2

かけ算

> 4 / 2

わり算

> 3^2

２乗，べき乗

> sqrt(9)

平方根

> 9^(1/2) ½乗（平方根）

お役立ちサイト

• R‐Tips

http://cse.naro.affrc.go.jp/takezawa/r‐tips/r.html

• Rjp Wiki

http://www.okada.jp.org/RWiki/

• そのほか多数のウェブサイトがあります

• わからないことは，キーワードや関数名を

使って，どんどんググってみよう！

オブジェクト

• オブジェクトとは

数値，数列（ベクトル），行列，データフレーム（表形式

のデータ）などを名前で操作できる

> a <‐ 3.14 > b <‐ c( 1, 2, 3, 4, 5, 6) > x <‐ matrix(nrow = 2, ncol = 3) > data <‐ read.csv(“HandSize.csv”, header = T) 

• 要素の指定

b[3]

x[2, 3]

data$HL

• オブジェクトと数値の演算，オブジェクトどうしの

演算をうまく使えば，効率良く計算処理できる

データの読み書き

• CSV（コンマ区切りテキスト）ファイルを使うの

が便利

– エクセルからの書き出し

• 『名前をつけて保存』CSV（カンマ区切り）（*.csv）

– Rでの読み込み

• read.csv(“ファイル名.csv”, header = T)

– Rでの書き出し

• write.csv(オブジェクト名，“ファイル名.csv”）

部分データの取り出し

• [,]を使った行操作

> d[d$HL==“F”,]

# dの中からd$HL==“F”がTRUEである行を取り出した部分データ

> d[d$HL==“F”,]$HL

# 同上の部分データの中のHL値

• subest()を使った操作

> d.f <‐ subset(d, HL == “F”)

# dのうちHL==“H”である部分データをd.fに代入

> d.f$HL

作図…plot(), boxplot(), histogram()

• plot(HL, data = d)

– データ番号順にHLをプロットする

• boxplot(HL, data = d)

– 箱ヒゲ図を描く

• plot(HL ~ Sex, data = d)

– 性別にHLのボックスプロットを描く（Sexがカテゴリカル

変数なので）

• plot(HL ~ BL, data = d)

– 散布図（両方とも連続変数なので）

• histogram(d$HL, breaks = 数列)

– breaksに数列を指定すると階級幅を

データの代表値aと残差平方和

• データX

₁

, X

₂

, ...X

_n

を１つの数値 a で代表させるこ

とを考える

• データと a との食い違い（残差）

– X

₁

– a, X

₂

– a, ...X

_n

– a

• n個のデータ全体としてどれだけ食い違っている

か（残差の総量）？

– Σ（X

_i

– a)

# 残差にプラスとマイナスがあると打ち

消し合うので良くない

– Σ|X

_i

– a|

# 絶対値は計算が面倒

– Σ（X

_i

– a)

2

_{# 2乗和だと何かと便利（これを残差平}

方和と呼ぶ）

平均は残差平方和を最小にする

データの代表値

• 代表値 a に対する残差平方和を a の関数f(a)と考

える（aの値を変えると残差平方和は増減する）

f(a)  = (X

₁

– a)

2

+ (X

2

– a)

2

+ ...(X

n

– a)

2

= X

₁2

+ X

22

+ ...X

n2

– 2a(X

1

+ X

2

+ ...X

n

) + na

2

• 残差を最小にする a は，集団の良い代表値であろ

う

• f(a)の値を最小にする a は，f(a)をaで微分した導

関数f’(a) = 0として求められる

–2(X

₁

+ X

₂

+ ...X

_n

) + 2a = 0

a = (X

₁

+ X

₂

+ ...X

_n

) /n

（平方完成でも同じように確認できる）

正規分布 N（μ, σ

2

）

• データの確率分布

データのバラツキ具合を確率的

に表現するもの

• 正規分布

自然界で観察される連続変数

（体長など）がしばしば従うシンプ

ルな分布

平均μのまわりに左右対称の釣

鐘型にデータがバラつく（分散σ

2

）

素直で性質の良い分布なので，

データが正規分布に従うと考える

と色々計算しやすい

統計モデルって？

• データの変動を数学的に説明するもの

• 例）

手長は雌雄によって違う（分散分析）

手長は体長によって違う（直線回帰）

• 線形モデル lm(応答変数∼説明変数)

残差が分散の等しい正規分布に従うと仮定する

モデル

１つの関数lm()を使って分散分析，直線回帰，共

分散分析などができる

m0）全平均モデル（nullモデル）

全平均を計算するだけでもモデル （Null model） 全平均μが集団を代表する値である μの推定値の回りの変動は ランダム誤差である （正規分布に従う） （データ数 − パラメータ数） …自由度 …残差平方和 残差分散（残差平均平方）：



 2 ) (x  …残差標準誤差

√

残差分散

m0）全平均モデル

P値って何？

lm(formula = HL ~ 1, data = d) Residuals:

Min 1Q Median 3Q Max -1.64762 -1.02262 0.00238 0.65238 2.85238 Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 17.6476 0.1683 104.8 <2e-16 *** ---全平均 残差標準偏差 √残差平方 P値 ｔ値 帰無仮説の下で統計検定量（t値）がここまで 大きな値を取る確率 帰無仮説って何？ _{ここでは全平均 = 0} じゃあ t値って何？

統計検定の手順

① 帰無仮説と対立仮説を設定する  帰無仮説 確認したいこととは逆の仮説（HLの平均はゼロ）  対立仮説 確認したい仮説（HLの平均はゼロじゃない） ② 帰無仮説が正しいと仮定して，統計検定量（この場合t値） を計算する ③ 帰無仮説が正しい時，統計検定量がそこまで外れた値を取 る確率（p値）を求める ④ 有意水準αを設定する（通常α = 0.05 = 1/20） ⑤ p値が有意水準より小さければ，帰無仮説を捨てて（棄却 し），対立仮説を採用する ⑥ 『有意水準α でHLの平均はゼロと有意に異なる』といった結 論が得られる

有意水準 α

• 帰無仮説を棄却するp値の基準値

• たとえばp = 0.05で帰無仮説を棄却した場合

– 帰無仮説が正しいとき，その現象が起る確率は1/20以下 – 帰無仮説の下でそんな希な現象が起ったと考えるより， そもそも帰無仮説が間違っていたと判断する方が妥当で あろう...と判断する

• 有意水準を大きくすると，（本当は差がないのにある

と）ウソを言う確率が高くなる（Type Iエラー）

• 有意水準を小さくし過ぎると，違いがあっても検出し

にくくなる（検出力の低下）

m0）全平均モデル

P値って何？

lm(formula = HL ~ 1, data = d) Residuals:

Min 1Q Median 3Q Max -1.64762 -1.02262 0.00238 0.65238 2.85238 Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 17.6476 0.1683 104.8 <2e-16 *** ---残差標準偏差 √残差平方 全平均 ｔ値 P値 帰無仮説の下で統計検定量（t値）がここまで 大きな値を取る確率 帰無仮説って何？ _{ここでは全平均 = 0} じゃあ t値って何？

ｔ 値とは？

(平均値−予測値）÷ （xの分散/データ数）

> t.cal <‐ (mean(d$HL) – 0) / sqrt(var(d$HL)/n)   自由度＝n‐1のt分布に従う（nはデータ数） 25 d$HLの平均は0と違うかどうか？ t=104.8, df=41, p-value=0 平均値の違いをデータのバラツキ（分散）とデータ数を考慮した 上で評価する尺度（検定統計量） 自由度41のt分布

m1）ANOVA（分散分析）モデル

平均値が雌雄によって異なると説明するモデル lm(HL ~ Sex, data = d)

全平均値だけで説明するモデル lm(HL ~ 1, data = d)

m1）ANOVA（分散分析）モデル

▶ メスの平均値は０と有意に異なる

▶ 雌雄の平均差がゼロ（帰無仮説）と考えたとき，_{t値がこんなに} 大きくなる確率は非常に小さい＝雌雄の平均は有意に異なる

lm(formula = HL ~ Sex, data = d) Residuals:

Min 1Q Median 3Q Max -1.5944 -0.4944 -0.1160 0.6375 1.9056 Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 16.9375 0.1463 115.778 < 2e-16 *** SexM 1.6569 0.2235 7.415 4.99e-09 *** ---メスの平均がゼ ロと異なるかど うかのp値 メス（アルファ ベットで早い方） の平均値 オス平均値−メス平均値 雌雄の平均値の差が ゼロと異なるかどうかのp値

...ということで

• メスとオスの手長の平均値は異なることが確

認されました

• lm(HL ~ Sex, data = d)でやっているのは２標本

のt検定と同じことです

> var.test(BL.f, BL.m)

• 分散が等しいことを確認

> t.test(BL.f, BL.m, ,var.equal = TRUE)

• BL.fとBL.mの平均値が等しいか検定 • t値，p値はlm()と全く同じになります

しかし...

• オスとメスでは体の大きさが違うはず

m2）直線回帰モデル

▶ 手長は体長と直線関係にある ▶ _{Y切片はゼロと有意に異なるとは言えない} ▶ 傾きはゼロと有意に異なる lm(formula = HL ~ BL, data = d) Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) -1.33514 1.83554 -0.727 0.471 BL 0.11501 0.01111 10.354 7.02e-13 *** ---Y切片 Y切片がゼロと異なる かどうかのp値 傾き 傾きがゼロであ るかどうかのp値

m3）ANCOVA（共分散分析）モデル

▶ メスの切片はゼロと有意に異ならない ▶ 傾きはゼロと有意に異なる

▶ 雌雄の切片は有意に異ならない

lm(formula = HL ~ BL + Sex, data = d) Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 2.72263 2.71358 1.003 0.3219 BL 0.08902 0.01698 5.243 5.8e-06 *** SexM 0.54183 0.27437 1.975 0.0554 . ---切片の雌雄差 _{切片の雌雄差が} ゼロであるかどう かのp値 0.05に近く微妙

m4）ANCOVAモデル（交互作用あり）

▶ 切片はゼロと有意に異ならない ▶ メスの傾きはゼロと有意に異なる

▶ オスの切片はメスの切片と有意に異ならない ▶ オスの傾きはメスの傾きと有意に異ならない

lm(formula = HL ~ BL + Sex + BL:Sex, data = d) Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 4.72776 3.44481 1.372 0.17798 BL 0.07646 0.02156 3.546 0.00106 ** SexM -5.01788 5.87759 -0.854 0.39860 BL:SexM 0.03320 0.03506 0.947 0.34965

---モデル式と直線関係式

モデル式 関係式（メス） 関係式（オス） HL ~ 1 ̂ HL ~ Sex ̂ ̂ +  HL ~ BL ̂ +  BL HL ~ BL + Sex ̂ +  BL ̂ +  ) +  BL HL ~ BL + Sex    + BL:Sex ̂ +  BL ̂   +  ) +  BL

a

Intercept

b BL

a

_m

SexM

b

_m

BL:SexM

summary係数表出力

以上の結果から...

１）手長だけを雌雄比較すると

▶

メス ＜ オス （有意差あり）

２）体長と手長の関係を考慮すると

▶

手長∼体長 （手長は体長に比例する）

３）上の直線関係を考慮した上で

▶

同一体長で雌雄比較すると雌雄に有意差はない

▶

_{ただし，p値は微妙なので，もう少しデータを増や}

すことが望ましい

150 160 170 180 16 17 18 19 20 BL HL male female Model 0: HL ~ 1 150 160 170 180 16 17 18 19 20 BL HL Model 1: HL ~ Sex 150 160 170 180 16 17 18 19 20 BL HL Model 2: HL ~BL 150 160 170 180 16 17 18 19 20 BL HL Model 3: HL ~ BL + Sex 150 160 170 180 16 17 18 19 20 BL HL

Model 4: HW ~ BL + Sex + BL:Sex

今までに検討したモデル

統計検定は，単純な モデルを帰無仮説とし て，複雑なモデル（対 立仮説）の妥当性を 逆説的に検証する

赤池情報量規準AIC

☞ 『真に正しいモデル』からの距離を反映する尺度 ☞ AIC値が小さい方がデータの変動を合理的に説明する 良いモデル（差が2以内だと微妙） ☞ 絶対値に意味はなく相対評価 （包含関係にないモデルも比較できる）

最大対数尤度

パラメータ数

モデルのデータへのあては まり 過剰説明に対するペナ ルティー データ数が違うと比較できないので注意！（e.g. 欠損データ）

AICによるモデル選択

AICでは Model 3がベスト Model 4も捨てがたい でも係数では雌雄差 が有意じゃないという 微妙な結果に

統計検定とAICによる

モデル選択

BL*Sex BL +Sex BL Sex 1 AIC

BL*Sex − 0.3496 0.105 2.544e‐05 8.194e‐12 75.8 BL+Sex − − 0.0554  5.802e‐06 1.441e‐12  74.8 BL − − − N/A 7.017e‐13 76.8 Sex − − − − 4.994e‐09  95.2 1 − − − − − 129.5 共分散分析 交互作用あり 共分散分析 直線回帰. 分散分析 全平均 ☞ この例では検定とモデル選択の結果が必ずしも一致せず， スッキリしない結果になりました（そういうこともあります） ☞ スッキリさせるには，データを増やす，測定誤差を小さくする

まとめ

• ExcelとR

☞Excelはデータの入力チェックなどに便利

☞Rに馴れるて来るとExcelの良くない面が見えてくる

• 統計モデルとは

☞データの変動を単純化して数学的に説明するもの

☞推定値のまわりの残差平方和が小さくなるようにパ

ラメータを推定する

• 統計検定とは

☞帰無仮説と対立仮説を用意し，ある有意水準で帰無

仮説を捨てられるかどうか検討する（実験系向き）

• モデル選択とは

☞複数のモデルを比較し，データの変動を合理的に説

明するものを選び出す（フィールド科学向き）

参考書

• 舟尾暢男：データ解析環境「R」―定番フ リーソフトの基本操作からグラフィックス、 統計解析まで ，工学社 • 久保拓也：データ解析のための統計モ デリング入門――一般化線形モデル・ 階層ベイズモデル・MCMC (確率と情報 の科学)，岩波書店 • Crawley: 統計学:Rを用いた入門書，共 立出版 • グラフェン：一般線形モデルによる生物 科学のための現代統計学―あなたの 実験をどのように解析するか，共立出 版

尤度比検定の制約事項

• 比較するモデルは包含関係（

ネスト

）になっ

ていなければならない

• 何回も検定を繰り返すと

多重比較

の問題を

生じる

• 帰無仮説が棄却されなかった場合，帰無仮

説が正しいとは判断できない

• カイ2乗近似が良くないことがあるので，

ブートストラップ

を使う方が望ましい

尤度比検定 anova(h0, h1)

• 2つのモデルh

₀

とh

₁

の最大対数尤度の差を2倍した

ものは，自由度＝パラメータ数差のカイ2乗分布に

従う

• 正規誤差の下では，F検定と同値

• カイ2乗近似は必ずしもよろしくないので，ブートスト

ラップ（パラメトリック，ノンパラメトリック）を使う方が

望ましいようだ

p.value = 1 – pchisq(2*(MLL.h1$value – MLL.h0$value), h1のパラメータ数 – h0のパラメータ数)