Help Time Series Model - Sequential Population Analysis

[Allan Debertin]

Hi,

I am attempting to code a Sequential Population Model using STAN, based on WinBUGS code from Millar & Meyer (2000) "Bayesian state-space modeling of age-structured data: fitting a model is just the beginning".

Here is the WinBUGS/JAGS code:

cat("model{

######Vague priors on first row and column#######

for(y in 1:(Y+1)) {LogRecruits[y]~dnorm(10,.01);

    LogN[1,y]<-LogRecruits[y]; }

    for(a in 2:(A+1)){LogInitials[a-1]~dnorm(10,.01);

    LogN[a,1]<-LogInitials[a-1]; }

    

    ########Distribution for data, I[a,y]######

    for (a in 1:A) {

    for (y in 1:Y) { MeanLogI[a,y]<-Logq[a]+(LogN[a,y]+10*LogN[a+1,y+1])/12;

    Prec[a,y]<-PrecLogI/a;

    I[a,y]~dlnorm(MeanLogI[a,y],Prec[a,y]); } }

    

    #####Prior on mortalities, hierarchical######NB: log(0.2)=-1.61

    for(a in 1:A){ MeanLogAgeM[a]~dnorm(-1.21,1);}

    for(y in 1:Y){ MeanLogDeltaM[y]~dnorm(0,1);}

    PrecLogM~dgamma(1,1);

    for(a in 1:A) {

    for(y in 1:Y) { MeanLogM[a,y]<-MeanLogAgeM[a]+MeanLogDeltaM[y];

    M[a,y]~dlnorm(MeanLogM[a,y],PrecLogM); } }

    

    #####Cohort equation#####

    for(a in 2:(A+1)) {

    for(y in 2:(Y+1)) {

    LogN[a,y]<-log(max(exp(LogN[a-1,y-1])*exp(-M[a-1,y-1])-C[a-1,y-1]*exp(-0.5*M[a-1,y-1]),10)); } }

    

    

    ######Prior on q's######

    MeanLogq<--10; PrecLogq<-1;

    for (a in 1:A){ Logq[a]~dnorm(MeanLogq,PrecLogq); }

    

    ######Prior on q's, quadratic######

    #gamma~dnorm(-8,1); 

    #alpha~dunif(0.0001,0.02); 

    #beta~dnorm(5,0.25);

    #for (a in 1:A){ MeanLogq[a]<-gamma-alpha*(a-beta)*(a-beta);

    #                Logq[a]~dnorm(MeanLogq[a],16); }

    

    ######Prior on precision of log(I)######

    PrecLogI~dgamma(10,1);

    

    ######Vector of 1983 numbers for monitoring#####

    for(a in 1:A) {N2001s[a]<-exp(LogN[a,6])/1000; }

    }", 

    

    file="agestructure.jag")

Below is the code that works in STAN, but isn't exactly the same, see highlighted part.: 

FYI:  A = 8, Y = 25. for the generated quantities. I could change the size to read real N_log[A+1,N+1]

data {

    int<lower=0> A;

    int<lower=0> Y;

    matrix[A,Y] I;

    matrix[A,Y] C;

    }

    parameters {

    real<lower=-50,upper=50> LogN[A,Y];

    real<lower=-50, upper=50> Logq[A];

    real<lower=0, upper=2> M[A,Y];

    real<lower=-50, upper=50> MeanLogAgeM[A];

    real<lower=-50, upper=50> MeanLogDeltaM[Y]; 

    real<lower=0, upper=10>  PrecLogM;

    real<lower=0, upper=50> PrecLogI;

    real<lower=0, upper=30> LogInitials[9];

    real<lower=0, upper=30> LogRecruits[26];

    }

    model {

   

    PrecLogI ~ gamma(10,1);

       for(a in 1:A){ MeanLogAgeM[a]~normal(-1.21,1);}

       for(y in 1:Y){ MeanLogDeltaM[y]~normal(0,1);}

        PrecLogM ~ gamma(1,1);

       for(a in 2:(9)){LogInitials[a-1]~normal(10,.01);}

       for(y in 1:(26)){LogRecruits[y]~normal(10,.01);}

    //prior on q's

        {real MeanLogq;

         real PrecLogq;     

          MeanLogq <- -10;

          PrecLogq <- 1;

          for (a in 1:A) {

          Logq[a] ~ normal(MeanLogq, PrecLogq);

          }

         }

    // state equations

    // **************

     //mortality equation  

     {real MeanLogM[A,Y];

 

       for(a in 1:A) {

        for(y in 1:Y) { 

          MeanLogM[a,y] <- MeanLogAgeM[a] + MeanLogDeltaM[y];

           M[a,y] ~ lognormal(MeanLogM[a,y],PrecLogM); 

        } 

       }

      }

         

    // observation equation

    // ********************

    //Distribution for data, I[a,y]

      { real MeanLogI[A,Y];

        real Prec[A,Y];

         for (a in 1:A) {

              for (y in 1:Y) { 

              MeanLogI[a,y]<- Logq[a]+ (LogN[a,y]);

                     Prec[a,y]<-PrecLogI/a;

                     I[a,y] ~ lognormal(MeanLogI[a,y],Prec[a,y]); 

        } 

      }

    }

    }

