---
title: "Advanced Mixed-Effects Demand Modeling"
author: "Brent Kaplan"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Advanced Mixed-Effects Demand Modeling}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",   # R console style comments
  dev = "png",
  dpi = 144,
  fig.width = 7,
  fig.height = 5,
  warning = FALSE,  # Suppress warnings for cleaner output
  message = FALSE,  # Suppress messages for cleaner output
  cache = TRUE,     # Cache for faster rebuilds
  cache.path = "mixed-demand-advanced_cache/",
  autodep = TRUE
)

library(beezdemand)
library(dplyr)
library(ggplot2)

data("apt", package = "beezdemand", envir = environment())
data("ko", package = "beezdemand", envir = environment())

# Set `BEEZDEMAND_VIGNETTE_MODE=full` to run the slower / less stable examples.
render_fast <- tolower(Sys.getenv("BEEZDEMAND_VIGNETTE_MODE", "fast")) != "full"

cache_key_object <- function(x) {
  tmp <- tempfile(fileext = ".rds")
  saveRDS(x, tmp)
  on.exit(unlink(tmp), add = TRUE)
  unname(tools::md5sum(tmp))
}

# Invalidate cache when the package data changes.
knitr::opts_chunk$set(
  cache.extra = list(
    beezdemand_version = as.character(utils::packageVersion("beezdemand")),
    apt = cache_key_object(apt),
    ko = cache_key_object(ko)
  )
)

# Optional fits used later (kept NULL in fast renders).
fit_two_factors_add <- NULL
fit_two_factors_int <- NULL
fit_simplified_example <- NULL
fit_collapsed_same <- NULL
fit_collapsed_asymmetric <- NULL
fit_drug_slope_random <- NULL
fit_age_add <- NULL
fit_age_rhs <- NULL
```

# Introduction

This vignette covers advanced topics for mixed-effects nonlinear demand
modeling with `beezdemand`. It assumes you are already familiar with the basics
covered in `vignette("mixed-demand")`, including fitting models with
`fit_demand_mixed()`, inspecting results with `tidy()` / `glance()` /
`augment()`, and basic plotting.

Topics covered here include:

- **Multi-factor models** (additive and interaction)
- **Collapsing factor levels** (asymmetric specifications for Q0 and alpha)
- **Estimated marginal means** and **pairwise comparisons**
- **Multi-factor visualization** with faceting
- **Complex random effects structures**
- **Continuous covariates** and `fixed_rhs`
- **Trends** with `emtrends`

```{r shared-setup}
quick_nlme_control <- nlme::nlmeControl(
  msMaxIter = 100,
  niterEM = 20,
  maxIter = 100, # Low iterations for speed
  pnlsTol = 0.1,
  tolerance = 1e-4, # Looser tolerance
  opt = "nlminb",
  msVerbose = FALSE
)

# Prepare data subsets used throughout
ko_alf <- ko[ko$drug == "Alfentanil", ]

# Fit the one-factor model used as a baseline/fallback in later sections
fit_one_factor_dose <- fit_demand_mixed(
  data = ko_alf,
  y_var = "y_ll4",
  x_var = "x",
  id_var = "monkey",
  factors = "dose",
  equation_form = "zben",
  nlme_control = quick_nlme_control,
  start_value_method = "heuristic"
)
```

## Model with Two Factors (Additive)

Here, $Q_{0}$ and $\alpha$ vary by drug and dose additively.
For this, we need more data than just ko_alf, so we use `ko`.
Note: With complex models and small sample sizes, convergence can be challenging.
The start_value_method = "pooled_nls" is often more robust for complex models.

```{r two_factor_add_note, echo=FALSE, results='asis'}
if (render_fast) {
  cat(
    "*This example is computationally intensive and is not run during standard vignette building. ",
    "To run it, set `BEEZDEMAND_VIGNETTE_MODE=full` as an environment variable before building. ",
    "When successful, the output shows fixed-effect estimates for each drug and dose combination ",
    "on Q0 and alpha, plus random-effects variance components.*\n"
  )
}
```

```{r two_factor_add, eval = !render_fast, include = !render_fast}
# This model might be slow or fail to converge with the small 'ko' and quick_nlme_control
# In a real vignette, use a dataset known to converge or increase iterations.
# For demonstration, we'll proceed.
fit_two_factors_add <- try(
  fit_demand_mixed(
    data = ko,
    y_var = "y_ll4",
    x_var = "x",
    id_var = "monkey",
    factors = c("drug", "dose"),
    factor_interaction = FALSE, # Additive effects
    equation_form = "zben",
    nlme_control = quick_nlme_control,
    start_value_method = "pooled_nls"
  ),
  silent = TRUE
)

