d <- read.csv('https://stats.dip.jp/01_ds/data/load/load_tp_y2022_h12_tokyo.csv')
library(DT)
datatable(d)

n <- nrow(d)
ii.tr <- sample(1:n, size = floor(0.8*n))
d.tr <- d[ ii.tr, ]
d.te <- d[-ii.tr, ]
fit <- lm(gw ~ tp, data = d.tr)
tp_hat <- predict(fit, newdata = d.te)
RMSE <- sqrt(mean((d.te$tp - tp_hat)^2))
RMSE

## [1] 17.97735

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(tp_hat, d.te$tp))

##     MBE   MAE   MAPE  RMSE
## 1 16.99 16.99 132.44 17.98

fit <- lm(gw ~ poly(tp, 2, raw = TRUE), data = d)
library(sjPlot)

## Learn more about sjPlot with 'browseVignettes("sjPlot")'.

tab_model(fit)

	gw
Predictors	Estimates	CI	p
(Intercept)	57.11	55.11 – 59.10	<0.001
tp [1st degree]	-2.85	-3.07 – -2.62	<0.001
tp [2nd degree]	0.08	0.07 – 0.08	<0.001
Observations	364
R2 / R2 adjusted	0.693 / 0.691

tp.p <- seq(0, max(d$gw), 5)
conf <- predict(fit, newdata = data.frame(tp = tp.p), 
                interval = 'confidence')
pred <- predict(fit, newdata = data.frame(tp = tp.p), 
                interval = 'prediction')
COL <- c(rgb(255,   30,  42,  105, max = 255), 
         rgb( 50,   10, 225,  105, max = 255), 
         rgb(  50, 155,  110,  105, max = 255), 
         rgb(140, 140, 40,  105, max = 255), 
         rgb(180, 180, 140,  105, max = 255)) 
matplot(x = d$tp, y = d$gw,
        type = 'p', pch = 16, col = COL[1], ylim = c(min(d$gw)-5, max(d$gw)),
        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(tp.p, conf[, 'lwr'], conf[, 'upr'], col = COL[4])
gray.area(tp.p, pred[, 'lwr'], pred[, 'upr'], col = COL[5])
matlines(x = tp.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％予測区間'))