Dear All,

I would appreciate your advice on a mechanistic model of neighborhood dynamics. My goal is to estimate the intensity of competition among trees. I have 45 sites, containing all possible combinations of four tree species. In each site, 6 trees of each species have been cored, meaning that I have observations of growth on up to 24 “focal” trees per site. In addition, I have measured the diameters of all trees within 10 m of every focal tree. For these “non-focal” trees, I have no record of growth, but only their final sizes. I estimate growth rates using a nonlinear model of growth that peaks early, and slows as trees grow. I reduce the rate of growth by calculating a neighborhood crowding index “NCI”, which increases with the number of neighbors and their size, and decreases with their distance to the target tree. It is easy to calculate NCI for every tree in the final year, because at that time, I have observed sizes for every tree. However, to take advantage of growth data from the cores, I also want it for previous years. Therefore, I estimate the size of every non-focal tree in each year of the past, using the current values of the growth model parameters (*Max*, *Dopt* and *K*), under the assumption that the focal trees are representative of the non-focal trees. Then, I re-calculate an inferred NCI for every tree in every year.

The issue is that the model runs very slowly. Even when only examining a single site, with 6 focal individuals and 62 neighbors, it can take 10 minutes for 100 iterations, which is far too few for convergence. And I’d like to run it for the full dataset, using 45 sites, and many hundreds of trees.

To reduce the number of loops, I have attempted to re-write the code, looping over individuals, rather than observations. However, this requires exponentiating row_vectors pulled from the matrix *DBH_all_n*, which is not allowed in stan. Also, please note that the exponentiation is essential. Eventually, I plan to estimate the effects of distance to neighbors and the the effect of the size of neighbors. The exponents 2.5 and 1.8 are temporary placeholders.

So, I would appreciate any advice you can offer about how I could speed up the performance of the model which is pasted here below.

thanks, and best wishes,

–tim

