R语言基于mgcv包进行广义可加模型及交互作用演示（2）

广义可加模型(generalized additive models，GAMs)是广义线性模型和可加模型的结合，由 Hastie T 和 Tibshirani Ｒ于1986 年首先提出，其不要求应变量与自变量满足线性关系，适用于非线性数据的研究。既往我们已经介绍了使用VGAM包进行广义可加模型，今天我们继续来介绍使用mgcv包进行广义可加模型，这里我们使用我们原来的臭氧的数据（公众号回复：臭氧可以获得数据），主要描述的是臭氧浓度和大气一些相关指标的情况，因为有些数据是非线性的，使用Logistic回归不合适我们先导入数据看一下

library(foreign)
library(mgcv)
be <- read.spss("E:/r/test/ozone.sav",
                use.value.labels=F, to.data.frame=T)
names(be)

数据中有七个变量，ozon每日臭氧水平为结局变量，Inversion base height（ibh）反转基准高度，Pressure gradient (mm Hg) 压力梯度（mm Hg），Visibility (miles) 能见度（英里），Temperature (degrees F) 温度（华氏度），Day of the year日期，vh我也不知道是什么，反正就是一参数，这里所有的变量都是连续的。
我们先来进行单因素分析，假设我们想了解臭氧浓度和反转基准高度的关系

mgam<-gam(ozon~dpg+s(ibh),data=bc,family= poisson (),model=T)
summary(mgam)

截距项Intercept的Z=101.64，P小于0.05，表明截距项有统计学意义

假设检验结果表明平滑函数s(ibh)有统计意义，（X=687.5，P小于0.05），表明 ibh与
因变量臭氧氧浓度具有某种非线性关系。

伪判定系数 Ｒ－sq. ( adj) 显 示 了 模 型 的 解 释 能 力 为R-sq.(adj) =0.459，Deviance explained 表明模型可解释的偏差为51.6%。
可以进一步绘预测变量的制残差图

plot(mgam)

由上图我们可以指导，高度过低和过高对臭氧浓度都存在影响。同理可以做出其他指标的单因素分析。
做完单因素分析后我们开始做多因素分析，我们先建立好模型方程

mgam1<-gam(ozon~s(vh)+s(ibh)+s(dpg)+s(vis)+s(temp),data=be,family=poisson(),model=T)
summary(mgam1)

这里可以看到vh指标大于0.05了，接下来我们对模型进行共线性分析

由此可以看到vh和temp3项都超过0.5了，有可能存在共线性，我们把它们剔除掉，重新进行拟合

mgam2<-gam(ozon~s(ibh)+s(dpg)+s(vis),data=be,family=poisson(),model=T)
summary(mgam2)

从上图表明ibh，dpg，vis具有统计学意义，进一步绘图

plot(mgam2,pages=1)

可以通过三维立体图来显示变量的影响，下图可知ibh，dpg可能存在一定的交互

vis.gam(mgam2)

广义可加模型要求平滑函数效应是可加的，如果两变量间具有交互效应，如果ibh和dpg具有交互效应，则可以通过如下拟合

mgam3<-gam(ozon~ti(ibh)+ti(dpg)+ti(ibh,dpg),family=poisson(),data=be)
summary(mgam3)

由上表可知，ibh和dpg存在交互效应，我们可以进一步绘制交互可视化图形

vis.gam(mgam3,c("ibh","dpg"),theta=30,ticktype="detailed",xlab="高度",ylab ="压力") #模型可视化

由上图可知，压力会对高度存在交互影响，在一定的高度下，过大或者过小的压力都会对臭氧浓度造成影响。
参考文献
1.mgcv包解释文件
2.向伟. 广义可加模型在出生缺陷影响因素分析中的应用及R语言实现过程[J]. 中国妇幼保健,v.29(29):4711-4715.