直線を描いたら曲線を描きたい

　これまでの様々な線形回帰では、あくまでも直線関係＝線形関係（厳密には正しくないが）で表されるデータに関して解析を行ってきた。しかし、実際の世界ではそんな関係のデータは少数だろう。例えば、植物の高さと葉面積の関係を調べているとしよう。高さは1次元だが、面積は2次元のデータだ。それゆえ、高さを説明変数として葉面積を被説明変数とすると、2次式の関係になる可能性が高い。そうなってくると、もはやこれまでのただの線形回帰ではデータの様子を正しく推定することは不可能となるだろう。直線を描けるようになったら、曲線関係を描きたくなる、そういうものである。これに対する解決方策は3つある。1つ目は今回紹介する、非線形な関係をモデル式に取りこんで最小二乗法をおこなう回帰、2つ目は一般化線形モデル、3つ目は一般化加法モデルである。一般化線形モデルは最尤推定、一般化加法モデルは最小二乗推定なので、今回紹介する方法は一般化加法モデルにどちらかといえば近い。

非線形最小二乗法

　今回紹介する方法は、例え、説明変数が非線形に被説明変数に影響を与えているとしても、発想は線形回帰と同じく、「回帰曲線」と被説明変数の偏差の二乗和が最小になるように「回帰曲線」を推定する方法である。線形回帰の時の仮定と同じく、この偏差は平均0の正規分布に従い、説明変数によらず分散が一定であることを仮定している。どんな関係式が非線形になるかを、特に生物学において登場する機会の多い例を交えながら紹介しよう。

・多項式回帰

　以下のように、モデルが説明変数の多項式で表される場合、多項式回帰polynomial regressionと呼ぶ。単回帰は多項式回帰の特殊な例ということになる。例えば、植物高を説明変数としたとき、葉面積の総和は2次関数の関係になることが期待される。

一般の多項式のグラフを思い浮かべてもらうとわかるように、多項式回帰をした場合、曲がったデータの関係を表現することが可能である。2乗項を含めば一山or一谷の関係、3乗項を含めば一つの山と一つの谷を持った関係……、というように次数の大きな項を含むほど複雑な関係を表現できる。ただし、どれくらい大きな次数まで含むかは、モデルを組む時にあらかじめ決めておく必要がある。回帰係数は計算機が計算してくれるパラメータだが、このような人間側が決めなければいけないパラメータをハイパーパラメータhyperparameterと呼ぶ（その点で、回帰するモデルを選択する過程もハイパーパラメータである）。大きな次数を含むほど複雑な関係を表現できる一方で、解釈はより難しくなるし、手元のデータを表現することに重きが置かれて、新しいデータに適応できないようなモデルが作られる（これを過学習overfittingという）問題も付きまとう。こういうときにはモデル選択（例1、例2）によってどの次数まで含むかを検討する必要が出てくる。また、各項x、x^2、x^3……は互いに強く相関するので、多重共線性の問題が生じうることも念頭に置いておく必要がある。1次ではなく、2次や3次の関係が強い影響を及ぼしうることが理論的に明白な場合は2次の項、3次の項だけを説明変数に加えることも考えよう。

　では具体的な解析例を紹介する。以下のようなデータを解析してみよう。

------------------------------------------------------

library(plotn)

x1 <- c(15.9, 14.8, 2.3, 21.5, 19.8, 2.3, 13.4, 23.1,

        -8.2, 24, 23, -5.9, -1.3, 28.9, 1.4, -3.6, 12.1,

        9.9, -3.4, -2.8, 17.9, -4.6, -5.2, 29, 25.5, 2.8, 

        29.1, 4.6, -2.9, 8.7, 24.1, 25.7, 3.3, -5.1, 13.2,

        -2.7, 21.2, 29.2, 6.3, 22.2, 1.5, 16.3, 5.4, 25.4,

        6.7, 5.3, 0.9, 14.1, 8.6, -9, 14.9, 29, 13.7, -2.2,

        8.2, 24.9, 27.2, 25.7, 8.8, 7.3, 2, -4, 25.7, 25.1,

        18.9, -2.7, 15.7, 26.7, 19, 8.6, -4.1, 26.5, 11.2, 

        8.8, 0.5, -0.1, -7.6, 7.2, -5.7, -8.7, 1.3, 4.2, 

        -7.8, 7.2, 27.3, -3.1, 27.8, -8.4, 24.3, -0.2, -3.9,

        6, 5.6, -0.3, 1.8, 12.8, 16.9, 9.4, -0.2, 3.6)

