佐賀大学農学部 生物生産学科 動物生産学分野 和田研究室

非線形モデルの当てはめ

最小自乗法

推定されたモデルから得られた予測値と実際のデータとの差を残差と言う。 残差の自乗和を最小にするように、モデルのパラメータの推定を行う手法を 最小自乗法という。通常の回帰分析、重回帰分析などの計算式は、この 最小自乗法によって得られたものである。

線形モデル

実は回帰分析と分散分析を統一的に理解することができる。 どちらのモデルも線形であり、計算手法も最小自乗法を用いているからだ。 回帰分析と分散分析を統一的に理解することによって、より複雑で 新しい統計手法への理解が可能となる。

また、これによって分散分析法を実験計画法から切り離すことが可能と なり、畜産分野で頻繁に見られる不揃いデータの分散分析が容易に なった(最小自乗分散分析法)。さらに、実験計画法から切り離された ことで、家畜育種分野では後述するBLUP法などのフィールドデータを 用いた分析法へと発展した。

ただし、すべての問題が最小自乗法で解けるわけではない。BLUP法などの 混合モデルでは後述する最尤法によらなければならない。

(参考)線形モデルの最小自乗解
直線回帰モデル
y = a + xb + e
ここで、yは目的変数ベクトル、xは説明変数ベクトル、 eは誤差ベクトル、aとbはパラメータで、aは定数項、bは回帰係数 です。

このモデルを次のように書き直します。

y = 1a + xb + e

ここで1は全ての要素が1のベクトルです。さて、1aとxb をひとつにしてしまいましょう。

y = Xc + e

X= ( (1, x1), (1, x2), ... ), c= ( a, b ) ですね。

残差は e = y - Xc ですから、残差自乗和は

= (y - Xc)'(y - Xc)

= y'y - 2c'X'y + c'(X'X)c

微分して=0とおけば残差自乗和を最小値にするパラメータの値が求められます。

-2 X'y + 2 X'Xc = 0

X'Xc = X'y

この式を正規方程式と呼び、この連立1次方程式をcについて 解けば、パラメータの推定値を求めることができる。 c^ = (X'X)-1X'y

非線形モデル

モデルを各パラメータで偏微分したときに、一次偏導関数がパラメータ自身を含むモデルを非線形モデルと呼ぶ。 例えば、logistic曲線などの成長曲線モデルが非線形モデルである。

非線形モデルを最小自乗法で当てはめる時には、多くの場合、計算式を 解析的に求めることができない。そこで、多くの反復解法が考案されて きた。通常はニュートン法とその改良手法を使えば十分であるが、非線形性 の強い場合にはマルカート法やPowell法を使う必要がある。また、偏微分 すら求められない場合には、数値微分法をもとにしたアルゴリズムを 使用しなければならない。

例題

対照閉鎖群集団(RR)の58世代目において、生時から15週齢までの 平均体重のデータについて、統計言語Rを用いてBrodyの 成長曲線を当てはめてみる。

Excelで次のデータを入力する。 このデータを「ファイルの種類」を「テキスト(タブ区切り)」にして growth.txtと名前をつけて保存する。

Rを起動してgrowth.txtを読み込む。

> gdata <- read.table("growth.txt")
> gdata
   Age Weight
0    0    6.7
1    1   17.0
2    2   34.5
3    3   56.5
4    4   80.3
5    5   94.6
6    6  100.7
7    7  103.7
8    8  104.9
9    9  103.5
10  10  104.8
11  11  105.7
12  12  106.6
13  13  108.0
14  14  107.9
15  15  107.3
> 
非線形回帰パッケージnlsとgdataを使用するという宣言をする。
> library(nls)
> attach(gdata)
パラメータの初期値をA=100.0, B=1.0, k=0.1として計算する。
> result <- nls( Weight ? A*(1.0-B*exp(-k*Age)), data=gdata, start = list( A=100.0, B=1.0, k=0.1))
結果のサマリーを表示する。
> summary(result)

Formula: Weight ? A * (1 - B * exp(-k * Age))

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
A 112.63730    3.85214  29.240 3.02e-13 ***
B   1.03431    0.05552  18.631 9.27e-11 ***
k   0.28909    0.03906   7.401 5.19e-06 ***
---
Signif. codes:  0  `***'  0.001  `**'  0.01  `*'  0.05  `.'  0.1  ` '  1 

Residual standard error: 7.194 on 13 degrees of freedom

Correlation of Parameter Estimates:
        A      B
B -0.3012       
k -0.7871 0.5903

成熟体重(A)は112.64、成熟速度(k)は0.2891、積分定数は1.0343と推定 されました。 寄与率を計算してみます。
> predict <- Weight + resid(result)
> cor( Weight, predict)*cor( Weight, predict)
[1] 0.9682295
> 


家畜家禽の発育様相の改良

最終更新年月日 2004年6月4日

ホームページ

ywada@cc.saga-u.ac.jp