R語言廣義相加模型 (GAMs)分析預測CO2時間序列數據|附代碼數據

2023-02-23   tecdat拓端

原標題:R語言廣義相加模型 (GAMs)分析預測CO2時間序列數據|附代碼數據

全文下載連結:http://tecdat.cn/?p=20904

最近我們被客戶要求撰寫關於廣義相加模型 (GAMs)的研究報告,包括一些圖形和統計輸出。

環境科學中的許多數據不適合簡單的線性模型,最好用廣義相加模型(GAM)來描述

這基本上就是具有 光滑函數的廣義線性模型(GLM)的擴展 。當然,當您使用光滑項擬合模型時,可能會發生許多複雜的事情,但是您只需要了解基本原理即可。

理論

讓我們從高斯線性模型的方程開始 :

GAM中發生的變化是存在光滑項:

這僅意味著對線性預測變量的貢獻現在是函數f。從概念上講,這與使用二次項()或三次項()作為預測變量沒什麼不同。

在這裡,我們將重點放在樣條曲線上。在過去,它可能類似於分段線性函數。

例如,您可以在模型中包含線性項和光滑項的組合

或者我們可以擬合廣義分布和隨機效應

一個簡單的例子

讓我們嘗試一個簡單的例子。首先,讓我們創建一個數據框,並創建一些具有明顯非線性趨勢的模擬數據,並比較一些模型對該數據的擬合程度。

x <- seq(0, pi * 2, 0.1)

sin_x <- sin(x)

y <- sin_x + rnorm(n = length(x), mean = 0, sd = sd(sin_x / 2))

Sample <- data.frame(y,x)

library(ggplot2)

ggplot(Sample, aes(x, y)) + geom_point()

嘗試擬合普通的線性模型:

lm_y <- lm(y ~ x, data = Sample)

並使用geom_smooth in 繪製帶有數據的擬合線 ggplot

ggplot(Sample, aes(x, y)) + geom_point() + geom_smooth(method = lm)

查看圖或 summary(lm_y),您可能會認為模型擬合得很好,但請查看殘差圖

plot(lm_y, which = 1)

顯然,殘差未均勻分布在x的值上,因此我們需要考慮一個更好的模型。

點擊標題查閱往期內容

【視頻】廣義相加模型(GAM)在電力負荷預測中的應用

左右滑動查看更多

01

02

03

04

運行分析

在R中運行GAM。

要運行GAM,我們使用:

gam_y <- gam(y ~ s(x), method = "REML")

要提取擬合值,我們可以predict :

predict(gam_y, data.frame(x = x_new))

但是對於簡單的模型,我們還可以利用中的 method = 參數來 geom_smooth指定模型公式。

您可以看到該模型更適合數據,檢查診斷信息。

check.gam 快速簡便地查看殘差圖。

gam.check(gam_y)

##

## Method: REML Optimizer: outer newton

## full convergence after 6 iterations.

## Gradient range [-2.37327e-09,1.17425e-09]

## (score 44.14634 & scale 0.174973).

## Hessian positive definite, eigenvalue range [1.75327,30.69703].

## 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 5.76 1.19 0.9

對模型對象使用summary將為您提供光滑項(以及任何參數項)的意義,以及解釋的方差。在這個例子中,非常合適。「edf」是估計的自由度——本質上,數量越大,擬合模型就越搖擺。大約為1的值趨向於接近線性項。

##

## Family: gaussian

## Link function: identity

##

## Formula:

## y ~ s(x)

##

## Parametric coefficients:

## Estimate Std. Error t value Pr(>|t|)

## (Intercept) -0.01608 0.05270 -0.305 0.761

##

## Approximate significance of smooth terms:

## edf Ref.df F p-value

## s(x) 5.76 6.915 23.38 <2e-16 ***

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##

## R-sq.(adj) = 0.722 Deviance explained = 74.8%

## -REML = 44.146 Scale est. = 0.17497 n = 63

光滑函數項

如上所述,我們將重點介紹樣條曲線,因為樣條曲線是最常實現的光滑函數(非常快速且穩定)。那麼,當我們指定s(x)時實際發生了什麼 ?

好吧,這就是我們說要把y擬合為x個函數集的線性函數的地方。默認輸入為薄板回歸樣條-您可能會看到的常見樣條是三次回歸樣條。三次回歸樣條曲線具有 我們在談論樣條曲線時想到的傳統 結點–在這種情況下,它們均勻分布在協變量範圍內。

基函數

我們將從擬合模型開始,記住光滑項是一些函數的和,

首先,我們提取_基本函數_集 (即光滑項的bj(xj)部分)。然後我們可以畫出第一和第二基函數。

model_matrix <- predict(gam_y, type = "lpmatrix")

plot(y ~ x)

現在,讓我們繪製所有基函數的圖,然後再將其添加到GAM(y_pred)的預測中。

matplot(x, model_matrix[,-1], type = "l", lty = 2, add = T)

