散文網(wǎng) » 生活 »日常 » 拓端tecdat|R語(yǔ)言廣義相加模型 (GAMs)分析預(yù)測(cè)CO2時(shí)間序列數(shù)據(jù)

拓端tecdat|R語(yǔ)言廣義相加模型 (GAMs)分析預(yù)測(cè)CO2時(shí)間序列數(shù)據(jù)

2021-07-15 11:01 作者:拓端tecdat 0人讀過 | 我要投稿

原文 |?http://tecdat.cn/?p=20904?
來源 |?拓端數(shù)據(jù)部落公眾號(hào)

環(huán)境科學(xué)中的許多數(shù)據(jù)不適合簡(jiǎn)單的線性模型，最好用廣義相加模型（GAM）來描述。

這基本上就是具有?光滑函數(shù)的廣義線性模型（GLM）的擴(kuò)展?。當(dāng)然，當(dāng)您使用光滑項(xiàng)擬合模型時(shí)，可能會(huì)發(fā)生許多復(fù)雜的事情，但是您只需要了解基本原理即可。

理論

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

GAM中發(fā)生的變化是存在光滑項(xiàng)：

這僅意味著對(duì)線性預(yù)測(cè)變量的貢獻(xiàn)現(xiàn)在是函數(shù)f。從概念上講，這與使用二次項(xiàng)（）或三次項(xiàng)（）作為預(yù)測(cè)變量沒什么不同。

在這里，我們將重點(diǎn)放在樣條曲線上。在過去，它可能類似于分段線性函數(shù)。

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

或者我們可以擬合廣義分布和隨機(jī)效應(yīng)

一個(gè)簡(jiǎn)單的例子

讓我們嘗試一個(gè)簡(jiǎn)單的例子。首先，讓我們創(chuàng)建一個(gè)數(shù)據(jù)框，并創(chuàng)建一些具有明顯非線性趨勢(shì)的模擬數(shù)據(jù)，并比較一些模型對(duì)該數(shù)據(jù)的擬合程度。


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?繪制帶有數(shù)據(jù)的擬合線?ggplot

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

查看圖或?summary(lm_y)，您可能會(huì)認(rèn)為模型擬合得很好，但請(qǐng)查看殘差圖

plot(lm_y, which = 1)

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

運(yùn)行分析

在R中運(yùn)行GAM。

要運(yùn)行GAM，我們使用：

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

要提取擬合值，我們可以predict??：

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

但是對(duì)于簡(jiǎn)單的模型，我們還可以利用中的?method =?參數(shù)來?geom_smooth指定模型公式。

您可以看到該模型更適合數(shù)據(jù)，檢查診斷信息。

check.gam?快速簡(jiǎn)便地查看殘差圖。

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

對(duì)模型對(duì)象使用summary將為您提供光滑項(xiàng)（以及任何參數(shù)項(xiàng)）的意義，以及解釋的方差。在這個(gè)例子中，非常合適?！癳df”是估計(jì)的自由度——本質(zhì)上，數(shù)量越大，擬合模型就越搖擺。大約為1的值趨向于接近線性項(xiàng)。


##
## 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

光滑函數(shù)項(xiàng)

如上所述，我們將重點(diǎn)介紹樣條曲線，因?yàn)闃訔l曲線是最常實(shí)現(xiàn)的光滑函數(shù)（非?？焖偾曳€(wěn)定）。那么，當(dāng)我們指定s(x)時(shí)實(shí)際發(fā)生了什么?？

好吧，這就是我們說要把y擬合為x個(gè)函數(shù)集的線性函數(shù)的地方。默認(rèn)輸入為薄板回歸樣條-您可能會(huì)看到的常見樣條是三次回歸樣條。三次回歸樣條曲線具有?我們?cè)谡務(wù)摌訔l曲線時(shí)想到的傳統(tǒng)?–在這種情況下，它們均勻分布在協(xié)變量范圍內(nèi)。

基函數(shù)

我們將從擬合模型開始，記住光滑項(xiàng)是一些函數(shù)的和，

首先，我們提取集??（即光滑項(xiàng)的bj（xj）部分）。然后我們可以畫出第一和第二基函數(shù)。


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

現(xiàn)在，讓我們繪制所有基函數(shù)的圖，然后再將其添加到GAM（y_pred）的預(yù)測(cè)中。


matplot(x, model_matrix[,-1], type = "l", lty = 2, add = T)
lines(y_pred ~ x_new, col = "red", lwd = 2)

現(xiàn)在，最容易想到這樣-每條虛線都代表一個(gè)函數(shù)（bj），據(jù)此?gam?估算系數(shù)（βj），將它們相加即可得出對(duì)應(yīng)的f（x）的貢獻(xiàn)（即先前的等式）。對(duì)于此示例而言，它很好且簡(jiǎn)單，因?yàn)槲覀儍H根據(jù)光滑項(xiàng)對(duì)y進(jìn)行建模，因此它是相當(dāng)相關(guān)的。順便說一句，您也可以只使用?plot.gam?繪制光滑項(xiàng)。

