library(tidyverse)
library(psych)
library(mediation)
library(knitr)
library(kableExtra)
setwd("/Users/anjalisingh/Downloads/")

raw <- read_csv(
  "#294906_+'2(form;+assigned)+x+2(severity;+assigned)'_June+5,+2026_08.32.csv",
  skip = 0
)

df <- raw %>%
  filter(!is.na(participantId) & participantId != "") %>%
  filter(Deb == 1, E_AC == -3, R_AC == 3) %>%
  mutate(across(c(E1, E2, E3, R1, R2, R3), as.numeric)) %>%
  mutate(
    E_AVG  = (E1 + E2 + E3) / 3,
    R_AVG  = (R1 + R2 + R3) / 3,
    Cond2_d = ifelse(Cond2 == "Pill", 1, 0),
    Cond1_d = ifelse(Cond1 == "S",    1, 0),
    MC_D   = as.numeric(MC_D),
    MC_F   = as.numeric(MC_F),
    Sev    = as.numeric(Sev),
    FP1    = as.numeric(FP1),
    FP2    = as.numeric(FP2),
    FP_AVG = (FP1 + FP2) / 2
  )

Sample & Scale Reliability

Sample Size

cat("Total rows after cleaning:", nrow(df))
## Total rows after cleaning: 483
df %>%
  count(Cond1, Cond2) %>%
  rename(Severity = Cond1, Format = Cond2, N = n) %>%
  kable(caption = "Cell sizes after attention check exclusions") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Cell sizes after attention check exclusions
Severity Format N
NS Gum 125
NS Pill 114
S Gum 112
S Pill 132

Cronbach’s Alpha

alpha_E <- psych::alpha(df[, c("E1", "E2", "E3")])
alpha_R <- psych::alpha(df[, c("R1", "R2", "R3")])

tibble(
  Scale    = c("Efficacy (E1-E3)", "Risk (R1-R3)"),
  `raw alpha`    = c(round(alpha_E$total$raw_alpha, 3),
                     round(alpha_R$total$raw_alpha, 3)),
  `std alpha`    = c(round(alpha_E$total$std.alpha, 3),
                     round(alpha_R$total$std.alpha, 3)),
  `avg r`        = c(round(alpha_E$total$average_r, 3),
                     round(alpha_R$total$average_r, 3))
) %>%
  kable(caption = "Scale reliability") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Scale reliability
Scale raw alpha std alpha avg r
Efficacy (E1-E3) 0.974 0.974 0.925
Risk (R1-R3) 0.936 0.937 0.831

Manipulation Checks

Format MC (Drug-like)

t_d <- t.test(MC_D ~ Cond2, data = df)
cat(sprintf("MC_D: t(%.1f) = %.3f, p < .001\n", t_d$parameter, t_d$statistic))
## MC_D: t(454.3) = -17.675, p < .001
plot_mc_d <- df %>%
  group_by(Cond2) %>%
  summarise(M = mean(MC_D, na.rm = TRUE),
            SE = sd(MC_D, na.rm = TRUE) / sqrt(n()), .groups = "drop")

ggplot(plot_mc_d, aes(x = Cond2, y = M, fill = Cond2)) +
  geom_bar(stat = "identity", width = 0.5) +
  geom_errorbar(aes(ymin = M - SE, ymax = M + SE), width = 0.15) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  scale_fill_manual(values = c("Gum" = "#E69F00", "Pill" = "#0072B2")) +
  scale_x_discrete(labels = c("Gum" = "Gummy", "Pill" = "Pill")) +
  labs(title = "Drug-like Perception by Format",
       subtitle = paste0("t(", round(t_d$parameter, 1), ") = ",
                         round(t_d$statistic, 2), ", p < .001"),
       x = "Format", y = "Mean Drug-like Rating (MC_D)") +
  theme_classic(base_size = 13) +
  theme(legend.position = "none")

Format MC (Food-like)

