Generate Stan code and data objects for mvgam models
Arguments
- object
An object of class
mvgamormvgam_prefit, returned from a call tomvgam- ...
ignored
Value
Either a character string containing the fully commented Stan code to fit a mvgam model or a named list containing the data objects needed to fit the model in Stan.
Examples
simdat <- sim_mvgam()
mod <- mvgam(y ~ s(season) +
s(time, by = series),
family = poisson(),
data = simdat$data_train,
run_model = FALSE)
# View Stan model code
stancode(mod)
#> // Stan model code generated by package mvgam
#> data {
#> int<lower=0> total_obs; // total number of observations
#> int<lower=0> n; // number of timepoints per series
#> int<lower=0> n_sp; // number of smoothing parameters
#> int<lower=0> n_series; // number of series
#> int<lower=0> num_basis; // total number of basis coefficients
#> vector[num_basis] zero; // prior locations for basis coefficients
#> matrix[total_obs, num_basis] X; // mgcv GAM design matrix
#> array[n, n_series] int<lower=0> ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?)
#> matrix[9, 18] S1; // mgcv smooth penalty matrix S1
#> matrix[9, 18] S2; // mgcv smooth penalty matrix S2
#> matrix[9, 18] S3; // mgcv smooth penalty matrix S3
#> matrix[9, 18] S4; // mgcv smooth penalty matrix S4
#> int<lower=0> n_nonmissing; // number of nonmissing observations
#> array[n_nonmissing] int<lower=0> flat_ys; // flattened nonmissing observations
#> matrix[n_nonmissing, num_basis] flat_xs; // X values for nonmissing observations
#> array[n_nonmissing] int<lower=0> obs_ind; // indices of nonmissing observations
#> }
#> parameters {
#> // raw basis coefficients
#> vector[num_basis] b_raw;
#>
#> // smoothing parameters
#> vector<lower=0>[n_sp] lambda;
#> }
#> transformed parameters {
#> // basis coefficients
#> vector[num_basis] b;
#> b[1 : num_basis] = b_raw[1 : num_basis];
#> }
#> model {
#> // prior for (Intercept)...
#> b_raw[1] ~ student_t(3, 0, 2.5);
#>
#> // prior for s(season)...
#> b_raw[2 : 10] ~ multi_normal_prec(zero[2 : 10],
#> S1[1 : 9, 1 : 9] * lambda[1]
#> + S1[1 : 9, 10 : 18] * lambda[2]);
#>
#> // prior for s(time):seriesseries_1...
#> b_raw[11 : 19] ~ multi_normal_prec(zero[11 : 19],
#> S2[1 : 9, 1 : 9] * lambda[3]
#> + S2[1 : 9, 10 : 18] * lambda[4]);
#>
#> // prior for s(time):seriesseries_2...
#> b_raw[20 : 28] ~ multi_normal_prec(zero[20 : 28],
#> S3[1 : 9, 1 : 9] * lambda[5]
#> + S3[1 : 9, 10 : 18] * lambda[6]);
#>
#> // prior for s(time):seriesseries_3...
#> b_raw[29 : 37] ~ multi_normal_prec(zero[29 : 37],
#> S4[1 : 9, 1 : 9] * lambda[7]
#> + S4[1 : 9, 10 : 18] * lambda[8]);
#>
#> // priors for smoothing parameters
#> lambda ~ normal(5, 30);
#> {
#> // likelihood functions
#> flat_ys ~ poisson_log_glm(flat_xs, 0.0, b);
#> }
#> }
#> generated quantities {
#> vector[total_obs] eta;
#> matrix[n, n_series] mus;
#> vector[n_sp] rho;
#> array[n, n_series] int ypred;
#> rho = log(lambda);
#>
#> // posterior predictions
#> eta = X * b;
#> for (s in 1 : n_series) {
#> mus[1 : n, s] = eta[ytimes[1 : n, s]];
#> ypred[1 : n, s] = poisson_log_rng(mus[1 : n, s]);
#> }
#> }
#>
#>
# View Stan model data
sdata <- standata(mod)
str(sdata)
#> List of 21
#> $ y : num [1:75, 1:3] 1 1 1 0 3 0 0 0 0 0 ...
#> $ n : int 75
#> $ X : num [1:225, 1:37] 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:37] "X.Intercept." "V2" "V3" "V4" ...
#> $ S1 : num [1:9, 1:18] 3.758 -0.625 -1.816 0.191 2.357 ...
#> $ zero : num [1:37] 0 0 0 0 0 0 0 0 0 0 ...
#> $ S2 : num [1:9, 1:18] 8.555 -1.22 4.352 0.822 -5.594 ...
#> $ S3 : num [1:9, 1:18] 8.555 -1.22 4.352 0.822 -5.594 ...
#> $ S4 : num [1:9, 1:18] 8.555 -1.22 4.352 0.822 -5.594 ...
#> $ p_coefs : Named num 0
#> ..- attr(*, "names")= chr "(Intercept)"
#> $ p_taus : num 0.988
#> $ ytimes : int [1:75, 1:3] 1 4 7 10 13 16 19 22 25 28 ...
#> $ n_series : int 3
#> $ sp : Named num [1:8] 0.368 0.368 0.368 0.368 0.368 ...
#> ..- attr(*, "names")= chr [1:8] "s(season)1" "s(season)2" "s(time):seriesseries_11" "s(time):seriesseries_12" ...
#> $ y_observed : num [1:75, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
#> $ total_obs : int 225
#> $ num_basis : int 37
#> $ n_sp : num 8
#> $ n_nonmissing: int 225
#> $ obs_ind : int [1:225] 1 2 3 4 5 6 7 8 9 10 ...
#> $ flat_ys : num [1:225] 1 1 1 0 3 0 0 0 0 0 ...
#> $ flat_xs : num [1:225, 1:37] 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:37] "X.Intercept." "V2" "V3" "V4" ...
#> - attr(*, "trend_model")= chr "None"