Make stancode for the jsdm model
jsdm_stancode.RdThis function returns the Stan code used to fit the model as specified by the data list, family and method.
Usage
jsdm_stancode(
method,
family,
prior = jsdm_prior(),
log_lik = TRUE,
site_intercept = "none",
beta_param = "cor",
zi_param = "constant"
)
# S3 method for class 'jsdmstan_model'
print(x, ...)Arguments
- method
The method, one of
"gllvm"or"mglmm"- family
is the response family, must be one of
"gaussian","neg_binomial","poisson","binomial","bernoulli","zi_poisson", or"zi_neg_binomial". Regular expression matching is supported.- prior
The prior, given as the result of a call to
jsdm_prior()- log_lik
Whether the log likelihood should be calculated in the generated quantities (by default
TRUE), required for loo- site_intercept
Whether a site intercept should be included, potential values
"none"(no site intercept),"grouped"(a site intercept with hierarchical grouping) or"ungrouped"(site intercept with no grouping)- beta_param
The parameterisation of the environmental covariate effects, by default
"cor". See details for further information.- zi_param
For the zero-inflated families, whether the zero-inflation parameter is a species-specific constant (default,
"constant"), or varies by environmental covariates ("covariate").- x
The jsdm_stancode object
- ...
Currently unused
Details
Environmental covariate effects ("betas") can be parameterised in two
ways. With the "cor" parameterisation all covariate effects are assumed
to be constrained by a correlation matrix between the covariates. With the
"unstruct" parameterisation all covariate effects are assumed to draw
from a simple distribution with no correlation structure. Both parameterisations
can be modified using the prior object.
Examples
jsdm_stancode(family = "gaussian", method = "gllvm")
#> //Generated by jsdmstan
#> functions{
#>
#>
#> }
#> data{
#> int<lower=1> N; // Number of sites
#> int<lower=1> S; // Number of species
#> int<lower=1> D; // Number of latent dimensions
#> int<lower=0> K; // Number of predictor variables
#> matrix[N, K] X; // Predictor matrix
#> real Y[N,S]; //Species matrix
#> }
#> transformed data{
#>
#> // Ensures identifiability of the model - no rotation of factors
#> int<lower=1> M;
#> M = D * (S - D) + choose(D, 2) + D;
#>
#> }
#> parameters{
#>
#> //betas are hierarchical with covariance model
#> vector<lower=0>[K] sigmas_preds;
#> matrix[K, S] z_preds;
#> // covariance matrix on betas by predictors
#> corr_matrix[K] cor_preds;
#> // Factor parameters
#> vector[M] L; // Non-zero factor loadings
#> real<lower=0> sigma_L; // variance of species loadings
#> // Latent variables
#> matrix[D, N] LV_uncor; // Per-site latent variable
#> real<lower=0> sigma[S]; // Gaussian parameters
#> }
#> transformed parameters{
#>
#> // covariance matrix on betas by preds
#> matrix[K, S] betas;
#>
#> // Construct factor loading matrix
#> matrix[S, D] Lambda_uncor;
#> // Constraints to allow identifiability of loadings
#> for (i in 1:(D-1)) {
#> for (j in (i+1):(D)){
#> Lambda_uncor[i,j] = 0;
#> }
#> }
#> {
#> int index;
#> index = 0;
#> for (j in 1:D) {
#> for (i in j:S) {
#> index = index + 1;
#> Lambda_uncor[i, j] = L[index];
#> }
#> }
#> }
#>
#> betas = diag_pre_multiply(sigmas_preds, cor_preds) * z_preds;
#>
#> }
#> model{
#>
#> matrix[N,S] mu;
#>
#> // model
#> matrix[N, S] LV_sum = ((Lambda_uncor * sigma_L) * LV_uncor)';
#> mu = (X * betas) + LV_sum;
#>
#>
#>
#> // Species parameter priors
#> sigmas_preds ~ normal(0,1) ;
#> to_vector(z_preds) ~ normal(0,1) ;
#> // covariance matrix priors
#> cor_preds ~ lkj_corr(1) ;
#>
#> // Factor priors
#> to_vector(LV_uncor) ~ normal(0,1) ;
#> L ~ normal(0,1) ;
#> sigma_L ~ normal(0,1) ; // Variance of factor loadings
#>
#> //Standard deviation parameters
#> sigma ~ normal(0,1) ;
#>
#>
#> for(i in 1:N) Y[i,] ~ normal(mu[i,], sigma);
#>
#> }
#> generated quantities{
#>
#> // Calculate linear predictor, y_rep, log likelihoods for LOO
#> matrix[N, S] log_lik;
#>
#> // Sign correct factor loadings and factors
#> matrix[D, N] LV;
#> matrix[S, D] Lambda;
#> for(d in 1:D){
#> if(Lambda_uncor[d,d] < 0){
#> Lambda[,d] = -1 * Lambda_uncor[,d];
#> LV[d,] = -1 * LV_uncor[d,];
#> } else {
#> Lambda[,d] = Lambda_uncor[,d];
#> LV[d,] = LV_uncor[d,];
#> }
#> }
#> {
#> matrix[N, S] linpred;
#> linpred = (X * betas) + ((Lambda_uncor * sigma_L) * LV_uncor)' ;
#>
#> for(i in 1:N) {
#> for(j in 1:S) { log_lik[i, j] = normal_lpdf(Y[i, j] | linpred[i, j], sigma[j]);
#> }
#> }
#> }
#>
#> }
#>
jsdm_stancode(family = "poisson", method = "mglmm")
#> //Generated by jsdmstan
#> functions{
#>
#>
#> }
#> data{
#> int<lower=1> N; // Number of sites
#> int<lower=1> S; // Number of species
#>
#> int<lower=0> K; // Number of predictor variables
#> matrix[N, K] X; // Predictor matrix
#> int<lower=0> Y[N,S]; //Species matrix
#> }
#> transformed data{
#>
#> }
#> parameters{
#>
#> //betas are hierarchical with covariance model
#> vector<lower=0>[K] sigmas_preds;
#> matrix[K, S] z_preds;
#> // covariance matrix on betas by predictors
#> corr_matrix[K] cor_preds;
#> // species covariances
#> vector<lower=0>[S] sigmas_species;
#> matrix[S, N] z_species;
#> corr_matrix[S] cor_species;
#> }
#> transformed parameters{
#>
#> // covariance matrix on betas by preds
#> matrix[K, S] betas;
#>
#> matrix[N, S] u;
#> u = (diag_pre_multiply(sigmas_species, cor_species) * z_species)';
#>
#> betas = diag_pre_multiply(sigmas_preds, cor_preds) * z_preds;
#>
#> }
#> model{
#>
#> matrix[N,S] mu;
#>
#> // model
#> mu = (X * betas) + u;
#>
#>
#>
#> // Species parameter priors
#> sigmas_preds ~ normal(0,1) ;
#> to_vector(z_preds) ~ normal(0,1) ;
#> // covariance matrix priors
#> cor_preds ~ lkj_corr(1) ;
#>
#> // Species parameter priors
#> sigmas_species ~ normal(0,1) ;
#> to_vector(z_species) ~ normal(0,1) ;
#> cor_species ~ lkj_corr(1) ;
#>
#>
#> for(i in 1:N) Y[i,] ~ poisson_log(mu[i,]);
#>
#> }
#> generated quantities{
#>
#> // Calculate linear predictor, y_rep, log likelihoods for LOO
#> matrix[N, S] log_lik;
#>
#> {
#> matrix[N, S] linpred;
#> linpred = (X * betas) + u ;
#>
#> for(i in 1:N) {
#> for(j in 1:S) { log_lik[i, j] = poisson_log_lpmf(Y[i, j] | linpred[i, j]);
#> }
#> }
#> }
#>
#> }
#>