if (!inherits(fit_two_factors_add, "try-error")) {
  print(fit_two_factors_add)
} else {
  cat(
    "Two-factor additive model failed to converge with current settings for the vignette.\n"
  )
  fit_two_factors_add <- NULL # Ensure it's NULL for later checks
}
```

### Model with Two Factors and Interaction

This allows the effect of dose to be different for each drug (and vice-versa).

```{r two_factor_int_note, echo=FALSE, results='asis'}
if (render_fast) {
  cat(
    "*This example is computationally intensive and is not run during standard vignette building. ",
    "The interaction model allows the effect of dose to differ by drug (and vice versa), ",
    "producing drug:dose interaction terms in the fixed-effects table.*\n"
  )
}
```

```{r two_factor_int, eval = !render_fast, include = !render_fast}
fit_two_factors_int <- try(
  fit_demand_mixed(
    data = ko,
    y_var = "y_ll4",
    x_var = "x",
    id_var = "monkey",
    factors = c("drug", "dose"),
    factor_interaction = TRUE, # Interaction effect
    equation_form = "zben",
    nlme_control = quick_nlme_control,
    start_value_method = "pooled_nls"
  ),
  silent = TRUE
)

if (!inherits(fit_two_factors_int, "try-error")) {
  print(fit_two_factors_int)
} else {
  cat(
    "Two-factor interaction model failed to converge with current settings for the vignette.\n"
  )
  fit_two_factors_int <- NULL
}
```

### Using Different equation_form and y_var

The simplified equation form expects y_var to be raw consumption.

```{r simplified_equation_note, echo=FALSE, results='asis'}
if (render_fast) {
  cat(
    "*This example uses `equation_form = \"simplified\"` with raw consumption values (no LL4 ",
    "transformation needed). The simplified equation handles zeros natively. When this example ",
    "runs, the output resembles the `zben` model output but with `y` as the dependent variable.*\n"
  )
}
```

```{r simplified_equation, eval = !render_fast, include = !render_fast}
better_nlme_control <- nlme::nlmeControl(
  msMaxIter = 100, # Max iterations for nlmeStruct optimization
  niterEM = 50, # Max EM iterations
  maxIter = 100, # Max iterations for the optim call for fixed effects
  pnlsTol = 1e-4, # Tolerance for PNLS convergence
  tolerance = 1e-7, # Tolerance for outer optimization loop
  opt = "nlminb",
  msVerbose = TRUE # Set to TRUE to see optimization steps
)

fit_simplified_example <- try(
  fit_demand_mixed(
    data = ko_alf, # Using simpler Alfentanil data
    y_var = "y", # Raw consumption
    x_var = "x",
    id_var = "monkey",
    factors = "dose",
    equation_form = "simplified",
    nlme_control = quick_nlme_control,
    start_value_method = "pooled_nls"
  ),
  silent = TRUE
)

if (!inherits(fit_simplified_example, "try-error")) {
  print(fit_simplified_example)
} else {
  cat(
    "Simplified equation model failed to converge with current settings for the vignette.\n"
  )
  fit_simplified_example <- NULL
}
```

## Collapsing Factor Levels

Sometimes you might want to group levels of a factor, and you may want different
groupings for $Q_{0}$ versus $\alpha$. The `collapse_levels` argument allows you
to specify separate collapsing schemes for each parameter.

**Structure:**
```r
collapse_levels = list(
  Q0 = list(factor_name = list(new_level = c("old_level1", "old_level2"), ...)),
  alpha = list(factor_name = list(new_level = c("old_level1", ...), ...))
)
```

Either `Q0` or `alpha` (or both) can be omitted to keep original levels for that parameter.

### Example: Same collapsing for both parameters

```{r}
# Ensure levels to collapse are present in ko_alf$dose
# levels(ko_alf$dose) are "0.001", "0.003", "3e-04"
fit_collapsed_same <- try(
  fit_demand_mixed(
    data = ko_alf,
    y_var = "y_ll4",
    x_var = "x",
    id_var = "monkey",
    factors = "dose",
    collapse_levels = list(
      Q0 = list(
        dose = list(low_doses = c("3e-04", "0.001"), high_dose = "0.003")
      ),
      alpha = list(
        dose = list(low_doses = c("3e-04", "0.001"), high_dose = "0.003")
      )
    ),
    equation_form = "zben",
    nlme_control = quick_nlme_control
  ),
  silent = TRUE
)

