4 Main Analyses
4.1 Analysis 1: Effect of Reappraisal on Emotion Ratings
This analysis tests whether the reappraisal intervention reduces negative emotions, replicating the basic treatment effect from the original study.
4.1.1 Normality Check
# Kolmogorov-Smirnov test for normality
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# Visual check
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 Linear Mixed Model: Basic Treatment Effect
# Model 1: Basic effect of condition on emotion ratings
# Following the original paper's approach with random intercepts
model_basic <- lmer(
emotion_rating ~ condition +
(1 | group_id) +
(1 | participant_id) +
(1 | trial),
data = df,
REML = TRUE
)
# Model summary
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 Model with Intervention Proportion
# Model 2: Effect of intervention proportion on emotion ratings
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 Analysis 2: Dose-Response Relationship
This analysis examines the relationship between the proportion of treated participants and emotion reduction, testing different functional forms (linear, quadratic, cubic, exponential).
4.2.1 Model Comparison: Linear vs. Quadratic vs. Cubic vs. Exponential
# Prepare data for model comparison
# Center proportion for polynomial models
df <- df %>%
mutate(
prop_centered = scale(intervention_proportion, scale = FALSE),
prop_scaled = scale(intervention_proportion)
)
# Linear model
model_linear <- lmer(
emotion_rating ~ intervention_proportion + condition +
(1 | group_id) + (1 | participant_id) + (1 | trial),
data = df, REML = FALSE
)
# Quadratic model
model_quadratic <- lmer(
emotion_rating ~ poly(intervention_proportion, 2, raw = TRUE) + condition +
(1 | group_id) + (1 | participant_id) + (1 | trial),
data = df, REML = FALSE
)
# Cubic model
model_cubic <- lmer(
emotion_rating ~ poly(intervention_proportion, 3, raw = TRUE) + condition +
(1 | group_id) + (1 | participant_id) + (1 | trial),
data = df, REML = FALSE
)
# Exponential transformation
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
)
# Model comparison using 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 Best Fitting Model Summary
# Select best model based on AIC (typically exponential as in the original paper)
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# Use the quadratic model (commonly the best fit, similar to original study)
cat("=== Quadratic Model Summary ===\n")## === 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 Visualize Dose-Response Relationship
# Get predictions from each model
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
# Predictions (using fixed effects only for visualization)
# We'll use aggregate data for cleaner visualization
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"
)
# Plot with smoothed lines
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 Analysis 3: Separate Analysis for Treated and Non-Treated Participants
Following the original paper’s approach, we analyze treated and non-treated participants separately.
4.3.1 Non-Treated Participants Only
# Filter non-treated participants
df_nontreated <- df %>% filter(condition == "Non-treated")
cat("=== Analysis of Non-Treated Participants ===\n")## === Analysis of Non-Treated Participants ===## N observations: 2240# Model for non-treated participants
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# Effect of proportion on non-treated participants' emotions
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 Treated Participants Only
# Filter treated participants
df_treated <- df %>% filter(condition == "Treated")
cat("=== Analysis of Treated Participants ===\n")## === Analysis of Treated Participants ===## N observations: 960# Model for treated participants
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 Visualize Separate Effects
# Create separate visualizations
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"
)
# Plot for non-treated
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"
)
# Plot for treated
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"
)
# Combine plots
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 Analysis 4: Semantic Similarity Analysis (Replicating Semantic Projection Analysis)
This analysis tests whether non-treated participants adopt reappraisal language from treated participants, providing evidence for the social learning mechanism.
4.4.1 Model: Effect of Proportion on Semantic Similarity
# Full model with interaction
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 Similarity by Condition and Proportion
# Summary statistics
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 |
# Plot
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 Semantic Similarity Over Trials
# Calculate means by trial, condition, and proportion
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"
)
# Plot similarity over trials by proportion
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 Analysis 5: Effect of Closeness
# Model with closeness interaction
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# Summary by closeness
closeness_summary <- df %>%
group_by(closeness, condition, proportion_factor) %>%
summarise(
mean_emotion = mean(emotion_rating),
se_emotion = sd(emotion_rating) / sqrt(n()),
.groups = "drop"
)
# Plot
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 Analysis 6: Robust Mixed Effects Models
Given potential violations of normality, we also fit robust mixed-effects models.
# Robust mixed model using 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 |
