In brentkaplan/beezdemand: Behavioral Economic Easy Demand

knitr::opts_chunk$set(
  message = FALSE,
  collapse = TRUE,
  comment = NA,
  fig.width = 6,
  fig.height = 4
)

library(beezdemand)
library(dplyr)
library(purrr)
library(tidyr)
library(ggplot2)

load("../data/etm.RData")

Data Structure

glimpse(etm)

Typical columns: - id: participant identifier

• x: alternative product price

• y: consumption level

• target: condition type (e.g., "alt")

• group: product category

Checking Unsystematic Data

etm |>
  dplyr::filter(group %in% "E-Cigarettes" & id %in% 1)


unsys_one <- etm |>
  filter(group %in% "E-Cigarettes" & id %in% 1) |>
  check_unsystematic_cp()

unsys_one

unsys_all <- etm |>
  group_by(id, group) |>
  nest() |>
  mutate(unsys = map(data, check_unsystematic_cp)) |>
  unnest(unsys)

unsys_all

summary(unsys_all)

Nonlinear Model Fitting

Two Stage

fit_one <- etm |>
  dplyr::filter(group %in% "E-Cigarettes" & id %in% 1) |>
  fit_cp_nls(
    equation = "exponentiated",
    return_all = TRUE
  )

summary(fit_one)

plot(fit_one, x_trans = "log10")

fit_all <- etm |>
  group_by(id, group) |>
  nest() |>
  mutate(
    unsys = map(data, check_unsystematic_cp),
    fit = map(data, fit_cp_nls, equation = "exponentiated", return_all = TRUE),
    summary = map(fit, summary),
    plot = map(fit, plot, x_trans = "log10"),
    glance = map(fit, glance),
    tidy = map(fit, tidy)
  )

fit_all

fit_all$summary[[2]]

rand_id <- sample(unique(etm$id), 1)
print(rand_id)
fit_all$summary[[rand_id]]
fit_all$plot[[rand_id]]

Fit to Group (pooled by group)

fit_pooled <- etm |>
  group_by(group) |>
  nest() |>
  mutate(
    unsys = map(data, check_unsystematic_cp),
    fit = map(data, fit_cp_nls, equation = "exponentiated", return_all = TRUE),
    summary = map(fit, summary),
    plot = map(fit, plot, x_trans = "log10"),
    glance = map(fit, glance),
    tidy = map(fit, tidy)
  )

fit_pooled

fit_pooled |>
  dplyr::filter(group %in% "E-Cigarettes") |>
  pull(summary)

fit_pooled |>
  dplyr::filter(group %in% "E-Cigarettes") |>
  pull(plot)

Fit to Group (mean)

fit_mean <- etm |>
  group_by(group, x) |>
  summarise(
    y = mean(y)
  ) |>
  ungroup() |>
  group_by(group) |>
  nest() |>
  mutate(
    unsys = map(data, check_unsystematic_cp),
    fit = map(data, fit_cp_nls, equation = "exponentiated", return_all = TRUE),
    summary = map(fit, summary),
    plot = map(fit, plot, x_trans = "log10"),
    glance = map(fit, glance),
    tidy = map(fit, tidy)
  )

fit_mean

fit_mean |>
    unnest(cols = c(glance, tidy)) |>
    select(
        group,
        term,
        estimate
    ) |>
        ggplot(aes(x = group, y = estimate, group = term)) +
        geom_bar(stat = "identity") +
        facet_wrap(~term)

fit_mean |>
  dplyr::filter(group %in% "E-Cigarettes") |>
  pull(summary)

fit_mean |>
  dplyr::filter(group %in% "E-Cigarettes") |>
  pull(plot)

Linear Model Fitting

fit_one_linear <- etm |>
  dplyr::filter(group %in% "E-Cigarettes" & id %in% 1) |>
  fit_cp_linear(
    type = "fixed",
    log10x = TRUE,
    return_all = TRUE
  )

summary(fit_one_linear)
plot(fit_one_linear, x_trans = "log10")

Linear Mixed-Effects Model

fit_mixed <- fit_cp_linear(
  etm,
  type = "mixed",
  log10x = TRUE,
  group_effects = "interaction",
  random_slope = FALSE,
  return_all = TRUE
)

summary(fit_mixed)

# plot fixed effects only
plot(fit_mixed, x_trans = "log10", pred_type = "fixed")

# plot random effects only
plot(fit_mixed, x_trans = "log10", pred_type = "random")

# plot both fixed and random effects
plot(fit_mixed, x_trans = "log10", pred_type = "all")

Extracting Model Coefficients

glance(fit_one)
tidy(fit_one)

extract_coefficients(fit_mixed)

Post-hoc Estimated Marginal Means and Comparisons

cp_posthoc_slopes(fit_mixed)
cp_posthoc_intercepts(fit_mixed)

brentkaplan/beezdemand documentation built on June 11, 2025, 3:50 a.m.