y <- c(-26.8, -24.1, 0.7, -33, -26.1, -1.2, -27.3, -40.8,

       43, -36.8, 4.2, 36.4, 16.2, 1.4, 32.7, 16, -18.5, 

       -29.6, 7.9, 12.3, -25.3, 5.1, 27.3, -22.5, 11.9, -25.9, 

       16.9, -40.9, 46.1, -33, -14.3, -20.9, -48.1, 32.2, -34.6, 

       2.8, -43.7, -12, -35.5, -68.4, -12.6, -60.1, -7, -47.1,

       -24.7, -0.7, -4.4, -49.7, -22, 101.9, -55.2, 21.6, -42.2, 

       27.2, -34.2, -15.3, -10.1, 1.9, -42.2, -19.6, 17.4, 46.8, 

       -26.5, 10.3, -20.9, 3.8, -37.4, -5.1, -59.9, -41.2, 13.3, 

       -47.6, -41, -33.2, 8.8, 12.4, 34.2, -79.5, 69.3, 41.4, 17.9, 

       -25.9, 70.4, -39.8, -8.1, 40.3, -21.1, 39.7, 1.6, 14.8, 

       -0.5, -38, -28.9, 8.6, -14.8, -36.1, -24.1, -64.4, 9, 10)

d <- data.frame(x1 = x1, y = y)

head(d)

##     x1     y

## 1 15.9 -26.8

## 2 14.8 -24.1

## 3  2.3   0.7

## 4 21.5 -33.0

## 5 19.8 -26.1

## 6  2.3  -1.2

plotn(y ~ x1, data = d)#図1の描画

------------------------------------------------------

データはどうやら15付近を最小値とする曲線関係にあるようだ。2次以上のどんな項で解析してもいいのだが、一番シンプルな2次項で回帰することを考える。また、2次関数の性質から、最大・最小値がx = 0まわりに存在しないのであれば、1次の項を加えて軸をずらしたほうが良いだろう。以下のように解析が可能である。なお、多項式回帰であればlm関数を使って解析できる。非線形の項で紹介したものの、多項式回帰は線形回帰とみなしてもよい。

------------------------------------------------------

fit1 <- lm(y ~ x1 + I(x1^2), data = d) #多項式回帰、2次以上の項はI()で囲む必要がある。もし、そうしないと単回帰になる

fit2 <- summary(fit1)

fit2

## 

## Call:

## lm(formula = y ~ x1 + I(x1^2), data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -50.943 -13.128  -0.644  12.167  38.737 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  1.27278    2.35769    0.54    0.591    

## x1          -5.48166    0.37538  -14.60   <2e-16 ***

## I(x1^2)      0.18593    0.01585   11.73   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 18.08 on 97 degrees of freedom

## Multiple R-squared:  0.701,  Adjusted R-squared:  0.6948 

## F-statistic: 113.7 on 2 and 97 DF,  p-value: < 2.2e-16

------------------------------------------------------

切片は1.27、1次項の係数は-5.48、2次項は0.19と推定された。二次関数の性質を思い出せば、軸の位置は-(-5.48/(0.19 * 2)) = 14.42付近であることがわかる。うまく推定されていそうな雰囲気がある。では、回帰曲線をグラフに描き込んでみよう。以下のようにできる。

------------------------------------------------------