generated quantities {

  

  //cohort equation

  real N_log[9,26];

  for(y in 1:(26)) {

    N_log[1,y]<-LogRecruits[y];

  }

  for(a in 2:(9)){

    N_log[a,1]<-LogInitials[a];

  }

  for(a in 2:(9)) {

    for(y in 2:(26)) {

      N_log[a,y]<-log(fmax(exp(LogN[a-1,y-1])*exp(-M[a-1,y-1])-C[a-1,y-1]*exp(-0.5*M[a-1,y-1]),10));

    } 

  }

}

 I want to change the highlighted line to read below, since the survey data I am working with also occurs late in the year and the index likely has to do more with the next year in the time step.

  ########Distribution for data, I[a,y]######

  for (a in 1:A) {

    for (y in 1:Y) { MeanLogI[a,y]<-Logq[a]+(LogN[a,y]+11*LogN[a+1,y+1])/12;

                     Prec[a,y]<-PrecLogI/a;

                     I[a,y] ~ lognormal(MeanLogI[a,y],Prec[a,y]); } }

With the change in the code that I highlighted, I get the error message that:

INDEX OPERATOR [] OUT OF BOUNDS; index=26; lower bound=1; upper bound=25; index position=2; LogN

Which makes sense, because A = 8 and Y=25, but I want to include those following LogN estimates into the analysis.

I think I need to either need to code for "Missing data" or using slicing commands such as the "block" function, although I haven't found any examples using the "block" function.

Any help would be great.

Thanks,

Allan Debertin

[Bob Carpenter]

I'm not sure what you're asking. No matter what you do, you won't be able
to index beyond the size of array you have. The LogN[a+1, y+1] expression runs
over the bounds of LogN.

What data are you thinking of as missing? Or in terms that can
be translated straight into Stan, what density are you after for
the a+1 and y+1 case?

> Any help would be great.

I'd strongly suggest picking a formatting style for your code
and sticking to it. There's a suggestion in the back of the manual.

The second index doesn't need to be parenthesized --- "for (a in 2:9)" works.

You can make everything much more efficient by using vector
forms of the functions. You can also use local variables in Stan
rather than having to make every local into an array as in BUGS. So this