t_f <- t.test(MC_F ~ Cond2, data = df)
cat(sprintf("MC_F: t(%.1f) = %.3f, p < .001\n", t_f$parameter, t_f$statistic))
## MC_F: t(378.5) = 27.965, p < .001
plot_mc_f <- df %>%
  group_by(Cond2) %>%
  summarise(M = mean(MC_F, na.rm = TRUE),
            SE = sd(MC_F, na.rm = TRUE) / sqrt(n()), .groups = "drop")

ggplot(plot_mc_f, aes(x = Cond2, y = M, fill = Cond2)) +
  geom_bar(stat = "identity", width = 0.5) +
  geom_errorbar(aes(ymin = M - SE, ymax = M + SE), width = 0.15) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  scale_fill_manual(values = c("Gum" = "#E69F00", "Pill" = "#0072B2")) +
  scale_x_discrete(labels = c("Gum" = "Gummy", "Pill" = "Pill")) +
  labs(title = "Food-like Perception by Format",
       subtitle = paste0("t(", round(t_f$parameter, 1), ") = ",
                         round(t_f$statistic, 2), ", p < .001"),
       x = "Format", y = "Mean Food-like Rating (MC_F)") +
  theme_classic(base_size = 13) +
  theme(legend.position = "none")

Severity MC

t_s <- t.test(Sev ~ Cond1, data = df)
cat(sprintf("Sev: t(%.1f) = %.3f, p < .001\n", t_s$parameter, t_s$statistic))
## Sev: t(468.3) = -27.781, p < .001
plot_sev <- df %>%
  group_by(Cond1) %>%
  summarise(M = mean(Sev, na.rm = TRUE),
            SE = sd(Sev, na.rm = TRUE) / sqrt(n()), .groups = "drop")

ggplot(plot_sev, aes(x = Cond1, y = M, fill = Cond1)) +
  geom_bar(stat = "identity", width = 0.5) +
  geom_errorbar(aes(ymin = M - SE, ymax = M + SE), width = 0.15) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  scale_fill_manual(values = c("NS" = "#CC79A7", "S" = "#009E73")) +
  scale_x_discrete(labels = c("NS" = "Non-Severe", "S" = "Severe")) +
  labs(title = "Perceived Severity by Condition",
       subtitle = paste0("t(", round(t_s$parameter, 1), ") = ",
                         round(t_s$statistic, 2), ", p < .001"),
       x = "Severity Condition", y = "Mean Perceived Severity (Sev)") +
  theme_classic(base_size = 13) +
  theme(legend.position = "none")


Main Effects: DVs

Risk Perception (R_AVG)

mod_R2 <- lm(R_AVG ~ Cond1_d * Cond2_d, data = df)
summary(mod_R2)
## 
## Call:
## lm(formula = R_AVG ~ Cond1_d * Cond2_d, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2753 -0.9419 -0.1607  1.0581  4.0853 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.0853     0.1299  -8.355  7.1e-16 ***
## Cond1_d           0.5794     0.1890   3.066  0.00229 ** 
## Cond2_d           0.6146     0.1881   3.268  0.00116 ** 
## Cond1_d:Cond2_d   0.1666     0.2649   0.629  0.52966    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.452 on 479 degrees of freedom
## Multiple R-squared:  0.1063, Adjusted R-squared:  0.1007 
## F-statistic: 18.99 on 3 and 479 DF,  p-value: 1.186e-11
plot_df <- df %>%
  group_by(Cond1, Cond2) %>%
  summarise(M = mean(R_AVG, na.rm = TRUE),
            SE = sd(R_AVG, na.rm = TRUE) / sqrt(n()), .groups = "drop")

ggplot(plot_df, aes(x = Cond1, y = M, color = Cond2, group = Cond2)) +
  geom_point(size = 3) +
  geom_line(linewidth = 1) +
  geom_errorbar(aes(ymin = M - SE, ymax = M + SE), width = 0.1) +
  scale_color_manual(values = c("Gum" = "#E69F00", "Pill" = "#0072B2"),
                     labels = c("Gum" = "Gummy", "Pill" = "Pill")) +
  scale_x_discrete(labels = c("NS" = "Non-Severe", "S" = "Severe")) +
  labs(title = "Effect of Severity and Format on Risk Perception",
       x = "Condition Severity", y = "Mean Risk Perception (R_AVG)",
       color = "Format") +
  theme_classic(base_size = 13) +
  theme(legend.position = "top")

