Edited at

R言語 - 回帰モデルを使った分析

More than 3 years have passed since last update.


はじめに

 この Tips では、回帰分析を紹介しています。サンプルの分析では、ランニングのトレーニング目標を考えるために、過去のトレーニング量に対して、回帰分析としていくつかのモデル式を適用し、モデルの良さを評価した上で、今後のランニング目標を決定しています。モデルは、累乗、指数、漸近指数、ロジスティック曲線、ゴンペルツ曲線を当てはめて、モデルの評価には、AIC(赤池情報量規準)を使用しました。

モデル名
モデル式

累乗モデル
$y=ax^b$

指数モデル
$y=ab^x$

漸近指数モデル
$y=Asymp+(R0-Asymp)exp(-exp(lrc)*x)$

ロジスティック曲線モデル
$y=Asym/(1+exp((xmid-x)/scal))$

ゴンペルツ曲線モデル
$y=Asym*exp(-b2*b3^x)$


データの作成


サンプルデータ

# 過去のトレーニング期間のデータを作成

timeFrame <- c( "Sep-12", "Oct-12", "Nov-12", "Dec-12", "Jan-13", "Feb-13", "Mar-13", "Apr-13", "May-13", "Jun-13", "Jul-13", "Aug-13", "Sep-13", "Oct-13", "Nov-13", "Dec-13", "Jan-14", "Feb-14", "Mar-14", "Apr-14", "May-14", "Jun-14", "Jul-14", "Aug-14", "Sep-14", "Oct-14", "Nov-14", "Dec-14", "Jan-15", "Feb-15", "Mar-15", "Apr-15", "May-15", "Jun-15", "Jul-15", "Aug-15", "Sep-15", "Oct-15", "Nov-15", "Dec-15" )
# 過去のトレーニング期間の回数データを作成
trainingSequence <- as.numeric(c(1:length(timeFrame)))
# 過去のトレーニング量(走行距離(Km))のデータを作成
runDistance <- as.numeric(c( "65", "33", "31.4", "55.6", "47", "77", "100.3", "50.35", "77.6", "80", "81.5", "53", "126.5", "85.1", "96", "128.8", "75.7", "80.6", "91.1", "83.1", "51", "63.5", "64.2", "63.1", "74.8", "86", "136", "116", "125", "105", "105", "93", "95.5", "85", "88", "100", "105.5", "89", "92", "117" ))
# トレーニングログをデータフレームとして作成
trainingLog <- data.frame( trainingSequence, timeFrame, runDistance )


累乗モデルへの当てはめ


R

# nlsファンクションを使って、累乗モデルに当てはめる

powerModel <- nls( trainingLog$runDistance~a*trainingLog$trainingSequence^b, start=list( a=1, b=1 ))
# 分析結果の確認
summary ( powerModel )


分析結果例

Formula: trainingLog$runDistance ~ a * trainingLog$trainingSequence^b

Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
a 44.22670 7.97657 5.545 2.4e-06 ***
b 0.22765 0.05831 3.904 0.000375 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 20.92 on 38 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: 6.662e-06

 累乗モデルからパラメータaとbの推定値を得られたので、次の式を推定することができました。


累乗モデル

トレーニング量=44.22670*トレーニング回数^{0.22765}



指数モデルへの当てはめ


R

# nlsファンクションを使って、指数モデルに当てはめる

exponentModel <- nls( trainingLog$runDistance~a*b^trainingLog$trainingSequence, start=list( a=1, b=1 ))
# 分析結果の確認
summary ( exponentModel )


分析結果例

Formula: trainingLog$runDistance ~ a * b^trainingLog$trainingSequence

Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
a 63.278770 6.040147 10.48 9.22e-13 ***
b 1.013546 0.003649 277.72 < 2e-16 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 21.82 on 38 degrees of freedom
Number of iterations to convergence: 8
Achieved convergence tolerance: 2.917e-06

 指数モデルからパラメータaとbの推定値を得られたので、次の式を推定することができました。


指数モデル

トレーニング量=63.27877*1.013546^{トレーニング回数}



漸近指数モデルへの当てはめ


R

# nlsファンクションを使って、漸近指数モデルに当てはめる

asymptoteModel <- nls( trainingLog$runDistance~SSasymp( trainingLog$trainingSequence, Asymp, R0, lrc ))
# 分析結果の確認
summary ( asymptoteModel )


分析結果例

Formula: trainingLog$runDistance ~ SSasymp(trainingLog$trainingSequence, Asymp, R0, lrc)

Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
Asymp 99.1614 8.8335 11.226 1.79e-13 ***
R0 38.1248 16.4724 2.314 0.02630 *
lrc -2.3442 0.6049 -3.875 0.00042 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 21.11 on 37 degrees of freedom
Number of iterations to convergence: 0
Achieved convergence tolerance: 2.524e-06

 漸近指数モデルからパラメータ Asymp、R0、lrc を得られたので、次の式を推定することができました。


漸近指数モデル

