重复测量资料在临床数据中非常普遍,常用重复测量的方差分析进行统计分析,但是经常面临的问题有:

①临床资料又常常含有缺失值,例如采用某新药治疗疾病,分别在治疗前,治疗后1月,治疗后3月测量Y指标,但由于病人依从性等原因,导致治疗3月后缺失几例数据。

②Y不满足正态性、方差齐性,且样本量不是很大。

怎么办?推荐分析神器之一:混合效应模型。本文结合文献,分享基于R语言实现混合效应分析的方法,主要采用nlme包中lme函数。

主要内容:

1.可视化不同组Hb随时间的变化趋势

2.时间作为分类变量,构建混合效应模型

3.时间作为连续变量,构建混合效应模型

4.模型1和模型2对比和选择

5.模型残差检验

 文献分享 

这篇文章是2021年发表在Neuroimage上,影响因子是5.8,作者观察了4个时间点,通过重复测量三个连续性指标,构建混合效应模型研究正常睡眠和睡眠不足对大脑微观结构的影响。《Sleep and sleep deprivation difffferentially alter white matter microstructure: A mixed model design utilising advanced diffffusion modelling》。

混合效应回归分析 混合效应检验_拟合

混合效应回归分析 混合效应检验_数据_02

作者观察的时间点time分别是早上、晚上、第二天早上。观察的指标Y分别是INVF(神经突内体积分数),exMD(神经突外平均扩散率)和exRD(桡神经突外扩散性)。

上述折线图的横轴是时间点,纵轴是边际估计均值,每个时间点上加了误差线,这个图直观的看出两组人群各指标随时间增加的变化。

混合效应回归分析 混合效应检验_拟合_03

混合效应分析结果显示:

随着时间的增加,INVF 和 exRD 有显著的变化(p<0.05),exMD没有显著变化。而exMD在TP1到TP3和TP1到TP4 的时间点中存在显著的交互作用,而INVF 和 exRD没有交互作用,可以参考原文翻译。

Results from the linear mixed model revealed significant main effects of time in INVF and exRD, but not exMD (see Table 2). Further, there were significant group × time interaction effects in exMD from TP1 to TP3 and from TP1 to TP4, but not INVF and exRD. 

 加载R包和数据 

本案例数据来自外部数据集,共计22名患者,分为组1和组2,测量的指标是血红蛋白浓度Hb,测量的时间点分别是t1,t2,t3,t4。数据概况如下表:

数据结构:自变量X是分组变量,Y指标是4个时间点重复测量Hb浓度。

研究思路:1:Hb随t(时间)的变化趋势是什么?2:组1和组2相比,Hb随t的变化趋势是否不同?

混合效应回归分析 混合效应检验_开发语言_04

1 加载数据

library(nlme); ##混合效应模型library(ggplot2)##作图library(emmeans)##边际估计dt <- read.csv('dt.csv')col = c("#006400","#ef1eef","purple")  ##制作一个颜色变量

注:首先需要将数据集转化为混合效应模型可分析的结构(俗称“宽数据”变为“长数据”),这部分需要先处理数据,关于数据宽变长,后面有机会再给大家整下。这里加载进来的数据是已经处理好。

2 可视化Hb随时间的变化趋势

下面是图一的代码ggplot代码实现ggplot(data=dt,       aes(time,Hb,color=group,group = as.factor(id)))+  geom_point(size=1,fill="White")+  geom_line(          linetype=1,          alpha=0.7,          size=1)+  scale_colour_manual(values=col)+  scale_fill_manual(values=col)+  xlab("time")+  ylab("Hb")+  ylim(0,80)+  theme_classic()+  theme(legend.position="top")

混合效应回归分析 混合效应检验_r语言_05

混合效应回归分析 混合效应检验_混合效应回归分析_06

结果解读:

图一是每个观察对象的Hb随时间的变化趋势,用两种颜色区别两组。图中每条线表示一个患者,纵坐标是Hb,横坐标是时间点。

图二是两组的Hb的估计边际均值随时间变化趋势。横坐标是时间点,纵坐标是估计边际均值,其中这个”均值”跟普通均值稍微有点差异,可以简单理解为均值。每个点上的误差线表示估计均值的标准误SE。

从以上两个图可以看出:Hb随时间呈线性上升的趋势,组2上升斜率低于组1,红色的线在绿色线下方。

3 时间作为分类变量,考察时间点和分组的交互效应

fit1<-lme(Hb~factor(time)*group,random=list(id=~1),          method="ML",data=dt)summary(fit1)

lme的参数说明:

formula:混合效应模型表达式,类似线性回归和二元逻辑回归,*表示分别考虑time,group,以及他们的交互作用。

