广义相加混合效应模型,是混合效应和相加模型的结合，不仅可以引进随机效应，还可以对重复测量的X(自变量)、t及其他协变量使用曲线拟合，可以满足上述的分析要求。本文主要采用mgcv包中gamm函数分析。

文献分享

文献1

这篇文章是基于MIMIC-III数据库2021年5月份发表在Frontiers in Medicine（IF=5.091分）上的文献。收集ARDS患者入院后7天内重复测量的NLR值，结局变量：30天死亡率和院内死亡率；作者①首先采用多因素Logistic回归，并拟合光滑曲线探索基线NLR与结局变量之间的相关性；②进一步采用广义相加混合模型（GAMM）用于分析死亡者和存活者之间NLR差异随时间的早期变化。The Association Between the Baseline and the Change inNeutrophil-to-Lymphocyte Ratio and Short-Term Mortality in Patients With Acute Respiratory Distress Syndrome

image.png

image.png

上图是：拟合平滑曲线探索基线NLR水平与结局变量之间的关系的结果，关于这个图是分析和绘制。

image.png

image.png

上图是：采用GAMM模型分析存活者与死亡者之间的NLR差异随时间（入ICU后第一周）的变化，图片中红色表示死亡组，蓝色表示存活组，显示两组患者的NLR趋势明显不同，两组之间的差异随时间而增加。

文献2

这篇文章是2021年发表European Radiology的一个关于新冠的研究，新冠患者按症状分成重症组和轻微组，入院3到30天内，采用CT多次重复扫描肺部肺部感染量和感染总量比，观察两者随时间的变化和影响。Generalized additive mixed model to evaluate the association between total pulmonary infection volume and volume ratio, and clinical types, in patients with COVID-19 pneumonia: a propensity score analysis.

image.png

作者采用线性混合效应模型研究时间和症状对感染量和感染总量比的影响，结果见表四。同时采用GAMM模型研究感染量和总量比随时间的变化趋势，如上图所示。红色是重症组，蓝色是轻微组，结果：

图a是总肺部感染量随时间的影响，同时受临床类型的显着影响（交互作用 p =0.01）。轻微组随着时间增加而减少，而重症组患者肺部感染总量，随着时间增加每天增加14.66 cm3（95% CI：3.92 至 25.40）。

图b 肺部感染总量比随时间的影响，同时受到临床类型的显着影响（p表示相互作用= 0.01）。轻微组随着时间增加而减少，而重症组肺部感染总量比，随着时间增加每天增加 0.45% (95% CI: 0.13 to 0.77).

本案例数据来自外部数据集，根据5个城市每年的某传染病感染率与癌症死亡率的数据，分析这两者之间的关系，其中Infec_Disease表示某传染病感染率，Cancer_Mortality表示癌症死亡率，具体数据概况如下表：

数据特点：①某传染病感染率与癌症死亡率，均是重复测量的数据集，都随时间变化；②各个城市提供的数据年份不完全相同，如城市3的资料年份是1969-2008，共计32个数据；城市5的资料的年份是1955-2008，共计38个数据。

分析思路：1传染病感染率(X)和肿瘤死亡率(Y)分别随t(时间)的变化趋势是什么？(请看上篇推文)2.传染病感染率(X)和肿瘤死亡率(Y)有没有联系?（本篇解决）3.有没有滞后效应，有没有非线性关系？

在分析X与Y有没有联系的时候，需要考虑以下两个问题：

①X,Y都随时间变化，因此在分析X与Y的关系时，要调整t的影响，如何调整呢？调整的不好，t的混杂因素作用并没有得到有效控制，得出来的结果就不可靠。②Y随X变化是否是线性变化关系？有没有分段或者阈值效应呢？如果不是线性变化，而是直接用直线拟合，就不能发现他们之间得关系。

接下来带着这两个问题，往下走：

image.png

1数据加载

library(mgcv)  gamm
library(dplyr) 数据转换
dt <- read.csv('cancer.csv')
col<- c("#FF5151",'blue')##制作一个颜色变量

2GAMM模型构建

