Autoregressive model for distribution?


I have a set of vegetation plots organized by treatments and blocks. Within each plots, we harvested some individuals of each species and measured their height different years after installation. I am interested in modelling how distribution of height vary in time for each species.

I have already made a hierarchical joint model of mean and dispersion with the general form :

  for(n in 1:N){
   mu[n] = exp( beta_0 + beta_m[m[n]] + beta_b[b[n]] + beta_s[s[n]] + 
    beta_1s[s[n]] * t[n] );
   phi[n] = exp( gamma_0 + gamma_m[m[n]] + gamma_b[b[n]] + gamma_s[s[n]] + 
   gamma_1s[s[n]] * t[n] );

shape = mu .* mu ./ phi;
rate = mu ./ phi;
y ~ gamma(shape, rate); 

But because we are dealing with time, I would like to implement an autoregressive model describing how species s’ height distribution in treatment m and block b at time t depend upon the same distribution at time t-1.

I was thinking about first estimating mu and phi for each case, and then model their dependencies, but I have hard time finding how to code that.

Does someone have an idea??

Thank you very much!

I thought about something like this, ind being the index describing the combination of factor (treatment + block + species), and t being the vector of time describing when each observation have been measured :

  for(n in 1:N){
    shape[ind[n],t[n]] = mu[ind[n], t[n]] .* mu[ind[n], t[n]] ./ phi[ind[n], t[n]]
    rate[ind[n],t[n]] = mu[ind[n],t[n]] ./ phi[ind[n],t[n]];

    y[n] ~ gamma{shape[ind[n],t[n]], rate[ind[n],t[n]]);
for(ti in 1:T){
mu[,ti] ~ gamma(alpha + beta*mu[,ti-1], phi_1);
phi[,ti] ~ gamma(gamma_0 + gamma*mu[,ti-1], phi_2);

How many unique time points do you have?

I have 7 species in 5 treatements, repeated twice (2 blocks).

So I have 7 * 5 * 2 = 70 combinations measured during four years = 280 time points.

You can put together time series models like this just about any way you want that makes sense for your problem.

Whenever I see something like this I get worried:

for(n in 1:N){
    shape[ind[n],t[n]] = mu[ind[n], t[n]] .* ...

because you want to define each shape[i, j] exactly once and it’s not clear that things are all getting defined or aren’t being mulitply defined here.

The times-series parameterizations can be tricky because of the correlations among the elements of mu and phi in the posterior.

You can save some computation computing rate first and re-using:

rate = mu ./ phi;
shape = mu .* rate;

With multi-indexing and elementwise division you can do this all at once (if the values are all matrice):

rate[ind, t] = mu[ind, t] ./ phi[ind, t];

but see the point above to make sure this is really the right thing to do in your problem.