关于算法:R语言广义相加模型-GAMs分析预测CO2时间序列数据

51次阅读

共计 4262 个字符,预计需要花费 11 分钟才能阅读完成。

原文链接:http://tecdat.cn/?p=20904 

环境迷信中的许多数据不适宜简略的线性模型,最好用狭义相加模型(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 的值上,因而咱们须要思考一个更好的模型。

运行剖析

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


最受欢迎的见解

1. 用 SPSS 预计 HLM 档次线性模型模型

2.R 语言线性判别分析(LDA),二次判别分析(QDA)和正则判别分析(RDA)

3. 基于 R 语言的 lmer 混合线性回归模型

4.R 语言 Gibbs 抽样的贝叶斯简略线性回归仿真剖析

5. 在 r 语言中应用 GAM(狭义相加模型)进行电力负荷工夫序列剖析

6. 应用 SAS,Stata,HLM,R,SPSS 和 Mplus 的分层线性模型 HLM

7.R 语言中的岭回归、套索回归、主成分回归:线性模型抉择和正则化

8.R 语言用线性回归模型预测空气质量臭氧数据

9.R 语言分层线性模型案例

正文完
 0