非線形モデルの当てはめ

最小自乗法

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

線形モデル

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

また、これによって分散分析法を実験計画法から切り離すことが可能となり、畜産分野で頻繁に見られる不揃いデータの分散分析が容易になった(最小自乗分散分析法)。さらに、実験計画法から切り離されたことで、家畜育種分野では後述する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