原文連結：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不是正態分佈的，則使用廣義線性模型_（Nelder＆Wedderburn，1972）_，其中y通過連結函式進行變換，但再次假設f（y）和x線性相關。如果不是這種情況，並且關係在x的範圍內變化，則可能不是最合適的。我們在這裡有一些選擇：

• 我們可以使用線性擬合，但是如果這樣做的話，我們會在資料的某些部分上面或者下面。

• 我們可以分為幾類。我在下面的圖中使用了三個，這是一個合理的選擇。同樣，我們可能處於資料某些部分之下或之上，而在類別之間的邊界附近似乎是準確的。例如，如果x = 49時，與x = 50相比，y是否有很大不同？

• 我們可以使用多項式之類的變換。下面，我使用三次多項式，因此模型適合：。這些的組合使函式可以光滑地近似變化。這是一個很好的選擇，但可能會極端波動，並可能在資料中引起相關性，從而降低擬合度。

4樣條曲線

多項式的進一步細化是擬合“分段”多項式，我們在資料範圍內將多項式鏈在一起以描述形狀。“樣條線”是分段多項式，以繪圖員用來繪製曲線的工具命名。物理樣條曲線是一種柔性條，可以彎曲成形，並由砝碼固定。在構造數學樣條曲線時，我們有多項式函式，二階導數連續，固定在“結”點上。

下面是一個ggplot2物件，該物件的geom_smooth的公式包含ns函式中的“自然三次樣條” 。這種樣條曲線為“三次”

，並且使用10個結

5光滑函式

樣條曲線可以是光滑的或“搖擺的”，這可以通過改變節點數（k）或使用光滑懲罰γ來控制。如果我們增加結的數目，它將更“搖擺”。這可能會更接近資料，而且誤差也會更小，但我們開始“過度擬合”關係，並擬合我們資料中的噪聲。當我們結合光滑懲罰時，我們會懲罰模型中的複雜度，這有助於減少過度擬合。

6廣義相加模型（GAM）

廣義加性模型（GAM）（Hastie，1984）使用光滑函式（如樣條曲線）作為迴歸模型中的預測因子。

這些模型是嚴格可加的，這意味著我們不能像正常回歸那樣使用互動項，但是我們可以通過重新引數化作為一個更光滑的模型來實現同樣的效果。事實並非如此，但本質上，我們正轉向一種模型，如：

摘自Wood _（2017）_的GAM的更正式示例是：

其中：

• μi≡E（Yi），Y的期望

• Yi〜EF（μi，ϕi），Yi是一個響應變數，根據均值μi和形狀引數ϕ的指數族分佈。

• Ai是任何嚴格引數化模型分量的模型矩陣的一行，其中θ為對應的引數向量。

• fi是協變數xk的光滑函式，其中k是每個函式的基礎。

如果您要建立迴歸模型，但懷疑光滑擬合會做得更好，那麼GAM是一個不錯的選擇。它們適合於非線性或有噪聲的資料。

7 gam擬合

那麼，如何為上述S型資料建立 GAM模型？

在這裡，我將使用三次樣條迴歸：

