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

原文出处：拓端数据部落公众号

最近我们被客户要求撰写关于GAMs的研究报告，包括一些图形和统计输出。

我们在研究工作中使用广义加性模型（GAMs）。mgcv软件包是一套优秀的软件，可以为非常大的数据集指定、拟合和可视化GAMs。

这篇文章介绍一下广义加性模型（GAMs）目前可以实现的功能。

我们需要加载mgcv

library('mgcv')

受欢迎的例子数据集

dat中的数据在GAM相关的研究中得到了很好的研究，包含了一些协变量--标记为x0到x3--这些协变量在不同程度上与因变量有非线性关系。

我们想通过使用样条来逼近协变量和因变量之间的真实关系来尝试拟合这些关系。为了拟合一个加性模型，我们使用

 gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), dat,  "REML")

mgcv提供了一个summary()方法，用来提取关于拟合GAM的信息。

check()函数，用于检查模型中的每个光滑函数是否使用了足够数量的基函数。你可能没有直接使用check()——会输出其他诊断结果，也会产生四个模型诊断图。

绘制光滑函数图

为了将估计的GAMs可视化，mgcv提供了plot.gam()方法和vis.gam()函数，从对象中产生类似ggplot2的图。为了使GAM模型中的四个估计光滑函数可视化，我们将使用

plot(mod)

结果是绘制mod GAM中每一个光滑函数。

使用plot函数在绘图设备上绘制多个面板，并将各个绘图排成一行。

提取光滑函数数据

用于处理mod中表示的基础光滑函数，如果你想提取用于构建该图的大部分数据，你可以使用smooth()函数。

smooth(mod, "x1")

诊断图

由check()产生的诊断图

check(mod)

结果是一个包含四个诊断图的数组，包括模型残差的Q-Q图（左上）和直方图（左下），残差与线性预测器的图（右上），以及观察值与拟合值的图。

这四张图中的每一张都是通过用户可访问的函数生成的，函数实现了一个特定的图。例如，qqplot(mod)产生上图左上方的Q-Q图。

qqplot(mod)

qqplot(mod)的结果是一个残差的Q-Q图，其中的参考量值是通过模拟拟合模型的数据而得到。

还可以处理目前可用的许多更专业的光滑函数。例如，二维光滑函数。


plot(mod)

二维光滑函数的默认绘制方式是使用plot()。

和因子光滑函数交互项，相当于光滑曲线的随机斜率和截距，被画在一个面板上，颜色被用来区分不同的随机光滑函数。


## 模拟数据
f0 <- function(x) 2 * sin(pi * x)
f1 <- function(x, a=2, b=-1) exp(a * x)+b
f2 <- function(x) 0.2 * x^11 * (10 * (1 - x))^6 + 10 * (10 * x)^3 * (1 - x)^10
f <- f0(x0) + f1(x1, a[fac], b[fac]) + f2(x2)
fac <- factor(fac)
y <- f + rnorm(n) * 2


plot(mod)

含有因子-光滑函数交互项的更复杂的GAM的结果，bs = 'fs'。

还能做什么？

可以处理mgcv可以估计的大多数光滑函数，包括带有因子和连续副变量的按变量光滑函数、随机效应光滑函数（bs = 're'）、二维张量积光滑函数，以及带有参数项的模型。

参考文献

Augustin, N. H., Sauleau, E.-A., and Wood, S. N. (2012). On quantile quantile plots for generalized linear models. Computational statistics & data analysis 56, 2404–2409. doi:10.1016/j.csda.2012.01.026.