090614 R 講習会 in 名大 生命農学研究科: R で GLM をやってみよう

森林生態生理学研究分野 玉木 一郎

で GLM をやってみよう〜glm 関数の利用と手計算での最尤推定〜
玉木 一郎∗

2009/06/14
概要： 本テキストでは，まず一般化線形モデル（Generalized Linear Model: GLM）の概要について説明を行う。続いて，応答変数が二項分布とガンマ分布に従う場合の 2 つの例題を用いて，実際に glm 関数を用いたパラメータ推定，予測，信頼区間の推定を行う。さらに，同例題のデータを用いて直接尤度関数を記述し，optim 関数を用いた最尤法によるパラメータ推定を行うことで，GLM と最尤推定への理解を深めることを目的とした。

研究で用いられるようになりました。普通の林学の講義

目次で取り扱われることはまず無いかと思いますが，論文や学会発表などで名前を聞いたことのあるという人も多い1

はじめに: GLM とは？

122

2.1
2.2
2.3
2.42.5ドブソン（2008）2) P. 144 の二項分布の例
データを眺める . . . . . . . . . . . . . . . . . . . . . .
統計モデリングと glm 関数によるパラメータ推定 . . .
予測値とその 95% 信頼区間の推定 . . . . . . . . . . .
パラメトリック・ブートストラップ . . . . . . . . . . .
optim 関数による最尤推定 . . . . . . . . . . . . . . . .

43.1
3.2
3.3
3.4
3.5

McCullagh & Nelder（1989）5) P. 300 のガンマ分布の例
データを眺める . . . . . . . . . . . . . . . . . . . . . .
統計モデリングと glm 関数によるパラメータ推定 . . .
予測値とその 95% 信頼区間の推定 . . . . . . . . . . .
パラメトリック・ブートストラップ . . . . . . . . . . .
optim 関数による最尤推定 . . . . . . . . . . . . . . . .

44455測子，応答変数の期待値と線形予測子を結びつけている

5関数 g を連結関数と呼びます。連結関数は恒等関数や対4.1
4.2
4.3終わりに他に勉強すると良いこと . . . . . . . . . . . . . . . . .
参考書・サイト . . . . . . . . . . . . . . . . . . . . . .
R で統計をする際に . . . . . . . . . . . . . . . . . . .

34のではないかと思います。応答変数を y，説明変数を xとすると，一般化線形モデルは以下の式で示すことがで22233きます*3 。

g[E(yi )] = ηi = β0 + β1 xi1 + · · · + β j xi j + · · · + βm xim
yi ∼ f (yi | θ)
E は期待値を，β は係数を，添字 i はデータの番号を，
j はパラメータの番号を示します。ここで，η を線形予

666数関数，ロジット関数など様々な形をとることができます。どれを用いるかはモデリングの仮定やモデルのあてはまりを吟味した上で決めれば良いかと思います。特に

1 はじめに: GLM とは？

問題がなければ，応答変数が従うと仮定した確率分布にまずはじめに，このテキストでは，ふるまいに注目し対応する正準連結関数を使用するのが良いでしょう。∼ている変数のことを応答変数，これらを説明するために（チルダ）は左辺が右辺の確率分布に従うことを意味し用いる変数を説明変数と呼び*1 ，応答変数は何らかの確ます。 f は y が従う確率分布の確率関数（離散分布の場率分布に従うランダムな変数，説明変数はばらつきのな合）もしくは確率密度関数（連続分布の場合）を，θ はい固定された値として扱います。

そのパラメータを示します。一般化線形モデルでは，応GLM は Generalized Linear Model の頭文字を並べた

答変数はその型に応じて，正規分布や二項分布，ポアソもので，日本語では一般化線形モデルと呼ばれるもの

ン分布などの指数型分布族に含まれる確率分布に従うとです。普通の線形モデル（重回帰や分散分析）やロジス仮定します*4 。

ティック回帰，対数線形回帰などの応答変数の従う確率

一般化線形モデルでは最尤法でパラメータ推定を行い分布が異なる線形モデルを同一の枠組みの中で取り扱おます。先ほどの確率（密度）関数 f をデータが与えられうというものです*2 。近年 R が普及するにつれ，多くのた時のパラメータ θ の関数（L(θ | yi )）としてみると，n∗ 名古屋大学大学院 生命農学研究科 森林生態生理学研究分野・日帰と分散分析までを総合的にとらえたものです。
*3 いわゆる誤差項は無いことに注意。