random:设置随机效应,这里的随机效应是不同患者的随机截距效应。

method:模型参数的估计方法,有两种REML和ML。这里选择ML。

data:是数据集

4 模型1结果解读

混合效应回归分析 混合效应检验_r语言_07

模型1结果解释:

fit1返回的结果很多,

第一:最上面是模型的AIC,BIC和极大似然值;

第二:模型的随机效应,随机效应用方差表示;14.62是随机截距的误差项,6.5是残差

第三:模型的固定效应,也是我们最关注的核心分析结果。

第二列Value是回归系数,代表的是各组各时间点的均值。下面依次罗列。

Intercept:     表示组1的t1时间点的Hb的均值,是26.56;

factor(time)2,表示组1的 t2相比于t1,Hb的差值,是10.19;

factor(time)3,表示组1的 t3相比于t1,Hb的差值,是13.44;

factor(time)4,表示组1的 t4相比于t1,Hb的差值,是15.95;

group2,        表示t1时刻,组2相比于组1的Hb的差值 是-0.43;

factor(time)2:group2,交互项,表示组2的t2与t1两时点差与组1相比,差是-3.58。

factor(time)3:group2,交互项,表示组2的t3与t1两时点差与组1相比,差是-1.61。

factor(time)4:group2,交互项,表示组2的t4与t1两时点差与组1相比,差是-1.06。

根据以上的模型显示:随着时间的增加,Hb逐步增加,而组2比组1增长慢。根据模型的P值,可知,Hb随时间增加,其变化是显著的。而由交互检验P>0.05,显示组1和组2随时间的变化趋势基本一致,没有统计学差异,二者没有交互作用。

5 时间作为连续变量,考察时间点和分组的交互效应

fit2<-lme(Hb~time*factor(group),random=list(id=~1),          method="ML",data=dt)summary(fit2)

6 模型2的结果解读

混合效应回归分析 混合效应检验_拟合_08

模型2结果解读

第一:同上;

第二:同上;

第三:模型的固定效应,也是我们最关注的核心分析结果。

第二列 Value是回归系数:

Intercept:  表示组1随时间拟合直线的截距,表示组1的t1时间点Hb的均值,是23.68;

time,表示组1随时间拟合直线的斜率,表示组1的t2相比于t1,t3相比于t2的Hb的差值,是5.11;

factor(group)2,表示t1时刻,组2相比于组1的Hb的差值 是-1.68;

time:factor(group)2,交互项,表示组2与组1相比,两条拟合直线斜率的差是-0.12。

根据以上的模型显示:随时间增加,Hb变化是显著的。而由交互检验显示,P>0.05,显示组1和组2随时间的变化趋势基本一致,二者没有交互作用。

结果和模型1完全一致。

7 模型1和模型2对比和选择

两个模型似然比检验anova(fit1,fit2)

混合效应回归分析 混合效应检验_r语言_09

结果解释:

模型对比:

模型1和模型2的分析结果完全一致,而模型2更加简洁的给出了两组的Hb随t变化的速度和差别,不过模型2是假使Hb随t呈线性变化,关于选择是否恰当,可以采用似然比检验模型1和模型2,如上,P=0.3941,差异不显著,说明简化的模型2可以代替模型1,本次我们选择模型2作为最终的模型,也就意味这Hb随时间变化是线性的。

如果差异显著,P<0.05,则说明不能用用直线关系代替,可以考虑引入t的平方或者三次方拟合。

模型最终结果分析:

以上的分析已经回答出一开始的两个问题:

1.Hb随着时间增加而直线上升;

2.随着时间的增加,组1的Hb增长率是5.11个单位,组2的增长率是4.99(计算5.11-0.12),但组1和组2增长的差别不显著(P=0.9282),即二者增长率一致。

8 模型残差检验

对模型2进行残差检验Shapiro-Francia法检验sf.test(fit2$residuals[,"id"])  P=0.1128Kolmogorov-Smirnov法检验正态lillie.test(fit2$residuals[,"id"])P=0.3038可视化残差分布hist(fit2$residuals[,1], plot=T, prob=T,       main="Histogram of residuals",xlab="Residuals")lines(density(fit2$residuals[,1],na.rm=TRUE))plot(fit2)

混合效应回归分析 混合效应检验_拟合_10

混合效应回归分析 混合效应检验_拟合_11

结果解释:

结果显示两种正态检验方法的P均大于0.05,说明模型的残差符合正态分布,图1是残差随机分布在直线两边,直方图显示的残差成对称的左右分布,结果均说明,模型残差符合正态分布,模型的构建是合理的。