R による回帰分析

単回帰分析では独立変数 x を一つだけ用いて、従属変数 y を説明する分析手法である。次のように数式で表すことができる。

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

数式中の β0 と β1 を求めることで、y を x で説明するというものである。ここでは、R を利用して、単回帰分析の方法を示す。サンプルデータとして 1987 年から 2010 年の二酸化炭素濃度のデータを利用する。(ryourico2.csv

# データの読み込み
x <- read.csv("https://stat.biopapyrus.net/data/ryourico2.csv", header = T)

# 読み込んだデータを確認
head(x)
##   year co2ppm
## 1 1987  351.2
## 2 1988  354.0
## 3 1989  355.9
## 4 1990  356.7
## 5 1991  358.1
## 6 1992  358.4

読み込んだデータに対して単回帰分析を行う。lm 関数を利用する。また、year を独立変数とし、co2ppm を従属変数変数として、年が経つ毎に二酸化炭素の濃度がどのような振る舞いをするのかを分析する。

result <- lm(co2ppm ~ year, data = x)

co2ppm ~ year と与えることで、R は co2ppm = β0 + β1 year」と解釈して解析を進める。また、co2ppmyear はデータ x の列名であるから、data = x として与える。

解析結果の確認は summary 関数を利用する。

summary(result)
## Call:
## lm(formula = co2ppm ~ year, data = x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7920 -0.9308 -0.1964  0.8649  2.0894 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.190e+03  7.186e+01  -44.39   <2e-16 ***
## year         1.782e+00  3.596e-02   49.55   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
## 
## Residual standard error: 1.219 on 22 degrees of freedom
## Multiple R-squared: 0.9911,	Adjusted R-squared: 0.9907 
## F-statistic:  2455 on 1 and 22 DF,  p-value: < 2.2e-16 

(Intercept) の Estimate 値が β0 の値であり、year の Estimate 値が β1 の値である。すなわち、回帰式は「co2ppm = 1.782*year - 3.190」のようになっている。

また、Pr(>|t|) の列について、これはそれぞれ「β0 = 0」、「β1 = 0」を検定した p-value である。p-value がゼロに近づけば、β0 ≠ 0、β1 ≠ 0 であることを意味する。つまり、p-value がゼロに近づけば、求められた β0 は β1 は回帰式で意味のある係数、寄与率の高い係数として扱うことができる。逆に言えば、p-value が 1 に近づけば、それに該当する β0 または β1 を省略しても結果があまり影響を受けない。

解析結果をグラフにする時、plot 関数および abline 関数を用いる。

plot(x)           # すべてのデータを散布図で描く
abline(result)    # 回帰分析の結果を読み込んで回帰直線を描く
Rのlm関数で回帰分析を行った例

また、predict 関数を利用することで信頼区間や予測区間をプロットすることできる。

# 信頼区間を求める
conf <- predict(result, interval="confidence")
head(conf)
##        fit      lwr      upr
## 1 350.2470 349.2460 351.2480
## 2 352.0288 351.0909 352.9667
## 3 353.8106 352.9340 354.6872
## 4 355.5923 354.7749 356.4098
## 5 357.3741 356.6130 358.1353
## 6 359.1559 358.4478 359.8640


# 予測区間を求める
pred <- predict(result, interval="prediction")
head(pred)
##        fit      lwr      upr
## 1 350.2470 347.5271 352.9669
## 2 352.0288 349.3315 354.7261
## 3 353.8106 351.1340 356.4871
## 4 355.5923 352.9345 358.2502
## 5 357.3741 354.7331 360.0151
## 6 359.1559 356.5297 361.7821


# プロット
plot(x)           #散布図
abline(result)    #回帰直線

# 信頼区間(赤・実線)
lines(x[, 1], conf[, 2], col = "red", lty = 1)  #下界
lines(x[, 1], conf[, 3], col = "red", lty = 1)  #上界

# 予測区間(青・ダッシュ線)
lines(x[, 1], pred[, 2], col = "blue", lty = 2) #下界
lines(x[, 1], pred[, 3], col = "blue", lty = 2) #上界
Rのlm関数で回帰分析を行った例