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.