当解释变量与效应变量间关系不明确时,通常可以使用广义相加模型来检测比变量间是否具有非线性关系。

广义相加模型通过光滑样条函数 、核函数或者局部回归光滑函数,对变量进行拟合。GAM采用模型中的每个预测变量并将其分成多个部分(由'结'​​分隔),然后将多项式函数分别拟合到每个部分。GAM的原理是最小化残差(拟合优度)同时最大化简约性(最低可能自由度)。回归模型中部分或全部的自变量采用平滑函数,降低线性设定带来的模型风险,对模型的假定不严,如不需要假定自变量线性相关于因变量(线性或非线性都可以),并且解决logistic回归当解释变量个数较多时容易引起维度灾难(Curse of dimensionality)。

R中的实现主要有两个包:gam和mgcv,前者为旧版本的包,后者为新版本。两个包的基本建模过程是相似的(两个函数都是gam ;要小心同时加载两个库),但幕后计算方法不同,优化和模型的参数也不同。当在同一数据集上使用相同的模型结构时,预计结果会略有不同。

具体解释可参见:https://blog.csdn.net/textboy/article/details/47277131

这里以mgcv包为例

1.构建基础模型

gam(y~s(x,k = , bs =)) / gam(y~te(x,k = , bs =))

k: sets up the dimensionality of the smoothing matrix for each term. Penalized regression smoothers. Using a substantially increased k to see if there is pattern in the residuals that could potentially be explained by increasing k. Default任意数字(normally 10 degree of freedom)。

bs: See smooth terms for the full list. tp – DEFAULT, thin plate regression spline,cr – penalized cubic regression spline三次样条, cs – shrinkage version of cr,cc – cyclic cubic regression spline, ps – P-spline,cp – cyclic p-spline, ad – adaptive smoothing, fs – factor smooth interaction.

s: smooth s(covariate, edf); te: tensor product smooth
 

install.packages('mgcv')
library(mgcv)
mode<-gam(y~s(x,df),family="",data)

2. 在基础模型构建后,最重要的是确定模型自由度。确定自由度有几种方式:

1)基于生物学知识和专家经验;2)赤池信息准则(AIC),根据AIC最小确定自由度;3)根据残差独立原则,最小化残差自相关确定自由度;4)广义交叉验证

以最小化残差自相关为例,可构建循环,设定自由度为i(从1到20),分别构建20个相应自由度的模型并计算出每个模型残差偏自相关绝对值的和stat,通过作图的方式将i与stat之间的关系显示出来。以下为示例代码:

stat<-NULL
mod<-list()
for (i in 1:20){
mod[[i]]<-gam(y~s(time,df=2*i)+as.factor(dow),family=possion,data=mydata)
tt<-sum(abs((pacf(mod[[i]]$residuals))$acf))
stat<-append(stat,tt)
}

上述代码中,列表对象mod中包含20个GAM模型,tt为每个模型残差偏自相关绝对值的和,stat则是由20个GAM模型的tt值构成的向量。

3.构建模型

final_mod<-gam(y~so2+s(time,df=14*2)+as.factor(dow),family=possion,data=mydata)
summary(final_mod)

4.模型比较与评价

在R中,GAM的比较与评价可以用aov(model1,model2,test="Chisq")或者AIC(model1,model2)来比较。通过比较选出最优模型后,可通过观察模型中非参数平滑函数的自由度改变对解释变量的影响大小来评判模型是否稳健,又称敏感性分析。

5.模型参数的解释

以下面例子为例:

library(mgcv)      #加载mgcv软件包,因为gam函数在这个包里
Data <- read.delim("Rice_insect.txt")     #读取txt数据,存到Data变量中
Data <- as.matrix(Data)     #转为矩阵形式
#查看Data数据:Data,查看第2列:Data[,2],第2行:Data[2,]<br>
Adult<-Data[,2]
Day<-Data[,3]
Precipitation<-Data[,4]<br>
result1 <- gam(log(Adult) ~ s(Day))     #此时,Adult为相应变量,Day为解释变量
summary(result1) 

结果为

Day的影响水平p-value=0.473,解释能力为14.3%,说明影响不明显。

其中,edf为自由度,理论上,当自由度接近1时,表示是线性关系;当自由度比1大,则表示为曲线关系。

result2 <- gam(log(Adult) ~ s(Precipitation))
summary(result2)

此时p-value为0.0774,说明该因子在P<0.1水平下影响显著。

参考文献:[1]张云权, 朱耀辉, 李存禄,等. 广义相加模型在R软件中的实现[J]. 中国卫生统计, 2015, 32(006):1073-1075.

[2]广义相加模型及其R实现 - 知乎 (zhihu.com)

Logo

旨在为数千万中国开发者提供一个无缝且高效的云端环境,以支持学习、使用和贡献开源项目。

更多推荐