R语言实现广义加模型

2020-04-14 15:02:19 浏览数 (2)

今天给大家介绍一个广义加模型(Generalized Additive Model),其是广义线性模型的扩展,其在线性联系函数的基础上增加了一个平滑函数。我们首先看下这个广义线性模型的定义:

首先是确定一个数据集属于什么分布,如正态分布,二项分布等,然后是线性组合,通过添加参数使得数据集满足一个线性方程。最后就是创建连接函数,也就是把我们的数据集生成一个线性的方程,如普通线性模型、对数线性模型等。广义加模型主要是通过对自变量引入平滑函数,降低线性设定带来的模型风险。接下来,我们看下实现模型需要的包:

代码语言:javascript复制
install.packages("mgcv")
install.packages("gamRR")

包安装好后,我们看下具体的函数,在mgcv中有两个函数都可以实现gam,首先我们看下gam函数:

其中主要的参数:

Formula 模型的形式。其中和GLM一样可以对各自变量平滑处理的几个形式:

s(x1,x2,...,k=12,fx=FALSE,bs="tp",by=z,id=1)。其中bs主要涉及以下几种:tp – DEFAULT, thin plateregression 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 smoothinteraction;K代表的自由度,也就是需要处理的次数;fix指的是不是需要增加处罚。

te(x,z,bs=c("tp","tp"),m=c(2,3),k=c(5,10))对多个变量的平滑处理。

另外ti,t2这些都是遵循splines的参数列表。

Family 指的数据的分布形式。

Data 在这里data需要以列表形式的数据格式输入。

Method 平滑参数估计方法。“GCV.Cp”用于未知的标度参数,而Mallows’Cp/UBRE/AIC用于已知的标度。也可以是“GACV.Cp”,就是用GACV代替GCV,但是等价于“GCV.Cp”。“REML”用于REML估计,包括未知的规模,“P-REML”用于REML估计,但使用的是Pearson估计的规模。“ML”和“P-ML”相似,但使用极大似然代替REML。除了指数族之外,“REML”是默认值,唯一的其他选项是“ML”。

其它的参数我们不去赘述了。

那么我们直接看实例:

代码语言:javascript复制
library(mgcv)
library(gamRR)
dat <- gamSim(1,100,dist="poisson",scale=.25)#为gam模拟样本数据
fit <- gam(y~s(x0) s(x1) s(x2) s(x3),family=poisson,dat,method="REML")
plot(fit,select=2)#通过select选择要绘制的变量

我们还可以利用gam.check(fit,pch=19)来看下我们模型的具体情况,也可以说做一个评估:

结果中我们可以看出通过k-index来判断是否K值太低。K-index小于1则代表k太低,同时如果edf越接近k那么越低。

此外,为了方便大数据量的计算,还引入了bam,其优点是内存占用比gam低得多,但是对于大型数据集,它也可以快得多。bam也可以在集群上通过调用parallel 包进行计算。

Bam和gam唯一的区别就是相当于bam对gam进行了并行化处理。

最后我们,还要引入另外一个包来计算广义加模型的相对风险比(RR),这个值在临床中是很常见的主要用来描述队列研究中分析暴露因素与发病的关联程度。这里呢,主要是用来描述针对某一个变量此模型的风险性,评估函数是gamRR(fit,ref,est,data,n.points,plot,ylim)。

我们直接看下实例:

代码语言:javascript复制
gamRR(
fit=fit,
ref=c(x0=dat$x0[1],x1=dat$x1[1],x2=dat$x2[1],x3=dat$x3[1]),
est="x1",
data=dat,
n.points=10,
 plot=TRUE,
ylim=NULL)

至此就模型的整个过程就结束了。

欢迎大家学习交流!

0 人点赞