I agree it is an indexing issue, but then if that is the issue the random intercept only model should work but doesn’t.
Here is the stan code, exact data and r code am using to check.
reduce_sum_intercetp.R (3.2 KB)
reduce_sum_logit_int.stan (1.2 KB)
sample_data.csv (67.8 KB)
This is the error I get:
Exact stan code I used:
functions {
real partial_sum(int[] y_slice,
int start, int end,
matrix x,
vector beta,
vector rbeta,
int[] ll) {
return bernoulli_logit_lpmf(y_slice[start:end] |(x[start:end,] * beta) + rbeta[ll[start:end]]);
}
}
data {
int<lower=0> N;//Number of observations
int<lower=1> K;//Number of predictors with non-random slope
int<lower=1> L;//Number of customers/groups
int<lower=0,upper=1> y[N];//Binary response variable
int<lower=1,upper=L> ll[N];//Number of observations in groups
matrix[N,K] x;
}
parameters {
real rbeta_mu; //mean of distribution of beta parameters
real<lower=0> rbeta_sigma; //variance of distribution of beta parameters
vector[L] beta_raw; //group-specific parameters beta
vector[K] beta;
}
transformed parameters {
vector[L] rbeta;
for (l in 1:L)
rbeta[l] = rbeta_mu + rbeta_sigma * beta_raw[l]; // coefficients on x
}
model {
rbeta_mu ~ normal(0,10);
rbeta_sigma ~ cauchy(0,5);
beta~normal(0,5);
for (l in 1:L)
beta_raw[l] ~ normal(0,1);
target += reduce_sum(partial_sum,y,1,x,beta,rbeta,ll);
}
and snap shot of model data: