超好理解的广义线性混合模型

超好理解的广义线性混合模型

引用

CSDN

1.

https://m.blog.csdn.net/llthxx/article/details/144106541

广义线性混合模型（GLMM）是统计学和机器学习领域中一个重要的模型，它结合了广义线性模型（GLM）和混合效应模型的特点，能够处理具有层级结构和非正态分布数据的复杂问题。本文将详细介绍GLMM的基本概念、组件以及在R语言中的实现方法，并通过具体的代码示例帮助读者更好地理解和应用这一模型。

广义线性混合模型介绍

1.广义线性模型（GLM）概述

在介绍GLMM之前，我们需要先回顾一下广义线性模型（GLM）的基本概念。广义线性模型是对经典线性回归模型的推广，适用于非正态分布的数据。GLM由以下三个主要组成部分构成：

• 随机成分：响应变量的分布属于指数分布族，如正态分布、二项分布、泊松分布等。
• 系统成分：线性预测子，表示自变量与响应变量之间的线性关系。
• 连接函数：将线性预测子与响应变量的期望值联系起来的函数，如对数函数、logit函数等。

2.混合效应模型

混合效应模型（Mixed Effects Model）是一种包含固定效应和随机效应的模型：

• 固定效应（Fixed Effects）：这些效应通常是我们关心的主要效应，比如实验中的处理效应、性别、年龄等，这些因素在整个样本中被认为是相同的。
• 随机效应（Random Effects）：这些效应通常用于捕捉数据中不同层级（如个体、组、学校等）之间的变异，随机效应在不同层次（如不同个体或组）之间变化，通常认为是从某个分布（如正态分布）中随机抽取的。

混合效应模型的一个重要特征是能够处理数据中的层级结构（例如，学生嵌套在班级中、病人嵌套在医院中等），并且考虑了这些层级之间的相关性。

3.广义线性混合模型（GLMM）

广义线性混合模型结合了广义线性模型和混合效应模型，它允许响应变量服从任何符合广义线性模型设定的分布（如二项分布、泊松分布等），并且能够考虑数据中的固定效应和随机效应。GLMM的基本结构可以描述为：

• 随机成分：响应变量的分布属于指数分布族。
• 系统成分：包含固定效应和随机效应的线性预测子。
• 连接函数：将线性预测子与响应变量的期望值联系起来的函数。

4.GLMM的组件

GLMM由以下几个重要组件构成：

• 固定效应：模型中需要估计的参数，通常是我们感兴趣的效应，如处理效应、协变量等。
• 随机效应：用于捕捉数据中不同层级之间的变异，通常假设服从正态分布。
• 连接函数：将线性预测子与响应变量的期望值联系起来的函数。
• 分布假设：响应变量的分布类型，如二项分布、泊松分布等。

代码实现

下面将通过具体的R代码示例，展示如何使用lme4包拟合GLMM，并进行模型评估和结果解释。

1.拟合混合效应模型：

library(lme4)

fit33 <- glmer(response.x ~ RT.mean + CURRENT_FIX_DURATION.mean +
               fixatype + sm + fixacount + (1 | sub), data = ic, family = "binomial")

• glmer()是来自lme4包的函数，用于拟合广义线性混合效应模型（Generalized Linear Mixed Effects Model）。
• response.x：因变量（响应变量），这是一个二分类变量，适用于二项分布（family = "binomial"）。
• 固定效应（RT.mean + CURRENT_FIX_DURATION.mean + fixatype + sm + fixacount）：
• RT.mean、CURRENT_FIX_DURATION.mean等是自变量，表示不同的解释变量（例如，反应时间、当前注视时长、注视类型等）。
• 随机效应：(1 | sub)表示随机截距效应，sub表示被试（subject），即模型假设每个被试有不同的截距。
• family = "binomial"指定模型为二项式回归模型，适用于二分类结果。

2.模型摘要：

summary(fit33)

summary(fit33)输出混合效应模型的详细统计信息，包括固定效应和随机效应的估计值、标准误、z值、p值等。

3.计算R²值：

library(MuMIn)
r2_values <- r.squaredGLMM(fit33)
print(r2_values)

• r.squaredGLMM()计算广义线性混合效应模型的R²值（解释的方差比例），它能够量化模型的拟合优度，反映模型对数据的解释能力。
• MuMIn包提供了多个功能来计算模型的AIC、R²等信息。

4.模型参数的置信区间：

confint(fit33, method = "Wald")

• confint(fit33, method = "Wald")计算Wald置信区间，Wald方法通过参数的标准误来计算每个估计值的置信区间。
• 另外的method = "boot"选项可以使用自助法（Bootstrap）计算置信区间，通过多次抽样来得到置信区间。

5.固定效应和随机效应：

fixef(fit33)
ranef(fit33)

• fixef(fit33)提取固定效应参数，即解释变量（如RT.mean、CURRENT_FIX_DURATION.mean等）的回归系数。
• ranef(fit33)提取随机效应，即每个被试的随机截距。

6.计算固定效应的odds比（Odd Ratios）：

exp(fixef(fit33))

