记录一下R做GAM画图（smooth spline）的一些情况

画图（smooth spline）的一些情况"/>

记录一下R做GAM画图（smooth spline）的一些情况

临时又有工作需要画GAM的图，网上的资料比较零散，尤其是本站的GAM的教程，多是从统计的角度去讲、少有专注于R程序实现的不说，一些本来是讲R的，居然还收费了……我真是#@￥%#@￥@#￥#
因此记录一下重新摸索的过程中，发现的一些坑和小经验。希望下次再需要摸索的时候，只看自己写的这玩意儿就够了。

广义相加模型（GAM）的一些简单定义

GAM：general additive model 广义相加模型，按照我浅薄的理解，就是把广义模型里面的某些自变量给设定成曲线项，用了个惩罚b样条（penalized B splines）去自动做非线性关系的拟合而不用人工再去手动一个一个一个尝试（什么是不是2次项啊？是不是3次项啊？先往哪边歪啊blabla……）
然后这玩意儿拟合出来的结果是个啥？能说明什么问题？我认为只能看一个大概的趋势。当然，线性模型一样可以看一个大概的趋势，因此这玩意儿主要也就是瞅一瞅有没有曲线的趋势。真要得到x与y具体是什么样的关联，得到或许能够体现一定问题的那个系数，还需要后面再根据gam的趋势，找个切点去做分段线性回归再去进行验证。

先看我之前最常用的gam图

简单的GAM用plot直接做（带predictive interval）

以截取的部分数据为例

包括5个变量location、year、x、y、corvar（协变量）

首先要拟合构建gam，需要用到mgcv包
顺便把tidyverse也加载了

library(mgcv)
library(tidyverse)

gam函数的参数最主要的前两个：formula和family
formula写明x和y，其中曲线项用s()表示
family指y的分布，这里假设是高斯分布
这个数据中，需要看的是x与y的非线性关系，将corvar作为线性的协变量调整，则代码为：

fit1 <- df %>%gam(formula = y ~ s(x) + corvar, family = "gaussian")plot(fit1)

查看图像得到下图：会自动给出predictive interval（旁边两条虚的曲线）以及在x轴上的样本量稠密度。

但当遇到类似我这次的数据与要求时，上述的图像就不太够用了。
这次的数据，第一个向量location表明，具有不同地区，因此有时候就需要根据地区来分组，既要看每一个组的非线性的趋势、又要看整体的大趋势、还要把它们放到一张图里……
这种时候一般人可能会先想到，我直接用管道函数连上group_by_，再画图行不行？
我试了下，搞不懂怎么用ggplot画出来gam的图，而plot()……好像不能连管道函数，
因此就需要用到下面的方法：

scatter.smooth() 同时需要展现各组趋势与整体趋势的smooth spline

这种时候先拟合、在做图的思路多少有些行不通了，就需要用到scatter.smooth或者是loess.smooth函数。这两个函数都是自带的。
首先尝试用scatter.smooth()来画图：
代码为：

scatter.smooth(x=df$x,y=df$y)

这个东西也不能用管道函数，因此没法直接分组，但神奇的是它会直接把点和线分别显示到图里：

那么剩下的，就是只要调整一下点，把不同location的点用不同的颜色或者形状来显示，并且附上图例，就算完成任务了。
幸运的是，在我搜到scatter.smooth的时候，同时搜到loess.smooth可以用来解决类似的问题。
这种思路就是用在plot()后面with()一个loess.smooth()的line()：
首先尝试一下把点和线分开设置颜色：

plot(y ~ x, data= df,col='green')
with(df,lines(loess.smooth(x,y),col='blue'))

得到图像为：

然后考虑把location的信息体现在图像上。
考虑到plot函数当中具有一些参数，其中col设定颜色、pch设定形状，并且都有对应的数字，
因此下一步需要把location的信息变成数字，并且统一。即新建一个向量numlocation：

mark1 <- df%>%count(location)
mark1$numlocation <- c(1:length(mark1$location))
mark1 <- mark1 %>% select(location,numlocation)
df <- left_join(df, mark1, by="location")