coefs <- coef(fit2)

xx <- seq(min(d$x1), max(d$x1), length = 200)

yy <- coefs[1] + coefs[2]*xx + coefs[3]*(xx)^2

plotn(y ~ x1, data = d)

overdraw(points(xx, yy, col = "red", type = "l")) #図2の描画

------------------------------------------------------

するとデータをうまく表現できた回帰曲線が推定できている。ちなみに、今後は曲線をRで描くことも増えるが、Rには滑らかな曲線を描くようなプログラムは一部の例外を除き存在しない。そこで、基本的には推定された曲線のxとyの値を細かく求めて、その座標を直線で結ぶことで曲線を表現する。上のスクリプトではxxとyyのように求めたものが該当する。xxを求めたときのseq関数のlength引数が、分割数の指定に該当し、小さい値にすれば粗い「曲線」、大きな値にすれば滑らかな「曲線」が描ける。

　ネタ晴らしをすれば上記のデータは以下のようにして生成した。

------------------------------------------------------

b1 <- 4

b2 <- -6

b3 <- 0.2

n <- 100

x1 <- runif(n, -10, 30)

y <- rnorm(n = n, mean = b1 + b2*x1 + b3*(x1)^2, sd = 20)

plotn(y ~ x1)#図3の描画

------------------------------------------------------

おおむねいい推定値になっていることがわかるだろう。

・指数回帰

　以下のように、モデルが説明変数の指数関数で表される場合、指数回帰exponential regressionと呼ぶ。eはネイピア数である。b1 = 0として解析されることも多いかもしれない。例えば、環境による制限がないときの時間経過に対する個体数の増加は指数関数になることが知られている。

b2は指数部の正負に影響する。b3は指数関数の急峻さを示すと考えればよいだろう。また、xが正か負の無限大になるとe部分は0に収束し、その時の収束値がb1である。e以外の底で考えてもよいが、以下のように変換公式を使ってeの式に帰着可能である。eは数学的に便利な性質をたくさん有するので、eを使うことが多い。

またxに足し引きする項ついていてもよいが、この場合も上記の式に帰着できる。

　では具体的な解析例を紹介する。以下のようなデータを解析してみよう。

------------------------------------------------------

x1 <- c(4, 6.3, 1.5, 24.4, 13.7, 9, -6.5, 19.1, -1.1, 12.3, 2.8, 

        -9.9, 14.8, -5.1, 10, 29, 3.9, 8.3, 26.8, 28.4, 2.6, 9.5,

        -6.7, 19.2, 16.9, 22.7, 15.6, 6.1, 1.1, 19.4, 21.9, 9, -9.5,

        9.4, 1.7, -3.9, 0.1, 15, 16.8, -6.2, 9.4, -3.9, 12.2, 2, 

        19.6, 29.3, 11.6, -8.5, 1, 22.6, 4.9, 13.2, 15.6, 5, 3.3,

        1.1, -3.7, 24.6, 8, -2.5, -5, 18.9, -5.9, 12.4, -0.4, 23.3,

        11, 10.2, 2.9, -7.4, -8.5, 27.2, 18.5, 7.5, 25, 20, 1.5, -7.1,

        23.7, 17.6, 0.6, 18.6, 8, 24.4, 14.5, -5.8, 21.9, 19.5, 17.9, 

        21.5, -8.3, 2.6, 23.4, 1.5, 17.2, -1.7, -9, 9.2, 2.3, 24.1)

