PROC MCMC help

[Jackie418]

I am trying to program a  poisson random effects model and having trouble.  I have seen the pump example but I have more categories of random effects  so looking for another example.

My model is representing counts as a unction of time (years) and random effects of site and year within a larger stratum:

log(lambda ijt) = Si + wj + gammait + betai (t-t0) + error

S is an intercept term for each stratum i and beta is a slope term for each stratum i, fixed effects.

w is a random effect which represents data from different sites within a stratum over time

gamma is an individual year effect within a straum (data are collected once a year)

Has anyone done anything like this in MCMC?

Thanks

[Funda_SAS]

I am trying to understand your model i.e. what is the diffence between wj and error -- do you have multiple observations coming from each site? Can you please share how your data looks like?

[Jackie418]

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;

[SteveDenham]

The code looks like it will do what you want.  Here is how I would go about adding in the year effect.  Subset the data to something like 4 or 5 years, and only one stratum, using again only 4 or 5 circles.  I would devise an appropriate hypothesis that allowed me to parametrize year as a single continuous variable--a linear trend for now.  Rewrite the w= equation to capture this linear trend.  For a first try, assume that the linear trend is identical over all circles, then step up to looking at different slopes per circle, adding enough complexity such that you still get credible results and address your scientific questions.  When you have a stable model that addresses what you want, increase the data.  I would add more strata first, then finally add more years, until you reach the capacity of your computer.  I say that because I think that will occur before you have all of the data included.

Steve Denham

[Jackie418]

Well I think my linear trend is captured in the term: c1[circle]*t

t is number of years since the start year.

I think this will give me the beta (slope) for each stratum that I a ultimately looking for.

But perhaps you mean that I should keep the intercept random by circle but not the slope?

In any case, I think the authors include a year effect so that one could make inferences about populations and trends within subsets of years.  I am still not sure how to do that. Maybe I just add another random intercept parameter...

Thanks again for your interest and help.

[SteveDenham]

I was wondering about "t" in that equation.  Given that it indexes year, you have a random slope and intercept model in hand.  I think that gives you what you want to compare to the published values.  If not, then I am missing something critical, and I hope I am not leading you down a path of futile work.

Steve Denham

[SteveDenham]

One thing to think very carefully about is that for almost all non-Gaussian distributions (for example, the Poisson) the specification of the model will NOT include a separable error term (see Chapter 1 in Stroup's Generalized Linear Mixed Models), unless you are modeling a specific overdispersion.  This may affect how you think about approaching this problem.

Steve Denham