Open main paper
Packages for this section
library(tidyverse)
library(latex2exp)
library(jsonlite)
library(partykit) # to predict the evtree from dictionary6 Integrated framework
Objectives
This section combines the premium spectrum, the local metrics, and the policyholder partitions into one monitoring table per scenario. The Excel workbooks summarize, in dollars and in people, where proxy vulnerability and commercial loading concentrate.
Data for this section
the_CAS_colors_full <- c('#F89708', '#FDB922', '#FED35D', '#95CDF8', '#205ED5', '#142345')
marg_dist <- readRDS('preds/the_empirical_proportions.rds')
dictionnary_leaves_trees <- readRDS('evtree/dictionnary_leaves_trees.rds')
group_grid_stats <- jsonlite::fromJSON('preds/group_grid_stats.json')
group_pop_stats <- jsonlite::fromJSON('preds/group_pop_stats.json')Functions for this section (long)
# Helper function to calculate summaries
calculate_summary <- function(column, variable_name, summary_functions) {
summary_df <- data.frame(Variable = character(), Statistic = character(), Value = numeric())
for (func in summary_functions) {
if (grepl("^\\d*\\.\\d+$", func)) { # Quantile (VaR)
quantile_value <- as.numeric(func)
value <- quantile(column, probs = quantile_value, na.rm = TRUE)
label <- paste0("VaR ", func) # Clear label
summary_df <- rbind(summary_df, data.frame(Variable = variable_name, Statistic = label, Value = value))
} else if (grepl("^TVaR \\d*\\.\\d+$", func)) { # Tail Value at Risk (TVaR)
quantile_value <- as.numeric(sub("TVaR ", "", func)) # Extract quantile
var_value <- quantile(column, probs = quantile_value, na.rm = TRUE)
tvar_value <- mean(column[column >= var_value], na.rm = TRUE)
label <- paste0("TVaR ", quantile_value)
summary_df <- rbind(summary_df, data.frame(Variable = variable_name, Statistic = label, Value = tvar_value))
} else if (func == "mean") { # Mean
value <- mean(column, na.rm = TRUE)
label <- "Mean"
summary_df <- rbind(summary_df, data.frame(Variable = variable_name, Statistic = label, Value = value))
} else if (func == "std") { # Standard deviation
value <- sd(column, na.rm = TRUE) # Sample standard deviation
label <- "Std"
summary_df <- rbind(summary_df, data.frame(Variable = variable_name, Statistic = label, Value = value))
} else if (grepl("^%\\d+$", func)) { # Percentage of specific value
target_value <- as.numeric(sub("%", "", func))
value <- sum(column == target_value, na.rm = TRUE) / sum(!is.na(column))
label <- paste0("% of ", target_value)
summary_df <- rbind(summary_df, data.frame(Variable = variable_name, Statistic = label, Value = value))
} else {
warning(paste("Unsupported summary function:", func, "for variable", variable_name))
}
}
return(summary_df)
}
summarize_demographics <- function(data, sensitive, variables, summary_functions = c('0.1', 'mean', '0.9', '%0')) {
# Validate input
if (!is.data.frame(data)) stop("Input `data` must be a data frame.")
if (!all(variables %in% colnames(data))) stop("Some specified variables are not in the dataset.")
if(sensitive %in% variables) variables <- setdiff(variables, sensitive)
# Initialize the summary output
summary_table <- data.frame(Variable = 'Observations', Statistic = 'Count', Value = nrow(data), stringsAsFactors = FALSE)
# Loop over specified variables
for (var in c(sensitive, variables)) {
if (is.numeric(data[[var]])) {
numeric_summary <- calculate_summary(data[[var]], var, summary_functions)
summary_table <- rbind(summary_table, numeric_summary)
} else if (is.factor(data[[var]]) || is.character(data[[var]])) {
# Categorical: Append (c) to variable name
categorical_var <- paste0(var, " (c)")
if(var == sensitive) categorical_var <- paste0('(Sensitive) ', categorical_var)
# Count occurrences for all levels
counts <- table(data[[var]], useNA = "no") # Exclude NA from counts
levels_percent <- counts / sum(counts) # Proportion
# Create summary table
category_summary <- data.frame(
Variable = rep(categorical_var, length(counts)),
Statistic = paste0("lvl % : ", names(counts)), # Clearer level labels
Value = as.numeric(levels_percent) # Proportions, no *100
)
if(length(names(counts)) == 2 & '1' %in% names(counts)){
category_summary <- category_summary[which(names(counts) == '1'),]
}
summary_table <- rbind(summary_table, category_summary)
} else {
# Skip unsupported types
warning(paste("Variable", var, "is not numeric or categorical. Skipping."))
}
}
rownames(summary_table) <- NULL
return(summary_table)
}
# Load required packages
if (!require("statmod")) install.packages("statmod", quiet = TRUE)
library(statmod)
summarize_loss_premium <- function(data, loss_col, premium_col, distribution = "normal") {
# Validate input
if (!is.data.frame(data)) stop("Input `data` must be a data frame.")
if (!(loss_col %in% colnames(data))) stop("Loss column not found in data.")
if (!(premium_col %in% colnames(data))) stop("Premium column not found in data.")
# Filter non-NA values
valid_data <- data[!is.na(data[[loss_col]]) & !is.na(data[[premium_col]]), ]
if (nrow(valid_data) == 0) stop("No valid data after removing NA values.")
# Initialize summary table
summary_table <- data.frame(Variable = character(), Statistic = character(), Value = numeric(), stringsAsFactors = FALSE)
# Add summaries for loss and premium
loss_summary <- calculate_summary(valid_data[[loss_col]], "Loss", summary_functions = c('%0', 'mean'))
severity_summary <- calculate_summary(valid_data[[loss_col]][valid_data[[loss_col]] > 0], "Severity", summary_functions = c('mean', '0.95', 'TVaR 0.95'))
premium_summary <- calculate_summary(valid_data[[premium_col]], "Premium", summary_functions = c('0.05', 'mean', '0.95'))
# Add MAE (Mean Absolute Error)
mae <- mean(abs(valid_data[[loss_col]] - valid_data[[premium_col]]), na.rm = TRUE)
perf_summary_temp <- data.frame(Variable = "Performance metrics of price", Statistic = "MAE", Value = mae)
# Calculate deviance
deviance <- switch(
tolower(distribution),
"normal" = {
value <- mean((valid_data[[loss_col]] - valid_data[[premium_col]])^2)
label <- "Normal"
list(value = value, label = label)
},
"poisson" = {
if (any(valid_data[[premium_col]] <= 0)) stop("Premiums must be > 0 for Poisson.")
value <- 2 * sum(valid_data[[loss_col]] * log(valid_data[[loss_col]] / valid_data[[premium_col]]) -
(valid_data[[loss_col]] - valid_data[[premium_col]]), na.rm = TRUE)
label <- "Poisson"
list(value = value, label = label)
},
"binary" = {
if (any(valid_data[[loss_col]] < 0 | valid_data[[loss_col]] > 1)) stop("Losses must be in [0,1] for binary.")
value <- -2 * sum(
valid_data[[loss_col]] * log(valid_data[[premium_col]]) +
(1 - valid_data[[loss_col]]) * log(1 - valid_data[[premium_col]]),
na.rm = TRUE
)
label <- "Binary"
list(value = value, label = label)
},
"tweedie" = {
if (!requireNamespace("statmod", quietly = TRUE)) stop("Install the 'statmod' package for Tweedie distribution.")
p <- statmod::tweedie.profile(valid_data[[loss_col]] ~ valid_data[[premium_col]])$p.max
value <- sum(
statmod::tweedie.dev(valid_data[[loss_col]], valid_data[[premium_col]], p = p),
na.rm = TRUE
)
label <- "Tweedie"
list(value = value, label = label)
},
stop("Unsupported distribution type. Choose from 'normal', 'poisson', 'binary', or 'tweedie'.")
)
perf_summary <- rbind(perf_summary_temp, data.frame(Variable = "Performance metrics of price",
Statistic = paste0('Dev. ', deviance$label),
Value = deviance$value))
# Combine all summaries
summary_table <- rbind(loss_summary, severity_summary, premium_summary, perf_summary)
rownames(summary_table) <- NULL
return(summary_table)
}
# Load required package for Wasserstein distance
if (!requireNamespace("transport", quietly = TRUE)) install.packages("transport", quiet = TRUE)
library(transport)
summarize_loss_n_price_per_group <- function(data, loss_col, premium_col, protected_group_col, pdx_col,
levels = c(0, 1)) {
# Validate input
required_cols <- c(loss_col, premium_col, protected_group_col, pdx_col)
if (!all(required_cols %in% colnames(data))) stop("Some required columns are not in the dataset.")
if (length(levels) != 2) stop("The 'levels' argument must specify exactly two levels for the protected group.")
# Filter non-NA values
valid_data <- data[complete.cases(data[, required_cols]), ]
if (nrow(valid_data) == 0) stop("No valid data after removing NA values.")
# Check which specified levels are present
present_levels <- levels[levels %in% unique(valid_data[[protected_group_col]])]
# Initialize the summary output
summary_table <- data.frame(Variable = character(), Statistic = character(), Value = numeric(), stringsAsFactors = FALSE)
# Premium statistics mean
for (group in levels) {
if (group %in% present_levels) {
group_data <- valid_data[valid_data[[protected_group_col]] == group, ]
premium_data <- group_data[[premium_col]]
# Premium summaries
premium_summary <- calculate_summary(premium_data,
paste0("Mean price"),
summary_functions = 'mean')
premium_summary$Statistic <- paste0('D=', group)
summary_table <- rbind(summary_table, premium_summary)
} else {
# Add NA rows for missing group
na_summary <- data.frame(Variable = "Mean price", Statistic = paste0('D=', group),
Value = NA, stringsAsFactors = FALSE)
summary_table <- rbind(summary_table, na_summary)
}
}
# Premium statistics Var0.95
for (group in levels) {
if (group %in% present_levels) {
group_data <- valid_data[valid_data[[protected_group_col]] == group, ]
premium_data <- group_data[[premium_col]]
loss_data <- group_data[[loss_col]]
# Premium summaries
premium_summary <- calculate_summary(premium_data,
paste0("VaR 0.95 of price"),
summary_functions = '0.95')
premium_summary$Statistic <- paste0('D=', group)
summary_table <- rbind(summary_table, premium_summary)
} else {
# Add NA rows for missing group
na_summary <- data.frame(Variable = "VaR 0.95 of price",
Statistic = paste0('D=', group),
Value = NA, stringsAsFactors = FALSE)
summary_table <- rbind(summary_table, na_summary)
}
}
# Loss Ratio
for (group in levels) {
if (group %in% present_levels) {
group_data <- valid_data[valid_data[[protected_group_col]] == group, ]
premium_data <- group_data[[premium_col]]
loss_data <- group_data[[loss_col]]
# Premium summaries
loss_ratio <- mean(loss_data, na.rm = TRUE) / mean(premium_data, na.rm = TRUE)
loss_ratio_summary <- data.frame(
Variable = "Loss Ratio",
Statistic = paste0("D=", group),
Value = loss_ratio
)
summary_table <- rbind(summary_table, loss_ratio_summary)
} else {
# Add NA rows for missing group
na_summary <- data.frame(Variable = "Mean price", Statistic = paste0('D=', group),
Value = NA, stringsAsFactors = FALSE)
summary_table <- rbind(summary_table, na_summary)
}
}
# Average propensity
propensity <- mean(ifelse(group_data[[protected_group_col]] == 1,
group_data[[pdx_col]],
1 - group_data[[pdx_col]]), na.rm = TRUE)
propensity_summary <- data.frame(
Variable = "P(D=1|X)",
Statistic = "Mean",
Value = propensity
)
summary_table <- rbind(summary_table, propensity_summary)
# Compute Wasserstein distance for the specified levels if both are present
if (all(levels %in% present_levels)) {
group_a <- valid_data[valid_data[[protected_group_col]] == levels[1], premium_col]
group_b <- valid_data[valid_data[[protected_group_col]] == levels[2], premium_col]
# Check if both groups have sufficient observations
if (length(group_a) > 5 && length(group_b) > 5) {
wd <- wasserstein1d(group_a, group_b)
} else {
wd <- NA # Assign NA if one or both groups have insufficient data
}
} else {
# If one or both levels are missing, output NA for Wasserstein
wd <- NA
}
# Add Wasserstein Distance to the Summary Table
wd_summary <- data.frame(
Variable = "Wasserstein Distance",
Statistic = paste(levels[1], "vs", levels[2]),
Value = wd
)
summary_table <- rbind(summary_table, wd_summary)
rownames(summary_table) <- NULL
return(summary_table)
}
summarize_fairness_spectrum <- function(data, fairness_cols, spectrum_labels = NULL, summary_functions = c('mean')) {
# Validate input
if (!is.data.frame(data)) stop("Input `data` must be a data frame.")
if (!all(fairness_cols %in% colnames(data))) stop("Some specified fairness columns are not in the dataset.")
if (!is.null(spectrum_labels) && length(spectrum_labels) != length(fairness_cols)) {
stop("The length of 'spectrum_labels' must match the length of 'fairness_cols'.")
}
# Use column names as default labels if spectrum_labels is not provided
if (is.null(spectrum_labels)) spectrum_labels <- fairness_cols
# Initialize the summary output
summary_table <- data.frame(Variable = character(), Statistic = character(), Value = numeric(), stringsAsFactors = FALSE)
# Loop over specified fairness columns in the order provided
for (i in seq_along(fairness_cols)) {
fairness_col <- fairness_cols[i]
spectrum_label <- spectrum_labels[i]
# Ensure the column is numeric
if (!is.numeric(data[[fairness_col]])) {
warning(paste("Fairness column", fairness_col, "is not numeric. Skipping."))
next
}
# Calculate summaries and append to the summary table
fairness_summary <- calculate_summary(data[[fairness_col]], spectrum_label, summary_functions)
summary_table <- rbind(summary_table, fairness_summary)
}
rownames(summary_table) <- NULL
return(summary_table)
}
summarize_local_measures <- function(data,
split_measures = 'proxy_vuln',
whole_pop_measures = c('risk_spread', 'fair_range', 'parity_cost'),
protected_group_col,
summary_functions = c('mean', 'TVaR 0.95')) {
# Validate input
required_cols <- c(split_measures, whole_pop_measures, protected_group_col)
if (!all(required_cols %in% colnames(data))) stop("Some required columns are not in the dataset.")
# Filter non-NA values
valid_data <- data[complete.cases(data[, required_cols]), ]
if (nrow(valid_data) == 0) stop("No valid data after removing NA values.")
# Initialize the summary output
summary_table <- data.frame(Variable = character(), Statistic = character(), Value = numeric(), stringsAsFactors = FALSE)
### split_measures
the_levels <- unique(valid_data[[protected_group_col]])
for (smry_fn in summary_functions) {
for (split_meas in split_measures) {
for (i in 1:(length(the_levels) + 1)) {
if(i == 1){
summary_all_data <- calculate_summary(valid_data[[split_meas]], 'temp', summary_functions = smry_fn)
summary_all_data$Variable <- paste0(summary_all_data$Statistic, ' of ', split_meas)
summary_all_data$Statistic <- 'All'
summary_table <- rbind(summary_table, summary_all_data)
} else {
j <- i - 1
partial_data <- valid_data[[split_meas]][valid_data[[protected_group_col]] == the_levels[j]]
summary_protected <- calculate_summary(partial_data, 'temp', summary_functions = smry_fn)
summary_protected$Variable <- paste0(summary_protected$Statistic, ' of ',split_meas)
summary_protected$Statistic <- paste0('D=', the_levels[j])
summary_table <- rbind(summary_table, summary_protected)
}
}
}
}
for (other_meas in whole_pop_measures) {
summary_all_data <- calculate_summary(valid_data[[other_meas]], other_meas, summary_functions = 'mean')
summary_table <- rbind(summary_table, summary_all_data)
}
rownames(summary_table) <- NULL
return(summary_table)
}
compute_sb0_sb1 <- function(data, mua_col, mub_col, d_col, pd) {
# Validate input
if (!is.data.frame(data)) stop("Input `data` must be a data frame.")
if (!(mua_col %in% colnames(data))) stop("mu_A column not found in the dataset.")
if (!(mub_col %in% colnames(data))) stop("mu_B column not found in the dataset.")
if (!(d_col %in% colnames(data))) stop("D column not found in the dataset.")
if (length(pd) != 2 || any(pd <= 0) || sum(pd) != 1) stop("PD must be a valid probability vector of length 2 summing to 1.")
# Extract the columns
mu_A <- data[[mua_col]]
mu_B <- data[[mub_col]]
D <- data[[d_col]]
# Compute mu_B0 and mu_B1
mu_B0 <- ifelse(D == 1,
(mu_A - pd[2] * mu_B) / pd[1], # Formula for mu_B0 when D = 1
mu_B) # mu_B0 = mu_B when D = 0
mu_B1 <- ifelse(D == 0,
(mu_A - pd[1] * mu_B) / pd[2], # Formula for mu_B1 when D = 0
mu_B) # mu_B1 = mu_B when D = 1
# Return the modified dataset with mu_B0 and mu_B1
return(
list(mu_B0, mu_B1)
)
}
implied_propensity <- function(data, premium_col, mua_col, mub_col, d_col, pd){
sb0sb1 <- compute_sb0_sb1(data, mua_col = mua_col, mub_col = mub_col, d_col = d_col, pd = pd)
return(
(data[[premium_col]] - sb0sb1[[1]])/(sb0sb1[[2]] - sb0sb1[[1]])
)
}
summarize_premium_spectrum <- function(data, premium_col,
mua_col,
mub_col,
group_col,
rank_col,
pd,
summary_functions = c('mean', 'TVaR 0.95')) {
# Validate input
required_cols <- c(premium_col, mua_col, mub_col, group_col, rank_col)
if (!all(required_cols %in% colnames(data))) stop("Some required columns are not in the dataset.")
# Compute mu_B0 and mu_B1
sb0sb1 <- compute_sb0_sb1(data, mua_col = mua_col, mub_col = mub_col, d_col = group_col, pd = pd)
mu_B0 <- sb0sb1[[1]]
mu_B1 <- sb0sb1[[2]]
# Compute Implied Propensity and Commercial Loading
implied_prop <- (data[[premium_col]] - mu_B0) / (mu_B1 - mu_B0)
commercial_loading <- data[[premium_col]] - data[[mua_col]]
# Initialize summary output
summary_table <- data.frame(Variable = character(), Statistic = character(), Value = numeric(), stringsAsFactors = FALSE)
the_levels <- unique(data[[group_col]])
## commercial loading
for (smry_fn in summary_functions) {
for (i in 1:(length(the_levels) + 1)) {
if(i == 1){
summary_all_data <- calculate_summary(commercial_loading, 'temp', summary_functions = smry_fn)
summary_all_data$Variable <- paste0(summary_all_data$Statistic, ' of comm. loading')
summary_all_data$Statistic <- 'All'
summary_table <- rbind(summary_table, summary_all_data)
} else {
j <- i - 1
partial_data <- commercial_loading[data[[group_col]] == the_levels[j]]
summary_protected <- calculate_summary(partial_data, 'temp', summary_functions = smry_fn)
summary_protected$Variable <- paste0(summary_protected$Statistic, ' of comm. loading')
summary_protected$Statistic <- paste0('D=', the_levels[j])
summary_table <- rbind(summary_table, summary_protected)
}
}
}
## propensity + rank
for (i in 1:(length(the_levels) + 1)) {
if(i == 1){
summary_all_data <- calculate_summary(implied_prop, 'temp', summary_functions = 'mean')
summary_all_data$Variable <- paste0(summary_all_data$Statistic, ' implied propensity')
summary_all_data$Statistic <- 'All'
summary_table <- rbind(summary_table, summary_all_data)
} else {
j <- i - 1
partial_data <- implied_prop[data[[group_col]] == the_levels[j]]
summary_protected <- calculate_summary(partial_data, 'temp', summary_functions = 'mean')
summary_protected$Variable <- paste0(summary_protected$Statistic, ' implied propensity')
summary_protected$Statistic <- paste0('D=', the_levels[j])
summary_table <- rbind(summary_table, summary_protected)
}
}
for (i in 1:(length(the_levels) + 1)) {
if(i == 1){
summary_all_data <- calculate_summary(data[[rank_col]], 'temp', summary_functions = 'mean')
summary_all_data$Variable <- paste0(summary_all_data$Statistic, ' rank among spectrum')
summary_all_data$Statistic <- 'All'
summary_table <- rbind(summary_table, summary_all_data)
} else {
j <- i - 1
partial_data <- data[[rank_col]][data[[group_col]] == the_levels[j]]
summary_protected <- calculate_summary(partial_data, 'temp', summary_functions = 'mean')
summary_protected$Variable <- paste0(summary_protected$Statistic, ' rank among spectrum')
summary_protected$Statistic <- paste0('D=', the_levels[j])
summary_table <- rbind(summary_table, summary_protected)
}
}
rownames(summary_table) <- NULL
return(summary_table)
}
summarize_vulnerability <- function(data, premium_col, mua_col, group_col, thresholds = NULL) {
# Validate input
required_cols <- c(premium_col, mua_col, group_col)
if (!all(required_cols %in% colnames(data))) stop("Some required columns are not in the dataset.")
# Compute implied commercial adjustment relative to mu_A
commercial_adjustment <- (data[[premium_col]] - data[[mua_col]]) / data[[mua_col]]
# Initialize summary output
summary_table <- data.frame(Variable = character(), Statistic = character(), Value = numeric(), stringsAsFactors = FALSE)
if(is.null(thresholds)){
thresholds <- round(c(0, quantile(commercial_adjustment, c(0.65, 0.75))), 3)
}
the_levels <- unique(data[[group_col]])
# Calculate % of people exceeding each threshold
for (threshold in thresholds) {
for (i in 1:(length(the_levels) + 1)) {
if(i == 1){
percentage <- mean(commercial_adjustment > threshold, na.rm = TRUE)
vulnerability_summary <- data.frame(
Variable = paste0("% of ppl. with (Comm. load)/(Aware)>", threshold * 100, "%"),
Statistic = paste("All"),
Value = percentage
)
summary_table <- rbind(summary_table, vulnerability_summary)
} else {
j <- i - 1
partial_data <- commercial_adjustment[data[[group_col]] == the_levels[j]]
percentage <- mean(partial_data > threshold, na.rm = TRUE)
vulnerability_summary <- data.frame(
Variable = paste0("% of ppl. with (Comm. load)/(Aware)>", threshold * 100, "%"),
Statistic = paste0('D=', the_levels[j]),
Value = percentage
)
summary_table <- rbind(summary_table, vulnerability_summary)
}
}
}
rownames(summary_table) <- NULL
return(summary_table)
}
summarize_fairness_table <- function(data, group_col, marg_dist,
sensitive_variable,
top_bottom_summary = c(3, 3),
demographics_variables = c("D", "X2", "X1")) {
# Validate input
if (!is.data.frame(data)) stop("Input `data` must be a data frame.")
if (!all(group_col %in% c("1", colnames(data))))
stop("All group_col values must be '1' or valid column names in the dataset.")
if (!sensitive_variable %in% colnames(data))
stop("The `sensitive_variable` must be a valid column name in the dataset.")
# Preprocess demographic variables to limit levels to max 4 (including "Other")
preprocess_factors <- function(data, variables) {
for (var in variables) {
if (var %in% colnames(data) && (is.factor(data[[var]]) || is.character(data[[var]]))) {
counts <- table(data[[var]], useNA = "no")
counts <- sort(counts, decreasing = TRUE)
if (length(counts) > 4) {
top_3_levels <- names(head(counts, 3))
data[[var]] <- factor(
ifelse(data[[var]] %in% top_3_levels, data[[var]], "Other"),
levels = c(top_3_levels, "Other")
)
} else {
data[[var]] <- factor(data[[var]], levels = names(counts))
}
}
}
return(data)
}
# Apply preprocessing
data <- preprocess_factors(data, demographics_variables)
# Fetch levels of the sensitive variable
sensitive_levels <- levels(factor(data[[sensitive_variable]]))
# Helper function to tag sections
add_section <- function(summary_table, section_name) {
summary_table$Section <- section_name
return(summary_table)
}
# List of summary functions
summary_functions <- list(
"Demographics" = function(subdata) summarize_demographics(subdata,
sensitive = sensitive_variable,
variables = demographics_variables,
summary_functions = c('mean', 'std')),
"Loss and commercial price" = function(subdata) summarize_loss_premium(subdata,
loss_col = "Y",
premium_col = "prem",
distribution = 'normal'),
"Loss and price per group" = function(subdata) summarize_loss_n_price_per_group(
subdata, loss_col = "Y", premium_col = "prem",
protected_group_col = sensitive_variable, pdx_col = "pdx",
levels = sensitive_levels
),
"Local Measures" = function(subdata) summarize_local_measures(
subdata, split_measures = c("proxy_vuln"),
whole_pop_measures = c("risk_spread", "fair_range", "parity_cost"),
protected_group_col = sensitive_variable
),
"Fairness Spectrum" = function(subdata) summarize_fairness_spectrum(subdata,
fairness_cols = c("mu_B", "mu_U", "mu_A", "mu_H", "mu_C"),
summary_functions = c('mean'),
spectrum_labels = c('Best-est.', 'Unaware', 'Aware',
'Hyperaware', 'Corrective')),
"Fairness implications of commercial price" = function(subdata) summarize_premium_spectrum(
subdata, premium_col = "prem", mua_col = "mu_A", mub_col = "mu_B", group_col = sensitive_variable, rank_col = "r", pd = marg_dist
),
"Vulnerability" = function(subdata) summarize_vulnerability(
subdata, premium_col = "prem", mua_col = "mu_A",
group_col = sensitive_variable,
thresholds = round(c(0, quantile((data$prem - data$mu_A)/data$mu_A, c(0.65, 0.75))), 3) )
)
# Initialize list to hold all summaries
all_summaries <- list()
# edit group_col. build multiple tables: full, full-middle, pruned
#dictionnary_leaves_trees
#dictionnary_leaves_trees_comm
the_filters <- list('1' = rep(1, nrow(data)),
'proxy' = dictionnary_proxy)
for (gr in setdiff(group_col, '1')) {
data[['g']] <- as.numeric(as.character(data[[gr]]))
data <- left_join(data, dictionnary_proxy$rpart$dict,
by = c('g' = 'pred'))
data[[paste0(gr, '_node_new')]] <- data$node_new
data[[paste0(gr, '_node_or')]] <- data$node_or
data$g <- NULL; data$node_new <- NULL; data$node_or <- NULL;
}
# Loop over each variable in group_col
for (col in group_col) {
if (col == "1") {
subdata <- data
summary_table <- do.call(
rbind,
lapply(names(summary_functions), function(section) {
add_section(summary_functions[[section]](subdata), section)
})
)
summary_table <- summary_table %>%
rename_with(~ paste0("All Data"), starts_with("Value"))
all_summaries[["All Data"]] <- summary_table
} else {
levels <- levels(factor(data[[col]]))
for (level in levels) {
subdata <- data[data[[col]] == level, ]
summary_table <- do.call(
rbind,
lapply(names(summary_functions), function(section) {
add_section(summary_functions[[section]](subdata), section)
})
)
summary_table <- summary_table %>%
rename_with(~ paste0(names(group_col)[which(group_col == col)], " - ", level), starts_with("Value"))
key <- paste0(col, " - ", level)
all_summaries[[key]] <- summary_table
}
}
}
# Combine summaries while preserving section order
section_order <- names(summary_functions)
final_table <- Reduce(function(x, y) dplyr::left_join(x, y, by = c("Section", "Variable", "Statistic")), all_summaries) %>%
dplyr::mutate(Section = factor(Section, levels = section_order)) %>% # Enforce order
dplyr::arrange(Section) # Sort by Section
# Ensure column order
ordered_col_names <- c("Section", "Variable", "Statistic", setdiff(names(final_table), c("Section", "Variable", "Statistic")))
final_table <- final_table[, ordered_col_names]
return(final_table)
}
summarize_fairness_groups <- function(data,
group_col,
marg_dist,
sensitive_variable,
dictionnary_proxy,
dictionnary_comm,
demographics_variables = c("D", "X2", "X1")) {
# Define a helper function to process each group type
process_group <- function(data, group_type, dict_proxy, dict_comm) {
# Map proxy vulnerable groups
data[[paste0("pg_", group_type)]] <- factor(
predict(dict_proxy$model, newdata = data, type = 'node') %>% unname,
levels = dict_proxy$dict$node_or,
labels = dict_proxy$dict$node_new
)
# Regroup proxy groups if needed
if (group_type == "regrouped") {
levels(data[[paste0("pg_", group_type)]]) <- recode_bottom_top_middle(levels(data[[paste0("pg_", group_type)]]))
}
# Map commercially loaded groups
data[[paste0("cg_", group_type)]] <- factor(
predict(dict_comm$model, newdata = data, type = 'node') %>% unname,
levels = dict_comm$dict$node_or,
labels = dict_comm$dict$node_new
)
# Regroup commercial groups if needed
if (group_type == "regrouped") {
levels(data[[paste0("cg_", group_type)]]) <- recode_bottom_top_middle(levels(data[[paste0("cg_", group_type)]]))
}
# Generate the fairness summary table
group_col = c(
"All data" = "1",
"Proxy vulnerable groups" = paste0("pg_", group_type),
"Commercially loaded groups" = paste0("cg_", group_type)
)
big_tab <- summarize_fairness_table(
data = data,
group_col = group_col,
sensitive_variable = sensitive_variable,
marg_dist = marg_dist
)
# Add a row for "Leaf composition"
to_add_to_big_tab_full <- rep(NA, ncol(big_tab))
names(to_add_to_big_tab_full) <- names(big_tab)
# Extract the node labels for proxy and commercial groups
the_nms_p <- names(big_tab)[grepl(names(group_col[2]), names(big_tab))]
the_nms_c <- names(big_tab)[grepl(names(group_col[3]), names(big_tab))]
the_nodes_new_p <- str_split(the_nms_p, "-", n = 2) %>%
lapply(function(pieces) pieces[2]) %>%
unlist()
the_nodes_new_c <- str_split(the_nms_c, "-", n = 2) %>%
lapply(function(pieces) pieces[2]) %>%
unlist()
the_nodes_or_p <- sapply(as.numeric(the_nodes_new_p), function(node_new_id) {
if (is.na(node_new_id)) return(NA)
dict_proxy$dict %>% filter(node_new == node_new_id) %>% pull(node_or)
}) %>% unlist()
the_nodes_or_c <- sapply(as.numeric(the_nodes_new_c), function(node_new_id) {
if (is.na(node_new_id)) return(NA)
dict_comm$dict %>% filter(node_new == node_new_id) %>% pull(node_or)
}) %>% unlist()
# Simplify paths for nodes
the_nodes_path_simp_p <- sapply(the_nodes_or_p, function(node_or) {
if (is.na(node_or)) return(NA)
dict_proxy$paths$simp_concat[[as.character(node_or)]]
})
the_nodes_path_simp_c <- sapply(the_nodes_or_c, function(node_or) {
if (is.na(node_or)) return(NA)
dict_comm$paths$simp_concat[[as.character(node_or)]]
})
# Add paths to the final row of the summary table
to_add_to_big_tab_full[the_nms_p] <- the_nodes_path_simp_p
to_add_to_big_tab_full[the_nms_c] <- the_nodes_path_simp_c
return(list('table' = big_tab,
'paths' = to_add_to_big_tab_full))
}
# Process each group type
suppressWarnings({
suppressMessages({
big_tab_f_full <-process_group(data, "full",
dictionnary_proxy$rpart,
dictionnary_comm$rpart)
big_tab_r_full <- process_group(data, "regrouped",
dictionnary_proxy$rpart,
dictionnary_comm$rpart)
big_tab_p_full <- process_group(data, "pruned",
dictionnary_proxy$rpart_prune,
dictionnary_comm$rpart_prune)
})
})
# Return all three summary tables in a list
return(list(
"All subgroups" = big_tab_f_full,
"Top-Bottom subgroups" = big_tab_r_full,
"Coarse groups" = big_tab_p_full
))
}
library(openxlsx)
export_to_excel <- function(summary_tables,
file_name = "summary.xlsx") {
# Create a new workbook
wb <- createWorkbook()
for (i in seq_along(summary_tables)) {
summary_table <- summary_tables[[i]]$table
sheet_name <- names(summary_tables)[i]
to_add <- summary_tables[[i]]$paths
# Ensure sheet name is a valid character
sheet_name <- enc2native(as.character(sheet_name))
# Ensure all column names in the workbook are properly encoded
colnames(summary_table) <- enc2native(as.character(colnames(summary_table)))
# Add a worksheet
addWorksheet(wb, sheet_name)
# Extract column names and split into groups and labels
original_colnames <- colnames(summary_table)
split_names <- strsplit(original_colnames, " - ", fixed = TRUE)
groups <- sapply(split_names, function(x) ifelse(length(x) > 1, x[1], "")) # Left part
labels <- sapply(split_names, function(x) ifelse(length(x) > 1, x[2], x[1])) # Right part
# Write the second header row (labels)
writeData(wb, sheet_name, t(labels), startRow = 2, colNames = FALSE)
# Write the first header row (groups) and merge cells for groups
for (group in unique(groups)) {
if (group != "") {
matching_cols <- which(groups == group) # Columns in the same group
mergeCells(wb, sheet_name, cols = matching_cols, rows = 1) # Merge cells for the group
writeData(wb, sheet_name, group, startCol = min(matching_cols),
startRow = 1, colNames = FALSE)
}
}
# Write the main data below the headers
writeData(wb, sheet_name, summary_table, startRow = 3, colNames = FALSE)
# Apply a new theme for the headers
header_style_1 <- createStyle(
fontSize = 14, fontColour = the_CAS_colors_full[6], fgFill = the_CAS_colors_full[4],
halign = "CENTER", textDecoration = "Bold", border = "Bottom", borderColour = "black", borderStyle = "thick"
)
header_style_2 <- createStyle(
fontSize = 12, fontColour = the_CAS_colors_full[6], fgFill = the_CAS_colors_full[4],
halign = "CENTER", textDecoration = "Bold"
)
addStyle(wb, sheet_name, style = header_style_1, rows = 1, cols = 1:ncol(summary_table), gridExpand = TRUE, stack = TRUE)
addStyle(wb, sheet_name, style = header_style_2, rows = 2, cols = 1:ncol(summary_table), gridExpand = TRUE, stack = TRUE)
# Merge vertically consecutive cells in the first two columns
for (col in 1:2) { # Only for Section and Variable columns
values <- summary_table[[col]]
start_row <- 3 # Data starts at row 3
for (i in seq_along(values)) {
if (i == 1 || values[i] != values[i - 1]) {
start_merge <- start_row + i - 1
}
if (i == length(values) || values[i] != values[i + 1]) {
end_merge <- start_row + i - 1
mergeCells(wb, sheet_name, cols = col, rows = start_merge:end_merge)
}
}
}
# Apply vertical alignment and rotation styles to the first two columns
vertical_align_style <- createStyle(halign = "CENTER", valign = "CENTER",
wrapText = TRUE)
rotated_style <- createStyle(
textRotation = 90,
halign = "CENTER",
valign = "CENTER",
wrapText = TRUE
)
# Apply rotation to the first column (Section)
addStyle(wb, sheet_name, style = rotated_style, rows = 3:(nrow(summary_table) + 2), cols = 1, gridExpand = TRUE, stack = TRUE)
# Apply vertical alignment to the second column (Variable)
addStyle(wb, sheet_name, style = vertical_align_style, rows = 3:(nrow(summary_table) + 2), cols = 2, gridExpand = TRUE, stack = TRUE)
# Add alternating row colors based on blocks of the 'Variable' column
variable_groups <- rle(as.character(summary_table$Variable)) # Group consecutive 'Variable' values
start_row <- 3 # Data starts at row 3
group_starts <- cumsum(c(start_row, head(variable_groups$lengths, -1)))
group_ends <- group_starts + variable_groups$lengths - 1
# Apply alternating styles to blocks
for (i in seq_along(group_starts)) {
group_rows <- group_starts[i]:group_ends[i]
if (i %% 2 == 0) { # Alternate colors
addStyle(
wb, sheet_name,
style = createStyle(fgFill = "grey85"), # Light gray for even groups
rows = group_rows,
cols = 3:ncol(summary_table),
gridExpand = TRUE
)
} else {
addStyle(
wb, sheet_name,
style = createStyle(fgFill = "white"), # Light gray for even groups
rows = group_rows,
cols = 3:ncol(summary_table),
gridExpand = TRUE
)
}
}
# Add vertical borders for column groups
for (group in unique(groups)) {
if (group != "") {
matching_cols <- which(groups == group)
left_border_style <- createStyle(border = "Left", borderColour = "black", borderStyle = "thick")
right_border_style <- createStyle(border = "Right", borderColour = "black", borderStyle = "thick")
addStyle(wb, sheet_name, style = left_border_style, rows = 1:(nrow(summary_table) + 2), cols = min(matching_cols), gridExpand = TRUE, stack = TRUE)
addStyle(wb, sheet_name, style = right_border_style, rows = 1:(nrow(summary_table) + 2), cols = max(matching_cols), gridExpand = TRUE, stack = TRUE)
}
}
# Add bold horizontal lines after groups of rows (based on the `Section` column)
unique_sections <- unique(summary_table$Section)
for (section in unique_sections) {
group_rows <- which(summary_table$Section == section) + 2 # Offset for header rows
last_row <- max(group_rows)
horizontal_line_style <- createStyle(border = "Bottom", borderColour = "black", borderStyle = "thick")
addStyle(wb, sheet_name, style = horizontal_line_style, rows = last_row, cols = 1:ncol(summary_table), gridExpand = TRUE, stack = TRUE)
}
# Apply number formatting for numeric columns (round to 2 decimal places for display)
numeric_cols <- which(sapply(summary_table, is.numeric)) # Identify numeric columns
if (length(numeric_cols) > 0) {
number_format_style <- createStyle(numFmt = "0.00")
addStyle(wb, sheet_name, style = number_format_style, rows = 3:(nrow(summary_table) + 2), cols = numeric_cols, gridExpand = TRUE, stack = TRUE)
}
# Apply 2-digit number formatting for 'Fairness Spectrum' rows
pct_rows <- which(grepl('%', summary_table$Variable) | grepl('%', summary_table$Statistic)) + 2
integer_rows <- which(grepl('Observations', summary_table$Variable)) + 2
digits3_rows <- which(grepl('Loss Ratio', summary_table$Variable) | grepl('Wasserstein', summary_table$Variable)) + 2
# Define styles for formatting
pct_style <- createStyle(numFmt = "0.00%") # Percentage with 2 decimal places
integer_style <- createStyle(numFmt = "0") # Integer format
digits3_style <- createStyle(numFmt = "0.000") # Decimal with 3 digits
# Apply percentage format
if (length(pct_rows) > 0) {
addStyle(
wb,
sheet = sheet_name,
style = pct_style,
rows = pct_rows,
cols = 4:ncol(summary_table), # Apply to numeric columns
gridExpand = TRUE,
stack = TRUE
)
}
# Apply integer format
if (length(integer_rows) > 0) {
addStyle(
wb,
sheet = sheet_name,
style = integer_style,
rows = integer_rows,
cols = 4:ncol(summary_table), # Apply to numeric columns
gridExpand = TRUE,
stack = TRUE
)
}
# Apply 3-decimal format
if (length(digits3_rows) > 0) {
addStyle(
wb,
sheet = sheet_name,
style = digits3_style,
rows = digits3_rows,
cols = 4:ncol(summary_table), # Apply to numeric columns
gridExpand = TRUE,
stack = TRUE
)
}
loss_rows <- which(summary_table$Variable == 'Loss Ratio')
premium_rows <- which(summary_table$Variable == 'Mean price')
proxy_rows <- which(summary_table$Variable == 'Mean of proxy_vuln')
cload_rows <- which(summary_table$Variable == 'Mean of comm. loading')
matched_rows <- c(loss_rows, premium_rows, proxy_rows, cload_rows) + 2
cols <- 5:(ncol(summary_table) + 1)
# Apply conditional formatting for each matched row
for (row in matched_rows) {
# Verify that cols and rows are numeric and valid
if (!is.numeric(cols) || !is.numeric(row)) {
stop("Columns and rows must be numeric.")
}
conditionalFormatting(
wb,
sheet = sheet_name,
cols = cols, rows = row,
type = 'colourScale',
style = c("#91CF60", "#FFFFBF", "#FC8D59")
)
}
# Auto-adjust column widths and set a smaller width for the first column
setColWidths(wb, sheet_name, cols = 1, widths = 5) # Smaller width for the Section column
setColWidths(wb, sheet_name, cols = 2, widths = 15) # Smaller width for the Section column
setColWidths(wb, sheet_name, cols = 3, widths = 8) # Smaller width for the Section column
setColWidths(wb, sheet_name, cols = 4:ncol(summary_table), widths = 8)
## Add custom content
# Write the custom row directly to the sheet
start_row <- nrow(summary_table) + 3
writeData(wb, sheet_name, t(to_add %>% unname),
startRow = start_row,
startCol = 1, colNames = FALSE)
# Create a style for the custom row
custom_style <- createStyle(
fgFill = "grey90", # Light grey background
fontColour = "grey50", # Dark grey text
wrapText = TRUE, # Enable text wrapping
valign = "TOP", # Center vertical alignment,
halign = 'LEFT',
textDecoration = "Bold"
)
# Apply the style to the custom row
addStyle(
wb, sheet = sheet_name, style = custom_style,
rows = start_row, cols = 1:(length(to_add)),
gridExpand = TRUE, stack = TRUE
)
# Adjust the row height for the new row (taller row for long text)
setRowHeights(wb, sheet_name, rows = start_row, heights = 60)
# Freeze the first 3 columns and the first 2 rows
freezePane(wb, sheet_name, firstActiveRow = 5, firstActiveCol = 4)
}
# Save the workbook to a file
saveWorkbook(wb, file_name, overwrite = TRUE)
message("Summary table successfully exported to: ", file_name)
}
library(dplyr)
# Generate a gradient of 5 colors between two given colors
generate_gradient <- function(start_color, end_color, n = 5) {
scales::gradient_n_pal(c(start_color, end_color))(seq(0, 1, length.out = n))
}
recode_bottom_top_middle <- function(levels, numeric = FALSE, n_bottom = 3, n_top = 3, factor_values = NULL) {
# Check if input is valid
if (is.null(levels) || length(levels) < 1) {
stop("Input levels must be a non-empty numeric vector.")
}
# Ensure levels are numeric
levels <- as.numeric(levels)
# If numeric = TRUE, factor_values must be provided
if (numeric && (is.null(factor_values) || length(factor_values) < 1)) {
stop("For numeric = TRUE, a valid 'factor_values' vector must be provided.")
}
# Handle cases where the total number of levels is insufficient for grouping
if (length(levels) <= (n_bottom + n_top + 1)) {
return(factor(levels, levels = as.character(levels))) # No middle group possible
}
# Define the bottom and top ranges dynamically based on the original order
bottom <- levels[1:n_bottom] # First `n_bottom` values
top <- levels[(length(levels) - n_top + 1):length(levels)] # Last `n_top` values
middle <- levels[!(levels %in% c(bottom, top))] # Remaining middle
# Dynamically create label for the middle group
if (length(middle) > 1) {
if (numeric) {
# Subset the factor_values for rows where the levels belong to the middle group
middle_values <- factor_values[factor_values %in% middle]
middle_mean <- round(mean(as.numeric(as.character(middle_values)), na.rm = TRUE), 3) # Compute mean
middle_label <- middle_mean
} else {
middle_label <- paste0(min(middle), "-", max(middle)) # Range as label
}
} else if (length(middle) == 1) {
middle_label <- as.character(middle) # Single value retains its original label
} else {
middle_label <- NULL
}
# Create a mapping for novel levels
recoded <- sapply(levels, function(x) {
if (x %in% bottom) return(as.character(x)) # Keep bottom levels as is
if (x %in% middle) return(as.character(middle_label)) # Group middle levels dynamically
if (x %in% top) return(as.character(x)) # Keep top levels as is
})
# Return as a factor with ordered levels maintaining the original order
return(factor(recoded, levels = unique(recoded), ordered = TRUE))
}Code to export the summary tables
if (!dir.exists('tables')) dir.create('tables')
for (scenario in names(group_pop_stats)) {
data <- group_pop_stats[[scenario]]$valid
# Extract dictionaries
dictionnary_proxy <- dictionnary_leaves_trees$proxy_vuln[[scenario]]
dictionnary_comm <- dictionnary_leaves_trees$comm_load[[scenario]]
export_to_excel(summary_tables = summarize_fairness_groups(
data = data %>% mutate(D = factor(D), X2 = factor(X2)),
marg_dist = marg_dist[[scenario]],
sensitive_variable = "D",
dictionnary_proxy = dictionnary_proxy,
dictionnary_comm = dictionnary_comm),
file_name = paste0('../docs/tables/summaries_', scenario, '.xlsx'))
}6.1 Download Disparity Tables
Note
You can download the Excel disparity tables for each scenario below.