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}}}\}\)
Cubic Utility	\(\beta_{\text{NC}} \cdot \text{NC} + \beta_{\text{NC}_{\text{sq}}} \cdot \text{NC}^2 + \beta_{\text{NC}_{\text{cu}}} \cdot \text{NC}^3\)	\(\{\beta_{\text{NC}}, \beta_{\text{NC}_{\text{sq}}}, \beta_{\text{NC}_{\text{cu}}}\}\)

# Load the necessary library
library(kableExtra)

systemfonts and textshaping have been compiled with different versions of Freetype. Because of this, textshaping will not use the font cache provided by systemfonts

# 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",
                 "Cubic 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,
                mxl_cubic$AIC)

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

ll_values <- c(mxl_access$LLout,
               mxl_access_log$LLout,
               mxl_access_linear$LLout,
               mxl_box_cox$LLout,
               mxl_access_log_linear$LLout,
               mxl_cubic$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	13898.94	14028.73	-6931.47
Box Cox Utility Function	128548.7	128775.6	-64251.35
Log-Linear Utility Function	128554.29	128830.52	-64249.15
Cubic Utility Function	128565.43	128940.32	-64244.72

Plots

### Quadratic
# Generalized version for Quadratic model calculations
wtp_quad <- function(x, y, beta_name, beta_sq_name) {
  y$estimate[[beta_name]] + 2 * y$estimate[[beta_sq_name]] * x
}
### Linear

wtp_lin <- function(y, name) {
  y$estimate[[name]]
}


### Log
wtp_log    <- function(x, y , name) {(y$estimate[[name]]/(x))} # advanced version




### Box-Cox



wtp_box_cox <- function(x, y, beta_name, lambda_name) {
  y$estimate[[beta_name]] * (x^(y$estimate[[lambda_name]] - 1))
}



### Log Linear
wtp_log_linear    <- function(x, y , name1, name2) {y$estimate[[name1]] + (y$estimate[[name2]]/(x))} # advanced version


### Cubic

wtp_cubic <- function(x, y, beta_name, beta_sq_name, beta_cu_name) {
  y$estimate[[beta_name]] + 2 * y$estimate[[beta_sq_name]] * x + 3 * y$estimate[[beta_cu_name]]*x^2
}

## generating data vectors for each line

data_vector <- c(1:100)

data_vec_wtp_hnv <- wtp_quad(data_vector, mxl_access, "beta_hnv", "beta_hnv_sq")
data_vec_wtp_pa <- wtp_quad(data_vector, mxl_access, "beta_pa", "beta_pa_sq")

data_vec_wtp_hnv_vis <- wtp_quad(data_vector, mxl_access, "hnv_visible", "hnv_visible_sq")
data_vec_wtp_pa_half <- wtp_quad(data_vector, mxl_access, "pa_half_access", "pa_half_access_sq")
data_vec_wtp_pa_full <- wtp_quad(data_vector, mxl_access, "pa_full_access", "pa_full_access_sq")




data_vec_wtp_hnv_lin <- wtp_lin(mxl_access_linear, "beta_hnv")
data_vec_wtp_pa_lin <- wtp_lin(mxl_access_linear, "beta_pa")

data_vec_wtp_hnv_lin_vis <- wtp_lin(mxl_access_linear, "hnv_visible")
data_vec_wtp_pa_lin_half <- wtp_lin(mxl_access_linear, "pa_half_access")
data_vec_wtp_pa_lin_full <- wtp_lin(mxl_access_linear, "pa_full_access")


data_vec_wtp_hnv_log <- wtp_log(data_vector, mxl_access_log, "beta_hnv")
data_vec_wtp_pa_log <- wtp_log(data_vector, mxl_access_log, "beta_pa")

data_vec_wtp_hnv_log_vis <- wtp_log(data_vector, mxl_access_log, "hnv_visible")
data_vec_wtp_pa_log_half <- wtp_log(data_vector, mxl_access_log, "pa_half_access")
data_vec_wtp_pa_log_full <- wtp_log(data_vector, mxl_access_log, "pa_full_access")

data_vec_wtp_hnv_box_cox <- wtp_box_cox(data_vector, mxl_box_cox, "beta_hnv", "lambda_hnv")
data_vec_wtp_pa_box_cox <- wtp_box_cox(data_vector, mxl_box_cox, "beta_pa", "lambda_pa")

data_vec_wtp_hnv_box_cox_vis <- wtp_box_cox(data_vector, mxl_box_cox, "hnv_visible", "lambda_hnv_vis")
data_vec_wtp_pa_box_cox_half <- wtp_box_cox(data_vector, mxl_box_cox, "pa_half_access", "lambda_pa_half")
data_vec_wtp_pa_box_cox_full <- wtp_box_cox(data_vector, mxl_box_cox, "pa_full_access", "lambda_pa_full")


data_vec_wtp_hnv_log_linear <- wtp_log_linear(data_vector, mxl_access_log_linear, "beta_hnv", "beta_hnv_log")
data_vec_wtp_pa_log_linear <- wtp_log_linear(data_vector, mxl_access_log_linear, "beta_pa", "beta_pa_log")


data_vec_wtp_hnv_log_linear_vis <- wtp_log_linear(data_vector, mxl_access_log_linear, "hnv_visible", "hnv_visible_log")
data_vec_wtp_pa_log_linear_half <- wtp_log_linear(data_vector, mxl_access_log_linear, "pa_half_access", "pa_half_access_log")
data_vec_wtp_pa_log_linear_full <- wtp_log_linear(data_vector, mxl_access_log_linear, "pa_full_access", "pa_full_access_log")

data_vec_wtp_hnv_cu <- wtp_cubic(data_vector, mxl_cubic, "beta_hnv", "beta_hnv_sq", "beta_hnv_cu")
data_vec_wtp_pa_cu <- wtp_cubic(data_vector, mxl_cubic, "beta_pa", "beta_pa_sq", "beta_pa_cu")

data_vec_wtp_hnv_cu_vis <- wtp_cubic(data_vector, mxl_cubic, "hnv_visible", "hnv_visible_sq", "hnv_visible_cu")
data_vec_wtp_pa_cu_half <- wtp_cubic(data_vector, mxl_cubic, "pa_half_access", "pa_half_access_sq", "pa_half_access_cu")
data_vec_wtp_pa_cu_full <- wtp_cubic(data_vector, mxl_cubic, "pa_full_access", "pa_full_access_sq", "pa_full_access_cu")

hnv_fun_all <- ggplot() +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv, col="Quadratic")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_lin, col="Linear")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_log, col="Log")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_box_cox, col="Box Cox")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_log_linear, col="Log Linear")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_cu, col="Cubic")) +
  # geom_jitter(aes(x = database$sq_hnv_area/1000, y = 0), height = 0.2, alpha = 0.6, color = "peru") +
  geom_density(aes(x = database$sq_hnv_area/1000, y = after_stat(scaled * 5)), fill = "peru", alpha = 0.5) +
  geom_vline(aes(xintercept = median(database$sq_hnv_area/1000)), linetype = "dashed", size = 0.3) +
  geom_hline(yintercept = 0, color = "black") +
  labs(color="Specification") +
  xlab("Status Quo HNV in 1000ha") +
  ylab("WTP (€/1000ha)") +
  coord_cartesian(ylim = c(-3, 20)) + theme_minimal()