if (
  !is.null(fit_collapsed_same) &&
    !inherits(fit_collapsed_same, "try-error") &&
    !is.null(fit_collapsed_same$model)
) {
  print(fit_collapsed_same)
  cat("\nQ0 params:", fit_collapsed_same$param_info$num_params_Q0, "\n")
  cat("alpha params:", fit_collapsed_same$param_info$num_params_alpha, "\n")
} else {
  cat("Collapsed levels model failed to converge.\n")
}
```

### Example: Different collapsing for Q0 and alpha

This is particularly useful when you want fine-grained distinctions for one
parameter but not the other. For instance, you might hypothesize that maximum
consumption ($Q_{0}$) varies by dose, but sensitivity ($\alpha$) does not.

```{r}
# Q0: keep 2 collapsed levels (low vs high)
# alpha: collapse all to 1 level (intercept only)
fit_collapsed_asymmetric <- try(
  fit_demand_mixed(
    data = ko_alf,
    y_var = "y_ll4",
    x_var = "x",
    id_var = "monkey",
    factors = "dose",
    collapse_levels = list(
      Q0 = list(
        dose = list(low_doses = c("3e-04", "0.001"), high_dose = "0.003")
      ),
      alpha = list(dose = list(all_doses = c("3e-04", "0.001", "0.003")))
    ),
    equation_form = "zben",
    nlme_control = quick_nlme_control
  ),
  silent = TRUE
)

if (
  !is.null(fit_collapsed_asymmetric) &&
    !inherits(fit_collapsed_asymmetric, "try-error") &&
    !is.null(fit_collapsed_asymmetric$model)
) {
  print(fit_collapsed_asymmetric)
  cat("\nQ0 params:", fit_collapsed_asymmetric$param_info$num_params_Q0, "\n")
  cat(
    "alpha params:",
    fit_collapsed_asymmetric$param_info$num_params_alpha,
    "\n"
  )
  cat(
    "\nQ0 formula:",
    fit_collapsed_asymmetric$formula_details$fixed_effects_formula_str_Q0,
    "\n"
  )
  cat(
    "alpha formula:",
    fit_collapsed_asymmetric$formula_details$fixed_effects_formula_str_alpha,
    "\n"
  )
} else {
  cat("Asymmetric collapsed levels model failed to converge.\n")
}
```

**Note:** When all levels of a factor are collapsed to a single level for a
parameter, that factor is automatically removed from the formula for that
parameter (it contributes only to the intercept). A message will inform you
when this occurs.

### EMMs and Comparisons with Collapsed Factors

When using `collapse_levels`, the `get_demand_param_emms()` and
`get_demand_comparisons()` functions automatically handle the asymmetric
factor structures. EMMs will be computed using the collapsed levels, and
comparisons will only be performed when there are multiple levels to compare.

**Differential Collapsing Output:** When Q0 and alpha have different numbers
of factor levels (e.g., Q0 keeps original `dose` levels while alpha collapses
to fewer groups), the EMM output will include separate factor columns:

- Original factor column (e.g., `dose`) for the uncollapsed parameter (Q0)
- Suffixed column (e.g., `dose_alpha`) for the collapsed parameter (alpha)

This results in a cross-joined table showing all combinations of Q0 levels
and alpha levels. For example, if dose has 5 original levels and alpha
collapses to 2 groups, the EMM table will have 10 rows (5 × 2).

```{r}
# Using the asymmetric collapsed model from above (if it converged)
if (
  !is.null(fit_collapsed_asymmetric) &&
    !inherits(fit_collapsed_asymmetric, "try-error") &&
    !is.null(fit_collapsed_asymmetric$model)
) {
  cat("--- EMMs with collapsed factors ---\n")
  # EMMs will show collapsed Q0 levels but single alpha value
  collapsed_emms <- get_demand_param_emms(
    fit_obj = fit_collapsed_asymmetric,
    factors_in_emm = "dose", # Use original factor name
    include_ev = TRUE
  )
  print(collapsed_emms)

  cat("\n--- Comparisons with collapsed factors ---\n")
  # Q0 will have comparisons (2 levels), alpha will have none (1 level)
  collapsed_comparisons <- get_demand_comparisons(
    fit_obj = fit_collapsed_asymmetric,
    compare_specs = ~dose, # Use original factor name
    params_to_compare = c("Q0", "alpha")
  )

  cat("\nQ0 contrasts (collapsed levels):\n")
  print(collapsed_comparisons$Q0$contrasts_log10)

  cat("\nalpha contrasts (intercept-only, no comparisons):\n")
  print(collapsed_comparisons$alpha$contrasts_log10)
} else {
  cat("Asymmetric collapsed model not available for EMM demonstration.\n")
}
```

## Visualizing Multi-Factor Models

### Plotting a Two-Factor Model with Faceting

Let's use fit_two_factors_add (if it converged: y_ll4 ~ drug + dose). We can
facet by one factor and color by another.

```{r plot_two_factor_faceting, cache = FALSE}
# Assuming fit_two_factors_add was successfully created earlier
# If it failed, this chunk won't produce a plot.
active_two_factor_fit <- if (
  !is.null(fit_two_factors_add) && !is.null(fit_two_factors_add$model)
) {
  fit_two_factors_add
} else {
  # Fallback to a simpler model if the two-factor one failed in the vignette
  if (!is.null(fit_one_factor_dose$model)) fit_one_factor_dose else NULL
}

