diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 00000000..62ad0084 Binary files /dev/null and b/.DS_Store differ diff --git a/NEON_manuscript/.DS_Store b/NEON_manuscript/.DS_Store new file mode 100644 index 00000000..eaf6690a Binary files /dev/null and b/NEON_manuscript/.DS_Store differ diff --git a/NEON_manuscript/Case studies/.DS_Store b/NEON_manuscript/Case studies/.DS_Store new file mode 100644 index 00000000..aa69cddc Binary files /dev/null and b/NEON_manuscript/Case studies/.DS_Store differ diff --git a/NEON_manuscript/Case studies/mvgam_case_study2_files/.DS_Store b/NEON_manuscript/Case studies/mvgam_case_study2_files/.DS_Store new file mode 100644 index 00000000..e7b32978 Binary files /dev/null and b/NEON_manuscript/Case studies/mvgam_case_study2_files/.DS_Store differ diff --git a/NEON_manuscript/Case studies/rsconnect/.DS_Store b/NEON_manuscript/Case studies/rsconnect/.DS_Store new file mode 100644 index 00000000..6ca3a5fd Binary files /dev/null and b/NEON_manuscript/Case studies/rsconnect/.DS_Store differ diff --git a/NEON_manuscript/Case studies/rsconnect/documents/.DS_Store b/NEON_manuscript/Case studies/rsconnect/documents/.DS_Store new file mode 100644 index 00000000..15248328 Binary files /dev/null and b/NEON_manuscript/Case studies/rsconnect/documents/.DS_Store differ diff --git a/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study1.Rmd/.DS_Store b/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study1.Rmd/.DS_Store new file mode 100644 index 00000000..b4aee0ce Binary files /dev/null and b/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study1.Rmd/.DS_Store differ diff --git a/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study1.Rmd/rpubs.com/.DS_Store b/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study1.Rmd/rpubs.com/.DS_Store new file mode 100644 index 00000000..9f1da32b Binary files /dev/null and b/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study1.Rmd/rpubs.com/.DS_Store differ diff --git a/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/.DS_Store b/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/.DS_Store new file mode 100644 index 00000000..b4aee0ce Binary files /dev/null and b/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/.DS_Store differ diff --git a/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/rpubs.com/.DS_Store b/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/rpubs.com/.DS_Store new file mode 100644 index 00000000..9f1da32b Binary files /dev/null and b/NEON_manuscript/Case studies/rsconnect/documents/mvgam_case_study2.Rmd/rpubs.com/.DS_Store differ diff --git a/NEON_manuscript/Figures/.DS_Store b/NEON_manuscript/Figures/.DS_Store new file mode 100644 index 00000000..f3a3eff1 Binary files /dev/null and b/NEON_manuscript/Figures/.DS_Store differ diff --git a/NEON_manuscript/Manuscript/.DS_Store b/NEON_manuscript/Manuscript/.DS_Store new file mode 100644 index 00000000..4854a53d Binary files /dev/null and b/NEON_manuscript/Manuscript/.DS_Store differ diff --git a/NEON_manuscript/Manuscript/Appendix_S5_IxodesJAGS_Hyp3.txt b/NEON_manuscript/Manuscript/Appendix_S5_IxodesJAGS_Hyp3.txt index a55a3db1..fbbff53d 100644 --- a/NEON_manuscript/Manuscript/Appendix_S5_IxodesJAGS_Hyp3.txt +++ b/NEON_manuscript/Manuscript/Appendix_S5_IxodesJAGS_Hyp3.txt @@ -1,157 +1,157 @@ -JAGS model code generated by package mvgam - -GAM formula: -y ~ s(siteID, bs = re) + s(cum_gdd, k = 5) + s(cum_gdd, siteID, k = 5, m = 1, bs = fs) + s(season, m = 2, bs = cc, k = 12) + s(season, series, m = 1, k = 4, bs = fs) - -Trend model: -RW - -Required data: -integer n; number of timepoints per series -integer n_series; number of series -matrix y; time-ordered observations of dimension n x n_series (missing values allowed) -matrix ytimes; time-ordered n x n_series matrix (which row in X belongs to each [time, series] observation?) -matrix X; mgcv GAM design matrix of dimension (n x n_series) x basis dimension -matrix S2; mgcv smooth penalty matrix S2 -matrix S4; mgcv smooth penalty matrix S4 -vector zero; prior basis coefficient locations vector of length ncol(X) -vector p_coefs; vector (length = 1) of prior Gaussian means for parametric effects -vector p_taus; vector (length = 1) of prior Gaussian precisions for parametric effects - - -#### Begin model #### -model { - -## GAM linear predictor -eta <- X %*% b - -## mean expectations -for (i in 1:n) { -for (s in 1:n_series) { -mu[i, s] <- exp(eta[ytimes[i, s]] + trend[i, s]) -} -} - -## latent factors evolve as time series with penalised precisions; -## the penalty terms force any un-needed factors to evolve as flat lines -for (j in 1:n_lv) { -LV[1, j] ~ dnorm(0, penalty[j]) -} - -for (i in 2:n) { -for (j in 1:n_lv){ -LV[i, j] ~ dnorm(LV[i - 1, j], penalty[j]) -} -} - -## shrinkage penalties for each factor's precision act to squeeze -## the entire factor toward a flat white noise process if supported by -## the data. The prior for individual factor penalties allows each factor to possibly -## have a relatively large penalty, which shrinks the prior for that factor's variance -## substantially. Penalties increase exponentially with the number of factors following -## Welty, Leah J., et al. Bayesian distributed lag models: estimating effects of particulate -## matter air pollution on daily mortality Biometrics 65.1 (2009): 282-291. -pi ~ dunif(0, n_lv) -X2 ~ dnorm(0, 1)T(0, ) - -# eta1 controls the baseline penalty -eta1 ~ dunif(-1, 1) - -# eta2 controls how quickly the penalties exponentially increase -eta2 ~ dunif(-1, 1) - -for (t in 1:n_lv) { -X1[t] ~ dnorm(0, 1)T(0, ) -l.dist[t] <- max(t, pi[]) -l.weight[t] <- exp(eta2[] * l.dist[t]) -l.var[t] <- exp(eta1[] * l.dist[t] / 2) * 1 -theta.prime[t] <- l.weight[t] * X1[t] + (1 - l.weight[t]) * X2[] -penalty[t] <- max(0.0001, theta.prime[t] * l.var[t]) -} - -## latent factor loadings: standard normal with identifiability constraints -## upper triangle of loading matrix set to zero -for (j in 1:(n_lv - 1)) { -for (j2 in (j + 1):n_lv) { -lv_coefs[j, j2] <- 0 -} -} - -## positive constraints on loading diagonals -for (j in 1:n_lv) { -lv_coefs[j, j] ~ dnorm(0, 1)T(0, 1); -} - -## lower diagonal free -for (j in 2:n_lv) { -for (j2 in 1:(j - 1)) { -lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1); -} -} - -## other elements also free -for (j in (n_lv + 1):n_series) { -for (j2 in 1:n_lv) { -lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1); -} -} - -## trend evolution depends on latent factors -for (i in 1:n) { -for (s in 1:n_series) { -trend[i, s] <- inprod(lv_coefs[s,], LV[i,]) -} -} - -## likelihood functions -for (i in 1:n) { -for (s in 1:n_series) { -y[i, s] ~ dnegbin(rate[i, s], r[s]) -rate[i, s] <- ifelse((r[s] / (r[s] + mu[i, s])) < min_eps, min_eps, -(r[s] / (r[s] + mu[i, s]))) -} -} - -## complexity penalising prior for the overdispersion parameter; -## where the likelihood reduces to a 'base' model (Poisson) unless -## the data support overdispersion -for (s in 1:n_series) { -r[s] <- 1 / r_raw[s] -r_raw[s] ~ dnorm(0, 0.1)T(0, ) -} - -## posterior predictions -for (i in 1:n) { -for (s in 1:n_series) { -ypred[i, s] ~ dnegbin(rate[i, s], r[s]) -} -} - -## GAM-specific priors -## parametric effect priors (regularised for identifiability) -for (i in 1:1) { b[i] ~ dnorm(p_coefs[i], p_taus[i]) } -## prior (non-centred) for s(siteID)... -for (i in 2:4) { -b_raw[i] ~ dnorm(0, 1) -b[i] <- mu_raw1 + b_raw[i] * sigma_raw1 -} -sigma_raw1 ~ dexp(1) -mu_raw1 ~ dnorm(0, 1) -## prior for s(cum_gdd)... -K2 <- S2[1:4,1:4] * lambda[2] + S2[1:4,5:8] * lambda[3] -b[5:8] ~ dmnorm(zero[5:8],K2) -## prior for s(cum_gdd,siteID)... -for (i in c(9:12,14:17,19:22)) { b[i] ~ dnorm(0, lambda[4]) } -for (i in c(13,18,23)) { b[i] ~ dnorm(0, lambda[5]) } -## prior for s(season)... -K4 <- S4[1:10,1:10] * lambda[6] -b[24:33] ~ dmnorm(zero[24:33],K4) -## prior for s(season,series)... -for (i in c(34:36,38:40,42:44,46:48,50:52,54:56,58:60,62:64)) { b[i] ~ dnorm(0, lambda[7]) } -for (i in c(37,41,45,49,53,57,61,65)) { b[i] ~ dnorm(0, lambda[8]) } -## smoothing parameter priors... -for (i in 1:8) { -lambda[i] ~ dexp(0.05)T(0.0001, ) -rho[i] <- log(lambda[i]) -} -} +JAGS model code generated by package mvgam + +GAM formula: +y ~ s(siteID, bs = re) + s(cum_gdd, k = 5) + s(cum_gdd, siteID, k = 5, m = 1, bs = fs) + s(season, m = 2, bs = cc, k = 12) + s(season, series, m = 1, k = 4, bs = fs) + +Trend model: +RW + +Required data: +integer n; number of timepoints per series +integer n_series; number of series +matrix y; time-ordered observations of dimension n x n_series (missing values allowed) +matrix ytimes; time-ordered n x n_series matrix (which row in X belongs to each [time, series] observation?) +matrix X; mgcv GAM design matrix of dimension (n x n_series) x basis dimension +matrix S2; mgcv smooth penalty matrix S2 +matrix S4; mgcv smooth penalty matrix S4 +vector zero; prior basis coefficient locations vector of length ncol(X) +vector p_coefs; vector (length = 1) of prior Gaussian means for parametric effects +vector p_taus; vector (length = 1) of prior Gaussian precisions for parametric effects + + +#### Begin model #### +model { + +## GAM linear predictor +eta <- X %*% b + +## mean expectations +for (i in 1:n) { +for (s in 1:n_series) { +mu[i, s] <- exp(eta[ytimes[i, s]] + trend[i, s]) +} +} + +## latent factors evolve as time series with penalised precisions; +## the penalty terms force any un-needed factors to evolve as flat lines +for (j in 1:n_lv) { +LV[1, j] ~ dnorm(0, penalty[j]) +} + +for (i in 2:n) { +for (j in 1:n_lv){ +LV[i, j] ~ dnorm(LV[i - 1, j], penalty[j]) +} +} + +## shrinkage penalties for each factor's precision act to squeeze +## the entire factor toward a flat white noise process if supported by +## the data. The prior for individual factor penalties allows each factor to possibly +## have a relatively large penalty, which shrinks the prior for that factor's variance +## substantially. Penalties increase exponentially with the number of factors following +## Welty, Leah J., et al. Bayesian distributed lag models: estimating effects of particulate +## matter air pollution on daily mortality Biometrics 65.1 (2009): 282-291. +pi ~ dunif(0, n_lv) +X2 ~ dnorm(0, 1)T(0, ) + +# eta1 controls the baseline penalty +eta1 ~ dunif(-1, 1) + +# eta2 controls how quickly the penalties exponentially increase +eta2 ~ dunif(-1, 1) + +for (t in 1:n_lv) { +X1[t] ~ dnorm(0, 1)T(0, ) +l.dist[t] <- max(t, pi[]) +l.weight[t] <- exp(eta2[] * l.dist[t]) +l.var[t] <- exp(eta1[] * l.dist[t] / 2) * 1 +theta.prime[t] <- l.weight[t] * X1[t] + (1 - l.weight[t]) * X2[] +penalty[t] <- max(0.0001, theta.prime[t] * l.var[t]) +} + +## latent factor loadings: standard normal with identifiability constraints +## upper triangle of loading matrix set to zero +for (j in 1:(n_lv - 1)) { +for (j2 in (j + 1):n_lv) { +lv_coefs[j, j2] <- 0 +} +} + +## positive constraints on loading diagonals +for (j in 1:n_lv) { +lv_coefs[j, j] ~ dnorm(0, 1)T(0, 1); +} + +## lower diagonal free +for (j in 2:n_lv) { +for (j2 in 1:(j - 1)) { +lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1); +} +} + +## other elements also free +for (j in (n_lv + 1):n_series) { +for (j2 in 1:n_lv) { +lv_coefs[j, j2] ~ dnorm(0, 1)T(-1, 1); +} +} + +## trend evolution depends on latent factors +for (i in 1:n) { +for (s in 1:n_series) { +trend[i, s] <- inprod(lv_coefs[s,], LV[i,]) +} +} + +## likelihood functions +for (i in 1:n) { +for (s in 1:n_series) { +y[i, s] ~ dnegbin(rate[i, s], r[s]) +rate[i, s] <- ifelse((r[s] / (r[s] + mu[i, s])) < min_eps, min_eps, +(r[s] / (r[s] + mu[i, s]))) +} +} + +## complexity penalising prior for the overdispersion parameter; +## where the likelihood reduces to a 'base' model (Poisson) unless +## the data support overdispersion +for (s in 1:n_series) { +r[s] <- 1 / r_raw[s] +r_raw[s] ~ dnorm(0, 0.1)T(0, ) +} + +## posterior predictions +for (i in 1:n) { +for (s in 1:n_series) { +ypred[i, s] ~ dnegbin(rate[i, s], r[s]) +} +} + +## GAM-specific priors +## parametric effect priors (regularised for identifiability) +for (i in 1:1) { b[i] ~ dnorm(p_coefs[i], p_taus[i]) } +## prior (non-centred) for s(siteID)... +for (i in 2:4) { +b_raw[i] ~ dnorm(0, 1) +b[i] <- mu_raw1 + b_raw[i] * sigma_raw1 +} +sigma_raw1 ~ dexp(1) +mu_raw1 ~ dnorm(0, 1) +## prior for s(cum_gdd)... +K2 <- S2[1:4,1:4] * lambda[2] + S2[1:4,5:8] * lambda[3] +b[5:8] ~ dmnorm(zero[5:8],K2) +## prior for s(cum_gdd,siteID)... +for (i in c(9:12,14:17,19:22)) { b[i] ~ dnorm(0, lambda[4]) } +for (i in c(13,18,23)) { b[i] ~ dnorm(0, lambda[5]) } +## prior for s(season)... +K4 <- S4[1:10,1:10] * lambda[6] +b[24:33] ~ dmnorm(zero[24:33],K4) +## prior for s(season,series)... +for (i in c(34:36,38:40,42:44,46:48,50:52,54:56,58:60,62:64)) { b[i] ~ dnorm(0, lambda[7]) } +for (i in c(37,41,45,49,53,57,61,65)) { b[i] ~ dnorm(0, lambda[8]) } +## smoothing parameter priors... +for (i in 1:8) { +lambda[i] ~ dexp(0.05)T(0.0001, ) +rho[i] <- log(lambda[i]) +} +} diff --git a/NEON_manuscript/next_todo.R b/NEON_manuscript/next_todo.R index 1bdc01fc..f9e89f1d 100644 --- a/NEON_manuscript/next_todo.R +++ b/NEON_manuscript/next_todo.R @@ -6,14 +6,263 @@ dat$true_corrs mod1 <- mvgam(formula = y ~ s(season, bs = 'cc') + s(series, bs = 're'), data_train = dat$data_train, - trend_model = 'AR3', + trend_model = 'AR1', family = 'poisson', use_lv = TRUE, n_lv = 2, use_stan = TRUE, run_model = T, burnin = 10) +mod1$model_file summary(mod1) +plot_mvgam_factors(mod1) +plot(mod1, type = 'residuals') +fake <- dat$data_test +fake$y <- NULL +plot(mod1, 'trend', data_test = fake) + + + +plot_mvgam_series = function(data_train, data_test, series, n_bins, + log_scale = FALSE){ + + 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 + if (n_plots > 12) pages <- 4 + if (pages != 0) { + ppp <- n_plots %/% pages + + if (n_plots %% pages != 0) { + ppp <- ppp + 1 + while (ppp * (pages - 1) >= n_plots) pages <- pages - 1 + } + + # Configure layout matrix + c <- r <- trunc(sqrt(ppp)) + 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)) + + } else { ppp <- 1; oldpar <- par()} + + all_ys <- lapply(seq_len(n_plots), function(x){ + if(log_scale){ + log( data.frame(y = data_train$y, + series = data_train$series, + time = data_train$time) %>% + dplyr::filter(series == levels(data_train$series)[x]) %>% + dplyr::pull(y) + 1) + } else { + data.frame(y = data_train$y, + series = data_train$series, + time = data_train$time) %>% + dplyr::filter(series == levels(data_train$series)[x]) %>% + dplyr::pull(y) + } + }) + + for(i in 1:n_plots){ + s_name <- levels(data_train$series)[i] + + truth <- data.frame(y = data_train$y, + series = data_train$series, + time = data_train$time) %>% + dplyr::filter(series == s_name) %>% + dplyr::select(time, y) %>% + dplyr::distinct() %>% + dplyr::arrange(time) %>% + dplyr::pull(y) + + if(log_scale){ + truth <- log(truth + 1) + ylab <- 'log(Y + 1)' + } else { + ylab <- 'Y' + } + + plot(1, type = "n", bty = 'L', + xlab = 'Time', + ylab = ylab, + ylim = range(unlist(all_ys)), + xlim = c(0, length(c(truth)))) + title(s_name, line = 0) + + if(n_plots > 1){ + for(x in 1:n_plots){ + lines(all_ys[[x]], lwd = 1.85, col = 'grey85') + } + } + + lines(x = 1:length(truth), y = truth, lwd = 3, col = "white") + lines(x = 1:length(truth), y = truth, lwd = 2.5, col = "#8F2727") + box(bty = 'L', lwd = 2) + + } + + } else { + s_name <- levels(data_train$series)[series] + truth <- data.frame(y = data_train$y, + series = data_train$series, + time = data_train$time) %>% + dplyr::filter(series == s_name) %>% + dplyr::select(time, y) %>% + dplyr::distinct() %>% + dplyr::arrange(time) %>% + dplyr::pull(y) + + layout(matrix(1:4, nrow = 2, byrow = TRUE)) + if(!missing(data_test)){ + test <- data.frame(y = data_test$y, + series = data_test$series, + time = data_test$time) %>% + dplyr::filter(series == s_name) %>% + dplyr::select(time, y) %>% + dplyr::distinct() %>% + dplyr::arrange(time) %>% + dplyr::pull(y) + + plot(1, type = "n", bty = 'L', + xlab = 'Time', + ylab = 'Y', + ylim = range(c(truth, test)), + xlim = c(0, length(c(truth, test)))) + title('Time series', line = 0) + + lines(x = 1:length(truth), y = truth, lwd = 2, col = "#8F2727") + lines(x = (length(truth)+1):length(c(truth, test)), y = test, lwd = 2, col = "black") + abline(v = length(truth)+1, col = '#FFFFFF60', lwd = 2.85) + abline(v = length(truth)+1, col = 'black', lwd = 2.5, lty = 'dashed') + box(bty = 'L', lwd = 2) + + if(missing(n_bins)){ + n_bins <- max(c(length(hist(c(truth, test), plot = F)$breaks), + 20)) + } + + hist(c(truth, test), border = "#8F2727", + lwd = 2, + freq = FALSE, + breaks = n_bins, + col = "#C79999", + ylab = 'Density', + xlab = 'Count', main = '') + title('Histogram', line = 0) + + acf(c(truth, test), + na.action = na.pass, bty = 'L', + lwd = 2.5, ci.col = 'black', col = "#8F2727", + main = '', ylab = 'Autocorrelation') + acf1 <- acf(c(truth, test), plot = F, + na.action = na.pass) + clim <- qnorm((1 + .95)/2)/sqrt(acf1$n.used) + abline(h = clim, col = '#FFFFFF', lwd = 2.85) + abline(h = clim, col = 'black', lwd = 2.5, lty = 'dashed') + abline(h = -clim, col = '#FFFFFF', lwd = 2.85) + abline(h = -clim, col = 'black', lwd = 2.5, lty = 'dashed') + box(bty = 'L', lwd = 2) + title('ACF', line = 0) + + ecdf_plotdat = function(vals, x){ + func <- ecdf(vals) + func(x) + } + + plot_x <- seq(min(c(truth, test), na.rm = T), + max(c(truth, test), na.rm = T)) + plot(1, type = "n", bty = 'L', + xlab = 'Count', + ylab = 'Empirical CDF', + xlim = c(min(plot_x), max(plot_x)), + ylim = c(0, 1)) + title('CDF', line = 0) + lines(x = plot_x, + y = ecdf_plotdat(c(truth, test), + plot_x), + col = "#8F2727", + lwd = 2.5) + box(bty = 'L', lwd = 2) + + } else { + plot(1, type = "n", bty = 'L', + xlab = 'Time', + ylab = 'Observations', + ylim = range(c(truth)), + xlim = c(0, length(c(truth)))) + title('Time series', line = 0) + + lines(x = 1:length(truth), y = truth, lwd = 2, col = "#8F2727") + box(bty = 'L', lwd = 2) + + if(missing(n_bins)){ + n_bins <- max(c(length(hist(c(truth), plot = F)$breaks), + 20)) + } + + hist(c(truth), border = "#8F2727", + lwd = 2, + freq = FALSE, + breaks = n_bins, + col = "#C79999", + ylab = 'Density', + xlab = 'Count', main = '') + title('Histogram', line = 0) + + + acf(c(truth), + na.action = na.pass, bty = 'L', + lwd = 2.5, ci.col = 'black', col = "#8F2727", + main = '', ylab = 'Autocorrelation') + acf1 <- acf(c(truth), plot = F, + na.action = na.pass) + clim <- qnorm((1 + .95)/2)/sqrt(acf1$n.used) + abline(h = clim, col = '#FFFFFF', lwd = 2.85) + abline(h = clim, col = 'black', lwd = 2.5, lty = 'dashed') + abline(h = -clim, col = '#FFFFFF', lwd = 2.85) + abline(h = -clim, col = 'black', lwd = 2.5, lty = 'dashed') + box(bty = 'L', lwd = 2) + title('ACF', line = 0) + + + ecdf_plotdat = function(vals, x){ + func <- ecdf(vals) + func(x) + } + + plot_x <- seq(min(truth, na.rm = T), + max(truth, na.rm = T)) + plot(1, type = "n", bty = 'L', + xlab = 'Count', + ylab = 'Empirical CDF', + xlim = c(min(plot_x), max(plot_x)), + ylim = c(0, 1)) + title('CDF', line = 0) + lines(x = plot_x, + y = ecdf_plotdat(truth, + plot_x), + col = "#8F2727", + lwd = 2.5) + box(bty = 'L', lwd = 2) + } + + layout(1) + } + +} + +dat <- sim_mvgam(T = 100, n_series = 4, n_lv = 2, + mu_obs = c(4, 6, 10, 14), trend_rel = 0.3, + seasonality = 'hierarchical') +plot_mvgam_series(data_train = dat$data_train, + n_bins = 20, + series = 'all', + log_scale = TRUE) # Good for testing model files without compiling stanc(model_code = mod1$model_file)$model_name diff --git a/R/add_stan_data.R b/R/add_stan_data.R index 1159ee57..90725caf 100644 --- a/R/add_stan_data.R +++ b/R/add_stan_data.R @@ -299,7 +299,7 @@ add_stan_data = function(jags_file, stan_file, use_lv = FALSE, sub("\\)", "", sub(".*\\:", "", stan_file[b_raw_indices[i]-1]))))), - ') {\nb[i] <- mu_raw[', i, '] + b_raw[i] * sigma_raw[',i, + ') {\nb[i] = mu_raw[', i, '] + b_raw[i] * sigma_raw[',i, '];\n}') min_beta[i] <- as.numeric(sub("for \\(i in ", "", sub("\\:.*", "", diff --git a/R/add_trend_lines.R b/R/add_trend_lines.R index 685d2c21..fe430e93 100644 --- a/R/add_trend_lines.R +++ b/R/add_trend_lines.R @@ -280,7 +280,7 @@ add_trend_lines = function(model_file, stan = FALSE, paste0('row_vector[num_basis] b_raw;\n\n// latent factor AR1 terms\nvector[n_lv] ar1;') model_file[grep('// dynamic factor estimates', model_file) + 6] <- - paste0('LV[2:n, j] ~ normal(ar1[j] * LV[1:(n - 1), j], sigma[s]);') + paste0('LV[2:n, j] ~ normal(ar1[j] * LV[1:(n - 1), j], 1);') model_file[grep('model \\{', model_file) + 2] <- paste0('\n// priors for AR parameters\nar1 ~ normal(0, 0.5);\n') diff --git a/R/plot.mvgam.R b/R/plot.mvgam.R index 08ef3bb2..fa85cf39 100644 --- a/R/plot.mvgam.R +++ b/R/plot.mvgam.R @@ -44,11 +44,11 @@ plot.mvgam = function(object, type = 'smooths', # Other errors and warnings will propagate from individual functions below if(type == 're'){ - plot_mvgam_randomeffects(object) + plot_mvgam_randomeffects(object, ...) } if(type == 'pterms'){ - plot_mvgam_pterms(object) + plot_mvgam_pterms(object, ...) } if(type == 'residuals'){ diff --git a/R/plot_mvgam_fc.R b/R/plot_mvgam_fc.R index 007f914c..3bfb9338 100644 --- a/R/plot_mvgam_fc.R +++ b/R/plot_mvgam_fc.R @@ -18,8 +18,10 @@ #'\code{realisations = TRUE}. Ignored otherwise #'@param hide_xlabels \code{logical}. If \code{TRUE}, no xlabels are printed to allow the user to add custom labels using #'\code{axis} from base \code{R} -#'@param ylab Optional \code{character} string specifying the y-axis label +#'@param xlab label for x axis. +#'@param ylab label for y axis. #'@param ylim Optional \code{vector} of y-axis limits (min, max) +#'@param ... further \code{\link[graphics]{par}} graphical parameters. #'@param return_forecasts \code{logical}. If \code{TRUE}, the function will plot the forecast #'as well as returning the forecast object (as a \code{matrix} of dimension \code{n_samples} x \code{horizon}) #'@details Posterior predictions are drawn from the fitted \code{mvgam} and used to calculate posterior @@ -30,8 +32,8 @@ #'@export plot_mvgam_fc = function(object, series = 1, data_test, realisations = FALSE, n_realisations = 15, - hide_xlabels = FALSE, ylab, ylim, - return_forecasts = FALSE){ + hide_xlabels = FALSE, xlab, ylab, ylim, + return_forecasts = FALSE, ...){ # Check arguments if(class(object) != 'mvgam'){ @@ -168,6 +170,10 @@ plot_mvgam_fc = function(object, series = 1, data_test, ylab <- paste0('Predicitons for ', levels(data_train$series)[series]) } + if(missing(xlab)){ + xlab <- 'Time' + } + pred_vals <- seq(1:length(preds_last)) if(hide_xlabels){ plot(1, type = "n", bty = 'L', @@ -175,13 +181,13 @@ plot_mvgam_fc = function(object, series = 1, data_test, xaxt = 'n', ylab = ylab, xlim = c(0, length(preds_last)), - ylim = ylim) + ylim = ylim, ...) } else { plot(1, type = "n", bty = 'L', - xlab = 'Time', + xlab = xlab, ylab = ylab, xlim = c(0, length(preds_last)), - ylim = ylim) + ylim = ylim, ...) } if(realisations){ @@ -223,20 +229,29 @@ plot_mvgam_fc = function(object, series = 1, data_test, time = data_test$time) } + # Show historical distribution in grey + last_train <- (NROW(data_train) / NCOL(object$ytimes)) + polygon(c(pred_vals[1:(NROW(data_train) / NCOL(object$ytimes))], + rev(pred_vals[1:(NROW(data_train) / NCOL(object$ytimes))])), + c(cred[1,1:(NROW(data_train) / NCOL(object$ytimes))], + rev(cred[9,1:(NROW(data_train) / NCOL(object$ytimes))])), + col = 'grey70', border = NA) + # Plot training and testing points points(dplyr::bind_rows(data_train, data_test) %>% dplyr::filter(series == s_name) %>% dplyr::select(time, y) %>% dplyr::distinct() %>% dplyr::arrange(time) %>% - dplyr::pull(y), pch = 16, col = "white", cex = 0.75) + dplyr::pull(y), pch = 16, col = "white", cex = 0.8) points(dplyr::bind_rows(data_train, data_test) %>% dplyr::filter(series == s_name) %>% dplyr::select(time, y) %>% dplyr::distinct() %>% dplyr::arrange(time) %>% dplyr::pull(y), pch = 16, col = "black", cex = 0.65) - abline(v = NROW(data_train) / NCOL(object$ytimes), lty = 'dashed') + abline(v = last_train, col = '#FFFFFF60', lwd = 2.85) + abline(v = last_train, col = 'black', lwd = 2.5, lty = 'dashed') # Calculate out of sample DRPS and print the score drps_score <- function(truth, fc, interval_width = 0.9){ @@ -306,13 +321,13 @@ plot_mvgam_fc = function(object, series = 1, data_test, dplyr::select(time, y) %>% dplyr::distinct() %>% dplyr::arrange(time) %>% - dplyr::pull(y),pch = 16, col = "white", cex = 0.75) + dplyr::pull(y),pch = 16, col = "white", cex = 0.8) points(data_train %>% dplyr::filter(series == s_name) %>% dplyr::select(time, y) %>% dplyr::distinct() %>% dplyr::arrange(time) %>% - dplyr::pull(y),pch = 16, col = "black", cex = 0.65 ) + dplyr::pull(y),pch = 16, col = "black", cex = 0.65) } if(return_forecasts){ diff --git a/R/plot_mvgam_resids.R b/R/plot_mvgam_resids.R index e7791fc7..3a0379f1 100644 --- a/R/plot_mvgam_resids.R +++ b/R/plot_mvgam_resids.R @@ -320,13 +320,14 @@ x <- sapply(1:length(idx), # Plot plot(median_preds[1:length(object$resids[[series]])], object$resids[[series]], - main = 'Resids vs Fitted Values', + bty = 'L', xlab = 'Fitted values', - ylab = 'Residuals', + 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)], @@ -358,8 +359,8 @@ for (k in 1:N) { y = c(resid_probs[k,5], resid_probs[k,5]), col = c_dark, lwd = 2) } -abline(h = 0, lty = 'dashed', col = 'white', lwd = 2.85) -abline(h = 0, lty = 'dashed', col = 'black', lwd = 2.5) +abline(h = 0, col = '#FFFFFF60', lwd = 2.85) +abline(h = 0, col = 'black', lwd = 2.5, lty = 'dashed') # Q-Q plot coords <- qqnorm(series_residuals[1,], plot.it = F) @@ -379,13 +380,15 @@ pred_vals <- coords$x[order(coords$x)] pred_vals <- pred_vals[complete.cases(cred[1,])] plot(x = pred_vals, y = cred[5,][complete.cases(cred[1,])], - main = 'Normal Q-Q Plot', + bty = 'L', xlab = 'Theoretical Quantiles', ylab = 'Sample Quantiles', pch = 16, col = 'white', cex = 1, - ylim = range(cred, na.rm = T)) + ylim = range(cred, na.rm = T), + tck = -0.04) +title('Normal Q-Q Plot', line = 0) 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) @@ -399,7 +402,7 @@ polygon(c(pred_vals, rev(pred_vals)), c(cred[4,][complete.cases(cred[1,])], rev(cred[6,][complete.cases(cred[1,])])), col = c_mid_highlight, border = NA) lines(pred_vals, cred[5,][complete.cases(cred[1,])], col = c_dark, lwd = 2.5) -qqline(cred[5,][complete.cases(cred[1,])], col = 'white', lwd = 3) +qqline(cred[5,][complete.cases(cred[1,])], col = '#FFFFFF60', lwd = 3) qqline(cred[5,][complete.cases(cred[1,])], col = 'black', lwd = 2.5) # ACF plot @@ -426,16 +429,16 @@ cred <- sapply(1:NCOL(resid_acf), probs = probs, na.rm = T)) cred <- cred[, -1] clim <- qnorm((1 + .95)/2)/sqrt(acf1$n.used) -plot(1, type = "n", +plot(1, type = "n", bty = 'L', xlab = 'Lag', - ylab = 'ACF', - main = 'ACF', + ylab = 'Autocorrelation', xlim = c(1, N-1), xaxt = 'n', ylim = range(c(cred, -clim - 0.05, clim + 0.05))) axis(1, at = seq(1, NCOL(cred), by = 2)) +title('ACF', line = 0) N <- N - 1 rect(xleft = x[seq(1, N*2, by = 2)], @@ -468,10 +471,10 @@ for (k in 1:N) { y = c(cred[5,k], cred[5,k]), col = c_dark, lwd = 2) } -abline(h = clim, lty = 'dashed', col = 'white', lwd = 2.85) -abline(h = clim, lty = 'dashed', col = 'black', lwd = 2.5) -abline(h = -clim, lty = 'dashed', col = 'white', lwd = 2.85) -abline(h = -clim, lty = 'dashed', col = 'black', lwd = 2.5) +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') # PACF plot pacf1 <- pacf(series_residuals[1,], plot = F, @@ -497,16 +500,16 @@ cred <- sapply(1:NCOL(resid_pacf), probs = probs, na.rm = T)) clim <- qnorm((1 + .95)/2)/sqrt(pacf1$n.used) -plot(1, type = "n", +plot(1, type = "n", bty = 'L', xlab = 'Lag', - ylab = 'Partial ACF', - main = 'pACF', + ylab = 'Autocorrelation', xlim = c(1, length(sorted_x)), xaxt = 'n', ylim = range(c(cred, -clim - 0.05, clim + 0.05))) axis(1, at = seq(1, NCOL(cred), by = 2)) +title('pACF', line = 0) rect(xleft = x[seq(1, N*2, by = 2)], xright = x[seq(2, N*2, by = 2)], @@ -538,11 +541,10 @@ for (k in 1:N) { y = c(cred[5,k], cred[5,k]), col = c_dark, lwd = 2) } -abline(h = clim, lty = 'dashed', col = 'white', lwd = 2.85) -abline(h = clim, lty = 'dashed', col = 'black', lwd = 2.5) -abline(h = -clim, lty = 'dashed', col = 'white', lwd = 2.85) -abline(h = -clim, lty = 'dashed', col = 'black', lwd = 2.5) - +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') layout(1) } diff --git a/R/plot_mvgam_trend.R b/R/plot_mvgam_trend.R index 8feb2f1f..894a3be9 100644 --- a/R/plot_mvgam_trend.R +++ b/R/plot_mvgam_trend.R @@ -12,10 +12,14 @@ #'\code{realisations = TRUE}. Ignored otherwise #'@param hide_xlabels \code{logical}. If \code{TRUE}, no xlabels are printed to allow the user to add custom labels using #'\code{axis} from base \code{R}. Ignored if \code{derivatives = TRUE} +#'@param xlab label for x axis. +#'@param ... further \code{\link[graphics]{par}} graphical parameters. #'@export plot_mvgam_trend = function(object, series = 1, data_test, realisations = FALSE, n_realisations = 15, - derivatives = FALSE, hide_xlabels = FALSE){ + derivatives = FALSE, hide_xlabels = FALSE, + xlab, + ...){ # Check arguments if(class(object) != 'mvgam'){ @@ -114,6 +118,10 @@ plot_mvgam_trend = function(object, series = 1, data_test, c_dark <- c("#8F2727") c_dark_highlight <- c("#7C0000") + if(missing(xlab)){ + xlab <- 'Time' + } + if(derivatives){ .pardefault <- par(no.readonly=T) par(.pardefault) @@ -122,10 +130,10 @@ plot_mvgam_trend = function(object, series = 1, data_test, mai = c(0.8, 0.8, 0.4, 0.4)) plot(1, type = "n", bty = 'L', - xlab = 'Time', + xlab = xlab, ylab = paste0('Estimated trend for ', levels(data_train$series)[series]), xlim = c(0, length(preds_last)), - ylim = range(cred)) + ylim = range(cred), ...) if(realisations){ for(i in 1:n_realisations){ @@ -157,9 +165,11 @@ plot_mvgam_trend = function(object, series = 1, data_test, if(!missing(data_test)){ if(class(data_train)[1] == 'list'){ - abline(v = length(data_train$y) / NCOL(object$ytimes), lty = 'dashed') + abline(v = length(data_train$y) / NCOL(object$ytimes), col = '#FFFFFF60', lwd = 2.85) + abline(v = length(data_train$y) / NCOL(object$ytimes), col = 'black', lwd = 2.5, lty = 'dashed') } else { - abline(v = NROW(data_train) / NCOL(object$ytimes), lty = 'dashed') + abline(v = NROW(data_train) / NCOL(object$ytimes), col = '#FFFFFF60', lwd = 2.85) + abline(v = NROW(data_train) / NCOL(object$ytimes), col = 'black', lwd = 2.5, lty = 'dashed') } } @@ -171,11 +181,11 @@ plot_mvgam_trend = function(object, series = 1, data_test, probs = probs, na.rm = T)) plot(1, type = "n", bty = "L", - xlab = 'Time', + xlab = xlab, ylab = '1st derivative', xlim = c(min(pred_vals), max(pred_vals)), ylim = c(min(cred, na.rm = T) - sd(first_derivs, na.rm = T), - max(cred, na.rm = T) + sd(first_derivs, na.rm = T))) + max(cred, na.rm = T) + sd(first_derivs, na.rm = T)), ...) if(realisations){ @@ -222,10 +232,10 @@ plot_mvgam_trend = function(object, series = 1, data_test, } else { plot(1, type = "n", bty = 'L', - xlab = 'Time', + xlab = xlab, ylab = paste0('Estimated trend for ', levels(data_train$series)[series]), xlim = c(0, length(preds_last)), - ylim = range(cred)) + ylim = range(cred), ...) } if(realisations){ @@ -258,9 +268,11 @@ plot_mvgam_trend = function(object, series = 1, data_test, if(!missing(data_test)){ if(class(data_train)[1] == 'list'){ - abline(v = length(data_train$y) / NCOL(object$ytimes), lty = 'dashed') + abline(v = length(data_train$y) / NCOL(object$ytimes), col = '#FFFFFF60', lwd = 2.85) + abline(v = length(data_train$y) / NCOL(object$ytimes), col = 'black', lwd = 2.5, lty = 'dashed') } else { - abline(v = NROW(data_train) / NCOL(object$ytimes), lty = 'dashed') + abline(v = NROW(data_train) / NCOL(object$ytimes), col = '#FFFFFF60', lwd = 2.85) + abline(v = NROW(data_train) / NCOL(object$ytimes), col = 'black', lwd = 2.5, lty = 'dashed') } } diff --git a/R/ppc.mvgam.R b/R/ppc.mvgam.R index 8fa1dadb..62448306 100644 --- a/R/ppc.mvgam.R +++ b/R/ppc.mvgam.R @@ -16,7 +16,11 @@ #'number of bins returned by a call to `hist` in base `R` #'@param legend_position The location may also be specified by setting x to a single keyword from the #'list "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" and "center". -#'This places the legend on the inside of the plot frame at the given location. +#'This places the legend on the inside of the plot frame at the given location. Or alternatively, +#'use "none" to hide the legend. +#'@param xlab label for x axis. +#'@param ylab label for y axis. +#'@param ... further \code{\link[graphics]{par}} graphical parameters. #'@details Posterior predictions are drawn from the fitted \code{mvgam} and compared against #'the empirical distribution of the observed data for a specified series to help evaluate the model's #'ability to generate unbiased predictions. For all plots apart from the 'rootogram', posterior predictions @@ -41,13 +45,15 @@ ppc <- function(x, what, ...){ #'@method ppc mvgam #'@export ppc.mvgam = function(object, data_test, series = 1, type = 'density', - n_bins, legend_position){ + n_bins, legend_position, xlab, ylab, ...){ # Check arguments type <- match.arg(arg = type, choices = c("rootogram", "mean", "hist", "density", "pit", "cdf", "prop_zero")) + optional_args <- list(...) + if(class(object) != 'mvgam'){ stop('argument "object" must be of class "mvgam"') } @@ -194,6 +200,15 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', pred_props <- pred_props[-which(pred_props < lower)] } obs_prop <- length(which(truths == 0)) / length(truths) + + if(missing(ylab)){ + ylab <- 'Density' + } + + if(missing(xlab)){ + xlab <- paste0('Predicted proportion of zeroes for ', levels(data_train$series)[series]) + } + hist(pred_props, lwd = 2, xlim = c(min(min(pred_props), min(obs_prop)), max(max(pred_props), max(obs_prop))), @@ -202,8 +217,9 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', length.out = 15), border = "#B97C7C", col = "#C79999", - ylab = 'Density', - xlab = paste0('Predicted proportion of zeroes for ', levels(data_train$series)[series])) + ylab = ylab, + xlab = xlab, + ...) abline(v = obs_prop, lwd = 3, col = 'white') abline(v = obs_prop, lwd = 2.5, col = 'black') box(bty = 'L', lwd = 2) @@ -212,14 +228,17 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', legend_position = 'topright' } - legend(legend_position, - legend = c(expression(hat(y)[propzero]), - expression(y[propzero])), - bg = 'white', - col = c(c_mid, - 'black'), - lty = 1, lwd = 2, - bty = 'n') + if(legend_position != 'none'){ + legend(legend_position, + legend = c(expression(hat(y)[propzero]), + expression(y[propzero])), + bg = 'white', + col = c(c_mid, + 'black'), + lty = 1, lwd = 2, + bty = 'n') + } + } if(type == 'rootogram'){ @@ -276,14 +295,23 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', repped_x[k] + min(diff(xpos))/2 else repped_x[k] - min(diff(xpos))/2) + if(missing(xlab)){ + xlab <- expression(y) + } + + if(missing(ylab)){ + ylab <- expression(sqrt(frequency)) + } + # Plot the rootogram plot(1, type = "n", bty = 'L', - xlab = expression(y), - ylab = expression(sqrt(frequency)), + xlab = xlab, + ylab = ylab, xlim = range(xpos), ylim = range(c(data$tyexp, data[,13], data[,5], - data$tyexp - data$ty))) + data$tyexp - data$ty)), + ...) rect(xleft = x[seq(1, N*2, by = 2)], xright = x[seq(2, N*2, by = 2)], ytop = data$tyexp, @@ -335,6 +363,15 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', pred_means <- pred_means[-which(pred_means < lower)] } obs_mean <- mean(truths) + + if(missing(ylab)){ + ylab <- 'Density' + } + + if(missing(xlab)){ + xlab <- paste0('Predicted mean for ', levels(data_train$series)[series]) + } + hist(pred_means, xlim = c(min(min(pred_means), min(obs_mean)), max(max(pred_means), max(obs_mean))), @@ -343,8 +380,9 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', breaks = seq(min(pred_means), max(pred_means), length.out = 20), border = "#B97C7C", col = "#C79999", - ylab = 'Density', - xlab = paste0('Predicted mean for ', levels(data_train$series)[series])) + ylab = ylab, + xlab = xlab, + ...) abline(v = obs_mean, lwd = 3, col = 'white') abline(v = obs_mean, lwd = 2.5, col = 'black') box(bty = 'L', lwd = 2) @@ -353,6 +391,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', legend_position = 'topright' } + if(legend_position != 'none'){ legend(legend_position, legend = c(expression(hat(mu)), expression(mu)), @@ -361,6 +400,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', 'black'), lty = 1, lwd = 2, bty = 'n') + } } @@ -384,11 +424,20 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', ymax <- max(c(max(cred), max(true_dens$y))) + if(missing(ylab)){ + ylab <- paste0('Predictive density for ', levels(data_train$series)[series]) + } + + if(missing(xlab)){ + xlab <- '' + } + plot(1, type = "n", bty = 'L', - xlab = '', - ylab = paste0('Predictive density for ', levels(data_train$series)[series]), + xlab = xlab, + ylab = ylab, xlim = c(min_x, max_x), - ylim = c(0, ymax)) + ylim = c(0, ymax), + ...) polygon(c(true_dens$x, rev(true_dens$x)), c(cred[1,], rev(cred[9,])), col = c_light, border = NA) @@ -407,6 +456,8 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', if(missing(legend_position)){ legend_position = 'topright' } + + if(legend_position != 'none'){ legend(legend_position, legend = c(expression(hat(y)), 'y'), @@ -416,6 +467,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', lty = 1, lwd = 2, bty = 'n') + } box(bty = 'L', lwd = 2) } @@ -436,21 +488,32 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', } bin_lims <- range(c(truths, as.vector(preds))) - delta <- diff(range(preds)) / n_bins - breaks <- seq(bin_lims[1], bin_lims[2] + delta, delta) + #delta <- diff(range(preds)) / n_bins + breaks <- seq(bin_lims[1], bin_lims[2], length.out = n_bins) xlim <- c(0, max(max(density(preds[1,])$x), max(density(truths)$x))) ylim <- c(0, max(c(max(hist(truths, breaks = breaks, plot = F)$density), max(hist(preds, breaks = breaks, plot = F)$density)))) + + if(missing(xlab)){ + xlab <- paste0('Count') + } + + if(missing(ylab)){ + ylab = '' + } + hist(preds, breaks=breaks, lwd = 2, main='', - xlab = paste0('Predictive histogram for ', levels(data_train$series)[series]), + xlab = xlab, + ylab = ylab, ylim=ylim, xlim=xlim, border = "#B97C7C", col = "#C79999", - freq = F) + freq = F, + ...) par(lwd=2) hist(truths, breaks=breaks, @@ -467,6 +530,8 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', if(missing(legend_position)){ legend_position = 'topright' } + + if(legend_position != 'none'){ legend(legend_position, legend = c(expression(hat(y)), 'y'), @@ -476,6 +541,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', lty = 1, lwd = 2, bty = 'n') + } } @@ -497,11 +563,20 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', function(n) quantile(pred_cdfs[,n], probs = probs)) + if(missing(ylab)){ + ylab = paste0('Predictive CDF for ', levels(data_train$series)[series]) + } + + if(missing(xlab)){ + xlab = '' + } + plot(1, type = "n", bty = 'L', - xlab = '', - ylab = paste0('Predictive CDF for ', levels(data_train$series)[series]), + xlab = xlab, + ylab = ylab, xlim = c(min(plot_x), max(plot_x)), - ylim = c(0, 1)) + ylim = c(0, 1), + ...) polygon(c(plot_x, rev(plot_x)), c(cred[1,], rev(cred[9,])), col = c_light, border = NA) @@ -528,6 +603,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', legend_position = 'bottomright' } + if(legend_position != 'none'){ legend(legend_position, legend = c(expression(hat(y)), 'y'), @@ -536,6 +612,7 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', 'black'), lty = 1, lwd = 2, bty = 'n') + } box(bty = 'L', lwd = 2) } @@ -563,8 +640,10 @@ ppc.mvgam = function(object, data_test, series = 1, type = 'density', barplot(pit_hist, lwd = 2, col = "#B97C7C", xlab = paste0('Predictive PIT for ', levels(data_train$series)[series]), - border = NA) - abline(h = 1, lty = 'dashed', lwd = 2) + border = NA, + ...) + abline(h = 1, col = '#FFFFFF60', lwd = 2.85) + abline(h = 1, col = 'black', lwd = 2.5, lty = 'dashed') box(bty = 'L', lwd = 2) } } diff --git a/R/sim_mvgam.R b/R/sim_mvgam.R index 3adcf36d..9d2db088 100644 --- a/R/sim_mvgam.R +++ b/R/sim_mvgam.R @@ -181,8 +181,6 @@ sim_mvgam = function(T = 100, trend_alphas <- runif(n_lv, 0.75, 1.2) trend_rhos <- runif(n_lv, 6, 12) - cat('trend_alphas = ', trend_alphas, '\n\n') - cat('trend_rhos = ', trend_rhos, '\n\n') # Generate latent GP trends trends <- do.call(cbind, lapply(seq_len(n_lv), function(lv){ set.seed(runif(1, 1, 100)) diff --git a/man/plot_mvgam_fc.Rd b/man/plot_mvgam_fc.Rd index e73e011f..a1af8d2b 100644 --- a/man/plot_mvgam_fc.Rd +++ b/man/plot_mvgam_fc.Rd @@ -11,9 +11,11 @@ plot_mvgam_fc( realisations = FALSE, n_realisations = 15, hide_xlabels = FALSE, + xlab, ylab, ylim, - return_forecasts = FALSE + return_forecasts = FALSE, + ... ) } \arguments{ @@ -42,12 +44,16 @@ empirical quantiles of the forecast distribution are shown} \item{hide_xlabels}{\code{logical}. If \code{TRUE}, no xlabels are printed to allow the user to add custom labels using \code{axis} from base \code{R}} -\item{ylab}{Optional \code{character} string specifying the y-axis label} +\item{xlab}{label for x axis.} + +\item{ylab}{label for y axis.} \item{ylim}{Optional \code{vector} of y-axis limits (min, max)} \item{return_forecasts}{\code{logical}. If \code{TRUE}, the function will plot the forecast as well as returning the forecast object (as a \code{matrix} of dimension \code{n_samples} x \code{horizon})} + +\item{...}{further \code{\link[graphics]{par}} graphical parameters.} } \value{ A base \code{R} graphics plot and an optional \code{matrix} of the forecast distribution diff --git a/man/plot_mvgam_trend.Rd b/man/plot_mvgam_trend.Rd index 87215f03..ee0c10bc 100644 --- a/man/plot_mvgam_trend.Rd +++ b/man/plot_mvgam_trend.Rd @@ -11,7 +11,9 @@ plot_mvgam_trend( realisations = FALSE, n_realisations = 15, derivatives = FALSE, - hide_xlabels = FALSE + hide_xlabels = FALSE, + xlab, + ... ) } \arguments{ @@ -34,6 +36,10 @@ estimated 1st derivative for the estimated trend} \item{hide_xlabels}{\code{logical}. If \code{TRUE}, no xlabels are printed to allow the user to add custom labels using \code{axis} from base \code{R}. Ignored if \code{derivatives = TRUE}} + +\item{xlab}{label for x axis.} + +\item{...}{further \code{\link[graphics]{par}} graphical parameters.} } \description{ Plot mvgam latent trend for a specified series diff --git a/man/ppc.mvgam.Rd b/man/ppc.mvgam.Rd index 27c8f3a7..72e264f7 100644 --- a/man/ppc.mvgam.Rd +++ b/man/ppc.mvgam.Rd @@ -4,7 +4,17 @@ \alias{ppc.mvgam} \title{Plot mvgam posterior predictive checks for a specified series} \usage{ -\method{ppc}{mvgam}(object, data_test, series = 1, type = "density", n_bins, legend_position) +\method{ppc}{mvgam}( + object, + data_test, + series = 1, + type = "density", + n_bins, + legend_position, + xlab, + ylab, + ... +) } \arguments{ \item{object}{\code{list} object returned from \code{mvgam}} @@ -28,7 +38,14 @@ number of bins returned by a call to \code{hist} in base \code{R}} \item{legend_position}{The location may also be specified by setting x to a single keyword from the list "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" and "center". -This places the legend on the inside of the plot frame at the given location.} +This places the legend on the inside of the plot frame at the given location. Or alternatively, +use "none" to hide the legend.} + +\item{xlab}{label for x axis.} + +\item{ylab}{label for y axis.} + +\item{...}{further \code{\link[graphics]{par}} graphical parameters.} } \value{ A base \code{R} graphics plot showing either a posterior rootogram (for \code{type == 'rootogram'}),