拓端資料tecdat|R語言廣義相加模型 (GAMs)分析預測CO2時間序列資料

環境科學中的許多資料不适合簡單的線性模型，最好用廣義相加模型（GAM）來描述。

這基本上就是具有 光滑函數的​​廣義線性模型（GLM）​​的擴充 。當然，當您使用光滑項拟合模型時，可能會發生許多複雜的事情，但是您隻需要了解基本原理即可。

理論

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

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

這僅意味着對線性預測變量的貢獻現在是函數f。從概念上講，這與使用二次項（

）或三次項（

）作為預測變量沒什麼不同。

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

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

或者我們可以拟合廣義分布和随機效應

一個簡單的例子

讓我們嘗試一個簡單的例子。首先，讓我們建立一個資料框，并建立一些具有明顯非線性趨勢的模拟資料，并比較一些模型對該資料的拟合程度。

1.  x <- seq(0, pi * 2, 0.1)
2.  sin_x <- sin(x)
3.  y <- sin_x + rnorm(n = length(x), mean = 0, sd = sd(sin_x / 2))
4.  Sample <- data.frame(y,x)
5.  library(ggplot2)
6.  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的值上，是以我們需要考慮一個更好的模型。

運作分析

在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)      

1.  ##
2.  ## Method: REML Optimizer: outer newton
3.  ## full convergence after 6 iterations.
4.  ## Gradient range [-2.37327e-09,1.17425e-09]
5.  ## (score 44.14634 & scale 0.174973).
6.  ## Hessian positive definite, eigenvalue range [1.75327,30.69703].
7.  ## Model rank = 10 / 10
8.  ##
9.  ## Basis dimension (k) checking results. Low p-value (k-index<1) may
10.  ## indicate that k is too low, especially if edf is close to k'.
11.  ##
12.  ## k' edf k-index p-value
13.  ## s(x) 9.00 5.76 1.19 0.9      

對模型對象使用summary将為您提供光滑項（以及任何參數項）的意義，以及解釋的方差。在這個例子中，非常合适。“edf”是估計的自由度——本質上，數量越大，拟合模型就越搖擺。大約為1的值趨向于接近線性項。

1.  ##
2.  ## Family: gaussian
3.  ## Link function: identity
4.  ##
5.  ## Formula:
6.  ## y ~ s(x)
7.  ##
8.  ## Parametric coefficients:
9.  ## Estimate Std. Error t value Pr(>|t|)
10.  ## (Intercept) -0.01608 0.05270 -0.305 0.761
11.  ##
12.  ## Approximate significance of smooth terms:
13.  ## edf Ref.df F p-value
14.  ## s(x) 5.76 6.915 23.38 <2e-16 ***
15.  ## ---
16.  ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
17.  ##
18.  ## R-sq.(adj) = 0.722 Deviance explained = 74.8%
19.  ## -REML = 44.146 Scale est. = 0.17497 n = 63      

光滑函數項

如上所述，我們将重點介紹樣條曲線，因為樣條曲線是最常實作的光滑函數（非常快速且穩定）。那麼，當我們指定​

​s(x)​

​時實際發生了什麼 ？

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

基函數

我們将從拟合模型開始，記住光滑項是一些函數的和，

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

1.  model_matrix <- predict(gam_y, type = "lpmatrix")
2.  plot(y ~ x)

現在，讓我們繪制所有基函數的圖，然後再将其添加到GAM（​

​y_pred​

​）的預測中。

1.  matplot(x, model_matrix[,-1], type = "l", lty = 2, add = T)
2.  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)      

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

1.  betas <- gam_y$coefficients
2.  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）-這是懲罰回歸樣條曲線的優點。但如果我們檢查一下模型，就會發現有些東西是混亂的。

1.  par(mfrow = c(2,2))
2.  gam.check(CO2_time)

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

周期性光滑項fintrannual（month）由基函數組成，與我們已經看到的相同，隻是樣條曲線的端點被限制為相等，這在模組化時是有意義的周期性（跨月/跨年）的變量。

現在，我們将看到 ​

​bs =​

​​ 用于選擇光滑器類型的​

​k =​

​ 參數和用于選擇結數的 參數，因為三次回歸樣條曲線具有固定的結數。我們使用12結，因為有12個月。

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

讓我們看一下拟合的光滑項：

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

1.  par(mfrow = c(2,2))
2.  gam.check(CO2_season_time)

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

1.  plot(CO2_season_time)

結果

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

您可以通過​

​plot​

​​ 在拟合的gam模型上調用函數來繪制局部效果 ，還可以檢視參數項，也可以使用 ​

​termplot​

​​ 函數。您可以​

​ggplot​

​​ 像本教程前面所述那樣使用 簡單的模型，但是對于更複雜的模型，最好知道如何使用​

​predict預測​

​資料 。

geom_line(aes(y = predicted_values)