Hello,

I am currently trying to fit a hierarchical variant of a specific type of network-based choice model (link for the interested, I am here only focusing on the second of the two sub-models: https://www.sociologicalscience.com/articles-v4-14-318/). The outcome is the choice of an interaction partner from a set of available target actors for a given relational event, where the sender in a given event can be considered fixed here. The number of available target actors can vary from 60 to 90 in my case, the total number of events is ca. 4500. The predictors consist of statistics regarding the history of past interactions linking the sender of a given interaction event with potential target actors (e.g. the number of past times the sender i chose a target j). The hierarchical model structure is over different types of interactions (in my case 3 different types of economic exchanges).

I have managed to build a working version based on a custom likelihood function, but it is relatively slow. I have tried to vectorize where possible but might have missed some opportunities. I also implemented a non-centered parametrization as suggested in the manual because I encountered divergend iterations, but the fitting process is actually quite a bit slower than with the centered parametrization, which seems odd. Gradient evaluation is in the ballpark of 0.15 to 0.2 seconds.

Here is the Stan model:

```
functions {
real choiceModel_lpmf(int receiver, int sender, vector beta, vector[] X, int A){
real lprob;
real denom = 0;
for (k in 1:A){
if (k == sender) continue;
denom += exp(dot_product(beta, X[k]));
}
lprob = dot_product(beta, X[receiver]) - log(denom);
return(lprob);
}
}
data {
int<lower=0> N_events; // number of events
int<lower=0> N_pred; // number of predictors / network statistics
int<lower=0> A[N_events]; // number actors present at each event
int<lower=1> N_L; // number of levels
int<lower=0,upper=N_L> L[N_events]; // levels for events
int<lower=1> receiver[N_events]; // receiver position at current event
int<lower=1> sender[N_events]; // sender position at current event
vector[N_pred] X[sum(A)]; // predictor matrix in unragged format
}
parameters {
vector[N_pred] mu;
vector<lower=0>[N_pred] sigma;
vector[N_pred] beta_raw[N_L];
}
transformed parameters {
vector[N_pred] beta[N_L];
for(l in 1:N_L)
beta[l] = mu + sigma .* beta_raw[l];
}
model {
int pos; // needed for segmenting the unragged array X
pos = 1;
// Priors
mu ~ normal(0,2);
sigma ~ exponential(1);
for(l in 1:N_L)
beta_raw[l] ~ std_normal();
// likelihood (evaluated only at the currently active actors)
for(n in 1:N_events){
receiver[n] ~ choiceModel(sender[n], beta[L[n]], segment(X, pos, A[n]), A[n]);
pos = pos + A[n];
}
}
generated quantities {
vector[N_events] log_lik;
int pos;
pos = 1;
for(n in 1:N_events){
log_lik[n] = choiceModel_lpmf(receiver[n] | sender[n], beta[L[n]], segment(X, pos, A[n]), A[n]);
pos = pos + A[n];
}
}
```

I would greatly appreciate any advice on how to improve performance or on more efficent coding / overall model structure, as I am relatively new to Stan and Bayesian Stats in general.

Thank you,

Jakob