Modeling Enodwment, Distance, and Income Effects in a Spatially Explicit Discrete Choice Experiment

Baseline utility functions

Function Type	Function \(f(\text{NC},\Phi)\)	Parameters (‘\(\Phi\)‘)
Linear	\(\beta_{\text{NC}} \cdot \text{NC}\)	\(\{\beta_{\text{NC}}\}\)
Quadratic Utility	\(\beta_{\text{NC}} \cdot \text{NC} + \beta_{\text{NC}_{\text{sq}}} \cdot \text{NC}^2\)	\(\{\beta_{\text{NC}}, \beta_{\text{NC}_{\text{sq}}}\}\)
Logarithmic	\(\beta_{\text{NC}} \cdot \log(\text{NC})\)	\(\{\beta_{\text{NC}}\}\)
Box-Cox	\(\beta_{\text{NC}} \cdot \frac{\text{NC}^\lambda - 1}{\lambda}\)	\(\{\beta_{\text{NC}}, \lambda\}\)
Log-Linear	\(\beta_{\text{NC}} \cdot \text{NC} + \beta_{\text{NC}_{\text{log}}} \cdot \log(\text{NC})\)	\(\{\beta_{\text{NC}}, \beta_{\text{NC}_{\text{log}}}\}\)

Code

# Load the necessary library
library(kableExtra)


Attaching package: 'kableExtra'

The following object is masked from 'package:dplyr':

    group_rows

Code

# Create a data frame with model names, AIC, BIC, and log-likelihood (LLout)
model_names <- c("Quadratic Utility Function",
                 "Log Utility Function",
                 "Linear Utility Function",
                 "Box Cox Utility Function",
                 "Log-Linear Utility Function")

# Extract AIC, BIC, and LLout values from the models
aic_values <- c(mxl_access$AIC,
                mxl_access_log$AIC,
                mxl_access_linear$AIC,
                mxl_box_cox$AIC,
                mxl_access_log_linear$AIC)

bic_values <- c(mxl_access$BIC,
                mxl_access_log$BIC,
                mxl_access_linear$BIC,
                mxl_box_cox$BIC,
                mxl_access_log_linear$BIC)

ll_values <- c(mxl_access$LLout,
               mxl_access_log$LLout,
               mxl_access_linear$LLout,
               mxl_box_cox$LLout,
               mxl_access_log_linear$LLout)

# Combine into a data frame
model_table <- data.frame(Model = model_names, AIC = aic_values, BIC = bic_values, LLout = ll_values)

# Round the values to 2 decimal places
model_table$AIC <- round(model_table$AIC, 2)
model_table$BIC <- round(model_table$BIC, 2)
model_table$LLout <- round(model_table$LLout, 2)

# Find the min and max for AIC, BIC, and LLout
min_aic <- min(model_table$AIC)
max_aic <- max(model_table$AIC)
min_bic <- min(model_table$BIC)
max_bic <- max(model_table$BIC)
min_ll <- min(model_table$LLout)
max_ll <- max(model_table$LLout)

# Use cell_spec to color only the cells with min/max values
model_table$AIC <- cell_spec(model_table$AIC, "html", 
                             color = ifelse(model_table$AIC == min_aic, "green", 
                                            ifelse(model_table$AIC == max_aic, "red", "black")))

model_table$BIC <- cell_spec(model_table$BIC, "html", 
                             color = ifelse(model_table$BIC == min_bic, "green", 
                                            ifelse(model_table$BIC == max_bic, "red", "black")))

model_table$LLout <- cell_spec(model_table$LLout, "html", 
                               color = ifelse(model_table$LLout == min_ll, "red", 
                                              ifelse(model_table$LLout == max_ll, "green", "black")))

# Create the table with kable and format the cells
model_table %>%
  kable("html", escape = F) %>%
  kable_styling(full_width = F)

Model	AIC	BIC	LLout
Quadratic Utility Function	128533.25	128809.48	-64238.63
Log Utility Function	128773.39	128950.97	-64368.7
Linear Utility Function	128574.99	128752.57	-64269.5
Box Cox Utility Function	128548.7	128775.6	-64251.35
Log-Linear Utility Function	128554.29	128830.52	-64249.15

