复现分析：情绪调节传染驱动群体间负面情绪的减少
复现研究
2026-01-19
1 引言
本R Markdown文档复现了以下论文的关键分析：
“Emotion regulation contagion drives reduction in negative intergroup emotions”（情绪调节传染驱动群体间负面情绪的减少），作者：Pinus等人（2025），发表于Nature Communications。
原研究探讨了情绪调节干预是否能在群体中传播（“情绪调节传染”），即使未接受干预的参与者也能减少负面情绪。本复现使用语义文本相似度数据来检验：
- 认知重评干预对情绪评分的影响
- 接受干预参与者比例与情绪结果之间的关系
- 通过语义分析证明情绪调节传染的存在
- 重评语言从干预组向非干预组的传播
2 设置与数据加载
2.1 加载所需R包
# 如果尚未安装，则安装所需包
packages <- c(
"readxl",
"tidyverse",
"lme4",
"lmerTest",
"performance",
"effectsize",
"emmeans",
"ggeffects",
"patchwork",
"scales",
"knitr",
"kableExtra",
"car",
"nortest",
"robustlmm"
)
# 安装缺失的包
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages, repos = "https://cloud.r-project.org")
# 加载包
invisible(lapply(packages, library, character.only = TRUE))
# 设置所有图形的主题
theme_set(theme_minimal(base_size = 12) +
theme(
panel.grid.minor = element_blank(),
legend.position = "bottom"
))2.2 加载数据
# 从桌面加载数据（MacOS路径）
# 根据您的系统调整路径
data_path <- "~/Desktop/Semantic Textual Similarity.xlsx"
# 读取Excel文件
df <- read_excel(data_path, sheet = "Sheet1")
# 显示数据结构
glimpse(df)## Rows: 3,200
## Columns: 9
## $ 被试编号 <dbl> 56, 56, 56, 56, 56, 56, 56, 56, 5…
## $ 组别 <chr> "1-1", "1-1", "1-1", "1-1", "1-1"…
## $ 干预比例 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ `亲密度（0=0；1=1）` <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ `实验条件（0=自然观看；认知重评=1）` <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ 试次 <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11…
## $ 文本 <chr> "一个人靠着墙闭眼坐在地上。", "让人感到害怕。", "被别人打…
## $ 语义相似度 <dbl> 0.31062492, 0.18147468, 0.2632054…
## $ 情绪评分 <dbl> 4, 2, 2, 2, 2, 1, 1, 2, 2, 3, 3, …2.3 数据清理与变量准备
# 将列名重命名为英文以便于分析
df <- df %>%
rename(
participant_id = `被试编号`,
group_id = `组别`,
intervention_proportion = `干预比例`,
closeness = `亲密度（0=0；1=1）`,
condition = `实验条件（0=自然观看；认知重评=1）`,
trial = `试次`,
text = `文本`,
semantic_similarity = `语义相似度`,
emotion_rating = `情绪评分`
)
# 转换为适当的数据类型
df <- df %>%
mutate(
participant_id = as.factor(participant_id),
group_id = as.factor(group_id),
condition = factor(condition, levels = c(0, 1),
labels = c("Non-treated", "Treated")),
closeness = factor(closeness, levels = c(0, 1),
labels = c("Low", "High")),
intervention_proportion = as.numeric(intervention_proportion),
trial = as.integer(trial),
semantic_similarity = as.numeric(semantic_similarity),
emotion_rating = as.numeric(emotion_rating)
)
# 创建百分比形式的比例变量以便于解释
df <- df %>%
mutate(
proportion_pct = intervention_proportion * 100,
proportion_factor = factor(proportion_pct,
levels = c(0, 20, 40, 60),
labels = c("0%", "20%", "40%", "60%"))
)
# 显示摘要
summary(df)## participant_id group_id intervention_proportion closeness
## 1 : 20 1-1 : 100 Min. :0.00 Low :1600
## 2 : 20 1-2 : 100 1st Qu.:0.15 High:1600
## 3 : 20 1-3 : 100 Median :0.30
## 4 : 20 1-4 : 100 Mean :0.30
## 5 : 20 2-1 : 100 3rd Qu.:0.45
## 6 : 20 2-2 : 100 Max. :0.60
## (Other):3080 (Other):2600
## condition trial text semantic_similarity
## Non-treated:2240 Min. : 1.00 Length:3200 Min. :0.0000
## Treated : 960 1st Qu.: 5.75 Class :character 1st Qu.:0.3676
## Median :10.50 Mode :character Median :0.4825
## Mean :10.50 Mean :0.4864
## 3rd Qu.:15.25 3rd Qu.:0.5945
## Max. :20.00 Max. :1.0000
##
## emotion_rating proportion_pct proportion_factor
## Min. :1.000 Min. : 0 0% :800
## 1st Qu.:2.000 1st Qu.:15 20%:800
## Median :4.000 Median :30 40%:800
## Mean :3.846 Mean :30 60%:800
## 3rd Qu.:5.000 3rd Qu.:45
## Max. :9.000 Max. :60
## 2.4 数据概览
## === 样本概览 ===## 总观测数: 3200## 独立参与者数: 160## 独立组数: 32## 每位参与者的试次数: 20##
## === 条件分布 ===##
## Non-treated Treated
## 2240 960##
## === 干预比例分布 ===##
## 0% 20% 40% 60%
## 800 800 800 800##
## === 条件 able 干预比例 ===##
## 0% 20% 40% 60%
## Non-treated 800 640 480 320
## Treated 0 160 320 4803 描述性统计
3.1 按条件和比例的汇总统计
# 按条件和比例的情绪评分
emotion_summary <- df %>%
group_by(proportion_factor, condition) %>%
summarise(
n = n(),
mean_emotion = mean(emotion_rating, na.rm = TRUE),
sd_emotion = sd(emotion_rating, na.rm = TRUE),
se_emotion = sd_emotion / sqrt(n),
mean_similarity = mean(semantic_similarity, na.rm = TRUE),
sd_similarity = sd(semantic_similarity, na.rm = TRUE),
se_similarity = sd_similarity / sqrt(n),
.groups = "drop"
)
kable(emotion_summary,
digits = 3,
caption = "按条件和干预比例的描述性统计",
col.names = c("比例", "条件", "样本量",
"情绪均值", "情绪标准差", "情绪标准误",
"相似度均值", "相似度标准差", "相似度标准误")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| 比例 |条 |
样本
| |||||||
|---|---|---|---|---|---|---|---|---|
| 0% | Non-treated | 800 | 2.421 | 1.331 | 0.047 | 0.349 | 0.134 | 0.005 |
| 20% | Non-treated | 640 | 3.698 | 1.778 | 0.070 | 0.476 | 0.169 | 0.007 |
| 20% | Treated | 160 | 4.144 | 1.711 | 0.135 | 0.572 | 0.136 | 0.011 |
| 40% | Non-treated | 480 | 4.077 | 1.354 | 0.062 | 0.514 | 0.171 | 0.008 |
| 40% | Treated | 320 | 4.509 | 1.438 | 0.080 | 0.590 | 0.146 | 0.008 |
| 60% | Non-treated | 320 | 4.613 | 1.554 | 0.087 | 0.514 | 0.155 | 0.009 |
| 60% | Treated | 480 | 5.135 | 1.568 | 0.072 | 0.586 | 0.145 | 0.007 |
3.2 可视化：按条件和比例的情绪评分
# 绘制情绪评分
p_emotion <- ggplot(emotion_summary,
aes(x = proportion_factor, y = mean_emotion,
color = condition, group = condition)) +
geom_line(linewidth = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = mean_emotion - se_emotion,
ymax = mean_emotion + se_emotion),
width = 0.1, linewidth = 0.8) +
scale_color_manual(values = c("Non-treated" = "#2E86AB", "Treated" = "#E94F37")) +
labs(
title = "Negative Emotion Ratings by Intervention Proportion",
x = "Proportion of Treated Participants",
y = "Negative Emotion Rating",
color = "Condition"
) +
theme(legend.position = "bottom")
print(p_emotion)
3.3 可视化：按条件和比例的语义相似度
# 绘制语义相似度
p_similarity <- ggplot(emotion_summary,
aes(x = proportion_factor, y = mean_similarity,
color = condition, group = condition)) +
geom_line(linewidth = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = mean_similarity - se_similarity,
ymax = mean_similarity + se_similarity),
width = 0.1, linewidth = 0.8) +
scale_color_manual(values = c("Non-treated" = "#2E86AB", "Treated" = "#E94F37")) +
labs(
title = "Semantic Similarity to Reappraisal by Intervention Proportion",
subtitle = "Higher values indicate greater similarity to pure reappraisal language",
x = "Proportion of Treated Participants",
y = "Similarity to Reappraisal",
color = "Condition"
) +
theme(legend.position = "bottom")
print(p_similarity)
4 主要分析
4.1 分析1：认知重评对情绪评分的影响
本分析检验认知重评干预是否能减少负面情绪，复现原研究中的基本干预效应。
4.1.1 正态性检验
# Kolmogorov-Smirnov正态性检验
ks_result <- ks.test(df$emotion_rating, "pnorm",
mean = mean(df$emotion_rating),
sd = sd(df$emotion_rating))
cat("Kolmogorov-Smirnov test for emotion ratings:\n")## Kolmogorov-Smirnov test for emotion ratings:## D = 0.1335 , p = < 2.22e-16## Anderson-Darling test:## A = 51.5679 , p = < 2.22e-16# 可视化检验
par(mfrow = c(1, 2))
hist(df$emotion_rating, breaks = 20, main = "Distribution of Emotion Ratings",
xlab = "Emotion Rating", col = "lightblue")
qqnorm(df$emotion_rating, main = "Q-Q Plot")
qqline(df$emotion_rating, col = "red")
4.1.2 线性混合模型：基本干预效应
# 模型1：条件对情绪评分的基本效应
# 按照原论文的方法使用随机截距
model_basic <- lmer(
emotion_rating ~ condition +
(1 | group_id) +
(1 | participant_id) +
(1 | trial),
data = df,
REML = TRUE
)
# 模型摘要
summary(model_basic)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: emotion_rating ~ condition + (1 | group_id) + (1 | participant_id) +
## (1 | trial)
## Data: df
##
## REML criterion at convergence: 10363.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9691 -0.6659 -0.0151 0.6515 3.8600
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_id (Intercept) 0.5597 0.7481
## group_id (Intercept) 1.0368 1.0182
## trial (Intercept) 0.1041 0.3226
## Residual 1.2820 1.1323
## Number of obs: 3200, groups: participant_id, 160; group_id, 32; trial, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.6711 0.2089 41.5482 17.577 < 2e-16 ***
## conditionTreated 0.5838 0.1535 134.2574 3.802 0.000217 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## conditnTrtd -0.221##
## === Effect Sizes ===effectsize::standardize_parameters(model_basic) %>%
kable(digits = 3) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Parameter | Std_Coefficient | CI | CI_low | CI_high |
|---|---|---|---|---|
| (Intercept) | -0.098 | 0.95 | -0.328 | 0.131 |
| conditionTreated | 0.327 | 0.95 | 0.158 | 0.496 |
4.1.3 包含干预比例的模型
# 模型2：干预比例对情绪评分的影响
model_proportion <- lmer(
emotion_rating ~ condition * intervention_proportion +
(1 | group_id) +
(1 | participant_id) +
(1 | trial),
data = df,
REML = TRUE
)
summary(model_proportion)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## emotion_rating ~ condition * intervention_proportion + (1 | group_id) +
## (1 | participant_id) + (1 | trial)
## Data: df
##
## REML criterion at convergence: 10337.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9891 -0.6628 -0.0166 0.6529 3.8556
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_id (Intercept) 0.5622 0.7498
## group_id (Intercept) 0.4548 0.6744
## trial (Intercept) 0.1040 0.3226
## Residual 1.2820 1.1323
## Number of obs: 3200, groups: participant_id, 160; group_id, 32; trial, 20
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 2.6434 0.2397 37.1943 11.030
## conditionTreated 0.5352 0.4370 142.4772 1.225
## intervention_proportion 3.5462 0.6444 38.3697 5.503
## conditionTreated:intervention_proportion -0.1544 0.9602 144.5449 -0.161
## Pr(>|t|)
## (Intercept) 2.77e-13 ***
## conditionTreated 0.223
## intervention_proportion 2.66e-06 ***
## conditionTreated:intervention_proportion 0.872
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cndtnT intrv_
## conditnTrtd -0.150
## intrvntn_pr -0.746 0.157
## cndtnTrtd:_ 0.160 -0.934 -0.2614.2 分析2：剂量-反应关系
本分析检验干预参与者比例与情绪减少之间的关系，测试不同的函数形式（线性、二次、三次、指数）。
4.2.1 模型比较：线性 vs. 二次 vs. 三次 vs. 指数
# 准备模型比较的数据
# 对多项式模型进行比例中心化
df <- df %>%
mutate(
prop_centered = scale(intervention_proportion, scale = FALSE),
prop_scaled = scale(intervention_proportion)
)
# 线性模型
model_linear <- lmer(
emotion_rating ~ intervention_proportion + condition +
(1 | group_id) + (1 | participant_id) + (1 | trial),
data = df, REML = FALSE
)
# 二次模型
model_quadratic <- lmer(
emotion_rating ~ poly(intervention_proportion, 2, raw = TRUE) + condition +
(1 | group_id) + (1 | participant_id) + (1 | trial),
data = df, REML = FALSE
)
# 三次模型
model_cubic <- lmer(
emotion_rating ~ poly(intervention_proportion, 3, raw = TRUE) + condition +
(1 | group_id) + (1 | participant_id) + (1 | trial),
data = df, REML = FALSE
)
# 指数变换
df$prop_exp <- exp(df$intervention_proportion) - 1
model_exponential <- lmer(
emotion_rating ~ prop_exp + condition +
(1 | group_id) + (1 | participant_id) + (1 | trial),
data = df, REML = FALSE
)
# 使用AIC进行模型比较
model_comparison <- data.frame(
Model = c("Linear", "Quadratic", "Cubic", "Exponential"),
AIC = c(AIC(model_linear), AIC(model_quadratic),
AIC(model_cubic), AIC(model_exponential)),
BIC = c(BIC(model_linear), BIC(model_quadratic),
BIC(model_cubic), BIC(model_exponential))
)
model_comparison <- model_comparison %>%
mutate(
delta_AIC = AIC - min(AIC),
delta_BIC = BIC - min(BIC)
) %>%
arrange(AIC)
kable(model_comparison,
digits = 2,
caption = "Model Comparison: Different Dose-Response Functions") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Model | AIC | BIC | delta_AIC | delta_BIC |
|---|---|---|---|---|
| Linear | 10350.31 | 10392.81 | 0.00 | 0.00 |
| Quadratic | 10350.53 | 10399.10 | 0.21 | 6.29 |
| Cubic | 10351.55 | 10406.19 | 1.24 | 13.38 |
| Exponential | 10351.90 | 10394.39 | 1.58 | 1.58 |
4.2.2 最佳拟合模型摘要
# 基于AIC选择最佳模型（通常如原论文所述为指数模型）
best_model_name <- model_comparison$Model[1]
cat("Best fitting model based on AIC:", best_model_name, "\n\n")## Best fitting model based on AIC: Linear## === Quadratic Model Summary ===## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: emotion_rating ~ poly(intervention_proportion, 2, raw = TRUE) +
## condition + (1 | group_id) + (1 | participant_id) + (1 | trial)
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 10350.5 10399.1 -5167.3 10334.5 3192
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9830 -0.6631 -0.0172 0.6536 3.8481
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_id (Intercept) 0.5523 0.7432
## group_id (Intercept) 0.3953 0.6287
## trial (Intercept) 0.1025 0.3202
## Residual 1.2820 1.1323
## Number of obs: 3200, groups: participant_id, 160; group_id, 32; trial, 20
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 2.4771 0.2583 36.9242
## poly(intervention_proportion, 2, raw = TRUE)1 6.1067 1.9986 32.3282
## poly(intervention_proportion, 2, raw = TRUE)2 -4.3125 3.1824 31.9399
## conditionTreated 0.4695 0.1552 127.9993
## t value Pr(>|t|)
## (Intercept) 9.591 1.44e-11 ***
## poly(intervention_proportion, 2, raw = TRUE)1 3.056 0.00448 **
## poly(intervention_proportion, 2, raw = TRUE)2 -1.355 0.18489
## conditionTreated 3.026 0.00300 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(_,2,r=TRUE)1 p(_,2,r=TRUE)2
## p(_,2,r=TRUE)1 -0.659
## p(_,2,r=TRUE)2 0.493 -0.955
## conditnTrtd 0.000 -0.078 0.0004.2.3 可视化剂量-反应关系
# 从每个模型获取预测值
prop_seq <- seq(0, 0.6, length.out = 100)
pred_df <- expand.grid(
intervention_proportion = prop_seq,
condition = factor(c("Non-treated", "Treated"),
levels = c("Non-treated", "Treated"))
)
pred_df$prop_exp <- exp(pred_df$intervention_proportion) - 1
# 预测（仅使用固定效应进行可视化）
# 使用汇总数据以获得更清晰的可视化
agg_data <- df %>%
group_by(intervention_proportion, condition) %>%
summarise(
mean_emotion = mean(emotion_rating),
se_emotion = sd(emotion_rating) / sqrt(n()),
n = n(),
.groups = "drop"
)
# 绘制带平滑线的图形
p_dose <- ggplot(agg_data, aes(x = intervention_proportion, y = mean_emotion,
color = condition)) +
geom_point(size = 4, alpha = 0.8) +
geom_errorbar(aes(ymin = mean_emotion - se_emotion,
ymax = mean_emotion + se_emotion),
width = 0.02, alpha = 0.8) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, alpha = 0.2) +
scale_color_manual(values = c("Non-treated" = "#2E86AB", "Treated" = "#E94F37")) +
scale_x_continuous(labels = scales::percent_format(),
breaks = c(0, 0.2, 0.4, 0.6)) +
labs(
title = "Dose-Response: Effect of Intervention Proportion on Negative Emotions",
subtitle = "Quadratic fit with 95% confidence bands",
x = "Proportion of Treated Participants",
y = "Negative Emotion Rating",
color = "Condition"
) +
theme(legend.position = "bottom")
print(p_dose)
4.3 分析3：干预组与非干预组的分别分析
按照原论文的方法，我们分别分析干预组和非干预组参与者。
4.3.1 仅非干预组参与者
# 筛选非干预组参与者
df_nontreated <- df %>% filter(condition == "Non-treated")
cat("=== Analysis of Non-Treated Participants ===\n")## === Analysis of Non-Treated Participants ===## N observations: 2240# 非干预组参与者的模型
model_nontreated_quad <- lmer(
emotion_rating ~ poly(intervention_proportion, 2, raw = TRUE) +
(1 | group_id) + (1 | participant_id) + (1 | trial),
data = df_nontreated,
REML = TRUE
)
summary(model_nontreated_quad)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: emotion_rating ~ poly(intervention_proportion, 2, raw = TRUE) +
## (1 | group_id) + (1 | participant_id) + (1 | trial)
## Data: df_nontreated
##
## REML criterion at convergence: 7137.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9556 -0.6591 -0.0196 0.6229 3.8832
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_id (Intercept) 0.4040 0.6356
## group_id (Intercept) 0.5852 0.7650
## trial (Intercept) 0.1061 0.3258
## Residual 1.2294 1.1088
## Number of obs: 2240, groups: participant_id, 112; group_id, 32; trial, 20
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 2.4707 0.2935 29.4139
## poly(intervention_proportion, 2, raw = TRUE)1 6.3219 2.3336 28.4518
## poly(intervention_proportion, 2, raw = TRUE)2 -4.7527 3.7848 30.2115
## t value Pr(>|t|)
## (Intercept) 8.419 2.51e-09 ***
## poly(intervention_proportion, 2, raw = TRUE)1 2.709 0.0113 *
## poly(intervention_proportion, 2, raw = TRUE)2 -1.256 0.2188
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(_,2,r=TRUE)1
## p(_,2,r=TRUE)1 -0.658
## p(_,2,r=TRUE)2 0.486 -0.956# 比例对非干预组参与者情绪的影响
cat("\n=== This tests whether increasing the proportion of treated participants\n")##
## === This tests whether increasing the proportion of treated participants## reduces emotions in NON-TREATED participants (emotion regulation contagion) ===4.3.2 仅干预组参与者
# 筛选干预组参与者
df_treated <- df %>% filter(condition == "Treated")
cat("=== Analysis of Treated Participants ===\n")## === Analysis of Treated Participants ===## N observations: 960# 干预组参与者的模型
model_treated_quad <- lmer(
emotion_rating ~ poly(intervention_proportion, 2, raw = TRUE) +
(1 | group_id) + (1 | participant_id) + (1 | trial),
data = df_treated,
REML = TRUE
)
summary(model_treated_quad)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: emotion_rating ~ poly(intervention_proportion, 2, raw = TRUE) +
## (1 | group_id) + (1 | participant_id) + (1 | trial)
## Data: df_treated
##
## REML criterion at convergence: 3196.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6209 -0.6841 0.0054 0.6712 3.3673
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_id (Intercept) 0.57391 0.7576
## group_id (Intercept) 0.42719 0.6536
## trial (Intercept) 0.09596 0.3098
## Residual 1.40815 1.1867
## Number of obs: 960, groups: participant_id, 48; group_id, 24; trial, 20
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 4.038 1.461 27.515
## poly(intervention_proportion, 2, raw = TRUE)1 -0.125 7.931 23.469
## poly(intervention_proportion, 2, raw = TRUE)2 3.255 9.592 21.509
## t value Pr(>|t|)
## (Intercept) 2.765 0.010 *
## poly(intervention_proportion, 2, raw = TRUE)1 -0.016 0.988
## poly(intervention_proportion, 2, raw = TRUE)2 0.339 0.738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(_,2,r=TRUE)1
## p(_,2,r=TRUE)1 -0.971
## p(_,2,r=TRUE)2 0.931 -0.9904.3.3 可视化分组效应
# 创建分组可视化
agg_nontreated <- df_nontreated %>%
group_by(intervention_proportion) %>%
summarise(
mean_emotion = mean(emotion_rating),
se_emotion = sd(emotion_rating) / sqrt(n()),
n = n(),
.groups = "drop"
)
agg_treated <- df_treated %>%
group_by(intervention_proportion) %>%
summarise(
mean_emotion = mean(emotion_rating),
se_emotion = sd(emotion_rating) / sqrt(n()),
n = n(),
.groups = "drop"
)
# 非干预组图
p_nontreated <- ggplot(agg_nontreated,
aes(x = intervention_proportion, y = mean_emotion)) +
geom_point(size = 4, color = "#2E86AB") +
geom_errorbar(aes(ymin = mean_emotion - se_emotion,
ymax = mean_emotion + se_emotion),
width = 0.02, color = "#2E86AB") +
geom_smooth(method = "lm", formula = y ~ poly(x, 2),
se = TRUE, color = "#2E86AB", fill = "#2E86AB", alpha = 0.2) +
scale_x_continuous(labels = scales::percent_format()) +
labs(
title = "Non-Treated Participants",
subtitle = "Evidence for emotion regulation contagion",
x = "Proportion of Treated Participants",
y = "Negative Emotion Rating"
)
# 干预组图
p_treated <- ggplot(agg_treated,
aes(x = intervention_proportion, y = mean_emotion)) +
geom_point(size = 4, color = "#E94F37") +
geom_errorbar(aes(ymin = mean_emotion - se_emotion,
ymax = mean_emotion + se_emotion),
width = 0.02, color = "#E94F37") +
geom_smooth(method = "lm", formula = y ~ poly(x, 2),
se = TRUE, color = "#E94F37", fill = "#E94F37", alpha = 0.2) +
scale_x_continuous(labels = scales::percent_format()) +
labs(
title = "Treated Participants",
subtitle = "Direct effect of reappraisal training",
x = "Proportion of Treated Participants",
y = "Negative Emotion Rating"
)
# 组合图形
p_nontreated + p_treated +
plot_annotation(
title = "Emotion Ratings by Condition: Evidence for Contagion Effect",
theme = theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"))
)
4.4 分析4：语义相似度分析（复现语义投影分析）
本分析检验非干预组参与者是否从干预组参与者那里习得重评语言，为社会学习机制提供证据。
4.4.1 模型：比例对语义相似度的影响
# 包含交互作用的完整模型
model_semantic <- lmer(
semantic_similarity ~ intervention_proportion * condition * trial +
(1 | group_id) + (1 | participant_id),
data = df,
REML = TRUE
)
summary(model_semantic)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: semantic_similarity ~ intervention_proportion * condition * trial +
## (1 | group_id) + (1 | participant_id)
## Data: df
##
## REML criterion at convergence: -3265.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8836 -0.6802 -0.0526 0.6305 3.5731
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_id (Intercept) 0.002201 0.04691
## group_id (Intercept) 0.002211 0.04702
## Residual 0.019225 0.13865
## Number of obs: 3200, groups: participant_id, 160; group_id, 32
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 3.422e-01 1.783e-02 4.618e+01
## intervention_proportion 2.254e-01 5.174e-02 6.353e+01
## conditionTreated 2.110e-01 4.172e-02 4.517e+02
## trial 3.646e-03 7.473e-04 3.036e+03
## intervention_proportion:conditionTreated -9.685e-02 9.091e-02 4.447e+02
## intervention_proportion:trial 5.083e-03 2.398e-03 3.036e+03
## conditionTreated:trial -7.568e-03 2.658e-03 3.036e+03
## intervention_proportion:conditionTreated:trial -2.443e-03 5.732e-03 3.036e+03
## t value Pr(>|t|)
## (Intercept) 19.193 < 2e-16 ***
## intervention_proportion 4.357 4.92e-05 ***
## conditionTreated 5.057 6.21e-07 ***
## trial 4.879 1.12e-06 ***
## intervention_proportion:conditionTreated -1.065 0.28733
## intervention_proportion:trial 2.120 0.03409 *
## conditionTreated:trial -2.847 0.00444 **
## intervention_proportion:conditionTreated:trial -0.426 0.67002
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrv_ cndtnT trial int_:T intr_: cndtT:
## intrvntn_pr -0.770
## conditnTrtd -0.189 0.170
## trial -0.440 0.357 0.188
## intrvntn_:T 0.204 -0.308 -0.926 -0.203
## intrvntn_p: 0.323 -0.487 -0.138 -0.733 0.277
## cndtnTrtd:t 0.124 -0.100 -0.669 -0.281 0.607 0.206
## intrvnt_:T: -0.135 0.204 0.613 0.307 -0.662 -0.418 -0.9174.4.2 按条件和比例的语义相似度
# 汇总统计
semantic_summary <- df %>%
group_by(proportion_factor, condition) %>%
summarise(
mean_similarity = mean(semantic_similarity, na.rm = TRUE),
se_similarity = sd(semantic_similarity, na.rm = TRUE) / sqrt(n()),
n = n(),
.groups = "drop"
)
kable(semantic_summary,
digits = 3,
caption = "Semantic Similarity to Reappraisal by Condition and Proportion") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| proportion_factor | condition | mean_similarity | se_similarity | n |
|---|---|---|---|---|
| 0% | Non-treated | 0.349 | 0.005 | 800 |
| 20% | Non-treated | 0.476 | 0.007 | 640 |
| 20% | Treated | 0.572 | 0.011 | 160 |
| 40% | Non-treated | 0.514 | 0.008 | 480 |
| 40% | Treated | 0.590 | 0.008 | 320 |
| 60% | Non-treated | 0.514 | 0.009 | 320 |
| 60% | Treated | 0.586 | 0.007 | 480 |
# 绘图
ggplot(semantic_summary, aes(x = proportion_factor, y = mean_similarity,
color = condition, group = condition)) +
geom_line(linewidth = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = mean_similarity - se_similarity,
ymax = mean_similarity + se_similarity),
width = 0.1) +
scale_color_manual(values = c("Non-treated" = "#2E86AB", "Treated" = "#E94F37")) +
labs(
title = "Semantic Similarity to Pure Reappraisal by Intervention Proportion",
subtitle = "Evidence for spread of reappraisal language",
x = "Proportion of Treated Participants",
y = "Similarity to Reappraisal",
color = "Condition"
) +
theme(legend.position = "bottom")
4.4.3 跨试次的语义相似度
# 计算按试次、条件和比例的均值
trial_summary <- df %>%
group_by(trial, proportion_factor, condition) %>%
summarise(
mean_similarity = mean(semantic_similarity, na.rm = TRUE),
se_similarity = sd(semantic_similarity, na.rm = TRUE) / sqrt(n()),
.groups = "drop"
)
# 按比例绘制跨试次的相似度
ggplot(trial_summary, aes(x = trial, y = mean_similarity,
color = condition, group = condition)) +
geom_line(linewidth = 0.8) +
geom_point(size = 1.5) +
geom_ribbon(aes(ymin = mean_similarity - se_similarity,
ymax = mean_similarity + se_similarity,
fill = condition), alpha = 0.2, color = NA) +
facet_wrap(~proportion_factor, ncol = 4) +
scale_color_manual(values = c("Non-treated" = "#2E86AB", "Treated" = "#E94F37")) +
scale_fill_manual(values = c("Non-treated" = "#2E86AB", "Treated" = "#E94F37")) +
labs(
title = "Semantic Similarity to Reappraisal Over Trials",
subtitle = "By intervention proportion condition (replicating Figure 4B from original paper)",
x = "Trial",
y = "Similarity to Reappraisal",
color = "Condition",
fill = "Condition"
) +
theme(legend.position = "bottom")
4.5 分析5：亲密度效应
# 包含亲密度交互作用的模型
model_closeness <- lmer(
emotion_rating ~ intervention_proportion * condition * closeness +
(1 | group_id) + (1 | participant_id) + (1 | trial),
data = df,
REML = TRUE
)
summary(model_closeness)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: emotion_rating ~ intervention_proportion * condition * closeness +
## (1 | group_id) + (1 | participant_id) + (1 | trial)
## Data: df
##
## REML criterion at convergence: 10311.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9495 -0.6600 -0.0149 0.6476 3.8663
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_id (Intercept) 0.4942 0.7030
## group_id (Intercept) 0.3927 0.6266
## trial (Intercept) 0.1040 0.3226
## Residual 1.2820 1.1323
## Number of obs: 3200, groups: participant_id, 160; group_id, 32; trial, 20
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 2.2334 0.3100
## intervention_proportion 3.2796 0.8510
## conditionTreated 1.6653 0.5830
## closenessHigh 0.8194 0.4263
## intervention_proportion:conditionTreated -1.4957 1.2808
## intervention_proportion:closenessHigh 0.5364 1.2035
## conditionTreated:closenessHigh -2.2531 0.8245
## intervention_proportion:conditionTreated:closenessHigh 2.6654 1.8113
## df t value
## (Intercept) 32.2525 7.205
## intervention_proportion 35.6971 3.854
## conditionTreated 140.3693 2.856
## closenessHigh 29.0220 1.922
## intervention_proportion:conditionTreated 142.3738 -1.168
## intervention_proportion:closenessHigh 35.6971 0.446
## conditionTreated:closenessHigh 140.3693 -2.733
## intervention_proportion:conditionTreated:closenessHigh 142.3738 1.472
## Pr(>|t|)
## (Intercept) 3.34e-08 ***
## intervention_proportion 0.000465 ***
## conditionTreated 0.004935 **
## closenessHigh 0.064485 .
## intervention_proportion:conditionTreated 0.244839
## intervention_proportion:closenessHigh 0.658533
## conditionTreated:closenessHigh 0.007090 **
## intervention_proportion:conditionTreated:closenessHigh 0.143343
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrv_ cndtnT clsnsH int_:T int_:H cndT:H
## intrvntn_pr -0.760
## conditnTrtd -0.154 0.158
## closenssHgh -0.688 0.553 0.112
## intrvntn_:T 0.165 -0.263 -0.934 -0.120
## intrvntn_:H 0.538 -0.707 -0.112 -0.782 0.186
## cndtnTrtd:H 0.109 -0.112 -0.707 -0.159 0.660 0.158
## intrvn_:T:H -0.117 0.186 0.660 0.170 -0.707 -0.263 -0.934# 按亲密度汇总
closeness_summary <- df %>%
group_by(closeness, condition, proportion_factor) %>%
summarise(
mean_emotion = mean(emotion_rating),
se_emotion = sd(emotion_rating) / sqrt(n()),
.groups = "drop"
)
# 绘图
ggplot(closeness_summary, aes(x = proportion_factor, y = mean_emotion,
color = condition, group = condition)) +
geom_line(linewidth = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = mean_emotion - se_emotion,
ymax = mean_emotion + se_emotion),
width = 0.1) +
facet_wrap(~closeness, labeller = labeller(closeness = c("Low" = "Low Closeness",
"High" = "High Closeness"))) +
scale_color_manual(values = c("Non-treated" = "#2E86AB", "Treated" = "#E94F37")) +
labs(
title = "Emotion Ratings by Closeness Condition",
x = "Proportion of Treated Participants",
y = "Negative Emotion Rating",
color = "Condition"
) +
theme(legend.position = "bottom")
4.6 分析6：稳健混合效应模型
鉴于可能存在正态性违背，我们还拟合了稳健混合效应模型。
# 使用robustlmm的稳健混合模型
model_robust <- rlmer(
emotion_rating ~ condition * intervention_proportion +
(1 | group_id) + (1 | participant_id),
data = df
)
summary(model_robust)## Robust linear mixed model fit by DAStau
## Formula: emotion_rating ~ condition * intervention_proportion + (1 | group_id) + (1 | participant_id)
## Data: df
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1111 -0.6662 0.0122 0.6565 3.8639
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant_id (Intercept) 0.4995 0.7068
## group_id (Intercept) 0.5068 0.7119
## Residual 1.3529 1.1631
## Number of obs: 3200, groups: participant_id, 160; group_id, 32
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.5407 0.2419 10.503
## conditionTreated 0.5279 0.4289 1.231
## intervention_proportion 3.7304 0.6770 5.510
## conditionTreated:intervention_proportion -0.1088 0.9431 -0.115
##
## Correlation of Fixed Effects:
## (Intr) cndtnT intrv_
## conditnTrtd -0.146
## intrvntn_pr -0.785 0.148
## cndtnTrtd:_ 0.156 -0.935 -0.244
##
## Robustness weights for the residuals:
## 2549 weights are ~= 1. The remaining 651 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.348 0.687 0.813 0.796 0.942 0.999
##
## Robustness weights for the random effects:
## 157 weights are ~= 1. The remaining 35 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.389 0.680 0.859 0.812 0.949 0.999
##
## Rho functions used for fitting:
## Residuals:
## eff: smoothed Huber (k = 1.345, s = 10)
## sig: smoothed Huber, Proposal 2 (k = 1.345, s = 10)
## Random Effects, variance component 1 (participant_id):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal 2 (k = 1.345, s = 10)
## Random Effects, variance component 2 (group_id):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal 2 (k = 1.345, s = 10)##
## === Comparison: Standard vs. Robust Estimates ===comparison_df <- data.frame(
Parameter = c("(Intercept)", "Condition (Treated)",
"Intervention Proportion", "Interaction"),
Standard = fixef(model_proportion),
Robust = fixef(model_robust)
)
kable(comparison_df, digits = 3) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Parameter | Standard | Robust | |
|---|---|---|---|
| (Intercept) | (Intercept) | 2.643 | 2.541 |
| conditionTreated | Condition (Treated) | 0.535 | 0.528 |
| intervention_proportion | Intervention Proportion | 3.546 | 3.730 |
| conditionTreated:intervention_proportion | Interaction | -0.154 | -0.109 |
5 汇总图（复现原论文图3）
# 面板A：主要结果（剂量-反应）
panel_a <- ggplot(agg_data, aes(x = intervention_proportion, y = mean_emotion,
color = condition)) +
geom_point(size = 4, alpha = 0.9) +
geom_errorbar(aes(ymin = mean_emotion - se_emotion,
ymax = mean_emotion + se_emotion),
width = 0.015, linewidth = 0.8) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2),
se = TRUE, alpha = 0.15, linewidth = 1.2) +
scale_color_manual(values = c("Non-treated" = "#2E86AB", "Treated" = "#E94F37")) +
scale_x_continuous(labels = scales::percent_format(accuracy = 1),
breaks = c(0, 0.2, 0.4, 0.6)) +
labs(
title = "A | Effect of Intervention Proportion on Negative Emotions",
x = "Proportion of Treated Individuals",
y = "Negative Emotion",
color = "Condition"
) +
theme(
legend.position = c(0.85, 0.85),
legend.background = element_rect(fill = "white", color = NA)
)
# 面板B：语义相似度
panel_b <- ggplot(semantic_summary, aes(x = proportion_factor, y = mean_similarity,
color = condition, group = condition)) +
geom_line(linewidth = 1.2) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = mean_similarity - se_similarity,
ymax = mean_similarity + se_similarity),
width = 0.15, linewidth = 0.8) +
scale_color_manual(values = c("Non-treated" = "#2E86AB", "Treated" = "#E94F37")) +
labs(
title = "B | Similarity to Reappraisal Language",
x = "Proportion of Treated Individuals",
y = "Similarity to Reappraisal",
color = "Condition"
) +
theme(
legend.position = c(0.85, 0.25),
legend.background = element_rect(fill = "white", color = NA)
)
# 组合面板
panel_a / panel_b +
plot_annotation(
title = "Results from the Main Study",
subtitle = "Replication of Pinus et al. (2025) Figure 3",
theme = theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray40")
)
)
6 结果摘要
## ============================================================## RESULTS SUMMARY: Emotion Regulation Contagion Replication## ============================================================## 1. BASIC TREATMENT EFFECT## ----------------------------------------basic_coef <- fixef(model_basic)
cat(sprintf(" Treated vs. Non-treated: beta = %.3f\n", basic_coef["conditionTreated"]))## Treated vs. Non-treated: beta = 0.584basic_ci <- confint(model_basic, parm = "conditionTreated", method = "Wald")
cat(sprintf(" 95%% CI: [%.3f, %.3f]\n\n", basic_ci[1], basic_ci[2]))## 95% CI: [0.283, 0.885]## 2. DOSE-RESPONSE RELATIONSHIP## ----------------------------------------## Best fitting model: Linearcat(sprintf(" AIC difference from best: Linear=%.2f, Quadratic=%.2f\n\n",
model_comparison$delta_AIC[model_comparison$Model == "Linear"],
model_comparison$delta_AIC[model_comparison$Model == "Quadratic"]))## AIC difference from best: Linear=0.00, Quadratic=0.21## 3. EFFECT ON NON-TREATED PARTICIPANTS (Contagion Effect)## ----------------------------------------nontreated_coef <- fixef(model_nontreated_quad)
cat(sprintf(" Linear effect of proportion: beta = %.3f\n", nontreated_coef[2]))## Linear effect of proportion: beta = 6.322## Quadratic effect: beta = -4.753## 4. SEMANTIC SIMILARITY (Evidence for Social Learning)## ----------------------------------------semantic_coef <- fixef(model_semantic)
cat(sprintf(" Effect of proportion on similarity: beta = %.3f\n",
semantic_coef["intervention_proportion"]))## Effect of proportion on similarity: beta = 0.225cat(sprintf(" Condition x Proportion interaction: beta = %.3f\n\n",
semantic_coef["intervention_proportion:conditionTreated"]))## Condition x Proportion interaction: beta = -0.097## ============================================================## CONCLUSION: Results provide evidence for emotion regulation## contagion, showing that non-treated participants show reduced## negative emotions and increased reappraisal language use when## the proportion of treated participants in their group increases.## ============================================================7 会话信息
## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS 26.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Asia/Shanghai
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] robustlmm_3.3-3 nortest_1.0-4 car_3.1-3 carData_3.0-5
## [5] kableExtra_1.4.0 knitr_1.50 scales_1.3.0 patchwork_1.3.2
## [9] ggeffects_2.3.2 emmeans_1.11.2-8 effectsize_1.0.1 performance_0.14.0
## [13] lmerTest_3.2-0 lme4_1.1-37 Matrix_1.7-1 lubridate_1.9.4
## [17] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.4
## [21] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.2
## [25] tidyverse_2.0.0 readxl_1.4.5
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.2 farver_2.1.2
## [4] fastmap_1.2.0 TH.data_1.1-4 bayestestR_0.16.0
## [7] digest_0.6.37 estimability_1.5.1 timechange_0.3.0
## [10] lifecycle_1.0.4 fastGHQuad_1.0.1 survival_3.7-0
## [13] magrittr_2.0.3 compiler_4.4.2 rlang_1.1.5
## [16] sass_0.4.10 tools_4.4.2 yaml_2.3.10
## [19] labeling_0.4.3 xml2_1.3.8 abind_1.4-8
## [22] multcomp_1.4-28 withr_3.0.2 numDeriv_2016.8-1.1
## [25] grid_4.4.2 datawizard_1.3.0 xtable_1.8-4
## [28] colorspace_2.1-1 MASS_7.3-61 insight_1.4.4
## [31] cli_3.6.4 mvtnorm_1.3-3 rmarkdown_2.29
## [34] reformulas_0.4.0 generics_0.1.3 robustbase_0.99-4-1
## [37] rstudioapi_0.17.1 tzdb_0.5.0 parameters_0.26.0
## [40] minqa_1.2.8 cachem_1.1.0 splines_4.4.2
## [43] cellranger_1.1.0 vctrs_0.6.5 boot_1.3-31
## [46] sandwich_3.1-1 jsonlite_2.0.0 hms_1.1.3
## [49] Formula_1.2-5 systemfonts_1.2.3 jquerylib_0.1.4
## [52] glue_1.8.0 DEoptimR_1.1-3-1 nloptr_2.2.1
## [55] codetools_0.2-20 stringi_1.8.7 gtable_0.3.6
## [58] munsell_0.5.1 pillar_1.10.1 htmltools_0.5.8.1
## [61] R6_2.6.1 textshaping_1.0.1 Rdpack_2.6.4
## [64] evaluate_1.0.3 lattice_0.22-6 rbibutils_2.3
## [67] bslib_0.9.0 Rcpp_1.0.14 svglite_2.2.1
## [70] nlme_3.1-166 mgcv_1.9-1 xfun_0.52
## [73] zoo_1.8-14 pkgconfig_2.0.3