hnv_fun_all_vis <- ggplot() +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_vis, col="Quadratic")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_lin_vis, col="Linear")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_log_vis, col="Log")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_box_cox_vis, col="Box Cox")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_log_linear_vis, col="Log Linear")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_hnv_cu_vis, col="Cubic")) +
  # geom_jitter(aes(x = database$sq_hnv_area/1000, y = 0), height = 0.2, alpha = 0.6, color = "peru") +
  geom_density(aes(x = database$sq_hnv_area/1000, y = after_stat(scaled * 5)), fill = "peru", alpha = 0.5) +
  geom_vline(aes(xintercept = median(database$sq_hnv_area/1000)), linetype = "dashed", size = 0.3) +
  geom_hline(yintercept = 0, color = "black") +
  labs(color="Specification") +
  xlab("Status Quo HNV in 1000ha") +
  ylab("WTP (€/1000ha)") +
  coord_cartesian(ylim = c(-3, 20)) + theme_minimal()

pa_fun_all <- ggplot() +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa, col="Quadratic")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_lin, col="Linear")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_log, col="Log")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_box_cox, col="Box Cox")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_log_linear, col="Log Linear")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_cu, col="Cubic")) +
  # geom_jitter(aes(x = database$sq_pa_area/1000, y = 0), height = 0.2, alpha = 0.6, color = "forestgreen") +
  geom_density(aes(x = database$sq_pa_area/1000, y = after_stat(scaled * 5)), fill = "forestgreen", alpha = 0.5) +
  geom_vline(aes(xintercept = median(database$sq_pa_area/1000)), linetype = "dashed", size = 0.3) +
  geom_hline(yintercept = 0, color = "black") +
  labs(color="Specification") +
  xlab("Status Quo PA in 1000ha") +
  ylab("WTP (€/1000ha)") +
  coord_cartesian(ylim = c(-5, 30)) + theme_minimal()

