---
title: "How to Use Cross-Price Demand Model Functions"
author: "Brent Kaplan"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{How to Use Cross-Price Demand Model Functions}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
    message = FALSE,
    warning = FALSE,
    include = TRUE,
    collapse = TRUE,
    comment = "#>",
    dev = "png",
    dpi = 144,
    fig.width = 7,
    fig.height = 5
)

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

# Load package datasets
data("etm", package = "beezdemand")
data("lowNicClean", package = "beezdemand")
data("cannabisCigarettes", package = "beezdemand")
data("ongoingETM", package = "beezdemand")

lnic <- lowNicClean |>
    mutate(
        target = if_else(commodity == "cigarettesAdj", "adjusting", "fixed"),
        group = "cigarettes"
    ) |>
    select(
        -commodity
    ) |>
    relocate(condition, .after = group)

can_cig <- cannabisCigarettes |>
    select(
        id,
        x,
        y,
        target = commodity
    )

ongoing_etm <- ongoingETM |>
    pivot_longer(
        -c(id, x, flavor),
        names_to = "target",
        values_to = "y",
        values_drop_na = TRUE
    ) |>
    select(id, x, y, target)
```

## Introduction

Cross-price demand analysis examines how consumption of one commodity changes
as the price of another commodity varies. This is central to understanding
economic relationships between goods:

- **Substitutes**: When the price of the target commodity increases and
  consumption of the alternative increases, the goods function as substitutes
  (e.g., e-cigarettes and combustible cigarettes).
- **Complements**: When the price of the target increases and consumption of the
  alternative *decreases*, the goods function as complements.
- **Independent**: When the price of one commodity does not meaningfully affect
  consumption of the other.

These relationships are quantified through cross-price elasticity parameters.
The `beezdemand` package provides nonlinear (`fit_cp_nls()`), linear
(`fit_cp_linear()`), and mixed-effects linear (`fit_cp_linear(type = "mixed")`)
approaches for cross-price modeling.

## Data Structure

```{r}
glimpse(etm)
```

Typical columns:
- `id`: participant identifier

- `x`: alternative product price

- `y`: consumption level

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

- `group`: product category

These are the default canonical column names. If your data uses different names,
pass `x_var`, `y_var`, `id_var`, `group_var`, or `target_var` arguments to the
fitting functions (e.g., `fit_cp_nls(data, x_var = "price", y_var = "qty")`).

## Complete Cross-Price Analysis Workflow

This section demonstrates a complete analysis workflow using data that contains multiple
experimental conditions: target consumption when the alternative is absent (`alone`),
target consumption when the alternative is present (`own`), and alternative consumption
as a function of target price (`alt`).

### Loading the cp Dataset

```{r load-cp}
# Load the cross-price example dataset
data("cp", package = "beezdemand")

# Examine structure
glimpse(cp)

# View conditions
table(cp$target)
```

The `cp` dataset contains:

- **alone**: Target commodity (cigarettes) consumption when alternative is *not* available
- **own**: Target commodity consumption when alternative *is* available
- **alt**: Alternative commodity (e-cigarettes) consumption as a function of target price

Note that `x` represents the **target commodity price** throughout all conditions.

### Step 1: Fit Target Demand (Alone Condition)

First, we fit a standard demand curve to the target commodity when the alternative is absent:

```{r fit-alone}
# Filter to alone condition
alone_data <- cp |>
    dplyr::filter(target == "alone")

# Fit demand curve (modern interface)
fit_alone <- fit_demand_fixed(
    data = alone_data,
    equation = "koff",
    k = 2
)

# View results
fit_alone
```

### Step 2: Fit Target Demand (Own Condition)

Next, fit the same demand model to target consumption when the alternative is present:

```{r fit-own}
# Filter to own condition
own_data <- cp |>
    dplyr::filter(target == "own")

# Fit demand curve
fit_own <- fit_demand_fixed(
    data = own_data,
    equation = "koff",
    k = 2
)

# View results
fit_own
```

### Step 3: Fit Cross-Price Model (Alt Condition)

Finally, fit the cross-price model to alternative consumption as a function of target price:

```{r fit-alt-cp}
# Filter to alt condition
alt_data <- cp |>
    dplyr::filter(target == "alt")

# Fit cross-price model
fit_alt <- fit_cp_nls(
    data = alt_data,
    equation = "exponentiated",
    return_all = TRUE
)

# View results
summary(fit_alt)
```

`fit_cp_nls()` uses a log10-parameterized optimizer internally (for numerical
stability), but `predict()` returns `y_pred` on the natural `y` scale.
For the `"exponential"` form, predictions may also include `y_pred_log10`.

### Comparing Results Across Conditions

```{r compare-conditions}
# Extract key parameters for each condition
coef_alone <- coef(fit_alone)
Q0_alone <- coef_alone$estimate[coef_alone$term == "q0"]
Alpha_alone <- coef_alone$estimate[coef_alone$term == "alpha"]

