diff --git a/NAMESPACE b/NAMESPACE index 50b88ca2..0cae7e37 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -35,6 +35,7 @@ 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) @@ -69,6 +70,7 @@ 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) @@ -78,7 +80,6 @@ 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) @@ -104,9 +105,6 @@ importFrom(brms,pstudent_t) importFrom(brms,qstudent_t) importFrom(brms,rstudent_t) importFrom(brms,student) -importFrom(ggplot2,theme_bw) -importFrom(ggplot2,theme_get) -importFrom(ggplot2,theme_set) importFrom(grDevices,devAskNewPage) importFrom(grDevices,hcl.colors) importFrom(grDevices,rgb) diff --git a/R/get_linear_predictors.R b/R/get_linear_predictors.R index 22d4bbf2..f95db526 100644 --- a/R/get_linear_predictors.R +++ b/R/get_linear_predictors.R @@ -79,7 +79,7 @@ trend_Xp_matrix = function(newdata, trend_map, series = 'all', trend_test <- newdata trend_indicators <- vector(length = length(trend_test$time)) - for(i in 1:length(trend_test$time)){ + for(i in 1:length(trend_test[[1]])){ trend_indicators[i] <- trend_map$trend[which(trend_map$series == trend_test$series[i])] } @@ -132,7 +132,7 @@ trend_Xp_matrix = function(newdata, trend_map, series = 'all', # Compute eigenfunctions test_eigenfunctions <- lapply(seq_along(gp_covariates), function(x){ - prep_eigenfunctions(data = newdata, + prep_eigenfunctions(data = trend_test, covariate = gp_covariates[x], by = by[x], level = level[x], diff --git a/R/gp.R b/R/gp.R index 99cc2ac3..db24d7d7 100644 --- a/R/gp.R +++ b/R/gp.R @@ -119,7 +119,8 @@ make_gp_additions = function(gp_details, data, L = L[x], mean = mean[x], scale = scale[x], - max_dist = max_dist[x]) + max_dist = max_dist[x], + initial_setup = TRUE) }) eigenfuncs <- rbind(eigenfuncs, @@ -167,40 +168,22 @@ make_gp_additions = function(gp_details, data, #' 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) + termlabs <- attr(terms(formula, keep.order = TRUE), 'term.labels') + return(grep('gp(', termlabs, fixed = TRUE)) } #' 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) + termlabs <- attr(terms(formula, keep.order = TRUE), 'term.labels') - termlabs <- attr(terms(formula), 'term.labels') - - # Drop these terms from the formula + # Replace the gp() terms with s() for constructing the initial model 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])){ @@ -221,17 +204,6 @@ gp_to_s <- function(formula){ } newformula <- reformulate(termlabs, rlang::f_lhs(formula)) - - # 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) } @@ -379,7 +351,8 @@ prep_eigenfunctions = function(data, mean = NA, max_dist = NA, scale = TRUE, - L){ + L, + initial_setup = FALSE){ # 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 @@ -410,11 +383,23 @@ prep_eigenfunctions = function(data, 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 + # for the observations where the by is a different level; + # the design matrix is always sorted by time and then by series + # in mvgam + if(initial_setup){ + sorted_by <- data.frame(time = data$time, + series = data$series, + byvar = data[[by]]) %>% + dplyr::arrange(time, series) %>% + dplyr::pull(byvar) + } else { + sorted_by <- data[[by]] + } + full_eigens <- matrix(0, nrow = length(data[[by]]), ncol = NCOL(eigenfunctions)) full_eigens[(1:length(data[[by]]))[ - data[[by]] == level],] <- eigenfunctions + sorted_by == level],] <- eigenfunctions eigenfunctions <- full_eigens } else { eigenfunctions <- eigenfunctions * data[[by]] @@ -435,7 +420,7 @@ prep_gp_covariate = function(data, k = 20){ # Get default gp param priors from a call to brms::get_prior() - def_gp_prior <- brms::get_prior(formula(paste0(response, ' ~ gp(', covariate, + def_gp_prior <- suppressWarnings(brms::get_prior(formula(paste0(response, ' ~ gp(', covariate, ifelse(is.na(by), ', ', paste0(', by = ', by, ', ')), 'k = ', k, @@ -443,7 +428,7 @@ prep_gp_covariate = function(data, scale, ', c = ', boundary, - ')')), data) + ')')), data)) def_gp_prior <- def_gp_prior[def_gp_prior$prior != '',] def_rho <- def_gp_prior$prior[min(which(def_gp_prior$class == 'lscale'))] if(def_rho == ''){ @@ -500,7 +485,8 @@ prep_gp_covariate = function(data, boundary = boundary, mean = NA, max_dist = covariate_max_dist, - scale = scale) + scale = scale, + initial_setup = TRUE) # Make attributes table byname <- ifelse(is.na(by), '', paste0(':', by)) diff --git a/R/marginaleffects.mvgam.R b/R/marginaleffects.mvgam.R index 3b8d8cd5..061f79e7 100644 --- a/R/marginaleffects.mvgam.R +++ b/R/marginaleffects.mvgam.R @@ -17,79 +17,6 @@ #' @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 ggplot2 theme_get theme_set theme_classic -#' @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 - def_theme <- theme_get() - theme_set(theme_classic(base_size = 12, base_family = 'serif')) - orig_col <- .Options$ggplot2.discrete.colour - orig_fill <- .Options$ggplot2.discrete.fill - orig_cont <- .Options$ggplot2.continuous.colour - options(ggplot2.discrete.colour = c("#B97C7C", - "#A25050", - "#8F2727", - "darkred", - "#630000", - "#300000", - "#170000"), - ggplot2.continuous.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, - ...) - theme_set(def_theme) - orig_col -> .Options$ggplot2.discrete.colour - orig_fill -> .Options$ggplot2.discrete.fill - orig_cont -> .Options$ggplot2.continuous.colour -} - #' Functions needed for working with marginaleffects #' @rdname mvgam_marginaleffects #' @export diff --git a/R/mvgam.R b/R/mvgam.R index 3f5669b3..b3bff3c5 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -1247,7 +1247,8 @@ mvgam = function(formula, dplyr::arrange(time, series) -> X # Matrix of indices in X that correspond to timepoints for each series - ytimes <- matrix(NA, nrow = length(unique(X$time)), ncol = length(unique(X$series))) + ytimes <- matrix(NA, nrow = length(unique(X$time)), + ncol = length(unique(X$series))) for(i in 1:length(unique(X$series))){ ytimes[,i] <- which(X$series == i) } @@ -1661,13 +1662,7 @@ mvgam = function(formula, } # Tidy the representation - clean_up <- vector() - for(x in 1:length(vectorised$model_file)){ - clean_up[x] <- vectorised$model_file[x-1] == "" & - vectorised$model_file[x] == "" - } - clean_up[is.na(clean_up)] <- FALSE - vectorised$model_file <- vectorised$model_file[!clean_up] + vectorised$model_file <- sanitise_modelfile(vectorised$model_file) if(requireNamespace('cmdstanr', quietly = TRUE)){ # Replace new syntax if this is an older version of Stan diff --git a/R/mvgam_setup.R b/R/mvgam_setup.R index 1f8fde02..a6b0572f 100644 --- a/R/mvgam_setup.R +++ b/R/mvgam_setup.R @@ -8,7 +8,7 @@ mvgam_setup <- function(formula, data = list(), na.action, drop.unused.levels = FALSE, - maxit = 40) { + maxit = 5) { if(missing(knots)){ # Initialise the GAM for a few iterations to get all necessary structures for diff --git a/R/plot.mvgam.R b/R/plot.mvgam.R index c889c1ab..1014e0ab 100644 --- a/R/plot.mvgam.R +++ b/R/plot.mvgam.R @@ -133,9 +133,17 @@ plot.mvgam = function(x, type = 'residuals', class = class(object2$mgcv_model$smooth[[x]])[1], mgcv_plottable = object2$mgcv_model$smooth[[x]]$plot.me) })) + + # Filter out any GP terms + if(!is.null(attr(object2$mgcv_model, 'gp_att_table'))){ + gp_names <- unlist(purrr::map(attr(object2$mgcv_model, 'gp_att_table'), 'name')) + smooth_labs %>% + 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", + call. = FALSE) smooth_labs$smooth_index <- 1:NROW(smooth_labs) - if(n_smooths == 0) stop("No terms to plot - nothing for plot.mvgam() to do.") # Leave out random effects and MRF smooths, and any others that are not # considered plottable by mgcv diff --git a/R/plot_effects.R b/R/plot_effects.R new file mode 100644 index 00000000..b56c3e28 --- /dev/null +++ b/R/plot_effects.R @@ -0,0 +1,53 @@ +#' 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(model, + 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 = model, + 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 548b5ffc..a7b3d85f 100644 --- a/R/plot_mvgam_smooth.R +++ b/R/plot_mvgam_smooth.R @@ -115,6 +115,17 @@ plot_mvgam_smooth = function(object, smooth_int <- smooth } + # Check whether this is actually a gp() term + if(!is.null(attr(object2$mgcv_model, 'gp_att_table'))){ + gp_names <- unlist(purrr::map(attr(object2$mgcv_model, 'gp_att_table'), 'name')) + 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', + call. = FALSE) + } + } + # Check whether this type of smooth is even plottable if(!object2$mgcv_model$smooth[[smooth_int]]$plot.me){ stop(paste0('unable to plot ', object2$mgcv_model$smooth[[smooth_int]]$label, @@ -251,8 +262,14 @@ plot_mvgam_smooth = function(object, } # Generate linear predictor matrix from fitted mgcv model - Xp <- obs_Xp_matrix(newdata = pred_dat, - mgcv_model = object2$mgcv_model) + if(trend_effects){ + Xp <- trend_Xp_matrix(newdata = pred_dat, + trend_map = object2$trend_map, + mgcv_model = object2$trend_mgcv_model) + } else { + Xp <- obs_Xp_matrix(newdata = pred_dat, + mgcv_model = object2$mgcv_model) + } # Zero out all other columns in Xp keeps <- object2$mgcv_model$smooth[[smooth_int]]$first.para: @@ -299,8 +316,14 @@ plot_mvgam_smooth = function(object, if(residuals){ # Need to predict from a reduced set that zeroes out all terms apart from the # smooth of interest - Xp2 <- obs_Xp_matrix(newdata = object2$obs_data, - mgcv_model = object2$mgcv_model) + if(trend_effects){ + Xp2 <- trend_Xp_matrix(newdata = object2$obs_data, + trend_map = object2$trend_map, + mgcv_model = object2$trend_mgcv_model) + } else { + Xp2 <- obs_Xp_matrix(newdata = object2$obs_data, + mgcv_model = object2$mgcv_model) + } if(!missing(newdata)){ stop('Partial residual plots not available when using newdata') diff --git a/R/sanitise_modelfile.R b/R/sanitise_modelfile.R new file mode 100644 index 00000000..ae7eadb1 --- /dev/null +++ b/R/sanitise_modelfile.R @@ -0,0 +1,31 @@ +#' Clean up a stan file +#' @noRd +sanitise_modelfile = function(model_file){ + + # Remove empty lines + clean_up <- vector() + for(x in 1:length(model_file)){ + clean_up[x] <- trimws(model_file[x]) == "" | + trimws(model_file[x]) == "NA" + } + clean_up[is.na(clean_up)] <- FALSE + model_file <- model_file[!clean_up] + + # Expand on backslashes to make model more readable + hashes <- vector() + hashes[1] <- FALSE + for(x in 2:length(model_file)){ + hashes[x] <- grepl('//', model_file[x], fixed = TRUE) & + trimws(model_file[x-1]) != "" & + (!grepl('{', model_file[x-1], fixed = TRUE) | + grepl('}', model_file[x-1], fixed = TRUE)) & + !grepl(';', model_file[x], fixed = TRUE) + } + + if(any(hashes)){ + model_file[hashes] <- paste0('\n', + model_file[hashes]) + } + model_file <- readLines(textConnection(model_file), n = -1) + return(model_file) +} diff --git a/R/stan_utils.R b/R/stan_utils.R index 8b78c303..9af1184d 100644 --- a/R/stan_utils.R +++ b/R/stan_utils.R @@ -9,7 +9,7 @@ code = function(object){ stop('argument "object" must be of class "mvgam" or "mvgam_prefit"') } - cat(object$model_file[!grepl('^\\s*$', object$model_file)], sep = '\n') + cat(object$model_file, sep = '\n') } #' @noRd @@ -2571,8 +2571,8 @@ add_trend_predictors = function(trend_formula, 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) + b_trend_lines <- gsub('l_gp', 'l_gp_trend', b_trend_lines) + b_trend_lines <- gsub('k_gp', 'k_gp_trend', 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'), @@ -2614,7 +2614,9 @@ add_trend_predictors = function(trend_formula, } # Check for gp() terms - if(any(grepl('l_gp', trend_model_file))){ + if(any(grepl('l_gp', trend_model_file)) & + any(grepl('k_gp', trend_model_file)) & + any(grepl('z_gp', trend_model_file))){ # Add spd_cov_exp_quad function from brms code if(any(grepl('functions {', model_file, fixed = TRUE))){ @@ -2665,10 +2667,10 @@ add_trend_predictors = function(trend_formula, } 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) + trend_model_file <- gsub('l_gp', 'l_gp_trend', trend_model_file) + trend_model_file <- gsub('k_gp', 'k_gp_trend', 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)) + names(idx_data) <- gsub('l_gp', 'l_gp_trend', names(idx_data)) model_data <- append(model_data, idx_data) l_lines <- grep('// approximate gp eigenvalues', trend_model_file, fixed = TRUE) @@ -2681,7 +2683,7 @@ add_trend_predictors = function(trend_formula, 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)) + names(idx_data) <- gsub('k_gp', 'k_gp_trend', names(idx_data)) model_data <- append(model_data, idx_data) k_lines <- grep('// basis functions for approximate gp', trend_model_file, fixed = TRUE) @@ -2698,7 +2700,7 @@ add_trend_predictors = function(trend_formula, fixed = TRUE) + 1 last <- end for(i in end:(end+50)){ - if(grepl('vector[trend_k_gp', trend_model_file[i], + if(grepl('vector[k_gp_trend', trend_model_file[i], fixed = TRUE)){ last <- i } else { @@ -2732,9 +2734,10 @@ add_trend_predictors = function(trend_formula, trend_model_file) - 1] if(any(grepl('normal(0, lambda', trend_model_file, fixed = TRUE))){ + idx_headers <- trend_model_file[grep('normal(0, lambda', + trend_model_file, fixed = TRUE)-1] spline_coef_headers <- c(spline_coef_headers, - trend_model_file[grep('normal(0, lambda', - trend_model_file, fixed = TRUE)-1]) + grep('//', idx_headers, value = TRUE)) } if(any(grepl('// prior for gp', trend_model_file))){ @@ -2880,7 +2883,6 @@ add_trend_predictors = function(trend_formula, } model_file <- readLines(textConnection(model_file), n = -1) - } # Add any parametric effect beta lines @@ -3141,6 +3143,18 @@ add_trend_predictors = function(trend_formula, model_file <- gsub('latent trend', 'latent state', model_file) + # Any final tidying for trend_level terms + model_file <- gsub('byseriestrend', 'bytrendtrend', model_file) + model_file <- gsub(':seriestrend', ':trendtrend', model_file) + + names(model_data) <- gsub('byseriestrend', 'bytrendtrend', names(model_data)) + names(model_data) <- gsub(':seriestrend', ':trendtrend', names(model_data)) + + names(trend_mvgam$mgcv_model$coefficients) <- + gsub('byseriestrend', 'bytrendtrend', names(trend_mvgam$mgcv_model$coefficients)) + names(trend_mvgam$mgcv_model$coefficients) <- + gsub(':seriestrend', ':trendtrend', names(trend_mvgam$mgcv_model$coefficients)) + return(list(model_file = model_file, model_data = model_data, trend_mgcv_model = trend_mvgam$mgcv_model, diff --git a/man/plot_predictions.mvgam.Rd b/man/plot_effects.mvgam.Rd similarity index 73% rename from man/plot_predictions.mvgam.Rd rename to man/plot_effects.mvgam.Rd index 77f37263..3d0f7833 100644 --- a/man/plot_predictions.mvgam.Rd +++ b/man/plot_effects.mvgam.Rd @@ -1,27 +1,37 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/marginaleffects.mvgam.R -\name{plot_predictions.mvgam} -\alias{plot_predictions.mvgam} +% 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_predictions.mvgam( +plot_effects(object, ...) + +\method{plot_effects}{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{...}{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{model}{Model object} \item{condition}{Conditional predictions @@ -53,22 +63,6 @@ acceptable values is returned in an error message. When \code{type} is \code{NUL 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. @@ -84,17 +78,6 @@ the \code{marginaleffects:::type_dictionary} dataframe.} \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 diff --git a/src/mvgam.dll b/src/mvgam.dll index 3f551697..57254383 100644 Binary files a/src/mvgam.dll and b/src/mvgam.dll differ diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf index 12ae3edb..aed6ffb2 100644 Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ diff --git a/tests/testthat/test-gp.R b/tests/testthat/test-gp.R new file mode 100644 index 00000000..df5947a3 --- /dev/null +++ b/tests/testthat/test-gp.R @@ -0,0 +1,124 @@ +context("gp") + +test_that("gp_to_s is working properly", { + # All true gp() terms should be changed to s() with k = k+1 +formula <- y ~ s(series) + gp(banana) + + infect:you + gp(hardcourt) + +expect_equal(attr(terms(mvgam:::gp_to_s(formula), keep.order = TRUE), + 'term.labels'), + attr(terms(formula(y ~ s(series) + + s(banana, k = 11) + + infect:you + + s(hardcourt, k = 11)), + keep.order = TRUE), + 'term.labels')) + +# Characters that match to 'gp' should not be changed +formula <- y ~ gp(starwars) + s(gp) +expect_equal(attr(terms(mvgam:::gp_to_s(formula), keep.order = TRUE), + 'term.labels'), + attr(terms(formula(y ~ s(starwars, k = 11) + s(gp)), + keep.order = TRUE), + 'term.labels')) +}) + +test_that("gp for observation models working properly", { + mod <- mvgam(y ~ s(series, bs = 're') + + gp(time, by = series) + + year:season, + data = gaus_data$data_train, + family = gaussian(), + run_model = FALSE) + + # Gp data structures should be in the model_data + expect_true("l_gp_time_byseriesseries_1" %in% names(mod$model_data)) + expect_true("b_idx_gp_time_byseriesseries_1" %in% names(mod$model_data)) + expect_true("k_gp_time_byseriesseries_1" %in% names(mod$model_data)) + + # These should match to the eigenvalues and eigenfunctions created by + # a similar brms call + brms_dat <- suppressWarnings(brms::make_standata(y ~ s(series, bs = 're') + + gp(time, by = series, k = 10, + c = 5/4) + + year:season, + data = gaus_data$data_train, + family = gaussian())) + # Eigenvalues should be identical + all.equal(as.vector(brms_dat$slambda_1_1), + mod$model_data$l_gp_time_byseriesseries_1) + + # Eigenfunctions will be nearly identical + row_s1 <- which(gaus_data$data_train$series == 'series_1' & + !is.na(gaus_data$data_train$y)) + col_s1 <- grep('gp(time):seriesseries_1', names(coef(mod$mgcv_model)), + fixed = TRUE) + + expect_true(identical(dim(brms_dat$Xgp_1_1), + dim(mod$model_data$X[row_s1,col_s1]))) + + expect_true(max(abs(brms_dat$Xgp_1_1 - + mod$model_data$X[row_s1,col_s1])) < 0.01) + + # The mgcv model formula should contain s() in place of gp() + expect_equal(attr(terms(mod$mgcv_model$formula, keep.order = TRUE), + 'term.labels'), + attr(terms(formula(y ~ s(time, by = series, k = 11) + + year:season + + s(series, bs = "re")), + keep.order = TRUE), + 'term.labels')) +}) + + +test_that("gp for process models working properly", { + mod <- mvgam(y ~ s(series, bs = 're'), + trend_formula = ~ + gp(time, by = trend) + + year:season, + data = beta_data$data_train, + family = betar(), + trend_model = 'AR1', + run_model = FALSE) + + # Model file should have prior lines for gp terms + expect_true(any(grepl('// prior for gp(time):trendtrend1_trend...', + mod$model_file, fixed = TRUE))) + + expect_true(any(grepl("b_trend[b_trend_idx_gp_time_bytrendtrend1] = sqrt(spd_cov_exp_quad(", + mod$model_file, fixed = TRUE))) + + # Gp data structures should be in the model_data + expect_true("l_gp_trend_time_bytrendtrend1" %in% names(mod$model_data)) + expect_true("b_trend_idx_gp_time_bytrendtrend1" %in% names(mod$model_data)) + expect_true("k_gp_trend_time_bytrendtrend1" %in% names(mod$model_data)) + + # These should match to the eigenvalues and eigenfunctions created by + # a similar brms call + brms_dat <- suppressWarnings(brms::make_standata(y ~ gp(time, by = series, k = 10, + c = 5/4), + data = beta_data$data_train, + family = gaussian())) + # Eigenvalues should be identical + expect_true(all.equal(as.vector(brms_dat$slambda_1_1), + mod$model_data$l_gp_trend_time_bytrendtrend1)) + + # Eigenfunctions will be nearly identical + row_s1 <- mod$model_data$ytimes_trend[,1] + col_s1 <- grep('gp(time):trendtrend1', names(coef(mod$trend_mgcv_model)), + fixed = TRUE) + + expect_true(identical(dim(brms_dat$Xgp_1_1), + dim(mod$model_data$X_trend[row_s1,col_s1]))) + + expect_true(max(abs(brms_dat$Xgp_1_1 - + mod$model_data$X_trend[row_s1,col_s1])) < 0.01) + + # The mgcv model formula should contain s() in place of gp() + expect_equal(attr(terms(mod$trend_mgcv_model$formula, keep.order = TRUE), + 'term.labels'), + attr(terms(formula(y ~ s(time, by = series, k = 11) + + year:season), + keep.order = TRUE), + 'term.labels')) +})