本学術振興会特別研究員*1 応答変数は他にも従属変数とか目的変数，説明変数は独立変数とか予測変数と呼ばれています。
*2 名前の良く似た一般線形モデル（General linear model）は重回*4 英語文献中で出てくる，いわゆる error distribution は誤差分布と訳されることもありますが，誤解をまねくので，このレジュメ中では応答変数の従う確率分布としました。

1

090614 R 講習会 in 名大 生命農学研究科: R で GLM をやってみよう

森林生態生理学研究分野 玉木 一郎

個のデータを観測した時の尤度関数 L は以下の式で示す

1.0●

ことができます。

●L(θ | y1 , · · · , yi , · · · , yn ) =n∏Mortality rate

0.8L(θ | yi )i

この尤度関数を最大にするようなパラメータ θ を最尤推

0.6
●0.4
●
●0.2

定量と呼びます。実際には L の対数をとったもの（対数●0.0尤度）を扱う方が便利な場合が多いです。

1.70

GLM を R で実行する上では，連結関数と応答変数の1.75

1.801.85

Gas concentration [log10CS2mgl−1]従う確率分布について，ある程度の知識があれば十分か
図1

ガス濃度とカブトムシの死亡率（死亡個体数 / サンプ
ル個体数）の関係。丸のサイズはサンプル個体数に比例する。

もしれませんが，尤度関数の立て方のあたりまで知っておくと理解が深まって，応用も利くと思います。

2 ドブソン（2008）2) P. 144 の二項分布の例
2.1

●

●

統計モデリングと glm 関数によるパラメータ推定2.2

ここで，データを眺める

この節では二硫化炭素ガスに 5 時間暴露されたカブト• カブトムシの死亡個体数（y）が死亡率（ p = y/N ）ムシの死亡個体数とガス濃度との関係を GLM で解析しとサンプル個体数（ N ）を上限値に持つ二項分布にてみます（表 1）。

従う。

• 死亡率（ p）はガス濃度（ x）が増加すると共に高ま

表 1 カブトムシの死亡データ

ガス濃度

実験に供したるが，0 から 1 の間の値をとる。
死亡数• 連結関数はロジット関数とする。

（log10 CS2 mgl−1 ） カブトムシの数と仮定すると，以下の式で表現することができます。

1.69075961.7242

60

131.7552

62

181.7842

56

28

1.81136352

11 + exp(−ηi )

1.83695953yi ∼ Binomial(pi , Ni )1.861062611.8839

6060

log(pi) = ηi = β0 + β1 xi1 − pipi =これを R の glm 関数では以下のように指定します。

family = binomial のデフォルトの連結関数はロジット関数になっています。詳しい指定方法については?family を参照して下さい。
¶

どうやらガス濃度が高くなるとカブトムシが死にやすくなりそうです。ただし，死亡数には上限があります。

> beetle <- cbind(dead.beetle, n.beetle - dead.beetle)> glm.beetle <- glm(beetle ˜ gas, family = binomial)> glm.beetle$coefficients(Intercept)
gas-60.71745
34.27033

ガス濃度を横軸に，死亡率（死亡個体数 / サンプル個体数）を縦軸にとり，図示してみましょう（図 1）。このような割り算データは分母にくる数のサイズが分かるようµに図示するのが誠実です（この例では処理間でサンプル2.3個体数がほとんどばらついていないのであんまり意味ないのですが...）。

¶¶> gas <- c(1.6907, 1.7242, 1.7552, 1.7842,+
1.8113, 1.8369, 1.8610, 1.8839)> n.beetle <- c(59, 60, 62, 56, 63, 59, 62, 60)> dead.beetle <- c(6, 13, 18, 28, 52, 53, 61, 60)> plot(gas, dead.beetle/n.beetle, cex = n.beetle/40)µ´

予測値とその 95% 信頼区間の推定

次にモデルによる予測値を図に重ねます（図 2）。

³

³

³

> x <- seq(min(gas), max(gas), by = 0.01)> eta.pred <- glm.beetle$coefficients["(Intercept)"] ++glm.beetle$coefficients["gas"]*x> p.pred <- 1/(1 + exp(-eta.pred))> lines(x, p.pred)µ´´続けて予測値の 95% 信頼区間（パラメータの 95% 信頼区間については後述）を計算し，図に重ねます。このと2

090614 R 講習会 in 名大 生命農学研究科: R で GLM をやってみよう

1.0●

Mortality rate