pa_fun_all_half <- ggplot() +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_half, col="Quadratic")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_lin_half, col="Linear")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_log_half, col="Log")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_box_cox_half, col="Box Cox")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_log_linear_half, col="Log Linear")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_cu_half, col="Cubic")) +
  # geom_jitter(aes(x = database$sq_pa_area/1000, y = 0), height = 0.2, alpha = 0.6, color = "forestgreen") +
  geom_density(aes(x = database$sq_pa_area/1000, y = after_stat(scaled * 5)), fill = "forestgreen", alpha = 0.5) +
  geom_vline(aes(xintercept = median(database$sq_pa_area/1000)), linetype = "dashed", size = 0.3) +
  geom_hline(yintercept = 0, color = "black") +
  labs(color="Specification") +
  xlab("Status Quo PA in 1000ha") +
  ylab("WTP (€/1000ha)") +
  coord_cartesian(ylim = c(-5, 30)) + theme_minimal()

pa_fun_all_full <- ggplot() +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_full, col="Quadratic")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_lin_full, col="Linear")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_log_full, col="Log")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_box_cox_full, col="Box Cox")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_log_linear_full, col="Log Linear")) +
  geom_line(aes(x=data_vector, y=data_vec_wtp_pa_cu_full, col="Cubic")) +
  # geom_jitter(aes(x = database$sq_pa_area/1000, y = 0), height = 0.2, alpha = 0.6, color = "forestgreen") +
  geom_density(aes(x = database$sq_pa_area/1000, y = after_stat(scaled * 5)), fill = "forestgreen", alpha = 0.5) +
  geom_vline(aes(xintercept = median(database$sq_pa_area/1000)), linetype = "dashed", size = 0.3) +
  geom_hline(yintercept = 0, color = "black") +
  labs(color="Specification") +
  xlab("Status Quo PA in 1000ha") +
  ylab("WTP (€/1000ha)") +
  coord_cartesian(ylim = c(-5, 30)) + theme_minimal()

hnv_fun_all

hnv_fun_all_vis

pa_fun_all

pa_fun_all_half

pa_fun_all_full

Plot WTP distributions