> { real MeanLogI[A,Y];
> real Prec[A,Y];
>
> for (a in 1:A) {
> for (y in 1:Y) {
> MeanLogI[a,y]<- Logq[a]+ (LogN[a,y]);
> Prec[a,y]<-PrecLogI/a;
> I[a,y] ~ lognormal(MeanLogI[a,y],Prec[a,y]);
> }

can be written more transparentlya nd efficiently using local variables as:

for (a in 1:A) {
for (y in 1:Y) {

real meanlog;
real prec;
prec <- PrecLogI / a;
meanlog <- Logq[a] + LogN[a,y];

I[a,y] ~ lognormal(MeanLogI[a,y], Prec[a,y]);
}
}

but its' even better if you just declare LogN as an array of
vectors (less efficiently as a matrix) and then write the entire
loop as:

vector<lower=-50,upper=50>[Y] LogN[A];

for (a in 1:A)

I[a] ~ lognormal(Logq[a] + LogN[a], PrecLogI / a);

We don't recommend the kind of truncated priors you're using, though
there are still some left in our manual. We have recommendations in
the regression chapter (and more in the problematic posteriors chapter
for special cases).

If you do use those truncated priors, just make sure your estimates aren't
butting up against a boundary.

- Bob

[Allan Debertin]

Hi Bob,

Thanks for the quick reply and advice on coding. I'll use the appropriate format in the future. I will check that the posterior distributions are not near the boundary.

I will try to rephrase my question in the coding. This is the part of the coding I would like to change.

# WinBUGS Code

 

  ########Distribution for data, I[a,y]######

   for (a in 1:A) {

   for (y in 1:Y) { MeanLogI[a,y]<-Logq[a]+(LogN[a,y]+10*LogN[a+1,y+1])/12;

   Prec[a,y]<-PrecLogI/a;

   I[a,y]~dlnorm(MeanLogI[a,y],Prec[a,y]); } }

   
 

//STAN CODE (my first attempt which is wrong)

   for (a in 1:A) {

   for (y in 1:Y) { MeanLogI[a,y]<-Logq[a]+(LogN[a,y]+10*LogN[a+1,y+1])/12;

                   Prec[a,y]<-PrecLogI/a;

                   I[a,y] ~ lognormal(MeanLogI[a,y],Prec[a,y]); } }

The expression runs in WinBUGS, because the order of the code does not matter, so I think it ignores LogN[a+1,y+1] (or uses NAs) for the last step. However, as you know, the order of code does matter in STAN.

To calculate  MeanLogI[a,y] I would to use Logq[a]+(LogN[a,y]+11*LogN[a+1,y+1])/12;, until the last available parameter point in the specified LogN matrix. 

For the last data point, I would only use Logq[a]+LogN[a,y].

I tried to code it this way(see below) to tackle it, but had  a different problems. I am sure there is a better approach using the "tail", "head" or "block" functions, but I'm not sure how to construct them given that there isn't many examples in the  STAN Manual.

    //Distribution for data, I[a,y]

  {

  real MeanLogI[A,Y];

  real Prec[A,Y];

    for (a in 1:A-1) {

      for (y in 1:Y-1) { 

        MeanLogI[a,y]<- Logq[a]+ (LogN[a,y]+11*LogN[a+1,y+1])/12;

      }

    }

    MeanLogI[8,25] <- Logq[8]+ LogN[8,25];

    for (a in 1:A) {

      for (y in 1:Y) { 

        Prec[a,y]<-PrecLogI/a;

        I[a,y] ~ lognormal(MeanLogI[a,y],Prec[a,y]); 

      }

    }

  }

## Error "Rejecting proposed initial value with zero density."

There are similar (more complex) models I would like to run in STAN that I have similar coding issues in terms of the last step of a time series, So if we can solve this, I can move on to those models. 

Many Thanks,

Allan

  

[Daniel Lee]

Hi Allan,

I'm still not sure what you're trying to do, so my comments might be off base.

Suggestions:

- debug by print()? Stan has a print() function that could be used to print out values so you're sure the data manipulation is done properly.

- vectorize? To me, it looks cleaner, but I can see where coming from BUGS/JAGS it may not (yet). For example, line 24, change: "for(a in 1:A){ MeanLogAgeM[a]~normal(-1.21,1);}" to "MeanLogAgeM ~ normal(-1.21, 1);"

- If you're still trying to get the model running, I would build up. To me, that means starting with just the priors and making sure that's ok. Then moving to the state equations, then finally the observation equations so that you can tell where things are going wrong.

For the observation equation part, the part that you've highlighted:

- what are MeanLogI[A, 1], ... MeanLogI[A, Y], and MeanLogI[1, Y], ..., MeanLogI[A, Y]? They don't seem to be set. Perhaps this is where the error is coming from.
- change: "I[a,y] ~ lognormal(MeanLogI[a,y],Prec[a,y]);" to "I[a, y] ~ lognormal(MeanLogI[a, y], PrecLogI / a);" and get rid of Prec.

For the parameter declarations:

- are all the bounds necessary? I don't see them in the BUGS/JAGS model and this changes the model by truncating the distributions to those bounds.

Daniel

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Allan Debertin]

Hi Daniel,

I'll try your suggestions and Bob's with the code I know works and get back to you. 

Thanks for letting me know the nuances more. Stay tuned.

[Bob Carpenter]

All you need to do is figure out what you want the model to do.
Once you know the density, we can show you how to code it in Stan.

- Bob

[Allan Debertin]

Hi Guys,

Using your advice I got the model to work. I also forgot to change Normal and Log-Normal distributions since they are different specifications between JAGS and STAN.

This model is running slow, so there is probably a better way to code it. See highlighted part which causes it to run slowly. If I use this:

//Distribution for data, I[a,y]

  for (a in 1:A) {     

     I[a] ~ lognormal(Logq[a]+(LogN[a],PrecLogI/a); 

   } 

 it runs much faster, but I want the auto-correlation between years to help estimate the index (except for the last year).

library(rstan)

sink("agestructure.stan")

cat("

  data {

    int<lower=0> A;

    int<lower=0> Y;

    real I[A,Y];

    real C[A,Y];

    }

  parameters {

    vector[Y] LogN[A];

    vector[A] Logq;

    vector<lower=0, upper=10>[Y] M[A];

    vector[A] MeanLogAgeM;

    vector[Y] MeanLogDeltaM; 

    real PrecLogM;

    real PrecLogI;

    vector[A] LogInitials;

    vector[Y+1] LogRecruits;

    }

  model {

    PrecLogI ~ gamma(10,1);

    PrecLogM ~ gamma(1,1);

    MeanLogAgeM ~ normal(-1.61,1);

    MeanLogDeltaM ~ normal(0,1);

    LogRecruits ~ normal(10,10);

    Logq ~ normal(-10, 1);

    

    for (a in 2:A+1){

      LogInitials[a-1] ~ normal(10,10);

    }

    //Distribution for data, I[a,y]

    {real MeanLogI[A,Y];

     MeanLogI[A,Y]<-Logq[8]+LogN[8,25];

     I[8,25] ~ lognormal(MeanLogI[8,25],PrecLogI/8); 

     for (a in 1:A-1) {

      for (y in 1:Y-1) {

        MeanLogI[a,y]<-Logq[a]+(LogN[a,y]+10*LogN[a+1,y+1])/12;

        I[a,y] ~ lognormal(MeanLogI[a,y],PrecLogI/a); 

      } 

     }

    }

    //mortality equation

    for(a in 1:A) {

        M[a] ~ lognormal(MeanLogAgeM[a]+MeanLogDeltaM,PrecLogM); 

    }

   

    }

generated quantities {

  

  //cohort equation

  real N_log[A+1,Y+1];

  for(y in 1:Y+1) {

    N_log[1,y]<-LogRecruits[y];

  }

  for(a in 2:A+1){

    N_log[a,1]<-LogInitials[a-1];

  }

  for(a in 2:A+1) {

    for(y in 2:Y+1) {

      N_log[a,y]<-log(fmax(exp(LogN[a-1,y-1])*exp(-M[a-1,y-1])-C[a-1,y-1]*exp(-0.5*M[a-1,y-1]),10));

    } 

  }

}

    ", fill=TRUE) 

sink()

Is there a better way to code the highlighted part?

Cheers,

Allan

[Daniel Lee]

Hi Allan,

What do you mean slow? Are the parameters mixing?

Are you getting reasonable inferences? There really isn't a point trying to make it fast if it's not correct.

Daniel

On Dec 16, 2014 3:52 PM, "Allan Debertin" <alland...@gmail.com> wrote:

Hi Guys,

Using your advice I got the model to work. I also forgot to change Normal and Log-Normal distributions since they are different specifications between JAGS and STAN.

This model is running slow, so there is probably a better way to code it. See highlighted part which causes it to run slowly. If I use this:

//Distribution for data, I[a,y]

{

  for (a in 1:A) {     

     MeanLogI[a,y]<-Logq[a]+(LogN[a,y]+10*LogN[a+1,y+1])/12;

     I[a] ~ lognormal(Logq[a]+(LogN[a],PrecLogI/a); 

   } 

 }

}

 it runs much faster, but I want the auto-correlation between years to help estimate the index (except for the last year).

library(rstan)

sink("agestructure.stan")

cat("

  data {

    int<lower=0> A;

    int<lower=0> Y;

    real I[A,Y];

    real C[A,Y];

    }

  parameters {

    vector[Y] LogN[A];

    vector[A] Logq;

    vector<lower=0, upper=10>[Y] M[A];

    vector[A] MeanLogAgeM;

    vector[Y] MeanLogDeltaM; 

    real PrecLogM;

    real PrecLogI;

    vector[A] LogInitials;

    vector[Y+1] LogRecruits;

    }

  model {

    PrecLogI ~ gamma(10,1);

    PrecLogM ~ gamma(1,1);

    MeanLogAgeM ~ normal(-1.61,1);

    MeanLogDeltaM ~ normal(0,1);

    LogRecruits ~ normal(10,10);

    Logq ~ normal(-10, 1);

    

    for (a in 2:A+1){

      LogInitials[a-1] ~ normal(10,10);

    }

    //Distribution for data, I[a,y]

    {real MeanLogI[A,Y];

     MeanLogI[A,Y]<-Logq[8]+LogN[8,25];

     I[8,25] ~ lognormal(MeanLogI[8,25],PrecLogI/8); 

     for (a in 1:A-1) {

      for (y in 1:Y-1) {

        MeanLogI[a,y]<-Logq[a]+(LogN[a,y]+10*LogN[a+1,y+1])/12;

        I[a,y] ~ lognormal(MeanLogI[a,y],PrecLogI/a); 

      } 

     }

    }

    //mortality equation

    for(a in 1:A) {

        M[a] ~ lognormal(MeanLogAgeM[a]+MeanLogDeltaM,PrecLogM); 

    }

   

    }

generated quantities {

  

  //cohort equation

  real N_log[A+1,Y+1];

  for(y in 1:Y+1) {

    N_log[1,y]<-LogRecruits[y];

  }

  for(a in 2:A+1){

    N_log[a,1]<-LogInitials[a-1];

  }

  for(a in 2:A+1) {

    for(y in 2:Y+1) {

      N_log[a,y]<-log(fmax(exp(LogN[a-1,y-1])*exp(-M[a-1,y-1])-C[a-1,y-1]*exp(-0.5*M[a-1,y-1]),10));

    } 

  }

}

    ", fill=TRUE) 

sink()

Is there a better way to code the highlighted part?

Cheers,

Allan

[Allan Debertin]

Nevermind, something still not quite right. I'll work on it some more.

If I replace the highlighted bit with   

for (a in 1:A) {     

     I[a] ~ lognormal(Logq[a]+(LogN[a],PrecLogI/a); 

   } 

It compiles and inferences are reasonable per se. 

However, if you look at the WinBUGS model posted with the first message. 

I basically want to know how to translate this into STAN.

    ########Distribution for data, I[a,y]######

    for (a in 1:A) {

    for (y in 1:Y) { MeanLogI[a,y]<-Logq[a]+(LogN[a,y]+10*LogN[a+1,y+1])/12;

    Prec[a,y]<-PrecLogI/a;

    I[a,y]~dlnorm(MeanLogI[a,y],Prec[a,y]); } }

The index (survey) is closer in time with LogN[Y+1] than LogN[Y], and is likely correlated between the two steps. I'm working with fish, so age class LogN[a] grows into Log[A+1] through that time step. 

[Bob Carpenter]

First translate it to a density, then translate it to Stan.
It doesn't make sense as a density in Stan as written because
we don't know what's intended by a statement including LogN[a+1,y+1]
when a+1 and y+1 are outside the bounds declared for LogN.

So that probably means first understanding what BUGS is doing
with your model. Then it should be straightforward to translate.

- Bob

[Luc Coffeng]

Just in case it hasn't been pointed out: Stan's normal and log-normal distributions use standard deviations (sd), and not precision (1/sd^2). I noticed that you are using (presumed) conjugate gamma priors for precision parameters, which in fact are interpreted as standard deviations in Stan. Although this may not be the cause of the current issues you are running into, it will cause results to differ from an analysis in BUGS/JAGS.

Luc

[Allan Debertin]

Thanks Luc, I'll look into it.