广 义 相 加 模 型 ( generalized additional model，GAM)是对传统广义线性模型的非参数拓展，可有效处理解释变量与效应变量间复杂的非线性关系。GAM 目前已广泛应用于空气污染流行病学研究中，主要用于分析空气污染或气象因素对人群健康事件(如发病、住院和死亡)的急性损害效应。
在基础模型构建之后，最重要的工作就是确定模型中非参数平滑函数的自由度 df。在广义相加模型中，由于平滑函数的自由度对模型的参数估计和模型稳定性有一定影响。因而，选择合适的自由度对模型构建有重大意义，通常根据以下评判准则进行设定:
①基于生物学知识和专家经验(包括敏感性分析)设置固定的自由度;
②赤池信息准则 (Akaike information criterion，AIC)，依据 AIC 最小选择自由度。
③依据残差独立原则，通过最小化模型残差自相关来选择自由度。实际工作中，我们根据基础模型残差的偏自相关(PACF)绝对值之和最小选取自由度。
④依据广义交叉验证 (generalized cross-valida-tion，GCV)预测污染物浓度的最佳模型(GCV-PM 10 )选择自由度，这种方法是最小化误差均方过程的一种简化。
今天我们来演示②赤池信息准则和③依据残差独立原则确定模型自由度，继续使用我们的美国芝加哥1987年至 2000年大气污染与死亡数据（公众号回复：芝加哥2，可以获得数据）做实例分析，我们先导入需要的R包和数据看看
library(foreign)
library(mgcv)
library(dlnm)
library(splines)
library(car)
be read.spss("E:/r/test/ozone.sav",
use.value.labels=F, to.data.frame=T)
names(be)
我们先来看看数据的构成，death：死亡人数 (per day)，pm10：大气污染物pm10的中位数值，pm25median，o3median：二氧化硫的中位数值，time：天数，这里就是我们的时间，tmpd：华氏温度，date：日期，dow为星期几
我们等会要用到dow这个变量，我们先对数据进行转换，把它变为周末和非周末，这样就把它变成二分类变量了
bc$dowrecode(bc$dow,"'Saturday'=1;
'Sunday'=1;else=0")
bc$dowfactor(bc$dow)
整理好数据后我们就可以开始分析了，我们先用赤池信息准则，也就是AIC来分析，我们要求出最小AIC的模型
假设我们的模型为
fitgam(death ~ ns(time, 7*14) + dow,family=poisson(), bc)
自由度df=7
我们想知道选7不是不是最适合模型的，我们可以通过查看AIC来判定
extractAIC(fit)
由此可知次模型的AIC值为40424.65，最笨的办法就是把1:20个自由度分别做出模型，然后进行比较AIC，然后取最小值AIC的模型。然而我们也可以通过一个小小的循环做出来
我们先建立3个空值，等下要用到
tmpNULL
ttNULL
statNULL
建立一个循环（如下），这个循环比较简单，不容易出错，第一句为方程，第二句为求AIC，第三局把AIC提取出来，第四局为把结果打印出来，但是在family这里我们要选择泊松分析，类泊松分析求不出AIC
for (i in 1:20) {
fitgam(death ~ ns(time, i*14)+dow,
family=poisson(), bc)
tmpextractAIC(fit)
tt(tmp[2])
print(tt)
}
在上面出来的AIC值中选一个最小的就可以了，如果我们的自由度选择很多，数据太多不好分辨，我们改为如下
tmpNULL
ttNULL
statNULL
for (i in 1:20) {
fitgam(death ~ ns(time, i*14)+dow,
family=poisson(), bc)
tmpextractAIC(fit)
tt(tmp[2])
statappend(stat,tt)
}
新生成的stat变量里包含了所有的AIC值，我们对它取最小值就可以了
stat
min(stat)
由此可知道，根据AIC法，自由度取17时AIC最小。
下面我们来介绍一下，
基础模型残差的偏自相关(PACF)绝对值之和最小选取自由度
主要就是求出模型残差的和
先建立3个空值
tmpNULL
ttNULL
statNULL
建立循环
for(i in 1:20){
fitgam(death ~ ns(time, i*14) + dow,
family=poisson(), bc)
ttsum(abs((pacf(fit$residuals))$acf))
statappend(stat,tt)
}
跑20个数据，需要等一点时间，根据各自电脑快慢速度时间也不同，结果最终也是存储在stat变量里面
stat
min(stat)
根据此法算出来df应该取6，根据目前专家经验 和学术论文报道时间的自由度应该是7左右，因此基础模型残差的偏自相关(PACF)绝对值之和算出来的自由度对于这种广义可加模型模型拟合天气污染的应用还是比较准确的。此过程中还生成了相关的残差图片，我就不一一放上来了
参考文献
mgcv包解释文件
tsModel包解释文件
杨敏娟, 潘小川. 北京市大气污染与居民心脑血管疾病死亡的时间序列分析[J]. 环境与健康杂志, 2008(4):4.
邓建明, 秦伯强, 王博雯. 广义可加模型在R中的快捷实现及蓝藻水华预测分析中的应用[J]. 生态学杂志, 2015, 34(3):8.
向伟. 广义可加模型在出生缺陷影响因素分析中的应用及R语言实现过程[J]. 中国妇幼保健, 2014, 29(29):5.
张云权, 朱耀辉, 李存禄,等. 广义相加模型在R软件中的实现[J]. 中国卫生统计, 2015, 32(6):3.
殷文军, 彭晓武, 宋世震,等. 广州市空气污染与城区居民心脑血管疾病死亡的时间序列分析[J]. 环境与健康杂志, 2012, 29(6):6.
路凤, 李亚伟, 李成橙,等. 时间序列分析在空气污染与健康领域的应用及其R软件实现[J]. 中国卫生统计, 2018, 35(4):4.
Original: https://blog.csdn.net/dege857/article/details/121760618
Author: 天桥下的卖艺者
Title: R语言mgcv包时间序列分析在空气污染与健康领域的应用（3）—模型自由度选择
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/697896/
转载文章受原作者版权保护。转载请注明原作者出处！