广义线性混合模型介绍
1. 广义线性模型（GLM）概述
在介绍GLMM之前，我们需要先回顾一下广义线性模型（GLM）的基本概念。广义线性模型是对经典线性回归模型的推广，适用于非正态分布的数据。GLM由以下三个主要组成部分构成：
2. 混合效应模型
混合效应模型（Mixed Effects Model）是一种包含固定效应和随机效应的模型：
- 固定效应（Fixed Effects）：这些效应通常是我们关心的主要效应，比如实验中的处理效应、性别、年龄等，这些因素在整个样本中被认为是相同的。
- 随机效应（Random Effects）：这些效应通常用于捕捉数据中不同层级（如个体、组、学校等）之间的变异，随机效应在不同层次（如不同个体或组）之间变化，通常认为是从某个分布（如正态分布）中随机抽取的。
混合效应模型的一个重要特征是能够处理数据中的层级结构（例如，学生嵌套在班级中、病人嵌套在医院中等），并且考虑了这些层级之间的相关性。
3. 广义线性混合模型（GLMM）
广义线性混合模型结合了广义线性模型和混合效应模型，它允许响应变量服从任何符合广义线性模型设定的分布（如二项分布、泊松分布等），并且能够考虑数据中的固定效应和随机效应。GLMM的基本结构可以描述为：
4. GLMM的组件
GLMM由以下几个重要组件构成：
代码
1. 拟合混合效应模型：
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' 为正类
