Multivariate Model

I am working on multivariate linear regression model in Stan I run the code, however, I have some questions:

  1. I run the first model without scale() function and second with scale() function, is there a criteria or a specific case we should use the scale() function?

  2. In multivariate linear regression model, there is two popular prior distribution can be applied either LKJ or Wishart (inv-wishart) which one is better to use? In the same time, LKJ(n) how to select the n value for example, n=1, or 2, or, 3,……?

  3. How to take a log for y’s inside the stan model?

  4. Which one the more efficient in the multivariate model is Cholesky_factor_corr or Cholesky_factor_cov?

  5. How to generate y_prediction for all y’s in the multivariate model?

Thank you

Ali

1 Like

Hi! Can you post the model code and model call? And if you can, a subset of the data or a simulated data set would be helpful.

1 Like

R code

dat <- within(list(), {
  y <- as.matrix(dplyr::select(df, F1, F2, F3, F4))
  X <- model.matrix(~ 0 + Total_Fines + TOC + Total_solids + Total_sand + Total_gravel,
                    data = df) 
  N <- nrow(y)
  K <- ncol(y)
  P <- ncol(X)
  alpha_loc <- rep(0, K)
  alpha_scale <- rep(10, K)
  beta_loc <- matrix(0, K, P)
  beta_scale <- matrix(2.5, K, P)
  Sigma_corr_shape <- 2
  Sigma_scale_scale <- 5
  
})
library(rstan)
fit.model <- stan("stanTest.stan", data = dat, iter = 20000, warmup = 5000, chains = 3)
Stan code
data {
  // multivariate outcome
  int N;
  int K;
  vector[K] y[N];
  // covariates
  int P;
  vector[P] X[N];
  
  // prior
  vector[K] alpha_loc;
  vector[K] alpha_scale;
  vector[P] beta_loc[K];
  vector[P] beta_scale[K];
  real Sigma_corr_shape;
  real Sigma_scale_scale;
}
parameters {
  // regression intercept
  vector[K] alpha;
  // regression coefficients
  vector[P] beta[K];
  // Cholesky factor of the correlation matrix
  cholesky_factor_corr[K] Sigma_corr_L;
  vector<lower=0>[K] Sigma_scale;
  // student-T degrees of freedom
  
}
transformed parameters {
  vector[K] mu[N];
  matrix[K, K] Sigma;
  // covariance matrix
  Sigma = diag_pre_multiply(Sigma_scale, Sigma_corr_L);
  for (i in 1:N) {
    for (k in 1:K) {
      mu[i, k] = alpha[k] + dot_product(X[i], beta[k]);
    }
  }
}
model {
  for (k in 1:K) {
    alpha[k] ~ normal(alpha_loc[k], alpha_scale[k]);
    beta[k] ~ normal(beta_loc[k], beta_scale[k]);
  }
  
  Sigma_scale ~ cauchy(0., Sigma_scale_scale);
  Sigma_corr_L ~ lkj_corr_cholesky(Sigma_corr_shape);
  y ~ multi_normal_cholesky(mu, Sigma);

}

|test.csv (2.9 KB)

1 Like