y <- c(12.8, 16.3, 10.2, 45.1, 1.9, -0.4, 9.8, 37.6, 17.2, 11.1, 13.2,

       12.5, 19.7, 11, 11.4, 150, 9.1, 4.8, 95.2, 118.4, 11.7, 15, 22.7, 

       35, 14.2, 46.2, 21.7, -8.2, 13.9, 12.6, 41.7, -20.1, 13.9, 29.5,

       9, -14.2, 20.3, 22.6, 12.9, -2, 2.9, 6.3, 14.5, -1.9, 19.5, 

       155.3, 16.2, -11.1, 8.7, 59.7, 14.4, 21.7, 17.5, 4.3, 9.6, -1.8,

       0.4, 50.5, 5.7, 15.8, -5.8, 9, 4.3, -2.7, 11, 44.6, 0.1, 8, 5.6,

       20.8, 13.3, 104.6, 32.8, 0.6, 67.9, 11.8, 2.9, 5.8, 66.1, 25, 8.4,

       10.9, -16.1, 43.8, 10.3, 8, 21.1, 32.7, 3.9, 40.7, 9.8, 7.8, 46.5,

       5.9, 18.4, -11.9, 21.1, 1.1, 22.7, 53)

d <- data.frame(x1 = x1, y = y)

head(d)

##     x1    y

## 1  4.0 12.8

## 2  6.3 16.3

## 3  1.5 10.2

## 4 24.4 45.1

## 5 13.7  1.9

## 6  9.0 -0.4

plotn(y ~ x1, data = d)#図4の描画

------------------------------------------------------

データはx1が大きくなると急速に大きな値になり、一方で小さくなれば何らかの収束値を持つような、曲線関係にあるようだ。このような関係の時は指数関数を疑うべきだ。ではこのようなモデルを持つ可能性があるときは、今までのlm関数では解析できない。その代わりnls関数を使い、モデルを明示して解析をする。

　同時にパラメータ探索時の初期値を設定する必要がある。モデルが複雑になったため、単回帰のように解析的（単純な式変形だけで）に最小二乗推定値を計算できない。そこで、nls関数はモデルとデータの間の偏差平方和をパラメータの関数として考え、この平方和が小さくなる方向に向かってパラメータを更新していくことで最終的なパラメータを探索している（数値的解）。大まかな発想としては、単回帰の時と同じで、平方和の関数をパラメータで偏微分して解析している（デフォルトの方法はガウス・ニュートン法、追記）。例えば、xが負の方向に行った時の収束値は0付近なので、b1の初期値は1、b2もb3も正の値でかつ極端に大きな値ではなさそうだから0.1を初期値としてみよう。これらの初期値はstart引数にlist形式で与える。

------------------------------------------------------

fit1 <- nls(y ~ b1 + b2*exp(b3*x1), data = d,

            start = list(b1 = 1, b2 = 0.1, b3 = 0.1))

#非線形回帰、dにx1とyが存在するのでこれらは変数、残りの文字はパラメータとして扱われる

#startにlist形式で各パラメータの初期値を与える

fit2 <- summary(fit1)

fit2

## 

## Formula: y ~ b1 + b2 * exp(b3 * x1)

## 

## Parameters:

##    Estimate Std. Error t value Pr(>|t|)    

## b1  6.71717    1.31880   5.093 1.73e-06 ***

## b2  0.23184    0.08944   2.592    0.011 *  

## b3  0.22043    0.01380  15.973  < 2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 9.866 on 97 degrees of freedom

## 

## Number of iterations to convergence: 9 

## Achieved convergence tolerance: 4.358e-06

------------------------------------------------------

上記のように推定するべきパラメータをb1～b3として、それぞれの推定値が求まっている。b1 = 6.71、b2 = 0.23、b3 = 0.22程度だ。非線形回帰を行う上での注意点として、初期値の値には細心の注意が必要である。というのも、特にガウス・ニュートン法では、平方差の最小（極小）値を与える正解パラメータ付近に初期値を設定しないとうまく収束しないことがある。例えば、グラフの形からb2やb3を負の値からスタートさせる意味はない。なので、データの様子からある程度の決め打ちが必要である。また、指数関数を含む時、その係数やべき指数のstartが0としてはならない。係数やべき指数が0以外なら指数関数であるが、0になると定数関数となり、関数の形が変わってしまう。すると平方和関数は0で微分不可能となり、ここではパラメータ更新できない。回帰曲線をグラフに描き込んでみよう。以下のようにできる。

