Last updated: October 2009. Kajiyama                    [ 目次に戻る ]

非線形回帰分析 A nonlinear regression analysis


  1. 非線形回帰分析の関数 nls( ) と非線形曲線の関数     nls=Nonlinear Least Squaresの略

    nls( 非線形モデル式, データフレーム名, start=初期値のリスト, control, algorithm, trace, subset, weights, na.action, model, lower, upper, ...)

    非線形回帰分析は単回帰分析と同じく,2変数の散布の曲線モデル Y = f( X ) の係数を推定する.

    Y に目的変数の測定値(データ)を与え,X に説明変数の測定値(データ) を与えたときに,2変数の差を最小とする最小二乗法で係数 a と b,c,・・・を推定する.

    ただしY = f( X ) の関数が線形ではないモデル式を使用する.  [非線形回帰分析のための関数リスト]

    累乗モデル Y= a Xb   [Y ~ a * X ^ b]  ただしx≠0.
    べき乗モデル
    
    指数モデル Y= a bX   [Y ~ a * b ^ X]
    
    漸近(ぜんきん)指数モデル Y= a bX+ c    [Y ~ a * b ^ X + c] 関数 ? SSasymp参照
    
    ロジスティク成長モデル Y= a /(1+b exp(-cX)) [Y ~ a /(1+b*exp(-c*X))] 関数 ? SSlogis参照
    
    ゴンペルツ成長モデル Y= a bexp(-c X)     [Y ~ a*b^exp(-c*X)]  関数 ? SSgompertz参照
    

    累乗モデル Y= a Xb の予測式の当てはめ

    指数モデル Y= a bX の予測式の当てはめ

    漸近指数モデル Y= a bX+ c の予測式の当てはめ

    ゴンペルツ成長モデル Y= a bexp(-c X)  の予測式の当てはめ

    予測式モデルの評価 赤池情報量規準 AIC


    参考 非線形関数とグラフ(1)     非線形関数とグラフ(2)     確率・統計用語(和英) 関数名が統一されていないので注意



  2. 「2007年度都道府県別総人口と順位」のデータの分布


    2007年度の47都道府県別の総人口(万人単位)と順位のデータをデータフレームにする.

    データの出典:政府統計局都道府県の指標 e-State

    ・「コンソール」 画面に,

    都道府県名 <- c( "北海道", "青森県", "岩手県", "宮城県", "秋田県", "山形県", "福島県", "茨城県", "栃木県", "群馬県", "埼玉県",
    "千葉県", "東京都", "神奈川県", "新潟県", "富山県", "石川県", "福井県", "山梨県", "長野県", "岐阜県", "静岡県", "愛知県", "三重県",
    "滋賀県", "京都府", "大阪府", "兵庫県", "奈良県", "和歌山県", "鳥取県", "島根県", "岡山県", "広島県", "山口県", "徳島県",
    "香川県", "愛媛県", "高知県", "福岡県", "佐賀県", "長崎県", "熊本県", "大分県", "宮崎県", "鹿児島県", "沖縄県" )


    # と記述し都道府県名を入力

    順位 <- c( 8, 29, 32, 15, 37, 34, 18, 11, 20, 19, 5, 6, 1, 2, 14, 38, 35, 43 , 41, 16, 17, 10, 4, 22, 30, 13, 3, 7, 28,
    39, 47, 46, 21, 12, 25, 44, 40, 27, 45, 9, 42, 26, 23, 33, 36, 24, 31 )


    # と記述し総人口の順位を入力

    総人口万 <- c( 557, 141, 136, 235, 112, 120, 207, 297, 201, 202, 709, 610, 1276, 888, 240, 111, 117, 82, 88, 218, 210, 380,
    736, 188, 140, 264, 881, 559, 141, 102, 60, 73, 195, 287, 147, 80, 101, 145, 78, 506, 86, 145, 183, 120, 114, 173, 137 )


    # と記述し総人口(万人単位)を入力

    # 次に,

    x <- data.frame( 都道府県名, 順位, 総人口万 )

    # と記述し3つの変数を結合して都道府県順の並びのデータフレームを作成

    x1 <- order( x$順位 ) # と記述し,順位変数の人口の順番で並べ替え

    人口分布2007 <- x[x1,] # と記述し,総人口の順番で並べ替えしたデータフレームを作成

    # 次に,

    plot (人口分布2007$順位, 人口分布2007$総人口万 ) # と記述し順位と総人口の散布図を作成,分布が曲線であることを確認する.


    > 都道府県名 <- c( "北海道", "青森県", "岩手県", "宮城県", "秋田県", "山形県", "福島県", "茨城県", "栃木県", "群馬県", "埼玉県", 
    "千葉県", "東京都", "神奈川県", "新潟県", "富山県", "石川県", "福井県", "山梨県", "長野県", "岐阜県", "静岡県", "愛知県", "三重県", 
    "滋賀県", "京都府", "大阪府", "兵庫県", "奈良県", "和歌山県", "鳥取県", "島根県", "岡山県", "広島県", "山口県", "徳島県", 
    "香川県", "愛媛県", "高知県", "福岡県", "佐賀県", "長崎県", "熊本県", "大分県", "宮崎県", "鹿児島県", "沖縄県" ) 
    
    > 順位 <- c( 8, 29, 32, 15, 37, 34, 18, 11, 20, 19, 5, 6, 1, 2, 14, 38, 35, 43 , 41, 16, 17, 10,  4, 22, 30, 13, 3, 7, 28,
    39, 47, 46, 21, 12, 25, 44, 40, 27, 45, 9, 42, 26, 23, 33, 36, 24, 31 ) 
    
    > 総人口万 <- c( 557, 141, 136, 235, 112, 120, 207, 297, 201, 202, 709, 610, 1276, 888, 240, 111, 117, 82, 88, 218, 210, 380, 
    736, 188, 140, 264, 881, 559, 141, 102, 60, 73, 195, 287, 147, 80, 101, 145, 78, 506, 86, 145, 183, 120, 114, 173, 137 ) 
    
    > x <- data.frame( 都道府県名, 順位, 総人口万 ) 
    
    > x1 <- order( x$順位 )    
    
    > 人口分布2007 <- x[x1,]    
    
    > str( 人口分布2007 ) 
    'data.frame':   47 obs. of  3 variables:
     $ 都道府県名: Factor w/ 47 levels "愛知県","愛媛県",..: 37 26 31 1 17 30 45 46 43 28 ...
     $ 順位      : num  1 2 3 4 5 6 7 8 9 10 ...
     $ 総人口万  : num  1276 888 881 736 709 ...
    
    > plot (人口分布2007$順位, 人口分布2007$総人口万 )
    

    図 人口分布2007



    結果: 散布図より,2つの変数は直線的な相関はないが,単調減少の相関関係を示している.

    この関係を非線形曲線の関数モデルで当てはめたい.

     ・累乗モデル Y= a Xb の予測式の当てはめ

     ・指数モデル Y= a bX の予測式の当てはめ

     ・漸近指数モデル Y= a bX+ c の予測式の当てはめ

     ・ゴンペルツ成長モデル Y= a bexp(-c X)  の予測式の当てはめ

     ・予測式モデルの評価 赤池情報量規準 AIC



  3. 「人口分布2007」への累乗モデル Y= a Xb 曲線の当てはめ


    (1) 非線形回帰分析

    累乗モデル Y= a Xb   Rでは Y ~ a * X ^ bを当てはめパラメータ a, b を推測する.

    ・「コンソール」 画面に,

    累乗モデル <- nls( 人口分布2007$総人口万 ~ a * 人口分布2007$順位 ^b, start=list( a=1, b=1 ) )

    # と記述し,非線形単回帰分析を行う

    summary ( 累乗モデル )

    # と記述し分析結果を表示する.

    > 累乗モデル <- nls( 人口分布2007$総人口万 ~ a * 人口分布2007$順位 ^b, start=list(a=1, b=1 ) ) 
    
    > summary( 累乗モデル ) 
    
    
    Formula: 人口分布2007$総人口万 ~ a * 人口分布2007$順位^b
    
    Parameters:
         Estimate    Std. Error   t value   Pr(>|t|)    
    a  1436.78886    62.00431     23.17    <2e-16 ***
    b    -0.62162     0.02446    -25.41    <2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    Residual standard error: 74 on 45 degrees of freedom
    
    Number of iterations to convergence: 14 
    Achieved convergence tolerance: 4.256e-06 
    

    結果

    累乗モデル Y= a Xb の係数を求めると,係数 a = 1436.78886 ,係数 b = -0.62162 である.

    よって累乗モデル Y= 1437 X-0.622 となる.


    (2) 散布図と累乗モデルの作図

    ・Rの「コンソール」画面に,

    plot (人口分布2007$順位, 人口分布2007$総人口万, main="Y=1437*X^(-0.622)" ) # と記述し測定値の散布図を作成

    par ( new=T ) # と記述しグラフの重ね図を指示

    lines(人口分布2007$順位, 1437*人口分布2007$順位^(-0.622), col = "blue" ) # と記述し,回帰曲線 Y= 1436 X-0.62 を作図


    > plot (人口分布2007$順位, 人口分布2007$総人口万, main="Y=1437*X^(-0.622)" )
    
    > par ( new=T )
    
    > lines(人口分布2007$順位, 1437*人口分布2007$順位^(-0.622), col = "blue" )
    


    (3) 結果と考察 ?

    青色は,回帰曲線 Y= a Xb =1437 X-0.622である.係数のP値はいずれも0.1%より小さく有意である.




  4. 「人口分布2007」への指数モデル Y= a bX 曲線の当てはめ


    (1) 非線形回帰分析

    累乗モデル Y= a bX   Rでは Y ~ a * b ^ Xを当てはめパラメータ a, b を推測する.

    ・「コンソール」 画面に,

    指数モデル <- nls(人口分布2007$総人口万~a*b^人口分布2007$順位, start=list( a=1, b=1 ) )

    # と記述し,非線形単回帰分析を行う

    summary ( 指数モデル )

    # と記述し分析結果を表示する.

    > 指数モデル <- nls(人口分布2007$総人口万~a*b^人口分布2007$順位, start=list( a=1, b=1 ) ) 
    
    > summary ( 指数モデル ) 
    
    Formula: 人口分布2007$総人口万 ~ a * b^人口分布2007$順位
    
    Parameters:
        Estimate   Std. Error t value  Pr(>|t|)    
    a  1.140e+03   4.579e+01   24.89   <2e-16 ***
    b  9.121e-01   4.566e-03  199.77   <2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    Residual standard error: 68.46 on 45 degrees of freedom
    
    Number of iterations to convergence: 12 
    Achieved convergence tolerance: 4.449e-06 
    

    結果

    指数モデル Y= a bX の係数を求めると,係数 a = 1.140e+03 ,係数 b = 9.121e-01 である.

    よって指数モデル Y= 1140*0.912 X となる.


    (2) 散布図と指数曲線の作図

    ・Rの「コンソール」画面に,

    plot (人口分布2007$順位, 人口分布2007$総人口万, main="Y=1140*0.912^X" ) # と記述し測定値の散布図を作成

    par ( new=T ) # と記述しグラフの重ね図を指示

    lines(人口分布2007$順位, 1140*0.912^人口分布2007$順位, col = "red" ) # と記述し指数モデル Y= 1140*0.912 X を作図

    lines(人口分布2007$順位, 1437*人口分布2007$順位^(-0.622), col = "blue" ) # と記述し累乗モデル Y= 1436 X-0.62 を作図

    > plot (人口分布2007$順位, 人口分布2007$総人口万, main="1140*0.912^X" )
    
    > par ( new=T )
    
    > lines(人口分布2007$順位, 1140*0.912^人口分布2007$順位, col = "red" )
    
    > lines(人口分布2007$順位, 1437*人口分布2007$順位^(-0.622), col = "blue" )
    


    (3) 結果と考察 ?

    赤色は,指数曲線 Y= a bX =1140*0.912^Xである.係数のP値はいずれも0.1%より小さく有意である.




  5. 「人口分布2007」への漸近(ぜんきん)指数モデル Y= a bX+ c 曲線の当てはめ


    (1) 非線形回帰分析

    漸近指数モデル a bX+ c   Rでは Y ~ (a * b ^ X + c) を当てはめパラメータ a,b,c を推測する.

    関数 SSasymp a=R0-Asym, b=exp(-exp(lrc)), c=Asymとし,
    
    Y ~ ( (R0-Asym) * exp(-exp(lrc)) ^ X + Asym) としてパラメータを自動的に推定する.
    
    書式 SSasymp(説明変数X, Asym, R0, lrc)
    


    ・「コンソール」 画面に,

    漸近指数モデル  <- nls( 人口分布2007$総人口万 ~ SSasymp( 人口分布2007$順位, Asymp, R0, lrc ) )

    # と記述し,非線形回帰分析を行う

    summary ( 漸近指数モデル )

    # と記述し分析結果を表示する.

    > 漸近指数モデル  <- nls( 人口分布2007$総人口万 ~ SSasymp( 人口分布2007$順位, Asymp, R0, lrc ) ) 
    
    > summary ( 漸近指数モデル ) 
    
    Formula: 人口分布2007$総人口万 ~ SSasymp(人口分布2007$順位, Asymp, R0, lrc)
    
    Parameters:
           Estimate    Std. Error t value  Pr(>|t|)    
    Asymp   105.138      8.552    12.29    7.92e-16 ***
    R0     1316.717     35.141    37.47     < 2e-16 ***
    lrc      -1.941      0.047   -41.29     < 2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    Residual standard error: 37.86 on 44 degrees of freedom
    
    Number of iterations to convergence: 0 
    Achieved convergence tolerance: 9.27e-07 
    

    結果

    漸近指数モデル Y= a bX+ c の係数a=R0-Asym, b=exp(-exp(lrc)), c=Asymを求めると,

    係数 Asym=105.138,R0=1316.717,lrc=-1.941 である.

    Asym <- 105.138
    R0 <- 1316.717
    lrc <- -1.941

    a=R0-Asym=1316.717-105.138=1211.579

    b=exp(-exp(lrc))= exp(-exp(-1.941))=0.8662685

    c=Asymp=105.138 である.

    よって漸近指数曲線モデル Y= 1212*0.866 X+ 105 となる.


    (2) 散布図と漸近指数モデルの作図

    ・Rの「コンソール」画面に,

    plot (人口分布2007$順位, 人口分布2007$総人口万, main="Y=1212*0.866^X+105" ) # と記述し測定値の散布図を作成

    par ( new=T ) # と記述しグラフの重ね図を指示

    lines(人口分布2007$順位, (1212*0.866^人口分布2007$順位+ 105), col = "green" ) # と記述し漸近指数モデルY= 1212*0.866 X+ 105 を作図

    lines(人口分布2007$順位, 1140*0.912^人口分布2007$順位, col = "red" ) # と記述し指数モデルY= 1140*0.912 X を作図

    lines(人口分布2007$順位, 1437*人口分布2007$順位^(-0.622), col = "blue" ) # と記述し累乗モデルY= 1436 X-0.62 を作図

    > plot (人口分布2007$順位, 人口分布2007$総人口万, main="1212*0.866^X+105" )
    
    > par ( new=T )
    
    > lines(人口分布2007$順位, (1212*0.866^人口分布2007$順位+ 105), col = "green" )  
    
    > lines(人口分布2007$順位, 1140*0.912^人口分布2007$順位, col = "red" )
    
    > lines(人口分布2007$順位, 1437*人口分布2007$順位^(-0.622), col = "blue" )
    


    (3) 結果と考察 ?

    緑色は,漸近指数モデル Y= a bX+ c =1212*0.866 X+ 105 である.係数のP値はいずれも0.1%より小さく有意である.





  6. 「人口分布2007」への ゴンペルツ成長モデル Y= a bexp(-c X) 曲線の当てはめ


    (1) 非線形回帰分析:ゴンペルツ成長モデル(Gompertz curve)

    ゴンペルツ成長モデル Y= a bexp(-c X)   Rでは Y ~ a*b^exp(-c*X) を当てはめパラメータ a,b,c を推測する.

    関数 SSgompertz a=Asym, b=exp(-b2), c=-log(b3)とし,
    
    Y~Asym*exp(-b2*b3^x) としてパラメータを自動的に推定する.
    
    書式 SSgompertz(説明変数X, Asym, b2, b3)
    


    ・「コンソール」 画面に,

    ゴンペルツモデル <- nls( 人口分布2007$総人口万 ~ SSgompertz( 人口分布2007$順位, Asym, b2, b3 ) )

    # と記述し,非線形回帰分析を行う

    summary ( ゴンペルツモデル )

    # と記述し分析結果を表示する.

    > ゴンペルツモデル <- nls( 人口分布2007$総人口万 ~ SSgompertz( 人口分布2007$順位, Asym, b2, b3 ) ) 
    
    > summary ( ゴンペルツモデル ) 
    
    Formula: 人口分布2007$総人口万 ~ SSgompertz(人口分布2007$順位, Asym, b2, 
        b3)
    
    Parameters:
           Estimate   Std. Error  t value   Pr(>|t|)    
    Asym  73.003319  10.305427    7.084   8.61e-09 ***
    b2    -2.954319   0.129888  -22.745    < 2e-16 ***
    b3     0.944866   0.004778  197.742    < 2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
    
    Residual standard error: 35 on 44 degrees of freedom
    
    Number of iterations to convergence: 0 
    Achieved convergence tolerance: 2.481e-07 
    

    結果

    ゴンペルツ成長モデル Y= Asym*exp(-b2*b3X) の係数 Asym,b2, b3を求めると,

    係数 Asym=73.003319,b2=-2.954319,b3=0.944866 である.

    よって関数 SSgompertz より Y= Asym*exp(-b2*b3X)= 73.0*exp(-(-2.95)*0.945X) ----(1) となる.

    元の式に戻すと,a=Asym=73.0,b=exp(-b2)=19.19,c=-log(b3)=0.0567 であるから,

    ゴンペルツ成長モデル Y= a bexp(-c X)= 73.0*19.19 exp(-(0.0567)*X) ------------(2) となる.

    式(1)(2)は同じである.


    (2) 散布図とゴンペルツ成長曲線の作図

    ・Rの「コンソール」画面に,

    plot (人口分布2007$順位, 人口分布2007$総人口万, main="Y=73.0* 19.19^exp(-(0.0567)*X" ) # と記述し測定値の散布図を作成

    par ( new=T ) # と記述しグラフの重ね図を指示

    lines(人口分布2007$順位, (73.0* 19.19^exp(-(0.0567)*人口分布2007$順位)), col = "yellow", lwd=2 ) # と記述し

    # ゴンペルツ曲線モデルY= 73.0*exp(-(-2.95)*0.945X) を作図


    lines(人口分布2007$順位, (1212*0.866^人口分布2007$順位+ 105), col = "green" )

    # と記述し漸近指数モデルY= 1212*0.866 X+ 105 を作図

    lines(人口分布2007$順位, 1140*0.912^人口分布2007$順位, col = "red" ) # と記述し指数モデルY= 1140*0.912 X を作図

    lines(人口分布2007$順位, 1437*人口分布2007$順位^(-0.622), col = "blue" ) # と記述し累乗モデルY= 1436 X-0.62 を作図

    > plot (人口分布2007$順位, 人口分布2007$総人口万, main="73.0*exp(-(-2.95)*0.945^x)" )  
    
    > par ( new=T )
    
    > lines(人口分布2007$順位, (73.0*exp(-(-2.95)*0.945^人口分布2007$順位)), col = "yellow", lwd=2 )  # ゴンペルツ成長モデル
    
    > lines(人口分布2007$順位, (1212*0.866^人口分布2007$順位+ 105), col = "green" ) # 漸近指数モデル 
    
    > lines(人口分布2007$順位, 1140*0.912^人口分布2007$順位, col = "red" ) # 指数モデル
    
    > lines(人口分布2007$順位, 1437*人口分布2007$順位^(-0.622), col = "blue" ) # 累乗モデル
    


    (3) 結果と考察 ?

    黄色は,ゴンペルツ成長モデル Y= a bexp(-c X)= 73.0*19.19 exp(-0.0567*X) である.係数のP値はいずれも0.1%より小さく有意である.



  7. 赤池情報量規準による「モデルの選択」


    赤池情報量規準は統計モデルの良さを評価するための指標で AIC とも呼ばれ,この指標が小さいほうが良いモデルと考えられる.

    ・「コンソール」 画面に,

    AIC( 累乗モデル ) # と記述し累乗モデルのAIC値を計算する.

    AIC( 指数モデル ) # と記述し指数モデルのAIC値を計算する.

    AIC( 漸近指数モデル ) # と記述し漸近指数モデルのAIC値を計算する.

    AIC( ゴンペルツモデル ) # と記述しゴンペルツ成長モデルのAIC値を計算する.

    > AIC( 累乗モデル ) 
    [1] 541.9217
    
    > AIC( 指数モデル ) 
    [1] 534.6032
    
    > AIC( 漸近指数モデル ) 
    [1] 479.8645
    
    > AIC( ゴンペルツモデル ) 
    [1] 472.4701
    

    結果

    4曲線と測定値の散布図から検討すると漸近指数モデルとゴンペルツ成長モデルがよく適合している.

    AICの値で確認すると,ゴンペルツモデルが小さい.

    ゴンペルツ成長モデル Y= 73.0*19.19 exp(-0.0567*X) を回帰式として選択する.




  8. 関数を調べるには,「コンソール」 画面に,? 関数名,コマンド名を入れる


    > ? SSasymp  # で関数SSasympを調べる.
    
    SSasymp(stats)      Asymptotic Regression Model
    
    This selfStart model evaluates the asymptotic regression function and its gradient. 
    It has an initial attribute that will evaluate initial estimates of 
    the parameters Asym, R0, and lrc for a given set of data. 
    
    Usage
    SSasymp(input, Asym, R0, lrc)
    
    Arguments
    input    a numeric vector of values at which to evaluate the model. 
    Asym     a numeric parameter representing the horizontal asymptote
             on the right side (very large values of input). 
    R0       a numeric parameter representing the response when input is zero. 
    lrc      a numeric parameter representing the natural logarithm of the rate constant. 
    
    Value
    a numeric vector of the same length as input. It is the value of the expression 
    
    Y=Asym+(R0-Asym)*exp(-exp(lrc)*input) 
    

非線形回帰分析                          [ 目次に戻る ]