I might be misunderstanding something, but let me clarify how my suggestion would apply to what I think youâ€™re describing. You have longitudinal count data, observed at two time points per case (N = number of cases; T = number of time points = 2). You have a model for the raw count data (N*2 observations), which you model as Poisson-distributed.

Y_1 \sim Poisson(\lambda_1) \\\
Y_2 \sim Poisson(\lambda_2)

You also have a model for the change (K) from time 1 to time 2 (N observations), which you model as Skellam-distributed.

(Y_1 - Y_2) = K \sim Skellam(\lambda_1, \lambda_2)

As you note, the log-likelihoods of the two models will not be comparable. The Poisson model has twice as many observations. But more importantly, the Skellam distribution is a collapsed version of the independent bivariate Poisson distributions. That is, the bivariate Poisson distribution assigns probabilities to each unique combination of Y_1 and Y_2. But the Skellam distribution assigns probabilities to *sets* of unique combinations that share the same difference. For example, the Skellam distribution assigns a probability to

PR(Y_1 - Y_2 = -1) = PR(Y_1 = 1, Y_2 = 2) + PR(Y_1 = 2, Y_2 = 3) + \text{ ...etc.}

whereas the Poisson distribution assigns a probability to

PR(Y_1 = 1, Y_2=2)

Since PR(Y_1 = 1, Y_2=2) is contained by PR(Y_1 - Y_2 = -1), the Skellam distribution will by definition generate a larger log-likelihood than the Poisson distribution, assuming the same data and parameter values.

In this case, the two models share the same parameters (\lambda_1 and \lambda_2), so you can easily calculate the log-likelihood of both models simultaneously. You can run the model twice, selecting one and then the other model for the purpose of inference, and the compare the appropriate log-likelihoods. See the rough, incomplete model code below as an example.

```
data{
int N;
int<lower = 0, upper = 1> P; // Use poisson model?
array[N,2] int<lower = 0> Y;
}
transformed data{
array[N] int K;
for(n in 1:N){
K[n] = Y[n,1] - Y[n,2];
}
}
parameters{
real<lower = 0> lambda1;
real<lower = 0> lambda2;
}
transformed parameters{
matrix[N,2] log_lik_poisson;
vector[N] log_lik_skellam;
for(n in 1:N){
log_lik_poisson[n,1] = poisson_lpmf(Y[n,1] | lambda1);
log_lik_poisson[n,2] = poisson_lpmf(Y[n,2] | lambda2);
log_lik_skellam[n] = skellam_lpmf(K[n] | lambda1, lambda2);
}
}
model{
if(P == 1){
target += log_lik_poisson;
} else{
target += log_lik_skellam
}
}
```