model {
# remember that k = number of claims + 1
for (i in 1:N) {
y[i] ~ dpois(lambda)
y.rep[i] ~ dpois(lambda)
for (k in 1:K) {
eq[i,k] <- equals(y[i], k-1) # needed to construct
# aggregate data
eq.rep[i,k] <- equals(y.rep[i], k-1)
}
}
for (k in 1:K) {
m[k] <- sum(eq[,k]) # aggregate counts
m.rep[k] <- sum(eq.rep[,k])
# log of expected counts
logE[k] <- log(N) - lambda + (k-1)*log(lambda)
- logfact(k-1)
# likelihood ratio statistic
LR[k] <- 2*m[k]*(log(m[k]+0.00001) - logE[k])
LR.rep[k] <- 2*m.rep[k]*(log(m.rep[k]+0.00001)
- logE[k])
}
G <- sum(LR[])
G.rep <- sum(LR.rep[])
P <- step(G - G.rep)
lambda ~ dgamma(0.5, 0.001) # Jeffreys prior
# more powerful invariant test for excess zero counts
T.rep <- m.rep[1]*m.rep[3]/(m.rep[2]*m.rep[2])
P.inv <- step(T - T.rep)
}

Data:
list(N=35,K=10,T=4,
y=c(0,0,0,0,0,0,0,0,0,0,
1,1,1,1,1,2,2,2,2,2,
2,2,2,2,2,3,3,3,3,3,
3,3,4,4,5))

Inits:
list(lambda=2)
          node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   G   7.96   1.465   0.0137   6.956   7.382   12.1   10001   10000
   G.rep   6.473   3.531   0.03121   1.668   5.747   15.09   10001   10000
   P   0.701   0.4578   0.004491   0.0   1.0   1.0   10001   10000
   P.inv   0.9895   0.1019   0.001009   1.0   1.0   1.0   10001   10000
   T.rep   0.7173   1.116   0.01024   0.08304   0.4628   2.72   10001   10000
   lambda   1.7   0.221   0.001996   1.291   1.693   2.161   10001   10000