纵向线性混合模型,线性混合模型和广义估计方程

背景

线性模型需要满足正态性、独立性和同方差性等假设，其中独立性是线性模型最重要的假设之一，独立性要求每一个数据点必须来自于不同的总体。但由于重复测量数据、区组数据以及空间相关数据不能满足独立性假设，因此常常利用线性混合效应模型对上述数据进行分析（如相关性分析）。

使用一般线性模型时，是需要满足以下3点假设的

正态性，因变量y符合正态分布独立性，不同类别y的观察值之间相互独立，相关系数为零方差齐性，不同类别y的方差相等

线性混合效应模型的例子

A. 假设你想研究一个小镇上（4个街区）的平均收入和教育水平的关系，不考虑街区区分选取1000个体样本。如果这样构建模型，忽视了样本之间的相关性，同一个街区的个体不是独立的，因此残差项会于街区关联。B. 假设你想研究焦虑y和尿酸，三酸甘油酯的水平，有1000个体，在24小时内测量了10次。如果这样构建模型，y的方差随着时间而变化，从而导致方差异质性，这可以通过残差方差异质性体现出来。在A，由于空间关联性导致样本非独立；在B，存在异质方差性。这两种情况都不适于使用常规的线性模型，但是都可以用混合线性模型。

以一个探索基因表达量与采样距离相关性为例。
该数据来自于“ntegrated transcriptomic analysis of distance-related field cancerization in rectal cancer patients ”。该研究的思路为，取距离肿瘤原发灶不同距离处样本，进行RNAseq测序，找到表达量与距肿瘤远近相关的基因，然后进行功能注释、生存分析等等。

本文以上述数据为依托，用lme4构建线性混合效应模型，找到表达量与距离相关的基因。

为什么采用线性混合效应模型
我们希望找到表达量与距离相关的基因，但因为来自同一个患者的表达量是相关的，所以这128个样本不满足普通线性回归的独立性假设，需要用到线性混合效应模型。

数据准备 # 下载GEO数据library(GEOquery)gset <- getGEO('GSE90627', destdir="./raw-data/", AnnotGPL = F, ## 注释文件 getGPL = F) # 下载注释文件gpl <- getGEO("GPL17077",destdir="./raw-data/")ids <- Table(gpl)# 获取表达矩阵a=gset[[1]]expr=exprs(a) # 探针id转换成gene symbolrownames(expr) <- ids[match(rownames(expr),ids$ID),7] # 获取表型信息pd=pData(a) pd <- pd[,c(1,2,37:40)]colnames(pd)[3:6] <- c("age","gender","patient","tissue")# 构建出patient和distance变量library(stringr)pd$distance <- str_split(pd$title,"_",simplify = T)[,2]pd$distance[97:128] <- rep("0cm",32)pd$patient <- str_split(pd$title,"_",simplify = T)[,3]pd$patient[97:128] <- str_split(pd$title[97:128],"_",simplify = T)[,2]str(pd) #查看数据结构save(expr,pd,file = "./Rdata/GSE90627_eSet.Rdata")# 文章的Supplementary Figure 1记录了每个患者的末端样本距离，复制到excel中，保存为sTable1.csvload("./Rdata/GSE90627_eSet.Rdata")# 读入取样末端距离数据NP <- read.csv("./raw-data/sTable1.csv",header = T, stringsAsFactors = F)# 将距离改为数字pd$distance2 <- substr(pd$distance,1,1)pd$distance2[1:32] <- NP$NP.cm. #前32个恰好是patient1-32的ProximalEnd样本pd$distance2 <- as.numeric(pd$distance2)pd$patient <- as.factor(pd$patient) 构建模型

我们先以基因“CXCL1”的表达量为例，基于线性混合效应模型计算基因表达量与距离的相关性。

首先，线性混合效应模型包括2个部分：

固定效应：被研究的变量，本例中是“距离”随机效应：希望被控制的影响因素，本例中是“患者”

常用的随机效应公式：

(1 | g) ：斜率固定，截距变化x + (x | g) ：斜率和截距都变化
其中，g 代表随机效应（因子型数据）；x 代表固定效应（数值型数据） df <- data.frame(CXCL1=expr["CXCL1",],pd[,c(5,7,8)])library(lme4)lmer.fit <- lmer(CXCL1~distance2+(distance2|patient),data = df)summary(lmer.fit)# summary 结果Linear mixed model fit by REML ['lmerMod']Formula: CXCL1 ~ distance2 + (1 | patient) Data: dfREML criterion at convergence: 456Scaled residuals: Min 1Q Median 3Q Max -1.91585 -0.89480 -0.05695 0.74833 2.34680 Random effects: Groups Name Variance Std.Dev. patient (Intercept) 0.01888 0.1374 Residual 1.94671 1.3952 Number of obs: 128, groups: patient, 32Fixed effects: Estimate Std. Error t value(Intercept) 11.55250 0.15984 72.275distance2 -0.20275 0.02039 -9.946Correlation of Fixed Effects: # 固定效应的相关性 (Intr)distance2 -0.618# 线性相关性coeff <- summary(lmer.fit)$coefficients[2,1] 模型可视化 library(ggeffects) library(ggplot2)pred.mm <- ggpredict(lmer.fit, terms = c("distance2")) ggplot(pred.mm) + geom_point(data = df,aes(x = distance2, y = CXCL1, colour = distance),position = "jitter") + geom_line(aes(x = x, y = predicted)) + # slope geom_ribbon(aes(x = x, ymin = predicted - std.error, ymax = predicted + std.error), fill = "lightgrey", alpha = 0.5) + # error band theme_minimal()

显著性检验

如何检验固定效应是否显著（即基因的表达量与距离的相关性是否显著）？

我们通过检验固定效应模型与去除固定效应模型的差异，判断固定相应的线性相关性。

注意：REML（residualmaximum likelihood），即残差最大似然，默认固定效应是正确的，lmer建模时默认REML=TRUE，所以在判断固定效应是否显著时，要将其设为FALSE，使用ML（最大似然）建模。

fit.full <- lmer(CXCL1~distance2+(distance2|patient),data = df,REML = F)fit.null <- lmer(CXCL1~(distance2|patient),data = df,REML = F)anova(fit.full,fit.null,test="LRT") #LRT代表极大似然检验# 结果Data: dfModels:fit.null: CXCL1 ~ (distance2 | patient)fit.full: CXCL1 ~ distance2 + (distance2 | patient) npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) fit.null 5 510.88 525.14 -250.44 500.88 fit.full 6 459.45 476.56 -223.73 447.45 53.428 1 2.683e-13 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1# pvaluepval <- anova(fit.full,fit.null,test="LRT")$`Pr(>Chisq)`[2]pval[1] 2.682786e-13