if (
  !is.null(active_two_factor_fit$model) &&
    !is.null(active_two_factor_fit$param_info$factors) &&
    length(active_two_factor_fit$param_info$factors) >= 1
) {
  # Determine factors for aesthetics based on what's in active_two_factor_fit
  color_factor <- if ("dose" %in% active_two_factor_fit$param_info$factors) {
    "dose"
  } else {
    NULL
  }
  facet_factor_name <- if (
    "drug" %in% active_two_factor_fit$param_info$factors
  ) {
    "drug"
  } else {
    # If only 'dose' is available from fit_one_factor_dose, we can't facet by 'drug'
    # So, maybe facet by dose instead, or don't facet.
    if (
      "dose" %in%
        active_two_factor_fit$param_info$factors &&
        is.null(color_factor)
    ) {
      color_factor <- "dose" # color by dose if not faceting by it
      NULL # No faceting
    } else {
      NULL
    }
  }

  facet_formula_plot <- if (!is.null(facet_factor_name)) {
    stats::as.formula(paste("~", facet_factor_name))
  } else {
    NULL
  }

  plot(
    active_two_factor_fit,
    inv_fun = ll4_inv,
    color_by = color_factor,
    # linetype_by = if("dose" %in% active_two_factor_fit$param_info$factors) "dose" else NULL, # Example
    facet_formula = facet_formula_plot,
    title = "Demand Curves (Population Fit)",
    observed_point_alpha = 0.5,
    ind_line_alpha = .5
  )
} else {
  cat(
    "A suitable two-factor or one-factor model object not available for this plotting example.\n"
  )
}
```

This example attempts to use fit_two_factors_add. If that model didn't converge
(common with minimal iterations for vignette speed), it falls back to
fit_one_factor_dose and adjusts aesthetics. The plot will show population lines,
colored by one factor and faceted by another (if two factors are available).

### Plotting a Model Fit with equation_form = "simplified"

If you fit a model using equation_form = "simplified" (which models raw y), the
inv_fun is typically identity because predictions are already on the natural scale.

```{r plot_simplified_model, cache = FALSE}
# Assuming fit_simplified_example converged earlier
if (
  !is.null(fit_simplified_example) && !is.null(fit_simplified_example$model)
) {
  plot(
    fit_simplified_example,
    inv_fun = identity, # Predictions are already on raw y scale
    color_by = "dose",
    shape_by = "dose",
    title = "Demand Model ('simplified' equation, Raw Y)"
  )
} else {
  cat("fit_simplified_example model object not available for plotting.\n")
}
```

Users can further customize the returned ggplot
object by adding more layers or theme adjustments. For instance, to add custom
axis limits or breaks:

```r
plot_object +
    ggplot2::scale_x_continuous(
        limits = c(0, 1000),
        breaks = c(0, 100, 500, 1000)
    )
```

## Analyzing Estimated Marginal Means (`get_demand_param_emms`)

This function helps interpret how factors affect $Q_{0}$ and $\alpha$, providing estimates on both log10 and natural scales, and optionally Essential Value (EV).

**Note on `collapse_levels`:** When a model is fit using `collapse_levels` with
asymmetric specifications (different collapsing for Q0 and alpha), the EMM functions
automatically use the appropriate collapsed factor names for each parameter. You
still specify the original factor name in `factors_in_emm`, and the output will
show the collapsed levels. If a parameter has only one level (intercept-only),
that parameter's values will be the same across all rows.

```{r}
# We'll use a model with factors.
# If fit_two_factors_add converged, use it. Otherwise, use fit_one_factor_dose.
# For the vignette, let's ensure we use one that is likely available.
# If fit_two_factors_add is NULL (failed to converge in example), this will use fit_one_factor_dose
emm_model_to_use <- if (
  !is.null(fit_two_factors_add) && !is.null(fit_two_factors_add$model)
) {
  fit_two_factors_add
} else if (!is.null(fit_one_factor_dose$model)) {
  fit_one_factor_dose
} else {
  NULL
}

