#Bayes 1 model with up to 4 joinpoints model{ for(i in 1:n){ Obs[i]~dpois(mean[i]) #Modelling of the expected cases for every period #the maximum and minimum bounds are imposed to avoid numerical overflow. #x1 and x2 are the initial and final years of study. log(mean[i])<-log(Pob[i])+alpha+beta0*(anyo[i]-(x1+x2)/2)+ min(max(delta[1]*beta[1]*BaseJoin[i,1]+ delta[2]*beta[2]*BaseJoin[i,2]+ delta[3]*beta[3]*BaseJoin[i,3]+ delta[4]*beta[4]*BaseJoin[i,4],-10),10) #estimated rate at every period rate[i]<-exp(alpha+beta0*(anyo[i]-(x1+x2)/2)+min(max(delta[1]*beta[1]*BaseJoin[i,1]+ delta[2]*beta[2]*BaseJoin[i,2]+delta[3]*beta[3]*BaseJoin[i,3]+delta[4]*beta[4]*BaseJoin[i,4],-10),10)) } #APC at every period for(i in 1:(n-1)){APC[i]<-100*(rate[i+1]-rate[i])/rate[i]} #Prior distributions for joinpoints' locations for(j in 1:4){tau[j]~dunif(tauinf,tausup)} tauinf<-x1+2 tausup<-x2-2 #This condition imposes an order restriction and a minimimum distance between joinpoints #one has to be given as data and necesarily has to be equal to 1. one~dbern(condicion) condicion<-step(tau[2]-tau[1]-2)*step(tau[3]-tau[2]-2)*step(tau[4]-tau[3]-2) #flat prior distributions for alpha and beta0 alpha~dflat() beta0~dflat() #Prior distribution for parameter beta beta[1:4]~dmnorm(mean.beta[1:4],prec.u[1:4,1:4]) for(i in 1:4){ mean.beta[i]<-0 for(k in 1:n){BaseJoin2[k,i]<-BaseJoin[k,i]*pow(w[k],1/2)} #prior precision matrix for beta. #The parameter g makes the prior distribution of beta to be a Cauchy distribution for(j in 1:4){ prec.u[i,j]<-g*inprod(BaseJoin2[,i],BaseJoin2[,j])*max(step(i-j)*step(j-i),delta[i]*delta[j]) } } #Prior distribution for g. g~dgamma(0.5,n2) n2<-n/2 #Parameter w for the prior precision of beta for(j in 1:n){ w[j]<-Pob[j]*exp(alpha+beta0*(anyo[j]-(x1+x2)/2)) } ######################### #Break-point definitions# ######################### for(j in 1:4){ for(i in 1:n){ BaseJoin[i,j]<-1+(anyo[i]-tau[j])*gamma[i,j] gamma[i,j]<-gamma1[j]+(gamma2[j]-gamma1[j])*step(anyo[i]-tau[j]) } gamma1[j]<-(d0*c2[j]-c0*d2[j])/(d2[j]*c1[j]-d1[j]*c2[j]) gamma2[j]<-(d1[j]*c0-d0*c1[j])/(d2[j]*c1[j]-d1[j]*c2[j]) c1[j]<-inprod(anyocent[,j],primerperiodo[,j]) c2[j]<-inprod(anyocent[,j],segundoperiodo[,j]) d1[j]<-inprod(anyocent2[,j],primerperiodo[,j]) d2[j]<-inprod(anyocent2[,j],segundoperiodo[,j]) for(i in 1:n){ anyocent[i,j]<-anyo[i]-tau[j] anyocent2[i,j]<-anyo[i]*anyocent[i,j] primerperiodo[i,j]<-step(-anyocent[i,j]) segundoperiodo[i,j]<-1-primerperiodo[i,j] } } c0<-n;d0<-sum(anyo[]) # Prior distributions for delta delta[1]~dbern(p1) delta[2]~dbern(p2) delta[3]~dbern(p3) delta[4]~dbern(p4) # Prior distributions for the probability of deltas being included in the model p1~dbeta(0.5,1.5) p2~dbeta(0.5,1.5) p3~dbeta(0.5,1.5) p4~dbeta(0.5,1.5) #Indicator variable pointing which of the models M_0,M_1,M_2,M_3 or M_4 is sampled ModelSampled<-sum(delta[]) }