Help with slow model

Modeling
1

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]);
  }
}
2

Try this. You can remove the loop for the priors as these are all vectorized. Next, NCI_n is actually just created from data so we can put that loop into the transformed data block and now this only runs once at the start of the program instead of for each mcmc iteration.

In the future when you want to estimate the exponents. I’d set each of those exponents to 1 in the transformed data block and take the NCI_n_lognumerator = log(num) and NCI_n_logdenominator = -log(den). Then in the model block do,

transformed data {
...
for(a in 0:(N_neis[serial_n]-1)){ 
NCI_n_lognumerator[n] = log(DBH_all_n[serials_of_neis[nei_pos_start_n + a], year]);
NCI_n_logdenominator[n] = -log(dists[nei_pos_start_n + a]);
}
}
...
model {
...
vector[N_obs] NCI_n = exp(alpha[1] * NCI_N_lognumerator + alpha[2] * NCI_N_lognumerator);
...
}

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
} 
transformed data{
  vector[N_obs] NCI_n = rep_vector(0.0, N_obs);
  
  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];
  for(a in 0:(N_neis[serial_n]-1)){ // loop over all neighbors of the focal tree in observation n
    NCI_n[n] += (DBH_all_n[serials_of_neis[nei_pos_start_n + a], year]^1.8) / (dists[nei_pos_start_n + a]^2.5);
  } 
}
}

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
    // 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 ~ 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 ~ 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 ~ 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 ~ 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 prediction_n = (Max[species_n] / (1.0 + NCI_eff[species_n] * NCI_n[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]);
  }
}
1 Like
3

Dear @spinkney,
Thanks for your quick reply and helpful suggestions. I have removed the loop over species for setting hyperparameters. I really appreciate the log and exponentiation workaround for exponentiating row_vectors. I should have thought of that myself. I will retain that in the model block, however, because I need to update NCI in every iteration of the model. The functions add_iter() and get_iter() are used to count the number of iterations performed in the model. Using them allows me to only update NCI once the draws from the posterior have begun to stabilise.

At this point, I would like to ask if you see any further ways that I can streamline or speed up the model. Any and all advice is welcome!
Thanks, and best wishes,
–tim

The model now looks like this:

 functions {
     void add_iter();
     int get_iter();
  }

data {
  // General information
  int<lower=0>  N_spp;           // Number of species in the dataset
  int<lower=0>  N_sites;         // Number of sites in the observations
  int<lower=0>  N_years;         // What's the maximum number of years for which a tree was observed?
  int<lower=0>  max_N_neis;      // How many neighbors does the individual with the most neighbros have? this defines the  sizes of several data matrices.
  int<lower=0>  N_all;           // How many individuals are there, in total?
  int<lower=0>  N_dists;         // How many distances are there, between trees and their neighbors?

  // Matrices
  vector<lower=0>[N_dists]dists;           // A vector of distances from trees to each of their neighbors
  int<lower=0>serials_of_neis[N_dists];    // What are the serial numbers of every neighbor of every tree? THis is organised so that every tree is a row, and every column contains the serial numbers of its neighbors.
  matrix<lower=0>[N_all,N_years]DBH_all;   // A matrix, representing the DBH of every tree in every year
  matrix<lower=0>[N_all,N_years]NCI;       // A matrix, representing the NCI of every tree in every year
  matrix<lower=0>[N_all,N_years]RW;        // A matrix, representing the ring width of every tree in every year

  // Vectors
  int<lower=0>  start_pos  [N_all];       // What is the starting position in the vector of neighbors of this tree?
  int<lower=0>  species    [N_all];       // What is the species of this tree?
  int<lower=0>  type       [N_all];       // What type of tree is this obs made on (focal = 1, or neighbor = 2)?
  int<lower=1>  site       [N_all];       // In which site did this tree grow?
  int<lower=0>  N_neis     [N_all];       // How many neighbors does this tree have?
  real<lower=1> marginality[N_all];       // How marginal is this tree?
}

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>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[N_spp];      // Maximum growth rate, one per species
  real<lower=0>Dopt[N_spp];     // Size, in cm, at which RW is greatest, one per species
  real<lower=0>K[N_spp];        // Kurtosis - breadth of the peak in growth, one per species
  real<lower=0,upper=2>NCI_eff;// Effect of NCI on growth. Cannot be negative, as that would mess with the interpretation
  real<lower=0,upper=3> size_eff;// effect of size on NCI, global - not varying among species
  real<lower=0,upper=3> dist_eff;// effect of distance on NCI, global - not varying among species
  row_vector    [N_sites]beta_site;// Parameter to capture variation in mortality among sites
  row_vector    [N_years]beta_year;// Parameter to capture variation in mortality among years
  row_vector    [N_all]  beta_tree;// Parameter to capture variation in growth among individual trees
  real<lower=0> sigma_site;     // Standard deviation among sites
  real<lower=0> sigma_year;     // Standard deviation among years
  real<lower=0> sigma_tree;     // Standard deviation among trees
  real<lower=0> overall_sigma;  // Error scale
}