首先采用曲线拟合传染病感染率(X)与t，回答分析思路的第二个问题。
fml <- "Cancer_Mortality~s(Infec_Disease,bs='cr')+s(Year,bs='cr')"
fit<-gamm(formula(fml),random=list(Country=~1),data=dt,family=gaussian)
fit

gamm的参数：

formula：GAMM模型表达式，s()表示是平滑项，平滑的方法可以通过bs参数自定义，其中“tp”代表薄板回归样条，“cr”代表三次回归样条。
random：设置随机效应，这里的随机效应是不同患者的随机截距效应。
family：链接函数，同glm函数，连续变量写gaussian，二分类写binomial
data：数据集，纵向数据。
其他参见帮助文档。

3绘制平滑曲线

绘制平滑曲线回答上述分析思路的第二个问题，传染病感染率(X)和肿瘤死亡率(Y)有没有联系?下面代码是调整曲线拟合t后，绘制X与Y曲线图，而调整曲线X，绘制t与Y的曲线图。

提取拟合值
pred<-predict.gam(fit$gam,type="terms",se.fit=TRUE)
mfit<-mfit+mean(fit$gam$fitted.values)-mean(mfit)d
at<-cbind(dt,mfit,sfit)rm(mfit,sfit)
计算拟合值95%CI置信区间
dat$y.low <- dat$mfit-1.96*dat$sfit
dat$y.upp <- dat$mfit+1.96*dat$sfit
attach(dat)
绘制曲线
co<-c(0.003, 19.2415, 55.15824 ,77.0227)  #前两个是x轴的范围，后两个y轴的范围。
plot(mfit~Infec_Disease,ylim=c(co[3],co[4]),xlim=c(co[1],co[2]),col=col[1],type="p", pch=20, ylab="", xlab="")
par(new=TRUE); 
plot(y.low~Infec_Disease,ylim=c(co[3],co[4]),xlim=c(co[1],co[2]),col=col[2], type="p", pch=16, lwd=2,ylab="", xlab="",cex=0.5)
par(new=TRUE);
plot(y.upp~Infec_Disease,ylim=c(co[3],co[4]),xlim=c(co[1],co[2]),col=col[2], type="p", pch=16, ylab='Cancer_Mortality', xlab='Infec_Disease',cex=0.5)

image.png

image.png

结果解读：

由上图一可以看出，调整曲线t后，可以看到X与Y成线性关系，横轴是Infec_Disease(X),纵轴是Cancer_Mortality(Y),随着Infec_Disease的增加，Cancer_Mortality呈线性下降。

由上图二可以看出，调整曲线X后，可以看到t与Y成非线性关系，横轴是t,纵轴是Cancer_Mortality(Y),随着t的增加，Cancer_Mortality呈非线性上升。

因此说明要分析Infec_Disease和Cancer_Mortality的关系时，需要调整t的影响，且要用非线性关系调整t的混杂作用。

接下来，我们采用直线分别拟Infec_Disease，曲线拟合t，重新构建GAMM模型。

4直线拟合GAMM

根据上述的分析结果，X与Y成线性，t与Y成曲线，因此我们采用直线分别拟Infec_Disease，曲线拟合t，重新构建GAMM模型，探索X与Y的关系。

fit<-gamm(formula(fml),random=list(Country=~1),
data=dt,family=gaussian)
summary(fit$gam) 查看相加部分结果
anova(fit$gam)  查看方差检验结果
summary(fit$lme)查看混合效应部分结果

结果解读：

一开始提到广义相加混合效应模型是混合模型和相加模型的结合，fit这个对象结果中同时包含gam和lme部分的结果，并且直线拟合传染病感染率(X)与癌症死亡率(Y)时，gam和lme的结果一致,主要看相加模型的线性回归项和混合效应模型的固定效应部分，这是回归的主要结果。

结果显示：
曲线调整t的混杂作用后，Y随X成反向变化，回归系数是-0.6437，传染病感染率每增加一个单位，癌症死亡率下降0.6437，那么如果用直线拟合t就不一定能充分调整t的混杂作用，就可能会影响对X作用的准确评估。这个地方还可以采用线性混合效应模型拟合t^2/t3调整t的作用。

其他关于R返回的结果见上述红框标注，详细可以参考前篇推文的介绍