lines(y_pred ~ x_new, col = "red", lwd = 2)

現在,最容易想到這樣-每條虛線都代表一個函數(bj),據此 gam 估算係數(βj),將它們相加即可得出對應的f(x)的貢獻(即先前的等式)。對於此示例而言,它很好且簡單,因為我們僅根據光滑項對y進行建模,因此它是相當相關的。順便說一句,您也可以只使用 plot.gam 繪製光滑項。

好的,現在讓我們更詳細地了解基函數的構造方式。您會看到函數的構造與因變量數據是分開的。為了證明這一點,我們將使用 smoothCon。

x_sin_smooth <- smoothCon(s(x), data = data.frame(x), absorb.cons = TRUE)

現在證明您可以從基本函數和估計係數到擬合的光滑項。再次注意,這裡簡化了,因為模型只是一個光滑項。如果您有更多的項,我們需要將線性預測模型中的所有項相加。

betas <- gam_y$coefficients

linear_pred <- model_matrix %*% betas

請看下面的圖,記住這 X 是基函數的矩陣。

通過 gam.models , smooth.terms 光滑模型類型的所有選項,基本函數的構造方式(懲罰等),我們可以指定的模型類型(隨機效應,線性函數,交互作用)。

真實例子

我們查看一些CO2數據,為數據擬合幾個GAM,以嘗試區分年度內和年度間趨勢。

首先加載數據 。

CO2 <- read.csv("co2.csv")

我們想首先查看年趨勢,因此讓我們將日期轉換為連續的時間變量(採用子集進行可視化)。

CO2$time <- as.integer(as.Date(CO2$Date, format = "%d/%m/%Y"))

我們來繪製它,並考慮一個平穩的時間項。

我們為這些數據擬合GAM

它擬合具有單個光滑時間項的模型。我們可以查看以下預測值:

plot(CO2_time)

請注意光滑項如何減少到「普通」線性項的(edf為1)-這是懲罰回歸樣條曲線的優點。但如果我們檢查一下模型,就會發現有些東西是混亂的。

par(mfrow = c(2,2))

gam.check(CO2_time)

殘差圖的上升和下降模式看起來很奇怪-顯然存在某種依賴關係結構(我們可能會猜測,這與年內波動有關)。讓我們再試一次,並引入一種稱為周期光滑項。

周期性光滑項fintrannual(month)由基函數組成,與我們已經看到的相同,只是樣條曲線的端點被約束為相等,這在建模時是有意義的周期性(跨月/跨年)的變量。

現在,我們將看到 bs = 用於選擇光滑器類型的k = 參數和用於選擇結數的 參數,因為三次回歸樣條曲線具有固定的結數。我們使用12結,因為有12個月。

s(month, bs = 'cc', k = 12) + s(time)

讓我們看一下擬合的光滑項:

從這兩個光滑項來看,我們可以看到,月度光滑項檢測到CO2濃度的月度上升和下降——從相對幅度(即月度波動與長期趨勢)來看,我們可以看出消除時間序列成分是多麼重要。讓我們看看現在的模型診斷是怎樣的:

par(mfrow = c(2,2))

gam.check(CO2_season_time)

好多了。讓我們看一下季節性因素如何與整個長期趨勢相對應。

plot(CO2_season_time)

結果

從本質上講,您可以將GAM的模型結果表示為任何其他線性模型,主要區別在於,對於光滑項,沒有單一係數可供推斷(即負、正、效應大小等)。因此,您需要依靠視覺上解釋光滑項(例如從對plot(gam_model)的調用)或根據預測值進行推斷。當然,你可以在模型中包含普通的線性項(無論是連續的還是分類的,甚至在方差分析類型的框架中),並像平常一樣從中進行推斷。事實上,GAM對於解釋一個非線性現象通常是有用的,這個非線性現象並不直接引起人們的興趣,但在推斷其他變量時需要加以解釋。

您可以通過plot 在擬合的gam模型上調用函數來繪製局部效果 ,還可以查看參數項,也可以使用 termplot 函數。您可以ggplot 像本教程前面所述那樣使用 簡單的模型,但是對於更複雜的模型,最好知道如何使用predict預測數據 。

geom_line(aes(y = predicted_values)

本文摘選 R語言廣義相加模型 (GAMs)分析預測CO2時間序列數據 ,點擊「閱讀原文」獲取全文完整資料。

點擊標題查閱往期內容

【視頻】廣義相加模型(GAM)在電力負荷預測中的應用

分位數回歸、GAM樣條曲線、指數平滑和SARIMA對電力負荷時間序列預測

實現廣義相加模型GAM和普通最小二乘(OLS)回歸

R語言非參數模型厘定保險費率:局部回歸、廣義相加模型GAM、樣條回歸

R語言廣義加性模型GAMs分析溫度、臭氧環境數據繪製偏回歸圖與偏殘差圖

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回歸