I am working on multivariate linear regression model in Stan I run the code, however, I have some questions:
-
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?
-
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,……?
-
How to take a log for y’s inside the stan model?
-
Which one the more efficient in the multivariate model is Cholesky_factor_corr or Cholesky_factor_cov?
-
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