diff --git a/.github/workflows/R-CMD-check-rstan.yaml b/.github/workflows/R-CMD-check-rstan.yaml index c1674136..1e68bbe1 100644 --- a/.github/workflows/R-CMD-check-rstan.yaml +++ b/.github/workflows/R-CMD-check-rstan.yaml @@ -23,8 +23,10 @@ jobs: - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} - {os: ubuntu-latest, r: 'release'} +# Use a CRAN-like environment to emulate CRAN submission checks env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + NOT_CRAN: false R_KEEP_PKG_SOURCE: yes steps: diff --git a/R/mvgam.R b/R/mvgam.R index 20303ae8..8b119e2f 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -396,13 +396,13 @@ #' #' # Fit the model using AR1 trends #' mod <- mvgam(y ~ s(season, bs = 'cc', k = 6), -#' trend_map = trend_map, -#' trend_model = AR(), -#' data = mod_data, -#' return_model_data = TRUE, -#' burnin = 300, -#' samples = 300, -#' chains = 2) +#' trend_map = trend_map, +#' trend_model = AR(), +#' data = mod_data, +#' return_model_data = TRUE, +#' burnin = 300, +#' samples = 300, +#' chains = 2) #' #' # The mapping matrix is now supplied as data to the model in the 'Z' element #' mod1$model_data$Z @@ -525,71 +525,6 @@ #' layout(1) #' #' -#' # Example of logistic growth with possible changepoints -#' # Simple logistic growth model -#' dNt = function(r, N, k){ -#' r * N * (k - N) -#' } -#' -#' # Iterate growth through time -#' Nt = function(r, N, t, k) { -#' for (i in 1:(t - 1)) { -#' -#' # population at next time step is current population + growth, -#' # but we introduce several 'shocks' as changepoints -#' if(i %in% c(5, 15, 25, 41, 45, 60, 80)){ -#' N[i + 1] <- max(1, N[i] + dNt(r + runif(1, -0.1, 0.1), -#' N[i], k)) -#' } else { -#' N[i + 1] <- max(1, N[i] + dNt(r, N[i], k)) -#' } -#' } -#' N -#' } -#' -#' # Simulate expected values -#' set.seed(11) -#' expected <- Nt(0.004, 2, 100, 30) -#' plot(expected, xlab = 'Time') -#' -#' # Take Poisson draws -#' y <- rpois(100, expected) -#' plot(y, xlab = 'Time') -#' -#' # Assemble data into dataframe and model. We set a -#' # fixed carrying capacity of 35 for this example, but note that -#' # this value is not required to be fixed at each timepoint -#' mod_data <- data.frame(y = y, -#' time = 1:100, -#' cap = 35, -#' series = as.factor('series_1')) -#' plot_mvgam_series(data = mod_data) -#' -#' # The intercept is nonidentifiable when using piecewise -#' # trends because the trend functions have their own offset -#' # parameters 'm'; it is recommended to always drop intercepts -#' # when using these trend models -#' mod <- mvgam(y ~ 0, -#' trend_model = PW(growth = 'logistic'), -#' family = poisson(), -#' data = mod_data, -#' burnin = 300, -#' samples = 300, -#' chains = 2) -#' summary(mod) -#' -#' # Plot the posterior hindcast -#' plot(mod, type = 'forecast') -#' -#' # View the changepoints with ggplot2 utilities -#' library(ggplot2) -#' mcmc_plot(mod, variable = 'delta_trend', -#' regex = TRUE) + -#' scale_y_discrete(labels = mod$trend_model$changepoints) + -#' labs(y = 'Potential changepoint', -#' x = 'Rate change') -#' -#' #' # Example showcasing how cbind() is needed for Binomial observations #' # Simulate two time series of Binomial trials #' trials <- sample(c(20:25), 50, replace = TRUE) diff --git a/R/piecewise_trends.R b/R/piecewise_trends.R index f5997fad..08ffc6d8 100644 --- a/R/piecewise_trends.R +++ b/R/piecewise_trends.R @@ -47,7 +47,70 @@ #' no missing values are allowed in `cap`. #' #' @rdname piecewise_trends +#' @examples +#' \dontrun{ +#' # Example of logistic growth with possible changepoints +#' # Simple logistic growth model +#' dNt = function(r, N, k){ +#' r * N * (k - N) +#' } #' +#' # Iterate growth through time +#' Nt = function(r, N, t, k) { +#' for (i in 1:(t - 1)) { +#' +#' # population at next time step is current population + growth, +#' # but we introduce several 'shocks' as changepoints +#' if(i %in% c(5, 15, 25, 41, 45, 60, 80)){ +#' N[i + 1] <- max(1, N[i] + dNt(r + runif(1, -0.1, 0.1), +#' N[i], k)) +#' } else { +#' N[i + 1] <- max(1, N[i] + dNt(r, N[i], k)) +#' } +#' } +#' N +#' } +#' +#' # Simulate expected values +#' set.seed(11) +#' expected <- Nt(0.004, 2, 100, 30) +#' plot(expected, xlab = 'Time') +#' +#' # Take Poisson draws +#' y <- rpois(100, expected) +#' plot(y, xlab = 'Time') +#' +#' # Assemble data into dataframe and model. We set a +#' # fixed carrying capacity of 35 for this example, but note that +#' # this value is not required to be fixed at each timepoint +#' mod_data <- data.frame(y = y, +#' time = 1:100, +#' cap = 35, +#' series = as.factor('series_1')) +#' plot_mvgam_series(data = mod_data) +#' +#' # The intercept is nonidentifiable when using piecewise +#' # trends because the trend functions have their own offset +#' # parameters 'm'; it is recommended to always drop intercepts +#' # when using these trend models +#' mod <- mvgam(y ~ 0, +#' trend_model = PW(growth = 'logistic'), +#' family = poisson(), +#' data = mod_data, +#' chains = 2) +#' summary(mod) +#' +#' # Plot the posterior hindcast +#' plot(mod, type = 'forecast') +#' +#' # View the changepoints with ggplot2 utilities +#' library(ggplot2) +#' mcmc_plot(mod, variable = 'delta_trend', +#' regex = TRUE) + +#' scale_y_discrete(labels = mod$trend_model$changepoints) + +#' labs(y = 'Potential changepoint', +#' x = 'Rate change') +#' } #' @export PW = function(n_changepoints = 10, changepoint_range = 0.8, diff --git a/man/mvgam.Rd b/man/mvgam.Rd index 99dfadf2..8f7aafd6 100644 --- a/man/mvgam.Rd +++ b/man/mvgam.Rd @@ -474,13 +474,13 @@ trend_map <- data.frame(series = unique(mod_data$series), # Fit the model using AR1 trends mod <- mvgam(y ~ s(season, bs = 'cc', k = 6), - trend_map = trend_map, - trend_model = AR(), - data = mod_data, - return_model_data = TRUE, - burnin = 300, - samples = 300, - chains = 2) + trend_map = trend_map, + trend_model = AR(), + data = mod_data, + return_model_data = TRUE, + burnin = 300, + samples = 300, + chains = 2) # The mapping matrix is now supplied as data to the model in the 'Z' element mod1$model_data$Z @@ -603,71 +603,6 @@ for(i in 2:50){ layout(1) -# Example of logistic growth with possible changepoints -# Simple logistic growth model -dNt = function(r, N, k){ - r * N * (k - N) -} - -# Iterate growth through time -Nt = function(r, N, t, k) { -for (i in 1:(t - 1)) { - - # population at next time step is current population + growth, - # but we introduce several 'shocks' as changepoints - if(i \%in\% c(5, 15, 25, 41, 45, 60, 80)){ - N[i + 1] <- max(1, N[i] + dNt(r + runif(1, -0.1, 0.1), - N[i], k)) - } else { - N[i + 1] <- max(1, N[i] + dNt(r, N[i], k)) - } - } - N -} - -# Simulate expected values -set.seed(11) -expected <- Nt(0.004, 2, 100, 30) -plot(expected, xlab = 'Time') - -# Take Poisson draws -y <- rpois(100, expected) -plot(y, xlab = 'Time') - -# Assemble data into dataframe and model. We set a -# fixed carrying capacity of 35 for this example, but note that -# this value is not required to be fixed at each timepoint -mod_data <- data.frame(y = y, - time = 1:100, - cap = 35, - series = as.factor('series_1')) -plot_mvgam_series(data = mod_data) - -# The intercept is nonidentifiable when using piecewise -# trends because the trend functions have their own offset -# parameters 'm'; it is recommended to always drop intercepts -# when using these trend models -mod <- mvgam(y ~ 0, - trend_model = PW(growth = 'logistic'), - family = poisson(), - data = mod_data, - burnin = 300, - samples = 300, - chains = 2) -summary(mod) - -# Plot the posterior hindcast -plot(mod, type = 'forecast') - -# View the changepoints with ggplot2 utilities -library(ggplot2) -mcmc_plot(mod, variable = 'delta_trend', - regex = TRUE) + -scale_y_discrete(labels = mod$trend_model$changepoints) + -labs(y = 'Potential changepoint', - x = 'Rate change') - - # Example showcasing how cbind() is needed for Binomial observations # Simulate two time series of Binomial trials trials <- sample(c(20:25), 50, replace = TRUE) diff --git a/man/piecewise_trends.Rd b/man/piecewise_trends.Rd index aba9914c..bb3bbfd0 100644 --- a/man/piecewise_trends.Rd +++ b/man/piecewise_trends.Rd @@ -65,6 +65,71 @@ transformed if you are using a \code{poisson()} or \code{nb()} family). It is th that you specify the \code{cap} values on the scale of your outcome. Note also that no missing values are allowed in \code{cap}. } +\examples{ +\dontrun{ +# Example of logistic growth with possible changepoints +# Simple logistic growth model +dNt = function(r, N, k){ + r * N * (k - N) +} + +# Iterate growth through time +Nt = function(r, N, t, k) { +for (i in 1:(t - 1)) { + + # population at next time step is current population + growth, + # but we introduce several 'shocks' as changepoints + if(i \%in\% c(5, 15, 25, 41, 45, 60, 80)){ + N[i + 1] <- max(1, N[i] + dNt(r + runif(1, -0.1, 0.1), + N[i], k)) + } else { + N[i + 1] <- max(1, N[i] + dNt(r, N[i], k)) + } + } + N +} + +# Simulate expected values +set.seed(11) +expected <- Nt(0.004, 2, 100, 30) +plot(expected, xlab = 'Time') + +# Take Poisson draws +y <- rpois(100, expected) +plot(y, xlab = 'Time') + +# Assemble data into dataframe and model. We set a +# fixed carrying capacity of 35 for this example, but note that +# this value is not required to be fixed at each timepoint +mod_data <- data.frame(y = y, + time = 1:100, + cap = 35, + series = as.factor('series_1')) +plot_mvgam_series(data = mod_data) + +# The intercept is nonidentifiable when using piecewise +# trends because the trend functions have their own offset +# parameters 'm'; it is recommended to always drop intercepts +# when using these trend models +mod <- mvgam(y ~ 0, + trend_model = PW(growth = 'logistic'), + family = poisson(), + data = mod_data, + chains = 2) +summary(mod) + +# Plot the posterior hindcast +plot(mod, type = 'forecast') + +# View the changepoints with ggplot2 utilities +library(ggplot2) +mcmc_plot(mod, variable = 'delta_trend', + regex = TRUE) + +scale_y_discrete(labels = mod$trend_model$changepoints) + +labs(y = 'Potential changepoint', + x = 'Rate change') +} +} \references{ Taylor, Sean J., and Benjamin Letham. "Forecasting at scale." The American Statistician 72.1 (2018): 37-45. } diff --git a/src/mvgam.dll b/src/mvgam.dll index c1becdd1..b908ba0c 100644 Binary files a/src/mvgam.dll and b/src/mvgam.dll differ