Brms gps

.github/CODE_OF_CONDUCT.md

@@ -0,0 +1,126 @@
+# Contributor Covenant Code of Conduct
+
+## Our Pledge
+
+We as members, contributors, and leaders pledge to make participation in our
+community a harassment-free experience for everyone, regardless of age, body
+size, visible or invisible disability, ethnicity, sex characteristics, gender
+identity and expression, level of experience, education, socio-economic status,
+nationality, personal appearance, race, caste, color, religion, or sexual
+identity and orientation.
+
+We pledge to act and interact in ways that contribute to an open, welcoming,
+diverse, inclusive, and healthy community.
+
+## Our Standards
+
+Examples of behavior that contributes to a positive environment for our
+community include:
+
+* Demonstrating empathy and kindness toward other people
+* Being respectful of differing opinions, viewpoints, and experiences
+* Giving and gracefully accepting constructive feedback
+* Accepting responsibility and apologizing to those affected by our mistakes,
+  and learning from the experience
+* Focusing on what is best not just for us as individuals, but for the overall
+  community
+
+Examples of unacceptable behavior include:
+
+* The use of sexualized language or imagery, and sexual attention or advances of
+  any kind
+* Trolling, insulting or derogatory comments, and personal or political attacks
+* Public or private harassment
+* Publishing others' private information, such as a physical or email address,
+  without their explicit permission
+* Other conduct which could reasonably be considered inappropriate in a
+  professional setting
+
+## Enforcement Responsibilities
+
+Community leaders are responsible for clarifying and enforcing our standards of
+acceptable behavior and will take appropriate and fair corrective action in
+response to any behavior that they deem inappropriate, threatening, offensive,
+or harmful.
+
+Community leaders have the right and responsibility to remove, edit, or reject
+comments, commits, code, wiki edits, issues, and other contributions that are
+not aligned to this Code of Conduct, and will communicate reasons for moderation
+decisions when appropriate.
+
+## Scope
+
+This Code of Conduct applies within all community spaces, and also applies when
+an individual is officially representing the community in public spaces.
+Examples of representing our community include using an official e-mail address,
+posting via an official social media account, or acting as an appointed
+representative at an online or offline event.
+
+## Enforcement
+
+Instances of abusive, harassing, or otherwise unacceptable behavior may be
+reported to the community leaders responsible for enforcement at [email protected]. 
+All complaints will be reviewed and investigated promptly and fairly.
+
+All community leaders are obligated to respect the privacy and security of the
+reporter of any incident.
+
+## Enforcement Guidelines
+
+Community leaders will follow these Community Impact Guidelines in determining
+the consequences for any action they deem in violation of this Code of Conduct:
+
+### 1. Correction
+
+**Community Impact**: Use of inappropriate language or other behavior deemed
+unprofessional or unwelcome in the community.
+
+**Consequence**: A private, written warning from community leaders, providing
+clarity around the nature of the violation and an explanation of why the
+behavior was inappropriate. A public apology may be requested.
+
+### 2. Warning
+
+**Community Impact**: A violation through a single incident or series of
+actions.
+
+**Consequence**: A warning with consequences for continued behavior. No
+interaction with the people involved, including unsolicited interaction with
+those enforcing the Code of Conduct, for a specified period of time. This
+includes avoiding interactions in community spaces as well as external channels
+like social media. Violating these terms may lead to a temporary or permanent
+ban.
+
+### 3. Temporary Ban
+
+**Community Impact**: A serious violation of community standards, including
+sustained inappropriate behavior.
+
+**Consequence**: A temporary ban from any sort of interaction or public
+communication with the community for a specified period of time. No public or
+private interaction with the people involved, including unsolicited interaction
+with those enforcing the Code of Conduct, is allowed during this period.
+Violating these terms may lead to a permanent ban.
+
+### 4. Permanent Ban
+
+**Community Impact**: Demonstrating a pattern of violation of community
+standards, including sustained inappropriate behavior, harassment of an
+individual, or aggression toward or disparagement of classes of individuals.
+
+**Consequence**: A permanent ban from any sort of public interaction within the
+community.
+
+## Attribution
+
+This Code of Conduct is adapted from the [Contributor Covenant][homepage],
+version 2.1, available at
+<https://www.contributor-covenant.org/version/2/1/code_of_conduct.html>.
+
+Community Impact Guidelines were inspired by
+[Mozilla's code of conduct enforcement ladder][https://github.com/mozilla/inclusion].
+
+For answers to common questions about this code of conduct, see the FAQ at
+<https://www.contributor-covenant.org/faq>. Translations are available at <https://www.contributor-covenant.org/translations>.
+
+[homepage]: https://www.contributor-covenant.org

.github/CONTRIBUTING.md

@@ -0,0 +1,51 @@
+# Contributing to mvgam
+
+This document outlines how to propose a change to mvgam.
+For a detailed discussion on contributing to this and other open source R packages, please see the [development contributing guide](https://rstd.io/tidy-contrib) and our [code review principles](https://code-review.tidyverse.org/).
+
+## Fixing typos
+
+You can fix typos, spelling mistakes, or grammatical errors in the documentation directly using the GitHub web interface, as long as the changes are made in the _source_ file. 
+This generally means you'll need to edit [roxygen2 comments](https://roxygen2.r-lib.org/articles/roxygen2.html) in an `.R`, not a `.Rd` file. 
+You can find the `.R` file that generates the `.Rd` by reading the comment in the first line.
+
+## Bigger changes
+
+If you want to make a bigger change, it's a good idea to first file an issue and make sure someone from the team agrees that it’s needed. 
+If you’ve found a bug, please file an issue that illustrates the bug with a minimal 
+[reprex](https://www.tidyverse.org/help/#reprex) (this will also help you write a unit test, if needed).
+See the tidyverse guide on [how to create a great issue](https://code-review.tidyverse.org/issues/) for more advice.
+
+### Pull request process
+
+*   Fork the package and clone onto your computer. If you haven't done this before, we recommend using `usethis::create_from_github("nicholasjclark/mvgam", fork = TRUE)`.
+
+*   Install all development dependencies with `devtools::install_dev_deps()`, and then make sure the package passes R CMD check by running `devtools::check()`. 
+    If R CMD check doesn't pass cleanly, it's a good idea to ask for help before continuing. 
+*   Create a Git branch for your pull request (PR). We recommend using `usethis::pr_init("brief-description-of-change")`.
+
+*   Make your changes, commit to git, and then create a PR by running `usethis::pr_push()`, and following the prompts in your browser.
+    The title of your PR should briefly describe the change.
+    The body of your PR should contain `Fixes #issue-number`.
+
+*  For user-facing changes, add a bullet to the top of `NEWS.md` (i.e. just below the first header). Follow the style described in <https://style.tidyverse.org/news.html>.
+
+### Code style
+
+*   New code should follow the tidyverse [style guide](https://style.tidyverse.org) where possible. 
+    You can use the [styler](https://CRAN.R-project.org/package=styler) package to apply these styles, but please don't restyle code that has nothing to do with your PR.  
+
+*  We use [roxygen2](https://cran.r-project.org/package=roxygen2), with [Markdown syntax](https://cran.r-project.org/web/packages/roxygen2/vignettes/rd-formatting.html), for documentation.  
+
+*  We use [testthat](https://cran.r-project.org/package=testthat) for unit tests. 
+   Contributions with test cases included are easier to accept.  
+
+## Code of Conduct
+
+Please note that the mvgam project is released with a
+[Contributor Code of Conduct](CODE_OF_CONDUCT.md). By contributing to this
+project you agree to abide by its terms.
+
+## Roadmap
+
+The mvgam package is in a stable state of development, with some degree of active subsequent development as envisioned by the primary authors.

NAMESPACE

@@ -173,6 +173,8 @@ importFrom(bayesplot,nuts_params)
 importFrom(bayesplot,pp_check)
 importFrom(brms,bernoulli)
 importFrom(brms,beta_binomial)
+importFrom(brms,brm)
+importFrom(brms,brmsterms)
 importFrom(brms,conditional_effects)
 importFrom(brms,dbeta_binomial)
 importFrom(brms,do_call)

NEWS.md

@@ -1,5 +1,6 @@
 # mvgam 1.1.4 (development version; not yet on CRAN)
 ## New functionalities
+* Added support for approximate `gp()` effects with more than one covariate and with different kernel functions (#79) 
 * Added function `jsdgam()` to estimate Joint Species Distribution Models in which both the latent factors and the observation model components can include any of mvgam's complex linear predictor effects. Also added a function `residual_cor()` to compute residual correlation, covariance and precision matrices from `jsdgam` models. See `?mvgam::jsdgam` and `?mvgam::residual_cor` for details
 * Added a `stability.mvgam()` method to compute stability metrics from models fit with Vector Autoregressive dynamics (#21 and #76)
 * Added functionality to estimate hierarchical error correlations when using multivariate latent process models and when the data are nested among levels of a relevant grouping factor (#75); see `?mvgam::AR` for an example

R/conditional_effects.R

@@ -98,8 +98,8 @@
 conditional_effects.mvgam = function(x,
                                      effects = NULL,
                                      type = 'response',
-                                     points = TRUE,
-                                     rug = TRUE,
+                                     points = FALSE,
+                                     rug = FALSE,
                                      ...){
 
   use_def_effects <- is.null(effects)

R/get_linear_predictors.R

@@ -1,4 +1,5 @@
 #' Function to prepare observation model linear predictor matrix
+#' @importFrom brms brmsterms
 #' @noRd
 obs_Xp_matrix = function(newdata, mgcv_model){
   suppressWarnings(Xp  <- try(predict(mgcv_model,
@@ -32,41 +33,41 @@ obs_Xp_matrix = function(newdata, mgcv_model){
   # 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
+
+    # Compute the gp() eigenfunctions for newdata using the supplied brms_mock object
+    # Requires a dataframe of all relevant variables for the gp effects
+    mock_terms <- brms::brmsterms(attr(mgcv_model, 'brms_mock')$formula)
+    terms_needed <- unique(all.vars(mock_terms$formula)[-1])
+    newdata_mock <- data.frame(newdata[[terms_needed[1]]])
+    if(length(terms_needed) > 1L){
+      for(i in 2:length(terms_needed)){
+        newdata_mock <- cbind(newdata_mock,
+                              data.frame(newdata[[terms_needed[i]]]))
+      }
+    }
+    colnames(newdata_mock) <- terms_needed
+    newdata_mock$.fake_gp_y <- rnorm(NROW(newdata_mock))
+    brms_mock_data <- brms::standata(attr(mgcv_model, 'brms_mock'),
+                                     newdata = newdata_mock,
+                                     internal = TRUE)
 
     # 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])
-    })
+    bys <- unlist(purrr::map(gp_att_table, 'by'), use.names = FALSE)
+    lvls <- unlist(purrr::map(gp_att_table, 'level'), use.names = FALSE)
+
+    # Extract eigenfunctions for each gp effect
+    eigenfuncs <- eigenfunc_list(stan_data = brms_mock_data,
+                                 mock_df = newdata_mock,
+                                 by = bys,
+                                 level = lvls)
 
     # 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]]
+      Xp[,c(starts[[i]]:ends[[i]])] <- eigenfuncs[[i]]
     }
   }
 
@@ -127,7 +128,7 @@ trend_Xp_matrix = function(newdata, trend_map, series = 'all',
                                     trend_map,
                                     forecast = forecast)
 
-  suppressWarnings(Xp_trend  <- try(predict(mgcv_model,
+  suppressWarnings(Xp_trend <- try(predict(mgcv_model,
                                             newdata = trend_test,
                                             type = 'lpmatrix'),
                                     silent = TRUE))
@@ -158,41 +159,41 @@ trend_Xp_matrix = function(newdata, trend_map, series = 'all',
   # 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
+
+    # Compute the gp() eigenfunctions for newdata using the supplied brms_mock object
+    # Requires a dataframe of all relevant variables for the gp effects
+    mock_terms <- brms::brmsterms(attr(mgcv_model, 'brms_mock')$formula)
+    terms_needed <- unique(all.vars(mock_terms$formula)[-1])
+    newdata_mock <- data.frame(trend_test[[terms_needed[1]]])
+    if(length(terms_needed) > 1L){
+      for(i in 2:length(terms_needed)){
+        newdata_mock <- cbind(newdata_mock,
+                              data.frame(trend_test[[terms_needed[i]]]))
+      }
+    }
+    colnames(newdata_mock) <- terms_needed
+    newdata_mock$.fake_gp_y <- rnorm(NROW(newdata_mock))
+    brms_mock_data <- brms::standata(attr(mgcv_model, 'brms_mock'),
+                                     newdata = newdata_mock,
+                                     internal = TRUE)
 
     # 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 = trend_test,
-                          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])
-    })
+    bys <- unlist(purrr::map(gp_att_table, 'by'), use.names = FALSE)
+    lvls <- unlist(purrr::map(gp_att_table, 'level'), use.names = FALSE)
+
+    # Extract eigenfunctions for each gp effect
+    eigenfuncs <- eigenfunc_list(stan_data = brms_mock_data,
+                                 mock_df = newdata_mock,
+                                 by = bys,
+                                 level = lvls)
 
     # 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]]
+      Xp_trend[,c(starts[[i]]:ends[[i]])] <- eigenfuncs[[i]]
     }
   }