Embedded Files
回帰分析とは、連続変数同士の関係を分析する手法の一つであり、ANOVAと並ぶ線形モデルの基本である。
ANOVAと異なるのは、説明変数がカテゴリ変数ではなく、連続変数なことである。
説明変数が1つ、目的変数も1つであるような回帰分析を、特に単回帰と呼ぶことがある。
要点
要点
変数間の因果関係に気をつけるべし（目的変数は結果、説明変数は原因）
回帰係数は説明変数がもつ影響の大きさを、決定係数はモデルの当てはまりを示す
野外観察データでは、因果関係の解釈に慎重になるべし
式変形が有効な場面がある
回帰分析の強みの一つは、定量的な結論を下せること
例題1
例題1
野外観察
野外観察
アゲハはミカン科植物を寄主植物として利用している（幼虫がミカン科植物を食べて育つ）。アゲハ成虫の個体群密度が、その場所のミカン科樹木の数によってどのように変化するか調べるため、15か所の公園でトランセクト調査を行い、kmあたり成虫個体数を記録した。また、公園ごとのミカン科樹木の数も調査した。その結果、以下のようなデータが得られた。
> ButterflyData
Site.id Butterfly.Abundance Tree.Abundance
1 1 13.36 12
2 2 8.41 9
3 3 15.64 11
4 4 18.38 15
5 5 4.39 5
6 6 10.73 8
7 7 6.98 8
8 8 7.99 8
9 9 12.33 14
10 10 16.54 11
11 11 18.30 14
12 12 20.18 13
13 13 6.70 6
14 14 14.07 14
15 15 13.33 10
散布図の作成（2つの変数のどちらをx軸に、どちらをy軸に取るか？）
散布図の作成（2つの変数のどちらをx軸に、どちらをy軸に取るか？）
まずは作図してデータの傾向を掴む必要があるのは、回帰分析の際でも同様である。
2つの連続変数間の関係は、散布図によって図示することができる。散布図を描くにあたって気をつけるべき重要なポイントは、2つの変数間の因果関係である。散布図の基本は「因果関係の原因をx軸に、結果をy軸に取る」ことである。今回の調査の目的は、「アゲハ成虫の個体群密度が、その場所のミカン科樹木の数によってどのように変化するか調べる」ことであり、ミカン科樹木の数→アゲハ成虫の個体群密度、という因果関係を暗黙に仮定している。したがって、ミカン科樹木の数をx軸に、アゲハ成虫の個体群密度をy軸に取るのが適切である。
この大前提を踏まえて、今回のデータの散布図を作成してみよう。
散布図は、plot関数によって描くことができる。boxplot関数と同じように、”y軸 ~ x軸”のようなモデル式をplot関数に入れればよい。
par(las=1, lwd=2, family="sans") #lasは軸ラベルの向きを指定、lwdは線の太さを指定、familyはフォントを指定する引数（sansならサンセリフ体になる）
plot(Butterfly.Abundance ~ Tree.Abundance, ButterflyData, pch=16, cex=2) #pchは散布図のシンボルを指定（16なら黒丸になる）、cexは散布図のシンボルの大きさを指定する引数
右図のような散布図を描くことができる。
寄主植物の数が増えるほど、チョウの個体群密度も直線的に増加していることが読み取れる。
回帰分析では、「説明変数と目的変数は直線関係にある」というかなり強い仮定をおく必要があるため、事前にその仮定がどの程度満たされているか確かめておくのは重要である。（直線関係でなければ回帰分析してはいけないということではなく、解釈に注意が必要になるということ）
単回帰の実践
単回帰の実践
回帰分析はlm関数によって行うことができる。
まずlm関数によって回帰モデルを作成し、次にそのモデルをsummary関数で囲むことで、結果の概要を出力できる。lm関数内では、モデル式"目的変数 ~ 説明変数"を指定する。
なお、ここでも重要なのは、2つの変数間の因果関係である。回帰分析は、説明変数が目的変数に与える影響を検証する手法（基本的には）なので、説明変数は因果関係の原因、目的変数は因果関係の結果と解釈される。したがって、今回の場合には、説明変数にはTree.Abundanceを、目的変数にはButterfly.Abundanceを指定する必要がある（plot関数内と同じである）。
m <- lm(Butterfly.Abundance ~ Tree.Abundance, ButterflyData) #モデル式、指定するデータ名の順に指定する
summary(m)
> m <- lm(Butterfly.Abundance ~ Tree.Abundance, ButterflyData) #モデル式、指定するデータ名の順に指定する
> summary(m)
Call:
lm(formula = Butterfly.Abundance ~ Tree.Abundance, data = ButterflyData)
Residuals:
Min 1Q Median 3Q Max
-4.7312 -1.6067 -0.0002 1.5638 4.4378
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.4048 2.4265 -0.579 0.573
Tree.Abundance 1.3190 0.2214 5.958 4.76e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.598 on 13 degrees of freedom
Multiple R-squared: 0.7319, Adjusted R-squared: 0.7113
F-statistic: 35.5 on 1 and 13 DF, p-value: 4.762e-05
出力のうちまず重要なのは、CoefficientsのTree.Abundance行である。Tree.AbundanceのEstimateは「ミカン科樹木の数の回帰係数（の推定値）」であり、「ミカン科樹木の数が1増えたときに、アゲハ成虫の個体群密度が変化する値」、言い換えれば「ミカン科樹木の数がアゲハ成虫の個体群密度に与える影響の大きさ」を表している。この係数の符号がプラスであれば正の影響を、マイナスであれば負の影響を表している。さらに、一番右のPr(>|t|)はP値のことで、隣のt value（t値）を使って、「回帰係数が有意にゼロと異なるか」を検定している。今回の場合、P値は0.001以下なので、回帰係数はゼロより有意に大きい、と解釈することができる。
今回の調査の問いは「アゲハ成虫の個体群密度が、その場所のミカン科樹木の数によってどのように変化するか」だった。回帰係数の結果を見る限り、この問いに対しては「ミカン科樹木の数はアゲハ成虫の個体群密度に対して有意な正の影響を与えた。具体的にはミカン科樹木の数が1増えると、アゲハ成虫の個体群密度は1.32増えた」と結論づけられる。
ただし、回帰分析でわかることはそれだけではない。もう一つ重要なのが、Multiple R-squaredである。これは、この線形モデルのR2値（決定係数）のことで、これは「総変動の何割をミカン科樹木の数が説明したか」を表す、モデルの当てはまりの指標の一つである。たとえミカン科樹木の数が有意な影響を持っていたとしても、それがアゲハ成虫個体群密度に最も影響する要因であるかはまた別の話である。もしR2値が0.1や0.2だったら、その説明変数はとても重要な要因であるとは言い難い。今回の場合、R2値は0.73なので、総変動の73パーセントをミカン科樹木の数が説明した、と解釈することができ、ミカン科樹木の数はある程度重要な要因だろうと考えることができる。
隣にあるAdjusted R-squaredは自由度調整済み決定係数と呼ばれるもので、説明変数の数が多いほど1に近づきやすくなる決定係数の性質を補正したものである。なお、以下のように作成したモデルをanova関数で囲むことによっても、変数の有意性を検証することができる。
> anova(m)
Analysis of Variance Table
Response: Butterfly.Abundance
Df Sum Sq Mean Sq F value Pr(>F)
Tree.Abundance 1 239.624 239.624 35.497 4.762e-05 ***
Residuals 13 87.758 6.751
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
単回帰の場合、anova関数によって行われるF検定は、回帰係数の有意性を検証するt検定（summary関数のアウトプット）と完全に同値である。
回帰直線の図示
回帰直線の図示
回帰分析の結果を図示する典型的な方法は、実データの散布図に回帰直線を上書きすることである。
回帰直線とは、説明変数の値によって目的変数を推定or予測するための直線である。
言葉で説明するよりも図を見た方が早いので、まずは描いてみる。
回帰直線を描く工程は以下である。
回帰線を描くために、説明変数だけを含む新しいデータを作成する。
上のデータを用いて、predict関数によって予測値とSEを計算させる。
計算した予測値とSEから、予測値±SE（95%信頼区間を求めたければ、予測値±SE×1.96）を計算し、格納する。
lines関数によって実データプロットに上書きする。
以上の工程は、以下のスクリプトによって実行できる。ここでは、回帰線とその95%信頼区間を描き加えている（SEの約1.96倍が95%信頼区間の上/下限）。
（'20 2/4追記：上記のSE×1.96は、データ数が大のときのみ有効かも。データ数大なら平均値の残差は近似的に正規分布に従うが、データ数小のときは平均値の残差（の1/SE）はt分布に従うので、SE×t(dferror, 0.975)（dferror=残差自由度）にしないと不正確。Rではse.fit*qt(dferror, 0.975)#dferrorは残差自由度 で幅を計算できるし、predict(m, new.data, interval="confidence")とすれば残差がt分布に従うと仮定した95%信頼区間を計算してくれる）
#1
new.data <- data.frame(
Tree.Abundance=min(ButterflyData$Tree.Abundance):max(ButterflyData$Tree.Abundance)) #minおよびmax関数は、数列内の最小値と最大値を返す
new.data
#2
pred <- predict(m, new.data, se.fit=TRUE) #5本から15本までのアゲハ個体群密度を予測
pred
#3
pred.mean <- pred$fit #予測値の格納
pred.li <- pred$fit-pred$se.fit*1.96 #信頼区間下限値の格納
pred.ui <- pred$fit+pred$se.fit*1.96 #信頼区間上限値の格納
##これらをそれぞれx座標、y座標の値としてプロットする
new.data$Tree.Abundance
pred.mean
pred.li
pred.ui
#4
par(las=1, lwd=2, family="sans")
plot(Butterfly.Abundance ~ Tree.Abundance, ButterflyData, pch=16, cex=2)
lines(new.data$Tree.Abundance, pred.mean, lty=1) #予測値の上書き
lines(new.data$Tree.Abundance, pred.li, lty=2) #信頼区間下限の上書き
lines(new.data$Tree.Abundance, pred.ui, lty=2) #信頼区間上限の上書き
> #1
> new.data <- data.frame(
+ Tree.Abundance=min(ButterflyData$Tree.Abundance):max(ButterflyData$Tree.Abundance))
> new.data
Tree.Abundance
1 5
2 6
3 7
4 8
5 9
6 10
7 11
8 12
9 13
10 14
11 15
>
> #2
> pred <- predict(m, new.data, se.fit=TRUE) #5本から15本までのアゲハ個体群密度を予測
> pred
$fit
1 2 3 4 5 6 7 8
5.190184 6.509187 7.828190 9.147193 10.466196 11.785198 13.104201 14.423204
9 10 11
15.742207 17.061210 18.380213
$se.fit
1 2 3 4 5 6 7 8
1.3966683 1.2071842 1.0304995 0.8744078 0.7518465 0.6811618 0.6787591 0.7452991
9 10 11
0.8650154 1.0193414 1.1949420
$df
[1] 13
$residual.scale
[1] 2.598192
>
> #3
> pred.mean <- pred$fit #予測値の格納
> pred.li <- pred$fit-pred$se.fit*1.96 #信頼区間下限値の格納
> pred.ui <- pred$fit+pred$se.fit*1.96 #信頼区間上限値の格納
>
> ##これらをそれぞれx座標、y座標の値としてプロットする
> new.data$Tree.Abundance
[1] 5 6 7 8 9 10 11 12 13 14 15
> pred.mean
1 2 3 4 5 6 7 8
5.190184 6.509187 7.828190 9.147193 10.466196 11.785198 13.104201 14.423204
9 10 11
15.742207 17.061210 18.380213
> pred.li
1 2 3 4 5 6 7 8
2.452714 4.143106 5.808411 7.433353 8.992576 10.450121 11.773833 12.962418
9 10 11
14.046777 15.063301 16.038127
> pred.ui
1 2 3 4 5 6 7 8
7.927654 8.875268 9.847969 10.861032 11.939815 13.120276 14.434569 15.883990
9 10 11
17.437637 19.059119 20.722299
>
> #4
> par(las=1, lwd=2, family="sans")
> plot(Butterfly.Abundance ~ Tree.Abundance, ButterflyData, pch=16, cex=2)
> lines(new.data$Tree.Abundance, pred.mean, lty=1) #予測値の上書き
> lines(new.data$Tree.Abundance, pred.li, lty=2) #信頼区間下限の上書き
> lines(new.data$Tree.Abundance, pred.ui, lty=2) #信頼区間上限の上書き
上記の工程はやや複雑なので、アウトプットを一度よく眺めて自分なりによく理解することを勧める。
いずれにしろ、以上のスクリプトによって、右のような図を描くことができる。樹木の数が8ならアゲハ個体群密度は約9、樹木が14ならアゲハは約17、というのが回帰直線の予測だということが読み取れる。
点線の区間は回帰直線の95%信頼区間で、これは真のアゲハ成虫密度（の平均値）の回帰直線による推定精度を表していると考えてよい※。回帰直線による推定は限られた値域のデータから行われているので、直線の両側の端は真ん中に比べてどうしても推定精度が悪くなる。したがって、真ん中に比べて両側の端が区間が広くなるのが、回帰直線の信頼区間の特徴である。
※正確には、同じ方法の野外調査と回帰分析を100回行った場合に、そのうちの95回はこの信頼区間内に真のアゲハ成虫密度の平均値が含まれると期待される、という区間。
野外調査における結果の解釈
野外調査における結果の解釈
t検定やANOVAに比べると、回帰分析は野外観察データで使う機会が多い手法だ（と思う）。理由はおそらく、連続変数を実験的に操作するのはカテゴリ変数に比べて労力がかかるからだろう（ミカン科樹木の数を6、8、10、12、、、と何段階にも操作するのは現実的ではない）。
ここで注意すべきなのは、一般的に、野外観察データの場合、変数間の因果関係を証明するのが困難だということである。もしかしたら、植生が豊か→ミカン科樹木が多い、植生が豊か→アゲハ成虫が多い、という因果関係が存在し、それが理由でミカン科樹木が多いとアゲハ成虫が多くなったのかもしれない。本当にミカン科樹木の数がアゲハ成虫の個体群密度に影響していたのか？という問いは、上述の単回帰では答えることのできない問いである。
今回の調査では、ミカン科樹木の数とアゲハ成虫の個体群密度以外に調査されている項目はなく、前者から後者への因果関係を立証する材料は乏しい。よって、「ミカン科樹木の数はアゲハ成虫の個体群密度に対して有意な正の影響を与えた」という結論はやや言い過ぎかもしれず、「両者の間には有意な正の関係があった」程度の結論の方が妥当かもしれない。野外観察データでは、因果関係の解釈は慎重に、というのが原則である。
回帰係数と回帰直線
回帰係数と回帰直線
単回帰では、回帰係数は2つ推定される。
1つは、説明変数の回帰係数（上で説明した「Tree.Abundanceの回帰係数」）であり、もう1つは、切片（intercept）の回帰係数である。
summaryの出力と右の図を見比べてみると、まず、Tree.Abundanceが0のときの回帰直線のy座標、すなわちこの直線のy切片が-1.40である。これは、Interceptの回帰係数-1.4048と一致している。
次に、Tree.Abundanceが0から10に増えた時のButterfly.Abundanceの増加分は11.79-(-1.40)=13.19となっている。よって、この直線の傾きは13.19/10=1.32である。これは、Tree.Abundanceの回帰係数1.3190と一致している。
このように、単回帰の2つの回帰係数は、回帰直線のy切片と傾きに対応している。
式で書けば以下のようになる。
Y = β0 + β1 × X
Y = β0 + β1 × X
ここで、β0は切片の回帰係数、β1は説明変数の回帰係数、Yは目的変数、Xは説明変数である。
残差と回帰分析の診断
残差と回帰分析の診断
残差とは、回帰直線と実データとの差分のことである。
例題1の野外観察データを例にとると、残差は右のような図で表される。
回帰分析では、残差の2乗の総和（残差の二乗和）が最も小さくなるような回帰係数を求めて回帰直線を引く（実はこの二乗和の値は、anova関数の出力にあるResidualsのSum Sqのことである）。
回帰分析で検定や母数推定を行う際には、残差の母集団について、1. 正規分布、2. 等分散（Tree.Abundanceの値が違ってもバラツキの程度が等しい）が仮定される。
このうち残差の正規性については、残差のヒストグラムを描くことで簡便にチェックできる。
m <- lm(Butterfly.Abundance ~ Tree.Abundance, ButterflyData)
hist(resid(m), main="", xlab="Residuals")
今回の場合、とりあえず一山型の形になっており、正規分布と極端に異なるようなことはなさそうに見える。
上記の例のように、回帰分析では、事前に正規性・等分散性の仮定を吟味するというよりは、まず線形モデルの当てはめを行ってから、残差を観察することで仮定がどの程度満たされていそうか診断することが多い。これは、ANOVAに比べて、生データを見ただけでは直感的に正規性・等分散性が満たされているかの判断が難しいことによると思う。
Rでは、作成した線形モデルを以下のスクリプトに入れることで、これらの仮定について診断する機能が搭載されている。
m <- lm(Butterfly.Abundance ~ Tree.Abundance, ButterflyData)
par(mfrow=c(2, 2)) #以下の図はpar(mfrow=c(2, 2), lwd=2, family="sans")と指定している
plot(m)
右のような図が出力される。内容は以下である。
Residuals vs Fitted: 説明変数と目的変数の関係が、真っ直ぐか（赤い線が水平であればあるほど良い）
Normal Q-Q: 正規性が満たされているか（白丸が点線の上に直線状に並ぶほど良い）
Scale-Location: 等分散性が満たされているか（赤い線が水平であればあるほど良い）
Residuals vs Leverage: 外れ値が存在するか（赤い点線より外側にある点が外れ値の恐れ）
最初の3つが特に重要である。
Normal Q-Qがちょっと理解しづらいのだが、こちらのページでとても分かりやすく解説されている。
今回の結果の場合、1の真っ直ぐさと2の正規性はある程度保たれていそうだが、3の等分散性については、目的変数の値が大きくなるほど残差が大きくなるような傾向がわずかに観察されている。
なお、診断の結果の良し悪しを判断する絶対的な基準は存在しないが、モデルの改善のやり方としては、変数変換やGLM、GLSによる異分散の仮定などが考えられるだろう。
イメージがつきやすいように、真っ直ぐでないデータ（nonlinear）、正規分布でないデータ（nonnormal）、等分散でないデータ（heteroscedastic）について、それぞれで回帰直線とResiduals vs Fitted、Normal Q-Q、Scale-Locationを比較してみよう。
まず、"Nonlinear"の場合でも、強引に真っ直ぐな回帰直線が引かれていることがわかる。しかし、回帰直線とデータとの差分は、直線上のy座標の値が小さいうちは下ぶれし、真ん中あたりでは上ぶれ、大きくなるとまた下ぶれしている。Residuals vs Fittedのプロットはこのような実データの傾向をよく反映している。
次に、"Nonnormal"の場合、回帰直線は実データにしっかり乗っているように見えるが、よく見ると直線付近にはほとんどデータ点が無く、直線よりやや離れた両側に点が集中している（要するに二山型データ）。このような場合、Normal Q-Qのプロットは右の図のように”途切れて”しまう。
最後に、"Heteroscedastic"の場合、回帰直線のy座標の値が大きくなるほど、データ点が直線の両側に大きく散らばるように（言い換えれば、残差が大きく）なっている。この傾向は、Scale-Locationにおける右肩上がりの赤線によく表れている。
このうち、Nonlinearについては、多項式（2乗の項などを含むお出る）回帰やGAMによる非線形回帰が有効そうだ。また、Heteroscedasticについては、lnなどの変数変換、GLSによる重みづけ回帰、GLMなどが有効だろう。Nonnormalについては統計で何とかするというよりは、考慮に入れるべき他の要因の見落としがないか疑うべきだろう（直線より上側と下側のデータでは、何かの要因が異なる可能性がある）。
さて、次の例では、因果関係が明確な操作実験のデータを取り上げ、もう少し複雑な解析を行ってみよう。
例題２
例題２
操作実験
操作実験
ある昆虫の幼虫の発育がどの程度の低温で停止するのか調べるため、幼虫を温度条件16, 18, 20, 22, 24, 26℃の6グループに分けて飼育し、孵化から蛹化までにかかった日数を記録した。その結果、以下のようなデータが得られた。
> DevData
LarvalID Temp DevTime
1 16_1 16 67
2 16_2 16 51
3 16_3 16 53
4 16_4 16 63
5 16_5 16 42
6 18_1 18 39
7 18_2 18 32
8 18_3 18 28
9 18_4 18 36
10 18_5 18 35
11 20_1 20 25
12 20_2 20 27
13 20_3 20 25
14 20_4 20 29
15 20_5 20 24
16 22_1 22 25
17 22_2 22 23
18 22_3 22 25
19 22_4 22 23
20 22_5 22 23
21 24_1 24 20
22 24_2 24 23
23 24_3 24 23
24 24_4 24 18
25 24_5 24 23
26 26_1 26 17
27 26_2 26 19
28 26_3 26 18
29 26_4 26 18
30 26_5 26 20
真っ直ぐでないデータ
真っ直ぐでないデータ
par(lwd=2, family="sans", las=1)
plot(DevTime ~ Temp, DevData)
今回のデータの場合、温度は実験者が操作しており、温度→日数という因果関係が明らかである。
温度をx軸に、日数をy軸にとってデータをプロットすると、右のように下に凸の真っ直ぐでない関係を見て取ることができる。
一番低い温度である16℃でも、発育が停止してはいないようである。
DevTime.m <- lm(DevTime~Temp, DevData)
par(mfrow=c(1, 2), lwd=2, family="sans", las=1)
plot(DevTime~Temp, DevData)
abline(DevTime.m) #abline関数は、事前に作成した線形モデルを指定することで、簡便に回帰直線を上書きできる関数
plot(DevTime.m, which=1) #which=1と指定することで、Residuals vs Fittedのみ呼び出すことができる
このデータに対して単純に単回帰を行うと、右のような回帰直線が描かれるが、やはり回帰直線はデータにはイマイチ当てはまっていない。
おまけに、この直線を引いたところで、この昆虫の発育がどの程度の低温で止まるのかは全くわからない（直線を左に伸ばし続けても、有限の発育日数の値の範囲から逸脱することはない）。
式変形により、直線関係を示す形に持ち込む
式変形により、直線関係を示す形に持ち込む
このようなデータに対処する方法の一つは、適切な式変形によって、直線関係を示す形に持ち込むことである。
実は、一般的に昆虫の発育は以下のような反比例の式によって表されることが経験的にわかっている。
ここで、Dは発育日数、Tは温度である。T0は”発育零点”と呼ばれる定数で、TがT0であるとき、発育日数は無限大となり、発育が停止する温度と解釈することができる。
この式は、両辺の逆数を取るという非常に単純な式変形によって、以下のような直線の式にすぐに持ち込むことができる（re-parameterizationの入門になりそう）。
温度と発育日数の逆数が直線関係にあることは、以下のように温度と発育日数の逆数を散布図でプロットすることですぐに確かめることができる。ここで、関数Iは、式変形の加減乗除をモデル式における加減乗除から区別するための関数である。
par(lwd=2, family="sans", las=1)
plot(I(1/DevTime)~Temp, DevData)
したがって、発育日数の逆数に対して温度による単回帰を行えば、−T0 / Kと1 / Kの値を求めることができる。
これを行うためには、以下のようにモデル式の左辺をI(1/DevTime)と入力すればよい。
DevRate.m <- lm(I(1/DevTime) ~ Temp, DevData)
summary(DevRate.m)
> DevRate.m <- lm(I(1/DevTime) ~ Temp, DevData)
> summary(DevRate.m)
Call:
lm(formula = I(1/DevTime) ~ Temp, data = DevData)
Residuals:
Min 1Q Median 3Q Max
-0.0067518 -0.0027259 0.0002303 0.0020799 0.0073184
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0320722 0.0045652 -7.025 1.22e-07 ***
Temp 0.0033593 0.0002146 15.656 2.23e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.004014 on 28 degrees of freedom
Multiple R-squared: 0.8975, Adjusted R-squared: 0.8938
F-statistic: 245.1 on 1 and 28 DF, p-value: 2.228e-15
ここで推定されるInterceptの回帰係数（β0）とTempの回帰係数（β1）が、−T0 / Kおよび1 / Kと対応している。
すなわち、
なので、
となる。
後は、推定した回帰係数の値を上式に代入すれば、発育が停止する温度である発育零点T0を求めることができる。
回帰係数は、関数coefで線形モデルを囲むことで呼び出すことができる。
co <- coef(DevRate.m)
co
-co[1]/co[2]
> co <- coef(DevRate.m)
> co
(Intercept) Temp
-0.032072242 0.003359339
> -co[1]/co[2]
(Intercept)
9.54719
以上から、発育零点は約9.55℃と求められる。したがって、「この昆虫の発育が止まるのは、約9.55℃である」と結論することができる。
以下のようなスクリプトにより、単回帰によって求めた反比例関係が実際にデータに沿っているか目で確認することができる。
y <- function(x) k*(1/(x-t))
k <- 1/co[2]
t <- -co[1]/co[2]
par(lwd=2, family="sans", las=1)
plot(DevTime~Temp, DevData, xlim=c(5, 30), ylim=c(0, 180))
curve(y, add=TRUE, xlim=c(9.55, 30))
abline(v=t, lty=2, lwd=1)
右の図のように、曲線はしっかり実データに沿っているようである。さらに、9.55℃付近の垂直線に沿って曲線が漸近的に増加している様子が観察できる。したがって、この昆虫の発育が9.55℃付近で停止するのはある程度正しそうである。
以上のように、回帰分析の強みの一つは、回帰係数の値を適切に解釈することで、差が有る無いだけではない、定性的な結論を下せることである。
'20 1/12
Google Sites
Report abuse