database_wtp <- database %>%
  mutate(
    sq_pa_area_1000  = sq_pa_area / 1000,
    sq_hnv_area_1000 = sq_hnv_area / 1000,

    # PA WTPs
    wtp_pa_quadratic   = wtp_quad(sq_pa_area_1000,  mxl_access,            "beta_pa",        "beta_pa_sq"),
    wtp_pa_linear      = wtp_lin(mxl_access_linear,                       "beta_pa"),
    wtp_pa_log         = wtp_log(sq_pa_area_1000,   mxl_access_log,        "beta_pa"),
    wtp_pa_box_cox     = wtp_box_cox(sq_pa_area_1000, mxl_box_cox,         "beta_pa",        "lambda_pa"),
    wtp_pa_log_linear  = wtp_log_linear(sq_pa_area_1000, mxl_access_log_linear, "beta_pa", "beta_pa_log"),

    # HNV WTPs
    wtp_hnv_quadratic  = wtp_quad(sq_hnv_area_1000, mxl_access,            "beta_hnv",       "beta_hnv_sq"),
    wtp_hnv_linear     = wtp_lin(mxl_access_linear,                       "beta_hnv"),
    wtp_hnv_log        = wtp_log(sq_hnv_area_1000,  mxl_access_log,        "beta_hnv"),
    wtp_hnv_box_cox    = wtp_box_cox(sq_hnv_area_1000, mxl_box_cox,        "beta_hnv",       "lambda_hnv"),
    wtp_hnv_log_linear = wtp_log_linear(sq_hnv_area_1000, mxl_access_log_linear, "beta_hnv", "beta_hnv_log"),
    
    wtp_pa_full_quadratic   = wtp_quad(sq_pa_area_1000, mxl_access, "pa_full_access", "pa_full_access_sq"),
    wtp_pa_full_linear      = wtp_lin(mxl_access_linear, "pa_full_access"),
    wtp_pa_full_log         = wtp_log(sq_pa_area_1000, mxl_access_log, "pa_full_access"),
    wtp_pa_full_box_cox     = wtp_box_cox(sq_pa_area_1000, mxl_box_cox, "pa_full_access", "lambda_pa_full"),
    wtp_pa_full_log_linear  = wtp_log_linear(sq_pa_area_1000, mxl_access_log_linear, "pa_full_access", "pa_full_access_log"),

    wtp_pa_half_quadratic   = wtp_quad(sq_pa_area_1000, mxl_access, "pa_half_access", "pa_half_access_sq"),
    wtp_pa_half_linear      = wtp_lin(mxl_access_linear, "pa_half_access"),
    wtp_pa_half_log         = wtp_log(sq_pa_area_1000, mxl_access_log, "pa_half_access"),
    wtp_pa_half_box_cox     = wtp_box_cox(sq_pa_area_1000, mxl_box_cox, "pa_half_access", "lambda_pa_half"),
    wtp_pa_half_log_linear  = wtp_log_linear(sq_pa_area_1000, mxl_access_log_linear, "pa_half_access", "pa_half_access_log"),

    wtp_hnv_vis_quadratic   = wtp_quad(sq_hnv_area_1000, mxl_access, "hnv_visible", "hnv_visible_sq"),
    wtp_hnv_vis_linear      = wtp_lin(mxl_access_linear, "hnv_visible"),
    wtp_hnv_vis_log         = wtp_log(sq_hnv_area_1000, mxl_access_log, "hnv_visible"),
    wtp_hnv_vis_box_cox     = wtp_box_cox(sq_hnv_area_1000, mxl_box_cox, "hnv_visible", "lambda_hnv_vis"),
    wtp_hnv_vis_log_linear  = wtp_log_linear(sq_hnv_area_1000, mxl_access_log_linear, "hnv_visible", "hnv_visible_log")
  )



wtp_long <- database_wtp %>%
  select(starts_with("wtp_pa_"), starts_with("wtp_hnv_")) %>%
  pivot_longer(
    cols = everything(),
    names_to = "variable",
    values_to = "wtp"
  ) %>%
  mutate(
    type = case_when(
      grepl("^wtp_pa_", variable)  ~ "PA",
      grepl("^wtp_hnv_", variable) ~ "HNV"
    ),
    attribute = case_when(
      grepl("^wtp_pa_full_", variable) ~ "PA Full",
      grepl("^wtp_pa_half_", variable) ~ "PA Half",
      grepl("^wtp_pa_", variable)      ~ "PA",
      grepl("^wtp_hnv_vis_", variable) ~ "HNV Visible",
      grepl("^wtp_hnv_", variable)     ~ "HNV"
    ),
    functional_form = case_when(
      grepl("_log_linear$", variable) ~ "Log-Linear",
      grepl("_box_cox$", variable)    ~ "Box-Cox",
      grepl("_quadratic$", variable)  ~ "Quadratic",
      grepl("_linear$", variable)     ~ "Linear",
      grepl("_log$", variable)        ~ "Log"
    )
  )

