原文鏈接：http://tecdat.cn/?p=9024
用GAM進行建模時間序列
我已經准備了一個文件，其中包含四個用電時間序列以進行分析。數據操作將由data.table程序包完成
將提及的智能電表數據讀到data.table。
DT <- as.data.table(read_feather("DT_4_ind"))
使用GAM回歸模型。將工作日的字符轉換為整數，並使用recode包中的函數car重新編碼工作日，以適應一周中出現的情況：1.星期一，…，7星期日。
DT[, week_num := as.integer(car::recode(week,
"'Monday'='1';'Tuesday'='2';'Wednesday'='3';'Thursday'='4';
'Friday'='5';'Saturday'='6';'Sunday'='7'"))]將信息存儲在日期變量中，以簡化工作。
n_type <- unique(DT[, type])
n_date <- unique(DT[, date])
n_weekdays <- unique(DT[, week])
period <- 48讓我們看一下用電量的一些數據並對其進行分析。
data_r <- DT[(type == n_type[1] & date %in% n_date[57:70])]
ggplot(data_r, aes(date_time, value)) +
geom_line() +
theme(panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.minor = element_line(colour = "grey90"),
panel.grid.major = element_line(colour = "grey90"),
panel.grid.major.x = element_line(colour = "grey90"),
axis.text = element_text(size = 10),
axis.title = element_text(size = 12, face = "bold")) +
labs(x = "Date", y = "Load (kW)")
在繪制的時間序列中可以看到兩個主要的季節性：每日和每周。我們在一天中有48個測量值，在一周中有7天，因此這將是我們用來對響應變量進行建模的自變量–電力負荷。
訓練我們的第一個GAM。通過平滑函數s對自變量建模，對於每日季節性，使用三次回歸樣條，對於每周季節性，使用P樣條。
gam_1 <- gam(Load ~ s(Daily, bs = "cr", k = period) +
s(Weekly, bs = "ps", k = 7),
data = matrix_gam,
family = gaussian)首先是其可視化功能。
layout(matrix(1:2, nrow = 1))
plot(gam_1, shade = TRUE)
我們在這里可以看到變量對電力負荷的影響。在左圖中，白天的負載峰值約為下午3點。在右邊的圖中，我們可以看到在周末消費量減少了。
讓我們使用summary函數對第一個模型進行診斷。
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Load ~ s(Daily, bs = "cr", k = period) + s(Weekly, bs = "ps",
## k = 7)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2731.67 18.88 144.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Daily) 10.159 12.688 119.8 <2e-16 ***
## s(Weekly) 5.311 5.758 130.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.772 Deviance explained = 77.7%
## GCV = 2.4554e+05 Scale est. = 2.3953e+05 n = 672EDF：估計的自由度–可以像對給定變量進行平滑處理那樣來解釋（較高的EDF值表示更復雜的樣條曲線）。P值：給定變量對響應變量的統計顯着性，通過F檢驗進行檢驗（越低越好）。\（R ^ 2 \）–調整后的R平方（越高越好）。我們可以看到R-sq。（adj）值有點低...
讓我們繪制擬合值：
我們需要將兩個自變量的相互作用包括到模型中。
第一種交互類型對兩個變量都使用了一個平滑函數。
gam_2 <- gam(Load ~ s(Daily, Weekly),
data = matrix_gam,
family = gaussian)
summary(gam_2)$r.sq## [1] 0.9352108R平方值表明結果要好得多。
summary(gam_2)$s.table## edf Ref.df F p-value
## s(Daily,Weekly) 28.7008 28.99423 334.4754 0似乎也很好，p值為0，這意味着自變量很重要。擬合值圖：
現在，讓我們嘗試上述張量積交互。這可以通過function完成te，也可以定義基本函數。
## [1] 0.9268452與以前的模型相似gam_2。
summary(gam_3)$s.table## edf Ref.df F p-value
## te(Daily,Weekly) 23.65709 23.98741 354.5856 0非常相似的結果。讓我們看一下擬合值：
與gam_2模型相比，只有一點點差異，看起來te更合身。
## [1] 0.9727604summary(gam_4)$sp.criterion## GCV.Cp
## 34839.46summary(gam_4)$s.table## edf Ref.df F p-value
## te(Daily,Weekly) 119.4117 149.6528 160.2065 0我們可以在這里看到R方略有上升。
讓我們繪制擬合值：
這似乎比gam_3模型好得多。
## [1] 0.965618summary(gam_4_fx)$s.table## edf Ref.df F p-value
## te(Daily,Weekly) 335 335 57.25389 5.289648e-199我們可以看到R平方比模型gam_4低，這是因為我們過度擬合了模型。證明GCV程序（lambda和EDF的估計）工作正常。
因此，讓我們在案例（模型）中嘗試ti方法。
## [1] 0.9717469summary(gam_5)$sp.criterion## GCV.Cp
## 35772.35summary(gam_5)$s.table## edf Ref.df F p-value
## s(Daily) 22.583649 27.964970 444.19962 0
## s(Weekly) 5.914531 5.995934 1014.72482 0
## ti(Daily,Weekly) 85.310314 110.828814 41.22288 0然后使用t2。
## [1] 0.9738273summary(gam_6)$sp.criterion## GCV.Cp
## 32230.68summary(gam_6)$s.table## edf Ref.df F p-value
## t2(Daily,Weekly) 98.12005 120.2345 86.70754 0我還打印了最后三個模型的GCV得分值，這也是在一組擬合模型中選擇最佳模型的良好標准。我們可以看到，對於t2相應模型gam_6，GCV值最低。
在統計中廣泛使用的其他模型選擇標准是AIC（Akaike信息准則）。讓我們看看三個模型：
AIC(gam_4, gam_5, gam_6)## df AIC
## gam_4 121.4117 8912.611
## gam_5 115.8085 8932.746
## gam_6 100.1200 8868.628最低值在gam_6模型中。讓我們再次查看擬合值。
我們可以看到的模型的擬合值gam_4和gam_6非常相似。可以使用軟件包的更多可視化和模型診斷功能來比較這兩個模型。
第一個是function gam.check，它繪制了四個圖：殘差的QQ圖，線性預測變量與殘差，殘差的直方圖以及擬合值與響應的關系圖。讓我們為它們制作模型gam_4和gam_6。
gam.check(gam_4)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 7 iterations.
## The RMS GCV score gradiant at convergence was 0.2833304 .
## The Hessian was positive definite.
## The estimated model rank was 336 (maximum possible: 336)
## Model rank = 336 / 336
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## te(Daily,Weekly) 335.00 119.41 1.22 1gam.check(gam_6)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 9 iterations.
## The RMS GCV score gradiant at convergence was 0.05208856 .
## The Hessian was positive definite.
## The estimated model rank was 336 (maximum possible: 336)
## Model rank = 336 / 336
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## t2(Daily,Weekly) 335.00 98.12 1.18 1我們可以再次看到模型非常相似，只是在直方圖中可以看到一些差異。
layout(matrix(1:2, nrow = 1))
plot(gam_4, rug = FALSE, se = FALSE, n2 = 80, main = "gam n.4 with te()")
plot(gam_6, rug = FALSE, se = FALSE, n2 = 80, main = "gam n.6 with t2()")
該模型gam_6 有更多的“波浪形”的輪廓。因此，這意味着它對響應變量的適應性更高，而平滑因子更低。
vis.gam(gam_6, n.grid = 50, theta = 35, phi = 32, zlab = "",
ticktype = "detailed", color = "topo", main = "t2(D, W)")
我們可以看到最高峰值是Daily變量的值接近30（下午3點），而Weekly變量的值是1（星期一）。
vis.gam(gam_6, main = "t2(D, W)", plot.type = "contour",
color = "terrain", contour.col = "black", lwd = 2)
再次可以看到，電力負荷的最高值是星期一的下午3:00，直到星期四都非常相似，然后負荷在減少（周末）。