if (FALSE) {
# A short example to illustrate CAR(1) models
# Function to simulate CAR1 data with seasonality
sim_corcar1 = function(n = 125,
phi = 0.5,
sigma = 2,
sigma_obs = 0.75) {
# Sample irregularly spaced time intervals
time_dis <- c(1, runif(n - 1, 0, 5))
# Set up the latent dynamic process
x <- vector(length = n); x[1] <- -0.3
for (i in 2:n) {
# zero-distances will cause problems in sampling, so mvgam uses a
# minimum threshold; this simulation function emulates that process
if (time_dis[i] == 0) {
x[i] <- rnorm(
1,
mean = (phi^1e-3) * x[i - 1],
sd = sigma * (1 - phi^(2 * 1e-3)) / (1 - phi^2)
)
} else {
x[i] <- rnorm(
1,
mean = (phi^time_dis[i]) * x[i - 1],
sd = sigma * (1 - phi^(2 * time_dis[i])) / (1 - phi^2)
)
}
}
# Add 12-month seasonality
cov1 <- sin(2 * pi * (1:n) / 12)
cov2 <- cos(2 * pi * (1:n) / 12)
beta1 <- runif(1, 0.3, 0.7)
beta2 <- runif(1, 0.2, 0.5)
seasonality <- beta1 * cov1 + beta2 * cov2
# Take Gaussian observations with error and return
data.frame(
y = rnorm(n, mean = x + seasonality, sd = sigma_obs),
season = rep(1:12, 20)[1:n],
time = cumsum(time_dis)
)
}
# Sample two time series
dat <- rbind(
dplyr::bind_cols(
sim_corcar1(phi = 0.65, sigma_obs = 0.55),
data.frame(series = 'series1')
),
dplyr::bind_cols(
sim_corcar1(phi = 0.8, sigma_obs = 0.35),
data.frame(series = 'series2')
)
) %>%
dplyr::mutate(series = as.factor(series))
# mvgam with CAR(1) trends and series-level seasonal smooths
mod <- mvgam(
formula = y ~ -1,
trend_formula = ~ s(season, bs = 'cc', k = 5, by = trend),
trend_model = CAR(),
priors = c(
prior(exponential(3), class = sigma),
prior(beta(4, 4), class = sigma_obs)
),
data = dat,
family = gaussian(),
chains = 2,
silent = 2
)
# View usual summaries and plots
summary(mod)
conditional_effects(mod, type = 'expected')
plot(mod, type = 'trend', series = 1)
plot(mod, type = 'trend', series = 2)
plot(mod, type = 'residuals', series = 1)
plot(mod, type = 'residuals', series = 2)
mcmc_plot(
mod,
variable = 'ar1',
regex = TRUE,
type = 'hist'
)
# Now an example illustrating hierarchical dynamics
set.seed(123)
# Simulate three species monitored in three different regions
simdat1 <- sim_mvgam(
trend_model = VAR(cor = TRUE),
prop_trend = 0.95,
n_series = 3,
mu = c(1, 2, 3)
)
simdat2 <- sim_mvgam(
trend_model = VAR(cor = TRUE),
prop_trend = 0.95,
n_series = 3,
mu = c(1, 2, 3)
)
simdat3 <- sim_mvgam(
trend_model = VAR(cor = TRUE),
prop_trend = 0.95,
n_series = 3,
mu = c(1, 2, 3)
)
# Set up the data but DO NOT include 'series'
all_dat <- rbind(
simdat1$data_train %>%
dplyr::mutate(region = 'qld'),
simdat2$data_train %>%
dplyr::mutate(region = 'nsw'),
simdat3$data_train %>%
dplyr::mutate(region = 'vic')
) %>%
dplyr::mutate(
species = gsub('series', 'species', series),
species = as.factor(species),
region = as.factor(region)
) %>%
dplyr::arrange(series, time) %>%
dplyr::select(-series)
# Check priors for a hierarchical AR1 model
get_mvgam_priors(
formula = y ~ species,
trend_model = AR(gr = region, subgr = species),
data = all_dat
)
# Fit the model
mod <- mvgam(
formula = y ~ species,
trend_model = AR(gr = region, subgr = species),
data = all_dat,
chains = 2,
silent = 2
)
# Check standard outputs
summary(mod)
# Inspect posterior estimates for the correlation weighting parameter
mcmc_plot(mod, variable = 'alpha_cor', type = 'hist')
}
Run the code above in your browser using DataLab