精华内容
下载资源
问答
  • 线性混合模型

    2021-02-27 20:10:47
    但由于重复测量数据、区组数据以及空间相关数据不能满足独立性假设,因此常常利用线性混合效应模型对上述数据进行分析(如相关性分析)。 以一个探索基因表达量与采样距离相关性为例。 该数据来自于“ntegrated ...

    背景

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

    使用一般线性模型时,是需要满足以下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 symbol
    rownames(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.csv
    load("./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: df
    
    REML criterion at convergence: 456
    
    Scaled 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, 32
    
    Fixed effects:
                Estimate Std. Error t value
    (Intercept) 11.55250    0.15984  72.275
    distance2   -0.20275    0.02039  -9.946
    
    Correlation 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: df
    Models:
    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
    
    # pvalue
    pval <- anova(fit.full,fit.null,test="LRT")$`Pr(>Chisq)`[2]
    pval
    [1] 2.682786e-13
    
    展开全文
  • 广义线性混合模型pdf

    2018-11-09 10:06:55
    广义线性混合模型高清完整版,英文原版。清晰度很高,主要介绍线性混合模型
  • 随着软件包的进步,使用广义线性混合模型(GLMM)和线性混合模型(LMM)变得越来越容易。由于我们发现自己在工作中越来越多地使用这些模型,我们开发了一套R shiny工具来简化和加速与对象交互的lme4常见任务。

     原文链接:http://tecdat.cn/?p=3138

     

    随着软件包的进步,使用广义线性混合模型(GLMM)和线性混合模型(LMM)变得越来越容易。由于我们发现自己在工作中越来越多地使用这些模型,我们开发了一套R shiny工具来简化和加速与对象交互的lme4常见任务。 

     

    shiny的应用程序和演示

    演示此应用程序功能的最简单方法是使用Shiny应用程序,在此处启动一些指标以帮助探索模型。 

     

    在第一个选项卡上,该函数显示用户选择的数据的预测区间。该函数通过从固定效应和随机效应项的模拟分布中抽样并组合这些模拟估计来快速计算预测区间,以产生每个观察的预测分布。

    在下一个选项卡上,固定效应和组级效果的分布在置信区间图上显示。这些对于诊断非常有用,并提供了检查各种参数的相对大小的方法。

    在第三个标签上有一些方便的方法,显示效果的影响或程度predictInterval。对于每种情况,最多12个,在所选数据类型中,用户可以查看更改固定效应的影响。这允许用户比较变量之间的效果大小,以及相同数据之间的模型之间的效果大小。

    预测

    预测像这样。

     

    predict(m1, newdata
    #>        1        2        3        4        5        6        7        8
    #> 3.146336 3.165211 3.398499 3.114248 3.320686 3.252670 4.180896 3.845218
    #>        9       10
    #> 3.779336 3.331012
    

    预测lmglm

    predInte(m1, newdata = Eval[1:10, ], n.sims = 500, level = 0.9,
    #>         fit      lwr      upr
    #> 1  3.074148 1.112255 4.903116
    #> 2  3.243587 1.271725 5.200187
    #> 3  3.529055 1.409372 5.304214
    #> 4  3.072788 1.079944 5.142912
    #> 5  3.395598 1.268169 5.327549
    #> 6  3.262092 1.333713 5.304931
    

    预测区间较慢,因为它是模拟计算。

     

    可视化

    可视化检查对象的功能。最简单的是得到固定和随机效应参数的后验分布。

     

    head(Sim)
    #>          term        mean      median         sd
    #> 1 (Intercept)  3.22673524  3.22793168 0.01798444
    #> 2    service1 -0.07331857 -0.07482390 0.01304097
    #> 3   lectage.L -0.18419526 -0.18451731 0.01726253
    #> 4   lectage.Q  0.02287717  0.02187172 0.01328641
    #> 5   lectage.C -0.02282755 -0.02117014 0.01324410
    

    我们可以这样绘制:

    pltsim(sim(m1, n.sims = 100), level = 0.9, stat = 'median'
    

     

    我们还可以快速制作随机效应的图:

    head(Sims)
    #>   groupFctr groupID        term        mean      median        sd
    #> 1         s       1 (Intercept)  0.15317316  0.11665654 0.3255914
    #> 2         s       2 (Intercept) -0.08744824 -0.03964493 0.2940082
    #> 3         s       3 (Intercept)  0.29063126  0.30065450 0.2882751
    #> 4         s       4 (Intercept)  0.26176515  0.26428522 0.2972536
    #> 5         s       5 (Intercept)  0.06069458  0.06518977 0.3105805
    
    plotR((m1, n.sims = 100), stat = 'median', sd = TRUE

    有时,随机效应可能难以解释

     Rank(m1, groupFctr = "d")
    head(ranks)
    #>      d (Intercept) (Intercept)_var       ER pctER
    #> 1 1866   1.2553613     0.012755634 1123.806   100
    #> 2 1258   1.1674852     0.034291228 1115.766    99
    #> 3  240   1.0933372     0.008761218 1115.090    99
    #> 4   79   1.0998653     0.023095979 1112.315    99
    #> 5  676   1.0169070     0.026562174 1101.553    98
    #> 6   66   0.9568607     0.008602823 1098.049    97

    效果模拟

     

    解释LMM和GLMM模型的结果很困难,尤其是不同参数对预测结果的相对影响。

    
    impact(m1, Eval[7, ], groupFctr = "d", breaks = 5,
    n.sims = 300, level = 0.9)
    
    #>   case bin   AvgFit     AvgFitSE nobs
    #> 1    1   1 2.787033 2.801368e-04  193
    #> 2    1   2 3.260565 5.389196e-05  240
    #> 3    1   3 3.561137 5.976653e-05  254
    #> 4    1   4 3.840941 6.266748e-05  265
    #> 5    1   5 4.235376 1.881360e-04  176

     

    结果表明yhat根据我们提供的newdata在组因子系数的大小方面,从第一个到第五个分位数的变化。

    ggplot(impSim, aes(x = factor(bin), y = AvgFit, ymin = AvgFit - 1.96*AvgFitSE,
    ymax = AvgFit + 1.96*AvgFitSE)) +

     

    非常感谢您阅读本文,有任何问题请在下面留言!


    最受欢迎的见解

    1.基于R语言的lmer混合线性回归模型

    2.R语言用Rshiny探索lme4广义线性混合模型(GLMM)和线性混合模型(LMM)

    3.R语言线性混合效应模型实战案例

    4.R语言线性混合效应模型实战案例2

    5.R语言线性混合效应模型实战案例

    6.线性混合效应模型Linear Mixed-Effects Models的部分折叠Gibbs采样

    7.R语言LME4混合效应模型研究教师的受欢迎程度

    8.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长

    9.使用SAS,Stata,HLM,R,SPSS和Mplus的分层线性模型HLM

    展开全文
  • 原文链接:R语言用Rshiny探索lme4广义线性混合模型(GLMM)和线性混合模型(LMM)​tecdat.cn随着lme4包装的进步,使用广义线性混合模型(GLMM)和线性混合模型(LMM)变得越来越容易。由于我们发现自己在工作中...

    原文链接:

    R语言用Rshiny探索lme4广义线性混合模型(GLMM)和线性混合模型(LMM)tecdat.cn
    6bf6181082885556f65ebd1acdd638dc.png

    随着lme4包装的进步,使用广义线性混合模型(GLMM)和线性混合模型(LMM)变得越来越容易。由于我们发现自己在工作中越来越多地使用这些模型,我们开发了一套R shiny工具来简化和加速与对象交互的lme4常见任务。

    shiny的应用程序和演示

    演示此应用程序功能的最简单方法是使用捆绑的Shiny应用程序,在此处启动一些指标以帮助探索模型。

    eccdd4b02238d0ea43e884b1c023eab7.png

    在第一个选项卡上,该函数显示用户选择的数据的预测间隔,这些预测间隔是使用predictInterval包中的函数计算的。该函数通过从固定效应和随机效应项的模拟分布中抽样并组合这些模拟估计来快速计算预测间隔,以产生每个观察的预测分布。

    52a6a2bcc55025659a6d457e022b6b9a.png

    在下一个选项卡上,固定效应和组级效果的分布在置信区间图上显示。这些对于诊断非常有用,并提供了检查各种参数的相对大小的方法。这个标签使得利用四个相关的功能merToolsFEsimplotFEsimREsimplotREsim它们可将自己的使用。

    11e5d7714aa2917bcd6690cef1d92acd.png

    在第三个标签上有一些方便的方法,通过利用的力量来显示效果的影响或程度predictInterval。对于每种情况,最多12个,在所选数据类型中,用户可以查看更改固定效果或其中一个分组级别术语的影响。使用该REimpact函数,如果所有其他情况保持相等,则使用模型的预测模拟每个案例,但是通过固定效应或随机效应项的分布来移动观察。这是在因变量的比例上绘制的,这允许用户比较变量之间的效果大小,以及相同数据之间的模型之间的效果大小。

    预测

    标准预测看起来像这样。

    <span style="color:#5c5c5c"><code>predict(m1, newdata = InstEval[1:10, ]) #> 1 2 3 4 5 6 7 8 #> 3.146336 3.165211 3.398499 3.114248 3.320686 3.252670 4.180896 3.845218 #> 9 10 #> 3.779336 3.331012</code></span>

    随着predictInterval我们获得更喜欢所生产的标准对象的预测lmglm

    <span style="color:#5c5c5c"><code>#predictInterval(m1, newdata = InstEval[1:10, ]) # all other parameters are optional predictInterval(m1, newdata = InstEval[1:10, ], n.sims = 500, level = 0.9, stat = 'median') #> fit lwr upr #> 1 3.074148 1.112255 4.903116 #> 2 3.243587 1.271725 5.200187 #> 3 3.529055 1.409372 5.304214 #> 4 3.072788 1.079944 5.142912 #> 5 3.395598 1.268169 5.327549 #> 6 3.262092 1.333713 5.304931 #> 7 4.215371 2.136654 6.078790 #> 8 3.816399 1.860071 5.769248 #> 9 3.811090 1.697161 5.775237 #> 10 3.337685 1.417322 5.341484</code></span>

    请注意,predictInterval它较慢,因为它是计算模拟。它还可以将所有模拟yhat值作为属性返回到预测对象本身。

    predictInterval大量使用包中的sim函数arm来绘制模型参数的分布。然后,它将这些模拟值组合起来,yhat为每个观测值创建分布。

    绘制

    merTools还提供了merMod可视化检查对象的功能。最简单的是得到固定和随机效应参数的后验分布。

    <span style="color:#5c5c5c"><code>feSims <- FEsim(m1, n.sims = 100) head(feSims) #> term mean median sd #> 1 (Intercept) 3.22673524 3.22793168 0.01798444 #> 2 service1 -0.07331857 -0.07482390 0.01304097 #> 3 lectage.L -0.18419526 -0.18451731 0.01726253 #> 4 lectage.Q 0.02287717 0.02187172 0.01328641 #> 5 lectage.C -0.02282755 -0.02117014 0.01324410 #> 6 lectage^4 -0.01940499 -0.02041036 0.01196718</code></span>

    我们也可以这样绘制:

    <span style="color:#5c5c5c"><code>plotFEsim(FEsim(m1, n.sims = 100), level = 0.9, stat = 'median', intercept = FALSE)</code></span>

    584b1814851b5b5fa9df78187386bb7f.png

    我们还可以快速制作随机效应的图:

    <span style="color:#5c5c5c"><code>reSims <- REsim(m1, n.sims = 100) head(reSims) #> groupFctr groupID term mean median sd #> 1 s 1 (Intercept) 0.15317316 0.11665654 0.3255914 #> 2 s 2 (Intercept) -0.08744824 -0.03964493 0.2940082 #> 3 s 3 (Intercept) 0.29063126 0.30065450 0.2882751 #> 4 s 4 (Intercept) 0.26176515 0.26428522 0.2972536 #> 5 s 5 (Intercept) 0.06069458 0.06518977 0.3105805 #> 6 s 6 (Intercept) 0.08055309 0.05872426 0.2182059</code></span>

    <span style="color:#5c5c5c"><code>plotREsim(REsim(m1, n.sims = 100), stat = 'median', sd = TRUE)</code></span>

    bfe495729b2052441e7daaa6a678f4b3.png

    有时,随机效应可能难以解释,并且并非所有这些都与零有意义地不同

    <span style="color:#5c5c5c"><code>ranks <- expectedRank(m1, groupFctr = "d") head(ranks) #> d (Intercept) (Intercept)_var ER pctER #> 1 1866 1.2553613 0.012755634 1123.806 100 #> 2 1258 1.1674852 0.034291228 1115.766 99 #> 3 240 1.0933372 0.008761218 1115.090 99 #> 4 79 1.0998653 0.023095979 1112.315 99 #> 5 676 1.0169070 0.026562174 1101.553 98 #> 6 66 0.9568607 0.008602823 1098.049 97</code></span>

    效果模拟

    解释LMM和GLMM模型的结果仍然很困难,尤其是不同参数对预测结果的相对影响。

    <span style="color:#5c5c5c"><code>impSim <- REimpact(m1, InstEval[7, ], groupFctr = "d", breaks = 5, n.sims = 300, level = 0.9) impSim #> case bin AvgFit AvgFitSE nobs #> 1 1 1 2.787033 2.801368e-04 193 #> 2 1 2 3.260565 5.389196e-05 240 #> 3 1 3 3.561137 5.976653e-05 254 #> 4 1 4 3.840941 6.266748e-05 265 #> 5 1 5 4.235376 1.881360e-04 176</code></span>

    结果REimpact表明yhat,根据我们提供的情况newdata,在组因子系数的大小方面,从第一个到第五个五分位数的变化。

    <span style="color:#5c5c5c"><code>library(ggplot2) ggplot(impSim, aes(x = factor(bin), y = AvgFit, ymin = AvgFit - 1.96*AvgFitSE, ymax = AvgFit + 1.96*AvgFitSE)) + geom_pointrange() + theme_bw() + labs(x = "Bin of `d` term", y = "Predicted Fit")</code></span>

    0bc9ef0ae5def4980bec03bec76719f5.png

    非常感谢您阅读本文,有任何问题请在下面留言!

    展开全文
  • 广义线性混合模型GLMM

    万次阅读 多人点赞 2019-03-09 22:46:44
    广义线性混合模型GLMM(Generalized Linear Mixed Model),是广义线性模型GLM 和线性混淆模型LMM 的扩展形式,于二十世纪九十年代被提出。GLMM因其借鉴了混合模型的思想,其在处理纵向数据(重复测量资料)时,被...

    广义线性混合模型GLMMGeneralized Linear Mixed Model),是广义线性模型GLM线性混淆模型LMM 的扩展形式,于二十世纪九十年代被提出。GLMM因其借鉴了混合模型的思想,其在处理纵向数据重复测量资料)时,被认为具有独特的优势。GLMM不仅擅长处理重复测量资料,还可以用于任何层次结构的数据(因为本质上又是多水平模型)。

    提到GLMM,有必要先介绍几个容易混淆的概念:GLM、LMM、MLM、GMM 和GEE。

     

    相关模型简介

    广义线性模型 GLM

    广义线性模型GLM,是大家经常接触的概念了,比如经典的Logistic模型。GLM是普通线性模型的扩展形式,由于普通线性回归的因变量必须服从正态分布,而实际问题中经常会遇到分类问题或计数问题的建模,GLM采用连接函数(Link Function),将因变量的分布进行了扩展,使得因变量只要服从指数分布族即可(如正态分布,二项分布,泊松分布,多项分布等)。

    GLM 可以分解为 Random Component、System Component 和 Link Function 三个部分。Random Component 为残差部分,取决于因变量的分布;System Component 为预测部分,又称 linear predictor,是拟合的关键;Link Function 为连接变化函数,用于将指数分布族转化成正态分布,或者说,对预测结果进行非线性映射(建立 linear predictor与 label 之间的变换关系),是LM成长为GLM的关键环节。

    需要强调的是,link function 是从 label 映射到 linear predictor的过程,link function的反函数称为响应函数 response function。响应函数 把 linear predictor 直接映射到了预测目标 label。较常见的link function如 logit函数(又称log-odds);较常用的响应函数如 logistic(又称sigmoid,是二分类中的相应函数)和 softmax(是sigmoid的扩展形式,用于多分类问题),这两个都是 logit 的反函数。

    Logistic为例,如下(本部分摘自:GLM(广义线性模型) 与 LR(逻辑回归) 详解):

    最后啰嗦一句,因变量为Bernoulli Distribution也就是对二分类问题建模,因变量为Binomial Distribution也就是对多分类问题建模,因变量为Poisson Distribution也就是对计数问题建模(注意区分计数问题和多分类问题)。

    本文讲得比较简略,有两篇博客对GLM总结得比较棒,给出链接如下,值得一读:

    GLM(广义线性模型) 与 LR(逻辑回归) 详解

    广义线性模型GLM

     

    线性混合模型 LMM

    本部分参考自:《高级医学统计学》和 Wiki: Mixed_model

    线性混合模型LMM,又称混合线性模型MLM、混合模型MM、多水平模型MLM、随机系数模型RCM、等级线性模型HLM 等。首先看一下 Wiki上对混合模型MM的介绍:A mixed model (or more precisely mixed error-component model) is a statistical model containing both fixed effects and random effects. (注意:fixed在这里译为固定,不同于mixed混合)

    混合模型擅长于处理纵向数据(重复测量数据)和有缺失的数据,并且往往优于ANOVA等方法。

    在混合模型中,需要区分两个概念:random effectsrandom errors

    以矩阵定义混合模型,可以写成:

          y=X\beta+Z\gamma +\epsilon

    y 是观测值的向量,服从多元正态分布,且平均值可以表示为 E(y)=X\beta

    \beta 是固定因子的效应值(与X对应的固定效应参数向量)

    \gamma 是随机因子的效应值,服从多元正态分布,且平均值为 E(\gamma)=0,它的方差为 Var(\gamma)=G

    \epsilon 是残差的向量矩阵,它的平均值为 E(\epsilon )=0,它的方差为 Var(\epsilon )=R

    X 为固定效应自变量的设计矩阵(可包括连续性变量和分类变量,甚至可包含交互项或二次项等),Z 为随机效应变量构造的设计矩阵。

    [ 注意:切勿将固定效应狭义理解为主要变量,而应该是所有可能的解释变量(如分组变量和时间变量),包括这些变量之间的交互项。而随机效应则是假定的随机效应部分(这部分的意义应当从多水平模型的角度来理解了) ]

    该模型为固定效应 X\beta 和随机效应 Z\gamma 的混合,且固定效应和随机效应均与响应变量为线性关系,因此称为线性混合模型。

    注意:当满足球形检验时,重复测量资料的线性混合效应模型可退化为一般线性模型。

     

    混合模型的假定为 \gamma\sim N(0,G)\epsilon \sim N(0,R),其中 Cov(\gamma,\epsilon )=0,即两者的协方差为0(二者互相独立)。可以给出Henderson's "mixed model equations" (MME):

          \begin{pmatrix} X'R^{-1}X & X'R^{-1}Z \\ Z'R^{-1}X & Z'R^{-1}Z + G^{-1} \end{pmatrix} \begin{pmatrix} \hat{\boldsymbol{\beta}} \\ \hat{\boldsymbol{u}} \end{pmatrix} = \begin{pmatrix} X'R^{-1}\boldsymbol{y} \\ Z'R^{-1}\boldsymbol{y} \end{pmatrix}

    The solutions to the MME, \hat{\beta} and \hat{u} are best linear unbiased estimates (BLUE) and predictors (BLUP) for \beta and u(此处的u指的就是\gamma,有的版本习惯使用 u 来替代 \gamma 字符), respectively. 拟合混合模型还可以使用 EM 算法。

    工具包:R (nlme包中的lme方法,或 lme4 包中的lmer方法), Python (statsmodels包)。

     

    多水平模型 MLM

    本部分参考自:《高级医学统计学》

    多水平模型其实和线性混合模型LMM是等价的,只是理解的角度不同而已。MLM是从模型组建的多个水平来理解,关注构建过程;LMM则仅关注模型构建的结果(固定效应部分+随机效应部分)。多水平模型可以分层表述,然后整合成一个公式(即等价于LMM的公式)。下面以两水平模型为例,进行解读。

    一个包含“2个水平1的解释变量(x和z)和1个水平2的解释变量(w)”的两水平模型可以表述为:

          \\ y_{ij}=\beta _{0j}+\alpha _1x_{1ij}+\beta _{1j}z_{1ij}+e_{ij} \\ \beta _{0j}=\gamma _{00}+\gamma _{01}w_{1j}+u_{0j} \\ \beta _{1j}=\gamma _{10}+\gamma _{11}w_{1j}+u_{1j}

    其中,i=1,2,...,N (N是总样本量),j=1,2,...,J (J是水平2的解释变量的w的取值个数,假定w为分类变量)。则 y_{ij} 表示在变量w的第 j 种取值的情况中的第 i  个个体的结局测量值。第1水平方程(第1个等式)中,截距\beta _{0j} 带有下标 j,表示其值随 w 的取值变化而变化;系数\beta _{1j} 带有下标 j,表示变量z_{1ij} 对 y_{ij} 的效应随 w 的取值变化而变化;而系数\alpha _1 不带有下标 j,表示变量x_{1ij} 对 y_{ij} 的效应不随 w 的取值变化而变化。在两个第2水平方程(第2、3个等式)中,第1水平的回归系数变成了因变量。关于其他参数如e和u的规则,此处跳过(感兴趣的可查阅统计书《高级医学统计学》)。

    从概念上来讲,该模型的建立是从顶向下的,先进行第1水平的参数计算(通过枚举 j 来获得 j 组回归系数\beta _{0j}\beta _{1j});然后使用估计的回归系数进行第2水平的参数计算,生成多个第2水平的方程。这种计算步骤是传统的计算方法,现在的计算其实是同步进行的。

    如果将两个第2水平的方程代入到第1水平的方程中,可以得到:

          y_{ij}=(\gamma _{00}+\gamma _{01}w_{1j}+\alpha _1x_{1j}+\gamma _{10}z_{1ij}+\gamma _{11}w_{1j}z_{1ij})+(u_{0j}+u_{1j}z_{1ij}+e_{ij})

    这是一个组合模型,该式右边分为两部分,第一个括号部分是各个解释变量及其交互项产生的效应,第二个括号部分是复合残差结构。第一部分便可对应为LMM中提到的固定效应部分,第二部分可对应为LMM中提到的随机效应部分(包括纯粹残差项)。

     

    更一般地,两水平模型可表述为:

          \\ y_{ij}=\beta _{0j}+ \sum_{p=1}^P\alpha _px_{pij}+\sum_{q=1}^Q\beta _{qj}x_{qij}+e_{ij} \\ \beta _{0j}=\gamma _{00}+\sum_{m=1}^M\gamma _{0m}w_{mj}+u_{0j} \\ ... \\ \beta _{Qj}=\gamma _{Q0}+\sum_{m=1}^M\gamma _{Qm}w_{mj}+u_{Qj}

    将Q个第2水平的方程代入到第1水平的方程中,可以得到:

          \\ y_{ij}=(\gamma _{00}+\sum_{m=1}^M\gamma _{0m}w_{mj}+\sum_{p=1}^P\alpha _px_{pij}+\sum_{q=1}^Q\gamma _qz_{qij}+\sum_{q=1}^Q\sum_{m=1}^M\gamma _{qm}w_{mj}z_{qij}) \\+(u_{0j}+\sum_{q=1}^Qz_{qij}u_{qj}+e_{ij})

    该组合模型由两部分组成:固定效应部分(第一个括号中)和随机效应部分(第二个括号中)。

    MLM的参数估计十分复杂,模型构建的步骤也比较繁琐,此处都不进行讲解。

     

    高斯混合模型 GMM

    高斯混合模型GMM(Gaussian Mixed Model)指的是多个高斯分布函数的线性组合,理论上GMM可以拟合出任意类型的分布,通常用于解决同一集合下的数据包含多个不同的分布的情况(或者是同一类分布但参数不一样,或者是不同类型的分布,比如正态分布和伯努利分布)。

    设有随机变量X,则混合高斯模型可以用下式表示:

          

     为混合模型中的第k 个分量(component)。比如有两个聚类,可以用两个二维高斯分布来表示,那么分量数K=2 \pi_k 是混合系数(mixture coefficient),且满足:

          \sum _{k=1}^K\pi_k=1, 0\leq \pi_k\leq 1

    实际上,可以认为\pi_k就是每个分量的权重。

    GMM常用于聚类。如果要从 GMM 的分布中随机地取一个点的话,实际上可以分为两步:首先随机地在这 K 个 Component 之中选一个,每个 Component 被选中的概率实际上就是它的系数\pi_k,选中 Component 之后,再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布,转化为已知的问题。

    可以用EM算法估计GMM参数。

    此处介绍较为简略,有一篇博客讲得比较好,值得参考:高斯混合模型(GMM)及其EM算法的理解

    另外,还有个概念叫广义矩方法,也简称GMM,其与GEE密切相关,可参考:广义矩方法(GMM)和广义估计方程(GEE)到底有什么区别

     

    广义估计方程 GEE

    广义估计方程(generalized estimating equation, GEE)用于估计广义线性模型的参数(其中线性模型的结果之间可能存在未知的相关性)。于1986年由Liang和Zeger首次提出,是在广义线性模型和重复测量数据中,运用准似然估计方法估计参数的一种用于分析相关性数据的回归模型。

    详细介绍请参考这篇博客:广义估计方程GEE

     

    广义线性混合模型 GLMM

    广义线性混合模型GLMM,可以看做是线性混合模型LMM的扩展形式,使得因变量不再要求满足正态分布;也可以看作是GLM的扩展形式,使得可以同时包含固定效应随机效应

    回顾一下,LMM模型的一般形式为:

          y=X\beta+Z\gamma +\epsilon

    y 是N*1的向量,表示观测值;X是N*p的矩阵,表示固定效应自变量;\beta 是p*1的向量,表示固定效应参数向量;Z是N*q的矩阵,表示随机效应变量;\gamma 是q*1的向量(\gamma 在某些版本中也写成u),表示随机因子的效应值;\epsilon 是N*1的向量,表示残差(随机误差)。

          

    GLMM在此基础上做了一些改动。令 linear predictor, \eta, 表示固定效应和随机效应的组合(随机误差不包含在内),即:

          \eta =X\beta+Z\gamma

    令g(⋅)表示link function,用来连接 linear predictor 和 label,h(⋅)为g(⋅)的反函数,即response function。则有:

          g(E(y))=\eta , E(y)=H(\eta )=u, 因此:y=h(\eta )+\epsilon

    此处的 link function 和 response function 的示例,请直接参考GLM中的介绍(但此处会额外接触到几个概念:带随机效应的Logistic回归中的 probability density function 或简称PDF,和带随机效应的Poisson回归中的probability mass function 或简称PMF)。结果的解读,和GLM中的解读类似,细微的差别仅在于随机效应部分的解读。

    借鉴知乎上的一个理解:

    举个例子,我们认为疗效可能与服药时间相关,但是这个相关并不是简简单单的疗效随着服药时间的变化而改变。更可能的是疗效的随机波动的程度与服药时间有关。比如说,在早上10:00的时候,所有人基本上都处于半饱状态,此时吃药,相同剂量药物效果都差不多。但在中午的时候,有的人还没吃饭, 有的人吃过饭了,有的人喝了酒,结果酒精和药物起了反应,有的人喝了醋,醋又和药物起了另一种反应。显然,中午吃药会导致药物疗效的随机误差非常大。这种疗效的随机误差(而非疗效本身)随着时间的变化而变化,并呈一定分布的情况,必须用广义线性混合模型了。对于固定效应来说,参数的含义是,自变量每变化一个单位,应变量平均变化多少。而对于随机效应而言,参数是服从正态分布的一个随机变量,也就是说对于两个不同的自变量的值,对应变量的影响不一定是相同的。

     

    一篇文献以一个案例对以上几种模型进行了比较,值得一读:GEE、GLMM和MLM分析卫生重复测量资料的效果比较

     

    参考资料

    万崇华等. 高级医学统计学. 科学出版社.

    Wiki: Generalized_linear_model

    Wiki: Mixed_model

    Wiki: Generalized_linear_mixed_model

    Introduction to generalized linear mixed models

    GLM(广义线性模型) 与 LR(逻辑回归) 详解

    广义线性模型GLM

    广义估计方程GEE

    高斯混合模型(GMM)及其EM算法的理解

    混合模型初探

    周婷,兰蓝,邱建青,杜春霖,李晓松,张韬.GEE、GLMM和MLM分析卫生重复测量资料的效果比较[J].现代预防医学,2017,44(16):2881-2885+2899

    展开全文
  • 广义线性混合模型

    千次阅读 2019-02-16 18:13:43
    1.线性模型和线性混合模型区别 线性模型的表达式为:pitch~age+ε. 即两部分:固定项age和误差项ε。 广义线性混合模型表达式为:pitch~age+(1|subject)+ε 三部分:固定项age,随机项(1|subject)和误差项ε。 为...
  • 43、广义线性混合模型-2020.09.14 43、广义线性混合模型-2020.09.14 43、广义线性混合模型-2020.09.14 43、广义线性混合模型-2020.09.14 43、广义线性混合模型-2020.09.14
  • 线性混合模型中方差分量的经验Bayes检验
  • 而GLMM (generalized linear mixed model)是广义线性混合模型。 广义线性模型GLM很简单,举个例子,药物的疗效和服用药物的剂量有关。这个相关性可能是多种多样的,可能是简单线性关系(发烧时吃一片...
  • 一般线性模型混合线性模型 生命科学的数学统计和机器学习 (Mathematical Statistics and Machine Learning for Life Sciences) This is the seventeenth article from my column Mathematical Statistics and ...
  • Pavilion非线性混合模型在塞克LLDPE装置上的应用pdf,Pavilion非线性混合模型在塞克LLDPE装置上的应用
  • 一般线性模型混合线性模型 生命科学的数学统计和机器学习 (Mathematical Statistics and Machine Learning for Life Sciences) This is the eighteenth article from the column Mathematical Statistics and ...
  • 线性和广义线性混合模型及其统计诊断 [费宇 编著] 2013年版
  • 欢迎关注”生信修炼手册”!Linear Mixde Model, 简称LMM, 称之为线性混合模型。从名字也可以看出,这个模型和一般线性模型有着很深的渊源。线性混合模型是在一般线性模型的...
  • 该R包设计为使用通用线性混合模型(glmm)对基因表达进行建模。这使我们可以包括随机效果和固定效果。对于包的目的,我们使用glmer功能从这符合一个GLMM包。 该软件包特别关注不同React组或治疗组之间基因表达随时间...
  • 本文阐述了如何用极大似然估计(ML)和约束极大似然估计(REML)对线性混合模型(多层线性模型)进行参数估计,给出了对数似然函数,梯度和黑塞矩阵的计算过程,以及优化矩阵求逆的Woodbury恒等式,并且用python实现...
  • 进行数据分析时,会发现有时候一个模型中的变量之间可能具有相关性(correlation),比如面积和长度就具有高度的相关性,如果同时对这些参数建模,就存在共线性问题,所以一般是只针对其中一个参数建模。而这种...
  • 基于线性混合模型/神经网络模型的高光谱混合像元分解算法比较-以徐州地区Hyperion影像为例,蔡燕,赵银娣,混合像元的普遍存在使得传统的基于像元级的遥感分类和面积测量很难达到使用要求的精度,所以混合像元的...
  • 该文件介绍了一种新的线性混合模型参数的估计方法即谱分解估计
  • 线性混合模型与普通的线性模型不同的地方是除了有固定效应外还有随机效应。 ___________________________________________________________________________________ 一、线性混合模型理论 由两个部分来决定...
  • 1.线性混合模型(linear mixed model): h(x,y)代表传感器在像素(x,y)坐标处采集的特征,该特征可以认为是一个N维(N-d)的向量,N为光谱波段的范围;也可以认为是端元向量ei(纯像素或像素只代表图像中的一个...
  • scikit-learn包装器,用于R中的广义线性混合模型方法 这是一个轻量级的包装器,可以通过R通过python拟合广义线性多元多级模型。它可以通过brms调用轻松使贝叶斯模型拟合: : 它具有足够的灵活性以扩展到其他基于R...
  • 包含一组R函数和脚本,用于基于引导方案或Wald-z来计算百分数置信区间,用于使用经典方法和鲁棒方法估算的线性混合模型。 该脚本允许平衡的纵向数据设计。 实现了两种引导程序方案:野生引导程序和参数引导程序。 ...
  • 原文链接:如何解决线性混合模型中畸形拟合(Singular fit)的问题​tecdat.cn假设我们有一个模型mod <- Y ~ X*Condition + (X*Condition|subject)# Y = logit variable # X = continuous variable # Condition = ...
  • 线性模型是我们最常见到的、最理想的数学模型,基本的线性模型是数据科学...然而现实生活中的线性问题,很大几率不适用于基本的线性模型,需要使用线性混合模型来描述。Tensorflow edward提供对这类问题的解决方案。
  • iMap4:iMap4-使用线性混合模型的眼动数据(例如,注视图)的空间映射

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,554
精华内容 621
关键字:

线性混合模型