Rで非線形回帰分析＆ggplot2で描画：対数近似の例（R）


  1. Part 1：データの整形（tidyr::gathar()）
  2. Part 2：非線形回帰分析（nls()＆glm()で対数近似・パラメータ推定）
  3. Part 3：ggplot2で描画（複数のデータセットを重ね描き）


Part 2：非線形回帰分析（nls()＆glm()で対数近似・パラメータ推定）

まずは平均値のデータでカーブフィッティングを行います。
ggplot2では非線形近似曲線を描いてくれるので、グラフの描画には実際には必要ありません。
しかし被験者ごと、あるいは群平均ごとに主観的等価点（PSE）や丁度可知差異（JND）を算出する場合には、モデルを作ってパラメータ（係数）の推定を行う必要があります。

非線形回帰分析には、nls()が使えます。
nlsは非線形最小二乗法（nonlinear least squere）でパラメータ推定を行います。
未知のパラメータa, bを求める際、初期値（start）を設定する必要があります。
Rによる非線形最小二乗法

fit <- nls(rate ~ a * log(duration) + b, start = c(a = 1, b = 1), data = d2)
> fit
Nonlinear regression model
  model: rate ~ a * log(duration) + b
   data: d2
     a      b 
0.5136 0.3853 
 residual sum-of-squares: 0.5382

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.706e-08

今回の例では、a = 0.5136、b = 0.3853と出力されます。
値を取り出す場合は、

a <- summary(fit)$coefficients[2]
b <- summary(fit)$coefficients[1]

同じように、一般化線形モデルでもパラメータ推定を行うことができます。

fit2 < -glm(rate ~ log(duration), family = gaussian, data=d2)
> fit

Call:  glm(formula = rate ~ log(duration), data = d2)

Coefficients:
  (Intercept)  log(duration)  
       0.3853         0.5136  

Degrees of Freedom: 13 Total (i.e. Null);  12 Residual
Null Deviance:	    4.524 
Residual Deviance: 0.5382 	AIC: 0.11

a <- fit2$coefficients[2]
b <- fit2$coefficients[1] # fit2$coefficients["(Intercept)"] でもOK。

なお、MatlabのCurve fitting toolboxのfittypeでも同じ結果を得ることができます。
Matlab ドキュメンテーション
これでパラメータ推定ができました。
Part 3へ続く。

↑ PAGE TOP