gam(Y~s(X,bs="cr")

上面的設定意味著：

• s()指定光滑器。還有其他選項，但是s是一個很好的預設選項

• bs=“cr”告訴它使用三次迴歸樣條（'basis'）。

• s函式計算出要使用的預設結數，但是您可以將其更改為k=10，例如10個結。

8模型輸出：

檢視模型摘要：

1. ##
2. ##Family:gaussian
3. ##Linkfunction:identity
5. ##Parametriccoefficients:
6. ##EstimateStd.ErrortvaluePr(>|t|)
7. ##(Intercept)43.96590.830552.94<2e-16***
8. ##---
9. ##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1
10. ##
11. ##Approximatesignificanceofsmoothterms:
12. ##edfRef.dfFp-value
13. ##s(X)6.0877.143296.3<2e-16***
14. ##---
15. ##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1
16. ##
17. ##R-sq.(adj)=0.876Devianceexplained=87.9%
18. ##GCV=211.94Scaleest.=206.93n=300

• 顯示了我們截距的模型係數，所有非光滑引數將在此處顯示

• 每個光滑項的總體含義如下。

• 這是基於“有效自由度”（edf）的，因為我們使用的樣條函式可以擴充套件為許多引數，但我們也在懲罰它們並減少它們的影響。

9檢查模型：

該gam.check()函式可用於檢視殘差圖，但它也可以測試光滑器以檢視是否有足夠的結來描述資料。但是如果p值很低，則需要更多的結。

1. ##
2. ##Method:GCVOptimizer:magic
3. ##Smoothingparameterselectionconvergedafter4iterations.
4. ##TheRMSGCVscoregradientatconvergencewas1.107369e-05.
5. ##TheHessianwaspositivedefinite.
6. ##Modelrank=10/10
7. ##
8. ##Basisdimension(k)checkingresults.Lowp-value(k-index<1)may
9. ##indicatethatkistoolow,especiallyifedfisclosetok'.
10. ##
11. ##k'edfk-indexp-value
12. ##s(X)9.006.091.10.97

10它比線性模型好嗎？

讓我們對比具有相同資料的普通線性迴歸模型：

anova(my\_lm,my\_gam)

1. ##AnalysisofVarianceTable
2. ##
3. ##Model1:Y~X
4. ##Model2:Y~s(X,bs="cr")
5. ##Res.DfRSSDfSumofSqFPr(>F)
6. ##1298.0088154
7. ##2292.91606135.08732754026.161<2.2e-16***
8. ##---
9. ##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1

我們的方差分析函式在這裡執行了f檢驗，我們的GAM模型明顯優於線性迴歸。

11小結

所以，我們看了什麼是迴歸模型，我們是如何解釋一個變數y和另一個變數x的。其中一個基本假設是線性關係，但情況並非總是這樣。當關系在x的範圍內變化時，我們可以使用函式來改變這個形狀。一個很好的方法是在“結”點處將光滑曲線連結在一起，我們稱之為“樣條曲線”

我們可以在常規迴歸中使用這些樣條曲線，但是如果我們在GAM的背景中使用它們，我們同時估計了迴歸模型以及如何使我們的模型更光滑。

上面的示例顯示了基於樣條的GAM，其擬合度比線性迴歸模型好得多。

12用GAM進行建模用電負荷時間序列

我已經準備了一個檔案，其中包含四個用電時間序列來進行分析。資料操作將由data.table程式包完成。

將提及的智慧電錶資料讀到data.table。

DT<-as.data.table(read\_feather("ind"))

使用GAM迴歸模型。將工作日的字元轉換為整數，並使用recode包中的函式重新編碼工作日：1.星期一，…，7星期日。

1. DT\[,week_num:=as.integer(car::recode(week,
2. "'Monday'='1';'Tuesday'='2';'Wednesday'='3';'Thursday'='4';
3. 'Friday'='5';'Saturday'='6';'Sunday'='7'"))\]

將資訊儲存在日期變數中，以簡化工作。

1. n_type<-unique(DT\[,type\])
2. n_date<-unique(DT\[,date\])
3. n_weekdays<-unique(DT\[,week\])
4. period<-48

讓我們看一下用電量的一些資料並對其進行分析。

1. data\_r<-DT\[(type==n\_type\[1\]&date%in%n_date\[57:70\])\]
3. ggplot(data\_r,aes(date\_time,value))+
4. geom_line()+
5. theme(panel.border=element_blank(),
6. panel.background=element_blank(),
7. panel.grid.minor=element_line(colour="grey90"),
8. panel.grid.major=element_line(colour="grey90"),
9. panel.grid.major.x=element_line(colour="grey90"),
10. axis.text=element_text(size=10),
11. axis.title=element_text(size=12,face="bold"))+
12. labs(x="Date",y="Load(kW)")

在繪製的時間序列中可以看到兩個主要的季節性：每日和每週。我們在一天中有48個測量值，在一週中有7天，因此這將是我們用來對因變數–電力負荷進行建模的自變數。

訓練我們的第一個GAM。通過平滑函式s對自變數建模，對於每日季節性，使用三次樣條迴歸，對於每週季節性，使用P樣條。

1. gam_1<-gam(Load~s(Daily,bs="cr",k=period)+
2. s(Weekly,bs="ps",k=7),
3. data=matrix_gam,
4. family=gaussian)

首先是視覺化。

1. layout(matrix(1:2,nrow=1))
2. plot(gam_1,shade=TRUE)

我們在這裡可以看到變數對電力負荷的影響。在左圖中，白天的負載峰值約為下午3點。在右邊的圖中，我們可以看到在週末負載量減少了。

讓我們使用summary函式對第一個模型進行診斷。

1. ##
2. ##Family:gaussian
3. ##Linkfunction:identity
4. ##
5. ##Formula:
6. ##Load~s(Daily,bs="cr",k=period)+s(Weekly,bs="ps",
7. ##k=7)
8. ##
9. ##Parametriccoefficients:
10. ##EstimateStd.ErrortvaluePr(>|t|)
11. ##(Intercept)2731.6718.88144.7<2e-16***
12. ##---
13. ##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1
14. ##
15. ##Approximatesignificanceofsmoothterms:
16. ##edfRef.dfFp-value
17. ##s(Daily)10.15912.688119.8<2e-16***
18. ##s(Weekly)5.3115.758130.3<2e-16***
19. ##---
20. ##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1
21. ##
22. ##R-sq.(adj)=0.772Devianceexplained=77.7%
23. ##GCV=2.4554e+05Scaleest.=2.3953e+05n=672

EDF：估計的自由度–可以像對給定變數進行平滑處理那樣來解釋（較高的EDF值表示更復雜的樣條曲線）。P值：給定變數對因變數的統計顯著性，通過F檢驗進行檢驗（越低越好）。調整後的R平方（越高越好）。我們可以看到R-sq.（adj）值有點低。

讓我們繪製擬合值：

我們需要將兩個自變數的互動作用包括到模型中。

第一種互動型別對兩個變數都使用了一個平滑函式。

1. gam_2<-gam(Load~s(Daily,Weekly),
4. summary(gam_2)$r.sq

##\[1\]0.9352108

R方值表明結果要好得多。

summary(gam_2)$s.table

1. ##edfRef.dfFp-value
2. ##s(Daily,Weekly)28.700828.99423334.47540

似乎也很好，p值為0，這意味著自變數很重要。擬合值圖：

現在，讓我們嘗試上述變數互動。這可以通過function完成te，也可以定義基本函式。

##\[1\]0.9268452

與以前的模型相似gam_2。

summary(gam_3)$s.table

1. ##edfRef.dfFp-value
2. ##te(Daily,Weekly)23.6570923.98741354.58560

非常相似的結果。讓我們看一下擬合值：

與gam_2模型相比，只有一點點差異，看起來te擬合更好。

##\[1\]0.9727604

summary(gam_4)$sp.criterion

1. ##GCV.Cp
2. ##34839.46

summary(gam_4)$s.table

1. ##edfRef.dfFp-value
2. ##te(Daily,Weekly)119.4117149.6528160.20650

我們可以在這裡看到R方略有上升。
讓我們繪製擬合值：

這似乎比gam_3模型好得多。

##\[1\]0.965618

summary(gam\_4\_fx)$s.table

1. ##edfRef.dfFp-value
2. ##te(Daily,Weekly)33533557.253895.289648e-199

我們可以看到R平方比模型gam_4低，這是因為我們過度擬合了模型。證明lambda和EDF的估計工作正常。

因此，讓我們在案例（模型）中嘗試ti方法。

##\[1\]0.9717469

summary(gam_5)$sp.criterion

1. ##GCV.Cp
2. ##35772.35

summary(gam_5)$s.table

1. ##edfRef.dfFp-value
2. ##s(Daily)22.58364927.964970444.199620
3. ##s(Weekly)5.9145315.9959341014.724820
4. ##ti(Daily,Weekly)85.310314110.82881441.222880

然後使用t2。

##\[1\]0.9738273

summary(gam_6)$sp.criterion

1. ##GCV.Cp
2. ##32230.68

summary(gam_6)$s.table

1. ##edfRef.dfFp-value
2. ##t2(Daily,Weekly)98.12005120.234586.707540

我還輸出了最後三個模型的GCV得分值，這也是在一組擬合模型中選擇最佳模型的良好標準。我們可以看到，對於t2相應模型gam_6，GCV值最低。

在統計中廣泛使用的其他模型選擇標準是AIC（Akaike資訊準則）。讓我們看看三個模型：

AIC(gam\_4,gam\_5,gam_6)

1. ##dfAIC
2. ##gam_4121.41178912.611
3. ##gam_5115.80858932.746
4. ##gam_6100.12008868.628

最低值在gam_6模型中。讓我們再次檢視擬合值。

我們可以看到的模型的擬合值gam_4和gam_6非常相似。可以使用軟體包的更多視覺化和模型診斷功能來比較這兩個模型。

第一個是functiongam.check，它繪製了四個圖：殘差的QQ圖，線性預測變數與殘差，殘差的直方圖以及擬合值與因變數的關係圖。讓我們診斷模型gam_4和gam_6。

gam.check(gam_4)

1. ##
2. ##Method:GCVOptimizer:magic
3. ##Smoothingparameterselectionconvergedafter7iterations.
4. ##TheRMSGCVscoregradiantatconvergencewas0.2833304.
5. ##TheHessianwaspositivedefinite.
6. ##Theestimatedmodelrankwas336(maximumpossible:336)
7. ##Modelrank=336/336
8. ##
9. ##Basisdimension(k)checkingresults.Lowp-value(k-index<1)may
10. ##indicatethatkistoolow,especiallyifedfisclosetok'.
11. ##
12. ##k'edfk-indexp-value
13. ##te(Daily,Weekly)335.00119.411.221

gam.check(gam_6)

1. ##
2. ##Method:GCVOptimizer:magic
3. ##Smoothingparameterselectionconvergedafter9iterations.
4. ##TheRMSGCVscoregradiantatconvergencewas0.05208856.
5. ##TheHessianwaspositivedefinite.
6. ##Theestimatedmodelrankwas336(maximumpossible:336)
7. ##Modelrank=336/336
8. ##
9. ##Basisdimension(k)checkingresults.Lowp-value(k-index<1)may
10. ##indicatethatkistoolow,especiallyifedfisclosetok'.
11. ##
12. ##k'edfk-indexp-value
13. ##t2(Daily,Weekly)335.0098.121.181

我們可以再次看到模型非常相似，只是在直方圖中可以看到一些差異。

1. layout(matrix(1:2,nrow=1))
2. plot(gam_4,rug=FALSE,se=FALSE,n2=80,main="gamn.4withte()")
3. plot(gam_6,rug=FALSE,se=FALSE,n2=80,main="gamn.6witht2()")

該模型gam_6有更多的“波浪形”的輪廓。因此，這意味著它對因變數的擬合度更高，而平滑因子更低。

1. vis.gam(gam_6,n.grid=50,theta=35,phi=32,zlab="",
2. ticktype="detailed",color="topo",main="t2(D,W)")

我們可以看到最高峰值是Daily變數的值接近30（下午3點），而Weekly變數的值是1（星期一）。

1. vis.gam(gam_6,main="t2(D,W)",plot.type="contour",
2. 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模型進行時間序列預測