ggplot(wtp_long, aes(x = wtp, color = functional_form, fill = functional_form)) +
  geom_density(alpha = 0.15, linewidth = 0.8) +
  facet_wrap(~ attribute, scales = "free") +
  labs(
    x = "WTP",
    y = "Density",
    color = "Functional form",
    fill = "Functional form"
  ) +
  theme_minimal() +
  coord_cartesian(xlim = c(-5, 30))

ggplot(wtp_long, aes(x = functional_form, y = wtp, fill = functional_form)) +
  geom_boxplot(outliers = FALSE) +
  facet_wrap(~ attribute, scales = "free") +
  labs(
    x = "Functional form",
    y = "WTP",
    fill = "Functional form"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

wtp_summary_table <- wtp_long %>%
  group_by(attribute, functional_form) %>%
  summarise(
    Min = min(wtp, na.rm = TRUE),
    Max = max(wtp, na.rm = TRUE),
    Mean = mean(wtp, na.rm = TRUE),
    Median = median(wtp, na.rm = TRUE),
    SD = sd(wtp, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    attribute = factor(attribute, levels = c("PA", "PA Half", "PA Full", "HNV", "HNV Visible")),
    functional_form = factor(functional_form, levels = c("Quadratic", "Linear", "Log", "Box-Cox", "Log-Linear"))
  ) %>%
  arrange(attribute, functional_form)

kable(wtp_summary_table, digits = 2)

attribute	functional_form	Min	Max	Mean	Median	SD
PA	Quadratic	3.40	9.10	7.90	7.97	0.80
PA	Linear	16.05	16.05	16.05	16.05	0.00
PA	Log	0.70	191.12	6.11	3.47	7.46
PA	Box-Cox	6.14	9.62	8.34	8.46	0.52
PA	Log-Linear	8.00	24.36	8.46	8.24	0.64
PA Half	Quadratic	16.14	18.17	16.57	16.54	0.28
PA Half	Linear	25.76	25.76	25.76	25.76	0.00
PA Half	Log	1.31	359.01	11.48	6.53	14.01
PA Half	Box-Cox	10.50	21.27	17.04	17.38	1.65
PA Half	Log-Linear	-13.92	17.66	16.76	17.20	1.24
PA Full	Quadratic	11.38	20.44	18.53	18.65	1.27
PA Full	Linear	25.56	25.56	25.56	25.56	0.00
PA Full	Log	1.62	444.84	14.22	8.09	17.36
PA Full	Box-Cox	17.47	22.89	19.07	18.87	0.74
PA Full	Log-Linear	16.98	80.76	18.79	17.91	2.50
HNV	Quadratic	-2.20	12.00	9.93	10.39	1.72
HNV	Linear	17.46	17.46	17.46	17.46	0.00
HNV	Log	0.73	244.68	8.57	6.27	8.32
HNV	Box-Cox	4.85	26.60	9.39	9.10	2.09
HNV	Log-Linear	5.73	149.01	10.34	8.99	4.89
HNV Visible	Quadratic	0.46	15.63	13.42	13.90	1.83
HNV Visible	Linear	19.92	19.92	19.92	19.92	0.00
HNV Visible	Log	1.04	347.50	12.17	8.91	11.82
HNV Visible	Box-Cox	7.72	36.05	14.02	13.65	2.82
HNV Visible	Log-Linear	8.71	194.43	14.68	12.93	6.34

wtp_table <- database_wtp %>%
  distinct(RID, .keep_all = TRUE) %>%
  summarise(across(starts_with("wtp_"), ~ sum(.x, na.rm = TRUE))) %>%
  pivot_longer(
    cols = everything(),
    names_to = "variable",
    values_to = "total_wtp"
  ) %>%
  mutate(
    functional_form = str_extract(variable, "(log_linear|box_cox|quadratic|linear|log)$"),
    functional_form = recode(
      functional_form,
      "quadratic" = "Quadratic",
      "linear" = "Linear",
      "log" = "Log",
      "box_cox" = "Box-Cox",
      "log_linear" = "Log-Linear"
    ),
    attribute_raw = str_remove(variable, "_(log_linear|box_cox|quadratic|linear|log)$"),
    attribute = recode(
      attribute_raw,
      "wtp_pa" = "PA",
      "wtp_pa_half" = "PA Half",
      "wtp_pa_full" = "PA Full",
      "wtp_hnv" = "HNV",
      "wtp_hnv_vis" = "HNV Visible"
    )
  ) %>%
  select(attribute, functional_form, total_wtp) %>%
  pivot_wider(
    names_from = functional_form,
    values_from = total_wtp
  ) %>%
  mutate(
    attribute = factor(attribute, levels = c("PA", "PA Half", "PA Full", "HNV", "HNV Visible"))
  ) %>%
  arrange(attribute) %>%
  rowwise() %>%
  mutate(
    `Max % deviation` = 100 * (
      max(c(Quadratic, Linear, Log, `Box-Cox`, `Log-Linear`), na.rm = TRUE) -
      min(c(Quadratic, Linear, Log, `Box-Cox`, `Log-Linear`), na.rm = TRUE)
    ) / abs(min(c(Quadratic, Linear, Log, `Box-Cox`, `Log-Linear`), na.rm = TRUE))
  ) %>%
  ungroup()

kable(wtp_table, digits = 2)

attribute	Quadratic	Linear	Log	Box-Cox	Log-Linear	Max % deviation
PA	112386.7	228267.5	86918.4	118642.5	120374.2	162.62
PA Half	235720.9	366394.6	163275.5	242361.7	238441.2	124.40
PA Full	263596.8	363589.3	202306.9	271247.1	267342.8	79.72
HNV	141303.7	248420.4	121910.8	133611.7	147076.5	103.77
HNV Visible	190880.0	283423.0	173142.5	199414.5	208815.6	63.69

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\}\)
One scaling parameter with radius logged	\((1 + \delta \cdot log(\text{D}))\left( \beta_{\text{NC}} \cdot \text{NC} + \beta_{\text{NC\_sq}} \cdot \text{NC}^2 \right)\)	\(\{\delta\}\)

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

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

