From 032eaa80cdd36be62fd3c3204c4afd3f44b21d88 Mon Sep 17 00:00:00 2001 From: Nicholas Clark Date: Fri, 6 Oct 2023 21:25:18 +1000 Subject: [PATCH] first gp() additions --- DESCRIPTION | 4 +- NAMESPACE | 5 + R/compute_edf.R | 7 +- R/dynamic.R | 18 + R/get_linear_predictors.R | 108 ++++-- R/get_mvgam_priors.R | 29 +- R/globals.R | 3 +- R/gp.R | 683 ++++++++++++++++++++++++++++++++++ R/marginaleffects.mvgam.R | 60 +++ R/mvgam.R | 81 +++- R/stan_utils.R | 129 ++++++- R/summary.mvgam.R | 224 +++++------ R/trends.R | 60 --- R/update_priors.R | 4 +- man/plot_predictions.mvgam.Rd | 107 ++++++ src/mvgam.dll | Bin 1051136 -> 1051136 bytes tests/testthat/Rplots.pdf | Bin 22864 -> 22882 bytes 17 files changed, 1316 insertions(+), 206 deletions(-) create mode 100644 R/gp.R create mode 100644 man/plot_predictions.mvgam.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 3cb84a50..c96db513 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: mvgam Title: Multivariate (Dynamic) Generalized Additive Models -Version: 1.0.5 -Date: 2023-06-03 +Version: 1.0.6 +Date: 2023-10-06 Authors@R: person("Nicholas J", "Clark", , "nicholas.j.clark1214@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-7131-3301")) diff --git a/NAMESPACE b/NAMESPACE index 3c8ad740..d5a027d1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -78,6 +78,7 @@ export(plot_mvgam_series) export(plot_mvgam_smooth) export(plot_mvgam_trend) export(plot_mvgam_uncertainty) +export(plot_predictions.mvgam) export(ppc) export(roll_eval_mvgam) export(score) @@ -131,6 +132,7 @@ importFrom(magrittr,"%>%") importFrom(marginaleffects,get_coef) importFrom(marginaleffects,get_predict) importFrom(marginaleffects,get_vcov) +importFrom(marginaleffects,plot_predictions) importFrom(marginaleffects,set_coef) importFrom(methods,new) importFrom(mgcv,bam) @@ -150,6 +152,7 @@ importFrom(posterior,as_draws_rvars) importFrom(posterior,rhat) importFrom(posterior,variables) importFrom(rlang,missing_arg) +importFrom(rlang,parse_expr) importFrom(rstantools,posterior_epred) importFrom(rstantools,posterior_linpred) importFrom(rstantools,posterior_predict) @@ -169,7 +172,9 @@ importFrom(stats,dlnorm) importFrom(stats,dnbinom) importFrom(stats,dnorm) importFrom(stats,dpois) +importFrom(stats,drop.terms) importFrom(stats,ecdf) +importFrom(stats,fitted) importFrom(stats,formula) importFrom(stats,frequency) importFrom(stats,gaussian) diff --git a/R/compute_edf.R b/R/compute_edf.R index 1bd7ae35..b802bc55 100644 --- a/R/compute_edf.R +++ b/R/compute_edf.R @@ -1,4 +1,5 @@ #' Compute approximate EDFs of smooths +#' @importFrom stats fitted #'@noRd compute_edf = function(mgcv_model, object, rho_names, sigma_raw_names){ @@ -91,7 +92,11 @@ compute_edf = function(mgcv_model, object, rho_names, sigma_raw_names){ (mgcv_model$off[i] + ncol(mgcv_model$S[[i]]) - 1) XWXS[ind, ind] <- XWXS[ind, ind] + mgcv_model$S[[i]] * lambda[i] } - edf <- diag(solve(XWXS, XWX)) + suppressWarnings(edf <- try(diag(solve(XWXS, XWX)), silent = TRUE)) + if(inherits(edf, 'try-error')){ + edf <- rep(1, length(coef(mgcv_model))) + names(edf) <- names(coef(mgcv_model)) + } mgcv_model$edf <- edf } diff --git a/R/dynamic.R b/R/dynamic.R index cd681b75..e0196341 100644 --- a/R/dynamic.R +++ b/R/dynamic.R @@ -145,6 +145,24 @@ interpret_mvgam = function(formula, N){ attr(newformula, '.Environment') <- attr(formula, '.Environment') + # Check if any terms use the gp wrapper, as mvgam cannot handle + # multivariate GPs yet + response <- terms.formula(newformula)[[2]] + tf <- terms.formula(newformula, specials = c("gp")) + which_gp <- attr(tf,"specials")$gp + if(length(which_gp) != 0L){ + gp_details <- vector(length = length(which_gp), + mode = 'list') + for(i in seq_along(which_gp)){ + gp_details[[i]] <- eval(parse(text = rownames(attr(tf, + "factors"))[which_gp[i]])) + } + if(any(unlist(lapply(purrr::map(gp_details, 'term'), length)) > 1)){ + stop('mvgam cannot yet handle multidimensional gps', + call. = FALSE) + } + } + # Check if any terms use the dynamic wrapper response <- terms.formula(newformula)[[2]] tf <- terms.formula(newformula, specials = c("dynamic")) diff --git a/R/get_linear_predictors.R b/R/get_linear_predictors.R index 248946bc..22d4bbf2 100644 --- a/R/get_linear_predictors.R +++ b/R/get_linear_predictors.R @@ -26,6 +26,47 @@ obs_Xp_matrix = function(newdata, mgcv_model){ type = 'lpmatrix')) } + # Check for any gp() terms and update the design matrix + # accordingly + if(!is.null(attr(mgcv_model, 'gp_att_table'))){ + # Compute the eigenfunctions from the supplied attribute table, + # and add them to the Xp matrix + + # Extract GP attributes + gp_att_table <- attr(mgcv_model, 'gp_att_table') + gp_covariates <- unlist(purrr::map(gp_att_table, 'covariate')) + by <- unlist(purrr::map(gp_att_table, 'by')) + level <- unlist(purrr::map(gp_att_table, 'level')) + k <- unlist(purrr::map(gp_att_table, 'k')) + scale <- unlist(purrr::map(gp_att_table, 'scale')) + mean <- unlist(purrr::map(gp_att_table, 'mean')) + max_dist <- unlist(purrr::map(gp_att_table, 'max_dist')) + boundary <- unlist(purrr::map(gp_att_table, 'boundary')) + L <- unlist(purrr::map(gp_att_table, 'L')) + + # Compute eigenfunctions + test_eigenfunctions <- lapply(seq_along(gp_covariates), function(x){ + prep_eigenfunctions(data = newdata, + covariate = gp_covariates[x], + by = by[x], + level = level[x], + k = k[x], + boundary = boundary[x], + L = L[x], + mean = mean[x], + scale = scale[x], + max_dist = max_dist[x]) + }) + + # Find indices to replace in the design matrix and replace with + # the computed eigenfunctions + starts <- purrr::map(gp_att_table, 'first_coef') + ends <- purrr::map(gp_att_table, 'last_coef') + for(i in seq_along(starts)){ + Xp[,c(starts[[i]]:ends[[i]])] <- test_eigenfunctions[[i]] + } + } + return(Xp) } @@ -46,31 +87,6 @@ trend_Xp_matrix = function(newdata, trend_map, series = 'all', trend_test$series <- trend_indicators trend_test$y <- NULL - # Because these are set up inherently as dynamic factor models, - # we ALWAYS need to forecast the full set of trends, regardless of - # which series (or set of series) is being forecast - # data.frame(series = trend_test$series, - # time = trend_test$time, - # row_num = 1:length(trend_test$time)) %>% - # dplyr::group_by(series, time) %>% - # dplyr::slice_head(n = 1) %>% - # dplyr::ungroup() %>% - # dplyr::arrange(time, series) %>% - # dplyr::pull(row_num) -> inds_keep - # - # if(inherits(newdata, 'list')){ - # trend_test <- lapply(trend_test, function(x){ - # if(is.matrix(x)){ - # matrix(x[inds_keep,], ncol = NCOL(x)) - # } else { - # x[inds_keep] - # } - # - # }) - # } else { - # trend_test <- trend_test[inds_keep, ] - # } - suppressWarnings(Xp_trend <- try(predict(mgcv_model, newdata = trend_test, type = 'lpmatrix'), @@ -95,6 +111,48 @@ trend_Xp_matrix = function(newdata, trend_map, series = 'all', newdata = testdat, type = 'lpmatrix')) } + + # Check for any gp() terms and update the design matrix + # accordingly + if(!is.null(attr(mgcv_model, 'gp_att_table'))){ + # Compute the eigenfunctions from the supplied attribute table, + # and add them to the Xp matrix + + # Extract GP attributes + gp_att_table <- attr(mgcv_model, 'gp_att_table') + gp_covariates <- unlist(purrr::map(gp_att_table, 'covariate')) + by <- unlist(purrr::map(gp_att_table, 'by')) + level <- unlist(purrr::map(gp_att_table, 'level')) + k <- unlist(purrr::map(gp_att_table, 'k')) + scale <- unlist(purrr::map(gp_att_table, 'scale')) + mean <- unlist(purrr::map(gp_att_table, 'mean')) + max_dist <- unlist(purrr::map(gp_att_table, 'max_dist')) + boundary <- unlist(purrr::map(gp_att_table, 'boundary')) + L <- unlist(purrr::map(gp_att_table, 'L')) + + # Compute eigenfunctions + test_eigenfunctions <- lapply(seq_along(gp_covariates), function(x){ + prep_eigenfunctions(data = newdata, + covariate = gp_covariates[x], + by = by[x], + level = level[x], + k = k[x], + boundary = boundary[x], + L = L[x], + mean = mean[x], + scale = scale[x], + max_dist = max_dist[x]) + }) + + # Find indices to replace in the design matrix and replace with + # the computed eigenfunctions + starts <- purrr::map(gp_att_table, 'first_coef') + ends <- purrr::map(gp_att_table, 'last_coef') + for(i in seq_along(starts)){ + Xp_trend[,c(starts[[i]]:ends[[i]])] <- test_eigenfunctions[[i]] + } + } + return(Xp_trend) } diff --git a/R/get_mvgam_priors.R b/R/get_mvgam_priors.R index d1b06ea1..3118dc1e 100644 --- a/R/get_mvgam_priors.R +++ b/R/get_mvgam_priors.R @@ -223,6 +223,31 @@ get_mvgam_priors = function(formula, # Ensure series and time variables are present data_train <- validate_series_time(data_train, name = 'data') + # Check for gp and mo terms in the formula + orig_formula <- gp_terms <- mo_terms <- NULL + if(any(grepl('gp(', attr(terms(formula), 'term.labels'), fixed = TRUE))){ + + # Check that there are no multidimensional gp terms + formula <- interpret_mvgam(formula, N = max(data_train$time)) + orig_formula <- formula + + # Keep intercept? + keep_intercept <- attr(terms(formula), 'intercept') == 1 + + # Indices of gp() terms in formula + gp_terms <- which_are_gp(formula) + + # Extract GP attributes + gp_details <- get_gp_attributes(formula) + + # Replace with s() terms so the correct terms are included + # in the model.frame + formula <- gp_to_s(formula) + if(!keep_intercept){ + formula <- update(formula, trend_y ~ . -1) + } + } + # Check for missing rhs in formula drop_obs_intercept <- FALSE if(length(attr(terms(formula), 'term.labels')) == 0 & @@ -240,7 +265,9 @@ get_mvgam_priors = function(formula, call. = FALSE) } } - orig_formula <- formula + if(is.null(orig_formula)){ + orig_formula <- formula + } # Validate observation formula formula <- interpret_mvgam(formula, N = max(data_train$time)) diff --git a/R/globals.R b/R/globals.R index bba0c982..7960d0f1 100644 --- a/R/globals.R +++ b/R/globals.R @@ -8,4 +8,5 @@ utils::globalVariables(c("y", "year", "smooth_vals", "smooth_num", "mod_call", "particles", "obs", "mgcv_model", "param_name", "outcome", "mgcv_plottable", "term", "data_test", "object", "row_num", "trends_test", - "trend", "trend_series", "trend_y", ".", "gam")) + "trend", "trend_series", "trend_y", ".", "gam", + "group", "mod", "row_id")) diff --git a/R/gp.R b/R/gp.R new file mode 100644 index 00000000..edff9fc9 --- /dev/null +++ b/R/gp.R @@ -0,0 +1,683 @@ +#' Make gp() attributes table and necessary stan lines +#' @noRd +make_gp_additions = function(gp_details, data, + newdata, + model_data, mgcv_model, + gp_terms){ + # Need to expand combination of GPs if any of the by variables + # is a factor; mvgam will drop unused levels automatically + by <- gp_details$by + gp_details$row_id <- 1:NROW(gp_details) + gp_covariates <- gp_details$gp_covariates + gp_details_orig <- gp_details + if(any(!is.na(by))){ + for(i in 1:length(by)){ + if(!is.na(by[i])){ + if(is.factor(data[[by[i]]])){ + nlevels <- length(levels(droplevels(data[[by[i]]]))) + new_details <- do.call(rbind, lapply(1:nlevels, function(x){ + gp_details_orig[i,] + })) + new_details$level <- levels(droplevels(data[[by[i]]])) + + rows_drop <- which(gp_details$gp_covariates == gp_covariates[i] & + gp_details$by == by[i]) + + gp_details <- gp_details[-rows_drop,] + gp_details <- rbind(gp_details, new_details) + } + } + } + } + + # Preserve ordering of terms + gp_details %>% + dplyr::arrange(row_id, level) %>% + dplyr::select(-row_id) -> gp_details + + # Prepare the GP objects and Stan data lines; + # Data updates need to happen BEFORE calling vectorise + # model_file modifications can happen AFTER drop_obs_intercept + # to ensure everything behaves predictably + gp_covariates <- gp_details$gp_covariates + scale <- gp_details$scale + boundary <- gp_details$boundary + by <- gp_details$by + level <- gp_details$level + + # Prep the covariate for GP modelling + gp_data <- lapply(seq_along(gp_covariates), function(x){ + + # Find the correct k to ensure that the total number of coefficients + # when using gp(k = k) is the same as when using s(k = k + 1) + smooth_terms <- unlist(purrr::map(mgcv_model$smooth, 'term')) + smooth_bys <- unlist(purrr::map(mgcv_model$smooth, 'by')) + if(any(smooth_bys == 'NA')){ + smooth_bys[smooth_bys == 'NA'] <- NA + } + term_k <- mgcv_model$smooth[[min(which(smooth_bys %in% by[x] & + smooth_terms == gp_covariates[x]))]]$df + + prep_gp_covariate(data = data, + covariate = gp_covariates[x], + by = by[x], + level = level[x], + scale = scale[x], + k = term_k, + boundary = boundary[x]) + }) + + # Consolidate Stan data objects and add to model_data + gp_stan_data <- do.call(c, purrr::map(gp_data, 'data_append')) + model_data <- append(model_data, gp_stan_data) + + # Consolidate attribute tables + gp_att_table <- purrr::map(gp_data, 'att_table') + + # Create updated design matrix by replacing the s() basis functions with + # the gp() eigenfunctions + coefs_replace <- list() + for(x in gp_terms){ + label <- attr(terms(formula(mgcv_model)), 'term.labels')[x] + s_attributes <- eval(rlang::parse_expr(label)) + if(s_attributes$by != 'NA'){ + coef_name <- paste0('s(', s_attributes$term, '):', s_attributes$by) + } else { + coef_name <- paste0('s(', s_attributes$term, ')') + } + which_replace <- grep(coef_name, names(coef(mgcv_model)), fixed = TRUE) + names(mgcv_model$coefficients)[which_replace] <- + gsub('s(', 'gp(', names(mgcv_model$coefficients)[which_replace], + fixed = TRUE) + coefs_replace[[x]] <- which_replace + } + + newX <- model_data$X + + # Training data eigenfunctions + eigenfuncs <- do.call(cbind, purrr::map(gp_data, 'eigenfunctions')) + + # Testing data eigenfunctions + if(!is.null(newdata)){ + gp_covariates <- unlist(purrr::map(gp_att_table, 'covariate')) + by <- unlist(purrr::map(gp_att_table, 'by')) + level <- unlist(purrr::map(gp_att_table, 'level')) + k <- unlist(purrr::map(gp_att_table, 'k')) + scale <- unlist(purrr::map(gp_att_table, 'scale')) + mean <- unlist(purrr::map(gp_att_table, 'mean')) + max_dist <- unlist(purrr::map(gp_att_table, 'max_dist')) + boundary <- unlist(purrr::map(gp_att_table, 'boundary')) + L <- unlist(purrr::map(gp_att_table, 'L')) + test_eigenfunctions <- lapply(seq_along(gp_covariates), function(x){ + prep_eigenfunctions(data = newdata, + covariate = gp_covariates[x], + by = by[x], + level = level[x], + k = k[x], + boundary = boundary[x], + L = L[x], + mean = mean[x], + scale = scale[x], + max_dist = max_dist[x]) + }) + + eigenfuncs <- rbind(eigenfuncs, + do.call(cbind, test_eigenfunctions)) + } + + newX[unlist(coefs_replace), ] <- t(eigenfuncs) + model_data$X <- newX + + # Consolidate Stan data lines + gp_stan_lines <- paste0(purrr::map(gp_data, 'data_lines'), collapse = '') + + # Add coefficient indices to attribute table and to Stan data + for(covariate in seq_along(gp_att_table)){ + coef_indices <- grep(gp_att_table[[covariate]]$name, + names(coef(mgcv_model)), fixed = TRUE) + gp_att_table[[covariate]]$first_coef <- min(coef_indices) + gp_att_table[[covariate]]$last_coef <- max(coef_indices) + + gp_names <- gp_att_table[[covariate]]$name + gp_names <- gsub('(', '_', gp_names, fixed = TRUE) + gp_names <- gsub(')', '_', gp_names, fixed = TRUE) + gp_names <- gsub(':', 'by', gp_names, fixed = TRUE) + + gp_stan_lines <- paste0(gp_stan_lines, + paste0('array[', gp_att_table[[covariate]]$k, + '] int b_idx_', + gp_names, + '; // gp basis coefficient indices\n')) + gp_idx_data <- list(coef_indices) + names(gp_idx_data) <- paste0('b_idx_', + gp_names) + model_data <- append(model_data, gp_idx_data) + } + + # Add the GP attribute table to the mgcv_model + attr(mgcv_model, 'gp_att_table') <- gp_att_table + + return(list(model_data = model_data, + mgcv_model = mgcv_model, + gp_stan_lines = gp_stan_lines, + gp_att_table = gp_att_table)) +} + +#' Which terms are gp() terms? +#' @noRd +which_are_gp = function(formula){ + tf <- terms.formula(formula, specials = c("gp")) + if(is.null(rlang::f_lhs(formula))){ + out <- attr(tf,"specials")$gp + } else { + out <- attr(tf,"specials")$gp - 1 + } + return(out) +} + +#' Convert gp() terms to s() terms for initial model construction#' +#' @importFrom stats drop.terms +#' @noRd +gp_to_s <- function(formula){ + # Extract details of gp() terms + gp_details <- get_gp_attributes(formula) + + # Drop these terms from the formula + which_gp <- which_are_gp(formula) + response <- rlang::f_lhs(formula) + + suppressWarnings(tt <- try(drop.terms(terms(formula), + which_gp, + keep.response = TRUE), + silent = TRUE)) + if(inherits(tt, 'try-error')){ + newformula <- as.formula(paste(response, '~ 1')) + } else { + tt <- drop.terms(terms(formula), which_gp, keep.response = TRUE) + newformula <- reformulate(attr(tt, "term.labels"), rlang::f_lhs(formula)) + } + + # Now replace the gp() terms with s() for constructing the initial model + s_terms <- vector() + for(i in 1:NROW(gp_details)){ + if(!is.na(gp_details$by[i])){ + s_terms[i] <- paste0('s(', + gp_details$gp_covariates[i], + ', by = ', + gp_details$by[i], + ', k = ', + gp_details$k[i] + 1, ')') + } else { + s_terms[i] <- paste0('s(', + gp_details$gp_covariates[i], + ', k = ', + gp_details$k[i] + 1, ')') + } + } + + if(length(attr(terms(newformula), 'term.labels')) == 0){ + rhs <- '1' + } else { + rhs <- attr(terms(newformula), 'term.labels') + } + + newformula <- as.formula(paste(response, '~', + paste(paste(rhs, + collapse = '+'), '+', + paste(s_terms, collapse = '+')))) + attr(newformula, '.Environment') <- attr(formula, '.Environment') + return(newformula) +} + + +#' Store attributes of the gp terms +#' @importFrom rlang parse_expr +#' @noRd +get_gp_attributes = function(formula){ + gp_terms <- rownames(attr(terms(formula), 'factors'))[ + grep('gp', rownames(attr(terms(formula), 'factors')))] + gp_attributes <- lapply(seq_along(gp_terms), function(x){ + eval(rlang::parse_expr(gp_terms[x])) + }) + + # Extract information necessary to construct the GP terms + gp_covariates <- unlist(purrr::map(gp_attributes, 'term')) + k <- unlist(purrr::map(gp_attributes, 'k')) + if(any(is.na(k))){ + k[is.na(k)] <- 10 + } + scale <- unlist(purrr::map(gp_attributes, 'scale')) + boundary <- unlist(purrr::map(gp_attributes, 'c')) + if(any(is.na(boundary))){ + boundary[is.na(boundary)] <- 5.0/4 + } + by <- unlist(purrr::map(gp_attributes, 'by')) + if(any(by == 'NA')){ + by[by == 'NA'] <- NA + } + + # Return as a data.frame + return(data.frame(gp_covariates, + k, + scale, + boundary, + by, + level = NA)) +} + + +#' Propagate a Hilbert basis GP +#' Evaluate Laplacian eigenfunction for a given GP basis function +#' @noRd +phi = function(boundary, m, centred_covariate) { + 1 / sqrt(boundary) * sin((m * pi)/(2 * boundary) * + (centred_covariate + boundary)) +} + +#' Evaluate eigenvalues for a given GP basis function +#' @noRd +lambda = function(boundary, m) { + ((m * pi)/(2 * boundary))^2 +} + +#' Spectral density squared exponential Gaussian Process kernel +#' @noRd +spd = function(alpha_gp, rho_gp, eigenvalues) { + (alpha_gp^2) * sqrt(2 * pi) * rho_gp * + exp(-0.5 * (rho_gp^2) * (eigenvalues^2)) +} + +#' @noRd +sim_hilbert_gp = function(alpha_gp, + rho_gp, + b_gp, + last_trends, + fc_times, + train_times, + mean_train_times){ + + num_gp_basis <- length(b_gp) + + # Get vector of eigenvalues of covariance matrix + eigenvalues <- vector() + for(m in 1:num_gp_basis){ + eigenvalues[m] <- lambda(boundary = (5.0/4) * + (max(train_times) - min(train_times)), + m = m) + } + + # Get vector of eigenfunctions + eigenfunctions <- matrix(NA, nrow = length(fc_times), + ncol = num_gp_basis) + for(m in 1:num_gp_basis){ + eigenfunctions[, m] <- phi(boundary = (5.0/4) * + (max(train_times) - min(train_times)), + m = m, + centred_covariate = fc_times - mean_train_times) + } + + # Compute diagonal of covariance matrix + diag_SPD <- sqrt(spd(alpha_gp = alpha_gp, + rho_gp = rho_gp, + sqrt(eigenvalues))) + + # Compute GP trend forecast + as.vector((diag_SPD * b_gp) %*% t(eigenfunctions)) +} + +#' Mean-center and scale the particular covariate of interest +#' so that the maximum Euclidean distance between any two points is 1 +#' @noRd +scale_cov <- function(data, covariate, by, level, scale = TRUE, + mean, max_dist){ + if(!is.na(by)){ + if(!is.na(level)){ + Xgp <- data[[covariate]][data[[by]] == level] + } + } else { + Xgp <- data[[covariate]] + } + + if(is.na(mean)){ + Xgp_mean <- mean(Xgp, na.rm = TRUE) + } else { + Xgp_mean <- mean + } + + if(is.na(max_dist)){ + Xgp_max_dist <- (abs(max(Xgp, na.rm = TRUE) - + min(Xgp, na.rm = TRUE))) + } else { + Xgp_max_dist <- max_dist + } + + if(scale){ + # Mean center and divide by max euclidean distance + (Xgp - Xgp_mean) / Xgp_max_dist + + } else { + # Just mean center + Xgp - Xgp_mean + } +} + +#' prep GP eigenfunctions +#' @noRd +prep_eigenfunctions = function(data, + covariate, + by = NA, + level = NA, + k, + boundary, + mean = NA, + max_dist = NA, + scale = TRUE, + L){ + + # Extract and scale covariate (scale set to FALSE if this is a prediction + # step so that we can scale by the original training covariate values supplied + # in mean and max_dist) + covariate_cent <- scale_cov(data = data, + covariate = covariate, + by = by, + level = level, + mean = mean, + max_dist = max_dist, + scale = scale) + + # Construct matrix of eigenfunctions + eigenfunctions <- matrix(NA, nrow = length(covariate_cent), + ncol = k) + if(missing(L)){ + L <- brms:::choose_L(covariate_cent, boundary) + } + + for(m in 1:k){ + # eigenfunctions[, m] <- phi(boundary = boundary * + # (max(covariate_cent) - + # min(covariate_cent)), + # m = m, + # centred_covariate = covariate_cent) + + eigenfunctions[, m] <- brms:::eigen_fun_cov_exp_quad(x = matrix(covariate_cent), + m = m, + L = L) + } + + # Multiply eigenfunctions by the 'by' variable if one is supplied + if(!is.na(by)){ + if(!is.na(level)){ + # no multiplying needed as this is a factor by variable, + # but we need to pad the eigenfunctions with zeros + # for the observations where the by is a different level + full_eigens <- matrix(0, nrow = length(data[[by]]), + ncol = NCOL(eigenfunctions)) + full_eigens[(1:length(data[[by]]))[ + data[[by]] == level],] <- eigenfunctions + eigenfunctions <- full_eigens + } else { + eigenfunctions <- eigenfunctions * data[[by]] + } + } + eigenfunctions +} + +#' Prep Hilbert Basis GP covariates +#' @noRd +prep_gp_covariate = function(data, + covariate, + by = NA, + level = NA, + scale = TRUE, + boundary = 5.0/4, + k = 20){ + + covariate_cent <- scale_cov(data = data, + covariate = covariate, + scale = scale, + by = by, + mean = NA, + max_dist = NA, + level = level) + + if(!is.na(by)){ + if(!is.na(level)){ + Xgp <- data[[covariate]][data[[by]] == level] + } + } else { + Xgp <- data[[covariate]] + } + + covariate_mean <- mean(Xgp, na.rm = TRUE) + covariate_max_dist <- ifelse(scale, + abs(max(Xgp, + na.rm = TRUE) - + min(Xgp, + na.rm = TRUE)), + 1) + # Check k + if(k > length(unique(covariate_cent))){ + warning('argument "k" > number of unique covariate values;\ndropping to the maximum allowed "k"', + call. = FALSE) + k <- length(unique(covariate_cent)) + } + + # Construct vector of eigenvalues for GP covariance matrix; the + # same eigenvalues are always used in prediction, so we only need to + # create them when prepping the data. They will need to be included in + # the Stan data list + L <- brms:::choose_L(covariate_cent, boundary) + eigenvalues <- vector() + for(m in 1:k){ + eigenvalues[m] <- sqrt(lambda(boundary = L, + m = m)) + } + + # Construct matrix of eigenfunctions; this will change depending on the values + # of the covariate, so it needs to be computed and included as data but also needs + # to be computed to make predictions + eigenfunctions <- prep_eigenfunctions(data = data, + covariate = covariate, + by = by, + level = level, + k = k, + boundary = boundary, + mean = NA, + max_dist = covariate_max_dist, + scale = scale) + + # Make attributes table + byname <- ifelse(is.na(by), '', paste0(':', by)) + covariate_name <- paste0('gp(', covariate, ')', byname) + if(!is.na(level)){ + covariate_name <- paste0(covariate_name, level) + } + att_table <- list(effect = 'gp', + name = covariate_name, + covariate = covariate, + by = by, + level = level, + k = k, + boundary = boundary, + L = L, + scale = scale, + mean = covariate_mean, + max_dist = covariate_max_dist, + eigenvalues = eigenvalues) + + # Items to add to Stan data + # Number of basis functions + covariate_name <- gsub('(', '_', covariate_name, fixed = TRUE) + covariate_name <- gsub(')', '_', covariate_name, fixed = TRUE) + covariate_name <- gsub(':', 'by', covariate_name, fixed = TRUE) + data_lines <- paste0('int k_', covariate_name, '; // basis functions for approximate gp\n') + append_dat <- list(k = k) + names(append_dat) <- paste0('k_', covariate_name, '') + + # Approximate GP eigenvalues + data_lines <- paste0(data_lines, paste0( + 'vector[', + 'k_', covariate_name, + '] l_', covariate_name, '; // approximate gp eigenvalues\n'), + collapse = '\n') + append_dat2 <- list(slambda = eigenvalues) + names(append_dat2) <- paste0('l_', covariate_name, '') + append_dat <- append(append_dat, append_dat2) + + # Return necessary objects in a list + list(att_table = att_table, + data_lines = data_lines, + data_append = append_dat, + eigenfunctions = eigenfunctions) +} + +add_gp_model_file = function(model_file, model_data, mgcv_model, gp_additions){ + + # Add data lines + model_file[grep('int ytimes[n, n_series];', + model_file, fixed = TRUE)] <- + paste0(model_file[grep('int ytimes[n, n_series];', + model_file, fixed = TRUE)], + '\n', + gp_additions$gp_stan_lines) + model_file <- readLines(textConnection(model_file), n = -1) + + # Replace the multi_normal_prec lines with spd_cov_exp_quad + gp_names <- unlist(purrr::map(attr(mgcv_model, 'gp_att_table'), 'name')) + gp_names_clean <- gsub('(', '_', gp_names, fixed = TRUE) + gp_names_clean <- gsub(')', '_', gp_names_clean, fixed = TRUE) + gp_names_clean <- gsub(':', 'by', gp_names_clean, fixed = TRUE) + s_to_remove <- list() + for(i in seq_along(gp_names)){ + s_name <- gsub('gp(', 's(', gp_names[i], fixed = TRUE) + to_replace <- grep(paste0('// prior for ', s_name, '...'), + model_file, fixed = TRUE) + 1 + pattern <- "S\\s*(.*?)\\s*\\[" + result <- regmatches(model_file[to_replace], + regexec(pattern, model_file[to_replace]))[[1]] + s_to_remove[[i]] <- unique(unlist(regmatches(result, + gregexpr("[[:digit:]]+", result)))) + + model_file[grep(paste0('// prior for ', s_name, '...'), + model_file, fixed = TRUE)] <- + gsub('s(', 'gp(', model_file[grep(paste0('// prior for ', s_name, '...'), + model_file, fixed = TRUE)], fixed = TRUE) + + model_file[to_replace] <- + paste0('z_', + gp_names_clean[i], + ' ~ std_normal();\n', + 'alpha_', + gp_names_clean[i], + ' ~ normal(0, 0.5);\n', + 'rho_', + gp_names_clean[i], + ' ~ inv_gamma(1.65, 5.97);\n', + 'b_raw[b_idx_', + gp_names_clean[i], + '] ~ std_normal();\n') + } + b_line <- max(grep('b[', model_file, fixed = TRUE)) + b_edits <- paste0('b[b_idx_', + gp_names_clean, + '] = sqrt(spd_cov_exp_quad(l_', + gp_names_clean, + ', alpha_', + gp_names_clean, + ', rho_', + gp_names_clean, + ')) .* z_', + gp_names_clean, + ';', + collapse = '\n') + model_file[b_line] <- paste0(model_file[b_line], + '\n', + b_edits) + model_file <- readLines(textConnection(model_file), n = -1) + + # Remove un-needed penalty matrices from the model file and the + # model data + for(i in seq_along(unique(unlist(s_to_remove)))){ + model_data[[paste0('S', + unique(unlist(s_to_remove))[i])]] <- NULL + model_file <- model_file[-grep(paste0('mgcv smooth penalty matrix S', + unique(unlist(s_to_remove))[i]), + model_file, fixed = TRUE)] + } + + # Add alpha, rho and z lines in parameters and model blocks + alpha_names <- paste(paste0('real alpha_', gp_names_clean, + ';'), + collapse = '\n') + rho_names <- paste(paste0('real rho_', gp_names_clean, + ';'), + collapse = '\n') + z_names <- paste(paste0('vector[k_', + gp_names_clean, + '] z_', + gp_names_clean, + ';'), + collapse = '\n') + model_file[grep("vector[num_basis] b_raw;", model_file, fixed = TRUE)] <- + paste0("vector[num_basis] b_raw;\n\n", + '// gp term sd parameters\n', + alpha_names, + '\n\n// gp term length scale parameters\n', + rho_names, + '\n\n// gp term latent variables\n', + z_names, + '\n') + model_file <- readLines(textConnection(model_file), n = -1) + + # Add spd_cov_exp_quad function from brms code + if(any(grepl('functions {', model_file, fixed = TRUE))){ + model_file[grep('functions {', model_file, fixed = TRUE)] <- + paste0('functions {\n', + '/* Spectral density function of a Gaussian process\n', + '* with squared exponential covariance kernel\n', + '* Args:\n', + '* l_gp: numeric eigenvalues of an SPD GP\n', + '* alpha_gp: marginal SD parameter\n', + '* rho_gp: length-scale parameter\n', + '* Returns:\n', + '* numeric values of the GP function evaluated at l_gp\n', + '*/\n', + 'vector spd_cov_exp_quad(data vector l_gp, real alpha_gp, real rho_gp) {\n', + 'int NB = size(l_gp);\n', + 'vector[NB] out;\n', + 'real constant = square(alpha_gp) * (sqrt(2 * pi()) * rho_gp);\n', + 'real neg_half_lscale2 = -0.5 * square(rho_gp);\n', + 'for (m in 1:NB) {\n', + 'out[m] = constant * exp(neg_half_lscale2 * square(l_gp[m]));\n', + '}\n', + 'return out;\n', + '}\n') + } else { + model_file[grep('Stan model code', model_file)] <- + paste0('// Stan model code generated by package mvgam\n', + 'functions {\n', + '/* Spectral density function of a Gaussian process\n', + '* with squared exponential covariance kernel\n', + '* Args:\n', + '* l_gp: numeric eigenvalues of an SPD GP\n', + '* alpha_gp: marginal SD parameter\n', + '* rho_gp: length-scale parameter\n', + '* Returns:\n', + '* numeric values of the GP function evaluated at l_gp\n', + '*/\n', + 'vector spd_cov_exp_quad(data vector l_gp, real alpha_gp, real rho_gp) {\n', + 'int NB = size(l_gp);\n', + 'vector[NB] out;\n', + 'real constant = square(alpha_gp) * (sqrt(2 * pi()) * rho_gp);\n', + 'real neg_half_lscale2 = -0.5 * square(rho_gp);\n', + 'for (m in 1:NB) {\n', + 'out[m] = constant * exp(neg_half_lscale2 * square(l_gp[m]));\n', + '}\n', + 'return out;\n', + '}\n}\n') + } + model_file <- readLines(textConnection(model_file), n = -1) + + return(list(model_file = model_file, + model_data = model_data)) +} diff --git a/R/marginaleffects.mvgam.R b/R/marginaleffects.mvgam.R index 061f79e7..2fb734da 100644 --- a/R/marginaleffects.mvgam.R +++ b/R/marginaleffects.mvgam.R @@ -17,6 +17,66 @@ #' @author Nicholas J Clark NULL +#' Effect plot as implemented in \pkg{marginaleffects} +#' +#' Convenient way to call marginal or conditional effect plotting functions +#' implemented in the \pkg{marginaleffects} package +#' @importFrom marginaleffects plot_predictions +#' @importFrom bayesplot color_scheme_set color_scheme_get +#' @inheritParams marginaleffects::plot_predictions +#' @return A \code{\link[ggplot2:ggplot]{ggplot}} object +#' that can be further customized using the \pkg{ggplot2} package, +#' or a `data.frame` (if `draw=FALSE`) +#' +#' @export +plot_predictions.mvgam = function(model, + condition = NULL, + by = NULL, + newdata = NULL, + type = NULL, + vcov = NULL, + conf_level = 0.95, + wts = NULL, + transform = NULL, + points = 0, + rug = FALSE, + gray = FALSE, + draw = TRUE, + ...){ + # Set red colour scheme + orig_col <- .Options$ggplot2.discrete.colour + orig_fill <- .Options$ggplot2.discrete.fill + options(ggplot2.discrete.colour = c("#B97C7C", + "#A25050", + "#8F2727", + "darkred", + "#630000", + "#300000", + "#170000"), + ggplot2.discrete.fill = c("#B97C7C", + "#A25050", + "#8F2727", + "darkred", + "#630000", + "#300000", + "#170000")) + plot_predictions(model = model, + condition = condition, + by = by, + newdata = newdata, + type = type, + vcov = vcov, + conf_level = conf_level, + wts = wts, + transform = transform, + points = points, + rug = rug, + gray = gray, + draw = draw, + ...) + +} + #' Functions needed for working with marginaleffects #' @rdname mvgam_marginaleffects #' @export diff --git a/R/mvgam.R b/R/mvgam.R index 7e898317..3f5669b3 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -616,6 +616,31 @@ mvgam = function(formula, validate_pos_integer(samples) validate_pos_integer(thin) + # Check for gp and mo terms in the formula + orig_formula <- gp_terms <- mo_terms <- NULL + if(any(grepl('gp(', attr(terms(formula), 'term.labels'), fixed = TRUE))){ + + # Check that there are no multidimensional gp terms + formula <- interpret_mvgam(formula, N = max(data_train$time)) + orig_formula <- formula + + # Keep intercept? + keep_intercept <- attr(terms(formula), 'intercept') == 1 + + # Indices of gp() terms in formula + gp_terms <- which_are_gp(formula) + + # Extract GP attributes + gp_details <- get_gp_attributes(formula) + + # Replace with s() terms so the correct terms are included + # in the model.frame + formula <- gp_to_s(formula) + if(!keep_intercept){ + formula <- update(formula, trend_y ~ . -1) + } + } + # Check for missing rhs in formula drop_obs_intercept <- FALSE if(length(attr(terms(formula), 'term.labels')) == 0 & @@ -633,7 +658,9 @@ mvgam = function(formula, call. = FALSE) } } - orig_formula <- formula + if(is.null(orig_formula)){ + orig_formula <- formula + } # Check for brmspriors if(!missing("data")){ @@ -1475,6 +1502,27 @@ mvgam = function(formula, inits <- family_inits(family = family_char, trend_model, smooths_included, model_data) + # Include any GP term updates + if(!is.null(gp_terms)){ + gp_additions <- make_gp_additions(gp_details = gp_details, + data = data_train, + newdata = data_test, + model_data = stan_objects$model_data, + mgcv_model = ss_gam, + gp_terms = gp_terms) + stan_objects$model_data <- gp_additions$model_data + ss_gam <- gp_additions$mgcv_model + + gp_names <- unlist(purrr::map(gp_additions$gp_att_table, 'name')) + gp_names <- gsub('gp(', 's(', gp_names, fixed = TRUE) + rhos_change <- list() + for(i in seq_along(gp_names)){ + rhos_change[[i]] <- grep(gp_names[i], rho_names, fixed = TRUE) + } + rho_names[c(unique(unlist(rhos_change)))] <- gsub('s(', 'gp(', + rho_names[c(unique(unlist(rhos_change)))], + fixed = TRUE) + } # Vectorise likelihoods vectorised <- vectorise_stan_lik(model_file = stan_objects$stan_file, @@ -1585,6 +1633,33 @@ mvgam = function(formula, vectorised$model_file <- readLines(textConnection(vectorised$model_file), n = -1) } + # Remaining model file updates for any GP terms + if(!is.null(gp_terms)){ + final_gp_updates <- add_gp_model_file(model_file = vectorised$model_file, + model_data = vectorised$model_data, + mgcv_model = ss_gam, + gp_additions = gp_additions) + vectorised$model_file <- final_gp_updates$model_file + vectorised$model_data <- final_gp_updates$model_data + } + + # Update monitor pars for any GP terms + if(any(grepl('real alpha_gp', + vectorised$model_file, fixed = TRUE)) & + !lfo){ + alpha_params <- trimws(gsub(';', '', gsub('real ', + '', + grep('real alpha_gp', + vectorised$model_file, fixed = TRUE, + value = TRUE), fixed = TRUE))) + rho_params <- trimws(gsub(';', '', gsub('real ', + '', + grep('real rho_gp', + vectorised$model_file, fixed = TRUE, + value = TRUE), fixed = TRUE))) + param <- c(param, alpha_params, rho_params) + } + # Tidy the representation clean_up <- vector() for(x in 1:length(vectorised$model_file)){ @@ -1674,7 +1749,7 @@ mvgam = function(formula, # outside of mvgam #### if(!run_model){ unlink('base_gam.txt') - output <- structure(list(call = formula, + output <- structure(list(call = orig_formula, trend_call = if(!missing(trend_formula)){ trend_formula } else { @@ -2133,7 +2208,7 @@ mvgam = function(formula, } #### Return the output as class mvgam #### - output <- structure(list(call = formula, + output <- structure(list(call = orig_formula, trend_call = if(!missing(trend_formula)){ trend_formula } else { diff --git a/R/stan_utils.R b/R/stan_utils.R index 052b9e58..8b78c303 100644 --- a/R/stan_utils.R +++ b/R/stan_utils.R @@ -2570,6 +2570,9 @@ add_trend_predictors = function(trend_formula, b_trend_lines <- gsub('\\bb\\b', 'b_trend', b_trend_lines) b_trend_lines <- gsub('raw', 'raw_trend', b_trend_lines) b_trend_lines <- gsub('num_basis', 'num_basis_trend', b_trend_lines) + b_trend_lines <- gsub('idx', 'trend_idx', b_trend_lines) + b_trend_lines <- gsub('l_gp', 'trend_l_gp', b_trend_lines) + b_trend_lines <- gsub('k_gp', 'trend_k_gp', b_trend_lines) model_file[grep("// derived latent states", model_file, fixed = TRUE)] <- paste0('// process model basis coefficients\n', paste(b_trend_lines, collapse = '\n'), @@ -2589,10 +2592,11 @@ add_trend_predictors = function(trend_formula, # Add any multinormal smooth lines if(any(grepl('multi_normal_prec', trend_model_file)) | - any(grepl('// priors for smoothing parameters', trend_model_file))){ + any(grepl('// priors for smoothing parameters', trend_model_file)) | + any(grepl('// prior for gp', trend_model_file))){ trend_smooths_included <- TRUE - # Replace any noncontiguous indices from trend model so names aren't + # Replace any indices from trend model so names aren't # conflicting with any possible indices in the observation model if(any(grepl('idx', trend_model_file))){ trend_model_file <- gsub('idx', 'trend_idx', trend_model_file) @@ -2600,7 +2604,8 @@ add_trend_predictors = function(trend_formula, names(idx_data) <- gsub('idx', 'trend_idx', names(idx_data)) model_data <- append(model_data, idx_data) - idx_lines <- grep('int trend_idx', trend_model_file) + idx_lines <- c(grep('int trend_idx', trend_model_file), + grep('// gp basis coefficient indices', trend_model_file)) model_file[min(grep('data {', model_file, fixed = TRUE))] <- paste0('data {\n', paste(trend_model_file[idx_lines], @@ -2608,6 +2613,107 @@ add_trend_predictors = function(trend_formula, model_file <- readLines(textConnection(model_file), n = -1) } + # Check for gp() terms + if(any(grepl('l_gp', trend_model_file))){ + + # Add spd_cov_exp_quad function from brms code + if(any(grepl('functions {', model_file, fixed = TRUE))){ + model_file[grep('functions {', model_file, fixed = TRUE)] <- + paste0('functions {\n', + '/* Spectral density function of a Gaussian process\n', + '* with squared exponential covariance kernel\n', + '* Args:\n', + '* l_gp: numeric eigenvalues of an SPD GP\n', + '* alpha_gp: marginal SD parameter\n', + '* rho_gp: length-scale parameter\n', + '* Returns:\n', + '* numeric values of the GP function evaluated at l_gp\n', + '*/\n', + 'vector spd_cov_exp_quad(data vector l_gp, real alpha_gp, real rho_gp) {\n', + 'int NB = size(l_gp);\n', + 'vector[NB] out;\n', + 'real constant = square(alpha_gp) * (sqrt(2 * pi()) * rho_gp);\n', + 'real neg_half_lscale2 = -0.5 * square(rho_gp);\n', + 'for (m in 1:NB) {\n', + 'out[m] = constant * exp(neg_half_lscale2 * square(l_gp[m]));\n', + '}\n', + 'return out;\n', + '}\n') + } else { + model_file[grep('Stan model code', model_file)] <- + paste0('// Stan model code generated by package mvgam\n', + 'functions {\n', + '/* Spectral density function of a Gaussian process\n', + '* with squared exponential covariance kernel\n', + '* Args:\n', + '* l_gp: numeric eigenvalues of an SPD GP\n', + '* alpha_gp: marginal SD parameter\n', + '* rho_gp: length-scale parameter\n', + '* Returns:\n', + '* numeric values of the GP function evaluated at l_gp\n', + '*/\n', + 'vector spd_cov_exp_quad(data vector l_gp, real alpha_gp, real rho_gp) {\n', + 'int NB = size(l_gp);\n', + 'vector[NB] out;\n', + 'real constant = square(alpha_gp) * (sqrt(2 * pi()) * rho_gp);\n', + 'real neg_half_lscale2 = -0.5 * square(rho_gp);\n', + 'for (m in 1:NB) {\n', + 'out[m] = constant * exp(neg_half_lscale2 * square(l_gp[m]));\n', + '}\n', + 'return out;\n', + '}\n}\n') + } + model_file <- readLines(textConnection(model_file), n = -1) + + trend_model_file <- gsub('l_gp', 'trend_l_gp', trend_model_file) + trend_model_file <- gsub('k_gp', 'trend_k_gp', trend_model_file) + idx_data <- trend_mvgam$model_data[grep('l_gp', names(trend_mvgam$model_data))] + names(idx_data) <- gsub('l_gp', 'trend_l_gp', names(idx_data)) + model_data <- append(model_data, idx_data) + + l_lines <- grep('// approximate gp eigenvalues', trend_model_file, fixed = TRUE) + model_file[min(grep('data {', model_file, fixed = TRUE))] <- + paste0('data {\n', + paste(trend_model_file[l_lines], + collapse = '\n')) + model_file <- readLines(textConnection(model_file), n = -1) + } + + if(any(grepl('k_gp', trend_model_file))){ + idx_data <- trend_mvgam$model_data[grep('k_gp', names(trend_mvgam$model_data))] + names(idx_data) <- gsub('k_gp', 'trend_k_gp', names(idx_data)) + model_data <- append(model_data, idx_data) + + k_lines <- grep('// basis functions for approximate gp', trend_model_file, fixed = TRUE) + model_file[min(grep('data {', model_file, fixed = TRUE))] <- + paste0('data {\n', + paste(trend_model_file[k_lines], + collapse = '\n')) + model_file <- readLines(textConnection(model_file), n = -1) + + # Update the parameters block with gp params + start <- grep("// gp term sd parameters", trend_model_file, + fixed = TRUE) + end <- grep("// gp term latent variables", trend_model_file, + fixed = TRUE) + 1 + last <- end + for(i in end:(end+50)){ + if(grepl('vector[trend_k_gp', trend_model_file[i], + fixed = TRUE)){ + last <- i + } else { + break + } + } + gp_params <- paste(trend_model_file[start:last], + collapse = '\n') + + model_file[min(grep('parameters {', model_file, fixed = TRUE))] <- + paste0('parameters {\n', + gp_params) + model_file <- readLines(textConnection(model_file), n = -1) + } + if(any(grepl("int n_sp; // number of smoothing parameters", model_file, fixed = TRUE))){ model_file[grep("int n_sp; // number of smoothing parameters", @@ -2630,6 +2736,12 @@ add_trend_predictors = function(trend_formula, trend_model_file[grep('normal(0, lambda', trend_model_file, fixed = TRUE)-1]) } + + if(any(grepl('// prior for gp', trend_model_file))){ + spline_coef_headers <- c(spline_coef_headers, + trend_model_file[grep('// prior for gp', + trend_model_file, fixed = TRUE)]) + } spline_coef_headers <- gsub('...', '_trend...', spline_coef_headers, fixed = TRUE) @@ -2646,6 +2758,16 @@ add_trend_predictors = function(trend_formula, } } + if(any(grepl('// prior for gp', trend_model_file))){ + starts <- grep('// prior for gp', trend_model_file, fixed = TRUE) + 1 + ends <- grep('// prior for gp', trend_model_file, fixed = TRUE) + 4 + for(i in seq_along(starts)){ + spline_coef_lines <- c(spline_coef_lines, + paste(trend_model_file[starts[i]:ends[i]], + collapse = '\n')) + } + } + spline_coef_lines <- gsub('_raw', '_raw_trend', spline_coef_lines) spline_coef_lines <- gsub('lambda', 'lambda_trend', spline_coef_lines) spline_coef_lines <- gsub('zero', 'zero_trend', spline_coef_lines) @@ -2757,7 +2879,6 @@ add_trend_predictors = function(trend_formula, "// posterior predictions") } - model_file <- readLines(textConnection(model_file), n = -1) } diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R index 1070b735..4c843f01 100644 --- a/R/summary.mvgam.R +++ b/R/summary.mvgam.R @@ -180,84 +180,82 @@ if(object$family == 'Gamma'){ } } -smooth_labs <- do.call(rbind, lapply(seq_along(object$mgcv_model$smooth), function(x){ - data.frame(label = object$mgcv_model$smooth[[x]]$label, - term = paste(object$mgcv_model$smooth[[x]]$term, collapse = ','), - class = class(object$mgcv_model$smooth[[x]])[1]) -})) - -if(any(!is.na(object$sp_names)) & !all(smooth_labs$class == 'random.effect')){ - if(!is.null(object$trend_call)){ - cat("\nGAM smoothing parameter (rho) estimates:\n") - } else { - cat("\nGAM observation smoothing parameter (rho) estimates:\n") - } - rho_coefs <- mcmc_summary(object$model_output, 'rho', - digits = digits, - variational = object$algorithm %in% c('fullrank', - 'meanfield'))[,c(3:7)] - rownames(rho_coefs) <- paste0(object$sp_names, '_rho') +if(all(is.na(object$sp_names))){ + +} else { + if(any(unlist(purrr::map(object$mgcv_model$smooth, inherits, 'random.effect')))){ + re_labs <- unlist(lapply(purrr::map(object$mgcv_model$smooth, 'term'), + paste, collapse = ','))[ + unlist(purrr::map(object$mgcv_model$smooth, inherits, 'random.effect'))] + re_labs <- gsub('series', 'trend', re_labs) + re_sds <- mcmc_summary(object$model_output, 'sigma_raw', + ISB = TRUE, digits = digits, + variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] - # Don't print random effect lambdas as they follow the prior distribution - if(any(smooth_labs$class == 'random.effect')){ - re_smooths <- smooth_labs %>% - dplyr::filter(class == 'random.effect') %>% - dplyr::pull(label) + re_mus <- mcmc_summary(object$model_output, 'mu_raw', + ISB = TRUE, digits = digits, + variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] - print(rho_coefs[!rownames(rho_coefs) %in% re_smooths,]) + rownames(re_sds) <- paste0('sd(',re_labs,')') + rownames(re_mus) <- paste0('mean(',re_labs,')') - } else { - print(rho_coefs) + if(!is.null(object$trend_call)){ + cat("\nGAM observation model group-level estimates:\n") + } else { + cat("\nGAM group-level estimates:\n") + } + print(rbind(re_mus, re_sds)) } } -if(any(!is.na(object$sp_names))){ - cat("\nApproximate significance of GAM observation smooths:\n") - suppressWarnings(printCoefmat(summary(object$mgcv_model)$s.table[, c(1,3,4), drop = FALSE], - digits = min(3, digits + 1), - signif.stars = getOption("show.signif.stars"), - has.Pvalue = TRUE, na.print = "NA", - cs.ind = 1)) -} - -if(any(smooth_labs$class == 'random.effect')){ - re_smooths <- smooth_labs %>% - dplyr::filter(class == 'random.effect') %>% - dplyr::pull(term) - - if(object$fit_engine == 'jags'){ - re_sds <- mcmc_summary(object$model_output, - paste0('sigma_raw', - seq_along(re_smooths)), - digits = digits, - variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] +if(!is.null(attr(object$mgcv_model, 'gp_att_table'))){ + gp_names <- unlist(purrr::map(attr(object$mgcv_model, 'gp_att_table'), 'name')) + alpha_params <- gsub(':', 'by', gsub(')', '_', + gsub('(', '_', paste0('alpha_', gp_names), + fixed = TRUE), fixed = TRUE)) + alpha_summary <- mcmc_summary(object$model_output, alpha_params, + ISB = TRUE, digits = digits, + variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] + + rownames(alpha_summary) <- paste0('alpha_', gp_names) + + rho_params <- gsub(':', 'by', gsub(')', '_', + gsub('(', '_', paste0('rho_', gp_names), + fixed = TRUE), fixed = TRUE)) + rho_summary <- mcmc_summary(object$model_output, rho_params, + ISB = TRUE, digits = digits, + variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] + rownames(rho_summary) <- paste0('rho_', gp_names) - re_mus <- mcmc_summary(object$model_output, - paste0('mu_raw', - seq_along(re_smooths)), - digits = digits, - variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] + if(!is.null(object$trend_call)){ + cat("\nGAM observation model gp term marginal deviation (alpha) and length scale (rho) estimates:\n") } else { - re_sds <- mcmc_summary(object$model_output, 'sigma_raw', - ISB = TRUE, - digits = digits, - variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] - - re_mus <- mcmc_summary(object$model_output, 'mu_raw', - ISB = TRUE, - digits = digits, - variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] + cat("\nGAM gp term marginal deviation (alpha) and length scale (rho) estimates:\n") } + print(rbind(alpha_summary, rho_summary)) +} - rownames(re_sds) <- paste0('sd(',re_smooths,')') - rownames(re_mus) <- paste0('mean(',re_smooths,')') +if(any(!is.na(object$sp_names))){ + gam_sig_table <- summary(object$mgcv_model)$s.table[, c(1,3,4), drop = FALSE] + if(!is.null(attr(object$mgcv_model, 'gp_att_table'))){ + gp_names <- unlist(purrr::map(attr(object$mgcv_model, 'gp_att_table'), 'name')) + if(all(rownames(gam_sig_table) %in% gsub('gp(', 's(', gp_names, fixed = TRUE))){ - if(!is.null(object$trend_call)){ - cat("\nGAM observation model group-level estimates:\n") - print(rbind(re_mus, re_sds)) - } else { - cat("\nGAM group-level estimates:\n") - print(rbind(re_mus, re_sds)) + } else { + gam_sig_table <- gam_sig_table[!rownames(gam_sig_table) %in% + gsub('gp(', 's(', gp_names, fixed = TRUE),,drop = FALSE] + + if(!is.null(object$trend_call)){ + cat("\nApproximate significance of GAM observation smooths:\n") + } else { + cat("\nApproximate significance of GAM observation smooths:\n") + } + suppressWarnings(printCoefmat(gam_sig_table, + digits = min(3, digits + 1), + signif.stars = getOption("show.signif.stars"), + has.Pvalue = TRUE, na.print = "NA", + cs.ind = 1)) + } } } @@ -576,54 +574,68 @@ if(!is.null(object$trend_call)){ if(all(is.na(object$trend_sp_names))){ } else { - to_print <- vector(length = length(object$trend_mgcv_model$smooth)) - for(i in 1:length(object$trend_mgcv_model$smooth)){ - if(inherits(object$trend_mgcv_model$smooth[[i]], 'random.effect')){ - to_print[i] <- FALSE - } else { - to_print[i] <- TRUE - } - } + if(any(unlist(purrr::map(object$trend_mgcv_model$smooth, inherits, 'random.effect')))){ + re_labs <- unlist(lapply(purrr::map(object$trend_mgcv_model$smooth, 'term'), + paste, collapse = ','))[ + unlist(purrr::map(object$trend_mgcv_model$smooth, inherits, 'random.effect'))] + re_labs <- gsub('series', 'trend', re_labs) + re_sds <- mcmc_summary(object$model_output, 'sigma_raw_trend', + ISB = TRUE, digits = digits, + variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] + + re_mus <- mcmc_summary(object$model_output, 'mu_raw_trend', + ISB = TRUE, digits = digits, + variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] + + rownames(re_sds) <- paste0('sd(',re_labs,')_trend') + rownames(re_mus) <- paste0('mean(',re_labs,')_trend') + + cat("\nGAM process model group-level estimates:\n") + print(rbind(re_mus, re_sds)) + } + } + + if(!is.null(attr(object$trend_mgcv_model, 'gp_att_table'))){ + gp_names <- unlist(purrr::map(attr(object$trend_mgcv_model, 'gp_att_table'), 'name')) + alpha_params <- gsub(':', 'by', gsub(')', '_', + gsub('(', '_', paste0('alpha_', gp_names), + fixed = TRUE), fixed = TRUE)) + alpha_summary <- mcmc_summary(object$model_output, alpha_params, + ISB = TRUE, digits = digits, + variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] - if(all(to_print == FALSE)){ + rownames(alpha_summary) <- gsub('series', 'trend', paste0('alpha_', gp_names)) + + rho_params <- gsub(':', 'by', gsub(')', '_', + gsub('(', '_', paste0('rho_', gp_names), + fixed = TRUE), fixed = TRUE)) + rho_summary <- mcmc_summary(object$model_output, rho_params, + ISB = TRUE, digits = digits, + variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] + rownames(rho_summary) <- gsub('series', 'trend', paste0('rho_', gp_names)) + + cat("\nGAM process model gp term marginal deviation (alpha) and length scale (rho) estimates:\n") + print(rbind(alpha_summary, rho_summary)) + } + + if(any(!is.na(object$trend_sp_names))){ + gam_sig_table <- summary(object$trend_mgcv_model)$s.table[, c(1,3,4), drop = FALSE] + if(!is.null(attr(object$trend_mgcv_model, 'gp_att_table'))){ + gp_names <- unlist(purrr::map(attr(object$trend_mgcv_model, 'gp_att_table'), 'name')) + if(all(rownames(gam_sig_table) %in% gsub('gp(', 's(', gp_names, fixed = TRUE))){ } else { - cat("\nGAM process smoothing parameter (rho) estimates:\n") - rho_coefs <- mcmc_summary(object$model_output, 'rho_trend', - digits = digits, - variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] - rownames(rho_coefs) <- paste0(object$trend_sp_names, '_rho_trend') - print(rho_coefs[to_print,]) + gam_sig_table <- gam_sig_table[!rownames(gam_sig_table) %in% + gsub('gp(', 's(', gp_names, fixed = TRUE),,drop = FALSE] cat("\nApproximate significance of GAM process smooths:\n") - suppressWarnings(printCoefmat(summary(object$trend_mgcv_model)$s.table[,c(1,3,4), drop = FALSE], + suppressWarnings(printCoefmat(gam_sig_table, digits = min(3, digits + 1), signif.stars = getOption("show.signif.stars"), has.Pvalue = TRUE, na.print = "NA", cs.ind = 1)) - - } - - if(any(to_print == FALSE)){ - re_labs <- unlist(lapply(purrr::map(object$trend_mgcv_model$smooth, 'term'), - paste, collapse = ','))[ - unlist(purrr::map(object$trend_mgcv_model$smooth, inherits, 'random.effect'))] - re_labs <- gsub('series', 'trend', re_labs) - re_sds <- mcmc_summary(object$model_output, 'sigma_raw_trend', - ISB = TRUE, digits = digits, - variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] - - re_mus <- mcmc_summary(object$model_output, 'mu_raw_trend', - ISB = TRUE, digits = digits, - variational = object$algorithm %in% c('fullrank', 'meanfield'))[,c(3:7)] - - rownames(re_sds) <- paste0('sd(',re_labs,')_trend') - rownames(re_mus) <- paste0('mean(',re_labs,')_trend') - - cat("\nGAM process model group-level estimates:\n") - print(rbind(re_mus, re_sds)) - } + } } } @@ -649,8 +661,6 @@ if(object$fit_engine == 'jags'){ } else { cat('\nRhat looks reasonable for all parameters\n') } - - } } diff --git a/R/trends.R b/R/trends.R index ba5e7aff..32377390 100644 --- a/R/trends.R +++ b/R/trends.R @@ -71,66 +71,6 @@ sim_gp = function(last_trends, h, rho_gp, alpha_gp){ h)) } -#' Propagate a Hilbert basis GP -#' Evaluate Laplacian eigenfunction for a given GP basis function -#' @noRd -phi = function(boundary, m, centred_covariate) { - 1 / sqrt(boundary) * sin((m * pi)/(2 * boundary) * - (centred_covariate + boundary)) -} - -#' Evaluate eigenvalues for a given GP basis function -#' @noRd -lambda = function(boundary, m) { - ((m * pi)/(2 * boundary))^2 -} - -#' Spectral density squared exponential Gaussian Process kernel -#' @noRd -spd = function(alpha_gp, rho_gp, eigenvalues) { - (alpha_gp^2) * sqrt(2 * pi) * rho_gp * - exp(-0.5 * (rho_gp^2) * (eigenvalues^2)) -} - -#' @noRd -sim_hilbert_gp = function(alpha_gp, - rho_gp, - b_gp, - last_trends, - fc_times, - train_times, - mean_train_times){ - - num_gp_basis <- length(b_gp) - - # Get vector of eigenvalues of covariance matrix - eigenvalues <- vector() - for(m in 1:num_gp_basis){ - eigenvalues[m] <- lambda(boundary = (5.0/4) * - (max(train_times) - min(train_times)), - m = m) - } - - # Get vector of eigenfunctions - eigenfunctions <- matrix(NA, nrow = length(fc_times), - ncol = num_gp_basis) - for(m in 1:num_gp_basis){ - eigenfunctions[, m] <- phi(boundary = (5.0/4) * - (max(train_times) - min(train_times)), - m = m, - centred_covariate = fc_times - mean_train_times) - } - - # Compute diagonal of covariance matrix - diag_SPD <- sqrt(spd(alpha_gp = alpha_gp, - rho_gp = rho_gp, - sqrt(eigenvalues))) - - # Compute GP trend forecast - as.vector((diag_SPD * b_gp) %*% t(eigenfunctions)) -} - - #' AR3 simulation function #' @param last_trends Vector of trend estimates leading up to the current timepoint #' @param h \code{integer} specifying the forecast horizon diff --git a/R/update_priors.R b/R/update_priors.R index 5bb2745d..30845d86 100644 --- a/R/update_priors.R +++ b/R/update_priors.R @@ -39,8 +39,8 @@ update_priors = function(model_file, model_file, fixed = TRUE))){ # Updating parametric effects - if(any(grepl(paste0(priors$param_name[i], '...'), model_file, fixed = TRUE))){ - header_line <- grep(paste0(priors$param_name[i], '...'), model_file, fixed = TRUE) + if(any(grepl(paste0('// prior for ', priors$param_name[i], '...'), model_file, fixed = TRUE))){ + header_line <- grep(paste0('// prior for ', priors$param_name[i], '...'), model_file, fixed = TRUE) newprior <- paste(trimws(strsplit(priors$prior[i], "[~]")[[1]][2])) model_file[header_line + 1] <- paste(trimws(strsplit(model_file[header_line + 1], "[~]")[[1]][1]), '~', diff --git a/man/plot_predictions.mvgam.Rd b/man/plot_predictions.mvgam.Rd new file mode 100644 index 00000000..77f37263 --- /dev/null +++ b/man/plot_predictions.mvgam.Rd @@ -0,0 +1,107 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/marginaleffects.mvgam.R +\name{plot_predictions.mvgam} +\alias{plot_predictions.mvgam} +\title{Effect plot as implemented in \pkg{marginaleffects}} +\usage{ +plot_predictions.mvgam( + model, + condition = NULL, + by = NULL, + newdata = NULL, + type = NULL, + vcov = NULL, + conf_level = 0.95, + wts = NULL, + transform = NULL, + points = 0, + rug = FALSE, + gray = FALSE, + draw = TRUE, + ... +) +} +\arguments{ +\item{model}{Model object} + +\item{condition}{Conditional predictions +\itemize{ +\item Character vector (max length 3): Names of the predictors to display. +\item Named list (max length 3): List names correspond to predictors. List elements can be: +\itemize{ +\item Numeric vector +\item Function which returns a numeric vector or a set of unique categorical values +\item Shortcut strings for common reference values: "minmax", "quartile", "threenum" +} +\item 1: x-axis. 2: color/shape. 3: facets. +\item Numeric variables in positions 2 and 3 are summarized by Tukey's five numbers \code{?stats::fivenum} +}} + +\item{by}{Marginal predictions +\itemize{ +\item Character vector (max length 3): Names of the categorical predictors to marginalize across. +\item 1: x-axis. 2: color. 3: facets. +}} + +\item{newdata}{When \code{newdata} is \code{NULL}, the grid is determined by the \code{condition} argument. When \code{newdata} is not \code{NULL}, the argument behaves in the same way as in the \code{predictions()} function.} + +\item{type}{string indicates the type (scale) of the predictions used to +compute contrasts or slopes. This can differ based on the model +type, but will typically be a string such as: "response", "link", "probs", +or "zero". When an unsupported string is entered, the model-specific list of +acceptable values is returned in an error message. When \code{type} is \code{NULL}, the +default value is used. This default is the first model-related row in +the \code{marginaleffects:::type_dictionary} dataframe.} + +\item{vcov}{Type of uncertainty estimates to report (e.g., for robust standard errors). Acceptable values: +\itemize{ +\item FALSE: Do not compute standard errors. This can speed up computation considerably. +\item TRUE: Unit-level standard errors using the default \code{vcov(model)} variance-covariance matrix. +\item String which indicates the kind of uncertainty estimates to return. +\itemize{ +\item Heteroskedasticity-consistent: \code{"HC"}, \code{"HC0"}, \code{"HC1"}, \code{"HC2"}, \code{"HC3"}, \code{"HC4"}, \code{"HC4m"}, \code{"HC5"}. See \code{?sandwich::vcovHC} +\item Heteroskedasticity and autocorrelation consistent: \code{"HAC"} +\item Mixed-Models degrees of freedom: "satterthwaite", "kenward-roger" +\item Other: \code{"NeweyWest"}, \code{"KernHAC"}, \code{"OPG"}. See the \code{sandwich} package documentation. +} +\item One-sided formula which indicates the name of cluster variables (e.g., \code{~unit_id}). This formula is passed to the \code{cluster} argument of the \code{sandwich::vcovCL} function. +\item Square covariance matrix +\item Function which returns a covariance matrix (e.g., \code{stats::vcov(model)}) +}} + +\item{conf_level}{numeric value between 0 and 1. Confidence level to use to build a confidence interval.} + +\item{wts}{string or numeric: weights to use when computing average contrasts or slopes. These weights only affect the averaging in \verb{avg_*()} or with the \code{by} argument, and not the unit-level estimates themselves. Internally, estimates and weights are passed to the \code{weighted.mean()} function. +\itemize{ +\item string: column name of the weights variable in \code{newdata}. When supplying a column name to \code{wts}, it is recommended to supply the original data (including the weights variable) explicitly to \code{newdata}. +\item numeric: vector of length equal to the number of rows in the original data or in \code{newdata} (if supplied). +}} + +\item{transform}{A function applied to unit-level adjusted predictions and confidence intervals just before the function returns results. For bayesian models, this function is applied to individual draws from the posterior distribution, before computing summaries.} + +\item{points}{Number between 0 and 1 which controls the transparency of raw data points. 0 (default) does not display any points.} + +\item{rug}{TRUE displays tick marks on the axes to mark the distribution of raw data.} + +\item{gray}{FALSE grayscale or color plot} + +\item{draw}{\code{TRUE} returns a \code{ggplot2} plot. \code{FALSE} returns a \code{data.frame} of the underlying data.} + +\item{...}{Additional arguments are passed to the \code{predict()} method +supplied by the modeling package.These arguments are particularly useful +for mixed-effects or bayesian models (see the online vignettes on the +\code{marginaleffects} website). Available arguments can vary from model to +model, depending on the range of supported arguments by each modeling +package. See the "Model-Specific Arguments" section of the +\code{?marginaleffects} documentation for a non-exhaustive list of available +arguments.} +} +\value{ +A \code{\link[ggplot2:ggplot]{ggplot}} object +that can be further customized using the \pkg{ggplot2} package, +or a \code{data.frame} (if \code{draw=FALSE}) +} +\description{ +Convenient way to call marginal or conditional effect plotting functions +implemented in the \pkg{marginaleffects} package +} diff --git a/src/mvgam.dll b/src/mvgam.dll index b52e44d48f89ae694dc7b6e517c24eb7699d89fe..764c3c8312fe1cf513859d67453ab200d4b2b285 100644 GIT binary patch delta 77 zcmZpe;?OX~VL}J9`fK@#-M&maPdDCby~Wsii>dV%bL%aZ)?2Ktx7b>5vA5phXuZYR hdW);|7I*6{p4MBut+)8L-s0EV2-3QpeG~uBZ~$(QBUS(a delta 77 zcmZpe;?OX~VL}IU=t8lH-M&oB7LB)BZ!xysVrsp`+n*;mxA?U-g0yaD-^Bkj8~_$oA;SOw diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf index fdc1786ab531a2ae6ad99898057894e247f6d289..6eb07ea5bebee4d65437afe8732e7874c105673b 100644 GIT binary patch delta 8955 zcmZvh2|SeD`~MB$k*!Br%GzQnjD62iwrUbmn1r(LJ7d44Bq_?o*i}My#+H5077~Ur zW{}-rvW&6!zx8;2%lG&9@_Nm@?wRYHbD!(;d0*GL=Q|qm6AgJICnO~)DKDiYqaY`H z`J{xNyN%XkA6o$tEoEs*X_^22tkeuANG&UaRDoPlgT+F-4Dl0|YvL=kE8PeThC={O z4Uu7EWGv)eO$SENX zV1m?mCS;OqVp)*o)mtTRNpkrzLdJ9xEiH?S)l>}#es6Z3s7}oJ@<~do8qD8&uC5`i z`h7RB*7z#;Mx${5Nx_Un<5!1QzF((PNDcN^8}{Z{EkoBh@9XRmCq6FR&pX_g*^7v? zAYg*=RWVy5q~eSN_}&hYI05Wz6b}I1tZYys#_~;(`h%J~r`|xp7xGr6@~2{c@p0#% zTFz&gwk)9}_bRb4i2%n#A`1HZZzbWzmX>|~UY#F;xu&{<_YErDAosxm>w&$CHG;p{ z3wvGx>9U?=w1aKJ8Jp>Z!3*lR@;PaN(K7mf5EuEZ9T8V;uWwLR*9(5+hc#=Va9 z4ZuI-z}fXD;F?i*h>?gH+{=#*xlu4FJRyR;5WjrExq(luYg41pjC-B_MmGJ}@WfU? zm`&&-PTo)cO`V%fnu`mXfs=qIp?%aG&4Y4pA>rvKBp%pg#U$7^A{Y;&$i_WpJuMZ# z|C$PFj29c9W94*{o<{&2t9y}*2so=KWB=%$dS_ZNseA?d>Gd!NuvESRGtqPX;dk1) z1O;tX#+l`#4@lXtY>R1EQ8?7dv!JS4L*HD=*4%yi?t^!zv6o^Wxc9xBeQp1$k1x_B zik8*+4vn~kx7Sp1%NjqbUoxY`mu5So95q5JVR=5m9N!tqQGjdld;5;^l{LR>7H5nO zY`$=8J6~5{G}a8My{3*eL(qq0s-4|>efT2D@Bn&t;e7d(Y;xC~a5eJSchZ%b+s~jh zyHR7PEq%_67|Xxs$F8}TJxS|VNQ@X8VEz#F?d#}!DPCKG)8;cidH$@#P#-aQ$6I%} z)!=dl?pU7}who^GH}F-{aiOnrakOaFHNSz?j})`mv^m&zFR2eU?3k4ph*qa}*vz~Y zQpEz~C&oxC?DS?^8=FB3u|DOGck^-vcWtJZna8pV`Z*PM`6qnR&sW%U-j#A2&Te^P%;N1z z2-SQkOx%_HBhJV*F-KL3u?yf?uaHDi>JR%L) z{R}NX$lIeQM(@FC=B4`mimhl#Uf-?}m(_bCfV`-Ko8BJ>O^dAEl9y5V7$qlMZ*-C< zubp%Mx1@Dl)w!asP`M4;xss^(jJ-q{FYwFpP6bfG=WR9;Slc1hQFFMZJog~y?`eH+ zLWKM3*KZ+6x1>Eu#CM6aUp>bown#-CX|k{Ri>BSi#3~OLV)DMwB%0OaU&&=zrp3dN z0Qtqvr&?#u4%Di_NGo(o+&BL5UJ@KL3|aX$D8G36)o^27Xi8nV>5D>0Q{BG3SKr6t zSFtMdLW2z6H|NFbMGqW|W>I89`F(rO*Mr=qNKqyGBymMye&%?0`*NSRk1A@t6D&>q z`Y0uKYiN=3;YE-tDx ze+BcQew^>!y9PL7;s`qpR|RKNK{;o+ z0Q=Ge@qqa&J#BR=g@wxv&cdXJ)da`FVO&e{K=OzuTylXB@I@g9!9f}-L-4H@E~eQw zJWwHIzw`~0Ly*K`oA^r07>R0wDe9oEh48m318GhU@f zx)F`w%7@R^wf#^_d}3~BIm(9oc}^# z>ER$4*|EXXFVy36iBz!rpuB3485-!FQk+e&7b1!67RS$-eZ?(hB4WtxWND>_hV)A{ zEx|R2go9erOxp0P;c^Sxi>9us9;(hh9vNf{)QJ1uhpt^mll8tOSF`m#d_6Y>*Sqi< z=_k$ycxeXxm_s}C*Oy(QmkY9o3;7C?#0H(Vq6zT&*0j2^f|0j%Wksx`x|BYr?2*zO z^4y^ZN$lXsQ?H;pe9)juG{w0q9YKC>f~KTQ&mTUjzFJ@ALbFl1DeUHsSK%yB5p?ACxpe{L6UpwyjS|;6fVXYoahoOh~Urg1b3X zh*UsyErBgE2xS|+W$!3;{J_`LOEt|$^TtNmFmUSHNXkyf{1 zCp(Y17dks$J+?zy>-%|78Z`c(2XL!;_<|7N7LCWN6t^sd6W|;nHwhna6}b&S>NkEG zz74WBq&e_*zl5)Rm`H$oN5t2a>4+*J2%DEue-hw@&Y-k5%O0i>;0sCh z-qVH&LH4;EnBMxb*w`vze4Brb+J=>z9%=lncnyM2VNaJ-Fkt!8c*qRE_+Y+YO0FqG zU=Hwchy+K8MX^lZB1=l2bfZj=eJmlQAqzL>TFp1EvK@F z5JepBrw&$i?*~nyTD5T?FzvpVL8z8gU?M$SV6L8!&pba7i_VbBSWVLw@6lUKi`J!# zlcf180o^Y&Dr1@T>c7@}NaF$WKeu!;$;)ZBb#L)s>hD&VpzWV7rU%($e|KWW`JhE6 zls>kV6iSw=`Kvpiu}ZLpA}eb(3*wH*TatXhDrMQ-tuJyfvvY2t_`!DZ(qW62`>x8_ z{`7+n(7~BB1#dHDqu0ZfwZ02apBBdw0_4ekY)TvUTYX%lUw!DF(p1tfx`dRwDcw&W zj2PEh_t(F@LwTo?U*)D`d^uFIZ{UZ~zUb@h)gm6o#o4wrJL~Vs&k*U~vXoNTM5Rr# z0fo4Pbv;UXb7xM&Z^g46>}P2d}QX0dh2B` zcXCd3F^!Z`TM?)J|DoRV2z1R@2CzvV}|=v+L%+nWjwjJdpr_Tgbo3!o=RA=;-|Xwz+3ef1+uQ6@utre&lry?gM9t%>Oq?GS4-u@a%I9`BcfO?& zdy#%yIUART_izDHGt*KIydiZqYOt-oWKr@=np|Lg9HNL1KcK=ncet8XflZsMl|ikM zCC`TJ1ullh^^q?He~S#VkMaX&juxggShUxUVp&JZ#v<#G8qz`x@Q=LH-i^dy8>G9m z0N|Q57pa=+GEwhcKh7=lAo*LNqq{<;LRAM68-*h-3ip)HOS=fp_pym@rmE-@w5l85 zI9e>)aJMQV8{gX6RVy4A&SF3B<##)bA0!pf_h80V6ma^htlIp`{53H1J7oX)F66nI z2MqYP=V}m#RHa;2moFu9y;f(pf_U!<1LrvkOs*&UC=|}~%HZ{gcOKlF?c+NCaJ!O9 z?cw%;;vwcl7n>J9+72b7_6cbw@^E`@EnJ@gb^FUNTcP%x@;y)boZvMt4C)!3mP)Nh zqccTt__I1Hb(7QwnH~fi}IZ~`05L(h;1D{{1+Dl^iHeD{0XP_OtHj3n`vaREI&Y%RsCevg=06(XvU4% z^%0O7of(3<_~gNw2;-x@%fg<|w<0*GPw-X9cV)>7p?u}Y?mU;LXG*VyaENi)*2jn_ z(#LPDDxP7#&VARRg?Ic30w!YrZ1|-vyyQNlTW``B|M5vvCPt|GUjGcbD7%3UhW5@U zAlVm4g1n1>R0%{Fc~=e!T!3nrn6&2Z1VNj(e)K#xM_FZ}hUV>s-hWTG=p#h3a-T(i zGN@8@6%3>ch$tA zDjjwWEm>{~cTZk3qfAM-(4qw8@OQMkzL#O&(B?@mt-9rXbY0Kov@HW8TZyvgys2B0 zZh353%6<{``ysTkyvl!O3`v{soE}V@(R(!Sw4(&%N~^8Pk6MivVd$Dd*Lfp}Zo7Ic z=f={x##_PL=`g^fs4xmyfX#0H9;&@S-Ca&23*9!NXm-Y3$8X8PCi65ou%ScD-RX9? znLn}h5?e7pPXw;0q9*AYJ0+Vw>-!FQ@-d!VdF=#!63jF7Te`Ktp+AZ-Yol9*t^{ccXvT_wz$28=GEgSDeM- zVYNr&jk{21we8k(ACK?aD(#EQz*KC_HSD0vKT}tcwXbOn6HMpow7Pzj4sTqkh(1k_ zSPif2k8Tk7#vk=8Dv>7fB&+HQXMpblxLS;=-#O zHcMbDR*)hw$&c-nhL)<~qBx1`PTrSc-Ai8wEzm=|UzXc@IyAnx@C!TRc81=e;qF^m zlL9a)nO2PHJ5$;Gk%4x}mn-alaT{dwHu*MNmQ+JwCxIhvLcU(I5$D?apR1hGdO~&gvVr67M->bdVaJ}5FGY(y|R*rHx%P6H+DaI%j5L|LqaJ@23 zU+;1X+0L@y0mpL0qGC4LxJ+OX znKnyCVd`gGS}~SSUeHEtDgPLooUpSIfgQQ8id)iG3>y73^Yj0bRA$@iWat{3jdwh^ z`}q`DhlX1RowxOR5|5Ez+v*hf37}k#WdEZ^P}G z-le~xMN$&YK_{pjI49E=lRZ8%cb9AS@deLpoQEx~Gd078;>XX65{+E#-t(^%mb`Vr z?PQ~RK;NQVcd&FOkys2&zJ@)Y# z$}@ldzmzmy;#wGVl|%3^NA}kr{7R>w=YJ{pOiO7bt+pbgl==0)&{S2QOCagq9JrsC zQ$2_oO)XXorzJC|SB`(*Nl>U+-h5&AFnH`1H1$MvK|wMm(=Bt^5mefrm4Q?sr*hw@ zqCAxLWh)h6^dKR(m6eTq1^)8fWbPi*bmI!1d42v*F2vq8mGEMpD3NP)zg1U(B>d=t zNeRi_iocNiLJD{{X(W4ccCbC~)+t3lrHHL7UagZj9Z&e(9_vF9h?~zNo;|Ab1N@`;pim3+P~rP~r?m zO4S$Yk8w&s(W>Af$e?~fbA8(6w%bnADuf%b#6Gh_!+U z69GEgkjF?ve#+gs08no2gU}8eF^!)tRszNeIhsl17o-Z~Otj?hR;nF0&^-a4V ze8_x%S$ynYw%qPX$$B;4Wf+dgPpOwZ;V=|TyA`+_ARS8xci*Y&iCYRy=MmJqIPj)0 z*)XT7$UmLuPey?wI0(3Z4bmcTCfs`+bmm>|yz3=lLUF}k~*D~dA_H}Z=Z(~m@z4`hP5&@T)JMvbzvf{Y*FCx7 zI#b=(_?D)&z-ZXL=HpCwTIO=eG#AG8-LR{x$}v#KISO!Rmx7=W^k;mh5qr*EoG3xu zG?el=vbbh0EtW^o=4N}8>pR(kYp8Z5ry*Y9$^rSEp=c0kq-eZAkDKiwz$SRQa-p); zp7zwHNwCqQfWwc)3e)K^!L>v3*!R!pB@C(C?sxu78~I!tc*y3Od1*6z*~pz39eP7| zBrSBUX{!OX@V1Nj^D8LE*vE5wB;b|w%JMnfj$i8vS_ghq`Vydgs0szGx?`Ci1J#u9 z@K*yN2nt-Tp$VM2#{E8KP3@#exrvxPLXCAt1hlV92Mnsk01Y~S>sQTdbPiZ^Q!JoI zOuIr#h$aj-%*W-u)hw<LQ8{PpGBCGu259AGio|!vhE17pqabmU3JmJhzWe%LWV+ze{7#tXB zkDmWG7Thm!#3u}?+J;MVU^R;@rvEiN4PmV4hhMGeBKghh!ED!1sQhVuHF z#Ks==7^p6_O*>2FKB`|jER+V+?bn($s^iNhm}>Tm%+d`}kAc?N38tSJqJ?jwQ^{ z=>1DeG5dmJ|F#Gy@WH2aHKl$Avs}hXD7|=JN_L_!!w>HkJ11Te*$j zR^pW8LCO4t@Mtz;9A0~|J)$u25FtE+U_%=SHngfur8^9Y^300HoK0pwBV@1Lcj$q{ zctS33Gg8=Q+*6lp3k>xqxo&y5%%jnp9W;vaEFQOUerj@AE*KopY`mR4Mx*$Zjc4LY z7V(FL{_Q3yzAH;z4Mk4#(+6R11j$q}jyGwAhHM2-jbc|qaZf6N$Q;#XuLF^HukGfTQ#Y%NZMbmbLDM{?~+T3|6*xbfm+A9xO5Q>j%`X)O0R2>H4D1;`5p~eti_43Q(R*`?w^!EuKq90F* za+8_l(ied{O{Jb;)*cKqCjN{l!T2z1F~@{vBMwyFqTp1|lZP(2?V=;StAUwhP<$Tu zxI|uSY}i8O4u6&@y<)aAVMk%>Y0uGYrr{e5-4`-qYf>zM=eoF=$m|jFHFx@$W?pr0 z&}LnMUkUHhjK2jmnd2qSlxy{CrmCNrW=DTw^LQn=E$iCJp)pst%GP9a+?7_*bS#Ho zr2&WKT#diWt$ymB@T93qlhn8#GM3VFJz)q;`Ewxeip|Ou1he7@X$L*#@!z0&YRd)SGtLdr*4Ly~hlpKw?((n@U|z zosOtjR1irHvD<*~aMLS3H_h0$Q)pJRy{QUrL1|9_N^UZ0BqZ!Ip>oaqs2;C0N!1>7 zHA#u)ZZC*5%nvNf+%Tbw*$K_vKD>oo))b+kcN|7?YLx*!acVX{PsH#(%E2$5XdDf? zDMY+S9RrWzvX06$H}(`@Q1shlFZ?6NWz8ccXM|awAU~W-_UTR7CK~Ri&V#8i-w^C$ z?^Z;AI?{0i48dtD5xguJX05WUit`Ia80hEnA4>~_gnZS3lZG!j>528#;QPR@5z|Qb zfXqMIGPnzQ&Xqoe;>Zr4ExSBUHS;K@aW%F1LE08lVkDPeq zZJq$GSG((ej95%Hc(MeKT0-iw5g;BwI`9o2(V@=b^BJqem-!tuJo%u3F;n{V!q@$j z6{R2@T)^}KS>BdRQ-6`BV;lXVYGr?}l-;HbG&%+}I*3B2_SSRxgq}$w2+sy1F??Kd zyOUeuG|y{%5LrXB2&#$FG6?p%@G9~e+Mo2%$mi1{lCIo7HWb(V^VkIz2iGa_`-by40}$5e0vZX8{PQ_g3L5Am#Nxp z*OU%U@R)(>`EN=9N0_Wrv`A_$NnQBK7NU zZ+!sG{)@Sr9wb4cmR`PW%w+9ZfMxaEDy2ee`qcdN1!tqG>C)R^~!R%M1#pjE*C-tdg%t10K)ozGqAmDddCrzb(-5~06R6Pd$<05?A zTdst_Ih2pa%uCPv=^V=3RB}iaL#E*w^gh3iGS|9`0YGNEaCx|p-|$xP+R=F|Z#9_W z`u-fZj~PnGBg#A}p5A<=sPZAPn_#}-MFZvUU$2ZtZe7F$?XyBgl zeTvS(f~6vVdu4B6kE+85ulmGVBuVi`W`6TV?vEC2%QhRXA1ghat})9e#{hmhpR39SNBVb;qqf z7!Od#_|n8wF#k&Uir_60fOK?2V&d6>eK!Q*o$nv~woblj@3{@`Mo^wZ%6`}?W$oI( zVh^+SrLL26{hezIR{Q9UA4zxKG%gpO>4pN0f}o*)Fso_9$LR$pM?uS$A19>Q*pKc_ zn6*VXlKmjR`kY6KGw9$dxuf5MZjGwN~toF`I`=sldmUW&MfUsKcJ5ng`e2?lB9nau>P3-gjEa9Wprj)?)=VYNbX}WVrdcD{|`DLLQ z931QEvajcKamUeQzt{11_C&_6K|x42Nv_;?vHkD^nc;8+^_$vsN9S}|`>k8dY6|iK zl1G06ZUQI&ypd9pmY4fK*JP#sKi6axk$R?9cV%V%a+M@{9Y#qaNQe*-OazG%y)(+-NQ6{FkPwVcB3c+m9W{E3 zUdHH57^54V|H#e#|Nr^EdFFZM%syxDz23Ff`>wUmbP?sxBFacPAxTAPMQKSXNhzr- z3=+Do*0-K|*<290r6esWE%Sf(3Qci>l+w{kl)}pL&G+#&HPbL5MwhG_1m+FJ6}Hpg z1v$=GKIOgrgZYx?^vVP8MMLjSEb{xAYpeC!-rl<2#=wBq_GiDg=y)Mb%8BSHAr7Es z0=TkQODw21(V2^Pb*>&CaVsF^eH^N(Hn#}$PpKkZ*q}yKcM_GKSo^1RWitW6;P2KJ zJkJ-r2!6Cn4QsUs3h&_}zHe5T2<-L^*zOdlKA4vH_I`7yMx!Uk@NnOj%&P70XS(4% zp)_hyQxinCB@qe=faAID8cpDu+U>Gj)8pq}hdsUCFM_>xdwTOqO9w6n9@osZK z<%QfDM+$2D32H(6dNup=KF4dWz|Qg(F1coOYj33HaBT};3LKQ5Mgl|KLE9vsW0RmG zR?wzA-Gc%WdSGI&(cK@|Dr9}belW`;XqpjaX5Bf|ApDvCtGA>25k>3f^~yb*Dr;m1 z>rz_CIaaRE{JtW5fwx%?B)TfODs?{1tn2Pz_b2x?-sSf96fJGL)TtGYJNP-m%E$s2 z16E3-LU9LE^}6#1Kv#X>vD=h5@u|(!dBibYdGrbjqmocuI67^C!UQgS%hnxrv1B6m z+L+iP93Nti@m{U{j#bWr?KW4DTk@R0vae<35}hAwzgssj{UGQSV&RgP&y-r)>unck z^QjcM7q9!>WC$M8#E@$on3Dn>^dd z`}@5|3)6jT*QOsse>uy4-`G8c)`wPH^okha;*rft?=THd`q?m_ey%S9*(jvjYTqE7 zb81m6fu)d%85jV6-9?!3Bl3K>0j%@DcHJ^wiZ3oGKJIFtU$#=-xd(b=0!4U`VfM=4 zlq-AR{du995w8x1_Ds11r2JirEw4w(PKnqXIIk9mDzBZI@>E3%K z>e^rGr1ZmDVzraX=tGVEF=pZC(SlboE=%1=stIiY#33ztVD&oTaOQ#!n&05Q+;c<=k}&$-Y&I)ZvwCSQflv8GCQ zNYWQo62~W2T;J1gTfh3Q#&T{s_cA-Fk-yvl9Muj7>>J@zzOt|qW3`#gumYq^k{SGt z&I}%COnh?au3F5UtJW23o&AoUO;+6+nuL8c8K90el)Z#nE|@u zbG15z@n)5fYOhxn{NIl|O&727Ty=Q{Go~w+(ZglvNg5rm?>VpcTBc%7QEaT&v?{)< zmR4O(^hVS_?p4%T(N-J%x`iH3z$}X9yZfgPb|>b9!(Gg*7r818yZhs%eE4Hf{ z!@IV>FMG%vsw}~~jPgvpq!KE)m!lW1u?RAAwB751y>BUt_WU$Z6y=%!aj*}C{u(;s zLl7lh9CjlpXZQ+7^NHR7=!r|?uw3E1SW?np$D^K6dpFwvbl5`!88l|I%TmmAz3btp zH-0{kOh*Febr(?AhM;-gVx*^?#Ig?V3Lnfv2{`Udmi7?dbE zFpJ{sCCs3PFi(20=;i!bW;rf=D|D#-P{Z|GJH>t@7ER+x(gwz1Z;jUnl|;Wk3ou`Z zUwlC9BzN+>9(76{4<1DBIXI~r6jKJRcU8Q9Vy4Jn;WU4H9ELqqni)~voAgzHHNBRH z2sk^T=MyDRWQ5SlSdJfEvg=6G&u*+|u0%fP=atc3|Iv>dIxC_Zdmpr2UeZ+um?jTmJ=gYrMP$T(in^nFs`L0~@`B_SYs2QDJ#|2&k6~Zi z=L#B7_4g$RrH2ftuROMUtmbDcW#4DjovY|Gh%fb6UdqAe zeMmwJJ3f@HxLvLpy{CGlUQFsnfi4IZ|d0!Dyhypm-CN7clqO&y@7~KS2BhS0ArS)6F*~0eVC~gp2OVQaNgcZXuZj_!H{KB&$4A6=s8?$)crm zK{1E#4?!FIi^gKc%YD7N55;?sq=_@jwMep^LUl7rOfc}AlNqlWzb6_}y~|K)J6iZ5d@tap2%_F9NLaxlevy{8}l(dr^Lv8>Oj ztb}Pr`Lm9vyr7YF{$1A#z&jB_*CjRn@y}Xks^eu9Tw0;+Urk)v2EIDHpGB1xnCCyU zU-xee#p zT*QR(N=>7+R|w;Nyr+@w%O*h}EgUsv0(wxu3E%&!gqewO%{1@Gr{~@mREG!3fki z>+CX$G2L;y%fBzulhF0JX{f<^lSI_;tgA|T$?8ZPho#{^W6{ZhJO0l6o$e|Un$YTK zLunykx1yC0#JPg7MWm30K_~5Ym!n716QS<{s_oQEJ>?7S0esTJ_xWfN(SNj|m5hfD zwPU5yGh4^v$wqp;o*3EBs4{F)P|;n4c!QD0vS`muaE0-Pd8d2nHC{#GPZL`w?265! z)V?Y6cmogpMJc7yJ9V`SBX z7jJrYj4yIXmI7>+b^d}GcYON%$EEu4#mgBV4VT}asSkB!eBBT_5i#Sc8!38m`=TRU zaYkLKJQTdaFR>Hj*7}ms0vviUdRP|nU|HyA?5!{K%$&&rka=lFi|Dy73;QE)@i1>_LpVz^H+QYeRHSr zeeTPO{MDD97t%zF=XrEjv`0q{vUye9a=!yFzZ*LecfCvp6&gF_H5D(JIpaFe{CM#+ zrxc6F65Y)_H;s*n>)Zg%6B6(`0#?+sK9RgCr8Z@CfD({>Ycw^NLi+o*h~`eUknu*tK;)6 zed(cTlII&kYs#K7_Ja}P(sOL#=g+F>_&>qj#pJ9)F zCMOQ0#eTrox8nAu1&7W36?2IOPMh2X5E;r)@}>t_HDI_UQZTfFWPj~UU)qiHHjVsm zpE8wFwk2UOC|Cb-r$OBQiY^&<#GmEPLSvX0OkzMfhU3NFTqWlKCxOe3TgVUoG ztiaQ2Wun;&@T=NkLiHCJ`;GhYctZ``>;s9mXV{V0Q)RJDiM4K7Y>9`Q%*$_4#U;jN zfWDm`M{d5~Iy40@PCU2R)Mrj*7}gHY@Ol0hrH!{KuIr0=Os79uVtzm_=;O|+N~|cb zS#x>Cuxh8o;!H)bG>L+b>8Stkuc8e7Qk71&vc*UZi~wlyml65h{&z(%*OhnJbdj~* z=_+6P)g+BDO3lD;>BtQeeLGh}Px{%J6nX1AsV<0kL0#d$t#iKd{+^0iGn0r8{MA`~ zTlNbve=Yd$j}(Zn1V$9RuzK1Ba4$c;%I2KEzq*!h)e?0RPnOjQzRV1JhpuwbC(Gi5 zlKty#$z6=KF>6=u3sWd>Ilz&_#LTF=Tt5+ zz@A)6TfVl0E!9=sN^=@a2>fs>aP*!49t<-^e9NQe0+$6G*8VfV9a!QJ9}20NQhWUA z+#==+>;=zpWTU%9KXY5svwwZfNo(YtV(+@lDRExurN+jW7h!+<=x%q1?rTJOz&7dx z2LH182Q@qKht2bnHx)NV`kWF17Bxn4OSd#GE)R^G%var_Ul!P^(k4y6+>TUKRFfhtlzoSe`2YdVZJfmSnkT`9zvXQTCafZ%!1z5 zoex5K+d&LA*K_rP&`*+!GAa-?A~!~PiUwYA4*qsllMnVTIaOi*(p2~;hDvlCarOQ| z2(7XMXJ+a{bNB5>%3VKG3`kdZ8G7O`>a?BYW|7zrP}_<`f6T7ct;e+B;Hq+&C$Z|; z4@b)6HOb8jHW-9w?H7!6dXHTT0h4o;@yIasA3lDfm;?9jJ@?Pv|B?c^$;M@@6Pj+l z_`v@#-+en&A%N-(BLm4hqol5;vqk@Rj9nF;ESfy#pwqgqAc6qmoul=zA_-g1EHvNX z>kg*{TjsASjNg=xMLr)xGS4^7$|z6~^ms z@mbDbZq!}8eEl2kl@yi_R;cT#Dk0zLYU=M83|{GX)-!!m&90G5DKI ztRRWsQ~@xR`bH3`Zv+Je$?`Ubr|zLgTd{fT0MSFx4(i^#J^JzNf)Jj-lU1!!gdK6OM%PFHlYl7PNHUDs6O$#W(o zn2RTDk*kM?JYA|vSCe$I?p@-rt}~5|W`-&mYvcFXwI)(>Rc6ELw8pvqk(sFMIBX7X*sNdYmF&dHID zm1m++PzX?^#aVQHR1=@FWRLlm^$U$^LM0tFw^|JtQ|o@(ln*=C>0`k~_`pTP?t8!t zPL!FRTB;3>b-EMR*2=!{CmYDsqQeU;Eu9GQ1^=r*|9)n`GqI-msNT!P!M31D2^elp z7XkJ4L@%}vZ(>YpTPuJYFpFU@u!R}s=Ivv-O5nN zk6j+Ax4)?;t3wgS0DByzfp6>{PXkK=rww|1slB+W9ImE^WUDgGi^X!g8K+#KMyMPr zf7WuXD-9U8x#JDqjLQ08*J9N9wa?U8Sr3V61!W!R^ha1kD$Iy^tT5_)Yf@oo>-R*x z*$$e{6W1T)_ka^5W;eu%hV(fE!JYco?W0iHHPvlCqX=i{ABD}G7tU|Aj1!Twd9jL6=N9lFfro<&3x3B;Yw&nsL&nA}|RP!s|>NU|k^F0GEgQ zFV|awL6N=1!sGXi|Ek06lSx$rJJ-wI-bdlQ!#1dIL4E9koJ+hmulSn3gG6v~#QjrM z$a9z{_%8hkPx?p=l4jC3$j4FhHJ7^eN-dRu0lAfc+1iUbSMNEcXFO%D5be63z6>c0 z{|)*C^U|qumTrqQz`C9T58r;gmZiRzMX|vdh&{jUHH*u{w4Rbv2Hvn;$iK<4Ql*I> z*od{`5lrQ@PWCso*lbX&SYFIKdyu5P*)n^|P?OW~S_R~Qs0<>9hm!&-r{FNqPOzZq za7MHn14rOlJRz~lXiuOB}An`ffP>uV#>u|PmI!kninSGca>+;MYG${ zd=lk<36?E|id@7J>z$f$@al!R6NhE4yOjSjBO%1IxfM6t5cXv+d^t0BJHH3*ux*dd zZgeFSWxSnW(xOicX}|i!Mgy-o(l?0N_$Nnn(3Nw@)qR!N3vZ1j1)L5tg+H(W-0}=0 zIq0-|HQujM^h9;gG{^=UXxfhNE4mmJ=Ra{K1r~yH`}b60nW%+EvBgPqt;59|T@20n zIIR3Ge!1!j(Qy0hw21xNLB^49YkK`=DynZ!d@Z#9Ga+yRLS4(FIJR#tEyuDzc0(6T z{J^!hTg8pL9|3;_BsdU>?bm()O2$m%97E<;>ES0=gb6a)YJ!6-vNq zR+vS0;zP~mQvz#ODVjc3kbPGtMt1mAtVR%C=nLi(rpk_Bt@_X7oX zXzmiZStq#|DZ2V3gk`TcA#dVEeidma)-6fGO!tmS@a-m10L{Li zYvo%r<%MFsnJ0&jPV?eU@kZ0;S?t=AZ!WV;4$lcoDY#@MGIAl?&#zG+Qz<5%rzFx0 zt57%lQ#2QYU2YQC<+8dt1yLb&sgVP;tDk%?>guzPdIVcNzRqw|kEPk=5#ASf>^}<$ zW8=k@=ZYSJ%Y;+&_F+(9Isyq>@bB!^*<=c2Ax(3s_CiJR(Al$Y?_a(3FDl3x-St(w zV0}j=*y@qNb9!&s8!9q21jFif_fI~qQ1cRJF1s-ano2b=jDPE-eoe!hLkBT3S|!>j z#0U`$sgQ-%X$9i%eu+{x^Y+i2$Rt?pw<(b%z1l+OS}8kVsL*nYV~vuxlEuv_x1kw9 z3zjcRUiH>m(F&*1BwDwQp}0TY*gi$`{LX3K1o+vk^QH}z-DP>&wT1KfH|s1a=37^8 zbD0fh)4Anxxh;unk*O}U5UZ&qXaDFj{QF12j`NVFsBA4WZPGIV;YW01B~*#SYXd;? zq!lHYe{>)SUq1!wGqktx_~MGiFB;zyBa^QZ{E<^ zRQ|E6E@ckt8D+DHf$SBg5vrENAv3DDrGF7l#+gOhPIBq2)?*DEX6Rll((pRcLBcBS z*VlqzEm2^yG_C4A3>8&9hz_=5GD-haH~>0`5PXjsM8A_;^rH0J(rRC6B^6D-?F}RDcCA8_K!bHw{7f2SHV=9ZSvj=mEc=I0jV$GybyTf7DhP>Y)<#* zi`$j9D3Hw5Nc%87Sas1-t(5+R{#oaJh93-@U`q}r!5)|d>h!8wnUi@A3iKq(6P#J_ zfnE5GrJ@Wr(6eRdSzuK(@xfMAgYW<7y#Pop8jxC2_e~OlLefTTZb^Bj!he)|zJ$h|H@q7Gf+i^E+B^cg68$G1%3@1_~tU`?eEq7k@(~ zGw3Jy7+H97SblPNjI2lfEvK7_cV;T*1yIwZ6wuCN_8UkPLvU_8$KPFa5V*i3*XX~8 zp^R&li;rP9+ded?%`%?_urU!T0CvmU_h+Jt^O(4Tg_Z_$N(NY70wP042zoqQ1H}APcwvt%A}kh^e#Bz(27XRn7#3 zh`b94TdK8%NA!?EW{A->{y~#1^Ha$TlW`U%oc}g*C*d+H=)TDz;CT0UKotgsh@Jq%K#Y}S{iw09e^k1X5y`Pl+c__qh-gM=L{ zQSw@~HrW}@!nu!r+3`+tNnBlXyh2?4pb#dSO~Y$NFCW-4jdrGo%)1G>{yB0{H;3g{ zo)FV5S9>9KD!lSJ4eyK%q*}9*MS%0eZ&l(s#ruHbb-Pgs4?P6K;>PrcdnX_P5aljg zx8ynOoaGGViAWd?PR-PUOb|RFCH`e3pQN$*_G+yBzs#e9U_cNUGs>$X1v8&cO?9`` zO7Z`=wxbhBvrBDZ?*;n;Iz48xWG+(dwgZdX61d9xH~n zc?JQ>9Lpw#*uvElBi@1U9>B@ zjUEE$n~a-Z>iXR!f=b_5dpS=4HWmijjYj<+Z(+Y;z%{e|zr8Kkm7{3bqlOHEb)Y-6 z#R7ij4(tnPPT|wq!U1*o7(DJgnC$zOe_DsYV5`HkZn->exw8ATSD_&fpM-~0PS0JH z+g1*t1s6m2OLb*?olhmhzXO)|sgNOP6IH`AarmCv6@qG>GH9vb*i zAxkK3InS2Y(Y;o(z5|*bd#;l9OD+(o&_XmP6x3rRGbKuw0mUGWxat1YB?U6kNR zFs2T@DgPBr5XoMH2q)^(bwrn^%~j4ZFU=zZ7XyLcDBlY*xK}J)dQ>9~Ix9Ar?aGeM zP}K%03}+L<0Hy&(oIiHJv>4Q#Yr&Erqfg$A{3#Wc7_Qmrnno0Fn{19-+hc^`5jzQU6x3X>){3RnApp(zG7w_M&y6Ik5} z04M04C>~hXm_YUI0P$CS0Xsm}UqJ!_f7pxxWk9&rw2z$pD2tl65+Z?GpxJF7*G^qp z-!^3fOAvnYK}=(yKGjKkqWfJ$^uiEYCX+Kbj1vk_9oI z?8&SKs4W#iV4%u2R{Q>$2{uK@h}uGX5suL3ER|L(JBnD4(K_l3M^yYxJAKP5SBRkM z8!9z3)m*^32m<`$Acf~vbNe+4SppLc2VY z8Odl16H5bWx&L@v8vIE9*O-jde{Ll$EhGP*%SubjO36o`g=x~u%P|NFYUpV){15#6 BO*;Sp