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:
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.
glimpse(etm)
#> Rows: 240
#> Columns: 5
#> $ id <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
#> $ x <dbl> 2, 4, 8, 16, 32, 64, 2, 4, 8, 16, 32, 64, 2, 4, 8, 16, 32, 64, …
#> $ y <dbl> 0, 0, 0, 0, 0, 0, 1, 2, 2, 2, 2, 2, 3, 5, 5, 16, 17, 13, 0, 0, …
#> $ target <chr> "alt", "alt", "alt", "alt", "alt", "alt", "alt", "alt", "alt", …
#> $ group <chr> "Cigarettes", "Cigarettes", "Cigarettes", "Cigarettes", "Cigare…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")).
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).
# Load the cross-price example dataset
data("cp", package = "beezdemand")
# Examine structure
glimpse(cp)
#> Rows: 48
#> Columns: 5
#> $ id <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
#> $ x <dbl> 0.00, 0.01, 0.03, 0.06, 0.13, 0.25, 0.50, 1.00, 2.00, 4.00, 8.0…
#> $ y <dbl> 24.3430657, 22.8832117, 22.6058394, 21.0291971, 19.7810219, 15.…
#> $ target <chr> "alone", "alone", "alone", "alone", "alone", "alone", "alone", …
#> $ group <chr> "cigarettes", "cigarettes", "cigarettes", "cigarettes", "cigare…
# View conditions
table(cp$target)
#>
#> alone alt own
#> 16 16 16The cp dataset contains:
Note that x represents the target commodity
price throughout all conditions.
First, we fit a standard demand curve to the target commodity when the alternative is absent:
# 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
#>
#> Fixed-Effect Demand Model
#> ==========================
#>
#> Call:
#> fit_demand_fixed(data = alone_data, equation = "koff", k = 2)
#>
#> Equation: koff
#> k: fixed (2)
#> Subjects: 1 ( 1 converged, 0 failed)
#>
#> Use summary() for parameter summaries, tidy() for tidy output.Next, fit the same demand model to target consumption when the alternative is present:
# 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
#>
#> Fixed-Effect Demand Model
#> ==========================
#>
#> Call:
#> fit_demand_fixed(data = own_data, equation = "koff", k = 2)
#>
#> Equation: koff
#> k: fixed (2)
#> Subjects: 1 ( 1 converged, 0 failed)
#>
#> Use summary() for parameter summaries, tidy() for tidy output.Finally, fit the cross-price model to alternative consumption as a function of target price:
# 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)
#> Cross-Price Demand Model Summary
#> ================================
#>
#> Model Specification:
#> Equation type: exponentiated
#> Functional form: y ~ (10^log10_qalone) * 10^(I * exp(-(10^log10_beta) * x))
#> Fitting method: nls_multstart
#> Method details: Multiple starting values optimization with nls.multstart
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> log10_qalone 1.1720228 0.0025378 461.8250 < 2.2e-16 ***
#> I -0.3712609 0.0066765 -55.6075 < 2.2e-16 ***
#> log10_beta -0.1271127 0.0217331 -5.8488 5.696e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Confidence Intervals:
#> 2.5 % 97.5 %
#> log10_qalone 1.1665 1.17751
#> I -0.3857 -0.35684
#> log10_beta -0.1741 -0.08016
#>
#> Fit Statistics:
#> R-squared: 0.9973
#> AIC: 1.3
#> BIC: 4.39
#>
#> Parameter Interpretation (natural scale):
#> qalone (Q_alone): 14.86 - consumption at zero alternative price
#> I: -0.3713 - interaction parameter (substitution direction)
#> beta: 0.7463 - sensitivity parameter (sensitivity of relation to price)
#>
#> Optimizer parameters (log10 scale):
#> log10_qalone: 1.172
#> log10_beta: -0.1271fit_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.
# 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
#> Condition Q0_or_Qalone Alpha_or_I
#> 1 Alone (Target) 23.54818 0.01406589
#> 2 Own (Target) 21.19336 0.01562876
#> 3 Alt (Cross-Price) NA -0.37126092Interpretation:
Q0 between alone and own conditions shows how
the presence of an alternative affects baseline consumption of the
target commodityI parameter from the cross-price model indicates
whether the products are substitutes (I < 0) or complements (I >
0)# 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()etm |>
dplyr::filter(group %in% "E-Cigarettes" & id %in% 1)
#> # A tibble: 6 × 5
#> id x y target group
#> <dbl> <dbl> <dbl> <chr> <chr>
#> 1 1 2 3 alt E-Cigarettes
#> 2 1 4 5 alt E-Cigarettes
#> 3 1 8 5 alt E-Cigarettes
#> 4 1 16 16 alt E-Cigarettes
#> 5 1 32 17 alt E-Cigarettes
#> 6 1 64 13 alt E-Cigarettes
unsys_one <- etm |>
filter(group %in% "E-Cigarettes" & id %in% 1) |>
check_systematic_cp()
unsys_one$results
#> # A tibble: 1 × 15
#> id type trend_stat trend_threshold trend_direction trend_pass bounce_stat
#> <chr> <chr> <dbl> <dbl> <chr> <lgl> <dbl>
#> 1 1 cp NA 0.025 up NA 0.2
#> # ℹ 8 more variables: bounce_threshold <dbl>, bounce_direction <chr>,
#> # bounce_pass <lgl>, reversals <int>, reversals_pass <lgl>, returns <int>,
#> # n_positive <int>, systematic <lgl>
unsys_one_lnic <- lnic |>
filter(
target == "adjusting",
id == "R_00Q12ahGPKuESBT"
) |>
check_systematic_cp()
unsys_one_lnic$results
#> # A tibble: 1 × 15
#> id type trend_stat trend_threshold trend_direction trend_pass bounce_stat
#> <chr> <chr> <dbl> <dbl> <chr> <lgl> <dbl>
#> 1 R_00Q… cp NA 0.025 down NA 0
#> # ℹ 8 more variables: bounce_threshold <dbl>, bounce_direction <chr>,
#> # bounce_pass <lgl>, reversals <int>, reversals_pass <lgl>, returns <int>,
#> # n_positive <int>, systematic <lgl>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)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)"
)| group | n_subjects | pct_systematic |
|---|---|---|
| Cigarettes | 10 | 90 |
| Combustibles | 10 | 70 |
| E-Cigarettes | 10 | 80 |
| Non-Combustibles | 10 | 90 |
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)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)"
)| condition | n_subjects | pct_systematic |
|---|---|---|
| 100% | 67 | 97.0 |
| 2% | 57 | 91.2 |
| 2% NegFrame | 59 | 91.5 |
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)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)"
)| target | n_subjects | pct_systematic |
|---|---|---|
| cannabisFix | 99 | 67.7 |
| cigarettesFix | 99 | 77.8 |
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)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)"
)| target | n_subjects | pct_systematic |
|---|---|---|
| ECig | 47 | 83.0 |
| FixCig | 47 | 89.4 |
fit_one <- etm |>
dplyr::filter(group %in% "E-Cigarettes" & id %in% 1) |>
fit_cp_nls(
equation = "exponentiated",
return_all = TRUE
)
summary(fit_one)
#> Cross-Price Demand Model Summary
#> ================================
#>
#> Model Specification:
#> Equation type: exponentiated
#> Functional form: y ~ (10^log10_qalone) * 10^(I * exp(-(10^log10_beta) * x))
#> Fitting method: nls_multstart
#> Method details: Multiple starting values optimization with nls.multstart
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> log10_qalone 1.194135 0.060311 19.7995 0.0002815 ***
#> I -1.268376 0.837302 -1.5148 0.2270524
#> log10_beta -0.737728 0.265129 -2.7825 0.0688459 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Confidence Intervals:
#> 2.5 % 97.5 %
#> log10_qalone 1.002 1.386
#> I -3.933 1.396
#> log10_beta -1.581 0.106
#>
#> Fit Statistics:
#> R-squared: 0.8597
#> AIC: 34.06
#> BIC: 33.23
#>
#> Parameter Interpretation (natural scale):
#> qalone (Q_alone): 15.64 - consumption at zero alternative price
#> I: -1.268 - interaction parameter (substitution direction)
#> beta: 0.1829 - sensitivity parameter (sensitivity of relation to price)
#>
#> Optimizer parameters (log10 scale):
#> log10_qalone: 1.194
#> log10_beta: -0.7377
plot(fit_one, x_trans = "log10")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)
)# 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)"
)| id | group | term | estimate | std.error |
|---|---|---|---|---|
| 1 | Cigarettes | log10_qalone | -4.515900e+01 | 2.267248e+03 |
| 1 | Cigarettes | I | -2.979000e+00 | 2.266886e+03 |
| 1 | Cigarettes | log10_beta | -3.111000e+00 | 3.393980e+02 |
| 1 | Combustibles | log10_qalone | 3.010000e-01 | 0.000000e+00 |
| 1 | Combustibles | I | -5.212610e+02 | 3.023230e+02 |
| 1 | Combustibles | log10_beta | 5.720000e-01 | 3.400000e-02 |
| 1 | E-Cigarettes | log10_qalone | 1.194000e+00 | 6.000000e-02 |
| 1 | E-Cigarettes | I | -1.268000e+00 | 8.370000e-01 |
| 1 | E-Cigarettes | log10_beta | -7.380000e-01 | 2.650000e-01 |
| 1 | Non-Combustibles | log10_qalone | -4.521400e+01 | 4.216897e+04 |
| 1 | Non-Combustibles | I | -2.963000e+00 | 4.216867e+04 |
| 1 | Non-Combustibles | log10_beta | -3.742000e+00 | 6.217747e+03 |
| 2 | Cigarettes | log10_qalone | -3.078160e+02 | 1.915342e+06 |
| 2 | Cigarettes | I | 3.082850e+02 | 1.915342e+06 |
| 2 | Cigarettes | log10_beta | -4.310000e+00 | 2.700805e+03 |
| 2 | Combustibles | log10_qalone | -4.830000e-01 | 2.381000e+00 |
| 2 | Combustibles | I | 9.950000e-01 | 2.271000e+00 |
| 2 | Combustibles | log10_beta | -1.527000e+00 | 1.663000e+00 |
| 2 | E-Cigarettes | log10_qalone | -3.077080e+02 | 1.133752e+06 |
| 2 | E-Cigarettes | I | 3.082790e+02 | 1.133752e+06 |
| 2 | E-Cigarettes | log10_beta | -4.397000e+00 | 1.598611e+03 |
| 2 | Non-Combustibles | log10_qalone | -3.080510e+02 | 5.758916e+06 |
| 2 | Non-Combustibles | I | 3.082640e+02 | 5.758916e+06 |
| 2 | Non-Combustibles | log10_beta | -4.831000e+00 | 8.117098e+03 |
| 3 | Cigarettes | log10_qalone | 1.000000e+00 | 0.000000e+00 |
| 3 | Cigarettes | I | -4.017660e+02 | 2.325640e+02 |
| 3 | Cigarettes | log10_beta | -3.290000e-01 | 3.400000e-02 |
| 3 | Combustibles | log10_qalone | -4.481500e+01 | 3.020400e+03 |
| 3 | Combustibles | I | -1.428000e+00 | 3.020061e+03 |
| 3 | Combustibles | log10_beta | -3.171000e+00 | 9.393800e+02 |
| 3 | E-Cigarettes | log10_qalone | -4.501600e+01 | 5.324700e+01 |
| 3 | E-Cigarettes | I | -2.466000e+00 | 5.244800e+01 |
| 3 | E-Cigarettes | log10_beta | -2.244000e+00 | 1.166900e+01 |
| 3 | Non-Combustibles | log10_qalone | -4.510200e+01 | 2.997900e+01 |
| 3 | Non-Combustibles | I | -2.269000e+00 | 2.892100e+01 |
| 3 | Non-Combustibles | log10_beta | -2.094000e+00 | 7.884000e+00 |
| 4 | Cigarettes | log10_qalone | -4.395200e+01 | 5.640785e+03 |
| 4 | Cigarettes | I | -2.957000e+00 | 5.640440e+03 |
| 4 | Cigarettes | log10_beta | -3.308000e+00 | 8.422170e+02 |
| 4 | Combustibles | log10_qalone | 1.447000e+00 | 5.400000e-02 |
| 4 | Combustibles | I | -3.174723e+03 | 1.313616e+06 |
| 4 | Combustibles | log10_beta | 1.200000e-02 | 2.183300e+01 |
| 4 | E-Cigarettes | log10_qalone | -4.491600e+01 | 7.298490e+02 |
| 4 | E-Cigarettes | I | -2.872000e+00 | 7.294550e+02 |
| 4 | E-Cigarettes | log10_beta | -2.865000e+00 | 1.157390e+02 |
| 4 | Non-Combustibles | log10_qalone | -4.533800e+01 | 1.974971e+04 |
| 4 | Non-Combustibles | I | -1.332000e+00 | 1.974939e+04 |
| 4 | Non-Combustibles | log10_beta | -3.577000e+00 | 6.498072e+03 |
| 5 | Cigarettes | log10_qalone | 0.000000e+00 | 0.000000e+00 |
| 5 | Cigarettes | I | 0.000000e+00 | 0.000000e+00 |
| 5 | Cigarettes | log10_beta | -7.120000e+00 | 6.800000e-02 |
| 5 | Combustibles | log10_qalone | 1.448000e+00 | 0.000000e+00 |
| 5 | Combustibles | I | -6.810000e+00 | 1.150000e-01 |
| 5 | Combustibles | log10_beta | 1.800000e-02 | 3.000000e-03 |
| 5 | E-Cigarettes | log10_qalone | 9.030000e-01 | 0.000000e+00 |
| 5 | E-Cigarettes | I | -1.486004e+11 | 4.754863e+10 |
| 5 | E-Cigarettes | log10_beta | 1.070000e+00 | 6.000000e-03 |
| 5 | Non-Combustibles | log10_qalone | 1.724000e+00 | 0.000000e+00 |
| 5 | Non-Combustibles | I | -8.603600e+02 | 4.971870e+02 |
| 5 | Non-Combustibles | log10_beta | 7.000000e-02 | 2.700000e-02 |
| 6 | Cigarettes | log10_qalone | -4.504800e+01 | 8.240780e+02 |
| 6 | Cigarettes | I | -2.277000e+00 | 8.237050e+02 |
| 6 | Cigarettes | log10_beta | -2.892000e+00 | 1.642750e+02 |
| 6 | Combustibles | log10_qalone | -4.540600e+01 | 2.709070e+02 |
| 6 | Combustibles | I | -2.843000e+00 | 2.704450e+02 |
| 6 | Combustibles | log10_beta | -2.643000e+00 | 4.490300e+01 |
| 6 | E-Cigarettes | log10_qalone | 3.160000e-01 | 4.600000e-02 |
| 6 | E-Cigarettes | I | -4.890000e-01 | 1.680000e-01 |
| 6 | E-Cigarettes | log10_beta | -9.410000e-01 | 2.580000e-01 |
| 6 | Non-Combustibles | log10_qalone | -4.525000e+01 | 6.354410e+02 |
| 6 | Non-Combustibles | I | -2.619000e+00 | 6.350490e+02 |
| 6 | Non-Combustibles | log10_beta | -2.835000e+00 | 1.108870e+02 |
| 7 | Cigarettes | log10_qalone | 3.010000e-01 | 0.000000e+00 |
| 7 | Cigarettes | I | -2.154420e+02 | 1.342600e+01 |
| 7 | Cigarettes | log10_beta | -6.870000e-01 | 4.000000e-03 |
| 7 | Combustibles | log10_qalone | 1.082000e+00 | 5.900000e-01 |
| 7 | Combustibles | I | -5.340000e-01 | 4.990000e-01 |
| 7 | Combustibles | log10_beta | -1.639000e+00 | 1.047000e+00 |
| 7 | E-Cigarettes | log10_qalone | -4.466600e+01 | 4.022285e+03 |
| 7 | E-Cigarettes | I | -2.971000e+00 | 4.021935e+03 |
| 7 | E-Cigarettes | log10_beta | -3.235000e+00 | 5.996800e+02 |
| 7 | Non-Combustibles | log10_qalone | -4.536000e+01 | 9.208600e+03 |
| 7 | Non-Combustibles | I | -2.277000e+00 | 9.208267e+03 |
| 7 | Non-Combustibles | log10_beta | -3.413000e+00 | 1.779549e+03 |
| 8 | Cigarettes | log10_qalone | -4.510600e+01 | 4.086164e+04 |
| 8 | Cigarettes | I | -2.664000e+00 | 4.086134e+04 |
| 8 | Cigarettes | log10_beta | -3.735000e+00 | 6.701295e+03 |
| 8 | Combustibles | log10_qalone | 3.082040e+02 | 3.013474e+05 |
| 8 | Combustibles | I | -3.083770e+02 | 3.013469e+05 |
| 8 | Combustibles | log10_beta | -4.245000e+00 | 4.253780e+02 |
| 8 | E-Cigarettes | log10_qalone | 4.770000e-01 | 8.400000e-02 |
| 8 | E-Cigarettes | I | -8.805930e+02 | 4.930348e+05 |
| 8 | E-Cigarettes | log10_beta | -2.700000e-02 | 3.233800e+01 |
| 8 | Non-Combustibles | log10_qalone | 4.810000e-01 | 3.500000e-02 |
| 8 | Non-Combustibles | I | -2.251000e+00 | 7.430000e-01 |
| 8 | Non-Combustibles | log10_beta | -1.077000e+00 | 9.900000e-02 |
| 9 | Cigarettes | log10_qalone | -4.524400e+01 | 2.270866e+04 |
| 9 | Cigarettes | I | -2.910000e+00 | 2.270833e+04 |
| 9 | Cigarettes | log10_beta | -3.608000e+00 | 3.417431e+03 |
| 9 | Combustibles | log10_qalone | 3.010000e-01 | 0.000000e+00 |
| 9 | Combustibles | I | -1.630526e+11 | 5.225423e+10 |
| 9 | Combustibles | log10_beta | 7.700000e-01 | 6.000000e-03 |
| 9 | E-Cigarettes | log10_qalone | 6.520000e-01 | 1.170000e-01 |
| 9 | E-Cigarettes | I | -4.460000e-01 | 1.180000e-01 |
| 9 | E-Cigarettes | log10_beta | -1.387000e+00 | 3.800000e-01 |
| 9 | Non-Combustibles | log10_qalone | -4.419500e+01 | 4.062960e+02 |
| 9 | Non-Combustibles | I | -2.695000e+00 | 4.058770e+02 |
| 9 | Non-Combustibles | log10_beta | -2.736000e+00 | 6.985300e+01 |
| 10 | Cigarettes | log10_qalone | -4.541400e+01 | 7.800100e+01 |
| 10 | Cigarettes | I | -2.402000e+00 | 7.738400e+01 |
| 10 | Cigarettes | log10_beta | -2.352000e+00 | 1.664800e+01 |
| 10 | Combustibles | log10_qalone | 1.279000e+00 | 0.000000e+00 |
| 10 | Combustibles | I | -6.525461e+03 | 3.788643e+03 |
| 10 | Combustibles | log10_beta | 6.440000e-01 | 2.900000e-02 |
| 10 | E-Cigarettes | log10_qalone | 0.000000e+00 | 0.000000e+00 |
| 10 | E-Cigarettes | I | 0.000000e+00 | 0.000000e+00 |
| 10 | E-Cigarettes | log10_beta | -4.892000e+00 | 9.500000e-02 |
| 10 | Non-Combustibles | log10_qalone | 1.439000e+00 | 1.400000e-02 |
| 10 | Non-Combustibles | I | -8.492300e+01 | 1.419300e+02 |
| 10 | Non-Combustibles | log10_beta | 3.090000e-01 | 1.500000e-01 |
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)
)# 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"
)| group | term | estimate | std.error |
|---|---|---|---|
| Cigarettes | log10_qalone | 0.098 | 0.205 |
| Cigarettes | I | -0.868 | 1.251 |
| Cigarettes | log10_beta | -0.938 | 0.827 |
| Combustibles | log10_qalone | 0.969 | 0.098 |
| Combustibles | I | -0.661 | 0.683 |
| Combustibles | log10_beta | -0.743 | 0.579 |
| E-Cigarettes | log10_qalone | 0.540 | 0.105 |
| E-Cigarettes | I | -0.661 | 0.735 |
| E-Cigarettes | log10_beta | -0.742 | 0.622 |
| Non-Combustibles | log10_qalone | 0.924 | 0.134 |
| Non-Combustibles | I | -4.148 | 19.320 |
| Non-Combustibles | log10_beta | -0.271 | 0.904 |
# Show one plot example
fit_pooled |>
dplyr::filter(group == "E-Cigarettes") |>
dplyr::pull(plot) |>
pluck(1)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)
)# 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"
)| group | term | estimate | std.error |
|---|---|---|---|
| Cigarettes | log10_qalone | 0.098 | 0.074 |
| Cigarettes | I | -0.868 | 0.453 |
| Cigarettes | log10_beta | -0.938 | 0.299 |
| Combustibles | log10_qalone | 0.969 | 0.031 |
| Combustibles | I | -0.661 | 0.218 |
| Combustibles | log10_beta | -0.743 | 0.185 |
| E-Cigarettes | log10_qalone | 0.540 | 0.052 |
| E-Cigarettes | I | -0.661 | 0.367 |
| E-Cigarettes | log10_beta | -0.742 | 0.310 |
| Non-Combustibles | log10_qalone | 0.924 | 0.007 |
| Non-Combustibles | I | -4.148 | 0.992 |
| Non-Combustibles | log10_beta | -0.271 | 0.046 |
# 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)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)
#> Linear Cross-Price Demand Model Summary
#> =======================================
#>
#> Formula: y ~ log10(x)
#> Method: lm
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.13333 3.55773 0.0375 0.97190
#> log10(x) 9.20649 3.03472 3.0337 0.03864 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-squared: 0.697 Adjusted R-squared: 0.6213
plot(fit_one_linear, x_trans = "log10")fit_mixed <- fit_cp_linear(
etm,
type = "mixed",
log10x = TRUE,
group_effects = "interaction",
random_slope = FALSE,
return_all = TRUE
)
summary(fit_mixed)
#> Mixed-Effects Linear Cross-Price Demand Model Summary
#> ====================================================
#>
#> Formula: y ~ log10(x) * group + (1 | id)
#> Method: lmer
#>
#> Fixed Effects:
#> Estimate Std. Error t value
#> (Intercept) -0.10000 2.66378 -0.0375
#> log10(x) 0.80675 1.85305 0.4354
#> groupCombustibles 2.09333 3.07225 0.6814
#> groupE-Cigarettes 0.98667 3.07225 0.3212
#> groupNon-Combustibles 0.14667 3.07225 0.0477
#> log10(x):groupCombustibles 3.83445 2.62060 1.4632
#> log10(x):groupE-Cigarettes 0.78777 2.62060 0.3006
#> log10(x):groupNon-Combustibles 4.76459 2.62060 1.8181
#>
#> Random Effects:
#> Group Term Variance Std.Dev NA
#> id (Intercept) <NA> 23.76351 4.874783
#> Residual <NA> <NA> 54.45409 7.379301
#>
#> Model Fit:
#> R2 (marginal): 0.1121 [Fixed effects only]
#> R2 (conditional): 0.3819 [Fixed + random effects]
#> AIC: 1655
#> BIC: 1690
#>
#> Note: R2 values for mixed models are approximate.
# plot fixed effects only
plot(fit_mixed, x_trans = "log10", pred_type = "fixed")glance(fit_one)
#> # A tibble: 1 × 6
#> r.squared aic bic equation method transform
#> <dbl> <dbl> <dbl> <chr> <chr> <chr>
#> 1 0.860 34.1 33.2 exponentiated nls_multstart none
tidy(fit_one)
#> # A tibble: 3 × 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 log10_qalone 1.19 0.0603 19.8 0.000282
#> 2 I -1.27 0.837 -1.51 0.227
#> 3 log10_beta -0.738 0.265 -2.78 0.0688extract_coefficients(fit_mixed)
#> $fixed
#> (Intercept) log10(x)
#> -0.1000000 0.8067540
#> groupCombustibles groupE-Cigarettes
#> 2.0933333 0.9866667
#> groupNon-Combustibles log10(x):groupCombustibles
#> 0.1466667 3.8344541
#> log10(x):groupE-Cigarettes log10(x):groupNon-Combustibles
#> 0.7877715 4.7645940
#>
#> $random
#> $id
#> (Intercept)
#> 1 -1.0155373
#> 2 -2.0805203
#> 3 -2.6890820
#> 4 0.1255158
#> 5 11.0796263
#> 6 -3.3356788
#> 7 -2.3087309
#> 8 -2.4608713
#> 9 -2.7651522
#> 10 5.4504307
#>
#> with conditional variances for "id"
#>
#> $combined
#> $id
#> (Intercept) log10(x) groupCombustibles groupE-Cigarettes
#> 1 -1.11553732 0.806754 2.093333 0.9866667
#> 2 -2.18052028 0.806754 2.093333 0.9866667
#> 3 -2.78908198 0.806754 2.093333 0.9866667
#> 4 0.02551585 0.806754 2.093333 0.9866667
#> 5 10.97962630 0.806754 2.093333 0.9866667
#> 6 -3.43567877 0.806754 2.093333 0.9866667
#> 7 -2.40873092 0.806754 2.093333 0.9866667
#> 8 -2.56087134 0.806754 2.093333 0.9866667
#> 9 -2.86515219 0.806754 2.093333 0.9866667
#> 10 5.35043065 0.806754 2.093333 0.9866667
#> groupNon-Combustibles log10(x):groupCombustibles log10(x):groupE-Cigarettes
#> 1 0.1466667 3.834454 0.7877715
#> 2 0.1466667 3.834454 0.7877715
#> 3 0.1466667 3.834454 0.7877715
#> 4 0.1466667 3.834454 0.7877715
#> 5 0.1466667 3.834454 0.7877715
#> 6 0.1466667 3.834454 0.7877715
#> 7 0.1466667 3.834454 0.7877715
#> 8 0.1466667 3.834454 0.7877715
#> 9 0.1466667 3.834454 0.7877715
#> 10 0.1466667 3.834454 0.7877715
#> log10(x):groupNon-Combustibles
#> 1 4.764594
#> 2 4.764594
#> 3 4.764594
#> 4 4.764594
#> 5 4.764594
#> 6 4.764594
#> 7 4.764594
#> 8 4.764594
#> 9 4.764594
#> 10 4.764594
#>
#> attr(,"class")
#> [1] "coef.mer"cp_posthoc_slopes(fit_mixed)
#> Slope Estimates and Comparisons
#> ===============================
#>
#> Estimated Marginal Means:
#> group x.trend SE df lower.CL upper.CL
#> Cigarettes 0.01665965 0.03826583 223 -0.05874926 0.09206855
#> Combustibles 0.09584197 0.03826583 223 0.02043307 0.17125088
#> E-Cigarettes 0.03292730 0.03826583 223 -0.04248160 0.10833621
#> Non-Combustibles 0.11504957 0.03826583 223 0.03964066 0.19045847
#>
#> Degrees-of-freedom method: kenward-roger
#> Confidence level used: 0.95
#>
#> Significant interaction: No
#>
#> No significant interaction detected (alpha = 0.05 ). Pairwise slope comparisons not performed.
#> P-value adjustment method: tukey
cp_posthoc_intercepts(fit_mixed)
#> Intercept Estimates and Comparisons
#> ===================================
#>
#> Estimated Marginal Means:
#> group emmean SE df lower.CL upper.CL
#> Cigarettes -0.1000000 2.663776 59.69 -5.428913 5.228913
#> Combustibles 1.9933333 2.663776 59.69 -3.335580 7.322246
#> E-Cigarettes 0.8866667 2.663776 59.69 -4.442246 6.215580
#> Non-Combustibles 0.0466667 2.663776 59.69 -5.282246 5.375580
#>
#> Degrees-of-freedom method: kenward-roger
#> Confidence level used: 0.95
#>
#> Significant interaction: No
#>
#> No significant interaction detected (alpha = 0.05 ). Pairwise intercept comparisons not performed.
#> P-value adjustment method: tukeyvignette("beezdemand") – Getting started with
beezdemandvignette("model-selection") – Choosing the right model
class for your datavignette("group-comparisons") – Extra sum-of-squares
F-test for group comparisonsvignette("mixed-demand") – Mixed-effects nonlinear
demand modelsvignette("mixed-demand-advanced") – Advanced
mixed-effects topicsvignette("hurdle-demand-models") – Two-part hurdle
demand modelsvignette("migration-guide") – Migrating from
FitCurves()