coef_own <- coef(fit_own)
Q0_own <- coef_own$estimate[coef_own$term == "q0"]
Alpha_own <- coef_own$estimate[coef_own$term == "alpha"]

comparison <- data.frame(
    Condition = c("Alone (Target)", "Own (Target)", "Alt (Cross-Price)"),
    Q0_or_Qalone = c(
        Q0_alone,
        Q0_own,
        coef(fit_alt)["qalone"]
    ),
    Alpha_or_I = c(
        Alpha_alone,
        Alpha_own,
        coef(fit_alt)["I"]
    )
)

comparison
```

**Interpretation:**

- Comparing `Q0` between alone and own conditions shows how the presence of an alternative
  affects baseline consumption of the target commodity
- The `I` parameter from the cross-price model indicates whether the products are
  substitutes (I < 0) or complements (I > 0)

### Combined Visualization

```{r combined-plot, fig.width=8, fig.height=4}
# Create prediction data
x_seq <- seq(0.01, max(cp$x), length.out = 100)

# Get demand predictions for each condition
pred_alone <- predict(fit_alone, newdata = data.frame(x = x_seq))$.fitted
pred_own <- predict(fit_own, newdata = data.frame(x = x_seq))$.fitted

# Cross-price model predictions (always on the natural y scale)
pred_alt <- predict(fit_alt, newdata = data.frame(x = x_seq))$y_pred

# Combine into plot data
plot_data <- data.frame(
    x = rep(x_seq, 3),
    y = c(pred_alone, pred_own, pred_alt),
    Condition = rep(c("Target (Alone)", "Target (Own)", "Alternative (Cross-Price)"), each = length(x_seq))
)

# Plot
ggplot() +
    geom_point(data = cp, aes(x = x, y = y, color = target), alpha = 0.6) +
    geom_line(data = plot_data, aes(x = x, y = y, color = Condition), linewidth = 1) +
    scale_x_log10() +
    scale_y_log10() +
    labs(
        x = "Target Price",
        y = "Consumption",
        title = "Cross-Price Analysis: All Conditions",
        color = "Condition"
    ) +
    theme_bw()
```

## Checking Unsystematic Data

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


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

unsys_one$results

unsys_one_lnic <- lnic |>
    filter(
        target == "adjusting",
        id == "R_00Q12ahGPKuESBT"
    ) |>
    check_systematic_cp()

unsys_one_lnic$results
```

```{r compute-unsys-etm, results='hide'}
unsys_all <- etm |>
    group_by(id, group) |>
    nest() |>
    mutate(
        sys = map(data, check_systematic_cp),
        results = map(sys, ~ dplyr::select(.x$results, -id))
    ) |>
    select(-data, -sys) |>
    unnest(results)
```

```{r show-unsys-etm}
knitr::kable(
    unsys_all |>
        group_by(group) |>
        summarise(
            n_subjects = n(),
            pct_systematic = round(mean(systematic, na.rm = TRUE) * 100, 1),
            .groups = "drop"
        ),
    caption = "Systematicity check by product group (ETM dataset)"
)
```

### Demonstration from Rzeszutek et al. (2025)

#### Low Nicotine Study (Kaplan et al., 2018)

```{r compute-unsys-lnic, results='hide'}
unsys_all_lnic <- lnic |>
    filter(target == "fixed") |>
    group_by(id, condition) |>
    nest() |>
    mutate(
        sys = map(
            data,
            check_systematic_cp
        )
    ) |>
    mutate(results = map(sys, ~ dplyr::select(.x$results, -id))) |>
    select(-data, -sys) |>
    unnest(results) |>
    arrange(id)
```

```{r show-unsys-lnic}
knitr::kable(
    unsys_all_lnic |>
        group_by(condition) |>
        summarise(
            n_subjects = n(),
            pct_systematic = round(mean(systematic, na.rm = TRUE) * 100, 1),
            .groups = "drop"
        ),
    caption = "Systematicity check by condition (Low Nicotine study, Kaplan et al., 2018)"
)
```


#### Unpublished Cannabis and Cigarette Data

```{r compute-unsys-can-cig, results='hide'}
unsys_all_can_cig <- can_cig |>
    filter(target %in% c("cannabisFix", "cigarettesFix")) |>
    group_by(id, target) |>
    nest() |>
    mutate(
        sys = map(data, check_systematic_cp),
        results = map(sys, ~ dplyr::select(.x$results, -id))
    ) |>
    select(-data, -sys) |>
    unnest(results) |>
    arrange(id)
```

