拓端tecdat|【視頻】R語言廣義相加模型（GAM）在電力負荷預測中的應用
原文鏈接：http://tecdat.cn/?p=9024
原文出處：拓端數據部落公眾號
拓端tecdat：R語言廣義相加模型（GAM）在電力負荷預測中的應用
視頻：R語言廣義相加模型（GAM）在電力負荷預測中的應用
1導言
這篇文章探討了為什么使用廣義相加模型?是一個不錯的選擇。為此，我們首先需要看一下線性回歸，看看為什么在某些情況下它可能不是最佳選擇。
2回歸模型
假設我們有一些帶有兩個屬性Y和X的數據。如果它們是線性相關的，則它們可能看起來像這樣：
為了檢查這種關系，我們可以使用回歸模型。線性回歸是一種使用X來預測變量Y的方法。將其應用于我們的數據將預測成紅線的一組值：
這就是“直線方程式”。根據此等式，我們可以從直線在y軸上開始的位置（“截距”或α）開始描述，并且每個單位的x都增加了多少y（“斜率”），我們將它稱為x的系數，或稱為β）。還有一點自然的波動，如果沒有的話，所有的點都將是完美的。我們將此稱為“殘差”（?）。
數學上是：
或者，如果我們用實際數字代替，則會得到以下結果：
這篇文章通過考慮每個數據點和線之間的差異（“殘差）然后最小化這種差異來估算模型。
我們在線的上方和下方都有正誤差和負誤差，因此，通過對它們進行平方并最小化“平方和”，使它們對于估計都為正。這稱為“普通最小二乘法”或OLS。
3非線性關系如何？
因此，如果我們的數據看起來像這樣，我們該怎么辦：
我們剛剛看到的模型的關鍵假設之一是y和x線性相關。如果我們的y不是正態(tài)分布的，則使用廣義線性模型?_（Nelder＆Wedderburn，1972）_，其中y通過鏈接函數進行變換，但再次假設f（y）和x線性相關。如果不是這種情況，并且關系在x的范圍內變化，則可能不是最合適的。我們在這里有一些選擇：
我們可以使用線性擬合，但是如果這樣做的話，我們會在數據的某些部分上面或者下面。
我們可以分為幾類。我在下面的圖中使用了三個，這是一個合理的選擇。同樣，我們可能處于數據某些部分之下或之上，而在類別之間的邊界附近似乎是準確的。例如，如果x = 49時，與x = 50相比，y是否有很大不同？
我們可以使用多項式之類的變換。下面，我使用三次多項式，因此模型適合：
。這些的組合使函數可以光滑地近似變化。這是一個很好的選擇，但可能會極端波動，并可能在數據中引起相關性，從而降低擬合度。
4樣條曲線
多項式的進一步細化是擬合“分段”多項式，我們在數據范圍內將多項式鏈在一起以描述形狀?！皹訔l線”是分段多項式，以繪圖員用來繪制曲線的工具命名。物理樣條曲線是一種柔性條，可以彎曲成形，并由砝碼固定。在構造數學樣條曲線時，我們有多項式函數，二階導數連續(xù)，固定在“結”點上。
下面是一個ggplot2?對象，該?對象的?geom_smooth?的公式包含ns?函數中的“自然三次樣條”? 。這種樣條曲線為“三次”
，并且使用10個結
5光滑函數
樣條曲線可以是光滑的或“搖擺的”，這可以通過改變節(jié)點數（k）或使用光滑懲罰γ來控制。如果我們增加結的數目，它將更“搖擺”。這可能會更接近數據，而且誤差也會更小，但我們開始“過度擬合”關系，并擬合我們數據中的噪聲。當我們結合光滑懲罰時，我們會懲罰模型中的復雜度，這有助于減少過度擬合。
6廣義相加模型（GAM）
廣義加性模型（GAM）（Hastie，1984）使用光滑函數（如樣條曲線）作為回歸模型中的預測因子。
這些模型是嚴格可加的，這意味著我們不能像正常回歸那樣使用交互項，但是我們可以通過重新參數化作為一個更光滑的模型來實現(xiàn)同樣的效果。事實并非如此，但本質上，我們正轉向一種模型，如：
摘自Wood _（2017）_的GAM的更正式示例?是：
其中：
μi≡E（Yi），Y的期望
Yi?EF（μi，?i），Yi是一個響應變量，根據均值μi和形狀參數?的指數族分布。
Ai是任何嚴格參數化模型分量的模型矩陣的一行，其中θ為對應的參數向量。
fi是協(xié)變量xk的光滑函數，其中k是每個函數的基礎。
如果您要建立回歸模型，但懷疑光滑擬合會做得更好，那么GAM是一個不錯的選擇。它們適合于非線性或有噪聲的數據。
7 gam擬合
那么，如何?為上述S型數據建立 GAM模型？
在這里，我將使用三次樣條回歸?：
gam(Y?~?s(X,?bs="cr")
上面的設置意味著：
s()指定光滑器。還有其他選項，但是s是一個很好的默認選項bs=“cr”告訴它使用三次回歸樣條（'basis'）。s函數計算出要使用的默認結數，但是您可以將其更改為k=10，例如10個結。
8模型輸出：
查看模型摘要：
##?##?Family:?gaussian?##?Link?function:?identity?##?Parametric?coefficients:##?????????????Estimate?Std.?Error?t?value?Pr(>|t|)????##?(Intercept)??43.9659?????0.8305???52.94???<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(X)?6.087??7.143?296.3??<2e-16?***##?---##?Signif.?codes:??0?'***'?0.001?'**'?0.01?'*'?0.05?'.'?0.1?'?'?1##?##?R-sq.(adj)?=??0.876???Deviance?explained?=?87.9%##?GCV?=?211.94??Scale?est.?=?206.93????n?=?300
顯示了我們截距的模型系數，所有非光滑參數將在此處顯示
每個光滑項的總體含義如下。
這是基于“有效自由度”（edf）的，因為我們使用的樣條函數可以擴展為許多參數，但我們也在懲罰它們并減少它們的影響。
9檢查模型：
該?gam.check()?函數可用于查看殘差圖，但它也可以測試光滑器以查看是否有足夠的結來描述數據。但是如果p值很低，則需要更多的結。
##?##?Method:?GCV???Optimizer:?magic##?Smoothing?parameter?selection?converged?after?4?iterations.##?The?RMS?GCV?score?gradient?at?convergence?was?1.107369e-05?.##?The?Hessian?was?positive?definite.##?Model?rank?=??10?/?10?##?##?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##?s(X)?9.00?6.09?????1.1????0.97
10它比線性模型好嗎？
讓我們對比具有相同數據的普通線性回歸模型：
anova(my\_lm,?my\_gam)
##?Analysis?of?Variance?Table##?##?Model?1:?Y?~?X##?Model?2:?Y?~?s(X,?bs?=?"cr")##???Res.Df???RSS?????Df?Sum?of?Sq??????F????Pr(>F)????##?1?298.00?88154??????????????????????????????????????##?2?292.91?60613?5.0873?????27540?26.161?<?2.2e-16?***##?---##?Signif.?codes:??0?'***'?0.001?'**'?0.01?'*'?0.05?'.'?0.1?'?'?1
我們的方差分析函數在這里執(zhí)行了f檢驗，我們的GAM模型明顯優(yōu)于線性回歸。
11小結
所以，我們看了什么是回歸模型，我們是如何解釋一個變量y和另一個變量x的。其中一個基本假設是線性關系，但情況并非總是這樣。當關系在x的范圍內變化時，我們可以使用函數來改變這個形狀。一個很好的方法是在“結”點處將光滑曲線鏈接在一起，我們稱之為“樣條曲線”
我們可以在常規(guī)回歸中使用這些樣條曲線，但是如果我們在GAM的背景中使用它們，我們同時估計了回歸模型以及如何使我們的模型更光滑。
上面的示例顯示了基于樣條的GAM，其擬合度比線性回歸模型好得多。
12用GAM進行建模用電負荷時間序列
我已經準備了一個文件，其中包含四個用電時間序列來進行分析。數據操作將由data.table程序包完成。
將提及的智能電表數據讀到data.table。
DT?<-?as.data.table(read\_feather("ind"))
使用GAM回歸模型。將工作日的字符轉換為整數，并使用recode包中的函數重新編碼工作日：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)")
在繪制的時間序列中可以看到兩個主要的季節(jié)性：每日和每周。我們在一天中有48個測量值，在一周中有7天，因此這將是我們用來對因變量–電力負荷進行建模的自變量。
訓練我們的第一個GAM。通過平滑函數s對自變量建模，對于每日季節(jié)性，使用三次樣條回歸，對于每周季節(jié)性，使用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?=?672
EDF：估計的自由度–可以像對給定變量進行平滑處理那樣來解釋（較高的EDF值表示更復雜的樣條曲線）。P值：給定變量對因變量的統(tǒng)計顯著性，通過F檢驗進行檢驗（越低越好）。調整后的R平方（越高越好）。我們可以看到R-sq.（adj）值有點低。
讓我們繪制擬合值：
我們需要將兩個自變量的交互作用包括到模型中。
第一種交互類型對兩個變量都使用了一個平滑函數。
gam_2?<-?gam(Load?~?s(Daily,?Weekly),summary(gam_2)$r.sq
##?\[1\]?0.9352108
R方值表明結果要好得多。
summary(gam_2)$s.table
##?????????????????????edf???Ref.df????????F?p-value##?s(Daily,Weekly)?28.7008?28.99423?334.4754???????0
似乎也很好，p值為0，這意味著自變量很重要。擬合值圖：
現(xiàn)在，讓我們嘗試上述變量交互。這可以通過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.46
summary(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低，這是因為我們過度擬合了模型。證明lambda和EDF的估計工作正常。
因此，讓我們在案例（模型）中嘗試ti方法。
##?\[1\]?0.9717469summary(gam_5)$sp.criterion
##???GCV.Cp?##?35772.35
summary(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.68
summary(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值最低。
在統(tǒng)計中廣泛使用的其他模型選擇標準是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???????1
gam.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，直到星期四都非常相似，然后負荷在周末減少。
最受歡迎的見解
1.在python中使用lstm和pytorch進行時間序列預測
2.python中利用長短期記憶模型lstm進行時間序列預測分析
3.使用r語言進行時間序列（arima，指數平滑）分析
4.r語言多元copula-garch-模型時間序列預測
5.r語言copulas和金融時間序列案例
6.使用r語言隨機波動模型sv處理時間序列中的隨機波動
7.r語言時間序列tar閾值自回歸模型
8.r語言k-shape時間序列聚類方法對股票價格時間序列聚類
9.python3用arima模型進行時間序列預測