optim 関数による最尤推定

2.5●

n 個のデータを観測した時の尤度は L =

●

●

0.8

森林生態生理学研究分野 玉木 一郎∏nif (yi ) という掛け算の形をとります。このままでは f (y) が極端0.6に小さかったり大きかったりする場合にふれが大きくな●るため，対数尤度 LL = log(L) =0.4i log( f (yi )) について考えます。LL は単純に L の対数をとっただけなので L●●

0.2∑nが最大のとき，LL も最大値をとります。この LL を最●大にするようなパラメータを求めたいわけです。実際に

0.0
1.70

1.75

1.80

optim 関数は，与えられた関数を最小にするようなパラ

1.85

Gas concentration [log10CS2mgl−1]

メータを推定するので， NLL = −LL について考えます。
二項分布の確率関数は f (y) =N Cy py (1 − p)N−y なので，

図2

ガス濃度とカブトムシの死亡率（死亡個体数/サンプル個体数）の関係。丸のサイズはサンプル個体数に比例する。
実線は予測値を，破線は予測値の 95% 信頼区間を示す。

optim に与える関数を以下のように記述します。
¶> f.beetle1 <- function(parameters){+beta0 <- parameters[1]+beta1 <- parameters[2]+eta <- beta0 + beta1*gas+p <- 1/(1 + exp(-eta))+y <- dead.beetle+
N <- n.beetle+
NLL <- -sum(log(choose(N, y)*pˆy*(1 - p)ˆ(N - y)))+ }き，実際にはデータポイントごとに二項分布のサンプル数（N ）が異なりますが，今回はそれらの調和平均を丸めた値を用いて，N を固定して信頼区間を推定しました。
図を見ると信頼区間の幅は死亡率が 0 や 1 に近い時に狭µく，0.5 付近で最も広くなっています（物差をあてて見続いて適当な初期値を与えて optim を実行します。初ると分かりやすいです）。

¶´

期値は関数型を考えたり，グラフから適当に読み取るべ³

きなのですが，今回はズルして glm の推定結果に近い値> mean.N <- round(length(n.beetle)/sum(1/n.beetle))> p.lower <- qbinom(0.025, mean.N, p.pred)/mean.N> p.upper <- qbinom(0.975, mean.N, p.pred)/mean.N> lines(x, p.lower, lty = 2)
> lines(x, p.upper, lty = 2)µ2.4

パラメトリック・ブートストラップを入れておきました*5 。

¶

最初に行った glm.beetle$coefficients の値と見比たので，今度はパラメトリック・ブートストラップを用べて見てください。ほぼ等しいことが分かります。今回いた回帰係数の信頼区間の推定方法を紹介します。このは確率関数を直接記述しましたが，R には各種確率分布方法を簡単に説明すると，まず推定したパラメータを含の確率（密度）関数が最初から入っているので，信頼区むモデルを用いて乱数データを発生させ，その乱数デー間の計算をしたときと同様に，それらを利用することがタを応答変数にしてパラメータを再推定します。このシできます。

¶

ミュレーションを繰り返すことで，パラメータの信頼区間を推定します。

¶> n.sim <- 1000 # シミュレーション回数> beta <- data.frame(beta0 = rep(NA, n.sim),+
beta1 = rep(NA, n.sim))> eta <- glm.beetle$coefficients["(Intercept)"] ++glm.beetle$coefficients["gas"]*gas> p <- 1/(1 + exp(-eta))> N <- n.beetle> sapply(1:n.sim,+function(i){+y.random <- rbinom(length(p), N, p)+
beetle.random <- cbind(y.random, N - y.random)+
beta[i, ] <<- glm(beetle.random ˜ gas,+family = binomial)$coefficients+})> beta.ci <- apply(beta, 2, quantile,+
prob = c(0.025, 0.975))> beta.ci
beta0
beta1
2.5% -72.61395 29.59513
97.5% -52.42950 40.95304³> optim.beetle1 <- optim(c(-60, 30), fn = f.beetle1)´ >[1]optim.beetle1$par-60.73001 34.27737µ

先ほど，予測値の 95% 信頼区間の計算方法を紹介しµ³³´*5 できるかぎりそれっぽい値になるよう努力して下さい。初期値

3³> f.beetle2 <- function(parameters){+beta0 <- parameters[1]+beta1 <- parameters[2]+eta <- beta0 + beta1*gas+p <- 1/(1 + exp(-eta))+y <- dead.beetle+N <- n.beetle+
NLL <- -sum(dbinom(y, N, p, log = T))+ }> optim.beetle2 <- optim(c(-60, 30), fn = f.beetle2)µ´´