if (!is.null(emm_model_to_use)) {
  cat(
    "--- EMMs for model with factors:",
    paste(emm_model_to_use$param_info$factors, collapse = ", "),
    "---\n"
  )

  factors_for_emms <- emm_model_to_use$param_info$factors

  demand_emms_output <- get_demand_param_emms(
    fit_obj = emm_model_to_use,
    factors_in_emm = factors_for_emms, # Use factors from the model
    include_ev = TRUE
  )
  print(demand_emms_output)

  cat("\n--- EMMs for observed factor combinations only: ---\n")
  # This is useful if the EMM grid includes combinations not in data
  observed_demand_emms <- get_observed_demand_param_emms(
    fit_obj = emm_model_to_use,
    factors_in_emm = factors_for_emms,
    include_ev = TRUE
  )
  print(observed_demand_emms)
} else {
  cat(
    "No suitable model with factors converged for EMM analysis in the vignette.\n"
  )
}
```

## Performing Pairwise Comparisons (`get_demand_comparisons`)

Compare levels of factors for $Q_{0}$ and $\alpha$.

**Note on `collapse_levels`:** When using models fit with asymmetric
`collapse_levels`, comparisons are only performed for parameters that have
multiple levels. If a parameter was collapsed to a single level (intercept-only),
the comparisons for that parameter will be empty.

```{r}
# Using the same emm_model_to_use
if (
  !is.null(emm_model_to_use) && length(emm_model_to_use$param_info$factors) > 0
) {
  factors_present <- emm_model_to_use$param_info$factors

  if ("dose" %in% factors_present) {
    cat(
      "--- Pairwise comparisons for 'dose' (averaging over other factors if any): ---\n"
    )
    comparisons_dose <- get_demand_comparisons(
      fit_obj = emm_model_to_use,
      compare_specs = ~dose,
      contrast_type = "pairwise",
      adjust = "fdr",
      report_ratios = TRUE
    )
    print(comparisons_dose)
  }

  if (all(c("drug", "dose") %in% factors_present)) {
    cat(
      "\n--- Pairwise comparisons for 'drug' within each level of 'dose': ---\n"
    )
    # EMMs calculated over drug*dose, then contrast drug within each dose
    comparisons_drug_by_dose <- get_demand_comparisons(
      fit_obj = emm_model_to_use,
      compare_specs = ~ drug * dose,
      contrast_type = "pairwise",
      contrast_by = "dose", # Compare 'drug' levels, holding 'dose' constant
      adjust = "fdr",
      report_ratios = TRUE
    )
    print(comparisons_drug_by_dose)
  }
} else {
  cat(
    "No suitable model with factors converged for comparisons in the vignette.\n"
  )
}
```

## Advanced Topics

**Note**: This section covers advanced modeling techniques for experienced users.
These advanced topics require deeper understanding of mixed-effects models and
may be computationally intensive.

```{r advanced_topics_note, echo=FALSE, results='asis'}
if (render_fast) {
  cat(
    "**Performance Note**: The code examples in the sections below are computationally intensive ",
    "and are not evaluated during standard vignette building. The code is shown for reference. ",
    "To run these examples, set the environment variable `BEEZDEMAND_VIGNETTE_MODE=full` before ",
    "building vignettes.\n"
  )
}
```

### More Complex Random Effects Structures

This example shows how to specify a more complex random-effects structure (e.g.,
random slopes) for $Q_{0}$ and $\alpha$.

```{r complex_random_effects_fit, eval = !render_fast, include = !render_fast}
robust_nlme_control <- nlme::nlmeControl(
  msMaxIter = 100,
  niterEM = 50,
  maxIter = 100,
  pnlsTol = 1e-4,
  tolerance = 1e-7,
  opt = "nlminb",
  msVerbose = TRUE # Set to FALSE for less output once working
)

fit_drug_slope_random <- try(
  fit_demand_mixed(
    data = ko,
    y_var = "y_ll4",
    x_var = "x",
    id_var = "monkey", # This is your 'monkey' identifier
    factors = c("drug", "dose"),
    factor_interaction = FALSE, # Additive fixed effects as requested
    equation_form = "zben",
    random_effects = Q0 + alpha ~ drug, # Random slopes by drug
    covariance_structure = "pdDiag",
    nlme_control = robust_nlme_control,
    start_value_method = "pooled_nls" # Good for more complex models
  ),
  silent = TRUE
)

