車の燃費データを使って回帰分析する。

1 データ

単位はオリジナル（出典参照）から馴染みのあるものに変換した。

出典：【UCI Machine Learning Repository】Auto MPG Data Set

変数名 内容（単位）
km 燃費（km/L）
ncy 気筒数（本）
cc 排気量（cc）
hp 馬力（hp）
kg 車体重量（kg）
sec 加速性能：時速97km（60mile）に到達する時間（秒）
yr　 製造年（年）
type アメ車，日本車，欧州車の3区分
name 車名（英語） 
d <- read.csv('https://stats.dip.jp/01_ds/data/car_mileage.csv')
rownames(d) <- paste0('No.', 1:nrow(d))
colnames(d) <- c('km', 'ncy', 'cc', 'hp', 'kg', 'sec', 'yr', 'type', 'name')

library(DT)
datatable(d)

2 訓練データ・テストデータ

n <- nrow(d)

# 訓練データサイズの全サイズ80％とした。sample関数で非復元無作為抽出する。
ii.tr <- sample(1:n, size = floor(0.8*n))

d.tr <- d[ ii.tr, ] # 訓練データを抽出
d.te <- d[-ii.tr, ] # 試験データを抽出（負のインデックスの要素は除かれる）

3 フィッティング・予測

fit <- lm(km ~ kg + cc, data = d.tr) # 訓練データでフィッティング

kmhat <- predict(fit, newdata = d.te) # 試験データの説明変数を用いて予測

RMSE <- sqrt(mean((d.te$km -kmhat)^2)) # 予測精度指標（２乗平均平方根誤差）
RMSE
## [1] 2.008881

3.1 精度評価

get.accuracy <- function(yhat, y, digits = 2)
{
  d <- data.frame(MBE  = mean(yhat - y),
                  MAE  = mean(abs(yhat - y)),
                  MAPE = mean(abs((yhat - y) / y)) * 100,
                  RMSE = sqrt(mean((yhat - y)^2)))
  return(round(d, digits))
}

(a <- get.accuracy(d.te$km, kmhat))

4 信頼区間・予測区間

4.1 フィッティング

2次の多項式回帰モデルをlm関数を使用して訓練データにフィッティングさせる。

\[Y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \epsilon_i,\quad \epsilon_i \sim \mathbf{N}(0, \sigma^2)\]

fit <- lm(km ~ poly(kg, 2, raw = TRUE), data = d)
#summary(fit)
library(sjPlot)
tab_model(fit)
  km
Predictors Estimates CI p
(Intercept) 26.47 23.96 – 28.97 <0.001
kg [1st degree] -0.02 -0.02 – -0.01 <0.001
kg [2nd degree] 0.00 0.00 – 0.00 <0.001
Observations 392
R2 / R2 adjusted 0.715 / 0.714

4.2 予測

# 予測に使用する入力データ（車体重量）
kg.p <- seq(0, 3000, 100)

# 信頼区間
conf <- predict(fit, newdata = data.frame(kg = kg.p),
                interval = 'confidence')
# 予測区間
pred <- predict(fit, newdata = data.frame(kg = kg.p),
                interval = 'prediction')

4.3 グラフ

# カラーパレット
COL <- c(rgb(255,   0,   0,  105, max = 255), # 赤
         rgb(  0,   0, 255,  105, max = 255), # 青
         rgb(  0, 155,   0,  105, max = 255), # 緑
         rgb(140, 140, 140,  105, max = 255), # 暗灰
         rgb(180, 180, 180,  105, max = 255)) # 灰

[RGB_Color] https://www.rapidtables.com/web/color/RGB_Color.html

#cairo_pdf(file = 'auto-mpg_graph.pdf')
# 散布図
matplot (x = d$kg, y = d$km,
         type = 'p', pch = 16, col = COL[1], ylim = c(0, 20),
         main = '２次多項式回帰モデルを使った燃費予測', 
         xlab = '車体重量［kg］', ylab = '燃費［km］')

# ポリゴン描画（薄墨色の領域を描く）
gray.area <- function(x, lwr, upr, col)
{
  polygon(c(x, rev(x)), c(lwr, rev(upr)), col = col, border = NA)
}
gray.area(kg.p, conf[, 'lwr'], conf[, 'upr'], col = COL[4]) # 信頼区間
gray.area(kg.p, pred[, 'lwr'], pred[, 'upr'], col = COL[5]) # 予測区間

# 予測値曲線を描画
matlines(x = kg.p, y = conf[, 'fit'], col = COL[2], lwd = 3)

# 凡例（はんれい）
legend('topright',
       pch = c(16, NA, NA, NA), col = COL[-3],
       lty = c(NA, 1, NA, NA), lwd = c(NA, 3, NA, NA), 
       fill = c(NA, NA, COL[4], COL[5]), border = F,
       legend = c('燃費', '予測値', '95％信頼区間', '95％予測区間'))

#dev.off()

5 演習課題

12時の東京の気温（tp）と関東の電力需要（gw）を使用し回帰分析せよ。 信頼区間と予測区間も描画せよ。

説明変数 内容
date 日
hr 時
gw 電力需要［GW］ギガワット
tp 気温［℃ ］ 
d <- read.csv('https://stats.dip.jp/01_ds/data/load/load_tp_y2022_h12_tokyo.csv')

datatable(d)