によってうまくいかない場合があります。

090614 R 講習会 in 名大 生命農学研究科: R で GLM をやってみよう

統計モデリングと glm 関数によるパラメータ推定

3.2

3 McCullagh & Nelder（1989）5) P. 300 のガ
ンマ分布の例

ここで，• 凝固時間（y）は正の値を示す連続データなので，ガデータを眺める

3.1

森林生態生理学研究分野 玉木 一郎

ンマ分布に従う。

この節では血液の凝固時間（秒）と血漿濃度（%）の関• 血漿濃度の対数値（ x）に対する凝固時間の減少具合

係を GLM で解析してみます（表 2）。このデータは glmが双曲線的なので，連結関数には逆数関数を用いる。

関数のヘルプに入っているものですが，今回はそのうち• shape パラメータ（α）は常に一定。の lot1 だけ使います（そのため変数名を time に変更• scale パラメータ（ s）が線形予測子の変化に応じてしました）。

変化する*7 。

• ガンマ分布の平均は αs，分散は αs2 。

表 2 血液の凝固時間（秒）と血漿濃度（%）の関係

血漿濃度（%）

凝固時間（秒）と仮定すると，以下の式で表現することができます。

511810581542E(yi ) = αsi = 1/ηi2035yi ∼ Gamma(a, si )30274025

これを R の glm 関数では以下のように指定します。

6021family = Gamma のデフォルトの連結関数は逆数関数に

8019¶

10018ηi = β0 + β1 xiなっています。
> glm.time <- glm(time ˜ log(u), family = Gamma)> glm.time$coefficients(Intercept)
log(u)-0.01655438 0.01534311µ

血漿濃度が高くなるにつれて凝固時間は短くなるよ³´log(u) の回帰係数が正の値になっていますが，これはうです。横軸に血漿濃度の対数値*6 ，縦軸に凝固時間を連結関数が逆数関数のためです。この場合のように，回とって図示してみます（図 3）。

¶³帰係数だけ見て「説明変数は応答変数に対して正の影響> u <- c(5, 10, 15, 20, 30, 40, 60, 80, 100)
> time <- c(118, 58, 42, 35, 27, 25, 21, 19, 18)> plot(log(u), time, ylim = c(20, 140))µを及ぼしている」と判断しないで下さい。

3.3´

予測値とその 95% 信頼区間の推定

次にモデルによる予測値とその 95% 信頼区間を計算しますが，信頼区間を計算する際に，shape パラメータの最尤推定値が必要になります。そこで，まず shape パ

140120

ラメータを計算しておきます。ガンマ分布の GLM では●

shape パラメータは dispersion パラメータの逆数になTime [s]100

り ま す 。1/summary(glm.time)$dispersion で 計 算80できるのですが，glm 関数では繰り返し重み付け最小60●

自乗法でパラメータ推定を行っているので，この shape●40

パラメータの推定値は最尤推定値と異なる場合がありま

●
●

20
1.5

2.02.53.0●

3.5●

● ●

す*8 。というわけで，以下のようにして shape パラメー

4.04.5タ（a）の最尤推定値を計算します。

log(Plasma concentration) [log(%)]

図3*7 scale パラメータの代わりに rate パラメータ（1/s）を指定する

血漿濃度（%）の対数値と血液の凝固時間の関係こともできます。ガンマ分布を扱う際には，混乱を避けるために scale なのか rate なのか，どちらの話をしているのか明示しておく必要があります。
*8 といってもどちらの shape パラメータでも信頼区間の幅はそん
なに変らないのですが...。ちなみに回帰係数の推定値は最尤推定値とほとんど変わりません。

*6 いきなり対数をとっていますが，とらない場合も後で試してみて下さい。尤度がずいぶん低くなります。

4090614 R 講習会 in 名大 生命農学研究科: R で GLM をやってみよう¶

森林生態生理学研究分野 玉木 一郎³µ

同様に optim 関数で最尤推定します*9 。

¶´

この値を用いて予測値とその 95% 信頼区間を計算し，図示します（図 4）。