则数据框变成了：

将col和pch都设定为numlocation，另外再调整一下趋势线的粗细

plot(y ~ x, data= df,col=numlocation,pch=numlocation)with(df,lines(loess.smooth(x,y),col='blue',lwd=2))

剩下的就是图片直观、美观之类的工作了
一是图例需要有，然后是图例不能跟图像重合，三是字体
设置图例用legend函数，legend(x,y,legend,…)
其最关键的三个参数，第一个是图例的x轴位置、第二个是图例的y轴位置，第三个是传入图例的数据。

legend(0.9,1000, legend = mark1$location, xpd = TRUE,col = mark1$numlocation, pch = mark1$numlocation, cex = 0.7, bty = 'n')

此外xpd参数用来与par()函数相配合，防止图例盖在图像上面。
col参数和pch参数与图像当中，即plot函数中相对应，调整图例前面的小标记的颜色和形状。
cex参数为图例字体大小，
bty参数为图例外围一圈的框是否显示，默认为‘O’即显示，不显示为’n’。

避免图像与图例重合用par函数, par(mai=(bottom,left,top,right))
用来设定绘图区域，数字越大，相应区域的留白就越多。
与此同时，由于字体是整个图像保持一致，因此也可以在par这里设置。

Windows系统中查看字体用windowsFonts函数

windowsFonts()

显示包含如下字体，

$serif
[1] “TT Times New Roman”
$sans
[1] “TT Arial”
$mono
[1] “TT Courier New”

根据工作要求，选择新罗马，

最后把绘图区域设置、字体设置、图例、轴标题都加上，完整的代码为：

par(mai=c(1,1,1,2),family='serif') # 设定绘图区域和字体
plot(y~ x, data= df,col=numlocation,pch=numlocation,xlab='',ylab='')
#绘制散点图，调整点的形状和颜色，取消默认的x、y轴标题
with(df,lines(loess.smooth(x,y),col='blue',lwd=2)) # 拟合曲线、设定颜色和宽度
title(xlab = 'This is x', ylab = 'This is y')
# 重新写轴标题
legend(0.75,20, legend = mark1$location, xpd = TRUE, col = mark1$numlocation,  #图例的颜色pch = mark1$numlocation,  #图例的形状，与图中保持一致cex = 0.7, bty = 'n') #图例字号、图例外圈的方框变透明

得到的图像为：

ggplot 既要展现分组 还得有predictable interval的时候到了

仍然是以上面的数据为例

df %>% ggplot(aes(x,y))+geom_point()+geom_smooth(method = "gam")

通过aes设定x和y，通过geom_point画散点图、通过geom_smooth拟合曲线，就能得到初步的图：

可以看到曲线旁边的PI是有的，如果想要去掉，在geom_smooth里面把参数se设定为false即可

但此时，一来没有图例，二来原本计划分组展示的变量location，在散点图图像中体现出来的确实不同的location都是一样的黑色圆点。因此需要重新设定一下散点图的形状与颜色。在geom_smooth里传入颜色和形状的参数：

df %>% ggplot(aes(x,y))+geom_point(aes(shape=location,color=location))+geom_smooth(method = "gam")

另外需要注意的是，geom_point里面自带的颜色形状分类只到6个，当分类变量的取值>6时，就会出现分类不够的情况，需要另外+scale来对颜色和分类进行补充说明，例如当location有22个时：

后面的都分类都显示不出来 不只是图例 在散点图中也没有，因此需要补充上scale的代码
由于同时将location传入到shape和color参数，因此也需要分别设定，此时location的分类有22个
同时，将两列图例调整为一列，并改变x和y轴的名称

p1<-dataASMR %>% ggplot(aes(SDI,val))+geom_point(aes(shape=location,color=location))+scale_shape_manual(values = 1:22)+ #增加shape的分类scale_color_manual(values = 1:22)+ #增加color的分类geom_smooth(method = "gam")p1 <- p1+guides(shape=guide_legend(ncol = 1),color=guide_legend(ncol = 1))+labs(x='x',y='yy')p1