if (
  !inherits(fit_drug_slope_random, "try-error") &&
    !is.null(fit_drug_slope_random$model)
) {
  print(fit_drug_slope_random)

  cat("\nFixed Effects:\n")
  print(nlme::fixef(fit_drug_slope_random))

  cat("\nRandom Effects Structure (Standard Deviations and Correlations):\n")
  print(nlme::VarCorr(fit_drug_slope_random$model))

  cat(
    "\nRandom Effect Estimates for each Monkey (deviations from fixed effects):\n"
  )
  print(utils::head(nlme::ranef(fit_drug_slope_random)))

  cat("\nSubject-Specific Coefficients (Fixed + Random, log10 scale):\n")
  print(head(coef(fit_drug_slope_random))) # type="combined" is default

  all_params_emms <- get_demand_param_emms(
    fit_obj = fit_drug_slope_random,
    factors_in_emm = c("drug", "dose"), # Get all cells
    at = list(monkey = "A"), # Condition on monkey A
    include_ev = TRUE # Or FALSE if not needed
  )

  observed_params_emms <- get_observed_demand_param_emms(
    fit_obj = fit_drug_slope_random,
    factors_in_emm = c("drug", "dose"), # Get all cells
    at = list(monkey = "A"), # Condition on monkey A
    include_ev = TRUE # Or FALSE if not needed
  )
} else {
  cat("Random-slope model failed to converge with the current settings.\n")
  fit_drug_slope_random <- NULL
}
```

Now let's demonstrate how to extract individual-level predicted coefficients using the
`get_individual_coefficients()` function. This combines fixed and random effects to
calculate individual-level parameter estimates for each subject.

```{r complex_random_effects_coefs, eval = !render_fast, include = !render_fast}
if (!is.null(fit_drug_slope_random) && !is.null(fit_drug_slope_random$model)) {
  # Get individual coefficients in wide format (default)
  individual_coefs_wide <- get_individual_coefficients(
    fit_drug_slope_random,
    params = c("Q0", "alpha"),
    format = "wide"
  )

  cat("Individual coefficients (wide format, log10 scale):\n")
  print(individual_coefs_wide)

  # Get individual coefficients in long format
  individual_coefs_long <- get_individual_coefficients(
    fit_drug_slope_random,
    params = c("Q0", "alpha"),
    format = "long"
  )

  cat("\nIndividual coefficients (long format, log10 scale):\n")
  print(head(individual_coefs_long, 10))

  # Convert to natural scale for interpretation
  individual_coefs_natural <- individual_coefs_wide

  # Convert Q0 coefficients to natural scale
  q0_cols <- grep(
    "^estimated_Q0_",
    names(individual_coefs_natural),
    value = TRUE
  )
  for (col in q0_cols) {
    natural_col <- gsub("estimated_Q0_", "Q0_natural_", col)
    individual_coefs_natural[[natural_col]] <- 10^individual_coefs_natural[[
      col
    ]]
  }

  # Convert alpha coefficients to natural scale
  alpha_cols <- grep(
    "^estimated_alpha_",
    names(individual_coefs_natural),
    value = TRUE
  )
  for (col in alpha_cols) {
    natural_col <- gsub("estimated_alpha_", "alpha_natural_", col)
    individual_coefs_natural[[natural_col]] <- 10^individual_coefs_natural[[
      col
    ]]
  }

  cat("\nExample individual coefficients on natural scale:\n")
  print(individual_coefs_natural[
    1:3,
    c(
      "id",
      grep("Q0_natural_", names(individual_coefs_natural), value = TRUE)[1:2],
      grep("alpha_natural_", names(individual_coefs_natural), value = TRUE)[1:2]
    )
  ])

  # You can also get coefficients for specific parameters only
  q0_only <- get_individual_coefficients(
    fit_drug_slope_random,
    params = "Q0",
    format = "wide"
  )

  cat("\nQ0 coefficients only:\n")
  print(q0_only)
} else {
  cat("Random-slope model not available for coefficient extraction.\n")
}
```

### Continuous Covariates and `fixed_rhs`

In addition to factor predictors, you can include continuous covariates in the fixed-effects
linear models for `Q0` and `alpha`.

There are two convenient ways to do this:

1. Add continuous covariates additively using the `continuous_covariates` argument (no formula writing required).
2. Specify a full right-hand-side using `fixed_rhs` (formula), which gives you complete control,
   including interactions (e.g., factor-by-continuous).

Below we demonstrate both, using a made-up continuous covariate `age` assigned per subject (`monkey`).

```{r continuous_setup, eval = !render_fast, include = !render_fast}
set.seed(123)