------------------------------------------------------

coefs <- coef(fit2)

xx <- seq(min(d$x1), max(d$x1), length = 200)

yy <- coefs[1] + coefs[2]*exp(coefs[3]*xx) 

plotn(y ~ x1, data = d)#図5の描画

overdraw(points(xx, yy, col = "red", type = "l"))

------------------------------------------------------

するとデータをうまく表現できた回帰曲線が推定できている。ネタ晴らしをすれば上記のデータは以下のようにして生成した。

------------------------------------------------------

b1 <- 5

b2 <- 0.4

b3 <- 0.2

n <- 100

x1 <- runif(n, -10, 30)

y <- rnorm(n = n, mean = b1 + b2*exp(b3*x1), sd = 10)

plotn(y ~ x1)#図6の描画

------------------------------------------------------

おおむねいい推定値になっていることがわかるだろう。

・ロジスティック回帰

　以下のように、モデルが説明変数のロジスティック関数で表される場合、ロジスティック回帰logistic regressionと呼ぶ。eはネイピア数である。b2はこの関数の最小値、b3は最大値を規定する。b1は最大値と最小値の平均値での傾きを表し、b1 > 0なら右肩下がりで、xが大きくなれば最小値、小さくなれば最大値に近づく。b1 < 0なら右肩上がりになり、関係も逆になる。b4はこの関数が最大値、最小値の平均値に到達するxの値を規定する。この関数は座標（b4、(b3 - b2)/2）で点対称である。

例えば、環境による制限があるときの時間経過に対する個体数の増加はロジスティック関数になることが知られている。他にも、説明変数を除草剤濃度とし、除草剤処理しなかったときを生育量100%した時の応答もロジスティック関数のようになる。

　では具体的な解析例を紹介する。以下のようなデータを解析してみよう。

------------------------------------------------------

x1 <- c(24.4, 9.7, -8, 15.1, -1.7, -9.9, 22.7, 18.3, 16.8, 

        5.9, -9.5, 7.6, -1.5, -2.4, -5.8, -3.2, 14.8, 6.6, 

        5.7, 25.8, 20, 22.7, 26.1, 3.5, -6.8, -8, 27, 3.4, 

        4.2, 23.4, -4.5, -1.4, -5.3, -0.8, 1.5, 14.8, 4.7, 

        -7.4, 3.4, 6.5, 14.6, 18.3, 24.2, 14.7, 5.7, 9.5,

        19.3, -2.5, -1.2, -5.5, 15.9, -3.4, -2.6, -4.1, -4.6, 

        6.7, 9, 25.1, -7.1, -5.3, 9.6, 19.7, 7.1, 29.8, -6,

        16.4, 26.5, -4.1, -3.9, 1.9, 15, 9.4, 14.1, -9.3, 

        -2.3, -5.1, 15.9, -9.1, 15.3, 18.8, -2.2, -8.5,

        -7.5, 19.9, 19.8, 8.4, 4.2, 7.6, 29.7, 29.8, -7.9, 

        -4.7, 23.6, 20.1, 17.5, 12.3, -7.3, 9.2, -8.2, -3.7)

y <- c(14.7, 35.5, 82.6, 2.1, 98, 101.1, -4.4, -14.3, 5, 66.3,

       88.6, 52.5, 105.2, 84.7, 90.2, 101.9, 9.5, 73.7, 73.5, 

       -3.6, 10.4, 10.6, -3.1, 72.6, 109.4, 87.4, -8.5, 88.8, 

       86.9, 2.8, 107.9, 82.3, 79.6, 111.8, 96.1, -5.2, 86.5, 

       95, 75.4, 76.7, 16, -5.3, 0.2, 7.3, 58.4, 25.1, -5.1, 

       95.1, 91.7, 78.9, -8.5, 107.2, 95.2, 105.7, 86.1, 80.2,

       42.9, -11.1, 87.1, 94.4, 33, -2.6, 58.8, 10.7, 84.6, 

       15.2, -8.3, 88.8, 83.4, 99.3, 16.5, 45.8, 12.5, 75.8,

       104.3, 109.4, 14.8, 87.9, 10.7, -3.3, 96.4, 95.4, 93.8, 

       8.9, -4.7, 36.4, 90.9, 53.2, 15.1, -5.3, 87.5, 100.7, 

       6.6, 6.8, -4.3, 8.6, 101.7, 16.5, 111, 110.4)

