Spatial models for plant neighborhood dynamics in Stan
Cristina Barber, Andrii Zaiats, Cara Applestein and T.Trevor Caughlin
August 20, 2020
1. Introduction
Interactions between neighboring plants impact how plants grow, survive, and reproduce. Although interactions occur at the scale of individual plants, the consequences of plant-plant interactions shape population and community structure across terrestrial ecosystems. Plants tend to do worse in single-species neighborhoods than in many-species neighborhoods (Sortibrán et al. 2014), a plant level dynamic that helps explain how plant biodiversity is maintained across ecosystems, from montane deserts (Adler et al. 2010) to tropical rainforests (Comita et al. 2010). Managing plant neighborhoods, from thinning dense stands of trees (Cescatti and Piutti 1998, Lechuga et al. 2017) to planting species that will facilitate their neighbors (Gómez-Aparicio 2009), is a cornerstone of forestry, restoration, and agriculture. The importance of plant-plant interactions across basic and applied ecology points to the need for analyses to quantify how neighborhoods impact plant demography. Such analyses must account for space, as plants interact more with closer neighbors than with neighbors further away (Fig.1).
Figure 1: Spatial structure of a plant neighborhood. The seedling in the center of the plot experiences a range of plant-plant interactions, depending on the neighbor’s species identity, plant size, and physical distance.
The objective of this case study is to illustrate how to fit spatially-explicit models for plant-plant interactions in Stan. Using multiple examples of plant demographic rates, including growth, survival, and recruitment, we address the following objectives:
1. Develop models that account for non-linear relationships between interactions strength and pairwise distances between neighboring plants:
Explicitly estimating distance-decay parameters in neighborhood models enables greater flexibility and opportunities for novel inference. We use simulated data to demonstrate Stan’s ability to fit non-linear models for plant-plant interactions.
2. Efficient strategies for modeling pairwise interactions between many plants:
Spatially-explicit models require accounting for pairwise interactions between neighbors, leading to large matrices. Every individual element added to a pairwise distance matrix increases matrix size following a series of square numbers. However, not all plants represented in the pairwise matrix interact, so this kind of data accumulates many zeroes, which creates a sparse matrix. We will demonstrate how to efficiently handle sparse matrices in plant neighborhood models, using Stan’s segment function to eliminate zero-valued elements from computation.
3. For hierarchical neighborhood models, we explore the costs and benefits of the non-centered parameterization for model convergence:
Hierarchical structures, including random effects that represent spatial and temporal covariance, are common in ecological data. We show how plant-plant interaction models can incorporate random effects and then explore the consequences of centered and non-centered parameterization (Betancourt and Girolami 2013) for model-fitting.
2. Environment
This case study uses Stan version 2.19.1 and the following R packages.
library(rstan)
library(ggplot2)
library(bayesplot)
library(rethinking)
library(gridExtra)
library(kableExtra)3. Simulation: Effect of neighbor plants on growth
We start by simulating a spatially explicit data set of plant growth and then fitting a neighborhood model. We model the growth of stationary plants as a function of their intrinsic growth (growth in isolation) and their local neighborhood characteristics (i.e., neighbor size and proximity; Eq. 3).
\[\begin{split} growth\sim normal(\mu,\sigma ^2)\\ \ \mu_j = \alpha + \beta \textrm{s}_0 + a_3 \sum\limits_{j=1}^n \frac {\textrm{size}_{j}^{a_1}} {\textrm{distance}_{j}^2 a_2} \quad Eq. 3 \end{split}\]
Where \(\alpha\) is the intercept, \(\beta\) is the growth rate of plants depending on their initial size, \(\mathit{\textrm{s}_0}\) is the initial size of the target plant. Size, and distance are the matrices containing the size and distance of neighbors, and \(\mathit{a_1}\), \(\mathit{a_2}\) and \(\mathit{a_3}\) are parameters describing the relationship of growth with the neighbors’s size and distance.
3.1 Data simulation
The generative model includes two predictor variables (Fig.2); the size of plants at time t and t+1, and a spatial component; the x and y coordinates for each individual:
N=500
size_t0=rnorm(N,mean=10, sd=1)
size_t1=rnorm(N, mean=35,sd=1)
x=runif(N,min=0,max=14)
y=runif(N,min=0,max=14)
distmat=as.matrix(dist(cbind(x,y),method="euclidean"))
sizemat=matrix(rep(size_t0,500),nrow=500,ncol=500,byrow=T)
Figure 2: Map showing the simulated plant individuals and their initial size. The dots show where the plants are, and the size of the dot is proportional to the initial size of the plant.
Our generative model assumes that plant-plant interactions beyond a neighbor effect radius of 10 m are zero. We then simulate the effect of neighbors as a function of neighbor size (bigger neighbors have larger effects) and distance (closer neighbors have larger effects) using the following equation: growth=alpha+size_t0*beta+ a3*rowSums(exp(a1*log(sizemat4)-(distmat4^2*a2))). We assume growth is normally distributed, with a mean equal to the linear combination of an intercept (representing baseline growth), a size-dependent term, and the neighborhood effect term (Eq. 3).
#Here, the parameters are set to produce a biologically reasonable outcome:
alpha=0.8 # Intercept.
beta=0.2 # Size dependent growth (we hypothesize that bigger individuals grow faster).
a3=.035 # Neighbor effect parameter
a2=0.05 # Scaling parameter which determines the spatial reach of neighbor effect
a1=0.1 # Size dependent effect of the neighbor.
true_param=c(alpha,beta,a3,a2,a1,sd=0.01)
growth=rnorm(N,
mean=alpha+
size_t0*beta+
a3*rowSums(exp(a1*log(sizemat)-(distmat^2*a2))),
sd=0.01)The pairwise matrix representing plant-plant interactions can become large, so one solution is to assume that neighbors far from each other do not interact. This assumption establishes a neighbor effect radius, beyond which neighborhood effects are assumed to be negligible. The consequence of establishing a neighbor effect radius is that the effect of any individuals beyond the radius threshold is disregarded.
When we set the neighbor effect radius, we create a considerable amount of zero-valued matrix elements representing neighbor pairs with minimal interaction. One way to decrease computation time is to slice off zeros from the neighbor matrices as they do not contribute to the crowding term. The built-in segment function in Stan provides a way to eliminate zero-valued matrix elements from computation by storing non-zero elements and their position on the pairwise matrix in three separate vectors. The segment function requires creating a long vector of non-zero observations and two vectors indicating the position of non-zero elements within the matrix (Fig. 3).
Figure 3: Segment function requires the creation of a vector with no zeroes, an index of the non- zero values on each matrix row and an index of the position in the vector of the first values of the matrix rows. This diagram explains the composition of the non-zero values vector and the two vector index needed to use the segment function.
We compare Stan models fit with the segment function to the same models fit using the full matrix. To evaluate how choosing a neighborhood effect radius could introduce bias in parameter estimation, we compare the fit of six different models with neighbor effect radii of 5m, 10m, 15m, 20m. We note that the “true” effect radius, which is known in this case because the data were simulated, is 10 m. With real data, the decision of this radius should be based on biological knowledge, for example, the distance that roots extend from the base of a plant (Zaiats et al. 2020) or by testing at which neighbor effect radius the model has the best fit (Pacala and Silander 1985, Pacala and Silander 1987). For each neighbor effect radius, we set the effect of distant neighbors beyond the radius to zero.
distmat=ifelse(distmat>10,0,distmat) #here, the cut-off is at 10m radius. Because the neighbors further than the neighbor effect radius do not affect the response variable, we also set the size of plants in the pairwise matrix beyond the neighbor effect radius to zero.
sizemat=ifelse(distmat==0,0,sizemat)Then we create the different input databases for each of the neighbor effect radii.
# Combine the predictors in the list to create the model data base:
case_stanlist=list(N=N,growth=growth,size_t0=size_t0, distmat=distmat,sizemat=sizemat)3.2 Stan code
In our Stan code, neighbor interaction parameters are estimated in the transformed parameters block. The priors for the parameters alpha, beta, a1, a2 and, a3 are weakly informed priors following a normal(0,5) distribution. We use an exponential(1) distribution for sigma to prevent the model from exploring negative values for the standard deviation. The priors used in this example are relatively vague and could be tightened if prior information is available. The following code uses the pairwise matrix:
data{
int N; // number of individuals
vector [N] size_t0; // size of focal plants
vector [N] growth; // growth of focal plants, response
vector [N] sizemat[N]; // full size matrix
vector [N] distmat[N]; // full distance matrix
}
parameters{
real alpha; // Intercept
real beta; // Effect of target plant size on response
real sigma; // global SD
real<lower=0> a1; // size dependent effect of the neighbor
real a3; // net crowding effect
real<lower=0> a2; // spatial scale
}
transformed parameters{
vector[N] kernel;
vector[N] mu;
vector[N] smat[N];
vector[N] dmat[N];
for(i in 1:N){
for(j in 1:N){
smat[i,j]=sizemat[i,j]^a1;
dmat[i,j]=distmat[i,j]^2*a2;
}}
for(n in 1:N)
kernel[n]=sum(smat[n]./exp(dmat[n])); //we use ./ to perform a vectorial division
for(n in 1:N)
mu[n]=alpha+size_t0[n]*beta+a3*kernel[n];
}
model{
alpha~normal(0,5); // uninformed priors centered at zero;
beta~normal(0,5);
a1~normal(0,5);
a2~normal(0,5);
a3~normal(0,5);
sigma~exponential(1);
growth ~ normal(mu,sigma);
}To use the segment function, we have to create the three vectors described in Figure 3,including one vector that contains all the non-zero elements of the pairwise matrix and two index vectors that indicate the position of the non-zero values in the pairwise matrix.
# code to create the vector with no zeroes and the index vectors
n_nb=rep(NA,case_stanlist$N)#empty vector with length N (number of rows in the matrix)
#populate the vector, number of non zero observations in each row, i.e., group size index
for(i in 1:length(n_nb)) {
n_nb[i]=sum(case_stanlist$sizemat[i,]!=0)
}
case_stanlist$n_nb=n_nb #add to the list
a=as.vector(t(case_stanlist$sizemat)) #convert size matrix to a vector
b=as.vector(t(case_stanlist$distmat)) #convert distance matrix to a vector
size_observations=a[a!=0] #subset non-zero observations
dist_observations=b[a!=0] #subset non-zero observations
case_stanlist$size_vector=size_observations # add them to the list
case_stanlist$dist_vector=dist_observations
case_stanlist$obs=length(case_stanlist$size_vector) #total number of non-zero observations
pos=rep(NA,case_stanlist$N) # position index for each group
pos[1]=1
for(i in 2:(case_stanlist$N)){
pos[i]=pos[i-1]+case_stanlist$n_nb[i-1]
}
case_stanlist$pos=posThen to include the segment function in the Stan code we will substitute the size and distance pairwise matrices using the segment functions segment(size_vec, pos[n], n_nb[n]) and segment(dist_vec, pos[n], n_nb[n]) respectively. The code substituted by “…” is the same as in the previous code for the full matrix approach.
data{
int obs; // number of observation in the long vector
int N; // number of individuals
vector [N] size_t0; // size of focal plants
vector [N] growth; // growth of focal plants, response
int n_nb[N];
vector [obs] size_vector; // leng vector of neighbor sizes
vector [obs] dist_vector; // long vector of pairwise distances
int pos[N];
}
transformed data{
}
}
parameters{
...
}
transformed parameters{
vector[N] kernel;
vector[N] mu;
vector[obs] size_vec;
vector[obs] dist_vec;
for (i in 1:obs){
dist_vec[i]=dist_observations[i]^2;
size_vec[i]=size_observations[i]^a1;
}
for(n in 1:N)
kernel[n]=sum(exp(log(segment(size_vec, pos[n], n_nb[n]))-
(segment(dist_vec, pos[n], n_nb[n])*a2)));
for(n in 1:N)
mu[n]=alpha+size_t0[n]*beta+a3*kernel[n];
}
model{
...
}3.3 Output
We ran the code using the pairwise matrix and the segment function for the six neighbor effect radii and compared model fitting efficiency between them. We estimate the model fitting efficiency by dividing the elapsed time by the effective sample size (ESS) of 1000 iterations. Lower values represent more efficient sampling:
sum(get_elapsed_time(model)[,"sample"])/sum(summary(model)$summary[,'n_eff'] * 1000)The results show that the segment function is more efficient than the pairwise matrix for almost all neighbor effect radii. However, for the generative radius at 10m, both strategies are equally efficient. Using the segment function to fit neighbor interaction models is advantageous because the greater efficiency from the segment function compared with the matrix under different radii allows a faster exploration of models parametrized with different radii (Fig. 4).
Figure 4: Segment function is more efficient than the matrix approach for all neighbor effects, except 10 m. We note that 10 m was the “true” neighbor effect radius used to simulate data. The blue dashed line shows the change in efficiency using the matrix under different neighbor effect radii. The green line shows the change in efficiency using segment function under different neighbor effect radii. Smaller values of Time/ESS represent inreased efficiency and higher values represent decreased efficiency.
We can also check whether estimated parameter intervals capture the “true” values because the parameters values are known in this simulated example. Our results highlight the importance of choosing a neighbor effect radius to obtain accurate estimates that describes the plant-plant interactions (Fig. 5). The matrix and the segment function could recover the true parameters with a neighbor effect radius of 10 m. For the segment method, similar efficiency and parameter recovery was achieved with a neighbor effect radius of 15 m. Models fitted with other neighbor effect radii did not converge and were not able to recover the true parameters for both the matrix and the segment function. These results agree with Canham and Uriarte (2006). They determined that a slightly bigger radius than the mean radius can provide estimates with negligible biases. In comparison, a radius smaller than the mean radius provided more significant biases because it missed contributing individuals. Altogether, we would recommend erring on the side of overestimating the neighbor radius effect rather than having a too-small radius. Given this recommendation, we also suggest using the segment function as it is more efficient than the matrix for all radii tested (Fig.4).
Figure 5: The neighbor effect radius selection is essential to obtain unbiased estimates. The bias is similar for all the parameters and neighbors effect radii, so in this graph, we just show one of the parameters. The chosen parameter estimates are for a2, which estimates the rate of decrease of the effect of neighbors on growth with distance. The green shapes are the segment function parameters estimates. The blue shapes are the matrix parameter estimates. Each of the shapes corresponds to a different neighbor effect radius. The red line is the true parameter used in the simulation. 90% Credibility intervals (CI) showed in the figure. Note the CI for the 10, and 15 m radii are not visible.
Lastly, we check common diagnostic metrics to evaluate convergence at the generative radius of 10 m. There are no divergent transitions, and the ESS, and \(\hat{R}\) values do not reveal any sampling problems.
## [1] "Segment function at radius 10m:"## 0 of 8000 iterations ended with a divergence.| mean | se_mean | sd | 2.5% | 25% | 50% | 75% | 97.5% | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|---|---|---|
| alpha | 0.801 | 0.000 | 0.005 | 0.792 | 0.798 | 0.801 | 0.805 | 0.811 | 719.437 | 1.002 |
| beta | 0.200 | 0.000 | 0.000 | 0.199 | 0.199 | 0.200 | 0.200 | 0.201 | 886.849 | 1.001 |
| sigma | 0.010 | 0.000 | 0.000 | 0.009 | 0.010 | 0.010 | 0.010 | 0.011 | 7256.505 | 1.000 |
| a1 | 0.085 | 0.001 | 0.019 | 0.048 | 0.073 | 0.085 | 0.097 | 0.121 | 414.132 | 1.004 |
| a3 | 0.036 | 0.000 | 0.002 | 0.033 | 0.035 | 0.036 | 0.037 | 0.039 | 411.221 | 1.004 |
| a2 | 0.050 | 0.000 | 0.000 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 1915.463 | 1.001 |
4. Data example 1: Predicting recruitment of an invasive species based on parent tree location
4.1 Overview of the model
Our next example applies inverse modeling to estimate spatial patterns of invasive strangler fig tree recruitment. The focal tree species, Ficus microcarpa, is native to South Asia and is invading Florida. We constructed inverse models using spatial data on the abundance of seedling strangler figs in 30 x 30 m plots embedded in 300 m radius plots where the location and size of reproductive adult fig trees were recorded. The study design, described in Caughlin et al.(2012), resulted in fifty-two plots placed along a 250-km transect in Southwest Florida. The neighborhood model quantifies how seedling abundance declines with distance from adult trees (Eq.5).
\[\begin{split} recruitment \sim negative \, binomial(\mu,\phi)\\
\mu=a+b\sum \limits_{j=1}^n \frac{1}{c+\textrm{distance}_j}\textrm{CP} \quad Eq.5
\end{split}\]
The response variable in this example is recruitment, the number of seedlings within the monitored plots. Because our response variable is count data, we use the negative binomial distribution, parameterized in terms of a mean (\(\mu\)) and an overdispersion parameter(\(\phi\)). a is the intercept, b is the parameter describing the relationship between recruitment and the hyperbolic effect of distance, c is the rate of decay of the effect of distance on recruitment, and CP is the number of suitable recruitment sites (counts of host trees) in the 30 x 30 m plots.
The data on adult tree abundance are organized in a distance pairwise matrix with the number of rows equal to the total number of sites and the number of columns equal to the maximum number of adult trees observed at a site. Then the distance pairwise matrix is filled with the distances between seedling plots and parent trees. Once the matrix is filled, the distances beyond the neighbor effect radius are transformed to zero.
4.2 Stan Code
We use a very weakly informed prior normal(0,100) in a , b and j and an exponential(0.5) distribution for phy to prevent the model from exploring negative values for the shape parameters of the negative binomial. The original study did not include a log-link and exponenciated a , b and j to keep the parameters positive. To replicate the previous results, we ensure non-negative values for the mean of the negative binomial distribution by constraining parameters a , b and j to positive values using real<lower=0>, however we note that the log-link is the canonical link-function for the negative binomial distribution and probably a better way to ensure positive values for future studies. The Stan code for a model using the pairwise matrix is the following:
data {
int N; //number of available places for recruitment
int K; //number of parent trees
vector [K] dist[N]; //array containing N vectors with K distances
//of available places to parent trees
int x[N]; //count of seedlings
int CP[N]; //count of cabbage palms
vector [N] one; //vector of ones for the division
}
parameters {
real<lower=0> a;
real<lower=0> b;
real<lower=0> c;
real<lower=0> phy;
}
transformed parameters{
real mu[N];
real aa[N];
for(i in 1:N) {
aa[i] = sum(ones[i] ./(c+dist[i]));
}
for (n in 1:N){
mu[n]=(a+b*aa[n])*CP[n];
}
}
model{
a~normal(0,100);
b~normal(0,100);
c~normal(0,100);
phy~exponential(0.5);
x~neg_binomial_2(mu,phy);
}
}To use the segment function in this example, we consider that there are seedling plots that do not have any adult strangler fig trees nearby, indicating an empty neighborhood. An empty neighborhood will result in a zero in the n_nb vector, as in the second row of the predictor matrix in Figure 3. To account for these zeroes, we include an ifelse statement that prevents zero-valued elements of n_nb in the plant interactions kernel in the transformed parameters block:
data {
int N; //number of available places for recruitment
int K; //number of non-zero parent trees
vector [K] distrag; //vector containing all the non-zero distances
int x[N]; //count of seedlings
int CP[N]; //count of cabbage palms
int n_nb[N]; //vector giving the amount of non zero values
int pos [N]; //vector giving the position of non zero values
vector [N] one; //vector of ones for the division
}
parameters {
...
}
transformed parameters{
real mu[N];
real aa[N];
for(i in 1:N) {
aa[i] = sum(one[i] ./(c +segment(distrag,pos[i],n_nb[i])));
}
for (n in 1:N){
if (n_nb[n]==0){mu[n]=a;}
else{
mu[n]=(a+b*aa[n])*CP[n];
}
}
}
model{
...
}
}4.3 Output
The results of the model using matrix and segment strategies lead to similar parameter estimates (Fig. 6) and produce parameter estimates comparable to the frequentist approach used in the published paper (Fig. 7).
Figure 6: Parameters posterior density of the model describing seed dispersal of invasive strangler fig trees
Figure 7: Matrix and segment obtained similar estimates of the relationship between recruitment and the distance from parent tree. Curves show the relationship between recruitment and distance from parent tree parametrized using the sparse matrix, the segment function, and the frequentist maximum likelihood model fit using the mle4 package (Bates et al. 2015) from in the original study
While both the full matrix and the segment models show evidence of convergence, model efficiency is greater for the model fit with the segment function. Again, we estimate the model fitting efficiency by dividing the elapsed time by the effective sample size (ESS) of 1000 iterations. Fig_distmodelmatrix is the model fitted using the matrix, and Fig_distmodel is the model fitted using the segment function.
#Efficiency of the model fitted using the matrix:
sum(get_elapsed_time(Fig_distmodelmatrix)[,"sample"])/
sum(summary(Fig_distmodelmatrix)$summary[,'n_eff'] * 1000)## [1] 3.011334e-08#Efficiency of the model fitted using the segment function:
sum(get_elapsed_time(Fig_distmodel)[,"sample"])/
sum(summary(Fig_distmodel)$summary[,'n_eff'] * 1000)## [1] 8.417401e-095. Example 2: Seed germination probability
5.1 Overview of the model
Our second data example uses a Bayesian hierarchical framework to estimate how conspecific trees affect the seed to seedling transition in a tropical tree species from central Thailand. Plant neighborhoods considered in this model include seedling and adult tree neighbors, all of the same species as marked seeds placed on the forest floor in 1 \(m^2\) plots. The marked seeds comprised a seed addition experiment spanning a 5 km gradient of tree abundance. In this study, the neighbor effect radius was set at 10 m. The primary objective of the study was to quantify how the density of neighbors impacted the probability of seed germination. More information on the study can be found in Caughlin et al.(2015). We model the germination probability of seeds as proportional data using the binomial distribution (Eq. 6).
\[\begin{split} germination \sim binomial(\textrm{S,N}) \\ logit(\textrm{S})= \mu+b \textrm{Con.seedlings}+a \sum\limits_{j=1}^n \frac{\textrm{size}_j}{\textrm{distance}_j^{ger}}+e \quad Eq.6 \\ e_j \sim normal(0, \sigma) \\ \sigma \sim normal(0,1) \end{split}\]
The response variable in this example is the proportion of total added seeds within monitored plots. Input to the binomial distribution includes the total number of seeds added to each plot (N) and the number of successful germination events (S). Size and distance are the matrices containing the size and distance of neighbors. \(\mu\) is the baseline germination, b is the decrease in seed survival as a function of neighbor density, Con. Seedlings are the amount of conspecific seedling to represent the crowding effect, a is the parameter that represents the effect of neighbor size and distance on recruitment, ger is the distance decay of the effect of neighbor size and distance, and e is a random plot effect, to account for non-independence between seeds in the same plot.
The hierarchical structure of random effects, including the random plot effect in the seed germination model, can lead to model-fitting challenges. One common challenge in fitting these models is the high correlation between random effect parameters (“Neil’s funnel effect”), which causes the sampler to become “stuck” in the highly correlated area. When this happens, a non- centered parametrization can solve the pathology. However, if the non-centered parameterization is used, but Neil’s funnel does not occur, the non-centered parameterization can lead to an inefficient exploration of the probability surface. The relative performance of the centered vs. non-centered parameterizations depends on variance in the data as well as the sample size (Betancourt and Girolami 2013). In this case study, we compare a centered and a non-centered parameterization for the germination data, using a neighborhood model fit with the segment function in Stan.
5.2 Stan Code
For a centered parameterization using the segment function, we would add the random effect directly into the model g[n] = (..) + e[plots[n]]. We use a weakly informed prior normal(0,1) in a , b, mu and ger, for e we used a normal(0, sigma_plot) to center the random effects mean around zero.
data {
int N; //number of plots
int K; //number of non-zero parent trees
int M; //number of random levels
vector [K] sizeN; //matrix of neightbor size
vector [K] distN; //matrix of neightbor distances
int x[N]; //number of seedlings
int seeds[N]; //number of seeds
int am[N]; //vector giving the number of non zero values
int pos [N]; //vector giving the position of non zero values
int Cseedlings [N]; //number of conspecific seedlings
int plots[N]; //random effect of plots
}
parameters {
real a;
real b;
real<lower=0> ger;
real mu;
real e [M];
real sigma_plot;
}
transformed parameters{
real<lower=0, upper=1> s[N];
real g[N];
for (n in 1:N){
if (am[n]==0){g[n]=mu+ b*Cseedlings[n]+e[plots[n]];}
else{
g[n] = mu + b*Cseedlings[n] + a* sum(segment(sizeN,pos[n],am[n]) ./
exp(ger*log(segment(distN,pos[n],am[n])))) + e[plots[n]];
}
}
}
model{
a~normal(0,1);
b~normal(0,1);
ger~normal(0,1);
mu~nomal(0,1);
e~normal(0,sigma_plot);
sigma_plot~normal(0,1);
x~binomial_logit(seeds,g);
}
}If we would like to use a non-centered parameterization, we will add a few parameters for the random effect e in the transformed parameters block. The non-centered parametrization avoids the correlation between parameters using a linear model. The new parameters that we created are an intercept mu_plot, the random effect el[n], and a scaling parameter gamma_el. mu_plot is equivalent to the mean of the random effects, so we chose a prior normal(0,1) to maintain it around zero. el[n] is the random effect deviations from zero that follow a unit normal prior normal(0,1). The trick here is that the unit normal prior causes el[n] to be an orthogonal vector. As a consequence, the random effects parameters become uncorrelated and avoid Neil’s funnel effect. Because el[n] provides deviations at a scale of one, we have to multiply it by the scaling parameter gamma_el to bring the parameters to the model scale.
data{
...
}
parameters {
real a;
real b;
real<lower=0> ger;
real mu;
real mu_plot;
real el[M];
real gamma_el;
}
transformed parameters{
real<lower=0, upper=1> s[N];
real g[N];
real e [M];
for (n in 1:M){
e[n]= el[n]*gamma_el+mu_plot;
}
for (n in 1:N){
if (am[n]==0){g[n]=mu+ b*Cseedlings[n]+e[plots[n]];}
else{
g[n] = mu + b*Cseedlings[n] + a* sum(segment(sizeN,pos[n],am[n]) ./
exp(ger*log(segment(distN,pos[n],am[n])))) + e[plots[n]];
}
}
}
model{
a~normal(0,1);
b~normal(0,1);
ger~normal(0,1);
mu~nomal(0,1);
el~normal(0,1);
mu_plot~normal(0,1);
gamma_el~normal(0,1);
for (n in 1:N){
x~binomial_logit(seeds,g);
}
}5.3 Output
The results show that when using a centered parameterization, we obtain similar comparable parameter estimates to the study. Non-centered parameterization results in estimates further from the ones obtained in the study and with wider confidence intervals. This example illustrates how a non-centered parameterization can sometimes result in higher uncertainty than the centered parameterization.
Figure 8: In this case, the centered parameterization provided less uncertain results closer to the ones obtained in the original paper. The curves represent the relationship between recruitment and distance from the parent trees with the estimates obtained parametrizing the model using the centered and non-centered parameterization, and the Bayesian model fit using JAGS (Plummer 2003) from in the original study.
6. Final remarks
Neighbor interactions are important for plant population and community dynamics. However, implementing neighbor models presents several challenges, including large pairwise matrices for neighbor interactions and correlations between estimated parameters. We explored Stan’s segment function and concluded that the segment function provides a way to efficiently handle sparse matrices while providing near-identical estimates to models with the pairwise matrix as input. Our recommendation is to use the segment function as a default approach for handling sparse matrices that result from pairwise plant-plant interactions.
We have also demonstrated how hierarchical terms, including random effects to account for plants co-located within the same plot, can be incorporated into neighbor interaction models in Stan. Fitting hierarchical models can be challenging due to correlations between latent variables and their variance. Our case study shows how two options for handling pathologies associated with the correlation between parameters, centered and non-centered parameterization, can be incorporated into neighbor interaction models. In our case, the non-centered parameterization performed worse than the centered parameterization. However, we expect that the non-centered parametrization will perform better for some model formulations (Betancourt and Girolami 2013). Which parametrization should be used to fit the model is case-specific, and both parametrizations should be explored. In conclusion, our study demonstrates that Stan provides a robust Bayesian framework for modeling plant neighbor interactions, with the flexibility required to overcome inefficiencies associated with sparse matrices and incorporate hierarchical effects common in ecological data.
7. Acknowledgments
Support for this case study was provided by the National Science Foundation under grant #1415297 in the SBE program. We thank Jonah Gabry for his helpful comments on this case study and the participants in the StanCon2020 for the engaging discussions about this case study.