R语言实战多水平模型！（线性混合模型详细解读，初学者必备！）

医学和生信笔记

发布于 2026-03-17 17:35:36

1750

文章被收录于专栏：医学和生信笔记医学和生信笔记

关注公众号，发送R语言或python，可获取资料

💡专注R语言在🩺生物医学中的使用

设为“星标”，精彩不错过

本文目录：

理论知识
数据探索
空模型
添加1级水平的固定效应
添加2级水平的固定效应
具有随机斜率的MLM
具有交互效应的MLM
重复测量数据的MLM
广义混合效应模型
参考资料

unsetunset理论知识unsetunset

理论知识我就不献丑了，直接给大家推荐高手的解读！

关于多水平模型（multi-level models，MLM）的概念和理论知识，强烈推荐阅读冯国双老师的几篇文章，这是我目前见过的写的最通俗易懂的。多水平模型、混合模型、随机效应模型、固定效应模型、随机系数模型、方差成分模型等等，全都详细介绍到了：

我之前接触多水平模型很少，对一些概念很模糊，但是读完这些推文后，我感到豁然开朗、恍然大悟，强烈推荐大家仔细读一读，对于初学者帮助很大。

上面是理论知识，下面是R语言实战。

unsetunset数据探索unsetunset

数据可以从这个网站免费下载：https://www.learn-mlms.com/13-appendix.html。或者加入QQ群免费下载。

library(dplyr) # 数据操作
library(ggplot2) # 可视化
library(lme4) # 多水平模型
library(lmerTest) # 计算P值
library(performance) # 计算模型表现
data <- read.csv('datasets/heck2011.csv') # 加载数据

简单介绍下这个数据，一共有6871行，各个变量的含义如下，其中math是因变量，建模目的是用其他变量预测学生的数学成绩（math），或者叫探索和学生成绩有关系的变量：

schcode:学校id，一共有419所学校
Rid:每个学校内的学生id，每个学校内都是从1开始
id:学生id，从1到6871
female:是否是女性，1=是，0=否
ses:衡量学校内学生的社会经济地位构成的Z-score
femses:以性别（女性）为中心的测量学生社会经济地位的变量
math:学生数学成绩测试分数
ses_mean:以学校为单位衡量学生社会经济地位的均值
pro4yrc:每个学校中计划就读四年制大学的学生比例
public:学校类型，1=公立学校，0=私立学校

str(data)
## 'data.frame':    6871 obs. of  10 variables:
##  $ schcode : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Rid     : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ id      : int  6701 6702 6703 6704 6705 6706 6707 6708 6709 6710 ...
##  $ female  : int  1 1 1 0 0 0 0 1 0 1 ...
##  $ ses     : num  0.586 0.304 -0.544 -0.848 0.001 -0.106 -0.33 -0.891 0.207 -0.341 ...
##  $ femses  : num  0.586 0.304 -0.544 0 0 0 0 -0.891 0 -0.341 ...
##  $ math    : num  47.1 63.6 57.7 53.9 58 ...
##  $ ses_mean: num  -0.267 -0.267 -0.267 -0.267 -0.267 ...
##  $ pro4yrc : num  0.0833 0.0833 0.0833 0.0833 0.0833 ...
##  $ public  : int  0 0 0 0 0 0 0 0 0 0 ...

这个数据很明显是有层次的，学生是分别属于不同的学校的，所以这是一个两个水平的数据。学生是1级水平，学校是2级水平。

不同学校之间的水平是有差异的，为了演示先选择其中10个学校（一共有419个学校）：

data_sub <- data %>% 
  filter(schcode <= 10)

如果不考虑学校水平的差异，探索社会经济地位（ses）和成绩（math）之间的关系，可以画出如下的散点图和拟合线：

data_sub %>% 
  ggplot(mapping = aes(x = ses, y = math)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, fullrange = TRUE)
## `geom_smooth()` using formula = 'y ~ x'

如果考虑到学校之间的不同水平，按照学校分别拟合，会得到如下的图：

data_sub %>% 
  ggplot(mapping = aes(x = ses, y = math, colour = factor(schcode))) +
  geom_point() +
  geom_smooth(mapping = aes(group = schcode), method = "lm", se = FALSE, fullrange = TRUE) +
  labs(colour = "schcode")
