diff --git a/DESCRIPTION b/DESCRIPTION
index e1c0b94b..f99585ec 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -27,6 +27,7 @@ Imports:
posterior (>= 1.0.0),
loo (>= 2.3.1),
parallel,
+ collapse,
pbapply,
tweedie,
MASS,
diff --git a/R/as.data.frame.mvgam.R b/R/as.data.frame.mvgam.R
index e87ff1cc..9ccd26b9 100644
--- a/R/as.data.frame.mvgam.R
+++ b/R/as.data.frame.mvgam.R
@@ -54,8 +54,10 @@ as.data.frame.mvgam = function(x,
# Get a string of all possible variables to extract
all_vars <- variables(x)
- all_orig_vars <- unique(unlist(purrr::map(all_vars, 'orig_name')))
- all_alias_vars <- unique(unlist(purrr::map(all_vars, 'alias')))
+ all_orig_vars <- unlist(purrr::map(all_vars, 'orig_name'))
+ all_alias_vars <- unlist(purrr::map(all_vars, 'alias'))
+
+ all_orig_walias <- all_orig_vars[!is.na(all_alias_vars)]
all_alias_vars <- all_alias_vars[!is.na(all_alias_vars)]
# All possible var sets to extract
@@ -191,7 +193,7 @@ as.data.frame.mvgam = function(x,
names_to_use[[i]] <- NA
} else {
if(any(all_alias_vars == variable[i])){
- vars_to_extract[[i]] <- unname(unlist(purrr::map(all_vars, 'orig_name'))[
+ vars_to_extract[[i]] <- unname(all_orig_walias[
which(all_alias_vars == variable[i])])
names_to_use[[i]] <- variable[i]
} else {
diff --git a/R/globals.R b/R/globals.R
index 38cbceb8..33af3031 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -8,4 +8,4 @@ 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", "trend_series", "trend_y", "."))
diff --git a/R/plot.mvgam.R b/R/plot.mvgam.R
index 0165bb0e..69cebb65 100644
--- a/R/plot.mvgam.R
+++ b/R/plot.mvgam.R
@@ -160,8 +160,6 @@ plot.mvgam = function(x, type = 'residuals',
n_plots <- n_smooths
if (n_plots==0) stop("No suitable terms to plot - plot.mvgam() only handles smooths of 2 or fewer dimensions.")
pages <- 1
- .pardefault <- par(no.readonly=T)
- par(.pardefault)
if (n_plots > 4) pages <- 2
if (n_plots > 8) pages <- 3
@@ -179,7 +177,13 @@ plot.mvgam = function(x, type = 'residuals',
if (c<1) r <- c <- 1
if (c*r < ppp) c <- c + 1
if (c*r < ppp) r <- r + 1
- oldpar<-par(mfrow=c(r,c))
+
+ .pardefault <- par(no.readonly=T)
+ par(.pardefault)
+ oldpar<-par(mfrow=c(r,c),
+ mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(1.5, 0.5, 0))
} else { ppp<-1;oldpar<-par()}
diff --git a/R/plot_mvgam_factors.R b/R/plot_mvgam_factors.R
index 414c6c92..baa0061f 100644
--- a/R/plot_mvgam_factors.R
+++ b/R/plot_mvgam_factors.R
@@ -46,11 +46,20 @@ plot_mvgam_factors = function(object, plot = TRUE){
.pardefault <- par(no.readonly=T)
par(.pardefault)
if(object$n_lv <= 2){
- par(mfrow = c(1, 2))
+ par(mfrow = c(1, 2),
+ mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(1.5, 0.5, 0))
} else if(object$n_lv <= 4){
- par(mfrow = c(2, 2))
+ par(mfrow = c(2, 2),
+ mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(1.5, 0.5, 0))
} else {
- par(mfrow = c(3, 2))
+ par(mfrow = c(3, 2),
+ mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(1.5, 0.5, 0))
}
}
diff --git a/R/plot_mvgam_pterms.R b/R/plot_mvgam_pterms.R
index 745f5f2d..b59716df 100644
--- a/R/plot_mvgam_pterms.R
+++ b/R/plot_mvgam_pterms.R
@@ -55,11 +55,17 @@ if(length(pterms) > 0){
}
if(length(pterms) == 2){
- par(mfrow = c(2, 1))
+ par(mfrow = c(2, 1),
+ mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(1.5, 0.5, 0))
}
if(length(pterms) %in% c(3, 4)){
- par(mfrow = c(2, 2))
+ par(mfrow = c(2, 2),
+ mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(1.5, 0.5, 0))
}
for(i in 1:length(pterms)){
diff --git a/R/plot_mvgam_resids.R b/R/plot_mvgam_resids.R
index 53359b3d..b6684081 100644
--- a/R/plot_mvgam_resids.R
+++ b/R/plot_mvgam_resids.R
@@ -6,7 +6,6 @@
#'@importFrom stats complete.cases qqnorm qqline acf pacf na.pass
#'@param object \code{list} object returned from \code{mvgam}
#'@param series \code{integer} specifying which series in the set is to be plotted
-#'@param n_bins \code{integer} specifying the number of bins to use for binning fitted values
#'@param newdata Optional \code{dataframe} or \code{list} of test data containing at least 'series', 'y', and 'time'
#'in addition to any other variables included in the linear predictor of \code{formula}. If included, the
#'covariate information in \code{newdata} will be used to generate forecasts from the fitted model equations. If
@@ -25,7 +24,7 @@
#'Note, all plots use posterior medians of fitted values / residuals, so uncertainty is not represented.
#'@return A series of base \code{R} plots
#'@export
-plot_mvgam_resids = function(object, series = 1, n_bins = 15,
+plot_mvgam_resids = function(object, series = 1,
newdata, data_test){
# Check arguments
@@ -33,15 +32,7 @@ plot_mvgam_resids = function(object, series = 1, n_bins = 15,
stop('argument "object" must be of class "mvgam"')
}
- if(sign(series) != 1){
- stop('argument "series" must be a positive integer',
- call. = FALSE)
- } else {
- if(series%%1 != 0){
- stop('argument "series" must be a positive integer',
- call. = FALSE)
- }
- }
+ validate_pos_integer(series)
if(series > NCOL(object$ytimes)){
stop(paste0('object only contains data / predictions for ',
@@ -49,22 +40,12 @@ plot_mvgam_resids = function(object, series = 1, n_bins = 15,
call. = FALSE)
}
- if(sign(n_bins) != 1){
- stop('argument "n_bins" must be a positive integer',
- call. = FALSE)
- } else {
- if(n_bins%%1 != 0){
- stop('argument "n_bins" must be a positive integer',
- call. = FALSE)
- }
- }
-
if(!missing("newdata")){
data_test <- newdata
}
# Plotting colours
-probs = c(0.05, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.95)
+probs <- c(0.05, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.95)
c_light <- c("#DCBCBC")
c_light_highlight <- c("#C79999")
c_mid <- c("#B97C7C")
@@ -219,138 +200,61 @@ if(missing(data_test)){
}
}
-bottom_preds <- apply(preds, 2, function(x) quantile(x, 0.05, na.rm = TRUE))
-lower_preds <- apply(preds, 2, function(x) quantile(x, 0.1, na.rm = TRUE))
-median_preds <- apply(preds, 2, function(x) quantile(x, 0.5, na.rm = TRUE))
-upper_preds <- apply(preds, 2, function(x) quantile(x, 0.9, na.rm = TRUE))
-top_preds <- apply(preds, 2, function(x) quantile(x, 0.95, na.rm = TRUE))
-
# Graphical parameters
+.pardefault <- par(no.readonly=T)
+par(.pardefault)
layout(matrix(1:4, ncol = 2, nrow = 2, byrow = TRUE))
+oldpar <- par(mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(2, 0.5, 0))
-# Fitted vs redisuals plot
-n_fitted_bins = n_bins
-
-# Get x-axis values and bin if necessary to prevent overplotting
-sorted_x <- sort(unique(round(c(bottom_preds,
- lower_preds,
- median_preds,
- upper_preds,
- top_preds), 6)))
-sorted_x <- seq(from = min(sorted_x),
- to = max(sorted_x),
- length.out = length(sorted_x))
-
-if(length(sorted_x) > n_fitted_bins){
- sorted_x <- seq(min(sorted_x), max(sorted_x), length.out = n_fitted_bins)
- resid_probs <- do.call(rbind, lapply(seq_along(sorted_x), function(i){
- if(i == 1){
- quantile(series_residuals[which(preds <= sorted_x[i])], probs = probs, na.rm = TRUE)
- } else if(i == length(sorted_x)){
- quantile(series_residuals[which(preds >= sorted_x[i])], probs = probs, na.rm = TRUE)
- } else {
- quantile(series_residuals[which(preds > sorted_x[i - 1] &
- preds <= sorted_x[i])], probs = probs, na.rm = TRUE)
- }
- }))
-
-} else {
- resid_probs <- do.call(rbind, lapply(seq_along(sorted_x), function(i){
- if(i == 1){
- quantile(series_residuals[which(preds <= sorted_x[i])], probs = probs, na.rm = TRUE)
- } else if(i == length(sorted_x)){
- quantile(series_residuals[which(preds >= sorted_x[i])], probs = probs, na.rm = TRUE)
- } else {
- quantile(series_residuals[which(preds > sorted_x[i - 1] &
- preds <= sorted_x[i])], probs = probs, na.rm = TRUE)
- }
- }))
-}
-
-# Get polygon coordinates and plot
-N <- length(sorted_x)
-idx <- rep(1:N, each = 2)
-repped_x <- rep(sorted_x, each = 2)
-
-x <- sapply(1:length(idx),
- function(k) if(k %% 2 == 0)
- repped_x[k] + min(diff(sorted_x))/2 else
- repped_x[k] - min(diff(sorted_x))/2)
-
-# Plot
-# plot(median_preds[1:length(series_residuals)],
-# xlim = range(sorted_x),
-# series_residuals,
-# bty = 'L',
-# xlab = 'Fitted values',
-# ylab = 'DS residuals',
-# pch = 16,
-# col = 'white',
-# cex = 1,
-# ylim = range(resid_probs, na.rm = T))
-# title('Resids vs Fitted Values', line = 0)
-#
-# rect(xleft = x[seq(1, N*2, by = 2)],
-# xright = x[seq(2, N*2, by = 2)],
-# ytop = resid_probs[,9],
-# ybottom = resid_probs[,1],
-# col = c_light,
-# border = 'transparent')
-# rect(xleft = x[seq(1, N*2, by = 2)],
-# xright = x[seq(2, N*2, by = 2)],
-# ytop = resid_probs[,8],
-# ybottom = resid_probs[,2],
-# col = c_light_highlight,
-# border = 'transparent')
-# rect(xleft = x[seq(1, N*2, by = 2)],
-# xright = x[seq(2, N*2, by = 2)],
-# ytop = resid_probs[,7],
-# ybottom = resid_probs[,3],
-# col = c_mid,
-# border = 'transparent')
-# rect(xleft = x[seq(1, N*2, by = 2)],
-# xright = x[seq(2, N*2, by = 2)],
-# ytop = resid_probs[,6],
-# ybottom = resid_probs[,4],
-# col = c_mid_highlight,
-# border = 'transparent')
-#
-# for (k in 1:N) {
-# lines(x = c(x[seq(1, N*2, by = 2)][k],x[seq(2, N*2, by = 2)][k]),
-# y = c(resid_probs[k,5], resid_probs[k,5]),
-# col = c_dark, lwd = 2)
-# }
-
-series_residuals <- object$resids[[series]]
+# Extract expectation (fitted) values
preds <- hindcast(object, type = 'expected')$hindcasts[[series]]
+# Plot resids vs fitted
plot(1,
xlim = c(min(apply(preds, 2, function(x) quantile(x, 0.05, na.rm = TRUE))),
max(apply(preds, 2, function(x) quantile(x, 0.95, na.rm = TRUE)))),
bty = 'L',
- xlab = 'Fitted values',
- ylab = 'DS residuals',
+ xlab = '',
+ ylab = '',
pch = 16,
col = 'white',
cex = 1,
ylim = range(series_residuals, na.rm = T))
+
draws <- sample(1:dim(preds)[1],
- min(300, dim(preds)[1]),
+ min(200, dim(preds)[1]),
replace = FALSE)
for(i in draws){
points(x = preds[i,],
y = series_residuals[i,],
pch = 16,
- cex = 0.8,
- col = sample(c("#B97C7C60",
- "#A2505060",
- "#8F272760",
- "#7C000060"),
- 1))
+ cex = 0.7,
+ col = '#80808020')
}
title('Resids vs Fitted', line = 0)
-abline(h = 0, col = '#FFFFFF60', lwd = 2.85)
-abline(h = 0, col = 'black', lwd = 2.5, lty = 'dashed')
+title(ylab = "DS residuals", line = 1.5)
+title(xlab = "Fitted values", line = 1.5)
+
+# Plot quantile regression lines
+y <- apply(series_residuals, 2, function(x) quantile(x, 0.5, na.rm = TRUE))
+x <- apply(preds, 2, function(x) quantile(x, 0.5, na.rm = TRUE))
+predvals <- seq(min(preds),
+ max(preds),
+ length.out = 200)
+medmod <- gam(y ~ s(x))
+medpreds <- predict(medmod, newdata = data.frame(x = predvals),
+ type = 'response', se.fit = TRUE)
+
+polygon(c(predvals, rev(predvals)), c(medpreds$fit + 2*medpreds$se.fit,
+ rev(medpreds$fit - 2*medpreds$se.fit)),
+ col = "#7C000040",
+ border = NA)
+lines(x = predvals,
+ y = medpreds$fit,
+ col = "#7C000060",
+ lwd = 3)
# Q-Q plot
coords <- qqnorm(series_residuals[1,], plot.it = F)
@@ -375,14 +279,16 @@ pred_vals <- pred_vals[complete.cases(cred[1,])]
plot(x = pred_vals,
y = cred[5,][complete.cases(cred[1,])],
bty = 'L',
- xlab = 'Theoretical Quantiles',
- ylab = 'Sample Quantiles',
+ xlab = '',
+ ylab = '',
pch = 16,
col = 'white',
cex = 1,
ylim = range(cred, na.rm = T),
tck = -0.04)
title('Normal Q-Q Plot', line = 0)
+title(ylab = "Sample Quantiles", line = 1.5)
+title(xlab = "Theoretical Quantiles", line = 1.5)
polygon(c(pred_vals, rev(pred_vals)), c(cred[1,][complete.cases(cred[1,])],
rev(cred[9,][complete.cases(cred[1,])])),
col = c_light, border = NA)
@@ -424,8 +330,8 @@ cred <- sapply(1:NCOL(resid_acf),
cred <- cred[, -1]
clim <- qnorm((1 + .95)/2)/sqrt(acf1$n.used)
plot(1, type = "n", bty = 'L',
- xlab = 'Lag',
- ylab = 'Autocorrelation',
+ xlab = '',
+ ylab = '',
xlim = c(1, N-1),
xaxt = 'n',
ylim = range(c(cred,
@@ -433,6 +339,8 @@ plot(1, type = "n", bty = 'L',
clim + 0.05)))
axis(1, at = seq(1, NCOL(cred), by = 2))
title('ACF', line = 0)
+title(ylab = "Autocorrelation", line = 1.5)
+title(xlab = "Lag", line = 1.5)
N <- N - 1
rect(xleft = x[seq(1, N*2, by = 2)],
@@ -495,8 +403,8 @@ cred <- sapply(1:NCOL(resid_pacf),
clim <- qnorm((1 + .95)/2)/sqrt(pacf1$n.used)
plot(1, type = "n", bty = 'L',
- xlab = 'Lag',
- ylab = 'Autocorrelation',
+ xlab = '',
+ ylab = '',
xlim = c(1, length(sorted_x)),
xaxt = 'n',
ylim = range(c(cred,
@@ -504,6 +412,8 @@ plot(1, type = "n", bty = 'L',
clim + 0.05)))
axis(1, at = seq(1, NCOL(cred), by = 2))
title('pACF', line = 0)
+title(ylab = "Autocorrelation", line = 1.5)
+title(xlab = "Lag", line = 1.5)
rect(xleft = x[seq(1, N*2, by = 2)],
xright = x[seq(2, N*2, by = 2)],
@@ -539,6 +449,9 @@ abline(h = clim, col = '#FFFFFF60', lwd = 2.85)
abline(h = clim, col = 'black', lwd = 2.5, lty = 'dashed')
abline(h = -clim, col = '#FFFFFF60', lwd = 2.85)
abline(h = -clim, col = 'black', lwd = 2.5, lty = 'dashed')
+
+invisible()
+par(.pardefault)
layout(1)
}
diff --git a/R/plot_mvgam_series.R b/R/plot_mvgam_series.R
index 863bf873..bcd71ff7 100644
--- a/R/plot_mvgam_series.R
+++ b/R/plot_mvgam_series.R
@@ -156,8 +156,6 @@ plot_mvgam_series = function(object,
if(series == 'all'){
n_plots <- length(levels(data_train$series))
pages <- 1
- .pardefault <- par(no.readonly=T)
- par(.pardefault)
if (n_plots > 4) pages <- 2
if (n_plots > 8) pages <- 3
@@ -175,7 +173,13 @@ plot_mvgam_series = function(object,
if (c<1) r <- c <- 1
if (c*r < ppp) c <- c + 1
if (c*r < ppp) r <- r + 1
- oldpar<-par(mfrow=c(r,c))
+
+ .pardefault <- par(no.readonly=T)
+ par(.pardefault)
+ oldpar<-par(mfrow=c(r,c),
+ mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(1.5, 0.5, 0))
} else { ppp <- 1; oldpar <- par()}
@@ -266,6 +270,12 @@ plot_mvgam_series = function(object,
}
layout(matrix(1:4, nrow = 2, byrow = TRUE))
+ .pardefault <- par(no.readonly=T)
+ par(.pardefault)
+ oldpar <- par(mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(2, 0.5, 0))
+
if(!missing(data_test)){
test <- data.frame(y = data_test$y,
series = data_test$series,
@@ -277,11 +287,13 @@ plot_mvgam_series = function(object,
dplyr::pull(y)
plot(1, type = "n", bty = 'L',
- xlab = 'Time',
- ylab = ylab,
+ xlab = '',
+ ylab = '',
ylim = range(c(truth, test), na.rm = TRUE),
xlim = c(0, length(c(truth, test))))
title('Time series', line = 0)
+ title(ylab = ylab, line = 1.5)
+ title(xlab = "Time", line = 1.5)
if(lines){
lines(x = 1:length(truth), y = truth, lwd = 2, col = "#8F2727")
@@ -307,15 +319,17 @@ plot_mvgam_series = function(object,
freq = FALSE,
breaks = n_bins,
col = "#C79999",
- ylab = 'Density',
- xlab = paste0(y),
+ ylab = '',
+ xlab = '',
main = '')
title('Histogram', line = 0)
+ title(ylab = 'Density', line = 1.5)
+ title(xlab = paste0(y), line = 1.5)
acf(c(truth, test),
na.action = na.pass, bty = 'L',
lwd = 2.5, ci.col = 'black', col = "#8F2727",
- main = '', ylab = 'Autocorrelation')
+ main = '', ylab = '', xlab = '')
acf1 <- acf(c(truth, test), plot = F,
na.action = na.pass)
clim <- qnorm((1 + .95)/2)/sqrt(acf1$n.used)
@@ -325,6 +339,8 @@ plot_mvgam_series = function(object,
abline(h = -clim, col = 'black', lwd = 2.5, lty = 'dashed')
box(bty = 'L', lwd = 2)
title('ACF', line = 0)
+ title(ylab = 'Autocorrelation', line = 1.5)
+ title(xlab = 'Lag', line = 1.5)
ecdf_plotdat = function(vals, x){
func <- ecdf(vals)
@@ -336,11 +352,13 @@ plot_mvgam_series = function(object,
length.out = 100)
plot(1, type = "n", bty = 'L',
- xlab = paste0(y),
- ylab = 'Empirical CDF',
+ xlab = '',
+ ylab = '',
xlim = c(min(plot_x), max(plot_x)),
ylim = c(0, 1))
title('CDF', line = 0)
+ title(ylab = 'Empirical CDF', line = 1.5)
+ title(xlab = paste0(y), line = 1.5)
lines(x = plot_x,
y = ecdf_plotdat(c(truth, test),
plot_x),
@@ -350,11 +368,13 @@ plot_mvgam_series = function(object,
} else {
plot(1, type = "n", bty = 'L',
- xlab = 'Time',
- ylab = ylab,
+ xlab = '',
+ ylab = '',
ylim = range(c(truth), na.rm = TRUE),
xlim = c(0, length(c(truth))))
title('Time series', line = 0)
+ title(ylab = ylab, line = 1.5)
+ title(xlab = "Time", line = 1.5)
if(lines){
lines(x = 1:length(truth), y = truth, lwd = 2, col = "#8F2727")
@@ -374,16 +394,18 @@ plot_mvgam_series = function(object,
freq = FALSE,
breaks = n_bins,
col = "#C79999",
- ylab = 'Density',
- xlab = paste0(y),
+ ylab = '',
+ xlab = '',
main = '')
title('Histogram', line = 0)
+ title(ylab = 'Density', line = 1.5)
+ title(xlab = paste0(y), line = 1.5)
acf(c(truth),
na.action = na.pass, bty = 'L',
lwd = 2.5, ci.col = 'black', col = "#8F2727",
- main = '', ylab = 'Autocorrelation')
+ main = '', ylab = '', xlab = '')
acf1 <- acf(c(truth), plot = F,
na.action = na.pass)
clim <- qnorm((1 + .95)/2)/sqrt(acf1$n.used)
@@ -393,6 +415,8 @@ plot_mvgam_series = function(object,
abline(h = -clim, col = 'black', lwd = 2.5, lty = 'dashed')
box(bty = 'L', lwd = 2)
title('ACF', line = 0)
+ title(ylab = 'Autocorrelation', line = 1.5)
+ title(xlab = 'Lag', line = 1.5)
ecdf_plotdat = function(vals, x){
@@ -404,11 +428,13 @@ plot_mvgam_series = function(object,
to = max(truth, na.rm = T),
length.out = 100)
plot(1, type = "n", bty = 'L',
- xlab = paste0(y),
- ylab = 'Empirical CDF',
+ xlab = '',
+ ylab = '',
xlim = c(min(plot_x), max(plot_x)),
ylim = c(0, 1))
title('CDF', line = 0)
+ title(ylab = 'Empirical CDF', line = 1.5)
+ title(xlab = paste0(y), line = 1.5)
lines(x = plot_x,
y = ecdf_plotdat(truth,
plot_x),
@@ -418,5 +444,7 @@ plot_mvgam_series = function(object,
}
}
+ invisible()
+ par(.pardefault)
layout(1)
}
diff --git a/R/plot_mvgam_smooth.R b/R/plot_mvgam_smooth.R
index eecc1294..a3a5e0a5 100644
--- a/R/plot_mvgam_smooth.R
+++ b/R/plot_mvgam_smooth.R
@@ -361,8 +361,9 @@ plot_mvgam_smooth = function(object,
.pardefault <- par(no.readonly=T)
par(.pardefault)
par(mfrow = c(2, 1),
- mgp = c(2.5, 1, 0),
- mai = c(0.8, 0.8, 0.4, 0.4))
+ mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(1.5, 0.5, 0))
if(residuals){
plot(1, type = "n", bty = 'L',
diff --git a/R/plot_mvgam_trend.R b/R/plot_mvgam_trend.R
index f878fb0a..80657a2d 100644
--- a/R/plot_mvgam_trend.R
+++ b/R/plot_mvgam_trend.R
@@ -162,8 +162,9 @@ plot_mvgam_trend = function(object, series = 1, newdata, data_test,
.pardefault <- par(no.readonly=T)
par(.pardefault)
par(mfrow = c(2, 1),
- mgp = c(2.5, 1, 0),
- mai = c(0.8, 0.8, 0.4, 0.4))
+ mar=c(2.5, 2.3, 2, 2),
+ oma = c(1, 1, 0, 0),
+ mgp = c(1.5, 0.5, 0))
plot(1, type = "n", bty = 'L',
xlab = xlab,
diff --git a/README.md b/README.md
index 0d3f978a..37722f63 100644
--- a/README.md
+++ b/README.md
@@ -329,11 +329,11 @@ test_priors
#> 4 trend AR3 coefficient ar3 ~ std_normal();
#> 5 trend sd sigma ~ exponential(2);
#> example_change new_lowerbound new_upperbound
-#> 1 lambda ~ exponential(0.86); NA NA
-#> 2 ar1 ~ normal(0.36, 0.67); NA NA
-#> 3 ar2 ~ normal(0.81, 0.14); NA NA
-#> 4 ar3 ~ normal(-0.7, 0.91); NA NA
-#> 5 sigma ~ exponential(0.99); NA NA
+#> 1 lambda ~ exponential(0.06); NA NA
+#> 2 ar1 ~ normal(-0.74, 0.56); NA NA
+#> 3 ar2 ~ normal(-0.31, 0.74); NA NA
+#> 4 ar3 ~ normal(0.16, 0.52); NA NA
+#> 5 sigma ~ exponential(0.27); NA NA
```
Any of the above priors can be changed by modifying the `prior` column
@@ -361,17 +361,17 @@ lynx_mvgam <- mvgam(data = lynx_train,
#> Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup)
-#> Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup)
+#> Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup)
+#> Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 3 Iteration: 100 / 1000 [ 10%] (Warmup)
-#> Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 4 Iteration: 200 / 1000 [ 20%] (Warmup)
+#> Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup)
-#> Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 4 Iteration: 300 / 1000 [ 30%] (Warmup)
-#> Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup)
+#> Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 4 Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup)
@@ -379,40 +379,40 @@ lynx_mvgam <- mvgam(data = lynx_train,
#> Chain 3 Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 4 Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
-#> Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup)
-#> Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
+#> Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup)
+#> Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 3 Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 4 Iteration: 600 / 1000 [ 60%] (Sampling)
-#> Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling)
+#> Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 3 Iteration: 600 / 1000 [ 60%] (Sampling)
+#> Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 4 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 3 Iteration: 700 / 1000 [ 70%] (Sampling)
-#> Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling)
+#> Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 4 Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 3 Iteration: 800 / 1000 [ 80%] (Sampling)
-#> Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling)
-#> Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling)
+#> Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 4 Iteration: 900 / 1000 [ 90%] (Sampling)
+#> Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 3 Iteration: 900 / 1000 [ 90%] (Sampling)
-#> Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling)
-#> Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 1 finished in 26.0 seconds.
+#> Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
+#> Chain 2 finished in 23.8 seconds.
#> Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 4 finished in 26.5 seconds.
+#> Chain 4 finished in 25.0 seconds.
+#> Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
+#> Chain 1 finished in 25.2 seconds.
#> Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 3 finished in 27.2 seconds.
-#> Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 2 finished in 28.3 seconds.
+#> Chain 3 finished in 26.1 seconds.
#>
#> All 4 chains finished successfully.
-#> Mean chain execution time: 27.0 seconds.
-#> Total execution time: 28.4 seconds.
+#> Mean chain execution time: 25.0 seconds.
+#> Total execution time: 26.2 seconds.
```
Inspect the resulting model file, which is written in the `Stan`
@@ -601,29 +601,29 @@ summary(lynx_mvgam)
#> Fitted using Stan
#>
#> GAM coefficient (beta) estimates:
-#> 2.5% 50% 97.5% Rhat n.eff
-#> (Intercept) 6.77553825 6.798625000 6.82159025 1.00 2510
-#> s(season).1 -0.56987252 0.058678500 0.69993310 1.00 800
-#> s(season).2 -0.22794177 0.802511500 1.83384225 1.00 320
-#> s(season).3 -0.02880991 1.187900000 2.45974750 1.00 274
-#> s(season).4 -0.52695260 0.435967000 1.31273425 1.00 728
-#> s(season).5 -1.15109150 -0.114330500 0.98096142 1.01 369
-#> s(season).6 -1.01111050 0.004989565 1.09583825 1.01 425
-#> s(season).7 -0.67298565 0.368590500 1.36816900 1.00 562
-#> s(season).8 -0.99709282 0.185209500 1.79495750 1.00 257
-#> s(season).9 -1.20259625 -0.372599000 0.59993695 1.01 362
-#> s(season).10 -1.37199850 -0.691203500 -0.03623181 1.00 436
+#> 2.5% 50% 97.5% Rhat n.eff
+#> (Intercept) 6.774079500 6.79910500 6.82217125 1.00 2796
+#> s(season).1 -0.563900900 0.07171495 0.73274307 1.00 1073
+#> s(season).2 -0.215040125 0.75739600 1.78958300 1.00 379
+#> s(season).3 0.009095727 1.15399500 2.38966000 1.00 362
+#> s(season).4 -0.494926525 0.42365300 1.41556700 1.00 924
+#> s(season).5 -1.157187000 -0.09625150 0.96286345 1.00 446
+#> s(season).6 -1.031363750 0.07960010 1.14492800 1.00 672
+#> s(season).7 -0.693688575 0.38010600 1.44162925 1.00 598
+#> s(season).8 -1.033852250 0.16963850 1.70781550 1.00 309
+#> s(season).9 -1.182917750 -0.37433300 0.54939055 1.01 486
+#> s(season).10 -1.371112000 -0.69670700 -0.03767622 1.00 555
#>
#> GAM observation smoothing parameter (rho) estimates:
-#> 2.5% 50% 97.5% Rhat n.eff
-#> s(season) 2.082392 3.392625 4.264951 1 369
+#> 2.5% 50% 97.5% Rhat n.eff
+#> s(season) 2.164928 3.43265 4.224055 1 439
#>
#> Latent trend AR parameter estimates:
-#> 2.5% 50% 97.5% Rhat n.eff
-#> ar1[1] 0.7309469 1.1216950 1.42165350 1 527
-#> ar2[1] -0.8288280 -0.4066000 0.05082151 1 1431
-#> ar3[1] -0.4780412 -0.1362670 0.31519042 1 438
-#> sigma[1] 0.3964197 0.4939105 0.62026590 1 967
+#> 2.5% 50% 97.5% Rhat n.eff
+#> ar1[1] 0.7294796 1.120010 1.41617775 1.01 612
+#> ar2[1] -0.8173138 -0.392940 0.07362988 1.00 1341
+#> ar3[1] -0.4806702 -0.148709 0.25175495 1.00 539
+#> sigma[1] 0.4005763 0.494210 0.61646728 1.00 1146
#>
#> Stan MCMC diagnostics:
#> n_eff / iter looks reasonable for all parameters
@@ -731,7 +731,7 @@ plot(lynx_mvgam, type = 'forecast', newdata = lynx_test)
#> Out of sample DRPS:
- #> [1] 969.5303
+ #> [1] 1003.652
And the estimated latent trend component, again using the more flexible
`plot_mvgam_...()` option to show first derivatives of the estimated
@@ -828,57 +828,57 @@ lynx_mvgam_poor <- mvgam(data = lynx_train,
#> Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup)
+#> Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup)
+#> Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
-#> Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup)
-#> Chain 3 Iteration: 100 / 1000 [ 10%] (Warmup)
-#> Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup)
-#> Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
-#> Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
-#> Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup)
-#> Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup)
-#> Chain 3 Iteration: 400 / 1000 [ 40%] (Warmup)
-#> Chain 3 Iteration: 500 / 1000 [ 50%] (Warmup)
-#> Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
-#> Chain 3 Iteration: 600 / 1000 [ 60%] (Sampling)
+#> Chain 3 Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 4 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 4 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 4 Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 4 Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
-#> Chain 4 Iteration: 600 / 1000 [ 60%] (Sampling)
-#> Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling)
-#> Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling)
+#> Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
+#> Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling)
-#> Chain 3 Iteration: 700 / 1000 [ 70%] (Sampling)
-#> Chain 3 Iteration: 800 / 1000 [ 80%] (Sampling)
+#> Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup)
+#> Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup)
+#> Chain 3 Iteration: 400 / 1000 [ 40%] (Warmup)
+#> Chain 3 Iteration: 500 / 1000 [ 50%] (Warmup)
+#> Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
+#> Chain 4 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 4 Iteration: 700 / 1000 [ 70%] (Sampling)
+#> Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling)
+#> Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling)
+#> Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling)
+#> Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
+#> Chain 3 Iteration: 600 / 1000 [ 60%] (Sampling)
+#> Chain 3 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 4 Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 4 Iteration: 900 / 1000 [ 90%] (Sampling)
+#> Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
+#> Chain 2 finished in 1.1 seconds.
+#> Chain 4 finished in 1.1 seconds.
#> Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 1 finished in 1.1 seconds.
-#> Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling)
-#> Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
+#> Chain 3 Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 3 Iteration: 900 / 1000 [ 90%] (Sampling)
+#> Chain 1 finished in 1.2 seconds.
#> Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 2 finished in 1.2 seconds.
-#> Chain 3 finished in 1.2 seconds.
-#> Chain 4 finished in 1.2 seconds.
+#> Chain 3 finished in 1.3 seconds.
#>
#> All 4 chains finished successfully.
#> Mean chain execution time: 1.2 seconds.
-#> Total execution time: 1.3 seconds.
+#> Total execution time: 1.4 seconds.
```
The first approximator targets each model’s ability to simulate temporal
@@ -895,9 +895,9 @@ each competing model at the same set of timepoints.
``` r
compare_mvgams(lynx_mvgam, lynx_mvgam_poor, fc_horizon = 10)
#> RPS summaries per model (lower is better)
-#> Min. 1st Qu. Median Mean 3rd Qu. Max.
-#> Model 1 715.6817 983.1872 999.9234 1058.672 1176.271 1400.928
-#> Model 2 1346.9846 1622.1888 2219.2576 2391.485 3177.853 3730.506
+#> Min. 1st Qu. Median Mean 3rd Qu. Max.
+#> Model 1 709.4234 986.6266 1021.076 1069.682 1202.013 1427.319
+#> Model 2 1362.2546 1617.4447 2207.471 2388.815 3176.534 3731.437
#>
#> 90% interval coverages per model (closer to 0.9 is better)
#> Model 1 0.95
@@ -1082,53 +1082,53 @@ mod <- mvgam(y ~ s(season, bs = 'cc', k = 7) +
#> Chain 3 Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup)
-#> Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup)
-#> Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup)
+#> Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 4 Iteration: 200 / 1000 [ 20%] (Warmup)
-#> Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup)
+#> Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup)
-#> Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup)
+#> Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 4 Iteration: 300 / 1000 [ 30%] (Warmup)
+#> Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 3 Iteration: 400 / 1000 [ 40%] (Warmup)
#> Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup)
-#> Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 4 Iteration: 400 / 1000 [ 40%] (Warmup)
+#> Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 3 Iteration: 500 / 1000 [ 50%] (Warmup)
-#> Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup)
+#> Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 4 Iteration: 500 / 1000 [ 50%] (Warmup)
-#> Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
+#> Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 3 Iteration: 600 / 1000 [ 60%] (Sampling)
-#> Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 4 Iteration: 600 / 1000 [ 60%] (Sampling)
+#> Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 3 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 4 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 3 Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling)
-#> Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 4 Iteration: 800 / 1000 [ 80%] (Sampling)
+#> Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 3 Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 4 Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 2 finished in 9.5 seconds.
-#> Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 3 finished in 9.7 seconds.
#> Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling)
+#> Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
+#> Chain 2 finished in 10.0 seconds.
+#> Chain 3 finished in 10.0 seconds.
#> Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 4 finished in 10.1 seconds.
+#> Chain 4 finished in 10.4 seconds.
#> Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 1 finished in 12.0 seconds.
+#> Chain 1 finished in 12.3 seconds.
#>
#> All 4 chains finished successfully.
-#> Mean chain execution time: 10.3 seconds.
-#> Total execution time: 12.1 seconds.
+#> Mean chain execution time: 10.7 seconds.
+#> Total execution time: 12.4 seconds.
```
Inspect the summary to see that the posterior now also contains
@@ -1280,10 +1280,11 @@ mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE),
#> Running MCMC with 4 parallel chains...
#>
#> Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup)
+#> Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup)
-#> Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup)
+#> Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup)
@@ -1294,7 +1295,6 @@ mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE),
#> Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 3 Iteration: 400 / 1000 [ 40%] (Warmup)
-#> Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup)
#> Chain 4 Iteration: 200 / 1000 [ 20%] (Warmup)
#> Chain 4 Iteration: 300 / 1000 [ 30%] (Warmup)
#> Chain 4 Iteration: 400 / 1000 [ 40%] (Warmup)
@@ -1305,17 +1305,17 @@ mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE),
#> Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 4 Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
+#> Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup)
#> Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
-#> Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
-#> Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 3 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 4 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling)
+#> Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling)
#> Chain 3 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 4 Iteration: 700 / 1000 [ 70%] (Sampling)
-#> Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling)
+#> Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling)
#> Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 3 Iteration: 800 / 1000 [ 80%] (Sampling)
#> Chain 4 Iteration: 800 / 1000 [ 80%] (Sampling)
@@ -1326,15 +1326,15 @@ mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE),
#> Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling)
#> Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 1 finished in 1.3 seconds.
-#> Chain 3 finished in 1.3 seconds.
-#> Chain 4 finished in 1.3 seconds.
+#> Chain 1 finished in 1.2 seconds.
+#> Chain 3 finished in 1.2 seconds.
+#> Chain 4 finished in 1.2 seconds.
#> Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
-#> Chain 2 finished in 1.4 seconds.
+#> Chain 2 finished in 1.3 seconds.
#>
#> All 4 chains finished successfully.
-#> Mean chain execution time: 1.3 seconds.
-#> Total execution time: 1.5 seconds.
+#> Mean chain execution time: 1.2 seconds.
+#> Total execution time: 1.4 seconds.
```
Inspect the model summary, which now contains estimates for the
@@ -1343,7 +1343,7 @@ observation errors
``` r
summary(mod)
#> GAM formula:
-#> out ~ dynamic(temp, rho = 8, stationary = TRUE)
+#> out ~ s(time, by = temp, bs = "gp", m = c(-2, 8, 2), k = 27)
#>
#> Family:
#> gaussian
@@ -1444,7 +1444,7 @@ plot(mod, type = 'forecast', newdata = data_test)
#> Out of sample CRPS:
- #> [1] 1.768885
+ #> [1] 1.780129
There are many more extended uses for `mvgam` models, including the
ability to fit dynamic factor processes for analysing and forecasting
diff --git a/man/figures/README-unnamed-chunk-10-1.png b/man/figures/README-unnamed-chunk-10-1.png
index 3c028ccf..cce84b0b 100644
Binary files a/man/figures/README-unnamed-chunk-10-1.png and b/man/figures/README-unnamed-chunk-10-1.png differ
diff --git a/man/figures/README-unnamed-chunk-14-1.png b/man/figures/README-unnamed-chunk-14-1.png
index 94825220..9a4ba5f2 100644
Binary files a/man/figures/README-unnamed-chunk-14-1.png and b/man/figures/README-unnamed-chunk-14-1.png differ
diff --git a/man/figures/README-unnamed-chunk-15-1.png b/man/figures/README-unnamed-chunk-15-1.png
index 6397af79..80c8d07a 100644
Binary files a/man/figures/README-unnamed-chunk-15-1.png and b/man/figures/README-unnamed-chunk-15-1.png differ
diff --git a/man/figures/README-unnamed-chunk-19-1.png b/man/figures/README-unnamed-chunk-19-1.png
index 082a947c..733b45a7 100644
Binary files a/man/figures/README-unnamed-chunk-19-1.png and b/man/figures/README-unnamed-chunk-19-1.png differ
diff --git a/man/figures/README-unnamed-chunk-20-1.png b/man/figures/README-unnamed-chunk-20-1.png
index d905f579..482ca7c6 100644
Binary files a/man/figures/README-unnamed-chunk-20-1.png and b/man/figures/README-unnamed-chunk-20-1.png differ
diff --git a/man/figures/README-unnamed-chunk-21-1.png b/man/figures/README-unnamed-chunk-21-1.png
index 3016912a..8893932b 100644
Binary files a/man/figures/README-unnamed-chunk-21-1.png and b/man/figures/README-unnamed-chunk-21-1.png differ
diff --git a/man/figures/README-unnamed-chunk-22-1.png b/man/figures/README-unnamed-chunk-22-1.png
index 1f9c7e0f..51c3ec2e 100644
Binary files a/man/figures/README-unnamed-chunk-22-1.png and b/man/figures/README-unnamed-chunk-22-1.png differ
diff --git a/man/figures/README-unnamed-chunk-23-1.png b/man/figures/README-unnamed-chunk-23-1.png
index dd082876..dbb9347a 100644
Binary files a/man/figures/README-unnamed-chunk-23-1.png and b/man/figures/README-unnamed-chunk-23-1.png differ
diff --git a/man/figures/README-unnamed-chunk-25-1.png b/man/figures/README-unnamed-chunk-25-1.png
index 43522721..2826828e 100644
Binary files a/man/figures/README-unnamed-chunk-25-1.png and b/man/figures/README-unnamed-chunk-25-1.png differ
diff --git a/man/figures/README-unnamed-chunk-26-1.png b/man/figures/README-unnamed-chunk-26-1.png
index f95496ff..da925be4 100644
Binary files a/man/figures/README-unnamed-chunk-26-1.png and b/man/figures/README-unnamed-chunk-26-1.png differ
diff --git a/man/figures/README-unnamed-chunk-27-1.png b/man/figures/README-unnamed-chunk-27-1.png
index 835c4423..04b3004c 100644
Binary files a/man/figures/README-unnamed-chunk-27-1.png and b/man/figures/README-unnamed-chunk-27-1.png differ
diff --git a/man/figures/README-unnamed-chunk-28-1.png b/man/figures/README-unnamed-chunk-28-1.png
index 42172215..c5cff90c 100644
Binary files a/man/figures/README-unnamed-chunk-28-1.png and b/man/figures/README-unnamed-chunk-28-1.png differ
diff --git a/man/figures/README-unnamed-chunk-28-2.png b/man/figures/README-unnamed-chunk-28-2.png
index 4c16d4d6..5b086767 100644
Binary files a/man/figures/README-unnamed-chunk-28-2.png and b/man/figures/README-unnamed-chunk-28-2.png differ
diff --git a/man/figures/README-unnamed-chunk-29-1.png b/man/figures/README-unnamed-chunk-29-1.png
index 9f6edea8..43ab80d4 100644
Binary files a/man/figures/README-unnamed-chunk-29-1.png and b/man/figures/README-unnamed-chunk-29-1.png differ
diff --git a/man/figures/README-unnamed-chunk-30-1.png b/man/figures/README-unnamed-chunk-30-1.png
index 1d63f0c6..b800f966 100644
Binary files a/man/figures/README-unnamed-chunk-30-1.png and b/man/figures/README-unnamed-chunk-30-1.png differ
diff --git a/man/figures/README-unnamed-chunk-31-1.png b/man/figures/README-unnamed-chunk-31-1.png
index 5ee5483e..654f49c9 100644
Binary files a/man/figures/README-unnamed-chunk-31-1.png and b/man/figures/README-unnamed-chunk-31-1.png differ
diff --git a/man/figures/README-unnamed-chunk-32-1.png b/man/figures/README-unnamed-chunk-32-1.png
index f0433047..5222b2fb 100644
Binary files a/man/figures/README-unnamed-chunk-32-1.png and b/man/figures/README-unnamed-chunk-32-1.png differ
diff --git a/man/figures/README-unnamed-chunk-33-1.png b/man/figures/README-unnamed-chunk-33-1.png
index b712cdf8..9ccd84a1 100644
Binary files a/man/figures/README-unnamed-chunk-33-1.png and b/man/figures/README-unnamed-chunk-33-1.png differ
diff --git a/man/figures/README-unnamed-chunk-34-1.png b/man/figures/README-unnamed-chunk-34-1.png
index 463b2d77..ffd708e8 100644
Binary files a/man/figures/README-unnamed-chunk-34-1.png and b/man/figures/README-unnamed-chunk-34-1.png differ
diff --git a/man/figures/README-unnamed-chunk-35-1.png b/man/figures/README-unnamed-chunk-35-1.png
index 86574d44..d3e7a519 100644
Binary files a/man/figures/README-unnamed-chunk-35-1.png and b/man/figures/README-unnamed-chunk-35-1.png differ
diff --git a/man/figures/README-unnamed-chunk-36-1.png b/man/figures/README-unnamed-chunk-36-1.png
index 0c921d92..a4c528ca 100644
Binary files a/man/figures/README-unnamed-chunk-36-1.png and b/man/figures/README-unnamed-chunk-36-1.png differ
diff --git a/man/figures/README-unnamed-chunk-37-1.png b/man/figures/README-unnamed-chunk-37-1.png
index 2034970c..370572c4 100644
Binary files a/man/figures/README-unnamed-chunk-37-1.png and b/man/figures/README-unnamed-chunk-37-1.png differ
diff --git a/man/figures/README-unnamed-chunk-38-1.png b/man/figures/README-unnamed-chunk-38-1.png
index 11103161..175d64ae 100644
Binary files a/man/figures/README-unnamed-chunk-38-1.png and b/man/figures/README-unnamed-chunk-38-1.png differ
diff --git a/man/figures/README-unnamed-chunk-40-1.png b/man/figures/README-unnamed-chunk-40-1.png
index bf2055fd..40afa387 100644
Binary files a/man/figures/README-unnamed-chunk-40-1.png and b/man/figures/README-unnamed-chunk-40-1.png differ
diff --git a/man/figures/README-unnamed-chunk-40-2.png b/man/figures/README-unnamed-chunk-40-2.png
index 7ccc20c1..fbff6ff7 100644
Binary files a/man/figures/README-unnamed-chunk-40-2.png and b/man/figures/README-unnamed-chunk-40-2.png differ
diff --git a/man/figures/README-unnamed-chunk-40-3.png b/man/figures/README-unnamed-chunk-40-3.png
index eaf065ff..2dc854e7 100644
Binary files a/man/figures/README-unnamed-chunk-40-3.png and b/man/figures/README-unnamed-chunk-40-3.png differ
diff --git a/man/figures/README-unnamed-chunk-55-1.png b/man/figures/README-unnamed-chunk-55-1.png
index 7bd8d5ea..c844fd5b 100644
Binary files a/man/figures/README-unnamed-chunk-55-1.png and b/man/figures/README-unnamed-chunk-55-1.png differ
diff --git a/man/figures/README-unnamed-chunk-9-1.png b/man/figures/README-unnamed-chunk-9-1.png
index fda55589..89a1e59b 100644
Binary files a/man/figures/README-unnamed-chunk-9-1.png and b/man/figures/README-unnamed-chunk-9-1.png differ