トレーニング量=99.1614+(38.1248-99.1614)exp(-exp(-2.3442)トレーニング回数



ロジスティック曲線モデルへの当てはめ


R

# nlsファンクションを使って、ロジスティック曲線モデルを当てはめる

logisticsModel <- nls( trainingLog$runDistance ~ SSlogis(trainingLog$trainingSequence, Asym, xmid, scal ))
# 分析結果の確認
summary( logisticsModel )


分析結果例

Formula: trainingLog$runDistance ~ SSlogis(trainingLog$trainingSequence, Asym, xmid, scal)

Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
Asym 95.608 5.708 16.749 <2e-16 ***
xmid 1.854 2.396 0.774 0.4438
scal 5.774 2.923 1.976 0.0557 .
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 21.18 on 37 degrees of freedom
Number of iterations to convergence: 2
Achieved convergence tolerance: 8.752e-06

 ロジスティック曲線モデルからパラメータ Asym、xmid、scal を得られたので、次の式を推定することができました。


ロジスティック曲線

トレーニング量=95.608/( 1+exp(( 1.854 - トレーニング回数 )/5.774))



ゴンペルツ曲線モデルへの当てはめ


R

# nlsファンクションを使って、ゴンペルツ曲線モデルを当てはめる

GompertzModel <- nls( trainingLog$runDistance~SSgompertz( trainingLog$trainingSequence, Asym, b2, b3 ))
# 分析結果の確認
summary ( GompertzModel )


分析結果例

Formula: trainingLog$runDistance ~ SSgompertz(trainingLog$trainingSequence, Asym, b2, b3)

Parameters: # パラメータ推定値、標準誤差、t 値、p 値
Estimate Std. Error t value Pr(>|t|)
Asym 97.15226 6.86999 14.142 <2e-16 ***
b2 0.89681 0.33065 2.712 0.0101 *
b3 0.87785 0.06143 14.290 <2e-16 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 21.14 on 37 degrees of freedom
Number of iterations to convergence: 3
Achieved convergence tolerance: 4.49e-06

 ゴンペルツ曲線モデルからパラメータ Asym、b2、b3 を得られたので、次の式を推定することができました。


ゴンペルツ曲線モデル

トレーニング量=97.15226*exp(-0.89681*0.87785^{トレーニング回数})



グラフを確認する

 式を推定できたので、トレーニングログの散布図上に、それぞれの式をグラフ描画します。


R

# トレーニングログを描画

par( bg="darkgray" , col="white", col.axis="white", col.lab="white", col.main="white", col.sub="white")
plot (trainingLog$trainingSequence, trainingLog$runDistance, pch=19, axes = FALSE, main="Training Log", ylab="Distance(Km)", xlab="Date")
axis(side=1, at = trainingLog$trainingSequence, labels = trainingLog$timeFrame)
axis(side=2)
par( new=TRUE )
# 累乗モデルへの当てはめを描画
lines(trainingLog$trainingSequence, 44.22670*trainingLog$trainingSequence^(0.22765), col = "red", lwd=4 )
# 指数モデルへの当てはめを描画
lines(trainingLog$trainingSequence, 63.278770 * 1.013546^trainingLog$trainingSequence, col = "yellow", lwd=4 )
# 漸近指数モデルへの当てはめを描画
lines(trainingLog$trainingSequence,( 99.1614+(38.1248-99.1614)*exp(-exp(-2.3442)*trainingLog$trainingSequence) ) , col = "blue", lwd=4 )
# ロジスティック曲線モデルへの当てはめを描画
lines(trainingLog$trainingSequence,( 95.608/(1 + exp(( 1.854 - trainingLog$trainingSequence )/5.774)) ) , col = "green", lwd=4 )
# ゴンペルツ曲線モデルへの当てはめを描画
lines(trainingLog$trainingSequence,( 97.15226*exp(-0.89681*0.87785^trainingLog$trainingSequence )) , col = "pink", lwd=4)


AICでモデルを評価する

 最後に、AIC(赤池情報量規準)で統計モデルの良さを評価します。具体的には AIC ファンクションを実行して、相対的に値が最小のモデルを選択します。


R

# 累乗モデル(赤色のグラフ)

AIC(powerModel)
# 指数モデル(黄色のグラフ)
AIC(exponentModel)
# 漸近指数モデル(青色のグラフ)
AIC(asymptoteModel)
# ロジスティック成長モデル(緑色のグラフ)
AIC(logisticsModel)
# ゴンペルツ成長モデル(桃色のグラフ)
AIC(GompertzModel)


実行結果例

> # 累乗モデル(赤色のグラフ)

> AIC(powerModel)
[1] 360.7268
> # 指数モデル(黄色のグラフ)
> AIC(exponentModel)
[1] 364.08
> # 漸近指数モデル(青色のグラフ)
> AIC(asymptoteModel)
[1] 362.393
> # ロジスティック曲線モデル(緑色のグラフ)
> AIC(logisticsModel)
[1] 362.6532
> # ゴンペルツ曲線モデル(桃色のグラフ)
> AIC(GompertzModel)
[1] 362.5037


さいごに

 サンプルの分析では、グラフの確認とAICの評価から、累乗モデル(赤色のグラフ)による回帰モデルを選択して、今後のランニングのトレーニング目標を計画することにしました。具体的には次のような計画になりますが、年間の走行距離が約1272Kmという無理のない現実的なトレーニング目標になったと思います。


リファレンス

R 基本統計関数マニュアル

福岡大学 非線形回帰分析