全文下載鏈接：http://tecdat.cn/?p=9024

最近我們被客戶要求撰寫關于廣義相加模型（GAM）的研究報告，包括一些圖形和統(tǒng)計輸出。

這篇文章探討了為什么使用廣義相加模型?是一個不錯的選擇。為此，我們首先需要看一下線性回歸，看看為什么在某些情況下它可能不是最佳選擇。

回歸模型

假設我們有一些帶有兩個屬性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）使用光滑函數（如樣條曲線）作為回歸模型中的預測因子。

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

摘自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模型輸出：

查看模型摘要：

##?##?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點。在右邊的圖中，我們可以看到在周末負載量減少了。

點擊標題查閱往期內容

R語言對用電負荷時間序列數據進行K-medoids聚類建模和GAM回歸

左右滑動查看更多

01

02

03

04

讓我們使用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，這意味著自變量很重要。擬合值圖：

現在，讓我們嘗試上述變量交互。這可以通過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低，這是因為我們過度擬合了模型。證明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值最低。

在統(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???????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，直到星期四都非常相似，然后負荷在周末減少。

本文摘選?《?R語言廣義相加模型（GAM）在電力負荷預測中的應用?》?，點擊“閱讀原文”獲取全文完整資料。

點擊標題查閱往期內容

分位數回歸、GAM樣條曲線、指數平滑和SARIMA對電力負荷時間序列預測
實現廣義相加模型GAM和普通最小二乘(OLS)回歸
R語言非參數模型厘定保險費率：局部回歸、廣義相加模型GAM、樣條回歸
R語言廣義加性模型GAMs分析溫度、臭氧環(huán)境數據繪制偏回歸圖與偏殘差圖
R語言廣義相加（加性）模型（GAMs）與光滑函數可視化
R語言里的非線性模型：多項式回歸、局部樣條、平滑樣條、 廣義相加模型GAM分析
R語言用標準最小二乘OLS，廣義相加模型GAM?，樣條函數進行邏輯回歸LOGISTIC分類
R語言ISLR工資數據進行多項式回歸和樣條回歸分析
R語言中的多項式回歸、局部回歸、核平滑和平滑樣條回歸模型
R語言用泊松Poisson回歸、GAM樣條曲線模型預測騎自行車者的數量
R語言分位數回歸、GAM樣條曲線、指數平滑和SARIMA對電力負荷時間序列預測
R語言中的多項式回歸、B樣條曲線(B-spline Curves)回歸
R語言廣義相加模型 (GAMs)分析預測CO2時間序列數據
R語言中實現廣義相加模型GAM和普通最小二乘(OLS)回歸
在r語言中使用GAM（廣義相加模型）進行電力負荷時間序列分析
R語言用泊松Poisson回歸、GAM樣條曲線模型預測騎自行車者的數量
Python用廣義加性模型GAM進行時間序列分析
R語言廣義線性模型GLM、多項式回歸和廣義可加模型GAM預測泰坦尼克號幸存者
R語言中的廣義線性模型（GLM）和廣義相加模型（GAM）：多元（平滑）回歸分析保險資金投資組合信用風險敞口
R語言對用電負荷時間序列數據進行K-medoids聚類建模和GAM回歸
對用電負荷時間序列數據進行K-medoids聚類建模和GAM回歸