Ways to integrate distance

Interactions in quadratic utility function

Interaction Type	Function \(f(\text{NC},\Phi)\)	Radius Parameters (\(\Phi\))
Two linear interactions	\((\beta_{\text{NC}} + \delta \cdot \text{D}) \cdot \text{NC} + (\beta_{\text{NC\_sq}}+ \delta_{\text{sq}} \cdot \text{D}) \cdot \text{NC}^2\)	\(\{\delta,\delta_{\text{sq}}\}\)
Two log interactions	\((\beta_{\text{NC}} + \delta \cdot log(\text{D})) \cdot \text{NC} + (\beta_{\text{NC\_sq}}+ \delta_{\text{sq}} \cdot log(\text{D})) \cdot \text{NC}^2\)	\(\{\delta,\delta_{\text{sq}}\}\)
Radius also squared	\((\beta_{\text{NC}} + \delta_1 \cdot D + \delta_2 \cdot D^2)\cdot \text{NC} \;+\; (\beta_{\text{NC\_sq}} + \delta_{\text{sq},1}\cdot D + \delta_{\text{sq},2}\cdot D^2)\cdot \text{NC}^2\)	\(\{\delta_1,\delta_2, \delta_{\text{sq},1}\, \delta_{\text{sq},2}\}\)
One exponential decay	\(\exp(-\lambda \cdot \text{D})\left( \beta_{\text{NC}} \cdot \text{NC} + \beta_{\text{NC\_sq}} \cdot \text{NC}^2 \right)\)	\(\{\lambda\}\)

Interactions log utility function

Interaction Type	Function \(f(\text{NC},\Phi)\)	Radius Parameters (\(\Phi\))
One log interaction	\((\beta_{\text{NC}} + \delta \cdot log(D)) \cdot \log(\text{NC})\)	\(\{\delta\}\)

Plots

Two linear interactions

Heatmap plot

Code

plot_wtp_heatmap(
  model = radius_mean_model,
  attr = "pa_full",
  endowment_center = 11.7,
  radius_spec = "linear")

WTP functions for fixed radius

Code

pa <- plot_wtp_lines(
  model = radius_mean_model,
  attr = "pa",
  endowment_center = 11.7,
  radius_spec = "linear",
  median_line = 11
)

pa_half <- plot_wtp_lines(
  model = radius_mean_model,
  attr = "pa_half",
  endowment_center = 11.7,
  radius_spec = "linear",
  median_line = 11
)

pa_full <- plot_wtp_lines(
  model = radius_mean_model,
  attr = "pa_full",
  endowment_center = 11.7,
  radius_spec = "linear",
  median_line = 11
)

hnv <- plot_wtp_lines(
  model = radius_mean_model,
  attr = "hnv",
  endowment_center = 9.9,
  radius_spec = "linear",
  median_line = 9.9
)

hnv_vis <- plot_wtp_lines(
  model = radius_mean_model,
  attr = "hnv_vis",
  endowment_center = 9.9,
  radius_spec = "linear",
  median_line = 9.9
)

pa

Code

pa_half

Code

pa_full

Code

hnv

Code

hnv_vis

WTP functions for fixed endowment

Code

plot_wtp_distance_lines(
  model = radius_mean_model,
  attr = "pa",
  fixed_endowments = c(0, 10, 20, 40, 50),
  endowment_center = 11.7,
  radius_spec = "linear"
)

Code

plot_wtp_distance_lines(
  model = radius_mean_model,
  attr = "hnv",
  fixed_endowments = c(0, 10, 20, 40, 50),
  endowment_center = 9.9,
  radius_spec = "linear"
)

Two log interactions

WTP functions for fixed radius

Code

pa <- plot_wtp_lines(
  model = radius_mean_model_log,
  attr = "pa",
  endowment_center = 11.7,
  radius_spec = "log",
  median_line = 11
)

pa_half <- plot_wtp_lines(
  model = radius_mean_model_log,
  attr = "pa_half",
  endowment_center = 11.7,
  radius_spec = "log",
  median_line = 11
)