好的，現(xiàn)在讓我們更詳細(xì)地了解基函數(shù)的構(gòu)造方式。您會(huì)看到函數(shù)的構(gòu)造與因變量數(shù)據(jù)是分開的。為了證明這一點(diǎn)，我們將使用?smoothCon。

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

現(xiàn)在證明您可以從基本函數(shù)和估計(jì)系數(shù)到擬合的光滑項(xiàng)。再次注意，這里簡(jiǎn)化了，因?yàn)槟Ｐ椭皇且粋€(gè)光滑項(xiàng)。如果您有更多的項(xiàng)，我們需要將線性預(yù)測(cè)模型中的所有項(xiàng)相加。


betas <- gam_y$coefficients
linear_pred <- model_matrix %*% betas

請(qǐng)看下面的圖，記住這?X?是基函數(shù)的矩陣。

通過?gam.models?，?smooth.terms?光滑模型類型的所有選項(xiàng)，基本函數(shù)的構(gòu)造方式（懲罰等），我們可以指定的模型類型（隨機(jī)效應(yīng)，線性函數(shù)，交互作用）。

真實(shí)例子

我們查看一些CO2數(shù)據(jù)，為數(shù)據(jù)擬合幾個(gè)GAM，以嘗試區(qū)分年度內(nèi)和年度間趨勢(shì)。

首先加載數(shù)據(jù)?。

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

我們想首先查看年趨勢(shì)，因此讓我們將日期轉(zhuǎn)換為連續(xù)的時(shí)間變量（采用子集進(jìn)行可視化）。

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

我們來繪制它，并考慮一個(gè)平穩(wěn)的時(shí)間項(xiàng)。

我們?yōu)檫@些數(shù)據(jù)擬合GAM

它擬合具有單個(gè)光滑時(shí)間項(xiàng)的模型。我們可以查看以下預(yù)測(cè)值：

plot(CO2_time)

請(qǐng)注意光滑項(xiàng)如何減少到“普通”線性項(xiàng)的（edf為1）-這是懲罰回歸樣條曲線的優(yōu)點(diǎn)。但如果我們檢查一下模型，就會(huì)發(fā)現(xiàn)有些東西是混亂的。


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

殘差圖的上升和下降模式看起來很奇怪-顯然存在某種依賴關(guān)系結(jié)構(gòu)（我們可能會(huì)猜測(cè)，這與年內(nèi)波動(dòng)有關(guān)）。讓我們?cè)僭囈淮?，并引入一種稱為周期光滑項(xiàng)。

周期性光滑項(xiàng)fintrannual（month）由基函數(shù)組成，與我們已經(jīng)看到的相同，只是樣條曲線的端點(diǎn)被約束為相等，這在建模時(shí)是有意義的周期性（跨月/跨年）的變量。

現(xiàn)在，我們將看到?bs =?用于選擇光滑器類型的k =?參數(shù)和用于選擇結(jié)數(shù)的?參數(shù)，因?yàn)槿位貧w樣條曲線具有固定的結(jié)數(shù)。我們使用12結(jié)，因?yàn)橛?2個(gè)月。

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

讓我們看一下擬合的光滑項(xiàng)：

從這兩個(gè)光滑項(xiàng)來看，我們可以看到，月度光滑項(xiàng)檢測(cè)到CO2濃度的月度上升和下降——從相對(duì)幅度（即月度波動(dòng)與長(zhǎng)期趨勢(shì)）來看，我們可以看出消除時(shí)間序列成分是多么重要。讓我們看看現(xiàn)在的模型診斷是怎樣的：


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

好多了。讓我們看一下季節(jié)性因素如何與整個(gè)長(zhǎng)期趨勢(shì)相對(duì)應(yīng)。


plot(CO2_season_time)

結(jié)果

從本質(zhì)上講，您可以將GAM的模型結(jié)果表示為任何其他線性模型，主要區(qū)別在于，對(duì)于光滑項(xiàng)，沒有單一系數(shù)可供推斷（即負(fù)、正、效應(yīng)大小等）。因此，您需要依靠視覺上解釋光滑項(xiàng)（例如從對(duì)plot（gam_model）的調(diào)用）或根據(jù)預(yù)測(cè)值進(jìn)行推斷。當(dāng)然，你可以在模型中包含普通的線性項(xiàng)（無論是連續(xù)的還是分類的，甚至在方差分析類型的框架中），并像平常一樣從中進(jìn)行推斷。事實(shí)上，GAM對(duì)于解釋一個(gè)非線性現(xiàn)象通常是有用的，這個(gè)非線性現(xiàn)象并不直接引起人們的興趣，但在推斷其他變量時(shí)需要加以解釋。

您可以通過plot?在擬合的gam模型上調(diào)用函數(shù)來繪制局部效果?，還可以查看參數(shù)項(xiàng)，也可以使用?termplot?函數(shù)。您可以ggplot?像本教程前面所述那樣使用?簡(jiǎn)單的模型，但是對(duì)于更復(fù)雜的模型，最好知道如何使用predict預(yù)測(cè)數(shù)據(jù)?。