Efficacy Perception (E_AVG)

mod_E <- lm(E_AVG ~ Cond1_d * Cond2_d, data = df)
summary(mod_E)
## 
## Call:
## lm(formula = E_AVG ~ Cond1_d * Cond2_d, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9524 -0.7879  0.1228  1.0476  2.3653 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.6347     0.1130   5.615 3.33e-08 ***
## Cond1_d           0.3177     0.1644   1.932   0.0539 .  
## Cond2_d           0.2425     0.1636   1.482   0.1390    
## Cond1_d:Cond2_d  -0.4070     0.2305  -1.766   0.0781 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.264 on 479 degrees of freedom
## Multiple R-squared:  0.00867,    Adjusted R-squared:  0.002461 
## F-statistic: 1.396 on 3 and 479 DF,  p-value: 0.2431
plot_df_E <- df %>%
  group_by(Cond1, Cond2) %>%
  summarise(M = mean(E_AVG, na.rm = TRUE),
            SE = sd(E_AVG, na.rm = TRUE) / sqrt(n()), .groups = "drop")

ggplot(plot_df_E, aes(x = Cond1, y = M, color = Cond2, group = Cond2)) +
  geom_point(size = 3) +
  geom_line(linewidth = 1) +
  geom_errorbar(aes(ymin = M - SE, ymax = M + SE), width = 0.1) +
  scale_color_manual(values = c("Gum" = "#E69F00", "Pill" = "#0072B2"),
                     labels = c("Gum" = "Gummy", "Pill" = "Pill")) +
  scale_x_discrete(labels = c("NS" = "Non-Severe", "S" = "Severe")) +
  labs(title = "Effect of Severity and Format on Efficacy Perception",
       x = "Condition Severity", y = "Mean Efficacy Perception (E_AVG)",
       color = "Format") +
  theme_classic(base_size = 13) +
  theme(legend.position = "top")

Fair Price (FP_AVG)

mod_FP <- lm(FP_AVG ~ Cond1_d * Cond2_d, data = df)
summary(mod_FP)
## 
## Call:
## lm(formula = FP_AVG ~ Cond1_d * Cond2_d, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1491 -1.1491  0.0114  1.3897  3.4520 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      -0.4520     0.1531  -2.952  0.00331 **
## Cond1_d           0.5234     0.2227   2.350  0.01918 * 
## Cond2_d           0.6011     0.2217   2.711  0.00694 **
## Cond1_d:Cond2_d  -0.6839     0.3123  -2.190  0.02901 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.712 on 479 degrees of freedom
## Multiple R-squared:  0.01846,    Adjusted R-squared:  0.01232 
## F-statistic: 3.003 on 3 and 479 DF,  p-value: 0.03016
plot_df_FP <- df %>%
  group_by(Cond1, Cond2) %>%
  summarise(M = mean(FP_AVG, na.rm = TRUE),
            SE = sd(FP_AVG, na.rm = TRUE) / sqrt(n()), .groups = "drop")

ggplot(plot_df_FP, aes(x = Cond1, y = M, color = Cond2, group = Cond2)) +
  geom_point(size = 3) +
  geom_line(linewidth = 1) +
  geom_errorbar(aes(ymin = M - SE, ymax = M + SE), width = 0.1) +
  scale_color_manual(values = c("Gum" = "#E69F00", "Pill" = "#0072B2"),
                     labels = c("Gum" = "Gummy", "Pill" = "Pill")) +
  scale_x_discrete(labels = c("NS" = "Non-Severe", "S" = "Severe")) +
  labs(title = "Effect of Severity and Format on Fair Price",
       x = "Condition Severity", y = "Mean Fair Price (FP_AVG)",
       color = "Format") +
  theme_classic(base_size = 13) +
  theme(legend.position = "top")


Mediation: Risk Perception

MC_D -> R_AVG

med_mc  <- lm(MC_D  ~ Cond2_d,        data = df)
out_mc  <- lm(R_AVG ~ Cond2_d + MC_D, data = df)

