#' create_weight_grid
#' @description
#' function to create valid weight sets in long data.frame that can be
#' iterated through
#' @param x
#' @param wt_vals
#'
#' @returns
#' @export
#'
#' @examples
#' params_set_gte_1984 <- c(
#' "era5_land_soil_moisture_1m",
#' "cumu_era5_land_total_precipitation_sum",
#' "vhi",
#' "era5_land_snow_cover", # toying w/ adding this or not
#' "asi"
#' )
#' weight_combos <- create_weight_grid(params_set_gte_1984,wt_vals = c(0,seq(0.1, 1, by = 0.1)))
create_weight_grid <- function(x,wt_vals){
df <- map(x,\(xt){
tibble(
!!sym(xt) :=wt_vals
)
}
)|>
list_cbind()
df_expanded <- expand_grid(!!!df)
df_valid_combinations <- df_expanded[rowSums(df_expanded) == 1, ]
df_valid_combinations |>
mutate(
wt_id = row_number()
) |>
pivot_longer(
-wt_id, names_to = "parameter", values_to ="weight"
)
}
#' Title
#'
#' @param df
#' @param params_included
#' @param earliest_year
#'
#' @returns
#' @export
#'
#' @examples
#' df_env_model |>
#' normalize_to_z(params_included = l_params$gte1984,earliest_year = 1984)
normalize_to_z <- function(df,params_included, earliest_year){
df |>
filter(
year(yr_season)>=earliest_year
) |>
mutate(
pub_mo_label = as.character(pub_mo_label)
) |>
filter(
parameter %in% params_included,
pub_mo_label %in% c(month.abb[4:6])
) |>
group_by(
pub_mo_label, adm1_name, parameter
) |>
mutate(
zscore = scale(value,center=T,scale=T)[,1],
zscore = ifelse(parameter != "asi",zscore*-1,zscore)
) |>
ungroup()
}
#' Title
#'
#' @param df
#'
#' @returns
#' @export
#'
#' @examples
#' df_env_model |>
#' normalize_to_z(params_included = l_params$gte2000,earliest_year = 2001) |>
#' extract_truth_set()
extract_truth_set <- function(df){
df |>
filter(
month(pub_mo_date)==6,
parameter == "asi"
) |>
select(
yr_season, adm1_name, value,zscore_asi_Jun=zscore
) |>
distinct()
}
#' Title
#'
#' @param df
#' @param params_included
#'
#' @returns
#' @export
#'
#' @examples
summarise_z <- function(
df,
params_included,
weight_values
){
df_weight_grid <- create_weight_grid(x = params_included, wt_vals = weight_values)
df_filt <- df |>
filter(
parameter %in% params_included
)
split(df_filt,df_filt$pub_mo_label) |>
map(\(dft){
split(
df_weight_grid,
df_weight_grid$wt_id
) |>
map(
\(dft_w){
dft_weighted <- dft |>
left_join(
dft_w, by = "parameter"
)
dft_summarised <- dft_weighted |>
group_by(
yr_season ,
pub_mo_label,
adm1_name,
pub_mo_date,wt_id
) |>
summarise(
zscore = weighted.mean(zscore,w=weight,na.rm=T)
,.groups="drop"
) |>
mutate(
wt_set = list(dft_w)
)
}
) |>
list_rbind()
}
) |>
list_rbind()
}
#' Title
#'
#' @param df
#' @param params_included
#' @param earliest_year
#' @param rp
#'
#' @returns
#' @export
#'
#' @examples
weighted_classify <- function(df,
params_included,
earliest_year,
rp=3,
weight_values ){
df_normalized <- df |>
normalize_to_z(
params_included = params_included,
earliest_year = earliest_year
)
df_truth <- df_normalized |>
extract_truth_set()
df_z_weighted <- df_normalized |>
summarise_z(params_included = params_included,weight_values= weight_values)
df_z_weighted |>
left_join(
df_truth |>
select(-value)
) |>
utils$threshold_var(
var= "zscore",
by = c("pub_mo_label","adm1_name","wt_id"),
rp_threshold = rp
) |>
utils$threshold_var(
var= "zscore_asi_Jun",
by = c("pub_mo_label","adm1_name","wt_id"),
rp_threshold = rp
) |>
arrange(adm1_name, yr_season,wt_id) |>
ungroup()
}
single_indicator_performance <- function(df, parameter,earliest_year,rp){
df_normalized <- df |>
normalize_to_z(
params_included = parameter,
earliest_year = earliest_year
)
df_truth <- df_normalized |>
extract_truth_set()
df_classified <- df_normalized |>
left_join(
df_truth |>
select(-value)
) |>
utils$threshold_var(
var= "zscore",
by = c("pub_mo_label","adm1_name"),
rp_threshold = rp
) |>
utils$threshold_var(
var= "zscore_asi_Jun",
by = c("pub_mo_label","adm1_name"),
rp_threshold = rp
) |>
arrange(adm1_name, yr_season) |>
ungroup()
summarise_performance(df_classified,by = c("pub_mo_label","adm1_name"))
}
#' Title
#'
#' @param df
#' @param by
#'
#' @returns
#' @export
#'
#' @examples
summarise_performance <- function(
df,
by=c("pub_mo_label","adm1_name","wt_id","wt_set"),
parameter_subset = NULL
){
if(!is.null(parameter_subset)){
df |>
filter(
parameter %in% parameter_subset
)
}
df |>
mutate(
across(ends_with("_flag"),\(x) factor(x,levels = c("TRUE","FALSE")))
) |>
group_by(
across({{by}})
) |>
f_meas(zscore_asi_Jun_flag, zscore_flag, estimator= "binary",event_level = "first") |>
ungroup()
}
#' Title
#'
#' @param df
#'
#' @returns
#' @export
#'
#' @examples
top_performance_per_moment <- function(df,n=1){
# df= ldf_perf_all_models$gte1984,n=1
df_max <- df |>
group_by(pub_mo_label, wt_id) |>
summarise(
avg_estimate = mean(.estimate)
) |>
slice_max(
order_by = avg_estimate,
n= n
)|>
ungroup()
inner_join(
df_max,
select(df,
any_of(c("adm1_name","pub_mo_label", "wt_id", "wt_set",".estimate","asi_f1"))
),
by = c("pub_mo_label","wt_id")
)
}
#' Title
#'
#' @param df
#'
#' @returns
#' @export
#'
#' @examples
plot_optimal_compositions <- function(df, label_plot = F, pal){
df_labelled <- df|>
select(
pub_mo_label,adm1_name,wt_set,avg_estimate,.estimate
) |>
unnest(wt_set) |>
utils$label_parameters() |>
mutate(
p_label = glue(
"id: {wt_id}
f1: {scales::label_number(accuracy =0.001)(.estimate)}
avg f1: {scales::label_number(accuracy =0.001)(avg_estimate)}"
),
weight_pct_label = scales::percent(weight,accuracy = 1, trim = FALSE),
)
p <- df_labelled |>
group_by(pub_mo_label, adm1_name) |>
mutate(
pub_mo_facet = factor(pub_mo_label,levels= c("Apr","May","Jun")),
id = dense_rank(wt_id)
) |>
ungroup() |>
ggplot(
aes(x= id, y= weight,fill = parameter_label)
)+
geom_bar(
stat= "identity", color = "black"
)+
scale_fill_manual(values=pal)+
facet_grid(cols= vars(pub_mo_facet),
rows = vars(adm1_name)
,scales="free")+
labs(
x = "Different indicator weightings"
)+
scale_y_continuous(labels=scales::label_percent())+
theme(
legend.title = element_blank(),
axis.text.x = element_blank()
)
if(label_plot){
# df <- ldf_tops2$gte1984
df |>
mutate(
p_label = glue("{wt_id}: {scales::label_number(accuracy =0.01)(avg_estimate)}")
)
p <- p +
# composition plot
geom_text(
aes(x= id, y= weight, label = weight_pct_label),
position = position_stack(vjust = 0.5), color ="black"
)+
# wt plot
geom_label(aes(x= id, y= 1.1, label = p_label), color ="black",fill="beige",alpha=0.4)+
# theme(
# panel.spacing.y = unit(4, "lines")
# )+
expand_limits(y = 1.2)
}
p
}
#' Title
#' helper func to find weight ids that are simple and optimal
#' not a perfect solution, but gets us 90% there.
#' @param df
#'
#' @returns
#' @export
#'
#' @examples
plot_low_var_optimal <- function(df,
pal,
label_wt_id =T){
dfp <- df|>
select(
pub_mo_label, wt_set
) |>
unnest(wt_set) |>
group_by(
pub_mo_label,wt_id
) |>
mutate(
sd = sd(weight)
) |>
group_by(pub_mo_label) |>
slice_min(sd) |>
group_by(pub_mo_label) |>
mutate(
pub_mo_facet = factor(pub_mo_label,levels= c("Apr","May","Jun")),
id = dense_rank(wt_id)
) |>
ungroup() |>
utils$label_parameters()
p <- dfp |>
distinct() |>
mutate(
weight_pct_label = scales::percent(weight,accuracy = 1, trim = FALSE),
weight_id_label = glue("weight id: {wt_id}")
) |>
filter(pub_mo_label!= "Jun") |>
ggplot(
aes(
x= id,
y= weight,
fill = parameter_label
),
position = position_stack(vjust = 0.5)
)+
geom_bar(
stat= "identity", color = "black"
)+
geom_text(aes(x= id, y= weight, label = weight_pct_label),
position = position_stack(vjust = 0.5), color ="black")+
scale_fill_manual(values=pal) +
scale_y_continuous(labels =scales::label_percent())+
facet_grid(cols= vars(pub_mo_facet)
,scales="free")+
labs(
# title = "May publication weight combos that give avg 0.90 f1"
)
if(label_wt_id){
p <- p +
geom_text(
aes(
x= id+.25,
y= 0.5,
label = weight_id_label
), color ="black"
)
}
p
}
#' Title
#'
#' @param df
#' @param param_simple
#'
#' @returns
#' @export
#'
#' @examples
get_simple_weight_id <- function(df,param_simple){
df |>
select(wt_set) |>
unnest(wt_set) |>
distinct() |>
filter(
parameter == param_simple ,
weight == 1
) |>
pull(wt_id) |>
unique()
}
#' Title
#' @description
#' After plotting different optimal weight combinations you might want to compare specific
#' weighting compositions to a simpler model (i.e ASI). So this allows you to plug in the
#' weight ids for both april and may as well as `param_simple` which has always been ASI up to
#' this point in analysis. The output data.frame is interesting for plotting.
#'
#' @param df
#' @param wt_id_apr
#' @param wt_id_may
#' @param param_simple
#'
#' @returns
#' @export
#'
#' @examples
compare_to_simple_model <- function(df,wt_id_apr, wt_id_may,param_simple){
simple_model_id <-get_simple_weight_id(df= df, param_simple = param_simple)
df_chosen <- df |>
summarise_performance() |>
filter(
(pub_mo_label == "Apr" & wt_id==wt_id_apr)|
(pub_mo_label == "May" & wt_id==wt_id_may)|
wt_id == simple_model_id
) |>
filter(pub_mo_label != "Jun")
df_simple <- df_chosen |>
filter(
wt_id == simple_model_id
)
anti_join(df_chosen,df_simple) |>
left_join(
df_simple |>
select(
pub_mo_label,adm1_name,
estimate_simple =.estimate),
by = c("pub_mo_label","adm1_name")
)
}