d <- data.frame(x1 = x1, y = y)

head(d)

##     x1     y

## 1 24.4  14.7

## 2  9.7  35.5

## 3 -8.0  82.6

## 4 15.1   2.1

## 5 -1.7  98.0

## 6 -9.9 101.1

plotn(y ~ x1, data = d)#図7の描画

------------------------------------------------------

データはx1が大きいときや小さいときに極限として最大、最小値をもっている。こんな関係の時はロジスティック回帰を行うのが良いだろう。指数回帰の時と同様に、モデル式を明示して、初期値はstart引数にlist形式で与える。右肩下がりなので、b1には正の値でそこまで大きくない値を与える。b2は最小値を規定し、0付近なのでb2 = 0としてよいだろう。b3は最大値を規定し、100付近であるので、b3 = 100としてみる。b4はこのデータの関係において、最大値と最小値の平均になるときのxの値であるから、b4 = 5付近で探索させよう。

------------------------------------------------------

fit1 <- nls(y ~ b2 + (b3 - b2)/(1 + exp(b1*(x1 - b4))), data = d, 

            start = list(b1 = 0.1, b2 = 0, b3 = 100, b4 = 5)) #ロジスティック回帰

fit2 <- summary(fit1)

fit2

## 

## Formula: y ~ b2 + (b3 - b2)/(1 + exp(b1 * (x1 - b4)))

## 

## Parameters:

##    Estimate Std. Error t value Pr(>|t|)    

## b1  0.50585    0.06621   7.640 1.63e-11 ***

## b2  1.55082    1.73772   0.892    0.374    

## b3 95.05795    1.53366  61.981  < 2e-16 ***

## b4  8.14071    0.25489  31.938  < 2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 9.421 on 96 degrees of freedom

## 

## Number of iterations to convergence: 8 

## Achieved convergence tolerance: 2.452e-06

------------------------------------------------------

上記のように推定するべきパラメータをb1～b4として、それぞれの推定値が求まっている。b1 = 0.51、b2 = 1.55、b3 = 95.06、b4 = 8.14程度だ。指数回帰の時と同様、不適切な初期値はパラメータの収束を遅くさせたり、不可能にさせてしまう。特に、指数部分のb1は負にする意味はないだろう。回帰曲線をグラフに描き込んでみよう。以下のようにできる。

------------------------------------------------------

coefs <- coef(fit2)

xx <- seq(min(d$x1), max(d$x1), length = 200)

yy <- coefs[2] + (coefs[3] - coefs[2])/(1 + exp(coefs[1]*(xx - coefs[4])))

plotn(y ~ x1, data = d) #図8の描画

overdraw(points(xx, yy, col = "red", type = "l"))

------------------------------------------------------

するとデータをうまく表現できた回帰曲線が推定できている。ネタ晴らしをすれば上記のデータは以下のようにして生成した。

------------------------------------------------------

b1 <- 0.5

b2 <- 2

b3 <- 95

b4 <- 8

n <- 100

x1 <- runif(n, -10, 30)

y <- rnorm(n = n, mean = b2 + (b3 - b2)/(1 + exp(b1*(x1 - b4))), sd = 10)

plotn(y ~ x1)#図9の描画

------------------------------------------------------

おおむねいい推定値になっていることがわかるだろう。上記のような方法のほか、Rのdrcパッケージを使えば、比較的簡単にロジスティック回帰を実装できる。

