diff --git a/NAMESPACE b/NAMESPACE index 0cae7e37..28852917 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -35,7 +35,6 @@ S3method(pairs,mvgam) S3method(plot,mvgam) S3method(plot,mvgam_forecast) S3method(plot,mvgam_lfo) -S3method(plot_effects,mvgam) S3method(posterior_epred,mvgam) S3method(posterior_linpred,mvgam) S3method(posterior_predict,mvgam) @@ -70,7 +69,6 @@ export(pfilter_mvgam_fc) export(pfilter_mvgam_init) export(pfilter_mvgam_online) export(pfilter_mvgam_smooth) -export(plot_effects) export(plot_mvgam_factors) export(plot_mvgam_fc) export(plot_mvgam_pterms) @@ -133,7 +131,6 @@ 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) diff --git a/R/dynamic.R b/R/dynamic.R index 6bf8a347..c854af8e 100644 --- a/R/dynamic.R +++ b/R/dynamic.R @@ -233,9 +233,9 @@ interpret_mvgam = function(formula, N){ k <- term$k if(is.null(k)){ if(N > 8){ - k <- min(40, min(N, max(8, N))) + k <- min(40, min(N - 1, max(8, N - 1))) } else { - k <- N + k <- N - 1 } } diff --git a/R/gp.R b/R/gp.R index 7c5b45d3..974aa9b3 100644 --- a/R/gp.R +++ b/R/gp.R @@ -140,11 +140,7 @@ make_gp_additions = function(gp_details, data, 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_names <- clean_gpnames(gp_att_table[[covariate]]$name) gp_stan_lines <- paste0(gp_stan_lines, paste0('array[', gp_att_table[[covariate]]$k, '] int b_idx_', @@ -322,8 +318,7 @@ scale_cov <- function(data, covariate, by, level, scale = TRUE, } if(is.na(max_dist)){ - Xgp_max_dist <- (abs(max(Xgp, na.rm = TRUE) - - min(Xgp, na.rm = TRUE))) + Xgp_max_dist <- sqrt(max(brms:::diff_quad(Xgp))) } else { Xgp_max_dist <- max_dist } @@ -453,10 +448,7 @@ prep_gp_covariate = function(data, covariate_mean <- mean(Xgp, na.rm = TRUE) covariate_max_dist <- ifelse(scale, - abs(max(Xgp, - na.rm = TRUE) - - min(Xgp, - na.rm = TRUE)), + sqrt(max(brms:::diff_quad(Xgp))), 1) # Construct vector of eigenvalues for GP covariance matrix; the @@ -484,7 +476,8 @@ prep_gp_covariate = function(data, scale = scale, initial_setup = TRUE) - # Make attributes table + # Make attributes table using a cleaned version of the covariate + # name to ensure there are no illegal characters in the Stan code byname <- ifelse(is.na(by), '', paste0(':', by)) covariate_name <- paste0('gp(', covariate, ')', byname) if(!is.na(level)){ @@ -507,9 +500,7 @@ prep_gp_covariate = function(data, # 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) + covariate_name <- clean_gpnames(covariate_name) 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, '') @@ -531,6 +522,26 @@ prep_gp_covariate = function(data, eigenfunctions = eigenfunctions) } +#' Clean GP names so no illegal characters are used in Stan code +#' @noRd +clean_gpnames = function(gp_names){ + 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) + gp_names_clean <- gsub('.', '_', gp_names_clean, fixed = TRUE) + gp_names_clean <- gsub(']', '_', gp_names_clean, fixed = TRUE) + gp_names_clean <- gsub('[', '_', gp_names_clean, fixed = TRUE) + gp_names_clean <- gsub(';', '_', gp_names_clean, fixed = TRUE) + gp_names_clean <- gsub(':', '_', gp_names_clean, fixed = TRUE) + gp_names_clean <- gsub("'", "", gp_names_clean, fixed = TRUE) + gp_names_clean <- gsub("\"", "", gp_names_clean, fixed = TRUE) + gp_names_clean <- gsub("%", "percent", gp_names_clean, fixed = TRUE) + gp_names_clean <- gsub("[.]+", "_", gp_names_clean, fixed = TRUE) + gp_names_clean +} + +#' Update a Stan file with GP information +#' @noRd add_gp_model_file = function(model_file, model_data, mgcv_model, gp_additions){ rho_priors <- unlist(purrr::map(gp_additions$gp_att_table, 'def_rho')) @@ -547,9 +558,7 @@ add_gp_model_file = function(model_file, model_data, mgcv_model, gp_additions){ # 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) + gp_names_clean <- clean_gpnames(gp_names) s_to_remove <- list() for(i in seq_along(gp_names)){ s_name <- gsub('gp(', 's(', gp_names[i], fixed = TRUE) diff --git a/R/plot.mvgam.R b/R/plot.mvgam.R index 1014e0ab..c57e9e5f 100644 --- a/R/plot.mvgam.R +++ b/R/plot.mvgam.R @@ -141,7 +141,7 @@ plot.mvgam = function(x, type = 'residuals', dplyr::filter(!label %in% gsub('gp(', 's(', gp_names, fixed = TRUE)) -> smooth_labs } n_smooths <- NROW(smooth_labs) - if(n_smooths == 0) stop("No smooth terms to plot. Use plot_effects() to visualise other effects", + if(n_smooths == 0) stop("No smooth terms to plot. Use plot_predictions() to visualise other effects", call. = FALSE) smooth_labs$smooth_index <- 1:NROW(smooth_labs) diff --git a/R/plot_effects.R b/R/plot_effects.R deleted file mode 100644 index 07efa3c6..00000000 --- a/R/plot_effects.R +++ /dev/null @@ -1,53 +0,0 @@ -#' 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 -#' @name plot_effects.mvgam -#' @inheritParams marginaleffects::plot_predictions -#' @return A \code{\link[ggplot2:ggplot]{ggplot}} object -#' that can be further customized using the \pkg{ggplot2} package -#'@export -plot_effects <- function(object, ...){ - UseMethod("plot_effects", object) -} - -#' @name plot_effects.mvgam -#' @method plot_effects mvgam -#' @export -plot_effects.mvgam = function(object, - condition = NULL, - by = NULL, - newdata = NULL, - type = NULL, - conf_level = 0.95, - wts = NULL, - transform = NULL, - points = 0, - rug = FALSE, - ...){ - # Set colour scheme - col_scheme <- attr(bayesplot::color_scheme_get(), - 'scheme_name') - bayesplot::color_scheme_set('viridis') - - # Generate plot and reset colour scheme - out_plot <- plot_predictions(model = object, - condition = condition, - by = by, - newdata = newdata, - type = type, - vcov = NULL, - conf_level = conf_level, - wts = wts, - transform = transform, - points = points, - rug = rug, - gray = FALSE, - draw = TRUE, - ...) + bayesplot::bayesplot_theme_get() - color_scheme_set(col_scheme) - - # Return the plot - return(out_plot) -} diff --git a/R/plot_mvgam_smooth.R b/R/plot_mvgam_smooth.R index a7b3d85f..d983f631 100644 --- a/R/plot_mvgam_smooth.R +++ b/R/plot_mvgam_smooth.R @@ -121,7 +121,7 @@ plot_mvgam_smooth = function(object, if(any(grepl(object2$mgcv_model$smooth[[smooth_int]]$label, gsub('gp(', 's(', gp_names, fixed = TRUE), fixed = TRUE))){ - stop(smooth, ' is a gp() term. Use plot_effects() instead to visualise', + stop(smooth, ' is a gp() term. Use plot_predictions() instead to visualise', call. = FALSE) } } diff --git a/R/predict.mvgam.R b/R/predict.mvgam.R index c70824a7..5154a1ea 100644 --- a/R/predict.mvgam.R +++ b/R/predict.mvgam.R @@ -124,6 +124,9 @@ predict.mvgam = function(object, newdata, attr(all_linpreds, 'model.offset') <- 0 # Trend stationary predictions + if(!process_error){ + family_extracts <- list(sigma_obs = .Machine$double.eps) + } trend_predictions <- mvgam_predict(family = 'gaussian', Xp = all_linpreds, type = 'response', diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R index f5aef222..a8345a01 100644 --- a/R/summary.mvgam.R +++ b/R/summary.mvgam.R @@ -238,7 +238,8 @@ if(!is.null(attr(object$mgcv_model, 'gp_att_table'))){ 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')) + gp_names <- clean_gpnames(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))){ } else { @@ -596,7 +597,8 @@ if(!is.null(object$trend_call)){ } 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')) + gp_names <- clean_gpnames(unlist(purrr::map(attr(object$trend_mgcv_model, + 'gp_att_table'), 'name'))) alpha_params <- gsub('gp_', 'gp_trend_', gsub(':', 'by', gsub(')', '_', gsub('(', '_', paste0('alpha_', gp_names), fixed = TRUE), fixed = TRUE))) diff --git a/man/plot_effects.mvgam.Rd b/man/plot_effects.mvgam.Rd deleted file mode 100644 index b8979e2e..00000000 --- a/man/plot_effects.mvgam.Rd +++ /dev/null @@ -1,84 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plot_effects.R -\name{plot_effects.mvgam} -\alias{plot_effects.mvgam} -\alias{plot_effects} -\title{Effect plot as implemented in \pkg{marginaleffects}} -\usage{ -plot_effects(object, ...) - -\method{plot_effects}{mvgam}( - object, - condition = NULL, - by = NULL, - newdata = NULL, - type = NULL, - conf_level = 0.95, - wts = NULL, - transform = NULL, - points = 0, - rug = FALSE, - ... -) -} -\arguments{ -\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.} - -\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{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.} -} -\value{ -A \code{\link[ggplot2:ggplot]{ggplot}} object -that can be further customized using the \pkg{ggplot2} package -} -\description{ -Convenient way to call marginal or conditional effect plotting functions -implemented in the \pkg{marginaleffects} package -} diff --git a/src/RcppExports.o b/src/RcppExports.o index c1d41286..17845786 100644 Binary files a/src/RcppExports.o and b/src/RcppExports.o differ diff --git a/src/mvgam.dll b/src/mvgam.dll index 22a11451..30d04477 100644 Binary files a/src/mvgam.dll and b/src/mvgam.dll differ diff --git a/src/trend_funs.o b/src/trend_funs.o index ad153f68..8f575e59 100644 Binary files a/src/trend_funs.o and b/src/trend_funs.o differ diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf index 7e464b9d..cf61d416 100644 Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ diff --git a/tests/testthat/test-dynamic.R b/tests/testthat/test-dynamic.R index fa07394c..2a5aa437 100644 --- a/tests/testthat/test-dynamic.R +++ b/tests/testthat/test-dynamic.R @@ -32,10 +32,10 @@ test_that("dynamic to gp Hilbert is working properly", { 'gp(time, by = covariate, c = 5/4, k = 17, scale = TRUE)', fixed = TRUE) - # k will be fixed at N if N <= 8 + # k will be fixed at N-1 if N <= 8 expect_match(attr(terms(mvgam:::interpret_mvgam(formula = y ~ dynamic(covariate), N = 7)), 'term.labels'), - 'gp(time, by = covariate, c = 5/4, k = 7, scale = TRUE)', + 'gp(time, by = covariate, c = 5/4, k = 6, scale = TRUE)', fixed = TRUE) })