利用重抽样获取mgcv包的广义可加模型函数曲线的可信区间（3）

自采样目前广泛应用与统计学中，其原理很简单就是通过自身原始数据抽取一定量的样本（也就是取子集），通过对抽取的样本进行统计学分析，然后继续重新抽取样本进行分析，不断的重复这一过程N（大于500次以上）次，然后得到N个统计结果，然后进行区间分析，得到最终结果。

上一章我们介绍了BOOT重抽样获取C指数的可信区间，今天我们来介绍利用sample函数重抽样获取广义可加模型函数曲线的可信区间，这可是一个非常实用的技能，假设我们想了解某连续变量和结果之间的关系，可以使用mgcv包获得两者之间的曲线关系，但是mgcv不能做出95%可信区间，我们可以通过重抽样获取其可信区间。这可是很多的付费课程，付费软件的功能哦。

废话不多说了，我们先导入R包和数据

library(mgcv)
library(data.table)
library(ggplot2)
ac2<-read.csv("E:/r/test/xiyan.csv",sep=',',header=TRUE)

这是一个很简单的数据（公众号回复：吸烟数据，可以获得数据），age_w1是吸烟患者的年龄，cursmoke为患者是否吸烟。假设我们想了解年龄和吸烟状态存在的非线性关联。
建立模型

mgam.lr1<-gam(cursmoke~s(age_w1,k=3),family = binomial(link = "logit"),
               data = ac2)
summary(mgam.lr1)

P小于0.05表明存在非线性关系，绘图

plot(mgam.lr1,se=T)

该图表示年龄小的时候，吸烟状态变化很小，超过60岁以后变化很快，这样看不是很好理解，我们把它改成吸烟概率更加好理解，我们新建一个数据，把结果装进来

newdat<-as.data.table(expand.grid(
  age_w1=seq(
    from=min(ac2$age_w1,na.rm = T),
    to=max(ac2$age_w1,na.rm = T),
length.out = 1000)))

新建好数据后，在新数据上生成预测概率

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

绘图

ggplot(newdat,aes(age_w1,yhat))+
  geom_line()+
  xlab("year")+
  ylab("吸烟概率")

上图可知，随着患者年龄增加，吸烟概率降低，有不少粉丝问我，怎么把Y轴改成百分比显示，下面我们来改一下

ggplot(newdat,aes(age_w1,yhat))+
  geom_line()+
  scale_y_continuous(labels = scales::percent)+
  xlab("year")+
  ylab("吸烟概率")+
  coord_cartesian(xlim = range(newdat$age_w1),
                  ylim = c(0,.4),
                  expand = F)

OK，已经改好了，60岁的时候30%的吸烟概率，到了80岁只有10%了，接下来我们进行重抽样获取可信区间，我们打算重抽样500次，先建立一个空矩阵

nboot<-500
out<-matrix(NA_real_,ncol = nboot,nrow =nrow(newdat))

建立好以后就可以重抽样了，为了可重复，要设一个种子

set.seed(12345)

建立循环重抽样

这个过程有点费电脑，时间根据你的电脑速度来定，最近读了张敬信老师的向量化编程，深有启发，到时试试能不能用向量化编程来实现。

数据很大，上图只是一部分，求出数据后，剩下的就简单多了，我们看一下我们采样值和系统预测值是否存在偏差

mean(abs(newdat$yhat-rowMeans(out)))

由此可知偏差非常小，接下来就是生成可信区间数据

最后绘图

OK，到此全部完成，但是后台很多人私信问我BOOT重抽样怎么进行内部验证，看了今天内容应该就会了吧。留个问题给大家想想，假如我们用BOOT重抽样，上面的图应该怎么做？
原创不易，需要文章全部代码的朋友，请把本文章转发朋友圈集10个赞，截图发给我，嫌麻烦的也可以在这里下载：https://download.csdn.net/download/dege857/85404235

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