おわりに

　上記では代表的な非線形回帰を紹介した。上記の他にももっと複雑なモデルを組むこともでき、これらの活用によって解析の幅も広がるだろう。ただし、上記の再掲であるが、複雑なモデルは解釈が複雑である点、複雑なモデルほどパラメータの収束が難しい点、特に初期値設定に気を付けなければならない点などは常に念頭に置いておこう。また、上記のように初期依存性が高く、解析が難しい場合でも、時として線形モデルの別の拡張である一般化線形モデルの導入によって無理なく解析できる場合も多い。様々な選択肢をとれるように、解析道具を取り揃えていくとよいだろう。

追記：パラメータ最適値を求めるためのガウス・ニュートン法

　Rの非線形回帰において、パラメータを数値的に解くためにデフォルトで用いられているのはガウス・ニュートン法Gauss–Newton methodと呼ばれる繰り返し計算の方法である。その方法の詳細な計算には立ち入らないが、コンセプトだけは知っておいても損がない。特に、初期値がなぜ重要なのかを理解するにも役立つだろう。

　以下のように残差平方和関数R(b1, b2……)があったとする。f(xi)はパラメータb1、b2……で表される非線形モデル式である。

目標はこれまでの最小二乗法と同じく、パラメータの関数と見たときに、残差平方和関数の最小となるパラメータを求めることである。線形回帰では、f(xi)が変数の足し算だけで表されていたので、解析的に残差平方和を最小にするパラメータを求めることができた。しかし、非線形モデルでは解析的にパラメータを求められないので、近似計算などを利用しながら繰り返し計算して、最小値を実現するパラメータに近づけていく。

　Rの最小値を与えるパラメータ候補は、今までと同じくパラメータで偏微分した時の、偏微分係数∂R/∂b1 = ∂R/∂b2 =…… = 0となる点である。では、この偏微分係数 = 0となる点をどのように求めるかといえば、例えば、b1だけに注目すれば図形的に以下のような手順を踏めばよいことがわかる。

偏微分後の関数∂R/∂b1が0と交わる点付近に初期値b1startをとる。この点における関数∂R/∂b1の接線を求める。この接線が0と交わる点を次のパラメータ候補b1new1とする。次はこの点における関数∂R/∂b1の接線を求めて……という風に繰り返し計算をすれば、このb1は徐々にRを最小にする推定値^b1に近づくだろう。更新が進まなくなった時、それを推定値とすればよい。この方法はニュートン法Newton methodと呼ばれる。上記の計算からわかるように関数∂R/∂b1の接線を求めるということは、この関数の微分、すなわちRの2階微分が必要になる。1変数であれば難易度が低いが、多変数になると2階微分係数を行列として並べたもの、ヘッセ行列Hessian matrixが必要になり（この行列の行列式をヘシアンHessianと呼ぶ、Rで回帰分析するときは、この用語を覚えておこう。エラーの内容を読み解くのに必要になることがある）、さらにこのヘッセ行列の逆行列の計算が必要になる。これらの計算は非常に重たいので、ヘッセ行列を近似に基づいて計算を軽くしたものがガウス・ニュートン法である。

　上記の議論を見れば、適切な初期値を指定することがいかに大切なことかわかるだろう。例えば、Rの最小値の2つ右隣の極小値近辺にb1の初期値を設定したらどうなるだろう。上記の論法だけなら、このb1はRの最小値を実現するパラメータには近づかず、極小値に囚われてしまう。このような解を局所最適解local optimumと呼ぶ。対して、関数全体で見たときの最適値は大域最適解global optimumと呼ぶ。もちろん、実際の解析ではこのような局所最適解から逃れるための方法が実装されているのだが、もし、局所最適解付近の「谷」が深ければ、その谷から脱出できないこともある。ゆえに、初期値はなるべく、想定されうるパラメータ付近に設定するべきなのである。