广义可加模型 Python 主要模型参数 广义可加混合模型

收到粉丝投稿要求我做个这样的图，来自于2017年发表在Eur Urol期刊（SCI IF=17.9分）。The Impact of Local Treatment on Overall Survival in Patients with Metastatic Prostate Cancer on Diagnosis: A National Cancer Data Base Analysis. Eur Urol, 2017.

目前常用的曲线拟合主要是广义可加模型和RCS，既往我已经介绍了《利用重抽样获取mgcv包的广义可加模型函数曲线的可信区间（3）》，上图其实就是在这个的基础上进行分类的曲线拟合。

说的是什么意思呢？对于上图来说就是，男女就是两个变量对于随着年龄变化他们吸烟率的变化，我这张图两张曲线没有交点，第一张图有，表示横轴X演变到了一个点，不同类别的变量，概率发生了变化，原来红线高的，后面变低了。下面来演示一下，使用一个吸烟数据（公众号回复：吸烟数据2，可以获得数据），我们先把数据导入

ac1<-read.csv("E:/r/test/smoke2.csv",sep=',',header=TRUE)
library(VGAM)
library(data.table)
library(ggplot2)

这个数据很简单，其实就是上次的吸烟数据中加了性别，smoke_w1为吸烟状态，结局变量，sex是性别，age_w1是患者的年龄。

我们想了解不同性别随着年龄变化吸烟的情况，先把smoke_w1转成二分类的，就是是否吸烟

ac1$cursmoke<-as.integer(ac1$smoke_w1=="(1) Cur Smok")

现在可以正式分析了，我用的是vgam包来拟合广义可加模型，gam也是可以的，做法也是相同的

mgam.lr1<-vgam(cursmoke~sex+s(age_w1,df=3),family = binomialff(link = "logitlink"),
               data = ac1,model = T)

解析模型，P小于0.05

生成预测模型数据

在预测模型数据上生成概率

newdat$yhat<-predict(mgam.lr1,newdata=newdat,type="response")

生成了预测概率后就可以绘图了

ggplot(newdat,aes(age_w1,yhat,col=sex))+
  geom_line(size=2)+
  scale_color_viridis(discrete=T)+
  xlab("year")+
  ylab("吸烟概率")+
  coord_cartesian(xlim = range(ac1$age_w1),
                  ylim = c(0,0.5),
                  expand = F)

美化一下就可以用于发表了

我们还可以在这个的基础上进行重抽样获取它的可信区间，可以参考文章《利用重抽样获取mgcv包的广义可加模型函数曲线的可信区间（3）》的做法。如果两条曲线相交了，我们怎么找它的交点呢？通过预测值列表哪里找就可以了。

OK，本章结束觉得有用的话多多分享哟。