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