pa_full <- plot_wtp_lines(
  model = radius_mean_model_log,
  attr = "pa_full",
  endowment_center = 11.7,
  radius_spec = "log",
  median_line = 11
)

hnv <- plot_wtp_lines(
  model = radius_mean_model_log,
  attr = "hnv",
  endowment_center = 9.9,
  radius_spec = "log",
  median_line = 9.9
)

hnv_vis <- plot_wtp_lines(
  model = radius_mean_model_log,
  attr = "hnv_vis",
  endowment_center = 9.9,
  radius_spec = "log",
  median_line = 9.9
)

pa

Code

pa_half

Code

pa_full

Code

hnv

Code

hnv_vis

Exp model

WTP functions for fixed radius

Code

plot_wtp_lines_lambda(
  model = lambda_model,
  attr = "pa",
  endowment_center = 11.7,
  median_line = 11
)

WTP functions for fixed endowment

Code

plot_wtp_distance_lines_lambda(
  model = lambda_model,
  attr = "pa",
  fixed_endowments = c(0, 10, 20, 40, 50, 80),
  endowment_center = 11.7
)

Log model

WTP functions for fixed distance

Code

plot_wtp_lines_log(
  model = radius_log_model,
  attr = "pa",
  endowment_seq = seq(1, 100, by = 1),
  distance_values = c(0, 10, 20, 30, 50),
  median_line = 11
)

Code

plot_wtp_lines_log(
  model = radius_log_model,
  attr = "hnv",
  endowment_seq = seq(1, 100, by = 1),
  distance_values = c(0, 10, 20, 30, 50),
  median_line = 10
)

WTP functions for fixed endowment

Code

plot_wtp_distance_lines_log(
  model = radius_log_model,
  attr = "pa",
  fixed_endowments = c(1, 5, 10, 20, 40, 60)
)

Code

plot_wtp_distance_lines_log(
  model = radius_log_model,
  attr = "hnv",
  fixed_endowments = c(1, 5, 10, 20, 40, 60)
)

What do we need for the aggregation?

So far: WTP for all cells up to 15km is aggregated
WTP within cell only depends on endowment not on distance
Distance could enter in step 1 and/or 4/5 of aggregation?
Ideally enters in step 4/5, marginal cell WTP only dependent on endowment but cells that benefit from change within one cell should be selected based on distance
Close cells have a higher WTP than cells further away with WTP = 0 after some distance

Integrate Income

How to Intergrate Income?

Income can have effect on marginal utility of income, on utility derived from NC attributes or on both

Effect of income on NC attributes

Code

linear_inc <- apollo_loadModel("Results/MXL/Income_interactions/MXL wtp linear radius income")

Successfully loaded
  /home/sc.uni-leipzig.de/nc71qaxa/valugaps/Results/MXL/Income_interactions/MXL
  wtp linear radius income_model.rds

Code

plot_wtp_income <- function(model, param) {
  library(ggplot2)
  
  # get coefficients
  b <- if (!is.null(model$estimate)) model$estimate else model
  
  # mapping for your current model names
  name_map <- list(
    pa = list(base = "beta_pa", inc = "pa_inc", label = "PA"),
    hnv = list(base = "beta_hnv", inc = "hnv_inc", label = "HNV"),
    pa_half = list(base = "pa_half_access", inc = "pa_half_inc", label = "PA half access"),
    pa_full = list(base = "pa_full_access", inc = "pa_full_inc", label = "PA full access"),
    hnv_visible = list(base = "hnv_visible", inc = "hnv_vis_inc", label = "HNV visible")
  )
  
  if (!param %in% names(name_map)) {
    stop(paste0(
      "param must be one of: ",
      paste(names(name_map), collapse = ", ")
    ))
  }
  
  base_name <- name_map[[param]]$base
  inc_name  <- name_map[[param]]$inc
  lab       <- name_map[[param]]$label
  
  if (!base_name %in% names(b)) stop(paste("Missing coefficient:", base_name))
  if (!inc_name %in% names(b))  stop(paste("Missing coefficient:", inc_name))
  
  income_seq <- seq(-3000, 3000, by = 50)
  
  df <- data.frame(
    income = income_seq,
    wtp = unname(b[base_name]) + unname(b[inc_name]) * income_seq
  )
  
  ggplot(df, aes(x = income, y = wtp)) +
    geom_line(linewidth = 1.2) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
    geom_vline(xintercept = 0, linetype = "dotted", color = "grey50") +
    labs(
      title = paste("WTP over monthly income:", lab),
      x = "Monthly income (mean-centered)",
      y = "Willingness to pay"
    ) +
    theme_minimal(base_size = 13) +
    theme(
      plot.title = element_text(face = "bold"),
      panel.grid.minor = element_blank()
    )
}