```
data {
// General information
int<lower=0> Y; // Number of years in the observations
int<lower=0> S; // Number of sites in the observations
int<lower=0> SPP; // Number of species in the dataset
int<lower=0> N_years; // What's the maximum number of years for which a tree was observed?
int<lower=0> N_obs; // how many obervations were made, in total?
// Scalars with information about all individuals
int<lower=0> N_all; // How many individuals are there, in total?
int<lower=0> N_neis_tot; // How many neighbors are there, in total?
int<lower=0> N_dists; // The total number of distances between all pairs of trees.
// Vectors with information about every individual
int<lower=0> nei_pos_start[N_neis_tot]; // At what position in the vector do the neighbors for each tree begin?
int<lower=0> N_neis[N_neis_tot]; // How many neighbors are there for each non-focal tree?
int<lower=0> serials_of_neis[N_dists]; // What are the serial numbers of every neighbor of every tree?
real<lower=0> dists[N_dists]; // A vector, sorted by serial number, then by their neighbors, of distances from each tree to its neighbors
int<lower=0> type[N_obs]; // What type of tree is this obs made on (focal = 1, or neighbor = 2)?
int<lower=1> year[N_obs]; // In what year what the observation made?
int<lower=1> site[N_obs]; // In which site did this tree grow?
int<lower=0> serial[N_obs]; // What is the serial number of this tree?
int<lower=1> species[N_obs]; // what species is this tree?
int<lower=0> N_RW_foc; // Number of observations of RW on focal trees
int<lower=0> N_RW_nf; // Number of observations of RW on non-focal trees
real<lower=0> RW_foc[N_RW_foc]; // Response variable - Ring width (in mm) of focals
real<lower=0> RW_nf[N_RW_nf]; // Ring width (in mm) of non-focals
matrix<lower=0>[N_all,N_years]DBH_all; // A matrix, representing the DBH of every tree in every year
}
parameters {
real<lower=0>mu_Max; // Overall mean Max across the 4 species
real<lower=0>mu_Dopt; // Overall mean Dopt across the 4 species
real<lower=0>mu_K; // Overall mean K across the 4 species
real<lower=0>mu_NCI_eff; // Overall mean NCI_eff across the 4 species
real<lower=0>sigma_Max; // Standard deviation around the overall mean Max
real<lower=0>sigma_Dopt; // Standard deviation around the overall mean Dopt
real<lower=0>sigma_K; // Standard deviation around the overall mean K
real<lower=0>sigma_NCI; // Standard deviation around the overall mean NCI_eff
real<lower=0>Max[SPP]; // Maximum growth rate, one per species
real<lower=0>Dopt[SPP]; // Size, in cm, at which RW is greatest, one per species
real<lower=0>K[SPP]; // Kurtosis - breadth of the peak in growth, one per species
real<lower=0> NCI_eff[SPP]; // Effect of NCI on growth, one per species. Cannot be negative, as that would mess with the interpretation
real<lower=0> overall_sigma; // Error scale
}
model {
real preds[N_RW_foc];
matrix [N_all,N_years]RW_all_n = rep_matrix(0, N_all, N_years); // temporarily hold this iteration's predicted RW for non-focals, before cumsumming them together to produce a predicted DBH for that tree.
matrix [N_all,N_years]DBH_all_n = DBH_all; // make a temprary copy of DBH_all for each iteration, so that it can be updated with new estimates of tree size.
int counter = 1;
// Priors on hyperparameters
mu_Max ~ normal(4, 3); // Priors for max growth rate
mu_Dopt ~ normal(120,40); // Prior for optimal DBH
mu_K ~ normal(1.5, 1); // Prior for Kurtosis
mu_NCI_eff ~ normal(0.1, 1); // Prior for the NCI_effect
sigma_Max ~ normal(0, 5); // Prior for the max growth rate
sigma_Dopt ~ normal(0, 40); // Prior for the optimal DBH
sigma_K ~ normal(0, 5); // Prior for the Kurtosis
sigma_NCI ~ normal(0, 5); // Prior for the NCI effect
// Hyper parameters for species-level effects
for (s in 1:SPP){
// Max: Mean maximum growth rate
// Max[s] represents the maximal RW for each species.
// mu_Max represents maximal RW for the average species.
// sigma_Max represents the sd of maximal RW among species.
Max[s] ~ normal(mu_Max, sigma_Max);
// Dopt: Size, in mm, at which RW is greatest
// Dopt[s] Size, in mm, at which RW is greatest for each species.
// mu_Dopt Size, in mm, at which RW is greatest for the average species.
// sigma_Dopt represents the sd of DBH effects among species on RW.
Dopt[s] ~ normal(mu_Dopt, sigma_Dopt);
// K: Kurtosis - breadth of the peak in growth
// K[s] represents the breadth of the peak in growth for each species.
// mu_K represents the breadth of the peak in growth for the average species.
// sigma_K represents the sd of breadth of the peak in growth among species
K[s] ~ normal(mu_K, sigma_K);
// NCI_eff: Effect of crowding on RW, considering all neighbors to be functionally identical.
// NCI_eff represents the effect of crowding (NCI) on the RW of each species.
// mu_NCI_eff represents the effect of NCI on RW for the average species
// sigma_NCI represents the sd among species in the NCI-RW relationship
NCI_eff[s] ~ normal(mu_NCI_eff, sigma_NCI); //
}
// LOOP over all individuals, to estimate their growth patterns, and, if focal, to compare against observations
for(n in 1:N_obs){ // loop over all observations of non-focals, to predict their RW. THis will then be used in the next loop to compute NCI, updated each iteration.
int year_n = year[n];
int serial_n = serial[n];
int nei_pos_start_n = nei_pos_start[serial_n];
int species_n = species[n];
real NCI_n = 0.0;
real prediction_n;
for(a in 0:(N_neis[serial_n]-1)){ // loop over all neighbors of the focal tree in observation n
NCI_n += (DBH_all_n[serials_of_neis[nei_pos_start_n + a],year_n]^1.8) / (dists[nei_pos_start_n + a]^2.5);
}
prediction_n = (Max[species_n] / (1.0 + NCI_eff[species_n] * NCI_n)) *(exp(-1.0/2.0*(log(DBH_all_n[serial_n,year_n] / Dopt[species_n]) / K[species_n])^2.0));
if(type[n] == 1){
preds[counter] = prediction_n;
counter = 1 + counter;
} else {
RW_all_n[serial_n,year_n] = prediction_n;
}
}// end loop over observations
RW_foc ~ normal(preds, overall_sigma);
for(i in 1:rows(RW_all_n)){
DBH_all_n[i] = cumulative_sum(RW_all_n[i]);
}
}
```