Hello!

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!

Lucas

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.