This document reproduces all analyses reported in the paper across three studies (N = 549 total). Data are hosted on OSF at https://osf.io/tv7xb/.
library(dplyr)
library(ggplot2)
library(patchwork)
library(emmeans)
library(psych)
# ── Response-scale maps ──────────────────────────────────────────────────────
agree_map <- c(
"Strognly disagree" = 1, # typo present in raw data
"Strongly disagree" = 1,
"Somewhat disagree" = 2,
"Neither agree nor disagree" = 3,
"Somewhat agree" = 4,
"Strongly agree" = 5
)
likely_map <- c(
"Extremely unlikely" = 1,
"Somewhat unlikely" = 2,
"Neither likely nor unlikely" = 3,
"Somewhat likely" = 4,
"Extremely likely" = 5
)
wtp_map <- c(
"Nothing more" = 0,
"$0.50 more" = 0.5,
"$1.00 more" = 1,
"$1.50 more" = 1.5,
"$2.00 more" = 2,
"$2.50 more" = 2.5,
"$3.00 more" = 3,
"$3.50 more" = 3.5,
"$4.00 more" = 4,
"$4.50 more" = 4.5,
"$5.00+ more" = 5
)
yes_no_map <- c(
"Definitely no" = 1,
"Probably no" = 2,
"Not sure" = 3,
"Probably yes" = 4,
"Definitely yes" = 5
)
# ── Utility: replace "-99" with NA ──────────────────────────────────────────
na_99 <- function(x) {
x <- as.character(x)
x[x == "-99"] <- NA
x
}
as_num_na99 <- function(x) as.numeric(na_99(x))
# ── Shared Likert histogram theme & plot function ────────────────────────────
theme_apa_hist <- function(base_size = 14) {
theme_classic(base_size = base_size) %+replace%
theme(
panel.grid = element_blank(),
axis.line = element_line(color = "black", linewidth = 0.8),
plot.title = element_text(hjust = 0.5, face = "plain",
size = base_size + 1,
margin = margin(b = 10)),
axis.title.y = element_text(angle = 90, margin = margin(r = 10)),
axis.title.x = element_text(margin = margin(t = 8)),
plot.margin = margin(6, 8, 6, 8)
)
}
likert_labels <- c(
"Strongly\ndisagree", "Somewhat\ndisagree",
"Neither agree\nnor disagree", "Somewhat\nagree", "Strongly\nagree"
)
plot_likert_hist <- function(df, var, title, show_x = FALSE) {
df2 <- df |>
transmute(x = .data[[var]]) |>
filter(!is.na(x), x >= 1, x <= 5)
ggplot(df2, aes(x = x)) +
geom_bar(width = 1, color = "black", fill = "grey75", linewidth = 0.5) +
scale_x_continuous(
breaks = 1:5,
labels = if (show_x) likert_labels else NULL,
expand = c(0, 0)
) +
scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
labs(
title = title,
x = if (show_x) "Response" else NULL,
y = "Count"
) +
theme_apa_hist()
}
Goal: Test how often consumers infer that a GLP-1 “support” smoothie powder literally contains GLP-1 as an ingredient.
# ── Replace these URLs with your OSF download links ─────────────────────────
# Pattern: https://osf.io/FILEID/download
# S1 file is stored under https://osf.io/tv7xb/files/t7cyn — locate the
# individual file permalink on OSF and append /download.
# Example placeholder (update before publishing):
s1_url <- "https://osf.io/t7cyn/download"
dat1 <- read.csv(s1_url, comment.char = "#")
# Drop extra header rows and test responses
dat1 <- dat1[-c(1:10), ]
# Keep only participants who passed the attention check
dat1 <- subset(dat1, product.type == "Smoothie Supplement Powder")
x <- dat1$contains.glp1
blank <- is.na(x) | x == ""
has_glp1 <- grepl("\\bGLP-1\\b", x)
count_selected <- ifelse(blank, 0, lengths(strsplit(x, ",")))
N <- length(x)
tab1 <- data.frame(
Group = c(
"Selected GLP-1 as an ingredient",
"Selected GLP-1 as the *only* ingredient",
"Selected GLP-1 with other ingredients",
"Selected other ingredients but not GLP-1"
),
n = c(
sum(has_glp1),
sum(has_glp1 & count_selected == 1),
sum(has_glp1 & count_selected > 1),
sum(!has_glp1 & count_selected > 0)
)
)
tab1$`%` <- round(100 * tab1$n / N, 1)
knitr::kable(
tab1,
caption = "Table 1. Ingredient selections for the GLP-1 support product (Study 1, N = 117)."
)
| Group | n | % |
|---|---|---|
| Selected GLP-1 as an ingredient | 114 | 97.4 |
| Selected GLP-1 as the only ingredient | 73 | 62.4 |
| Selected GLP-1 with other ingredients | 41 | 35.0 |
| Selected other ingredients but not GLP-1 | 3 | 2.6 |
Because endorsement of GLP-1 as an ingredient was near ceiling (97.4% selected GLP-1), there was insufficient variability to support meaningful inferential tests. Results are limited to descriptive statistics.
Goal: Replicate the ingredient-inference effect using graded dose-response belief measures and three advertisement stimuli (N = 136).
s2_url <- "https://osf.io/m3xtc/download" # update with OSF permalink
dat2 <- read.csv(s2_url, comment.char = "#")
# Drop extra header rows; rows 52–53 are early participants who missed items
dat2 <- dat2[-c(1:11, 52:53), ]
dat2 <- subset(dat2, product.type == "Smoothie Supplement Powder")
# ── Recode core variables ────────────────────────────────────────────────────
dat2$effect.like.glp1 <- unname(agree_map[as.character(dat2$effect.like.glp1)])
dat2$similar.to.glp1 <- unname(agree_map[as.character(dat2$similar.to.glp1)])
dat2$works.like.glp1 <- unname(agree_map[as.character(dat2$works.like.glp1)])
# False-belief composite
dat2$glp1_belief_mean <- rowMeans(
dat2[, c("effect.like.glp1", "similar.to.glp1", "works.like.glp1")],
na.rm = TRUE
)
# Product evaluation composite
eval_items2 <- data.frame(
purchase1 = unname(likely_map[as.character(dat2$purchase1)]),
purchase2 = unname(agree_map[as.character(dat2$purchase2)]),
appealing = unname(agree_map[as.character(dat2$appealing)])
)
dat2$eval_mean <- rowMeans(eval_items2, na.rm = TRUE)
# WTP
dat2$wtp <- as.numeric(unname(wtp_map[as.character(dat2$wtp)]))
# Prescription-knowledge item
script_map <- c(
"Definitely no" = 1, "Probably no" = 2, "Not sure" = 3,
"Probably yes" = 4, "Definitely yes" = 5
)
dat2$script.required <- unname(script_map[as.character(dat2$script.required)])
# Demographics
dat2$age <- as.numeric(dat2$age)
dat2$bmi <- (as.numeric(dat2$weight) / (as.numeric(dat2$height_1)^2)) * 703
dat2$income <- as.numeric(factor(
na_99(dat2$income),
levels = c("$0 - $19,999","$20,000 - $39,999","$40,000 - $59,999",
"$60,000 - $89,999","$90,000 - $119,999",
"$120,000 - $149,999","$150,000+"),
ordered = TRUE
))
dat2$ses <- as.numeric(sub("^([0-9]+).*", "\\1", as.character(dat2$ses)))
dat2$gender_simple <- case_when(
dat2$gender %in% c("Female","Female,Non-binary / third gender") ~ "Female",
dat2$gender %in% c("Male","Male,Non-binary / third gender","Male,Transgender") ~ "Male",
TRUE ~ NA_character_
)
race2 <- dat2$race
race2[grepl(",", race2) | grepl("prefer not to say", race2, ignore.case = TRUE)] <-
"Other / Prefer not to say"
dat2$race_simple <- case_when(
race2 == "White or European American" ~ "White",
race2 == "Black or African American" ~ "Black",
race2 == "Asian" ~ "Asian",
is.na(race2) ~ NA_character_,
TRUE ~ "Other"
) |> factor(levels = c("White","Black","Asian","Other"))
taking2 <- unname(yes_no_map[as.character(dat2$taking.glp1)])
dat2$glp1_exposure <- case_when(
taking2 %in% c(4,5) ~ 1,
taking2 %in% c(1,2) ~ 0,
TRUE ~ NA_real_
)
dat2$condition <- as.factor(dat2$condition)
# False-belief items
cat("=== False-belief scale (3 items) ===\n")
## === False-belief scale (3 items) ===
psych::alpha(dat2[, c("effect.like.glp1","similar.to.glp1","works.like.glp1")])[1]
## $total
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd
## 0.9425831 0.9427783 0.9187885 0.8459634 16.47589 0.00860451 3.676471 1.108843
## median_r
## 0.8320618
# Inter-item correlations
cat("\nInter-item correlations:\n")
##
## Inter-item correlations:
round(cor(dat2[, c("effect.like.glp1","similar.to.glp1","works.like.glp1")],
use = "pairwise.complete.obs"), 2)
## effect.like.glp1 similar.to.glp1 works.like.glp1
## effect.like.glp1 1.00 0.88 0.83
## similar.to.glp1 0.88 1.00 0.83
## works.like.glp1 0.83 0.83 1.00
# Product evaluation items
cat("\n=== Product evaluation scale (3 items) ===\n")
##
## === Product evaluation scale (3 items) ===
psych::alpha(eval_items2)
##
## Reliability analysis
## Call: psych::alpha(x = eval_items2)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.97 0.97 0.95 0.91 29 0.005 2.4 1.4 0.91
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.95 0.97 0.97
## Duhachek 0.96 0.97 0.98
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## purchase1 0.95 0.95 0.91 0.91 20 0.0084 NA 0.91
## purchase2 0.95 0.95 0.91 0.91 20 0.0086 NA 0.91
## appealing 0.95 0.95 0.91 0.91 20 0.0085 NA 0.91
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## purchase1 136 0.97 0.97 0.94 0.93 2.2 1.3
## purchase2 135 0.97 0.97 0.94 0.93 2.5 1.4
## appealing 136 0.97 0.97 0.94 0.93 2.6 1.5
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## purchase1 0.43 0.19 0.16 0.17 0.05 0.00
## purchase2 0.38 0.13 0.16 0.25 0.08 0.01
## appealing 0.37 0.13 0.15 0.23 0.12 0.00
cat("\nInter-item correlations:\n")
##
## Inter-item correlations:
round(cor(eval_items2, use = "pairwise.complete.obs"), 2)
## purchase1 purchase2 appealing
## purchase1 1.00 0.91 0.91
## purchase2 0.91 1.00 0.91
## appealing 0.91 0.91 1.00
# Endorsement rates (somewhat/strongly agree)
agree_pct <- function(x) round(100 * mean(x >= 4, na.rm = TRUE), 1)
cat("Endorsement rates (somewhat/strongly agree):\n")
## Endorsement rates (somewhat/strongly agree):
cat(" Has GLP-1-like effects:", agree_pct(dat2$effect.like.glp1), "%\n")
## Has GLP-1-like effects: 73.3 %
cat(" Similar to GLP-1: ", agree_pct(dat2$similar.to.glp1), "%\n")
## Similar to GLP-1: 72.8 %
cat(" Works like GLP-1: ", agree_pct(dat2$works.like.glp1), "%\n")
## Works like GLP-1: 66.2 %
p1 <- plot_likert_hist(dat2, "effect.like.glp1", "Has GLP-1-like effects")
p2 <- plot_likert_hist(dat2, "similar.to.glp1", "Similar to GLP-1")
p3 <- plot_likert_hist(dat2, "works.like.glp1", "Works like GLP-1", show_x = TRUE)
p1 / p2 / p3
Figure 2. Distributions of endorsement that the product has GLP-1-like effects (Study 2).
# False-belief strength by condition
m_belief2 <- aov(glp1_belief_mean ~ condition, data = dat2)
cat("=== Beliefs by condition ===\n"); summary(m_belief2)
## === Beliefs by condition ===
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 0.79 0.3961 0.319 0.728
## Residuals 133 165.19 1.2421
# Product evaluations by condition
m_prefer2 <- aov(eval_mean ~ condition, data = dat2)
cat("\n=== Evaluations by condition ===\n"); summary(m_prefer2)
##
## === Evaluations by condition ===
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 3.07 1.536 0.838 0.435
## Residuals 133 243.64 1.832
# WTP by condition
m_wtp2 <- aov(wtp ~ condition, data = dat2)
cat("\n=== WTP by condition ===\n"); summary(m_wtp2)
##
## === WTP by condition ===
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 1.55 0.7758 0.339 0.713
## Residuals 133 304.71 2.2910
m_demog2 <- lm(
glp1_belief_mean ~ age + bmi + income + ses +
gender_simple + race_simple + glp1_exposure + condition,
data = dat2
)
summary(m_demog2)
##
## Call:
## lm(formula = glp1_belief_mean ~ age + bmi + income + ses + gender_simple +
## race_simple + glp1_exposure + condition, data = dat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8974 -0.3772 0.1976 0.7078 1.7750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.624106 0.540030 8.563 4.6e-14 ***
## age -0.007032 0.007660 -0.918 0.3604
## bmi -0.011276 0.009288 -1.214 0.2271
## income -0.115808 0.078010 -1.485 0.1403
## ses 0.005298 0.074303 0.071 0.9433
## gender_simpleMale 0.101440 0.201124 0.504 0.6149
## race_simpleBlack 0.608175 0.359617 1.691 0.0934 .
## race_simpleAsian -0.064022 0.361348 -0.177 0.8597
## race_simpleOther 0.265451 0.346526 0.766 0.4452
## glp1_exposure -0.467977 0.380182 -1.231 0.2208
## conditionmock 0.103269 0.249163 0.414 0.6793
## conditionsmoothie_king -0.093377 0.257641 -0.362 0.7177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.1 on 119 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.1107, Adjusted R-squared: 0.02851
## F-statistic: 1.347 on 11 and 119 DF, p-value: 0.2075
# Table 3: product evaluations ~ beliefs + demographics
m_eval2 <- lm(
eval_mean ~ glp1_belief_mean + age + bmi + income + ses +
gender_simple + race_simple + glp1_exposure + condition,
data = dat2
)
summary(m_eval2)
##
## Call:
## lm(formula = eval_mean ~ glp1_belief_mean + age + bmi + income +
## ses + gender_simple + race_simple + glp1_exposure + condition,
## data = dat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4918 -1.1428 -0.2576 1.0383 2.6879
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.649329 0.808749 0.803 0.42366
## glp1_belief_mean 0.353666 0.107990 3.275 0.00139 **
## age -0.011536 0.009055 -1.274 0.20518
## bmi 0.018661 0.011009 1.695 0.09270 .
## income -0.188597 0.092746 -2.033 0.04425 *
## ses 0.201981 0.087534 2.307 0.02277 *
## gender_simpleMale 0.083236 0.237184 0.351 0.72626
## race_simpleBlack 0.046515 0.428702 0.109 0.91378
## race_simpleAsian -0.485641 0.425736 -1.141 0.25630
## race_simpleOther -0.178304 0.409225 -0.436 0.66384
## glp1_exposure 0.390963 0.450710 0.867 0.38746
## conditionmock 0.334189 0.293735 1.138 0.25754
## conditionsmoothie_king 0.186300 0.303677 0.613 0.54074
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.296 on 118 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.1864, Adjusted R-squared: 0.1037
## F-statistic: 2.253 on 12 and 118 DF, p-value: 0.01324
# Table 4: WTP ~ beliefs + demographics
m_wtp2_full <- lm(
wtp ~ glp1_belief_mean + age + bmi + income + ses +
gender_simple + race_simple + glp1_exposure + condition,
data = dat2
)
summary(m_wtp2_full)
##
## Call:
## lm(formula = wtp ~ glp1_belief_mean + age + bmi + income + ses +
## gender_simple + race_simple + glp1_exposure + condition,
## data = dat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6375 -0.9488 -0.2903 0.5736 3.7960
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.308214 0.882607 -0.349 0.72755
## glp1_belief_mean 0.373319 0.117852 3.168 0.00196 **
## age -0.019426 0.009882 -1.966 0.05167 .
## bmi 0.002903 0.012014 0.242 0.80950
## income -0.130691 0.101216 -1.291 0.19915
## ses 0.292870 0.095528 3.066 0.00269 **
## gender_simpleMale -0.289948 0.258844 -1.120 0.26492
## race_simpleBlack 0.021751 0.467852 0.046 0.96300
## race_simpleAsian -0.648367 0.464616 -1.395 0.16549
## race_simpleOther -0.281454 0.446597 -0.630 0.52977
## glp1_exposure 0.802286 0.491870 1.631 0.10554
## conditionmock 0.188067 0.320560 0.587 0.55854
## conditionsmoothie_king 0.100154 0.331410 0.302 0.76303
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.415 on 118 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.2107, Adjusted R-squared: 0.1304
## F-statistic: 2.624 on 12 and 118 DF, p-value: 0.003859
cor.test(dat2$script.required, dat2$glp1_belief_mean)
##
## Pearson's product-moment correlation
##
## data: dat2$script.required and dat2$glp1_belief_mean
## t = -1.7237, df = 134, p-value = 0.08707
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.30797884 0.02158659
## sample estimates:
## cor
## -0.1472815
Goal: Replicate Study 2 in a larger, higher-powered sample and add psychological and social-media predictors (N = 363).
s3_url <- "https://osf.io/ypztv/download" # update with OSF permalink
dat3 <- read.csv(s3_url, comment.char = "#")
dat3 <- dat3[-c(1:10), ]
dat3 <- subset(dat3, product.type == "Smoothie Supplement Powder")
# ── Core DVs ─────────────────────────────────────────────────────────────────
dat3$effect.like.glp1 <- unname(agree_map[na_99(dat3$effect.like.glp1)])
dat3$similar.to.glp1 <- unname(agree_map[na_99(dat3$similar.to.glp1)])
dat3$works.like.glp1 <- unname(agree_map[na_99(dat3$works.like.glp1)])
dat3$glp1_belief_mean <- rowMeans(
dat3[, c("effect.like.glp1","similar.to.glp1","works.like.glp1")],
na.rm = TRUE
)
dat3$purchase1_num <- unname(likely_map[na_99(dat3$purchase1)])
dat3$purchase2_num <- unname(agree_map[na_99(dat3$purchase2)])
dat3$appealing_num <- unname(agree_map[na_99(dat3$appealing)])
dat3$prefer <- rowMeans(
dat3[, c("purchase1_num","purchase2_num","appealing_num")],
na.rm = TRUE
)
dat3$wtp <- as.numeric(unname(wtp_map[na_99(dat3$wtp)]))
# Prescription-knowledge item
dat3$script.required <- unname(script_map[as.character(dat3$script.required)])
# ── Demographics ─────────────────────────────────────────────────────────────
dat3$age <- as_num_na99(dat3$age)
dat3$bmi <- (as_num_na99(dat3$weight) / (as_num_na99(dat3$height_1)^2)) * 703
dat3$ses <- as.numeric(sub("^([0-9]+).*", "\\1", na_99(dat3$ses)))
dat3$income <- factor(
na_99(dat3$income),
levels = c("$0 - $19,999","$20,000 - $39,999","$40,000 - $59,999",
"$60,000 - $89,999","$90,000 - $119,999",
"$120,000 - $149,999","$150,000+"),
ordered = TRUE
) |> as.numeric()
dat3$gender_simple <- case_when(
na_99(dat3$gender) %in% c("Female","Female,Non-binary / third gender") ~ "Female",
na_99(dat3$gender) %in% c("Male","Male,Non-binary / third gender","Male,Transgender") ~ "Male",
TRUE ~ NA_character_
) |> factor(levels = c("Female","Male"))
race3 <- na_99(dat3$race)
race3[grepl(",", race3) | grepl("prefer not to say", race3, ignore.case = TRUE)] <-
"Other / Prefer not to say"
dat3$race_simple <- case_when(
race3 == "White or European American" ~ "White",
race3 == "Black or African American" ~ "Black",
race3 == "Asian" ~ "Asian",
is.na(race3) ~ NA_character_,
TRUE ~ "Other"
) |> factor(levels = c("White","Black","Asian","Other"))
dat3$taking.glp1_num <- unname(yes_no_map[na_99(dat3$taking.glp1)])
dat3$glp1_exposure <- case_when(
dat3$taking.glp1_num %in% c(4,5) ~ 1,
dat3$taking.glp1_num %in% c(1,2) ~ 0,
TRUE ~ NA_real_
)
dat3$politics <- recode(
na_99(dat3$politics),
"Very liberal" = "1",
"Liberal" = "2",
"Somewhat liberal" = "3",
"Neither liberal nor conservative"= "4",
"Somewhat conservative" = "5",
"Conservative" = "6",
"Very conservative" = "7",
.default = NA_character_
) |> as.numeric()
dat3$FL_16_DO <- factor(na_99(dat3$FL_16_DO))
# ── Health anxiety (6 SHAI items, columns 40–45) ─────────────────────────────
health_anxt <- as.data.frame(dat3[, 40:45])
score_item_0to3 <- function(x, levels_key) {
x <- na_99(x)
as.numeric(factor(x, levels = levels_key, ordered = TRUE)) - 1
}
s1h <- score_item_0to3(health_anxt[[1]], c("Never","Occasionally","Much of the time","Most of the time"))
s2h <- score_item_0to3(health_anxt[[2]], c("Without a problem","Most of the time",
"I try to resist thoughts of illness but am often unable to do so",
"Thoughts of illness are so strong I no longer try to resist them"))
s3h <- score_item_0to3(health_anxt[[3]], c("Not at all","Sometimes","Often","Always"))
s4h <- score_item_0to3(health_anxt[[4]], c("Never","Sometimes","Often","Always"))
s5h <- score_item_0to3(health_anxt[[5]], c("Very low","Fairly low","Moderate","High"))
s6h <- score_item_0to3(health_anxt[[6]], c("Do not worry about my health",
"Have a normal attitude about my health",
"Worry too much about my health","Am a hypochondriac"))
dat3$ha_mean <- rowMeans(cbind(s1h, s2h, s3h, s4h, s5h, s6h), na.rm = TRUE)
# ── Weight loss concerns (5 items, columns 46–52) ────────────────────────────
dat3$afraid_3lbs <- recode(dat3[,46],
"Very unafraid" = 1, "Unafraid" = 2,
"Neither afraid nor unafraid" = 3, "Afraid" = 4, "Very afraid" = 5)
dat3$weight_importance <- as.numeric(factor(dat3[,47],
levels = c("Not important at all","Barely important","Somewhat important",
"Important","Very important")))
dat3[,48][dat3[,48] == "-99"] <- NA
dat3$ever_dieted <- as.numeric(factor(dat3[,48],
levels = c("Definitely not","Probably not","Might or might not",
"Probably yes","Definitely yes")))
dat3$need.lose.weight <- ifelse(dat3[,51] == "Yes", 1, 0)
dat3$scale_matters <- as.numeric(factor(dat3[,52],
levels = c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
"Somewhat agree","Strongly agree")))
dat3$wlc.average <- rowMeans(cbind(
dat3$afraid_3lbs, dat3$weight_importance, dat3$ever_dieted,
dat3$need.lose.weight, dat3$scale_matters), na.rm = TRUE)
# ── Body surveillance (7 items, columns 53–59) ───────────────────────────────
bs_items <- names(dat3)[53:59]
likert7_key <- c("Strongly disagree"=1,"Disagree"=2,"Somewhat disagree"=3,
"Neither agree nor disagree"=4,"Somewhat agree"=5,
"Agree"=6,"Strongly agree"=7)
for (v in bs_items) {
x <- na_99(dat3[[v]]); x[x == "Not Applicable"] <- NA
dat3[[v]] <- unname(likert7_key[x])
}
rev_bs <- intersect(c("bs1","bs2","bs3","bs6","bs7"), names(dat3))
dat3[rev_bs] <- lapply(dat3[rev_bs], function(x) ifelse(is.na(x), NA, 8 - x))
dat3$body_surv_mean <- rowMeans(dat3[, bs_items], na.rm = TRUE)
# ── Body appreciation (10 items, columns 60–69) ──────────────────────────────
ba_items <- names(dat3)[60:69]
ba_key <- c("Never"=1,"Seldom"=2,"Sometimes"=3,"Often"=4,"Always"=5)
for (v in ba_items) {
x <- na_99(dat3[[v]]); x[x == "Not Applicable"] <- NA
dat3[[v]] <- unname(ba_key[x])
}
dat3$body_appreciation_mean <- rowMeans(dat3[, ba_items], na.rm = TRUE)
# ── Eating disorder risk (7 yes/no items, columns 70–76) ─────────────────────
ed_items <- names(dat3)[70:76]
yn_key <- c("No"=0,"Yes"=1)
for (v in ed_items) {
x <- na_99(dat3[[v]]); x[x == "Prefer not to say"] <- NA
dat3[[v]] <- unname(yn_key[x])
}
dat3$ed_risk_total <- rowSums(dat3[, ed_items], na.rm = TRUE)
cat("=== False-belief scale ===\n")
## === False-belief scale ===
psych::alpha(dat3[, c("effect.like.glp1","similar.to.glp1","works.like.glp1")])[1]
## $total
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd
## 0.9540386 0.9545878 0.938384 0.8751068 21.02053 0.004208446 3.617069 1.230947
## median_r
## 0.8806055
cat("\n=== Product evaluation scale ===\n")
##
## === Product evaluation scale ===
psych::alpha(dat3[, c("purchase1_num","purchase2_num","appealing_num")])[1]
## $total
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd
## 0.9669932 0.9680632 0.9563041 0.909942 30.31186 0.003014521 2.356122 1.368216
## median_r
## 0.9028927
cat("\n=== Health anxiety (SHAI) ===\n")
##
## === Health anxiety (SHAI) ===
psych::alpha(cbind(s1h, s2h, s3h, s4h, s5h, s6h))[1]
## $total
## raw_alpha std.alpha G6(smc) average_r S/N ase mean
## 0.8883078 0.8887597 0.8764536 0.5711083 7.989546 0.008653013 1.004736
## sd median_r
## 0.6093592 0.6029069
cat("\n=== Weight loss concerns ===\n")
##
## === Weight loss concerns ===
psych::alpha(cbind(dat3$afraid_3lbs, dat3$weight_importance,
dat3$ever_dieted, dat3$need.lose.weight,
dat3$scale_matters))[1]
## $total
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd
## 0.7540622 0.7833726 0.7620642 0.4196992 3.616221 0.01878833 2.884718 0.7933238
## median_r
## 0.3934151
cat("\n=== Body surveillance ===\n")
##
## === Body surveillance ===
psych::alpha(dat3[, bs_items])[1]
## $total
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd
## 0.8437645 0.842784 0.837652 0.433688 5.360677 0.01221992 4.079325 1.18841
## median_r
## 0.4128959
cat("\n=== Body appreciation ===\n")
##
## === Body appreciation ===
psych::alpha(dat3[, ba_items])[1]
## $total
## raw_alpha std.alpha G6(smc) average_r S/N ase mean
## 0.9619208 0.9619966 0.9623365 0.7168219 25.31346 0.002852702 3.333065
## sd median_r
## 0.9373581 0.7279823
cat("\n=== Eating disorder risk ===\n")
##
## === Eating disorder risk ===
psych::alpha(dat3[, ed_items])[1]
## $total
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd
## 0.7669352 0.7399519 0.743659 0.2890112 2.845443 0.01600541 0.2642601 0.2705937
## median_r
## 0.2931444
cat("Endorsement rates (somewhat/strongly agree):\n")
## Endorsement rates (somewhat/strongly agree):
cat(" Has GLP-1-like effects:", agree_pct(dat3$effect.like.glp1), "%\n")
## Has GLP-1-like effects: 67 %
cat(" Similar to GLP-1: ", agree_pct(dat3$similar.to.glp1), "%\n")
## Similar to GLP-1: 68.9 %
cat(" Works like GLP-1: ", agree_pct(dat3$works.like.glp1), "%\n")
## Works like GLP-1: 60.1 %
p1 <- plot_likert_hist(dat3, "effect.like.glp1", "Has GLP-1-like effects")
p2 <- plot_likert_hist(dat3, "similar.to.glp1", "Similar to GLP-1")
p3 <- plot_likert_hist(dat3, "works.like.glp1", "Works like GLP-1", show_x = TRUE)
p1 / p2 / p3
Figure 3. Distributions of endorsement that the product has GLP-1-like effects (Study 3).
m_belief3 <- aov(glp1_belief_mean ~ FL_16_DO, data = dat3)
cat("=== Beliefs by condition ===\n"); summary(m_belief3)
## === Beliefs by condition ===
## Df Sum Sq Mean Sq F value Pr(>F)
## FL_16_DO 2 19.7 9.827 6.684 0.00141 **
## Residuals 370 544.0 1.470
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nHolm-corrected pairwise contrasts:\n")
##
## Holm-corrected pairwise contrasts:
emmeans(m_belief3, pairwise ~ FL_16_DO, adjust = "holm")
## $emmeans
## FL_16_DO emmean SE df lower.CL upper.CL
## AmazonSmoothie 3.29 0.111 370 3.07 3.50
## Smoothie 3.73 0.106 370 3.52 3.94
## SmoothieKing 3.82 0.110 370 3.60 4.03
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## AmazonSmoothie - Smoothie -0.4453 0.153 370 -2.905 0.0078
## AmazonSmoothie - SmoothieKing -0.5312 0.156 370 -3.400 0.0022
## Smoothie - SmoothieKing -0.0859 0.152 370 -0.564 0.5731
##
## P value adjustment: holm method for 3 tests
m_prefer3 <- aov(prefer ~ FL_16_DO, data = dat3)
cat("\n=== Evaluations by condition ===\n"); summary(m_prefer3)
##
## === Evaluations by condition ===
## Df Sum Sq Mean Sq F value Pr(>F)
## FL_16_DO 2 9.2 4.618 2.486 0.0846 .
## Residuals 370 687.2 1.857
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m_wtp3 <- aov(wtp ~ FL_16_DO, data = dat3)
cat("\n=== WTP by condition ===\n"); summary(m_wtp3)
##
## === WTP by condition ===
## Df Sum Sq Mean Sq F value Pr(>F)
## FL_16_DO 2 7.0 3.476 1.576 0.208
## Residuals 370 815.9 2.205
demog_vars3 <- c("bmi","income","ses","age","gender_simple","race_simple",
"glp1_exposure","FL_16_DO")
psych_vars3 <- c("ha_mean","wlc.average","body_surv_mean",
"body_appreciation_mean","ed_risk_total","politics")
fml <- function(dv, rhs) as.formula(paste(dv, "~", paste(rhs, collapse = " + ")))
make_df <- function(dv, rhs) {
dat3 |> select(all_of(c(dv, rhs))) |> filter(complete.cases(pick(everything())))
}
df1 <- make_df("glp1_belief_mean", c(demog_vars3, psych_vars3))
message("Table 5 N: ", nrow(df1))
m1_s3 <- lm(fml("glp1_belief_mean", c(demog_vars3, psych_vars3)), data = df1)
summary(m1_s3)
##
## Call:
## lm(formula = fml("glp1_belief_mean", c(demog_vars3, psych_vars3)),
## data = df1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1386 -0.6162 0.1807 0.8567 2.5177
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.341447 0.572346 4.091 5.35e-05 ***
## bmi 0.008107 0.006049 1.340 0.181066
## income -0.036951 0.045738 -0.808 0.419714
## ses -0.042250 0.049081 -0.861 0.389929
## age -0.002025 0.005202 -0.389 0.697363
## gender_simpleMale -0.193812 0.132526 -1.462 0.144527
## race_simpleBlack 0.327344 0.204751 1.599 0.110792
## race_simpleAsian 0.350127 0.270356 1.295 0.196164
## race_simpleOther 0.430129 0.240495 1.789 0.074567 .
## glp1_exposure -0.611901 0.180255 -3.395 0.000767 ***
## FL_16_DOSmoothie 0.431067 0.151529 2.845 0.004709 **
## FL_16_DOSmoothieKing 0.465667 0.153938 3.025 0.002672 **
## ha_mean -0.051689 0.122650 -0.421 0.673696
## wlc.average 0.092939 0.108461 0.857 0.392101
## body_surv_mean -0.035231 0.064102 -0.550 0.582938
## body_appreciation_mean 0.200286 0.082902 2.416 0.016211 *
## ed_risk_total 0.041609 0.045327 0.918 0.359279
## politics 0.123817 0.034594 3.579 0.000394 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.16 on 346 degrees of freedom
## Multiple R-squared: 0.1608, Adjusted R-squared: 0.1196
## F-statistic: 3.9 on 17 and 346 DF, p-value: 4.992e-07
# Table 6: product evaluations ~ beliefs + demographics + psych
df2 <- make_df("prefer", c("glp1_belief_mean", demog_vars3, psych_vars3))
message("Table 6 N: ", nrow(df2))
m2_s3 <- lm(fml("prefer", c("glp1_belief_mean", demog_vars3, psych_vars3)), data = df2)
summary(m2_s3)
##
## Call:
## lm(formula = fml("prefer", c("glp1_belief_mean", demog_vars3,
## psych_vars3)), data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4716 -0.8689 -0.1019 0.8756 3.8496
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.122817 0.601357 -1.867 0.06273 .
## glp1_belief_mean 0.295015 0.055167 5.348 1.62e-07 ***
## bmi 0.008547 0.006224 1.373 0.17055
## income -0.014465 0.046979 -0.308 0.75834
## ses -0.033474 0.050419 -0.664 0.50719
## age 0.004316 0.005339 0.808 0.41946
## gender_simpleMale -0.081471 0.136413 -0.597 0.55074
## race_simpleBlack 0.193602 0.210883 0.918 0.35923
## race_simpleAsian 0.096632 0.278102 0.347 0.72845
## race_simpleOther -0.594601 0.247925 -2.398 0.01700 *
## glp1_exposure 0.372204 0.188026 1.980 0.04855 *
## FL_16_DOSmoothie 0.278354 0.157302 1.770 0.07769 .
## FL_16_DOSmoothieKing 0.262888 0.160041 1.643 0.10137
## ha_mean -0.108318 0.125891 -0.860 0.39016
## wlc.average 0.201815 0.111417 1.811 0.07096 .
## body_surv_mean 0.192970 0.065808 2.932 0.00359 **
## body_appreciation_mean 0.078692 0.085785 0.917 0.35962
## ed_risk_total 0.145039 0.046570 3.114 0.00200 **
## politics 0.053479 0.036150 1.479 0.13996
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.19 on 345 degrees of freedom
## Multiple R-squared: 0.2848, Adjusted R-squared: 0.2475
## F-statistic: 7.634 on 18 and 345 DF, p-value: < 2.2e-16
# Table 7: WTP ~ beliefs + demographics + psych
df3 <- make_df("wtp", c("glp1_belief_mean", demog_vars3, psych_vars3))
message("Table 7 N: ", nrow(df3))
m3_s3 <- lm(fml("wtp", c("glp1_belief_mean", demog_vars3, psych_vars3)), data = df3)
summary(m3_s3)
##
## Call:
## lm(formula = fml("wtp", c("glp1_belief_mean", demog_vars3, psych_vars3)),
## data = df3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5565 -0.9540 -0.3298 0.6195 4.6715
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.065819 0.704677 -2.932 0.0036 **
## glp1_belief_mean 0.347959 0.064645 5.383 1.36e-07 ***
## bmi 0.008376 0.007293 1.148 0.2516
## income 0.016361 0.055050 0.297 0.7665
## ses 0.024145 0.059082 0.409 0.6830
## age 0.006419 0.006256 1.026 0.3056
## gender_simpleMale 0.063314 0.159851 0.396 0.6923
## race_simpleBlack 0.435027 0.247115 1.760 0.0792 .
## race_simpleAsian -0.252520 0.325883 -0.775 0.4389
## race_simpleOther -0.392473 0.290521 -1.351 0.1776
## glp1_exposure -0.012206 0.220331 -0.055 0.9559
## FL_16_DOSmoothie 0.254049 0.184328 1.378 0.1690
## FL_16_DOSmoothieKing 0.007347 0.187538 0.039 0.9688
## ha_mean -0.192958 0.147521 -1.308 0.1917
## wlc.average 0.111755 0.130559 0.856 0.3926
## body_surv_mean 0.092216 0.077114 1.196 0.2326
## body_appreciation_mean 0.066048 0.100524 0.657 0.5116
## ed_risk_total 0.133020 0.054571 2.438 0.0153 *
## politics -0.009400 0.042361 -0.222 0.8245
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.394 on 345 degrees of freedom
## Multiple R-squared: 0.1756, Adjusted R-squared: 0.1326
## F-statistic: 4.084 on 18 and 345 DF, p-value: 9.262e-08
cor.test(dat3$script.required, dat3$glp1_belief_mean)
##
## Pearson's product-moment correlation
##
## data: dat3$script.required and dat3$glp1_belief_mean
## t = -2.0668, df = 371, p-value = 0.03944
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.206003181 -0.005206382
## sample estimates:
## cor
## -0.1066925
sessionInfo()
## R version 4.5.1 (2025-06-13)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.7.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## 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: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] psych_2.5.6 emmeans_1.11.2 patchwork_1.3.2 ggplot2_3.5.2
## [5] dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] sandwich_3.1-1 sass_0.4.10 generics_0.1.4 lattice_0.22-7
## [5] digest_0.6.37 magrittr_2.0.3 evaluate_1.0.4 grid_4.5.1
## [9] estimability_1.5.1 RColorBrewer_1.1-3 mvtnorm_1.3-3 fastmap_1.2.0
## [13] jsonlite_2.0.0 Matrix_1.7-3 survival_3.8-3 multcomp_1.4-28
## [17] scales_1.4.0 TH.data_1.1-4 codetools_0.2-20 jquerylib_0.1.4
## [21] mnormt_2.1.1 cli_3.6.5 rlang_1.1.6 splines_4.5.1
## [25] withr_3.0.2 cachem_1.1.0 yaml_2.3.10 tools_4.5.1
## [29] parallel_4.5.1 coda_0.19-4.1 vctrs_0.6.5 R6_2.6.1
## [33] zoo_1.8-14 lifecycle_1.0.4 MASS_7.3-65 pkgconfig_2.0.3
## [37] pillar_1.11.0 bslib_0.9.0 gtable_0.3.6 glue_1.8.0
## [41] xfun_0.56 tibble_3.3.0 tidyselect_1.2.1 rstudioapi_0.17.1
## [45] knitr_1.50 farver_2.1.2 xtable_1.8-4 htmltools_0.5.8.1
## [49] nlme_3.1-168 rmarkdown_2.29 labeling_0.4.3 compiler_4.5.1