## `geom_smooth()` using formula = 'y ~ x'

每个学校的ses的截距和斜率都是不一样的，所以不同学校对学生成绩的影响也是不同的，普通的多元线性回归没有考虑到这种差异，但是多水平模型可以发现这些差异并量化它们。

unsetunset空模型unsetunset

先给大家演示一下最简单的，也就是“只有随机截距的模型”（random intercept only model），又被称为空模型（null model）。在这种模型中，没有任何自变量来解释学生成绩的差异，模型只考虑了学校间的差异，这些差异是由随机截距来表示的。

null_model <- lmer(math ~ (1|schcode), data = data)
summary(null_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ (1 | schcode)
##    Data: data
## 
## REML criterion at convergence: 48877.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6336 -0.5732  0.1921  0.6115  5.2989 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schcode  (Intercept) 10.64    3.262   
##  Residual             66.55    8.158   
## Number of obs: 6871, groups:  schcode, 419
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  57.6742     0.1883 416.0655   306.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

在lme4和lmerTest中（lmerTest可以在结果中给出P值，lme4不会给出P值，这是这两个包的差别），拟合多水平模型的基本语法是：

# DV是因变量，IV是自变量
lmer(DV ~ 1 + IV1 + IV2 + ... + IVp + (random_effect1 + random_effect2 + ... + random_effect3 | grouping_variable), data = dataset)

该数据的因变量是math，波浪号右边的(1|schcode)表示分组变量schcode的随机截距，空模型没有其他自变量。

结果主要是两部分：

Random effects：随机效应的估计
Fixed effects：固定效应的估计

先看固定效应：截距的固定效应为57.67，表示所有学校的平均数学成绩，此时这个截距其实就是因变量的平均值，即：mean(data$math)=57.73391（微小差异，可以忽略）；

再看随机效应：不同学校之间成绩的方差（Variance）为10.64，不同学校之间的学生成绩会围绕总体均值（也就是此时的截距57.67）波动，其标准差（Std.Dev）是3.262，同一个学校内不同学生之间成绩的方差为66.55，标准差为8.158。

组内相关系数（intraclass correlation coefficient，ICC）是衡量多水平模型中群体间变异（即随机效应）与总变异（群体间变异+群体内变异）之比的一种指标。在多水平模型中，ICC用来评估因变量的变异有多少可以归因于群体间的差异，而有多少是群体内的个体差异。

若ICC为0，说明群体间的差异为0，说明数据不存在层次结构，在本例中也就是学校之间没有差异，可以不使用多水平模型，使用普通的多元线性回归即可。

根据模型输出，该例的总方差为10.64+66.55=77.19（var(data$math)），ICC=10.64/77.19=0.138；也就是说成绩总变异中的13.8%是学校不同导致的，剩下的86.2%的变异是由同一学校内学生之间的差异造成的。也就是说，大部分成绩变异来源于学校内部的差异，而学校之间的差异相对较小。

ICC也可以使用performance包自动计算：

performance::icc(null_model)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.138
##   Unadjusted ICC: 0.138

由于没有其他自变量，所以调整的ICC和未调整的ICC是一样的。

unsetunset添加1级水平的固定效应unsetunset

上面演示的空模型没有自变量，下面我们添加一个社会经济状况（ses）作为自变量，这个ses是在水平1单位（不同的学生，不同的学校属于水平2单位）间变化的，所以这个变量产生的效应属于水平1单位的固定效应。

像这种只有随机截距，没有随机斜率（但是有自变量，不是“空模型”）的多水平模型被称为随机截距模型，又叫方差成分模型，相关概念的理解请参考本文开头推荐的几篇文章。

用lme4拟合随机截距模型：

# RMEL表示限制性最大似然估计，这也是默认方法
ses_l1 <- lmer(math ~ ses + (1|schcode), data = data, REML = TRUE)
summary(ses_l1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ ses + (1 | schcode)
##    Data: data
## 
## REML criterion at convergence: 48215.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7733 -0.5540  0.1303  0.6469  5.6908 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schcode  (Intercept)  3.469   1.863   
##  Residual             62.807   7.925   
## Number of obs: 6871, groups:  schcode, 419
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   57.5960     0.1329  375.6989  433.36   <2e-16 ***
## ses            3.8739     0.1366 3914.6382   28.35   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## ses -0.025

先看固定效应：Intercept的估计值为57.5960，表示当所有自变量都为0的时候，因变量的预测值；ses的估计值为3.8739，表示ses每增加一个单位，数学成绩可以提高3.8739分；固定效应的解读和普通的多元线性回归没有差别。

加入自变量后，我们可以发现固定效应中截距的值发生了变化，因为此时的截距不再是因变量的平均值，它是一个基于当前模型的预测值。

再看随机效应：schcode的方差为3.469，残差的方差为62.807，该结果与上面的模型的结果解读是类似的。schcode作为分组变量，它的方差越大说明这个变量对因变量的影响越大，残差的方差表示模型无法解释的随机误差部分，这个值越大说明模型无法解释的随机误差越大。

未调整的ICC可以量化分层变量（这里是schcode）所能解释的方差（变异）比例。在模型中加入ses这个自变量后，未调整的ICC=0.046：

performance::icc(ses_l1)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.052
##   Unadjusted ICC: 0.046

说明在考虑社会经济地位的影响后，数学成绩差异中有4.6%是学校不同造成的。可以看到加入ses这个变量后，不同学校能够解释的变异减少了（13.8%到4.6%），这个也很好理解，因为有一部分变异被ses解释了，另外残差的变异也减小了（66.55到62.807），也是这个原因导致的。

多水平模型肯定是要比单水平模型的拟合程度更好的，因为它能够解释分组变量导致的变异，也就是让不能解释的变异更少了。下面我们拟合一个普通的多元线性回归，并比较一下两个模型。

f <- lm(math ~ ses, data = data)
compare_performance(f,ses_l1,metrics = "common") # 比较下两个模型
## # Comparison of Model Performance Indices
## 
## Name   |           Model |   AIC (weights) |   BIC (weights) |  RMSE |    R2 | R2 (adj.) | R2 (cond.) | R2 (marg.) |   ICC
## --------------------------------------------------------------------------------------------------------------------------
## f      |              lm | 48304.0 (<.001) | 48324.5 (<.001) | 8.131 | 0.143 |     0.143 |            |            |      
## ses_l1 | lmerModLmerTest | 48219.1 (>.999) | 48246.4 (>.999) | 7.810 |       |           |      0.167 |      0.121 | 0.052

结果表明，AIC、BIC、RMSE都是多水平模型更小，也就是模型表现更好。

多元线性回归中常用的提取结果的方法也都适用于多水平模型，比如计算可信区间等：

confint(ses_l1)
## Computing profile confidence intervals ...
##                 2.5 %    97.5 %
## .sig01       1.575429  2.144559
## .sigma       7.789560  8.063828
## (Intercept) 57.335234 57.856673
## ses          3.596455  4.152745

unsetunset添加2级水平的固定效应unsetunset

前面我们添加了ses作为1级水平的自变量，以解释学生在数学成绩方面的部分差异。下面我们再添加一个在水平2单位（不同的学校）间变化的自变量，该自变量在不同的学校间是不同的，但是在同一所学校内的所有学生中是相同的，符合条件的自变量有3个：

ses_mean:以学校为单位衡量学生社会经济地位的均值
pro4yrc:每个学校中计划就读四年制大学的学生比例
public:学校类型，1=公立学校，0=私立学校

我们选择public作为水平2单位的固定效应，这个模型依然是一个随机截距模型，或者叫方差成分模型：

ses_l1_public_l2 <- lmer(math ~ ses + public + (1|schcode), 
                         data = data, REML = TRUE)
summary(ses_l1_public_l2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ ses + public + (1 | schcode)
##    Data: data
## 
## REML criterion at convergence: 48216
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7718 -0.5541  0.1309  0.6477  5.6916 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schcode  (Intercept)  3.486   1.867   
##  Residual             62.807   7.925   
## Number of obs: 6871, groups:  schcode, 419
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   57.63143    0.25535  381.81733 225.693   <2e-16 ***
## ses            3.87338    0.13673 3928.37427  28.329   <2e-16 ***
## public        -0.04859    0.29862  385.93649  -0.163    0.871    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) ses   
## ses     0.013       
## public -0.854 -0.031

先看固定效应：Intercept的估计值为57.5960，表示当所有自变量都为0的时候，因变量的预测值；ses的估计值为3.87338，表示ses每增加一个单位，数学成绩可以提高3.87338分；public的估计值为-0.04859，说明相对于私立学校，公立学校在数学成绩上平均降低0.04859分。固定效应的解读和普通的多元线性回归没有差别。

再看随机效应：结果解读和上面一个模型的解读类似的，就不再重复了。schcode的方差变大了一些，残差的方差没有变化。

理论上如果新加入的自变量能够解释更多的因变量变异，那么残差的变异（方差）通常会减少，群体间（在本例中也就是学校间）的变异也会减少，因为这部分变异都被新加入的自变量解释了。但是很明显我们新加入的这个public自变量不太行，它几乎解释不了因变量的变异，从它的系数也可以看出来，只有-0.04859，比ses差远了。说明公立学校和私立学校对数学成绩影响很小（P值也表明这个变量没有统计学意义）。

我们也可以通过计算模型的一些指标看看这个自变量到底行不行，并且和前面的模型比较一下：

compare_performance(null_model,ses_l1,ses_l1_public_l2)
## # Comparison of Model Performance Indices
## 
## Name             |           Model |   AIC (weights) |  AICc (weights) |   BIC (weights) | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
## ------------------------------------------------------------------------------------------------------------------------------------------
## null_model       | lmerModLmerTest | 48881.8 (<.001) | 48881.8 (<.001) | 48902.3 (<.001) |      0.138 |      0.000 | 0.138 | 7.977 | 8.158
## ses_l1           | lmerModLmerTest | 48219.1 (0.729) | 48219.1 (0.729) | 48246.4 (0.988) |      0.167 |      0.121 | 0.052 | 7.810 | 7.925
## ses_l1_public_l2 | lmerModLmerTest | 48221.1 (0.271) | 48221.1 (0.271) | 48255.2 (0.012) |      0.167 |      0.121 | 0.053 | 7.810 | 7.925

结果发现，加入public后，AIC、BIC还变大了，R2没啥太大的变化，充分说明这个变量真的作用不大，可以不加。

unsetunset具有随机斜率的MLMunsetunset

前面我们选择了10个学校，以展示了不同学校间数学成绩（math）与社会经济状况（ses）之间的关系：

从图中可以看出，不同学校ses的截距和斜率值差异很大。例如，学校3的截距约为38，斜率较小且为正值，而学校8的截距约为55，斜率较大且为正值。

在ses_l1这个模型中，我们假定每个学校ses的斜率是相同的，只估计了随机截距的变异，但是从图中可以看出，其实每个学校ses的斜率都是不一样的。也就是说，我们只估计了ses的平均效应，却忽略了不同学校的ses是不同的。

下面我们在模型中添加一个随机斜率项来模拟学校间ses斜率的这种变异。

在lme4的语法中，只需要将想要估计的具有不同斜率的变量名字放在|前面即可：

# 估计ses的随机斜率和随机截距
ses_l1_random <- lmer(math ~ ses + (ses|schcode), 
                      data = data, REML = TRUE)
## boundary (singular) fit: see help('isSingular')
summary(ses_l1_random)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ ses + (ses | schcode)
##    Data: data
## 
## REML criterion at convergence: 48190.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8578 -0.5553  0.1290  0.6437  5.7098 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  schcode  (Intercept)  3.2042  1.7900        
##           ses          0.7794  0.8828   -1.00
##  Residual             62.5855  7.9111        
## Number of obs: 6871, groups:  schcode, 419
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   57.6959     0.1315  378.6378  438.78   <2e-16 ***
## ses            3.9602     0.1408 1450.7730   28.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## ses -0.284
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

上面的公式也可以写成math~ses+(1+ses|schcode)，1表示随机截距项，这是默认设置，也可以仅使用(ses|schcode)估计随机截距和随机斜率。如果你想从模型中排除随机截距，需要写成(0+ses|schcode)来覆盖默认设置。像这种既有随机截距又有随机斜率的模型又被称为随机系数模型。

先看固定效应：Intercept的估计值为57.6959，表示当所有自变量都为0的时候，因变量的预测值；ses的估计值为3.9602，表示ses每增加一个单位，数学成绩可以提高3.9602分。固定效应的解读和普通的多元线性回归没有差别。

再看随机效应：schcode（Intercept）的方差为3.2042，标准差是1.7900，它衡量的是不同学校之间截距的变异。中间的ses的方差为0.7794，标准差为0.8828，它衡量的不同学校之间斜率的变异，意思是不同学校的ses的斜率围绕总体平均斜率变化的方差为0.7794。残差的方差为62.5855，标准差为7.9111，它衡量的是模型无法解释的变异，可以发现在考虑了随机斜率后，残差的方差又变小了（62.807到62.5855）。Corr是-1表示随机截距和随机斜率的相关系数是-1，说明有些学校的截距越高，斜率就越低（结合上面的图看，是不是有这种趋势），也就是说：随着平均数学成绩的增加，ses与数学成绩之间的关系降低。

如果想查看每个学校的截距和斜率，可以使用ranef：

head(ranef(ses_l1_random))
## $schcode
##       (Intercept)           ses
## 1    0.9746642943 -0.4806908392
## 2    1.0450460989 -0.5154021638
## 3   -3.4842479301  1.7183824946
## 4    1.8810910566 -0.9277278791
## 省略......
## 418 -0.6233121254  0.3074088487
## 419 -1.2366360507  0.6098916565

unsetunset具有交互效应的MLMunsetunset

前面分别探索了ses和public对数学成绩的影响，如果要探索它们的交互效应，只需要像普通的单水平模型一样，将交互项添加到公式中即可，比如将ses和public的交互项添加到模型中：

# 添加交互效应
crosslevel_model <- lmer(math ~ ses + public + ses:public + (ses|schcode), 
                         data = data, REML = TRUE)
## boundary (singular) fit: see help('isSingular')
# 也可以写成math ~ ses*public + (ses|schcode)
summary(crosslevel_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ ses + public + ses:public + (ses | schcode)
##    Data: data
## 
## REML criterion at convergence: 48187.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8509 -0.5593  0.1294  0.6412  5.6998 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  schcode  (Intercept)  3.2144  1.7929        
##           ses          0.8013  0.8951   -1.00
##  Residual             62.5555  7.9092        
## Number of obs: 6871, groups:  schcode, 419
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   57.72440    0.25183  382.39815 229.216   <2e-16 ***
## ses            4.42383    0.27427 1283.55623  16.130   <2e-16 ***
## public        -0.02632    0.29472  387.41741  -0.089   0.9289    
## ses:public    -0.62520    0.31957 1363.95274  -1.956   0.0506 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) ses    public
## ses        -0.232              
## public     -0.852  0.197       
## ses:public  0.198 -0.858 -0.250
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

先看固定效应（由于public是分类变量，1表示公立学校，0表示私立学校，所以在进行回归分析时会自动进行哑变量编码，以0（也就是私立学校）为参考，这里涉及一个基础知识，即回归分析中的哑变量编码）：

截距（Intercept）是57.72440，即：当学校为是私立学校（public=0）且ses也为0时的平均数学预期成绩；
ses的估计值为4.42383，即：对于私立学校（public=0）来说，ses每增加一个单位，数学成绩会增加4.42383分；
public的估计值为-0.02632，即：公立学校（public=1）相比于私立学校（public=0），平均数学成绩会减少0.02632分；
ses:public的估计值为-0.62520，即：ses对数学成绩的影响在公立学校（public=1）比在私立学校（public=0）平均减少0.62520分。
借助这些系数，我们可以估算ses在公立学校（public=1）的预期斜率，即4.42383-0.62520=3.79863。因此公立学校中ses每增加一个单位，数学成绩会增加3.79863分，略小于私立学校（4.42383分）。

再看随机效应：（和上面ses_l1_random的结果解读基本类似，这里简单说一下）。schcode（Intercept）的方差为3.2144，标准差是1.7929，它衡量的是不同学校之间截距的变异。中间的ses的方差为0.8013，标准差为0.8951，它衡量的不同学校之间斜率的变异，意思是不同学校的ses的斜率围绕总体平均斜率变化的方差为0.8013。残差的方差为62.5555，标准差为7.9092，它衡量的是模型无法解释的变异。Corr是-1表示随机截距和随机斜率的相关系数是-1。

unsetunset重复测量数据的MLMunsetunset

重复测量数据的结构非常适合多水平模型，测量的时间点可以看做是1级水平，每个患者可以看做是2级水平，每个患者都包括多次测量数据，这是一个具有两个层次的结构。

下面用一个简单的例子进行演示。

使用某药治疗10个高血压患者，分别测量每个患者治疗前后的血压，请对该数据进行分析。

# 模拟数据
data12_1 <- data.frame(id=c(1:10,1:10),
                       stat=rep(c("治疗前","治疗后"),each=10),
                       bp=c(130,124,136,128,122,118,116,138,126,124,
                            114,110,126,116,102,100,98,122,108,106)
                       )
# 设置下因子水平，变成治疗后vs治疗前
data12_1$stat <- factor(data12_1$stat,levels = c("治疗前","治疗后"))
data12_1$id <- factor(data12_1$id)
str(data12_1)
## 'data.frame':    20 obs. of  3 variables:
##  $ id  : Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ stat: Factor w/ 2 levels "治疗前","治疗后": 1 1 1 1 1 1 1 1 1 1 ...
##  $ bp  : num  130 124 136 128 122 118 116 138 126 124 ...
data12_1 # 数据长这样
##    id   stat  bp
## 1   1 治疗前 130
## 2   2 治疗前 124
## 3   3 治疗前 136
## 4   4 治疗前 128
## 5   5 治疗前 122
## 6   6 治疗前 118
## 7   7 治疗前 116
## 8   8 治疗前 138
## 9   9 治疗前 126
## 10 10 治疗前 124
## 11  1 治疗后 114
## 12  2 治疗后 110
## 13  3 治疗后 126
## 14  4 治疗后 116
## 15  5 治疗后 102
## 16  6 治疗后 100
## 17  7 治疗后  98
## 18  8 治疗后 122
## 19  9 治疗后 108
## 20 10 治疗后 106

先画个图看看不同患者治疗前后的血压值。可以发现每个患者的斜率和截距都是不一样的：

library(ggplot2)
ggplot(data12_1, aes(stat,bp))+
  geom_line(aes(color=id,group = id))

下面就是建立多水平模型即可，这里为了省事，我们默认斜率是相等的，只考虑随机截距：

library(lmerTest)

f <- lmer(bp ~ stat + (1 | id), data = data12_1)
summary(f)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: bp ~ stat + (1 | id)
##    Data: data12_1
## 
## REML criterion at convergence: 113.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1422 -0.5311  0.1307  0.4070  1.5714 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 63.511   7.969   
##  Residual              4.889   2.211   
## Number of obs: 20, groups:  id, 10
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 126.2000     2.6153   9.6662   48.25 7.56e-13 ***
## stat治疗后  -16.0000     0.9888   9.0000  -16.18 5.83e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
## Warning in abbreviate(rn, minlength = 6): abbreviate used with non-ASCII chars
##            (Intr)
## stat治疗后 -0.189

结果解读不再重复。

unsetunset广义混合效应模型unsetunset

上面介绍的例子，都是假定因变量为连续分布，而在医学和公共卫生领域，许多应变量是离散型的，例如个体的健康状态可能与吸烟、饮酒、锻炼等日常生活方式有关。在离散型应变量的情形下，若数据具有层次结构特征，则最低水平的观察单位发生某事件的概率并不完全相互独立，故不再服从二项分布或Poisson分布，而服从超二项（extra-binomial）分布或超Poisson（extra-Poisson）分布。因此，在拟合这种类型的模型时，结局的聚集效应和离散型误差的复杂分布应考虑在模型中。假定在某试验中对某事件的测量为发生与不发生，即二分类的资料，若将其作为因变量，则在多水平框架内，处理这类资料的统计模型一般称为多水平广义线性模型。

使用方法和结果解读基本类似，以后遇到了再详细介绍。