Hi, I’m extremely new to using covariance matrices in Stan, and I’m running into some errors. My code is as follows:
// Define the beta-binomial regression model
data {
int<lower=0> N; // Number of adjectives
int<lower=0> y_comparative[N]; // Observed morphological comparative
int<lower=0> y_superlative[N]; // Observed morphological superlative
real FreqNonCentered_lemma[N]; // Lemma frequency (summed log freq of each)
real FreqNonCentered_comparative[N]; // Comparative frequency
real FreqNonCentered_superlative[N]; // Superlative frequency
int<lower=0> n_comparative[N]; // Total number of comparative forms
int<lower=0> n_superlative[N]; // Total number of superlative forms
int<lower=0> K; // Number of predictors
matrix[N, K] X; // Design matrix
matrix[2, N] mu; // Location parameter for bivariate_offsets
}
parameters {
vector[K] shared_betas; // Regression coefficients
real beta_freq;
real alpha;
matrix[2, N] bivariate_offsets;
real rho;
}
model {
vector[N] eta_comparative; // Probability of success for comparative
vector[N] eta_superlative; // Probability of success for superlative
vector[N] sigma_comparative; // Standard deviation for comparative
vector[N] sigma_superlative; // Standard deviation for superlative
matrix[2, 2] Sigmas[N]; // Covariance matrix for each observation
rho ~ uniform(-1, 1);
// Priors
shared_betas ~ normal(0, 10); // Prior for regression coefficients
alpha ~ normal(0, 10);
beta_freq ~ normal(0, 10);
// Calculate the standard deviations and covariance matrices
for (i in 1:N) {
sigma_comparative[i] = exp(alpha + FreqNonCentered_comparative[i] * beta_freq);
sigma_superlative[i] = exp(alpha + FreqNonCentered_superlative[i] * beta_freq);
Sigmas[i][1, 1] = square(sigma_comparative[i]);
Sigmas[i][2, 2] = square(sigma_superlative[i]);
Sigmas[i][1, 2] = sigma_comparative[i] * rho * sigma_superlative[i];
Sigmas[i][2, 1] = sigma_comparative[i] * rho * sigma_superlative[i];
bivariate_offsets[,i] ~ multi_normal(mu[,i], Sigmas[i]);
}
// Priors for bivariate_offsets
// Calculate the probability of success for each observation
eta_comparative = X * shared_betas + to_vector(bivariate_offsets[1, ]);
eta_superlative = X * shared_betas + to_vector(bivariate_offsets[2, ]);
// Likelihood
y_comparative ~ binomial_logit(n_comparative, eta_comparative);
y_superlative ~ binomial_logit(n_superlative, eta_superlative);
}
My data are
data.csv (26.8 KB)
.
My R code is here:
stan_model_file <- "./Models/bivariate_stan.stan"
# Read the data into R
N <- nrow(stan_data) # Number of observations
y_comparative <- round(stan_data$Count_COMP*stan_data$PercentMorph_COMP)+1 # Observed successes
y_superlative <- round(stan_data$Count_SUP*stan_data$PercentMorph_SUP)+1 # Observed successes
n_comparative <- stan_data$Count_COMP +1 # Total number of trials
n_superlative <- stan_data$Count_SUP +1 # Total number of trials
X <- stan_data %>%
select(mon_er,l_more,i_er,r_er,cluster_more,finalStress_more,initialCluster_est,sibfinal_most,rInitial_er,Log_Freq_HAL) %>%
as.matrix()
K <- ncol(stan_data %>% select(mon_er,l_more,i_er,r_er,cluster_more,finalStress_more,initialCluster_est,sibfinal_most,rInitial_er,Log_Freq_HAL)) # Number of predictors
FreqNonCentered_lemma <- log(stan_data$Count_COMP+stan_data$Count_SUP)
FreqNonCentered_superlative <- log(stan_data$Count_SUP)
FreqNonCentered_comparative <- log(stan_data$Count_COMP)
mu <- matrix(0, nrow = 2, ncol = N)
# a_ulam$Freq_noncentered
# Create a list of data for STAN
stan_data_list <- list(N = N,
y_comparative = y_comparative,
y_superlative = y_superlative,
n_comparative = n_comparative,
n_superlative = n_superlative,
X = X,
K = K,
FreqNonCentered_lemma = FreqNonCentered_lemma,
FreqNonCentered_comparative = FreqNonCentered_comparative,
FreqNonCentered_superlative = FreqNonCentered_superlative,
mu=mu
)
# Compile the STAN model
model <- cmdstan_model(stan_model_file)
# Fit the model to the data
fit <- model$sample(data = stan_data_list, cores = 4, chains = 4, threads = 2, init = .1)
I’m getting errors related to the covariance matrix not being positive definite and not symmetric, and I’m not quite sure how to proceed. I realize this is probably something fairly trivial for an experienced user, but I’m new to this aspect of stan and so not sure where to go from here. Any help or suggestions would be very much appreciated!