model {
  matrix [N_all,N_years]    DBH_all_n = DBH_all; // This iteration's DBH_all for every tree (rows) and years (cols)
  matrix [N_all,N_years]    pred;                // The predicted RW for each individual, before cumsumming them together  to produce a predicted DBH (for non-focal trees), or using them to compare against teh observed RWs (for focal trees)
  matrix [N_all, N_years]   NCI_n = NCI;        // This iteration's crowding index for a single tree (rows) and all years (cols)

  // 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
  sigma_Max  ~ normal(0,   5);    // Prior for the max growth rate
  sigma_Dopt ~ normal(0,  30);    // Prior for the optimal DBH
  sigma_K    ~ normal(0,   5);    // Prior for the Kurtosis
  sigma_year ~ normal(0,   2);    // Prior for the standard deviation among years
  sigma_site ~ normal(0,   2);    // Prior for the standard deviation among sites

  // Hyper parameters for species-level effects
  Max     ~ normal(mu_Max,     sigma_Max);
  Dopt    ~ normal(mu_Dopt,    sigma_Dopt);
  K       ~ normal(mu_K,       sigma_K);

  // LOOP over all individuals
  for(i in 1:N_all){
    int species_i = species[i];
    if(get_iter() > 20000){
      vector [N_neis[i]] dists_neis_i = log(dists[start_pos[i]:(start_pos[i] + N_neis[i]-1)])*dist_eff;   // The distances to the neighbors of tree i.
      matrix [N_neis[i], N_years]DBHs_neis_i = log(DBH_all_n[serials_of_neis[start_pos[i]:(start_pos[i] + N_neis[i]-1)]])*size_eff;   // The DBHs of the neighbors of tree i. Neighbors (rows) and years (cols).
      for(y in 1:N_years){
        NCI_n[i,y] = sum(NCI_eff * exp(col(DBHs_neis_i, y) - dists_neis_i)) * marginality[i];
      }
    }
    pred[i] = exp(-1.0 * NCI_n[i]) * Max[species_i] .* exp(-1.0/2.0*(log(DBH_all_n[i] / Dopt[species_i])/K[species_i]) .* (log(DBH_all_n[i] / Dopt[species_i]) / K[species_i])); 
    if(type[i] == 1){ // NON-FOCAL individuals
      if(get_iter() > 20000){
        DBH_all_n[i] = cumulative_sum(pred[i]);
        DBH_all_n[i] =  DBH_all_n[i] + DBH_all[i,N_years] - max(DBH_all_n); // this pushes the DBHs back up to their observed sizes at the end of teh study
      } else {          // FOCAL individuals
      target += normal_lpdf(RW[i] | (pred[i] + beta_year + beta_site[site[i]] + beta_tree[i]), overall_sigma); // Compare predictions and observations on a normal scale for individual n. 
    }
  }

  // Random effects for sites & years
  beta_site ~ normal(0, sigma_site);//sigma_site represents the standard deviation in RW among sites
  beta_year ~ normal(0, sigma_year);//sigma_year represents the standard deviation in RW among years
  beta_tree ~ normal(0, sigma_tree);//sigma_tree represents the standard deviation in RW among trees
}
generated quantities{
  add_iter(); // increment the counter after each iteration
}