¶³> x <- seq(1.6, 4.6, by = 0.1)
> eta.pred <- glm.time$coefficients["(Intercept)"] ++glm.time$coefficients["log(u)"]*x> time.pred <- 1/eta.pred> s.pred <- 1/eta.pred/a # scale パラメータ> time.lower <- qgamma(0.025, shape = a, scale = s.pred)> time.upper <- qgamma(0.975, shape = a, scale = s.pred)> lines(x, time.pred)> lines(x, time.lower, lty = 2)
> lines(x, time.upper, lty = 2)µ³

> f.time <- function(parameters){+beta0 <- parameters[1]+beta1 <- parameters[2]+a <- parameters[3] # shape パラメータ+eta <- beta0 + beta1*log(u)+s <- 1/eta/a # scale パラメータ+
NLL <- -sum(dgamma(time, shape = a,+scale = s, log = T))+ }> optim.time <- optim(c(-0.016, 0.015, 530), fn = f.time)> optim.time$par # beta0, beta1, shape パラメータ[1] -0.01655432
0.01534309 538.14832592µ´

4 終わりに´

今回は R の統計機能の紹介を担当してくれとのこと

で，R と言えば GLM かなと思い，自分が GLM の存在これまた少し分かりにくいので，物差をあてて見てほしを知った頃に知りたかったこと，調べるのが困難だったことを中心にまとめました*10 。それでもけっこうな量に140なってしまいました...。今回述べてはいないけど知って120●

おいて欲しいことは他にもあります*11 。キーワードを後

100

Time [s]optim 関数による最尤推定3.5> 1/summary(glm.time)$dispersion[1] 408.8208> library("MASS")> a <- gamma.shape(glm.time)$alpha # shape パラメータ> a[1] 538.1315

で述べておくので，自分で調べてみて下さい。

80

最近では GLM に関する書籍や情報も増えてきました

60が，これ一つだけを見ておけば大丈夫というものは特に●無いように思われます。自分がこれまで触れてきた中で

●

40●
●

20
1.5

2.0

2.5

3.0

3.5

●

●

● ●

4.04.5

良いと思うものをあげておくので参照してみて下さい。
いろいろな本を読んで，R で実験して自身の知識を高めて欲しいと思います。

log(Plasma concentration) [log(%)]

また，情報源によって表現が異なったり*12 ，勘違いし図 4 血漿濃度（%）の対数値と血液の凝固時間の関係。実線はモデルの予測値を，破線は予測値の 95% 信頼区間を示す。ているものも時にはあるので注意が必要です。本や論文中の数式も間違っていることも時にはあります。自分いのですが，x 軸の値が小さくなると（scale パラメータも極力それを無くすように記述してきましたが，文章やが大きくなる）と信頼区間の幅が大きくなるのが分かり

コードの誤り，不適切な表現等ありましたらご連絡下さます。

い*13 。

3.4

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

先ほどの二項分布の例と同様にパラメータの信頼区間を推定します。

¶> n.sim <- 1000
> beta <- data.frame(beta0 = rep(NA, n.sim),+
beta1 = rep(NA, n.sim))> eta <- glm.time$coefficients["(Intercept)"] ++glm.time$coefficients["log(u)"]*log(u)> s <- 1/eta/a # scale> sapply(1:n.sim,+function(i){+random.y <- rgamma(length(s),+shape = a, scale = s)+beta[i, ] <<- glm(random.y ˜ log(u),+family = Gamma)$coefficients+})> beta.ci <- apply(beta, 2, quantile,+
prob = c(0.025, 0.975))> beta.ci
beta0
beta1
2.5% -0.01809224 0.0146059397.5% -0.01483667 0.01601967µ³´5*9 警告が出ますが，ガンマ分布の shape パラメータや scale パラメータは正の値しかとらないので，パラメータ値の探索の過程でそれらが負になると警告が出ます。しかし，optim は気を取り直して計算を続けてくれているので気にしなくて良いです。method = "L-BFGS-B"にすると，パラメータの探索範囲を狭めることができますが，今回の例題ではデフォルトの"Nelder-Mead"の方がうまくパラメータを推定できるようです。気になる人は試してみて下さい。
*10 GLM を使ってみたくて使い始めた R ですが，実は自分の研究では glm 関数を使う機会がありませんでした...。
*11 多重検定とか多重共線性，検定 vs モデル選択，頻度主義者 vsベイジアンみたいに，面倒いし宗教的だったりして，自分ではあんまり触れたくないこともあります...。
*12 単に deviance と書いてあっても，−2LL のことを指さない場合もあります。
*13 連絡先（garageit(at)gmail.com）。あ，2 章と 3 章で変数名がかぶっているのがあるのは許して下さい...。