med_mc_result <- mediate(med_mc, out_mc,
                         treat = "Cond2_d", mediator = "MC_D",
                         robustSE = TRUE, sims = 1000)
summary(med_mc_result)
## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME           0.337484     0.096683     0.597080   0.006 ** 
## ADE            0.397329     0.038330     0.756321   0.030 *  
## Total Effect   0.734812     0.470606     0.997786  <2e-16 ***
## Prop. Mediated 0.455018     0.136482     0.930878   0.006 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 483 
## 
## 
## Simulations: 1000

MC_F -> R_AVG

med_mcf  <- lm(MC_F  ~ Cond2_d,        data = df)
out_mcf  <- lm(R_AVG ~ Cond2_d + MC_F, data = df)

med_mcf_result <- mediate(med_mcf, out_mcf,
                          treat = "Cond2_d", mediator = "MC_F",
                          robustSE = TRUE, sims = 1000)
summary(med_mcf_result)
## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                 Estimate 95% CI Lower 95% CI Upper p-value    
## ACME            0.424126     0.055982     0.764542   0.030 *  
## ADE             0.309939    -0.135331     0.765953   0.148    
## Total Effect    0.734065     0.460630     0.989548  <2e-16 ***
## Prop. Mediated  0.579140     0.075728     1.244433   0.030 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 483 
## 
## 
## Simulations: 1000

Mediation: Fair Price

MC_D -> FP_AVG

med_mc_fp  <- lm(MC_D   ~ Cond2_d,        data = df)
out_mc_fp  <- lm(FP_AVG ~ Cond2_d + MC_D, data = df)

med_mc_fp_result <- mediate(med_mc_fp, out_mc_fp,
                            treat = "Cond2_d", mediator = "MC_D",
                            robustSE = TRUE, sims = 1000)
summary(med_mc_fp_result)
## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                 Estimate 95% CI Lower 95% CI Upper p-value    
## ACME            0.516891     0.236344     0.782235  <2e-16 ***
## ADE            -0.252556    -0.638131     0.151543   0.226    
## Total Effect    0.264335    -0.037788     0.558633   0.082 .  
## Prop. Mediated  1.843023    -7.920690    14.677435   0.082 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 483 
## 
## 
## Simulations: 1000

MC_F -> FP_AVG

med_mcf_fp  <- lm(MC_F   ~ Cond2_d,        data = df)
out_mcf_fp  <- lm(FP_AVG ~ Cond2_d + MC_F, data = df)

med_mcf_fp_result <- mediate(med_mcf_fp, out_mcf_fp,
                             treat = "Cond2_d", mediator = "MC_F",
                             robustSE = TRUE, sims = 1000)
summary(med_mcf_fp_result)
## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                 Estimate 95% CI Lower 95% CI Upper p-value  
## ACME           -0.177419    -0.596004     0.234458   0.430  
## ADE             0.446678    -0.049502     0.962418   0.080 .
## Total Effect    0.269259    -0.039860     0.580711   0.084 .
## Prop. Mediated -0.584328    -7.997847     3.259502   0.490  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 483 
## 
## 
## Simulations: 1000

R_AVG -> FP_AVG

med_risk    <- lm(R_AVG  ~ Cond2_d,           data = df)
out_risk_fp <- lm(FP_AVG ~ Cond2_d + R_AVG,   data = df)

med_risk_fp_result <- mediate(med_risk, out_risk_fp,
                              treat = "Cond2_d", mediator = "R_AVG",
                              robustSE = TRUE, sims = 1000)
summary(med_risk_fp_result)
## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                 Estimate 95% CI Lower 95% CI Upper p-value  
## ACME           -0.098886    -0.192973    -0.012440   0.020 *
## ADE             0.371518     0.054123     0.702435   0.018 *
## Total Effect    0.272632    -0.027694     0.600062   0.086 .
## Prop. Mediated -0.330158    -3.048247     1.822511   0.106  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 483 
## 
## 
## Simulations: 1000

Report generated with R Markdown. Code visible via “Code” buttons.