WTP functions for fixed radius

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

pa_half

pa_full

hnv

hnv_vis

WTP functions for fixed endowment

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"
)

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

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

pa_half

pa_full

hnv

hnv_vis

Exp model

WTP functions for fixed radius

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

WTP functions for fixed endowment

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

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
)

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

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

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

One parameter scaling model

plot_wtp_distance_lines_dist_inc(
  model = dist_inc,
  attr = "pa",
  fixed_endowments = c(1, 5, 10, 20, 40, 80, 120),
  distance_seq = seq(1, 100, by = 0.25),
  ylim = c(-10, 30)
)

plot_wtp_distance_lines_dist_inc(
  model = dist_inc,
  attr = "pa_full",
  fixed_endowments = c(1, 5, 10, 20, 40, 80, 120),
  distance_seq = seq(1, 100, by = 0.25),
  ylim = c(-10, 50)
)

plot_wtp_lines_dist_inc(
  model = dist_inc,
  attr = "pa",
  endowment_seq = seq(1, 100, by = 1),
  distance_values = c(1, 5, 10, 20, 30, 50),
  median_line = 11,
  ylim = c(-10, 50)
)

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

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

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")

plot_wtp_income(linear_inc, "pa_half")

plot_wtp_income(linear_inc, "pa_full")

plot_wtp_income(linear_inc, "hnv")

plot_wtp_income(linear_inc, "hnv_visible")

Effect on income only on cost

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

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))

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

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

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

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