Rのlm関数で二酸化炭素濃度を増加勾配を分析する例

R による回帰分析

単回帰は、1 つの独立変数で 1 つ従属変数を説明したい場合に利用される分析手法である。そのため、単回帰は、独立変数 x と従属変数 y を用いて、次のように数式化することができる。

\[ y = \beta_{0} + \beta_{1}x \]

単回帰のモデル式中にあるパラメーター β0 と β1 を求めることで、1 次関数ができる。この 1 次関数を利用することで、観測データを正しく説明できるだけでなく、予測なども行えるようになる。

lm 関数を利用した単回帰分析

R では、lm 関数を利用することで、単回帰分析を行うことができる。ここで、R の anscombe データセットを利用して、lm 関数で単回帰してみる。anscombe データセットには、4 セットのデータが含まれている。

data(anscombe)
head(anscombe)
##   x1 x2 x3 x4   y1   y2    y3   y4
## 1 10 10 10  8 8.04 9.14  7.46 6.58
## 2  8  8  8  8 6.95 8.14  6.77 5.76
## 3 13 13 13  8 7.58 8.74 12.74 7.71
## 4  9  9  9  8 8.81 8.77  7.11 8.84
## 5 11 11 11  8 8.33 9.26  7.81 8.47
## 6 14 14 14  8 9.96 8.10  8.84 7.04

anscombe の 4 セットのデータのうち、まず 1 番目のセット(x1 と y1 列)に対して回帰分析をおこなってみる。

y <- anscombe$y1
x <- anscombe$x1

m <- lm(y ~ x)

summary(m)
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max
## -1.92127 -0.45577 -0.04136  0.70941  1.83882
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   3.0001     1.1247   2.667  0.02573 *
## x             0.5001     0.1179   4.241  0.00217 **
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared:  0.6665,	Adjusted R-squared:  0.6295
## F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217

回帰の結果には、回帰のモデル式に含まれているパラメーター β0 と β1 が計算される。それらの値は Coefficients 項の Estimate 列に表示されている。たとえば、この場合、切片(β0)が 3.0001 であり、傾き(β1)が 0.5001 となる。これにより、回帰直線のモデル式は、y = 0.5001x + 3.0001 と書くことができる。

Pr(>|t|) の列は、各パラメーターが 0 かどうかを検定した結果を p 値で示したものである。例えば、β1 の p 値が 1.0 のときは、β1 = 0 を意味する。このとき、モデル式は y = β0 となり、x と y の間に因果関係がない、と大まかに見積もることができる。逆に、β1 の p 値が 0.0 のときは、β1 ≠ 0 を意味する。そのため、モデル式は y = β0 + β1x となり、x と y の間に因果関係があると考えることができる。

解析結果を図にプロットするとき、abline 関数などを利用する。

plot(x, y, xlim = c(0, max(x)), ylim = c(0, max(y)))
abline(m)
R の lm 関数で回帰分析を行った結果

anscombe データセットには 4 セットのデータが含まれているので、この 4 セットのデータすべてに対して lm 関数で回帰分析をして、図示してみる。lm 関数のオプションの与え方は、上の方法とは異なるものの、実質には同じことが書かれている。

m1 <- lm(y1 ~ x1, data = anscombe)
m2 <- lm(y2 ~ x2, data = anscombe)
m3 <- lm(y3 ~ x3, data = anscombe)
m4 <- lm(y4 ~ x4, data = anscombe)


par(mfrow = c(2, 2))  # plot 4 figures
plot(anscombe$x1, anscombe$y1)
abline(m1)
plot(anscombe$x2, anscombe$y2)
abline(m2)
plot(anscombe$x3, anscombe$y3)
abline(m3)
plot(anscombe$x4, anscombe$y4)
abline(m4)

anscombe の 4 つのデータセットそれぞれに対して、回帰分析を行った結果は、次のように図示される。この図から見られるように、回帰分析は基本的に、全観測値からの誤差がもっとも小さくなるように回帰直線を計算しているだけである。データの分布によらず、回帰直線は必ず計算される。

R の lm 関数で anscombe の 4 つのデータセットに対して回帰分析を行った結果。

回帰直線が描けたらそれで終わりではない。anscombe データセットの場合、1 つ目のデータセット(x1, y1)は、回帰分析の結果はそのまま利用可能であるが、それ以外の場合は別の解析手法をとったほうがいい。例えば、2 つ目のデータセット(x2, y2)では、単回帰よりも 2 次曲線で近似曲線したほうがいいし、また、3 つ目と 4 つ目のデータセットでは、外れ値を取り除いてから単回帰したほうがいい。

信頼区間と予測区間

回帰分析で得られる回帰モデルに対して、信頼区間と予測区間を求めることができる。信頼区間は、母平均がある確率で入る区間のことをさす。これに対して、予測区間は、今後、新たに標本(データ)を取得するとき、その標本がある確率で入る区間のことをさす。

信頼区間

単回帰で計算された回帰直線から、信頼区間を取得するのに predict 関数を利用する。このとき、引数 interval="confidence" を指定する。デフォルトでは 95% 信頼区間が計算されるが、信頼区間の範囲を変更する場合は level 引数を調整する。信頼区間の計算結果は、3 列からなり、fit 列は x を回帰直線に代入して得られる値であり、lwr は信頼区間の下限で、upr は信頼区間の上限である。

y <- anscombe$y1
x <- anscombe$x1

m <- lm(y ~ x)
conf.interval <- predict(m, interval="confidence", level = 0.95)
head(conf.interval)
##         fit      lwr       upr
## 1  8.001000 7.116387  8.885613
## 2  7.000818 6.116205  7.885431
## 3  9.501273 8.141258 10.861287
## 4  7.500909 6.657464  8.344354
## 5  8.501091 7.503113  9.499069
## 6 10.001364 8.423422 11.579305

回帰分析の結果に、fit 列、lwr 列、および upr 列のデータを追加して図示する。

plot(x, y, xlim = c(0, max(x)), ylim = c(0, max(y)))
abline(m)

points(x, conf.interval[, 1], col = "orange")
lines(x[order(x)], conf.interval[order(x), 2], col = "blue")
lines(x[order(x)], conf.interval[order(x), 3], col = "darkgreen")
R の lm 関数で回帰分析を行った結果から求めた信頼区間。

予測区間

予測区間は、 predict 関数を利用して、引数 interval="prediction" を指定する。信頼区間の計算結果は、3 列からなり、fit 列は x を回帰直線に代入して得られる値であり、lwr は信頼区間の下限で、upr は信頼区間の上限である。

y <- anscombe$y1
x <- anscombe$x1

m <- lm(y ~ x)
conf.predict <- predict(m, interval="prediction", level = 0.95)
head(conf.predict)
##         fit      lwr       upr
## 1  8.001000 5.067072 10.934928
## 2  7.000818 4.066890  9.934747
## 3  9.501273 6.390801 12.611745
## 4  7.500909 4.579129 10.422689
## 5  8.501091 5.531014 11.471168
## 6 10.001364 6.789620 13.213107

回帰分析の結果に、fit 列、lwr 列、および upr 列のデータを追加して図示する。

plot(x, y, xlim = c(0, max(x)), ylim = c(0, max(y)))
abline(m)

points(x, conf.predict[, 1], col = "orange")
lines(x[order(x)], conf.interval[order(x), 2], col = "blue")
lines(x[order(x)], conf.interval[order(x), 3], col = "darkgreen")
R の lm 関数で回帰分析を行った結果から求めた予測区間。