#model
model
{
for (s in 1:N){
        for(t in 1:M){
                        y[s,t]~dpois(mu[s,t])
                        mu[s,t] <-e[s,t]*rr[s,t]
                        log(rr[s,t])<-a0+u[s]+v[s]+g[t]+psi[s,t]+j[t]+beta1*fs[s,t]+beta2*js[s,t]+beta3*qw[s,t]+beta4*qy[s,t]+beta5*rz[s,t]+beta6*qd[s,t]+beta7*ab[s,t]+beta8*ac[s,t]+beta9*ad[s,t]
}
        psi[s,1]~dnorm(0,taupsi)
        for(t in 2:M){
        psi[s,t]~dnorm(psi[s,t-1],taupsi)
}
        v[s]~dnorm(0,tauv)
}
        # CAR prior distribution for random effects:
        u[1:N]~car.normal(adj[],weights[],num[],tauu)
                for(x in 1:sumNumNeigh){
                        weights[x]<-1}
}
        g[t]~dnorm(0,taug)
}
                        j[1:M] ~ car.normal(adj[], weights[], num[], tauj)   
                for(t in 1:1) {               
                weights[t] <- 1;                        adj[t] <- t+1;                        num[t] <- 1
        }
        for(t in 2:(M-1)) {
                weights[2+(t-2)*2] <- 1;                        adj[2+(t-2)*2] <- t-1
                weights[3+(t-2)*2] <- 1;                        adj[3+(t-2)*2] <- t+1;                        num[t] <- 2
        }
        for(t in M:M) {
                weights[(T-2)*2 + 2] <- 1;                        adj[(T-2)*2 + 2] <- t-1;                        num[t] <- 1
}
taupsi<-pow(sdpsi,2)
taupsi~dgamma(0.5,0.0005)
tauu<-1/pow(au,2)
au<-pow(sdu,2)
au~dgamma(0.5,0.0005)
tauu~dgamma(0.5,0.0005)
tauv<-1/pow(muv,1)
muv<-pow(sdv,2)
muv~dgamma(0.5,0.0005)
taug<-1/pow(mug,1)
mug<-pow(sdg,2)
mug~dgamma(0.5,0.0005)
tau.err  ~ dgamma(0.01, 0.01)                                # measurement error precision
sigma.err <- 1 / sqrt(tau.err)
sigma2.err <- 1/tau.err
tauj  ~ dgamma(0.01, 0.01)                                # random walk precision
sigma <- 1 / sqrt(tau)
sigma2 <- 1/tau
for(t in 1:T) {    day[t] <- t   }                                # include this variable to use in time series (model fit) plot
a0~dflat()
a1~dnorm(0,0.0001)
beta1~dnorm(0.0,tau.beta1)
tau.beta1<-1/pow(sdbeta1,2)
sdbeta1~dunif(0,10)
beta2~dnorm(0.0,tau.beta2)
tau.beta2<-1/pow(sdbeta2,2)
sdbeta2~dunif(0,10)
beta3~dnorm(0.0,tau.beta3)
tau.beta3<-1/pow(sdbeta3,2)
sdbeta3~dunif(0,10)
beta4~dnorm(0.0,tau.beta4)
tau.beta4<-1/pow(sdbeta4,2)
sdbeta4~dunif(0,10)
beta5~dnorm(0.0,tau.beta5)
tau.beta5<-1/pow(sdbeta5,2)
sdbeta5~dunif(0,10)
beta6~dnorm(0.0,tau.beta6)
tau.beta6<-1/pow(sdbeta6,2)
sdbeta6~dunif(0,10)
beta7~dnorm(0.0,tau.beta7)
tau.beta7<-1/pow(sdbeta7,2)
sdbeta7~dunif(0,10)
beta8~dnorm(0.0,tau.beta8)
tau.beta8<-1/pow(sdbeta8,2)
sdbeta8~dunif(0,10)
beta9~dnorm(0.0,tau.beta9)
tau.beta9<-1/pow(sdbeta9,2)
sdbeta9~dunif(0,10)
#data
list(N=3, M=4, e=structure(.Data=c(1,2,3,1,2,3,1,2,3,1,2,3), .Dim=c(3,4)), y=structure( .Date=c(1,2,3,4,5,6,7,8,9,10,11,12), .Dim=c(3,4)), fs=structure(.Data=c(1,2,3,1,2,3,1,2,3,1,2,3), .Dim=c(3,4)), js=structure(.Data=c(1,2,3,1,2,3,1,2,3,1,2,3), .Dim = c(3,4)), qw=structure(.Data=c(1,2,3,1,2,3,1,2,3,1,2,3), .Dim = c(3,4)), qy=structure(.Data=c(1,2,3,1,2,3,1,2,3,1,2,3), .Dim = c(3,4)), rz=structure(.Data=c(1,2,3,1,2,3,1,2,3,1,2,3), .Dim = c(3,4)), qd=structure(.Data=c(1,2,3,1,2,3,1,2,3,1,2,3), .Dim = c(3,4)), ab=structure(.Data=c(1,2,3,1,2,3,1,2,3,1,2,3), .Dim = c(3,4)), ac=structure(.Data=c(1,2,3,1,2,3,1,2,3,1,2,3), .Dim = c(3,4)), ad=structure(.Data=c(1,2,3,1,2,3,1,2,3,1,2,3), .Dim = c(3,4)), num=c(3,5,3,5,8,5,5,8,5,3,5,3),
adj=c(5,4,2,
6,5,4,3,1,
6,5,2,
8,7,5,2,1,
9,8,7,6,4,3,2,1,
9,8,5,3,2,
11,10,8,5,4
12,11,10,9,7,6,5,4,
12,11,8,6,5,
11,8,7,
12,10,9,8,7,
11,9,8),
sumNumNeigh = 54))
不知道为什么发出来会有很多横向，我再发一次