単回帰

最終更新: 2017-03-15 00:45

単回帰

data.csv

X,Y
15,325
20,365
25,412
30,478
35,538
40,575
45,659
50,712
55,820
60,893

回帰分析を行う

Rのプロンプトに以下を入力する。

d <- read.csv('data.csv')
d.lm <- lm(Y ~ X, data=d)
summary(d.lm)

以下のような分析結果が出力される。

Call:
lm(formula = Y ~ X, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-34.145 -11.895  -6.809  18.491  32.291 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 106.0182    21.0606   5.034  0.00101 ** 
X            12.5782     0.5245  23.983 9.74e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 23.82 on 8 degrees of freedom
Multiple R-squared:  0.9863,	Adjusted R-squared:  0.9846 
F-statistic: 575.2 on 1 and 8 DF,  p-value: 9.737e-09

回帰直線を引く

上のコードに続けて以下を入力する。

plot(d)
abline(d.lm)

以下のような図が出力される(略)

信頼区間を表示する

ggplotを用いて信頼区間をプロットする。

library(ggplot2)

X_new <- data.frame(X=15:60)
conf_95 <- predict(d.lm, X_new, interval='confidence', level=0.95)
conf_95 <- data.frame(X_new, conf_95)
conf_50 <- predict(d.lm, X_new, interval='confidence', level=0.50)
conf_50 <- data.frame(X_new, conf_50)

p <- ggplot()
p <- p + theme_bw()
p <- p + geom_point(data=d, aes(x=X, y=Y), shape=1, size=3)
p <- p + geom_line(data=conf_50, aes(x=X, y=fit), size=1)
p <- p + geom_ribbon(data=conf_50, aes(x=X, ymin=lwr, ymax=upr), alpha=2/6)
p <- p + geom_ribbon(data=conf_95, aes(x=X, ymin=lwr, ymax=upr), alpha=1/6)
print(p)
# ggsave(p, file='confidence.png', dpi=300)

予測区間を表示する

同様に予測区間の表示も行う。

pred_95 <- predict(d.lm, X_new, interval='prediction', level=0.95)
pred_95 <- data.frame(X_new, pred_95)
pred_50 <- predict(d.lm, X_new, interval='prediction', level=0.50)
pred_50 <- data.frame(X_new, pred_50)

p <- ggplot()
p <- p + theme_bw()
p <- p + geom_point(data=d, aes(x=X, y=Y), shape=1, size=3)
p <- p + geom_line(data=pred_50, aes(x=X, y=fit), size=1)
p <- p + geom_ribbon(data=pred_50, aes(x=X, ymin=lwr, ymax=upr), alpha=2/6)
p <- p + geom_ribbon(data=pred_95, aes(x=X, ymin=lwr, ymax=upr), alpha=1/6)
print(p)
# ggsave(p, file='prediction.png', dpi=300)