• exp(fixef(fit33))计算固定效应的odds比。通过对固定效应系数（回归系数）取指数，得到自变量变动单位对应的odds比，这有助于理解每个自变量对因变量发生概率的影响。

7.过度离势检验：

deviance(fit33) / df.residual(fit33)

• 过度离势检验（Deviance / Residual Degrees of Freedom）是模型拟合优度的一个检验方法。
• Deviance是模型的拟合残差度量，反映模型的拟合效果。
• df.residual是残差自由度，表示数据中剩余的自由度。
• 如果结果接近 1，说明模型拟合合理，过度离势问题较小。

8.ROC曲线：

ic$pre <- predict(fit33, type = "response")
library(pROC)
modelroc <- roc(ic$response.x, ic$pre)
plot(modelroc, print.auc = TRUE, auc.polygon = TRUE, grid = c(0.1, 0.2),
     grid.col = c("green", "red"), max.auc.polygon = TRUE,
     auc.polygon.col = "skyblue", print.thres = TRUE)

• 预测概率：ic$pre <- predict(fit33, type = "response")使用模型fit33来预测响应变量的概率值（即response.x为1的概率）。
• ROC曲线绘制：使用pROC包中的roc()函数生成接收者操作特征曲线（ROC Curve），并使用plot()函数进行可视化，显示AUC（曲线下面积）以及最佳阈值等信息。
• print.auc = TRUE显示AUC值。
• auc.polygon = TRUE和max.auc.polygon = TRUE用于绘制AUC区域。
• grid = c(0.1, 0.2)和grid.col = c("green", "red")用于设置网格线。
• print.thres = TRUE显示最佳的分类阈值。

plot(ci(modelroc, of = "thresholds", thresholds = "best"))
ci(modelroc)

• ci(modelroc)计算并绘制ROC曲线的置信区间，特别是最佳阈值位置。

9.混淆矩阵：

ic$pred.cutoff <- ifelse(ic$pre >= 0.286, 1, 0)
confusionMatrix(data = factor(ic$pred.cutoff),
                reference = factor(ic$response.x),
                positive = "1")

• ic$pred.cutoff <- ifelse(ic$pre >= 0.286, 1, 0)根据预测的概率值（ic$pre），使用0.286作为分类阈值，将概率大于等于 0.286 的预测为1（正类），否则预测为0（负类）。
• confusionMatrix()计算并输出混淆矩阵，比较实际标签ic$response.x和预测标签ic$pred.cutoff之间的差异。
• 混淆矩阵提供了准确率（Accuracy）、敏感性（Sensitivity）、特异性（Specificity）等指标。
• positive = "1"指定1为正类。

10.性能指标：

• 准确率（Accuracy）：所有正确预测的比例。
• 敏感性（Sensitivity）：正类（1）被正确预测的比例，也叫召回率。
• 特异性（Specificity）：负类（0）被正确预测的比例。
• Kappa：反映分类一致性的统计量，接近1表示一致性较好。
• McNemar's Test：检验二分类预测的显著性，若p-value > 0.05，说明预测结果无显著差异。

总结

# 加载所需的库
library(lme4)       # 用于混合效应模型
library(MuMIn)      # 用于计算R²值
library(pROC)       # 用于绘制ROC曲线
library(caret)      # 用于混淆矩阵和其他分类评估

# 拟合广义线性混合效应模型 (GLMM)，响应变量是二分类的，使用binomial分布
fit33 <- glmer(response.x ~ RT.mean + CURRENT_FIX_DURATION.mean +
               fixatype + sm + fixacount + (1 | sub), data = ic, family = "binomial")

# 查看模型摘要
summary(fit33)

# 计算模型的R²值
r2_values <- r.squaredGLMM(fit33)
print(r2_values)

# 计算模型的置信区间 (Wald方法)
confint(fit33, method = "Wald")

# 提取固定效应和随机效应
fixef(fit33)  # 固定效应
ranef(fit33)  # 随机效应

# 计算固定效应的odds比
exp(fixef(fit33))  # 对固定效应取指数，计算odds比

# 计算过度离势检验 (deviance / residual df)，期望值接近1
deviance(fit33) / df.residual(fit33)

# 计算预测概率并进行ROC曲线绘制
ic$pre <- predict(fit33, type = "response")

# 绘制ROC曲线
modelroc <- roc(ic$response.x, ic$pre)
plot(modelroc, 
     print.auc = TRUE, 
     auc.polygon = TRUE, 
     grid = c(0.1, 0.2),
     grid.col = c("green", "red"), 
     max.auc.polygon = TRUE,
     auc.polygon.col = "skyblue", 
     print.thres = TRUE)  # 绘制ROC曲线，显示AUC值和最佳阈值

# 绘制最佳阈值位置
plot(ci(modelroc, of = "thresholds", thresholds = "best"))

# 计算并查看ROC的置信区间
ci(modelroc)

# 根据预测概率进行分类，并计算混淆矩阵
ic$pred.cutoff <- ifelse(ic$pre >= 0.286, 1, 0)  # 设置阈值为0.286

# 混淆矩阵分析
confusionMatrix(data = factor(ic$pred.cutoff),
                reference = factor(ic$response.x),
                positive = "1")  # '1' 为正类