Overview

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/.


Setup

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()
}

Study 1

Goal: Test how often consumers infer that a GLP-1 “support” smoothie powder literally contains GLP-1 as an ingredient.

Data import & cleaning

# ── 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")

Descriptive statistics (Table 1)

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)."
)
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.


Study 2

Goal: Replicate the ingredient-inference effect using graded dose-response belief measures and three advertisement stimuli (N = 136).

Data import & cleaning

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)

Scale reliability

# 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

Descriptive statistics & Figure 2

# 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).

Figure 2. Distributions of endorsement that the product has GLP-1-like effects (Study 2).

Condition differences (ANOVAs)

# 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

Demographic predictors of false beliefs (Table 2)

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

False beliefs and product evaluations (Tables 3 & 4)

# 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

Prescription knowledge and false beliefs

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

Study 3

Goal: Replicate Study 2 in a larger, higher-powered sample and add psychological and social-media predictors (N = 363).

Data import & cleaning

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)

Scale reliability

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

Descriptive statistics & Figure 3

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).

Figure 3. Distributions of endorsement that the product has GLP-1-like effects (Study 3).

Condition differences (ANOVAs)

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

Demographic and psychological predictors of false beliefs (Table 5)

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

False beliefs and product evaluations (Tables 6 & 7)

# 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

Prescription knowledge and false beliefs

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

Session Info

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