This is a very straightforward question that has frustrated me over the last two days and maybe I am just missing an obvious error. I want to estimate polling errors for past US elections from polling data with errors being correlated across states. Elections are indexed by I, states by S. Xi is the correlated component and I followed the instructions in the manual and have done this before but I get the “Log probability evaluates to log(0), i.e. negative infinity.” issue. Any thoughts where the issue is?

data {
int N;
int I;
int S;
vector[S] nu;
int i[N];
int s[N];
int y[N];
int n[N];
matrix[S, I] outcome;
}
transformed data {
matrix[S, I] logit_outcome;
logit_outcome = logit(outcome);
}
parameters {
real alpha;
vector[S] beta;
vector<lower = 0>[S] sigma;
matrix[S, I] raw_xi;
cholesky_factor_corr[S] L_Omega;
}
transformed parameters {
matrix[S, I] xi;
xi = (diag_pre_multiply(sigma, L_Omega) * raw_xi);
}
model {
vector[N] mu;
L_Omega ~ lkj_corr_cholesky(1);
sigma ~ student_t(nu, 0, 0.2);
alpha ~ normal(0, 0.2);
beta ~ normal(0, 0.2);
to_vector(raw_xi) ~ std_normal();
for (jj in 1:N){
mu[jj] = logit_outcome[s[jj], i[jj]] + alpha + beta[s[jj]] + xi[s[jj], i[jj]];
}
y ~ binomial_logit(n, mu);
}

data {
int N;
int I;
int S;
vector[S] nu;
int i[N];
int s[N];
int y[N];
int n[N];
matrix[S, I] outcome;
}
transformed data {
matrix[S, I] logit_outcome;
logit_outcome = logit(outcome);
}
parameters {
real alpha;
vector[S] beta;
vector<lower = 0>[S] sigma;
matrix[S, I] raw_xi;
cholesky_factor_corr[S] L_Omega;
}
transformed parameters {
matrix[S, I] xi;
xi = (diag_pre_multiply(sigma, L_Omega) * raw_xi);
}
model {
vector[N] mu;
L_Omega ~ lkj_corr_cholesky(1);
sigma ~ student_t(nu, 0, 0.2);
alpha ~ normal(0, 0.2);
beta ~ normal(0, 0.2);
to_vector(raw_xi) ~ std_normal();
for (jj in 1:N){
mu[jj] = logit_outcome[s[jj], i[jj]] + alpha + beta[s[jj]];
}
y ~ binomial_logit(n, mu);
}

This works just fine

data {
int N;
int I;
int S;
vector[S] nu;
int i[N];
int s[N];
int y[N];
int n[N];
matrix[S, I] outcome;
}
transformed data {
matrix[S, I] logit_outcome;
logit_outcome = logit(outcome);
}
parameters {
real alpha;
vector[S] beta;
}
model {
vector[N] mu;
alpha ~ normal(0, 0.2);
beta ~ normal(0, 0.2);
for (jj in 1:N){
mu[jj] = logit_outcome[s[jj], i[jj]] + alpha + beta[s[jj]];
}
y ~ binomial_logit(n, mu);
}