# Create a per-monkey continuous covariate: age
monkey_age <- unique(ko[, c("monkey")]) |>
  dplyr::mutate(age = round(runif(dplyr::n(), min = 18, max = 25), 1))

ko_age <- dplyr::left_join(ko, monkey_age, by = "monkey")
# Ensure a numeric version of dose is available for continuous modeling
ko_age <- ko_age |>
  dplyr::mutate(dose_num = suppressWarnings(as.numeric(as.character(dose))))

# Robust median for dose_num (fallbacks if all values are NA)
median_dose_num <- suppressWarnings(stats::median(
  ko_age$dose_num,
  na.rm = TRUE
))
if (is.na(median_dose_num) || !is.finite(median_dose_num)) {
  finite_vals <- ko_age$dose_num[is.finite(ko_age$dose_num)]
  if (length(finite_vals) > 0) {
    median_dose_num <- finite_vals[1]
  } else {
    median_dose_num <- 0
  }
}

# Sanity check
head(ko_age[, c("monkey", "age")])
```

#### A) Additive continuous covariate via `continuous_covariates`

Here we keep a single factor (`dose`) and add `age` as an additive continuous covariate. The RHS becomes
`~ 1 + dose + age` for both parameters.

```{r continuous_age_add, eval = !render_fast, include = !render_fast}
fit_age_add <- try(
  fit_demand_mixed(
    data = ko_age,
    y_var = "y_ll4",
    x_var = "x",
    id_var = "monkey",
    factors = "dose",
    continuous_covariates = c("age"),
    equation_form = "zben",
    nlme_control = quick_nlme_control,
    start_value_method = "pooled_nls"
  ),
  silent = TRUE
)

if (!inherits(fit_age_add, "try-error") && !is.null(fit_age_add$model)) {
  print(fit_age_add)

  # EMMs over dose at age = 21
  emms_age_add <- get_demand_param_emms(
    fit_obj = fit_age_add,
    factors_in_emm = c("dose"),
    at = list(age = 21),
    include_ev = TRUE
  )
  cat("\nEMMs (age = 21):\n")
  print(emms_age_add)

  # Pairwise comparisons for dose at age = 21 (log10 scale; ratios reported too)
  comps_age_add <- get_demand_comparisons(
    fit_obj = fit_age_add,
    compare_specs = ~dose,
    contrast_type = "pairwise",
    at = list(age = 21),
    adjust = "fdr",
    report_ratios = TRUE
  )
  cat("\nPairwise comparisons for dose (age = 21):\n")
  print(comps_age_add)

  # Plot population lines at age = 21 (natural scale on y)
  plot(
    fit_age_add,
    at = list(age = 21),
    inv_fun = ll4_inv,
    color_by = "dose",
    x_trans = "pseudo_log",
    y_trans = "pseudo_log",
    title = "Demand by Dose with Age as Covariate (Population Fit at age = 21)"
  )
} else {
  cat(
    "\nAdditive continuous covariate model failed to converge for the vignette settings.\n"
  )
}
```

#### B) fixed_rhs with a factor (drug) and a continuous covariate (dose)

Here we specify `drug` as a factor but treat `dose` as a continuous covariate
(`dose_num`, a numeric version of `dose`),
and include `age` as an additional continuous covariate. This uses a shared RHS
for `Q0` and `alpha`:
`~ 1 + drug + dose_num + age`.

```{r continuous_age_rhs, eval = !render_fast, include = !render_fast}
fit_age_rhs <- try(
  fit_demand_mixed(
    data = ko_age,
    y_var = "y_ll4",
    x_var = "x",
    id_var = "monkey",
    factors = c("drug"), # only 'drug' is treated as a factor
    fixed_rhs = "~ 1 + drug + dose_num + age",
    equation_form = "zben",
    nlme_control = quick_nlme_control,
    start_value_method = "pooled_nls"
  ),
  silent = TRUE
)

