d <- read.csv(‘https://stats.dip.jp/01_ds/data/milk_sales.csv’) str(d)

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,]#試験データを抽出(負のインデックスの要素は除かれる)

モデル1

fit1 <- lm(Sales ~ Price, data = d.tr)

モデル2

fit2 <- lm(Sales ~ Price + Display + Feature, data = d.tr)

モデル3

fit3 <- lm(Sales ~ Price + Display + Feature + DisFeat, data = d.tr)

library(sjPlot)

tab_model(fit1, fit2, fit3, show.aic = T, dv.labels = c(‘Model1’, ‘Model2’, ‘Model3’))

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)) }

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

(a <- get.accuracy(d.te$Sales, yhat))

同時プロモーション時の過去実績データ

is.both <- d\(Display & d\)Feature & d\(DisFeat d.both <- d[is.both, ] #d.both <- d[d\)DisFeat == 1, ] # こちらでも良い

予測用入力データ

d.p <- data.frame(Price = seq(min(d.both\(Price), max(d.both\)Price), 1), Display = 1, Feature = 1, DisFeat = 1)

信頼区間

conf <- predict(fit3, newdata = d.p, interval = ‘confidence’)

予測区間

pred <- predict(fit3, newdata = d.p, interval = ‘prediction’)

散布図

matplot(x = d.both\(Price, y = d.both\)Sales, ylim = c(60, 140), type = ‘p’, pch = 16, col = COL[1], main = ‘特別陳列とチラシ同時プロモ時の売上予測’, xlab = ‘価格[円]’, ylab = ‘売上[千円]’)

ポリゴン描画(薄墨色の領域を描く)

gray.area <- function(x, lwr, upr, col) { polygon(c(x, rev(x)), c(lwr, rev(upr)), col = col, border = NA) } gray.area(d.p\(Price, conf[, 'lwr'], conf[, 'upr'], col = COL[4]) # 信頼区間 gray.area(d.p\)Price, pred[, ‘lwr’], pred[, ‘upr’], col = COL[5]) # 予測区間

予測値曲線を描画

matlines(x = d.p$Price, y = conf[, ‘fit’], col = COL[2], lwd = 3) grid()

凡例(はんれい)

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%予測区間’))