```{r show-unsys-can-cig}
knitr::kable(
    unsys_all_can_cig |>
        group_by(target) |>
        summarise(
            n_subjects = n(),
            pct_systematic = round(mean(systematic, na.rm = TRUE) * 100, 1),
            .groups = "drop"
        ),
    caption = "Systematicity check by target (Cannabis/Cigarettes, unpublished data)"
)
```

#### In Progress Experimental Tobacco Marketplace Data

```{r compute-unsys-ongoing-etm, results='hide'}
unsys_all_ongoing_etm <- ongoing_etm |>
    # one person is doubled up
    distinct() |>
    filter(target %in% c("FixCig", "ECig")) |>
    group_by(id, target) |>
    nest() |>
    mutate(
        sys = map(data, check_systematic_cp),
        results = map(sys, ~ dplyr::select(.x$results, -id))
    ) |>
    select(-data, -sys) |>
    unnest(results)
```

```{r show-unsys-ongoing-etm}
knitr::kable(
    unsys_all_ongoing_etm |>
        group_by(target) |>
        summarise(
            n_subjects = n(),
            pct_systematic = round(mean(systematic, na.rm = TRUE) * 100, 1),
            .groups = "drop"
        ),
    caption = "Systematicity check by target (Ongoing ETM data)"
)
```

## Nonlinear Model Fitting

### Two Stage

```{r}
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")
```

```{r fit-individual, results='hide'}
fit_all <- etm |>
    group_by(id, group) |>
    nest() |>
    mutate(
        unsys = map(data, check_systematic_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)
    )
```

```{r show-individual-fits}
# Show parameter estimates for first 3 subjects only
knitr::kable(
    fit_all |>
        slice(1:3) |>
        unnest(tidy) |>
        select(id, group, term, estimate, std.error),
    digits = 3,
    caption = "Example parameter estimates (first 3 subjects)"
)

# Show one example plot
fit_all$plot[[2]]
```

### Fit to Group (pooled by group)

```{r fit-pooled, results='hide'}
fit_pooled <- etm |>
    group_by(group) |>
    nest() |>
    mutate(
        unsys = map(data, check_systematic_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)
    )
```

```{r show-pooled-results}
# Show tidy results instead of summary
knitr::kable(
    fit_pooled |>
        unnest(tidy) |>
        select(group, term, estimate, std.error),
    digits = 3,
    caption = "Pooled model parameter estimates by product group"
)

# Show one plot example
fit_pooled |>
    dplyr::filter(group == "E-Cigarettes") |>
    dplyr::pull(plot) |>
    pluck(1)
```

### Fit to Group (mean)

```{r fit-mean, results='hide'}
fit_mean <- etm |>
    group_by(group, x) |>
    summarise(
        y = mean(y),
        .groups = "drop"
    ) |>
    group_by(group) |>
    nest() |>
    mutate(
        unsys = map(data, check_systematic_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)
    )
```

```{r show-mean-results}
# Show tidy results
knitr::kable(
    fit_mean |>
        unnest(tidy) |>
        select(group, term, estimate, std.error),
    digits = 3,
    caption = "Mean model parameter estimates by product group"
)

# Show parameter estimates plot
fit_mean |>
    unnest(cols = c(glance, tidy)) |>
    select(
        group,
        term,
        estimate
    ) |>
    ggplot(aes(x = group, y = estimate, group = term)) +
    geom_bar(stat = "identity") +
    geom_point() +
    facet_wrap(~term, ncol = 1, scales = "free_y")

# Show one example plot
fit_mean |>
    dplyr::filter(group %in% "E-Cigarettes") |>
    dplyr::pull(plot) |>
    pluck(1)
```


## Linear Model Fitting

```{r}
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

```{r}
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

```{r}
glance(fit_one)
tidy(fit_one)
```

```{r}
extract_coefficients(fit_mixed)
```

## Post-hoc Estimated Marginal Means and Comparisons

```{r}
cp_posthoc_slopes(fit_mixed)
cp_posthoc_intercepts(fit_mixed)
```

## See Also

- `vignette("beezdemand")` -- Getting started with beezdemand
- `vignette("model-selection")` -- Choosing the right model class for your data
- `vignette("group-comparisons")` -- Extra sum-of-squares F-test for group comparisons
- `vignette("mixed-demand")` -- Mixed-effects nonlinear demand models
- `vignette("mixed-demand-advanced")` -- Advanced mixed-effects topics
- `vignette("hurdle-demand-models")` -- Two-part hurdle demand models
- `vignette("migration-guide")` -- Migrating from `FitCurves()`