plot_wtp_income(linear_inc, "pa")

Code

plot_wtp_income(linear_inc, "pa_half")

Code

plot_wtp_income(linear_inc, "pa_full")

Code

plot_wtp_income(linear_inc, "hnv")

Code

plot_wtp_income(linear_inc, "hnv_visible")

Effect on income only on cost

Code

pref_income <- apollo_loadModel("Results/MXL/MXL pref income")

Successfully loaded
  /home/sc.uni-leipzig.de/nc71qaxa/valugaps/Results/MXL/MXL pref
  income_model.rds

Code

plot_wtp_pref <- function(model,
                          attr = c("pa", "hnv", "pa_half_access", "pa_full_access", "hnv_visible"),
                          income_levels = c(0, 30, 60),
                          endowment = seq(0, 100, by = 1)) {
  
  
  b <- if (!is.null(model$estimate)) model$estimate else model
  attr <- match.arg(attr)
  
  # coefficient names in your estimated model
  coef_map <- list(
    pa = list(base = "beta_pa", sq = "beta_pa_sq", label = "PA"),
    hnv = list(base = "beta_hnv", sq = "beta_hnv_sq", label = "HNV"),
    pa_half_access = list(base = "pa_half_access", sq = "pa_half_access_sq", label = "PA half access"),
    pa_full_access = list(base = "pa_full_access", sq = "pa_full_access_sq", label = "PA full access"),
    hnv_visible = list(base = "hnv_visible", sq = "hnv_visible_sq", label = "HNV visible")
  )
  
  base_name <- coef_map[[attr]]$base
  sq_name   <- coef_map[[attr]]$sq
  lab       <- coef_map[[attr]]$label
  
  # checks
  needed <- c(base_name, sq_name, "beta_cost", "beta_income_cost")
  missing <- needed[!needed %in% names(b)]
  if (length(missing) > 0) {
    stop("Missing coefficients: ", paste(missing, collapse = ", "))
  }
  
  # create grid
  df <- expand.grid(
    endowment = endowment,
    income = income_levels
  )
  
  # exact formula as you wrote it
  df$wtp <- -(
    unname(b[base_name]) +
      2 * unname(b[sq_name]) * df$endowment
  ) / (
    -exp(unname(b["beta_cost"])) +
      unname(b["beta_income_cost"]) * df$income
  )
  
  df$income <- factor(df$income, levels = income_levels)
  
  ggplot(df, aes(x = endowment, y = wtp, color = income)) +
    geom_line(linewidth = 1.1) +
    geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
    labs(
      title = paste("WTP for", lab),
      subtitle = "WTP depending on endowment for different income levels",
      x = paste(lab, "endowment"),
      y = "WTP",
      color = "Income"
    ) +
    theme_minimal(base_size = 13) +
    theme(
      plot.title = element_text(face = "bold"),
      panel.grid.minor = element_blank(),
      legend.position = "bottom"
    )
}

plot_wtp_pref(pref_income, attr = "pa", income_levels = c(0, 30, 60), endowment = seq(0, 100, 1))

Code

plot_wtp_pref(pref_income, attr = "hnv", income_levels = c(0, 30, 60), endowment = seq(0, 100, 1))

Code

plot_wtp_pref(pref_income, attr = "pa_half_access", income_levels = c(0, 30, 60), endowment = seq(0, 100, 1))

Code

plot_wtp_pref(pref_income, attr = "pa_full_access", income_levels = c(0, 30, 60), endowment = seq(0, 100, 1))

Code

plot_wtp_pref(pref_income, attr = "hnv_visible", income_levels = c(0, 30, 60), endowment = seq(0, 100, 1))