読者です 読者をやめる 読者になる 読者になる

餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

Rのグラフに回帰分析の結果である数式を書き込む

R言語

データをプロットしたグラフに、推定された予測値や、推定されたパラメーターを数式ごと書き込みたいときはありませんか?
そう難しい話では無いのですが、幾つかコツがあるので例を作ってみました。

1. データセット

グラフを描くにはデータセットが要るので、今回は以下の数字を

# 奈良女子大学のページにあった実験観測データを拝借
t <- c(0,5,10,15,20,25,30) # 時間
Ts <- c(56.5, 32.1, 21.9, 17.0, 14.6, 11.9, 10.9) # 物体の温度
Tm <- 0 # 外気温(冷却材は氷でゼロ度だが、室温はもっと高い模様。なお、9.660583度を仮定すると、もっとも当てはまりが良い)

2. 推定を行う

推定を行い、予測値を計算してみましょう。モデルの制約によりpredict関数は使えません。

# ニュートンの冷却の法則を想定し、定数を推定する
dT <- (Ts - Tm)/(Ts[1]- Tm)
# OLSで推定。
r <- lm(log(dT) ~ t + 0)
# 予測値を計算
prdct <- Tm + (Ts[1] - Tm)*exp(coef(r)[1]*t)

3. プロットする座標の範囲を決める

グラフを重ねる前に縦横の座標の範囲を決めておかないと、ズレます*1

xlim <- c(0, max(t))
ylim <- c(0, 100)

4. データをプロットする

まずは観測値をプロットしてみます。

par(mar=c(5, 3.5, 1, 1), mgp=c(2, 1, 0))
plot(Ts ~ t, xlab="経過時間", ylab="", axes=FALSE, xlim=xlim, ylim=ylim)
x_at <- seq(xlim[1], xlim[2], 10)
y_at <- seq(ylim[1], ylim[2], 25)
x_labels <- sprintf("%d分", x_at)
y_labels <- sprintf("%d℃", y_at)
axis(1, at=x_at, labels=x_labels, pos=ylim[1])
axis(2, at=y_at, labels=y_labels, pos=xlim[1], las=2)

5. 推定値を重ね書きする

推定値をグラフに書き込みます。ここはRの紹介で良く触れられますね。

par(new=T)
plot(prdct ~ t, type="l", xlab="", ylab="", axes=FALSE, xlim=xlim, ylim=ylim, lty=2)

なお、こんな手間隙かけずにlines(t, prdct, lty=2)するべきだろと突っ込まれました。おっしゃる通り。
直線なのが気に入らない人は、plotの代わりに以下のようにcurveを使ってください。

par(new=T)
curve(Tm + (Ts[1] - Tm)*exp(coef(r)[1]*t), xname="t", xlab="", ylab="", axes=FALSE, xlim=xlim, ylim=ylim, lty=2)

6. 推定された式を凡例として表示する

まず、数式をsprintf関数で作って、parse関数で表現式に変えます。expression関数は、代入が手間なので使いません。使える数式の記号はhelp plotmathをすると一覧が見られます。

expr <- parse(text = sprintf("%.3f + (%.3f-%.3f)*e^{%1.3f*t}", Tm, Ts[1], Tm, coef(r)[1]))

次に、legend関数で凡例を表示します。推定値はlty、観測値はpchになりますが、使わない方には-1を入れていることに注意してください。

legend("topright", lty=c(2, -1), pch=c(-1, 1), legend=c(expr, "観測値"), bty="n")

7. データの出所を明記する

データを拝借しているので、出所を明記しておきましょう。

mtext("出所) 奈良女子大 - 「第11回の授業 (07/07/09)」の実験データから編集" , side=1, line=3.5, adj=-1)

8. 完成図

相関係数などを書いても良かったかも知れません。

*1:推定値のプロットにpar(new=T)とplot()を使うからズレるのであって、低水準関数のlines()を使えば座標の範囲を決めなくてもいいと言うツッコミがありましたが、ここを決めておかないと軸目盛などの位置もズレると言うか、あわせるのが困難になるので、やはり決めておいた方がいいです。