R语言mediation包进阶GLMM处理分类变量中介效应的完整解决方案1. 社会科学研究中的分类变量中介分析挑战在心理学、教育学等社会科学研究中我们经常遇到自变量、中介变量和因变量均为分类变量的情况。例如研究学校类型天主教/非天主教是否通过学生对学校的依恋程度喜欢/不喜欢影响校园暴力行为发生/未发生。这类研究设计需要特殊的中介效应分析方法。传统的中介效应分析通常假设连续型变量和线性关系但社会科学数据往往具有以下特征嵌套结构学生嵌套在班级/学校中分类变量二分或多分类的测量尺度非正态分布特别是当因变量为罕见事件时GLMM广义线性混合模型与mediation包的组合为解决这些问题提供了理想方案。这种组合能够处理分类变量考虑组内相关性提供准确的效应量估计2. 数据准备与变量编码2.1 示例数据集结构我们使用一个虚构的学生数据集包含以下关键变量变量名类型描述编码SCH_ID分类学校ID数值catholic二分是否天主教学校0非天主教,1天主教attachment二分学校依恋程度0不喜欢,1喜欢fight二分是否发生打架0否,1是gender二分学生性别0男,1女income有序家庭收入等级1-13级# 数据导入与查看 library(tidyverse) student - read_csv(student_data.csv) %% mutate(across(c(catholic, attachment, fight, gender), as.factor)) glimpse(student)2.2 分类变量编码原则对于中介分析中的分类变量需要特别注意编码方式二分变量建议使用0/1编码而非1/2编码确保结果解释时logit值方向正确例如relevel(factor(catholic), ref 0)多分类变量有序分类可视为连续或使用多项式模型无序分类必须创建虚拟变量注意当使用GLMM时分类预测变量的参考水平设置会直接影响结果解释。建议在分析前使用contrasts()函数检查对比矩阵。3. GLMM模型构建与验证3.1 两阶段模型设定中介分析需要建立两个模型中介变量模型预测中介变量(attachment)与自变量(catholic)的关系结果变量模型预测结果变量(fight)与自变量、中介变量的关系library(lme4) # 模型1中介变量模型 (M ~ X covariates) med_formula - attachment ~ catholic gender income (1|SCH_ID) med.fit - glmer(med_formula, family binomial(link logit), data student, control glmerControl(optimizer bobyqa)) # 模型2结果变量模型 (Y ~ X M X*M covariates) out_formula - fight ~ catholic * attachment gender income (1 attachment|SCH_ID) out.fit - glmer(out_formula, family binomial(link logit), data student, control glmerControl(optimizer bobyqa))3.2 模型诊断关键指标在运行模型后必须检查以下诊断指标收敛性确保模型已收敛检查optimizer警告尝试不同优化算法如bobyqa或Nelder_Mead奇异拟合检查随机效应方差是否接近0使用isSingular()函数可能需要简化随机效应结构多重共线性# 检查VIF值 car::vif(med.fit) car::vif(out.fit)模型比较使用ANOVA比较嵌套模型或通过AIC/BIC评估模型拟合4. 中介效应分析与结果解读4.1 mediation包实现library(mediation) set.seed(1234) # 保证结果可重复 med.out - mediate( med.fit, out.fit, treat catholic, # 处理变量 mediator attachment, # 中介变量 sims 1000, # 建议至少1000次模拟 boot TRUE, # 使用bootstrap boot.ci.type perc # 百分位置信区间 ) summary(med.out)4.2 结果解读框架典型输出包含三部分效应ACME (Average Causal Mediation Effect)中介效应间接效应表示X通过M影响Y的效应量ADE (Average Direct Effect)直接效应表示X直接影响Y的部分Total Effect总效应 ACME ADE对于分类变量效应量以概率尺度报告。解读时需注意效应方向系数符号表示效应方向统计显著性95%CI不包含0表示显著效应大小比较ACME与ADE的相对比例4.3 可视化呈现# 基础可视化 plot(med.out) # 使用ggplot2增强可视化 library(ggplot2) effects - data.frame( Effect c(ACME, ADE, Total), Estimate c(med.out$d0, med.out$z0, med.out$tau.coef), CI_lower c(med.out$d0.ci[1], med.out$z0.ci[1], med.out$tau.ci[1]), CI_upper c(med.out$d0.ci[2], med.out$z0.ci[2], med.out$tau.ci[2]) ) ggplot(effects, aes(x Effect, y Estimate)) geom_point(size 3) geom_errorbar(aes(ymin CI_lower, ymax CI_upper), width 0.1) labs(title 中介效应分解, y 效应量估计, x ) theme_minimal()5. 高级主题与问题排查5.1 常见问题解决方案模型不收敛增加迭代次数control glmerControl(optimizer bobyqa, optCtrl list(maxfun 2e5))标准化连续预测变量简化随机效应结构零膨胀问题考虑零膨胀模型使用brms包进行贝叶斯估计小样本校正使用Kenward-Roger或Satterthwaite自由度近似考虑贝叶斯方法提供更稳定的估计5.2 多类别中介分析当自变量或中介变量为多分类时# 对于三分类中介变量 med.fit - glmer(attachment_level ~ catholic (1|SCH_ID), family binomial(link logit), data student %% mutate(attachment_level case_when( attachment 0 ~ low, attachment 1 ~ medium, TRUE ~ high) %% ordered())) # 需要使用multinomial模型 library(nnet) med.fit - multinom(attachment_level ~ catholic, data student)5.3 敏感性分析评估未测量混杂因素的影响# 使用mediation包的sens参数 med.sens - medsens(med.out, rho.by 0.1, effect.type indirect) plot(med.sens)6. 完整案例代码与数据获取6.1 完整分析流程# 加载必要包 library(tidyverse) library(lme4) library(mediation) # 数据准备 student - read_csv(student_data.csv) %% mutate(across(c(catholic, attachment, fight, gender), ~ factor(.x, levels c(0, 1), labels c(No, Yes)))) # 模型构建 med.fit - glmer(attachment ~ catholic gender income (1|SCH_ID), family binomial, data student, control glmerControl(optimizer bobyqa)) out.fit - glmer(fight ~ catholic * attachment gender income (1|SCH_ID), family binomial, data student, control glmerControl(optimizer bobyqa)) # 中介分析 med.out - mediate(med.fit, out.fit, treat catholic, mediator attachment, sims 1000, boot TRUE) # 结果输出 summary(med.out) plot(med.out) # 效应量转换 effects - summary(med.out) logit_to_prob - function(x) exp(x)/(1exp(x)) data.frame( Effect c(ACME, ADE, Total), Logit_Scale c(effects$d0, effects$z0, effects$tau.coef), Probability_Scale logit_to_prob(c(effects$d0, effects$z0, effects$tau.coef)) )6.2 模拟数据生成如果无法获取真实数据可以使用以下代码生成模拟数据集set.seed(123) n_schools - 20 n_students - 200 sim_data - data.frame( SCH_ID rep(1:n_schools, each n_students/n_schools), catholic rbinom(n_students, 1, 0.5), gender rbinom(n_students, 1, 0.5), income sample(1:13, n_students, replace TRUE) ) %% group_by(SCH_ID) %% mutate( # 学校随机效应 school_effect rnorm(1, 0, 0.5), # 中介变量模型 logit_attachment -0.5 0.8*catholic 0.3*gender - 0.1*income school_effect, prob_attachment plogis(logit_attachment), attachment rbinom(n(), 1, prob_attachment), # 结果变量模型 logit_fight -1 0.5*catholic - 0.8*attachment 0.4*catholic*attachment 0.2*gender - 0.05*income school_effect, prob_fight plogis(logit_fight), fight rbinom(n(), 1, prob_fight) ) %% ungroup() write_csv(sim_data, student_sim_data.csv)
R语言mediation包实战：如何用GLMM处理分类变量的中介效应分析（附学生数据集）
R语言mediation包进阶GLMM处理分类变量中介效应的完整解决方案1. 社会科学研究中的分类变量中介分析挑战在心理学、教育学等社会科学研究中我们经常遇到自变量、中介变量和因变量均为分类变量的情况。例如研究学校类型天主教/非天主教是否通过学生对学校的依恋程度喜欢/不喜欢影响校园暴力行为发生/未发生。这类研究设计需要特殊的中介效应分析方法。传统的中介效应分析通常假设连续型变量和线性关系但社会科学数据往往具有以下特征嵌套结构学生嵌套在班级/学校中分类变量二分或多分类的测量尺度非正态分布特别是当因变量为罕见事件时GLMM广义线性混合模型与mediation包的组合为解决这些问题提供了理想方案。这种组合能够处理分类变量考虑组内相关性提供准确的效应量估计2. 数据准备与变量编码2.1 示例数据集结构我们使用一个虚构的学生数据集包含以下关键变量变量名类型描述编码SCH_ID分类学校ID数值catholic二分是否天主教学校0非天主教,1天主教attachment二分学校依恋程度0不喜欢,1喜欢fight二分是否发生打架0否,1是gender二分学生性别0男,1女income有序家庭收入等级1-13级# 数据导入与查看 library(tidyverse) student - read_csv(student_data.csv) %% mutate(across(c(catholic, attachment, fight, gender), as.factor)) glimpse(student)2.2 分类变量编码原则对于中介分析中的分类变量需要特别注意编码方式二分变量建议使用0/1编码而非1/2编码确保结果解释时logit值方向正确例如relevel(factor(catholic), ref 0)多分类变量有序分类可视为连续或使用多项式模型无序分类必须创建虚拟变量注意当使用GLMM时分类预测变量的参考水平设置会直接影响结果解释。建议在分析前使用contrasts()函数检查对比矩阵。3. GLMM模型构建与验证3.1 两阶段模型设定中介分析需要建立两个模型中介变量模型预测中介变量(attachment)与自变量(catholic)的关系结果变量模型预测结果变量(fight)与自变量、中介变量的关系library(lme4) # 模型1中介变量模型 (M ~ X covariates) med_formula - attachment ~ catholic gender income (1|SCH_ID) med.fit - glmer(med_formula, family binomial(link logit), data student, control glmerControl(optimizer bobyqa)) # 模型2结果变量模型 (Y ~ X M X*M covariates) out_formula - fight ~ catholic * attachment gender income (1 attachment|SCH_ID) out.fit - glmer(out_formula, family binomial(link logit), data student, control glmerControl(optimizer bobyqa))3.2 模型诊断关键指标在运行模型后必须检查以下诊断指标收敛性确保模型已收敛检查optimizer警告尝试不同优化算法如bobyqa或Nelder_Mead奇异拟合检查随机效应方差是否接近0使用isSingular()函数可能需要简化随机效应结构多重共线性# 检查VIF值 car::vif(med.fit) car::vif(out.fit)模型比较使用ANOVA比较嵌套模型或通过AIC/BIC评估模型拟合4. 中介效应分析与结果解读4.1 mediation包实现library(mediation) set.seed(1234) # 保证结果可重复 med.out - mediate( med.fit, out.fit, treat catholic, # 处理变量 mediator attachment, # 中介变量 sims 1000, # 建议至少1000次模拟 boot TRUE, # 使用bootstrap boot.ci.type perc # 百分位置信区间 ) summary(med.out)4.2 结果解读框架典型输出包含三部分效应ACME (Average Causal Mediation Effect)中介效应间接效应表示X通过M影响Y的效应量ADE (Average Direct Effect)直接效应表示X直接影响Y的部分Total Effect总效应 ACME ADE对于分类变量效应量以概率尺度报告。解读时需注意效应方向系数符号表示效应方向统计显著性95%CI不包含0表示显著效应大小比较ACME与ADE的相对比例4.3 可视化呈现# 基础可视化 plot(med.out) # 使用ggplot2增强可视化 library(ggplot2) effects - data.frame( Effect c(ACME, ADE, Total), Estimate c(med.out$d0, med.out$z0, med.out$tau.coef), CI_lower c(med.out$d0.ci[1], med.out$z0.ci[1], med.out$tau.ci[1]), CI_upper c(med.out$d0.ci[2], med.out$z0.ci[2], med.out$tau.ci[2]) ) ggplot(effects, aes(x Effect, y Estimate)) geom_point(size 3) geom_errorbar(aes(ymin CI_lower, ymax CI_upper), width 0.1) labs(title 中介效应分解, y 效应量估计, x ) theme_minimal()5. 高级主题与问题排查5.1 常见问题解决方案模型不收敛增加迭代次数control glmerControl(optimizer bobyqa, optCtrl list(maxfun 2e5))标准化连续预测变量简化随机效应结构零膨胀问题考虑零膨胀模型使用brms包进行贝叶斯估计小样本校正使用Kenward-Roger或Satterthwaite自由度近似考虑贝叶斯方法提供更稳定的估计5.2 多类别中介分析当自变量或中介变量为多分类时# 对于三分类中介变量 med.fit - glmer(attachment_level ~ catholic (1|SCH_ID), family binomial(link logit), data student %% mutate(attachment_level case_when( attachment 0 ~ low, attachment 1 ~ medium, TRUE ~ high) %% ordered())) # 需要使用multinomial模型 library(nnet) med.fit - multinom(attachment_level ~ catholic, data student)5.3 敏感性分析评估未测量混杂因素的影响# 使用mediation包的sens参数 med.sens - medsens(med.out, rho.by 0.1, effect.type indirect) plot(med.sens)6. 完整案例代码与数据获取6.1 完整分析流程# 加载必要包 library(tidyverse) library(lme4) library(mediation) # 数据准备 student - read_csv(student_data.csv) %% mutate(across(c(catholic, attachment, fight, gender), ~ factor(.x, levels c(0, 1), labels c(No, Yes)))) # 模型构建 med.fit - glmer(attachment ~ catholic gender income (1|SCH_ID), family binomial, data student, control glmerControl(optimizer bobyqa)) out.fit - glmer(fight ~ catholic * attachment gender income (1|SCH_ID), family binomial, data student, control glmerControl(optimizer bobyqa)) # 中介分析 med.out - mediate(med.fit, out.fit, treat catholic, mediator attachment, sims 1000, boot TRUE) # 结果输出 summary(med.out) plot(med.out) # 效应量转换 effects - summary(med.out) logit_to_prob - function(x) exp(x)/(1exp(x)) data.frame( Effect c(ACME, ADE, Total), Logit_Scale c(effects$d0, effects$z0, effects$tau.coef), Probability_Scale logit_to_prob(c(effects$d0, effects$z0, effects$tau.coef)) )6.2 模拟数据生成如果无法获取真实数据可以使用以下代码生成模拟数据集set.seed(123) n_schools - 20 n_students - 200 sim_data - data.frame( SCH_ID rep(1:n_schools, each n_students/n_schools), catholic rbinom(n_students, 1, 0.5), gender rbinom(n_students, 1, 0.5), income sample(1:13, n_students, replace TRUE) ) %% group_by(SCH_ID) %% mutate( # 学校随机效应 school_effect rnorm(1, 0, 0.5), # 中介变量模型 logit_attachment -0.5 0.8*catholic 0.3*gender - 0.1*income school_effect, prob_attachment plogis(logit_attachment), attachment rbinom(n(), 1, prob_attachment), # 结果变量模型 logit_fight -1 0.5*catholic - 0.8*attachment 0.4*catholic*attachment 0.2*gender - 0.05*income school_effect, prob_fight plogis(logit_fight), fight rbinom(n(), 1, prob_fight) ) %% ungroup() write_csv(sim_data, student_sim_data.csv)