090614 R 講習会 in 名大 生命農学研究科: R で GLM をやってみよう

4.1

森林生態生理学研究分野 玉木 一郎

他に勉強すると良いことません。統計学授業は一読しておくと良いと思

適当につらつらキーワードを並べておきます。

います。ここで offset 技も習得して下さい。

– RjpWiki6) : 一般的な関数について分からないこ• いろいろな確率分布の特徴（連続 or 離散，上限・下とは，まずはここで検索しましょう。

限，平均と分散の関係，関数型等々）。

• 回帰係数の検定。

4.3

R で統計をする際に• 情報量基準とモデル選択。

• なるべく新しいヴァージョンの R を使いましょう• 多重共線性。

（バグが修正されているので）。昔インストールした• 引数 offset の使い道。

古いのを使い続けないこと。

• 過分散（overdispersion）への対処。

• 検定とかする前に，必ず図を描いて，視覚的に確か• 一般化線形混合モデル（Generalized Linear Mixedめること（河崎君担当箇所参照）。

• 関数の使い方が分からない時，初めて使う時は必ずModel: GLMM）やランダム効果。
• ベイズ統計学。

関数の help を参照しましょう。

• マルコフ連鎖モンテカルロ（Marckov Chain Monte• まともなテキストエディタを使いましょう（安部君Carlo: MCMC）法。

担当箇所参照）。

• 統計学の歴史的経緯。
• ノンパラメトリック・ブートストラップ。
4.2

引用文献・サイト

参考書・サイト

1) Bolker BM, 2008. Ecological models and data in R.

個人的に役に立つと思う本とウェブサイトをあげておPrinceton University Press.

きます。

2) ドブソン AJ, 2008. 一般化線形モデル入門, 原著 第 2版（Dobson AJ, 2002. An introduction to generalized• 本linear models-2nd edition. Chapman & Hall. 田中豊・– Ecological models and data（Bolker 2008）1) : と森川敏彦・山中竹春・冨田誠 訳）. 共立出版.にかく最尤推定一筋な本です。すぐに役立つ知

3) 伊庭幸人・種村正美・大森裕浩・和合肇・佐藤整尚・

識は載っていませんが，自身の技術を底上げす高橋明彦, 2005. 計算統計 II マルコフ連鎖モンテカルることができます。

ロ法とその周辺. 岩波書店.

– 一般化線形モデル入門（ドブソン 2008）2) : 長

4) KuboWeb.

らく絶版で入手困難でしたが，最近第 2 版が翻訳されました。あまり系統立てて書かれていな

5) McCullagh P & Nelder JA, 1989. Generalized linear

い気がしますが，分かりやすいです。

models-2nd edition. Chapman & Hall.

– 計算統計 II（伊庭ら 2005）3) : GLM とは直接関

6) RjpWiki. http://www.okada.jp.org/RWiki/

係ありませんが，良い統計の本だと思います。

7) サルツブルグ D, 2006. 統計学を拓いた異才たち（Sals-– 統計学を拓いた異才たち（サルツブルグ 2006）7)

http://hosho.ees.hokudai.ac.jp/∼kubo/index-j.html

burg D, 2001. The lady tasting tesa. How statistics revo-: 統計の歴史を楽しく読めます。

lutionized science in the twentieth century. Henry Holt

– 統計モデル入門（丹後 2000）8) : 知りたいこと& Company. 竹内惠行・熊谷悦生 訳）日本経済新聞社.

はだいたい書いてある良い本だと思いますが，

8) 丹後俊郎, 2000. 統計モデル入門. 朝倉書店.

難解です。

9) ヴェナブルズ WN・リプリー BD, 2001. S-PLUS に– S-PLUS による統計解析（ヴェナブル・リプリー

よる統計解析, 原著 第 3 版（Venables WN & Ripley2001）9) : S-PLUS のコマンドが書いてありますBD, 1999 Modern applied statistics with S-PLUS-3rdが，R でも使えるので役に立ちます。難解です。

edition. Springer-Verlag. 伊藤幹夫・大津泰介・戸瀬信– 生物統計学入門（山田・北田 2004）10) : 確率分之・中東雅樹 訳）. シュプリンガー・フェアクラーク

布の話がきちんと書いてある良い入門書だと思東京.

います。

10) 山田作太郎・北田修一, 2004. 生物統計学入門. 成山• ウェブサイト堂書店.

– KuboWeb4) : 生態学者なら（でなくても）はずせ6