Sorry for the delay in reply. The data are Christmas bird counts over a period of 40 years (n=2000). The ultimate objective is to make inferences about population trends. I am attempting to replicate a published analysis so that is where the model comes from. The strata are 3 adjacent regions with some similar bird habitat. Within each stratum parties count birds in predefined circles of specific size. There are anywhere from 4 to 20 circles in a stratum in any given year. Some circles become inactive and some new circles enter in so there is a lot of missing data for the circles over time. Counts are overdispersed Poisson so that is why I started with the pump example. I've been thinking about this some more and I think each circle really has it's own intercept and slope (not just intercept) and so are random observations of the fixed stratum intercept and slope. That has changed my approach. There is another layer of complexity I left out. The effort at each circle in any given year varies. This relationship between effort and expected count has been supposedly worked out so you will see that additional term in the code below. I wrote this to test just modeling one of the strata. I haven't run it yet because it will take a long time on my small machine. But this will not give me a year effect. Not sure how to add that in. Any comments are welcome. Thanks. proc mcmc data=gull ntu=1000 thin=10 nmc=10000 seed=7893 nbi=4000 propcov=quanew diag=(mcse ess) outpost=postout monitor=(alpha beta H p s2 s30 s31) DIC; array c0[83]; array c1[83]; array llambda[2001]; parms alpha 0 beta 0 H 1 p 1; parms c0: 0; parms c1: 0; parms s2 1 ; parms s30 1 s31 1; parms llambda: 1; beginnodata; sd = sqrt(s2); endnodata; w =c0[circle]+(c1[circle]*t)+(H*(effort**p - 1)/p); if state_yr eq 1 then lp = lpdfnorm(llambda[state_yr], w, sd); else lp = lp + lpdfnorm(llambda[state_yr], w, sd); prior alpha: ~ normal(0, var = 10000); prior beta: ~ normal(0, var = 10000); prior H ~ normal(0, var = 10000); prior p ~ uniform(-4,4); prior c0 : ~ normal(alpha, var = s30); prior c1 : ~ normal(beta, var = s31); prior s2 ~ general(-log(s2)); prior llambda: ~ general(lp); hyperprior s30 ~ igamma(0.001, scale = 1000); hyperprior s31 ~ igamma(0.001, scale = 1000); lambda = exp(llambda[state_yr]); model count ~ poisson(lambda); run;
... View more