if (!inherits(fit_age_rhs, "try-error") && !is.null(fit_age_rhs$model)) {
  print(fit_age_rhs)

  # EMMs over drug (factor) at age = 21 and dose set to its median
  emms_age_rhs <- get_demand_param_emms(
    fit_obj = fit_age_rhs,
    factors_in_emm = c("drug"),
    at = list(age = 21, dose_num = median_dose_num),
    include_ev = TRUE
  )
  cat("\nEMMs (drug; age = 21, dose_num = median):\n")
  print(emms_age_rhs)

  # Pairwise comparisons: drug at age = 21 and dose set to its median
  comps_drug_by_dose_age <- get_demand_comparisons(
    fit_obj = fit_age_rhs,
    compare_specs = ~drug,
    contrast_type = "pairwise",
    at = list(age = 21, dose_num = median_dose_num),
    adjust = "fdr",
    report_ratios = TRUE
  )
  cat("\nPairwise comparisons of drug (age = 21, dose_num = median):\n")
  print(comps_drug_by_dose_age)

  # Plot population lines at age = 21 and dose_num = median, color by drug and facet by drug
  plot(
    fit_age_rhs,
    at = list(age = 21, dose_num = median_dose_num),
    inv_fun = ll4_inv,
    color_by = "drug",
    facet_formula = ~drug,
    x_trans = "pseudo_log",
    y_trans = "pseudo_log",
    title = "Demand by Drug with Dose Continuous (Population Fit at age = 21; dose_num = median)"
  )
} else {
  cat(
    "\nfixed_rhs model (drug factor, dose continuous) failed to converge for the vignette settings.\n"
  )
}
```

### Trends with `emtrends`

We can examine how the parameters change with continuous covariates using
`get_demand_param_trends()`, which wraps `emmeans::emtrends()` and returns
tidy results for Q0 and alpha trends on the log10 scale.

Below we compute trends with respect to `age` and `dose_num`, first overall and
then by `drug` (if the model includes `drug` as a factor).

```{r trends, eval = !render_fast, include = !render_fast}
# Choose a fitted model for trends (prefer fixed_rhs fit)
fit_for_trends <- NULL
if (
  exists("fit_age_rhs") && !is.null(fit_age_rhs) && !is.null(fit_age_rhs$model)
) {
  fit_for_trends <- fit_age_rhs
} else if (
  exists("fit_age_add") && !is.null(fit_age_add) && !is.null(fit_age_add$model)
) {
  fit_for_trends <- fit_age_add
}

if (!is.null(fit_for_trends)) {
  # Determine which of the candidate covariates are actually present in the model RHS
  # Use Q0 formula (or alpha if Q0 is NULL) - they may differ if collapse_levels was used
  rhs_formula_str <- fit_for_trends$formula_details$fixed_effects_formula_str_Q0
  if (
    is.null(rhs_formula_str) ||
      is.na(rhs_formula_str) ||
      !nzchar(rhs_formula_str)
  ) {
    rhs_formula_str <- fit_for_trends$formula_details$fixed_effects_formula_str_alpha
  }
  rhs_vars_tr <- tryCatch(
    all.vars(as.formula(rhs_formula_str)),
    error = function(e) character(0)
  )
  covariates_candidate <- c("age", "dose_num")
  covariates_in_model <- intersect(rhs_vars_tr, covariates_candidate)

  if (length(covariates_in_model) > 0) {
    cat("\n--- Trends (overall) at age = 21 and dose_num = median ---\n")
    trends_overall <- get_demand_param_trends(
      fit_obj = fit_for_trends,
      params = c("Q0", "alpha"),
      covariates = covariates_in_model,
      specs = ~1,
      at = list(age = 21, dose_num = median_dose_num)
    )
    print(trends_overall)

    # If 'drug' is a model factor, show trends by drug as well
    model_factors_tr <- fit_for_trends$param_info$factors
    if (!is.null(model_factors_tr) && ("drug" %in% model_factors_tr)) {
      cat("\n--- Trends by drug at age = 21 and dose_num = median ---\n")
      trends_by_drug <- get_demand_param_trends(
        fit_obj = fit_for_trends,
        params = c("Q0", "alpha"),
        covariates = covariates_in_model,
        specs = ~drug,
        at = list(age = 21, dose_num = median_dose_num)
      )
      print(trends_by_drug)
    } else {
      cat(
        "\nNote: The selected model does not include 'drug' as a factor; skipping by-drug trends.\n"
      )
    }
  } else {
    cat(
      "\nNo candidate continuous covariates ('age', 'dose_num') are present in the model RHS; skipping trends.\n"
    )
  }
} else {
  cat(
    "\nNo suitable fitted model available for computing trends in the vignette.\n"
  )
}
```

# See Also

- `vignette("mixed-demand")` -- Basic mixed-effects demand modeling
- `vignette("model-selection")` -- Choosing the right model class
- `vignette("fixed-demand")` -- Fixed-effect demand modeling
- `vignette("hurdle-demand-models")` -- Two-part hurdle demand models
- `vignette("cross-price-models")` -- Cross-price demand analysis
- `vignette("group-comparisons")` -- Group comparisons
- `vignette("migration-guide")` -- Migrating from `FitCurves()`
- `vignette("beezdemand")` -- Getting started with beezdemand