原文链接:http:///?p=20904环境科学中的许多数据不适合简单的线性模型,最好用广义相加模型(GAM)来描述。 这基本上就是具有 光滑函数的广义线性模型(GLM)的扩展 。当然,当您使用光滑项拟合模型时,可能会发生许多复杂的事情,但是您只需要了解基本原理即可。 理论让我们从高斯线性模型的方程开始 : GAM中发生的变化是存在光滑项: 这仅意味着对线性预测变量的贡献现在是函数f。从概念上讲,这与使用二次项()或三次项()作为预测变量没什么不同。 在这里,我们将重点放在样条曲线上。在过去,它可能类似于分段线性函数。 例如,您可以在模型中包含线性项和光滑项的组合 或者我们可以拟合广义分布和随机效应 一个简单的例子让我们尝试一个简单的例子。首先,让我们创建一个数据框,并创建一些具有明显非线性趋势的模拟数据,并比较一些模型对该数据的拟合程度。 x <- seq(0, pi * 2, 0.1) library(ggplot2) 尝试拟合普通的线性模型: lm_y <- lm(y ~ x, data = Sample) 并使用 ggplot(Sample, aes(x, y)) + geom_point() + geom_smooth(method = lm) 查看图或 plot(lm_y, which = 1) 显然,残差未均匀分布在x的值上,因此我们需要考虑一个更好的模型。 运行分析在R中运行GAM。 要运行GAM,我们使用: gam_y <- gam(y ~ s(x), method = "REML") 要提取拟合值,我们可以 predict(gam_y, data.frame(x = x_new)) 但是对于简单的模型,我们还可以利用中的 您可以看到该模型更适合数据,检查诊断信息。
gam.check(gam_y) ## 对模型对象使用summary将为您提供光滑项(以及任何参数项)的意义,以及解释的方差。在这个例子中,非常合适。“edf”是估计的自由度——本质上,数量越大,拟合模型就越摇摆。大约为1的值趋向于接近线性项。 ## 光滑函数项如上所述,我们将重点介绍样条曲线,因为样条曲线是最常实现的光滑函数(非常快速且稳定)。那么,当我们指定 好吧,这就是我们说要把y拟合为x个函数集的线性函数的地方。默认输入为薄板回归样条-您可能会看到的常见样条是三次回归样条。三次回归样条曲线具有 我们在谈论样条曲线时想到的传统 _结点_–在这种情况下,它们均匀分布在协变量范围内。 基函数我们将从拟合模型开始,记住光滑项是一些函数的和, 首先,我们提取_基本函数_集 (即光滑项的bj(xj)部分)。然后我们可以画出第一和第二基函数。 model_matrix <- predict(gam_y, type = "lpmatrix") 现在,让我们绘制所有基函数的图,然后再将其添加到GAM( matplot(x, model_matrix[,-1], type = "l", lty = 2, add = T) 现在,最容易想到这样-每条虚线都代表一个函数(bj),据此 好的,现在让我们更详细地了解基函数的构造方式。您会看到函数的构造与因变量数据是分开的。为了证明这一点,我们将使用 x_sin_smooth <- smoothCon(s(x), data = data.frame(x), absorb.cons = TRUE) 现在证明您可以从基本函数和估计系数到拟合的光滑项。再次注意,这里简化了,因为模型只是一个光滑项。如果您有更多的项,我们需要将线性预测模型中的所有项相加。 betas <- gam_y$coefficients 请看下面的图,记住这 通过 真实例子我们查看一些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)) 残差图的上升和下降模式看起来很奇怪-显然存在某种依赖关系结构(我们可能会猜测,这与年内波动有关)。让我们再试一次,并引入一种称为周期光滑项。 周期性光滑项fintrannual(month)由基函数组成,与我们已经看到的相同,只是样条曲线的端点被约束为相等,这在建模时是有意义的周期性(跨月/跨年)的变量。 现在,我们将看到 s(month, bs = 'cc', k = 12) + s(time) 让我们看一下拟合的光滑项: 从这两个光滑项来看,我们可以看到,月度光滑项检测到CO2浓度的月度上升和下降——从相对幅度(即月度波动与长期趋势)来看,我们可以看出消除时间序列成分是多么重要。让我们看看现在的模型诊断是怎样的: par(mfrow = c(2,2)) 好多了。让我们看一下季节性因素如何与整个长期趋势相对应。 plot(CO2_season_time) 结果从本质上讲,您可以将GAM的模型结果表示为任何其他线性模型,主要区别在于,对于光滑项,没有单一系数可供推断(即负、正、效应大小等)。因此,您需要依靠视觉上解释光滑项(例如从对plot(gam_model)的调用)或根据预测值进行推断。当然,你可以在模型中包含普通的线性项(无论是连续的还是分类的,甚至在方差分析类型的框架中),并像平常一样从中进行推断。事实上,GAM对于解释一个非线性现象通常是有用的,这个非线性现象并不直接引起人们的兴趣,但在推断其他变量时需要加以解释。 您可以通过 geom_line(aes(y = predicted_values) |
|