"Consider below the number of actuarial claims data for three groups of insurance policyholders,
year: 1 2 3 4 5
Grp1: 9 7 6 13 12
Grp2: 6 4 2 8 10
Grp3: 8 8 3 4 9
Run R and Jags to apply the following hierarchical model to analyze the data:
Yij ∼ Poisson(λij )
λij = Pijθj
θij ∼ Ga(α, β) Pij ∼ Ga(γ, δ)
α ∼ Ga(5, 5) γ ∼ U(0, 100)
β ∼ Ga(25, 1) δ ∼ U(0, 100),
where i = 1, 2, 3 and j = 1, . . . , 5.
What’s your conclusion about the group effect and the year effect?"
I have my model specifications drafted to be pulled into R using JAGS. My question is, how do I code in R to test for the effect of Group and the effect of Year separately? I've only ever used jags for one variable.
Here is my cookie-cutter JAGS code:
library(rjags)
forJags<-list( )
inits<-list( )
foo<jags.model(file="m2n4.bug",data = forJags,inits=inits)
out<-coda.samples(model=foo, variable.names = c( ), n.iter=50000,thin=5)
summary(out)
Here is my model:
model
{
for (i in 1:3,j in 1:5){
Y[i,j] ~ dpois(lambda[i,j])
lambda[i,j] = P[i,j]*theta[i,j]
theta[i,j] ~dgamma(alpha,beta)
P[i,j] ~ dgamma(gamma,delta)
}
alpha ~ dgamma(5,5)
beta ~ dgamma(25,1)
gamma ~ dunif(0,100)
delta ~ dunif(0,100)
}
Any input informing me of how to code such that I test for the effects separately would be huge.
Define your model as:
model
{
for (j in 1:5){P[j] ~ dgamma(gamma,delta)}
for (i in 1:3){
for(j in 1:5){
Y[i,j] ~ dpois(lambda[i,j])
lambda[i,j] = P[j]*theta[i,j]
theta[i,j] ~ dgamma(alpha,beta)
}
}
alpha ~ dgamma(5,5)
beta ~ dgamma(25,1)
gamma ~ dunif(0,100)
delta ~ dunif(0,100)
}
and then run:
library(rjags)
Y<-rbind(c(9, 7, 6, 13, 12),c( 6 ,4 ,2 ,8 ,10),c(8 ,8, 3, 4, 9))
forJags<-list('Y' = Y)
foo<-jags.model(file="m2n4.bug",data = forJags)
out<-coda.samples(model=foo, variable.names = c("theta","P"),n.iter=50000,thin=5)
summary(out)
You will then be able to see the separate effect for P (the year) and the single effect (theta)