Skip to content

Commit

Permalink
Merge pull request #92 from nicholasjclark/brms_gps
Browse files Browse the repository at this point in the history
update vignettes
nicholasjclark authored Nov 12, 2024
2 parents d723c7a + 8b08485 commit 5b3ff5c
Showing 103 changed files with 1,127 additions and 693 deletions.
32 changes: 19 additions & 13 deletions R/jsdgam.R
Original file line number Diff line number Diff line change
@@ -201,17 +201,16 @@
#' scale_color_viridis_c() +
#' theme_classic()
#'
#' # Inspect default priors for a joint species model with spatial factors
#' # Inspect default priors for a joint species model with three spatial factors
#' priors <- get_mvgam_priors(formula = count ~
#' # Environmental model includes species-level intercepts
#' # and random slopes for a linear effect of temperature
#' species +
#' # Environmental model includes random slopes for
#' # a linear effect of temperature
#' s(species, bs = 're', by = temperature),
#'
#' # Each factor estimates a different nonlinear spatial process, using
#' # 'by = trend' as in other mvgam State-Space models
#' factor_formula = ~ gp(lat, lon, k = 6, by = trend) - 1,
#' n_lv = 4,
#' n_lv = 3,
#'
#' # The data and grouping variables
#' data = dat,
@@ -223,20 +222,27 @@
#' head(priors)
#'
#' # Fit a JSDM that estimates hierarchical temperature responses
#' # and that uses four latent spatial factors
#' # and that uses three latent spatial factors
#' mod <- jsdgam(formula = count ~
#' # Environmental model includes species-level intercepts
#' # and random slopes for a linear effect of temperature
#' species +
#' # Environmental model includes random slopes for a
#' # linear effect of temperature
#' s(species, bs = 're', by = temperature),
#'
#' # Each factor estimates a different nonlinear spatial process, using
#' # 'by = trend' as in other mvgam State-Space models
#' factor_formula = ~ gp(lat, lon, k = 6, by = trend) - 1,
#' n_lv = 4,
#'
#' # Change default priors for fixed effect betas to standard normal
#' priors = prior(std_normal(), class = b),
#' n_lv = 3,
#'
#' # Change default priors for fixed random effect variances and
#' # factor P marginal deviations to standard normal
#' priors = c(prior(std_normal(),
#' class = sigma_raw),
#' prior(std_normal(),
#' class = `alpha_gp_trend(lat, lon):trendtrend1`),
#' prior(std_normal(),
#' class = `alpha_gp_trend(lat, lon):trendtrend2`),
#' prior(std_normal(),
#' class = `alpha_gp_trend(lat, lon):trendtrend3`)),
#'
#' # The data and the grouping variables
#' data = dat,
10 changes: 7 additions & 3 deletions R/plot_mvgam_smooth.R
Original file line number Diff line number Diff line change
@@ -168,13 +168,16 @@ plot_mvgam_smooth = function(object,
if(length(unlist(strsplit(smooth, ','))) >= 2L){

# Use default mgcv plotting for bivariate smooths as it is quicker
object2$mgcv_model <- relabel_gps(object2$mgcv_model)
if(inherits(object2$mgcv_model$smooth[[smooth_int]], 'tprs.smooth') |
inherits(object2$mgcv_model$smooth[[smooth_int]], 't2smooth') |
inherits(object2$mgcv_model$smooth[[smooth_int]], 'tensor.smooth')){
suppressWarnings(plot(object2$mgcv_model, select = smooth_int,
suppressWarnings(plot(object2$mgcv_model,
select = smooth_int,
residuals = residuals,
scheme = 2,
main = '', too.far = 0,
main = '',
too.far = 0,
contour.col = 'black',
hcolors = hcl.colors(25,
palette = 'Reds 2'),
@@ -183,7 +186,8 @@ plot_mvgam_smooth = function(object,
box(col = 'white')
box(bty = 'l', lwd = 2)
} else {
suppressWarnings(plot(object2$mgcv_model, select = smooth_int,
suppressWarnings(plot(object2$mgcv_model,
select = smooth_int,
residuals = residuals,
scheme = 2,
main = '', too.far = 0,
60 changes: 38 additions & 22 deletions R/stan_utils.R
Original file line number Diff line number Diff line change
@@ -3275,7 +3275,7 @@ add_trend_predictors = function(trend_formula,
#' @param quiet Logical (verbose or not?)
#' @details Utility function written by Michael Betancourt (https://betanalpha.github.io/)
#' @noRd
check_div <- function(fit, quiet=FALSE, sampler_params) {
check_div <- function(fit, quiet = FALSE, sampler_params) {
if(missing(sampler_params)){
sampler_params <- rstan::get_sampler_params(fit, inc_warmup=FALSE)
}
@@ -3298,7 +3298,7 @@ check_div <- function(fit, quiet=FALSE, sampler_params) {
#' @param quiet Logical (verbose or not?)
#' @details Utility function written by Michael Betancourt (https://betanalpha.github.io/)
#' @noRd
check_treedepth <- function(fit, max_depth = 10, quiet=FALSE,
check_treedepth <- function(fit, max_depth = 10, quiet = FALSE,
sampler_params) {
if(missing(sampler_params)){
sampler_params <- rstan::get_sampler_params(fit, inc_warmup=FALSE)
@@ -3325,7 +3325,7 @@ check_treedepth <- function(fit, max_depth = 10, quiet=FALSE,
#' @param quiet Logical (verbose or not?)
#' @details Utility function written by Michael Betancourt (https://betanalpha.github.io/)
#' @noRd
check_energy <- function(fit, quiet=FALSE, sampler_params) {
check_energy <- function(fit, quiet = FALSE, sampler_params) {
if(missing(sampler_params)){
sampler_params <- rstan::get_sampler_params(fit, inc_warmup=FALSE)
}
@@ -3359,29 +3359,37 @@ check_n_eff <- function(fit, quiet=FALSE, fit_summary, ignore_b_trend = FALSE) {
fit_summary <- rstan::summary(fit, probs = c(0.5))$summary
}

fit_summary <- fit_summary[-grep('ypred',
rownames(fit_summary)), ]

if(any(grep('LV', rownames(fit_summary)))){
fit_summary <- fit_summary[-grep('LV', rownames(fit_summary)), ]
fit_summary <- fit_summary[-grep('lv_coefs', rownames(fit_summary)), ]
fit_summary <- fit_summary[-grep('LV',
rownames(fit_summary)), ]
fit_summary <- fit_summary[-grep('lv_coefs',
rownames(fit_summary)), ]

if(any(grepl('L[', rownames(fit_summary), fixed = TRUE))){
fit_summary <- fit_summary[-grep('L[', rownames(fit_summary), fixed = TRUE), ]
fit_summary <- fit_summary[-grep('L[',
rownames(fit_summary), fixed = TRUE), ]
}
if(any(grepl('LV_raw[', rownames(fit_summary), fixed = TRUE))){
fit_summary <- fit_summary[-grep('LV_raw[', rownames(fit_summary), fixed = TRUE), ]
fit_summary <- fit_summary[-grep('LV_raw[',
rownames(fit_summary), fixed = TRUE), ]
}
}

if(ignore_b_trend){
if(any(grepl('b_trend[', rownames(fit_summary), fixed = TRUE))){
fit_summary <- fit_summary[-grep('b_trend[', rownames(fit_summary), fixed = TRUE), ]
if(any(grepl('_trend', rownames(fit_summary), fixed = TRUE))){
fit_summary <- fit_summary[-grep('_trend',
rownames(fit_summary), fixed = TRUE), ]
}

if(any(grepl('trend_mus[', rownames(fit_summary), fixed = TRUE))){
fit_summary <- fit_summary[-grep('trend_mus[', rownames(fit_summary), fixed = TRUE), ]
fit_summary <- fit_summary[-grep('trend_mus[',
rownames(fit_summary), fixed = TRUE), ]
}
}

N <- dim(fit_summary)[[1]]

iter <- dim(rstan::extract(fit)[[1]])[[1]]

neffs <- fit_summary[,'n_eff']
@@ -3406,34 +3414,42 @@ check_n_eff <- function(fit, quiet=FALSE, fit_summary, ignore_b_trend = FALSE) {
#' @param quiet Logical (verbose or not?)
#' @details Utility function written by Michael Betancourt (https://betanalpha.github.io/)
#' @noRd
check_rhat <- function(fit, quiet=FALSE, fit_summary, ignore_b_trend = FALSE) {
check_rhat <- function(fit, quiet = FALSE, fit_summary, ignore_b_trend = FALSE) {
if(missing(fit_summary)){
fit_summary <- rstan::summary(fit, probs = c(0.5))$summary
}

fit_summary <- fit_summary[-grep('ypred',
rownames(fit_summary)), ]

if(any(grep('LV', rownames(fit_summary)))){
fit_summary <- fit_summary[-grep('LV', rownames(fit_summary)), ]
fit_summary <- fit_summary[-grep('lv_coefs', rownames(fit_summary)), ]
fit_summary <- fit_summary[-grep('LV',
rownames(fit_summary)), ]
fit_summary <- fit_summary[-grep('lv_coefs',
rownames(fit_summary)), ]

if(any(grepl('L[', rownames(fit_summary), fixed = TRUE))){
fit_summary <- fit_summary[-grep('L[', rownames(fit_summary), fixed = TRUE), ]
fit_summary <- fit_summary[-grep('L[',
rownames(fit_summary), fixed = TRUE), ]
}
if(any(grepl('LV_raw[', rownames(fit_summary), fixed = TRUE))){
fit_summary <- fit_summary[-grep('LV_raw[', rownames(fit_summary), fixed = TRUE), ]
fit_summary <- fit_summary[-grep('LV_raw[',
rownames(fit_summary), fixed = TRUE), ]
}
}

if(ignore_b_trend){
if(any(grepl('b_trend[', rownames(fit_summary), fixed = TRUE))){
fit_summary <- fit_summary[-grep('b_trend[', rownames(fit_summary), fixed = TRUE), ]
if(any(grepl('_trend', rownames(fit_summary), fixed = TRUE))){
fit_summary <- fit_summary[-grep('_trend',
rownames(fit_summary), fixed = TRUE), ]
}

if(any(grepl('trend_mus[', rownames(fit_summary), fixed = TRUE))){
fit_summary <- fit_summary[-grep('trend_mus[', rownames(fit_summary), fixed = TRUE), ]
fit_summary <- fit_summary[-grep('trend_mus[',
rownames(fit_summary), fixed = TRUE), ]
}
}

N <- dim(fit_summary)[[1]]

no_warning <- TRUE
rhats <- fit_summary[,'Rhat']
if(max(rhats, na.rm = TRUE) > 1.05) no_warning <- FALSE
55 changes: 17 additions & 38 deletions R/summary.mvgam.R
Original file line number Diff line number Diff line change
@@ -929,45 +929,24 @@ gp_param_summary = function(object,

# Determine which parameters to extract
if(trend_effects){
alpha_params <- gsub('series',
'trend',
gsub('gp_',
'gp_trend_',
gsub(':',
'by',
gsub(')',
'_',
gsub('(', '_', paste0('alpha_', gp_names),
fixed = TRUE),
fixed = TRUE))))
rho_params <- gsub('series',
'trend',
gsub('gp_',
'gp_trend_',
gsub(':', 'by',
gsub(')',
'_',
gsub('(', '_', paste0('rho_', gp_names),
fixed = TRUE), fixed = TRUE))))

alpha_params <- gsub('gp_',
'gp_trend_',
gsub('series',
'trend',
paste0('alpha_', clean_gpnames(gp_names)),
fixed = TRUE),
fixed = TRUE)
rho_params <- gsub('gp_',
'gp_trend_',
gsub('series',
'trend',
paste0('rho_', clean_gpnames(gp_names)),
fixed = TRUE),
fixed = TRUE)
} else {
alpha_params <- gsub(':',
'by', gsub(')',
'_',
gsub('(',
'_',
paste0('alpha_',
clean_gpnames(gp_names)),
fixed = TRUE),
fixed = TRUE))
rho_params <- gsub(':',
'by', gsub(')',
'_',
gsub('(',
'_',
paste0('rho_',
clean_gpnames(gp_names)),
fixed = TRUE),
fixed = TRUE))
alpha_params <- paste0('alpha_', clean_gpnames(gp_names))
rho_params <- paste0('rho_', clean_gpnames(gp_names))
}

# Create summary tables
2 changes: 1 addition & 1 deletion README.Rmd
Original file line number Diff line number Diff line change
@@ -289,5 +289,5 @@ There are many more extended uses of `mvgam`, including the ability to fit hiera
This project is licensed under an `MIT` open source license

## Interested in contributing?
I'm actively seeking PhD students and other researchers to work in the areas of ecological forecasting, multivariate model evaluation and development of `mvgam`. Please reach out if you are interested (n.clark'at'uq.edu.au). Other contributions are also very welcome, but please see [The Contributor Instructions](https://github.com/nicholasjclark/mvgam/blob/main/.github/CONTRIBUTING.md) for general guidelines. Note that
I'm actively seeking PhD students and other researchers to work in the areas of ecological forecasting, multivariate model evaluation and development of `mvgam`. Please reach out if you are interested (n.clark'at'uq.edu.au). Other contributions are also very welcome, but please see [The Contributor Instructions](https://github.com/nicholasjclark/mvgam/blob/master/.github/CONTRIBUTING.md) for general guidelines. Note that
by participating in this project you agree to abide by the terms of its [Contributor Code of Conduct](https://dplyr.tidyverse.org/CODE_OF_CONDUCT).
38 changes: 19 additions & 19 deletions README.md
Original file line number Diff line number Diff line change
@@ -238,28 +238,28 @@ summary(lynx_mvgam)
#>
#> GAM coefficient (beta) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> (Intercept) 6.400 6.60 6.900 1.00 796
#> s(season).1 -0.650 -0.13 0.380 1.00 855
#> s(season).2 0.740 1.30 1.900 1.00 933
#> s(season).3 1.200 1.90 2.500 1.01 633
#> s(season).4 -0.046 0.55 1.100 1.00 834
#> s(season).5 -1.300 -0.69 -0.062 1.00 995
#> s(season).6 -1.300 -0.55 0.110 1.01 1006
#> s(season).7 0.021 0.72 1.400 1.00 757
#> s(season).8 0.620 1.40 2.100 1.00 763
#> s(season).9 -0.330 0.23 0.820 1.00 772
#> s(season).10 -1.300 -0.86 -0.400 1.00 1174
#> (Intercept) 6.400 6.60 6.900 1 942
#> s(season).1 -0.640 -0.13 0.400 1 1123
#> s(season).2 0.710 1.30 1.900 1 998
#> s(season).3 1.300 1.90 2.500 1 912
#> s(season).4 -0.045 0.52 1.200 1 856
#> s(season).5 -1.300 -0.70 -0.034 1 933
#> s(season).6 -1.200 -0.54 0.150 1 1147
#> s(season).7 0.062 0.73 1.500 1 928
#> s(season).8 0.610 1.40 2.100 1 1016
#> s(season).9 -0.370 0.21 0.820 1 936
#> s(season).10 -1.400 -0.87 -0.360 1 1117
#>
#> Approximate significance of GAM smooths:
#> edf Ref.df Chi.sq p-value
#> s(season) 9.99 10 46.5 <2e-16 ***
#> edf Ref.df Chi.sq p-value
#> s(season) 9.9 10 64.4 1.7e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Latent trend parameter AR estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> ar1[1] 0.60 0.83 0.98 1 608
#> sigma[1] 0.38 0.48 0.61 1 636
#> ar1[1] 0.60 0.83 0.98 1.01 643
#> sigma[1] 0.39 0.48 0.62 1.00 821
#>
#> Stan MCMC diagnostics:
#> n_eff / iter looks reasonable for all parameters
@@ -268,7 +268,7 @@ summary(lynx_mvgam)
#> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
#> E-FMI indicated no pathological behavior
#>
#> Samples were drawn using NUTS(diag_e) at Tue Nov 12 8:48:44 AM 2024.
#> Samples were drawn using NUTS(diag_e) at Tue Nov 12 10:11:54 AM 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split MCMC chains
#> (at convergence, Rhat = 1)
@@ -409,7 +409,7 @@ plot(lynx_mvgam, type = 'forecast', newdata = lynx_test)
<img src="man/figures/README-unnamed-chunk-21-1.png" alt="Plotting forecast distributions using mvgam in R" width="60%" style="display: block; margin: auto;" />

#> Out of sample DRPS:
#> 2380.298498
#> 2384.82381825

And the estimated latent trend component, again using the more flexible
`plot_mvgam_...()` option to show first derivatives of the estimated
@@ -568,7 +568,7 @@ summary(mod, include_betas = FALSE)
#> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
#> E-FMI indicated no pathological behavior
#>
#> Samples were drawn using NUTS(diag_e) at Tue Nov 12 8:49:29 AM 2024.
#> Samples were drawn using NUTS(diag_e) at Tue Nov 12 10:12:39 AM 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split MCMC chains
#> (at convergence, Rhat = 1)
@@ -608,7 +608,7 @@ areas of ecological forecasting, multivariate model evaluation and
development of `mvgam`. Please reach out if you are interested
(n.clark’at’uq.edu.au). Other contributions are also very welcome, but
please see [The Contributor
Instructions](https://github.com/nicholasjclark/mvgam/blob/main/.github/CONTRIBUTING.md)
Instructions](https://github.com/nicholasjclark/mvgam/blob/master/.github/CONTRIBUTING.md)
for general guidelines. Note that by participating in this project you
agree to abide by the terms of its [Contributor Code of
Conduct](https://dplyr.tidyverse.org/CODE_OF_CONDUCT).
Loading

0 comments on commit 5b3ff5c

Please sign in to comment.