回帰分析（予測）

d <- read.csv('https://stats.dip.jp/01_ds/data/load/load_tp_y2022_h12_tokyo.csv')
rownames(d) <- paste0('No.',1:nrow(d))
colnames(d) <- c('date','hr','gw','tp')

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)
gwhat <- predict(fit,newdata = d.te)

RMSE <- sqrt(mean((d.te$gw -gwhat)^2))
RMSE

## [1] 6.282425

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$gw,gwhat))

##    MBE  MAE  MAPE RMSE
## 1 0.24 5.27 14.46 6.28

fit <- lm(gw ~ poly(tp,2,raw = TRUE),data = d)
library(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,50,1)

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


matplot (x = d$tp,y = d$gw,type = 'p',pch = 16,col = COL[1],ylim = c(20,80),main = '２次多項式回帰モデルを使った需要予測',xlab = '気温［℃］',ylab = '電力需要［GW］')

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％予測区間'))