Learning how a tree branches out: A statistical modeling approach.

The increasingly large size of the graphical and numerical data sets collected with modern technologies requires constant update and upgrade of the statistical models, methods and procedures to be used for their analysis in order to optimize learning and maximize knowledge and understanding. This is the case for plant CT scanning (CT: computed tomography), including applications aimed at studying leaf canopies and the structural complexity of the branching patterns that support them in trees. Therefore, we first show after a brief review, how the CT scanning data can be leveraged by constructing an analytical representation of a tree branching structure where each branch is represented by a line segment in 3D and classified in a level of a hierarchy, starting with the trunk (level 1). Each segment, or branch, is characterized by four variables: (i) the position on its parent, (ii) its orientation, a unit vector in 3D, (iii) its length, and (iv) the number of offspring that it bears. The branching structure of a tree can then be investigated by calculating descriptive statistics on these four variables. A deeper analysis, based on statistical models aiming to explain how the characteristics of a branch are associated with those of its parents, is also presented. The branching patterns of three miniature trees that were CT scanned are used to showcase the statistical modeling framework, and the differences in their structural complexity are reflected in the results. Overall, the most important determinant of a tree structure appears to be the length of the branches attached to the trunk. This variable impacts the characteristics of all the other branches of the tree.

1 Introduction

Understanding the structural complexity of tree branching patterns, to explain light interception by leaf canopies and incorporate the information in the modeling of biological processes such as photosynthesis, has been the objective of numerous studies [1-6]. The process by which tree branches divide and subdivide starting from the trunk has been studied more particularly in relation to space occupancy and cover by the leaf canopy [4, 5]. Fractals have been used to quantify the complexity of tree crown architecture [4, 5], but other approaches such as probability models that rely less on some self-similarity assumption are worth investigating and assessing to model tree branching patterns provided the required data are available.

Since the early 2000s, modern technologies such as computed tomography (CT) scanning have been diverted from their original design (i.e., medical in the case of CT scanning), for applications with plants in general and trees in particular [7-9]. Tree crowns and plant leaf canopies [4, 5, 7] may have been CT scanned less often than root systems (for example, see [10, 11] and the review in [12]). This is likely because of the smaller size in general of a root system, the ‘hidden half’ of a plant, relative to the canopy, but both plant structures can be CT scanned within their respective limits. Thus, plant biologists are given access to representative and accurate 3D spatial data sets of a novel type to explore, once duly processed after they were collected non-invasively and non-destructively. Access to original 3D spatial data for crowns of trees of small size represents a tremendous opportunity for modelers, including statisticians and biomathematicians. Furthermore, information of this type about tree geometry, if extended, could enter predictive models for the abundance of tree epiphytes; see [13, 14].

Here, we characterize tree branching patterns with statistical models. We privilege an approach based on statistical modeling in our analysis, as an alternative to fractal geometry and Lindenmayer systems. Still, some similarity can be seen with the generation of 3D bush-like structures by a bracketed L-system [15, Fig 3.3]. We follow a systematic approach, i.e. no random sampling is performed, and develop and fit hierarchical models for key features of the branching pattern of three miniature conifers that had their crowns CT scanned in another study [5].

2 Computed tomography scanning of tree crowns

Regardless of the nature of the specimen, a basic principle of CT scanning technology is X-ray attenuation after some exposure time [16]. For a large number of “voxels” (i.e., 3D extension of 2D pixels), indirect measures of material density, called “CT numbers” (CTN), are computed from coefficients of X-ray attenuation measured by detectors. Accordingly, a CT scanner is calibrated so that CTN = −1000 for air and 0 for water. In a CT scanning session, the scale of observation refers to the size of the scanned specimen, a miniature tree in this work, while the scale of resolution is given by the dimensions of a voxel, and there is a difference of three orders of magnitude between the two scales; see Table 1 in [17]. The X-ray tube current and voltage were 50 mA and 120 kV for all the CT scanned miniature conifers from [5], and the X-ray doses that the specimens received in one exposure were far from being lethal [18].

The three miniature conifers used for examples here are two white spruces: one Picea glauca Pixie (height: 22.1 cm; diameter: 12.7 cm) and one Picea glauca Cy’s Wonder (height: 17.5 cm; diameter: 20.9 cm), plus one Norway spruce, Picea abies Tompa (height: 17.1 cm; diameter: 14.0 cm). These trees are presented in Table 2 and Fig 5 (row-column 2–4, row-column 2–3 and row-column 2–2, respectively) in [5]. Our work focusses on tree branching patterns, the CTN of branch voxels being in a specific range, but it is possible to display the tree crowns (including leaves), the leaf voxels having a CTN in a different range; see Fig 5 of [5]. The dimensions of a voxel depend on the CT scanning settings used for the tree. They are, respectively, 0.35 × 0.35 × 0.4 mm3, 0.47 × 0.47 × 0.4 mm3, and 0.29 × 0.29 × 0.4 mm3 for the ‘Pixie Tree’, ‘Wonder Tree’, and ‘Tompa Tree’ (surnames used hereafter). Scanning these three trees yielded about 145, 115, and 112 million CTNs (including air voxels with CTN of −1000), contained in 3D arrays made of 512 × 512 matrices. These are the raw data available for processing and analysis.

The skeletal branching patterns presented in Fig 1a, 1c and 1e below were prepared in a customised MATLAB (The MathWorks, Inc., Natick, MA, USA) graphical unit interface, by tracing branches using the 3D array of CT scanning data collected for the crown of each of the three miniature conifers [5]. Such tree branching patterns are said to be “skeletal” because they have a thickness of 1 voxel. These skeletons are at the basis of our statistical modeling, but are replaced by 3D branched structures in which each branch at each level is a line segment; see Fig 1b, 1d and 1f. This replacement is explained in detail in the next section.

Fig 1

(a)-(b) The Picea glauca Pixie specimen, alias ‘Pixie Tree’. (c)-(d) The Picea glauca Cy’s Wonder specimen, alias ‘Wonder Tree’. (e)-(f) The Picea abies Tompa specimen, alias ‘Tompa Tree’. (a) Skeletal branching pattern of Pixie Tree (obtained by CT scanning). (b) Analytical representation of the branching pattern for Pixie Tree. (c) Skeletal branching pattern of Wonder Tree (obtained by CT scanning). (d) Analytical representation of the branching pattern for Wonder Tree. (e) Skeletal branching pattern of Tompa Tree (obtained by CT scanning). (f) Analytical representation of the branching pattern for Tompa Tree.

3 Analytical representation of tree branching patterns as hierarchical sets of line segments

Following the graphical and quantitative analyses of CT scanning data reviewed in Section 2, a skeleton of tree branching pattern is produced where a branch is represented by a curved line in 3D. The goal of Section 3 is to explain how this skeleton can be approximated with a hierarchical set of line segments in 3D, called “analytical representation” of the tree branching pattern. This summarizes the large CT scanning data set in a small spreadsheet and will allow the characterization of a tree crown structure using the descriptive statistics and statistical models presented in Section 4.

In the construction of a line segment for a branch in 3D, the origin is identified as the point with 3D spatial coordinates in the 3 × 1 vector o = (x, y, z)⊤, at which the branch emanates from its parent, and the terminal point has coordinates e = (x, y, z)⊤. The direction v of a branch is calculated as v = (e − o)/‖e − o‖, with ‖e − o‖, the Euclidean distance between the end and origin of the branch, considered to be the length ℓ of the branch; v is a vector in S2, the unit sphere in 3D space. The position x of the branch relative to its parent is represented by the Euclidean distance between the origin of the branch o and that of its parent, say o, divided by the length of the parent branch, so the value of x belongs to the interval (0,1). For a given branch, n is the number of offspring branches that emanate from the branch considered. Thus, a branch in our analytical representation for a tree branching pattern is represented by the following four variables:

x the position of the branch relative to the parent (a non-negative real number smaller than 1);

v the branch orientation (a unit vector in 3D) is a 3 × 1 vector;

ℓ the branch length (a positive variable);

n the number of offspring, (a non-negative integer);

The results of the application of this analytical representation to the branching patterns of the three miniature conifers introduced in Section 2, as constructed from CT scanning data given in Fig 1a, 1c and 1e), are shown in Fig 1b, 1d and 1f; the similarity between the left and right panels of the same tree is noticeable. In the hierarchy, the trunk is the level 1 branch; the branches attached to the trunk are level 2 branches; level 3 branches originate from level 2 branches; and level 4 branches, from level 3 branches. Level 5 branches are very few and very small, and are not included in our data sets.

To develop a nested set of line segments, the ancestors must be identified for each terminal segment. A terminal branch (or terminal segment in the analytical representation) is one that has no offspring, that is, n = 0. All level 4 branches are terminal in our examples, since the data sets do not contain level 5 branches. The data sets have one row for each terminal branch, with values for the 4 variables (x, v, ℓ, n) for each ancestor and the terminal branch; see Table 1. The data matrix has four sets of columns, one for each level of the hierarchy. In Table 2 a subscript denotes the level of a variable, for instance x3 is the position of a level 3 branch on its level 2 parent. There are 7 columns per level: one for identifying the branch on its parent, id, one for the position x, three for the entries of v, one for the length ℓ, and another for the number of offspring n, in this order (level 4 has only 6 columns has its branches have no offspring). The number of rows in the data matrix is equal to the number of level 4 branches plus the numbers of level 2 and level 3 branches that do not bear offspring.

Table 1

The data entries for a level 4 branch of Pixie Tree.

id	x	Entries of v			ℓ	n
1	0	0	0	1	1	46
2	0.07	-0.69	0.12	0.72	0.32	8
1	0.34	-0.90	-0.38	0.21	0.10	3
2	0.98	-0.90	0.21	0.38	0.02	NA

Table 2

Description of the variables used in Section 4.

Level	Variable	Description
2	x 2	Position of a level 2 branch, standardized relative to the length of the trunk.
	c 2	Directional cosine between the level 2 branch and the trunk.
	ℓ 2	Length of the level 2 branch, standardized relative to the length of the trunk.
	n 2	Number of offspring emanating from a level 2 branch.
3	x 3	Position of a level 3 branch, standardized relative to the length of the level 2 parent branch from which it emanates.
	c 3	Directional cosine between the level 3 branch and the level 2 branch from which it emanates.
	ℓ 3	Length of the level 3 branch, standardized relative to the length of the trunk.
	n 3	Number of offspring emanating from a level 3 branch.
4	x 4	Position of a level 4 branch, normalized relative to the length of the level 3 parent branch from which it emanates.
	c 4	Directional cosine between the level 4 branch and the level 3 branch from which it emanates.
	ℓ 4	Length of the level 3 branch, standardized relative to the length of the trunk.

Table 1 presents the data for a level 4 branch of Pixie Tree; id identifies a branch on its parent at each level. It varies between 1 and the number of branches n of its parent. The first row gives information about the trunk. The trunk direction is the z–axis, as its v vector is (0, 0, 1)⊤. In the first row, n = 46 means that Pixie Tree has 46 level 2 branches. The second row gives information about a level 2 branch. Since id = 2, it is the second level 2 branch its position x2 = 0.07 means that its origin is the point (0, 0, 1)⊤0.07 on the trunk. Its direction is v2 = (−0.69, 0.12, 0.72)⊤ and its length is ℓ2 = 0.32. An analytical expression for the line segment of the level 2 branch in Table 1 is (0, 0, 1)⊤ 0.07 + (−0.69, 0.12, 0.72)⊤ 0.32 × z, for z ∈ (0, 1) where z are the evaluation points in the interval (0, 1). The origin of the level 3 branch in Table 1 corresponds to the value z = x3 = 0.34 on that segment. Analytical representations of the four segments are constructed recursively; that for the level 4 branch in Table 1 is

To facilitate comparisons between tree branching patterns, two standardizations are carried out. First, the length ℓ is divided by the trunk length so that the former is relative to the latter. Second, the direction of the trunk is the z–axis in the coordinate system where the directions v are recorded. This is done by multiplying all the unit vectors in the data set by a rotation matrix that maps the direction of the trunk to the vector (0, 0, 1)⊤. The analytical representations for the three trees investigated in this work, normalized in this way, are available as three text files in the S1 Data.

Table 2 presents the variables, ordered from level 2 to level 4, to which the statistical models presented in Section 4 will be fitted. The c-variable is the cosine between between the direction of an offspring and that of its parent. For instance . For the terminal branch presented in Table 1, its numerical value is The variable c measures the agreement between offspring and parent orientations. Its maximum value of 1 occurs when these two unit vectors are equal.

In Section 4, the variables of Table 2 are both dependent and independent variables. The determination of independent variables follows the hierarchy of Table 2. The explanatory variables for a given dependent variable are those above that variable in Table 2. Thus the model for x2 does not have any explanatory variable whereas that for ℓ4 can have up to 10 explanatory variables. This hierarchy corresponds, at least approximately, to the time at which these variables can be measured. The parents variables (x, c, ℓ, n) are, in general, determined before those of their offspring. Among the 4 variables for a branch, x can be measured before c and ℓ and n occur after x and c. The large number of candidate explanatory variables highlights the importance of variable selection that is treated in the next section.

A method for quantifying the canopy structure of large trees is presented in detail by [19]. It involves human beings actually climbing up trees to measure branch positions and lengths using a metric tape. Branches are organized by level in a hierarchy similar to ours, but without the subsequent analytical representation and statistical modeling. CT scanning, when possible, provides more precise measurements than in vivo data collection. Also, CT scanning data allows a complete reconstruction of the canopy, as illustrated in Fig 5 of [5].

4 Statistical models for the tree branching components in the hierarchy

This section presents the statistical models for the 4 components, x, v, ℓ, and n, of the tree branching at each level of the hierarchy, see Tables 1 and 2. All models have a “location parameter”, indexed by z, that depends on the relevant explanatory variables of Table 2. The first three models have, in addition, a shape parameter.

4.1 Modeling the position x of a branch

The relative position x of an offspring branch on its parent branch varies between 0 and 1. We propose a beta regression model, see [20] for x. Recall that the probability density function of the univariate beta distribution depends on two parameters, α, β > 0, and is given by where Γ(⋅) is the standard gamma function. The reparametrization in terms of the dispersion parameter ϕ = α + β > 0 and the expected value μ = α/(α + β) ∈ (0, 1) gives where ϕ > 0. In our beta regression model, the logit link function, allows to express the mean value in terms of explanatory variables where the columns of matrix Z are explanatory variables selected among those of Table 2 as discussed in Section 3.

4.2 Modeling the direction v of a branch

For v ∈ S2, the 3D unit sphere, we use the small circle model of [21] that depends on a unit vector u ∈ S2 and real parameters, θ and τ > 0. Its density is given by where the normalizing constant in the denominator can be expressed in terms of the standard normal cumulative distribution function Φ as

In our application of this model, u is the known direction of the parent of the branch with direction v and the parameters (θ, τ) are unknown. Note that the density 1 only depends on the cosine c = v⊤u in Table 2. The cosines c in Table 2 are sufficient statistics to fit this model. The maximum value of (1), corresponds to c = v⊤u = θ. Explanatory variables Z allow this most likely value of the cosine to depend on explanatory variables such as the branch position x and ancestor characteristics through the link function

We kept our model for the direction of branches simple, by using the identity link. Other link functions, such as a modified logit link function, were examined. They did not result in a noticeably better fit, as measured by the R2 defined in (4).

4.3 Modeling the length ℓ of a branch

The length ℓ of a branch is a positive variable. Its distribution may depend on the position x and the cosine c of the branch and on some ancestor characteristics. This is modeled using a Weibull regression model, see [22]. Recall that the Weibull distribution has density where η is a shape parameter and where the scale parameter σ > 0 may depend explanatory variables such as x, and c. Note that η = 1 gives the exponential distribution. We used the log link for σ: where β is a vector of regression parameter for σ.

4.4 Modeling the number of offspring n of a branch

As the random variable n takes non-negative integer values, we use a Poisson regression model where n is assumed to be a Poisson random variable with expectation λ > 0; see [23]. The parameter λ depends on branch characteristics contained in the vector Zλ, through

4.5 A unified approach to model fitting and model selection

The 3D line segment for a level 4 branch depends on eleven characteristics, (x2, v2, ℓ2, n2, x3, …, ℓ4), see Table 2. The models proposed in this section decompose the joint density for the eleven underlying random variables in terms of the marginal density for x2 times the conditional density of v2 given x2 times the conditional density of ℓ2 given x2 and v2 and so on. In each model, the conditioning variables are candidates explanatory variable that could enter the Z matrix for a particular dependent variable. Consider, for example, the cosines of the level 3 branches. According to the hierarchy of Table 2, the data set for this analysis is {(x2, c2, ℓ2, n2, x3, c3): i indexes level 3 branches}. When fitting model (1) to this data set we assume that given the ancestor variables {(x2, c2, ℓ2, n2, x3)} the cosines {c3} of the level 3 branches are statistically independent. Quadratic functions of the ancestor variables enter in the matrix Z of 3. Our proposal to learn the structure of the three trees considered in this work is to model the eleven characteristics, (x2, v2, ℓ2, n2, x3, …, ℓ4) by selecting, for each one, the most important explanatory variables.

R-packages are available to fit three of the four models proposed in this section. For the beta regression of Section 4.1, we used the R package betareg [24, 25]. To estimate the parameter of the Weibull regression model, we used the R package flexsurv [26] that uses the parametrization presented in Section 4.3. Finally, the Poisson regression was fitted using the R function glm [27]. We developed our own computer code to evaluate the log-likelihood of the Bingham and Mardia model presented in Section 4.2 and used the function nlminb [27] to maximize it.

The stepwise selection of the explanatory variables Z in each of the model is based on the Akaike Information Criterion (AIC) that is evaluated as where k is the number of parameters in the model and is the value of the log-likelihood function at the model parameter estimates.

The stepwise procedure starts with a baseline model where Z only has an intercept; its AIC is evaluated. Candidates models obtained by adding to the Z matrix one variable, allowable according to the hierarchy implied by Table 2, are fitted and their AIC are evaluated. The best candidate variable is the one corresponding to the model with the smallest AIC. If this AIC is smaller than the baseline AIC then this variable is added to the baseline model. The selection procedure is repeated: all variables not already in the model are tried and the AICs of the resulting model are compared. The procedure stops when adding any variable increases the baseline AIC. This first step selects linear explanatory variables. The results of this first step, for the three trees investigated here, are summarized in Table 4. The second step investigates quadratic terms obtained with the variables selected at step 1. A quadratic term is either a product of two variables or a variable squared. If q variables are selected at step 1 there are q(q + 1)/2 possible quadratic variables. This new pool of explanatory variables is investigated in a stepwise manner. The procedure is similar to that for step 1. Variables not already in the model are tried one at a time and the one that gives the smallest AIC is added to the baseline model if this results in a decrease of the AIC. The model fitting algorithm is described further in Algorithm 1 in Appendix A of S1 Appendix.

If the set of candidate explanatory variables is of size q, the total number of possible models is given by The first five integers in this sequence are 1, 3, 13, 95, 1337. The number of possible candidate models increases rapidly as the number of candidate variables increases. More information about this integer sequence can be found in Sloane in [28].

In addition to calculating the AIC for each model, the algorithm also provides the generalized R2 [29], which is defined as where is the value of the log-likelihood function of the full model and is the log-likelihood of the null model, containing only an intercept. This generalized R2 is derived from the likelihood test statistic used to test H0: θ = 0 vs. H1: θ ≠ 0 for some parameter θ [29, 30].

5 Results and discussion

The means and the standard deviations s in Table 3 provide interesting information about the trees’ branching strategy. The size of the data set n• for a given level i for a tree is the total number of offspring branches at the previous level, i − 1, for that tree. For instance the size n• = 146 for the level three branches of Pixie Tree is equal to the number of level 2 branches, 46, times 3.17, the mean for variable n2. Fig 2 gives a biplot representation of the 11 variables and of the three trees to help with the interpretation of the results. The x3 variable point is close to the tree point for Pixie, showing that the largest mean value for x3 is that of Pixie Tree; indeed Pixie Tree has the largest mean values for , where k = 2, 3, 4 represents the branch level in the hierarchical model and is larger than both and for all trees. The three variable points c2, c3, c4 are close to each other in the bottom right quadrant of the biplot and close to the tree point Wonder. This shows a positive association between that tend to decrease with k, also Wonder Tree has the largest mean cosines . For length, the distribution of the tree points relative to the variable points ℓ2, ℓ3 and ℓ4 reflects that Wonder Tree has the largest mean lengths at levels 2 and 4, and shares the largest mean length with Tompa Tree at level 3; also, the mean lengths decrease with k and the smallest mean lengths are for Pixie Tree whose shape is conical and distinct from the rounded form of the two others. The mean values and are negatively associated showing a different branching strategy for Pixie Tree, with , as compared to equal to respectively 22% and 47% for Wonder Tree and Tompa Tree respectively.

Table 3

Summary statistics of Pixie Tree, Wonder Tree and Tompa Tree (in this order, from top to bottom in each cell of the table) for the variables defined in Table 2.

Variable	m¯	s	n •	min	max	Tree
x 2	0.48	0.31	46	2 × 10−3	0.99	Pixie
	0.39	0.22	24	0.07	0.76	Wonder
	0.47	0.25	15	0.09	0.82	Tompa
c 2	0.66	0.25	46	0.05	0.96	Pixie
	0.78	0.14	24	0.45	1.00	Wonder
	0.58	0.19	15	0.27	0.92	Tompa
ℓ 2	0.19	0.13	46	0.02	0.43	Pixie
	0.37	0.18	24	0.15	0.79	Wonder
	0.32	0.18	15	0.06	0.59	Tompa
n 2	3.17	2.70	46	0	9	Pixie
	2.25	3.89	24	0	15	Wonder
	2	2.75	15	0	7	Tompa
x 3	0.59	0.26	146	0.10	1	Pixie
	0.46	0.20	54	0.04	0.82	Wonder
	0.41	0.27	30	0.05	0.80	Tompa
c 3	0.62	0.28	146	-0.49	1.00	Pixie
	0.67	0.25	54	-0.13	0.99	Wonder
	0.58	0.21	30	0.12	1	Tompa
ℓ 3	0.09	0.05	146	0.01	0.34	Pixie
	0.23	0.11	54	0.08	0.79	Wonder
	0.23	0.15	30	0.03	0.54	Tompa
n 3	0.37	1.20	146	0	9	Pixie
	0.50	1.68	54	0	11	Wonder
	0.93	1.36	30	0	4	Tompa
x 4	0.64	0.33	54	0.04	1	Pixie
	0.44	0.16	27	0.10	0.72	Wonder
	0.62	0.07	28	0.51	0.75	Tompa
c 4	0.60	0.30	54	-0.6	0.99	Pixie
	0.67	0.22	27	0.01	0.91	Wonder
	0.59	0.20	28	-0.25	0.92	Tompa
ℓ 4	0.08	0.06	54	0.02	0.39	Pixie
	0.21	0.10	27	0.09	0.47	Wonder
	0.09	0.04	28	0.03	0.21	Tompa

Fig 2

Biplot representation of the 11 variable means and of the three tree species constructed using the mean values reported in Table 3.

In Table 3, the standard deviation s provides a coarse measurement of the variability for each variable at different levels. To investigate whether and to which extent this variability could be explained, the models of Section 4 were fitted using the stepwise model selection procedure presented in Section 4.5. The explanatory variables that were found for each variable are reported in Table 4. A detailed presentation of the models selected for each of the 3 x 11 = 33 dependent variables is given in the Appendices B, C, D, and E. Table 4 presents the main effects selected by our stepwise selection algorithm, for each of the dependent variable. In total, the models for Pixie Tree and Tompa Tree have more explanatory variables (22 each) than those for Wonder Tree (14). The explanatory variable selected most often (i.e., in 11 models out of 24) is ℓ2, emphasizing that the length of the level 2 branch is a key determinant of the structure of the tree branching patterns. It is interesting to note the absence of relationships between offspring and parents variables, such as x and x. Indeed in 36 models fitted to levels 3 and 4 variables the parent variable is selected only 4 times (4/36), that is in 11% of the cases.

Table 4

Results of the stepwise model selection procedure.

Letters P, W and T represent Pixie Tree, Wonder Tree and Tompa Tree, respectively. A green (blue) shaded sub-cell means that the covariate has a positive (negative) estimated coefficient in the fitted model.

		Variable
		x 2	c 2	ℓ 2	n 2	x 3	c 3	ℓ 3	n 3	x 4	c 4
Model response	c 2	P
	ℓ 2	P	P
		W	W
		T	T
	n 2	P	P	P
		W	W
		T		T
	x 3			P
				T
	c 3		P	P
						W
	ℓ 3	P				P	P
					W	W	W
			T	T	T	T
	n 3	P		P	P	P	P	P
								W
				T	T			T
	x 4
							W
			T						T
	c 4		P	P						P
				W
						T	T
	ℓ 4	P		P	P
					W					W	W
		T	T			T	T	T	T	T

From Table 4, it appears that many of the level 2 variables are important predictors for almost all of the responses. Interestingly, this result applies to the level 4 response variables, whereas the level 4 branches are directly connected to the level 3 branches rather than the level 2 branches. This finding is further discussed below. Table 4 supports the expectation that longer branches tend to have more offspring as ℓ is selected as an explanatory variable for n, with a positive coefficient, in 5 of the 6 models for k = 2, 3.

To get a better understanding of the explanatory power of the ancestor variables, the R2 values for each of the models selected was calculated using (4) and is reported in Table 5.

Table 5

The R2 values for the selected modes of Table 4, including interaction effects presented in the supplementary material.

Model	R2 of the	R2 of the	R2 of the
Response	Pixie Tree	Wonder Tree	Tompa Tree
x 2	NA	NA	NA
x 3	0.16	0	0.16
x 4	0	0.29	0.24
c 2	0.15	0	0
c 3	0.08	0.04	0
c 4	0.52	0.18	0.45
ℓ 2	0.79	0.86	0.48
ℓ 3	0.31	0.60	0.83
ℓ 4	0.33	0.67	0.67
n 2	0.89	0.99	0.96
n 3	0.68	0.84	0.82

It is noticeable from Table 5 that certain models perform better than others. Position (x) and orientation (c) are not explained well, with the exceptions of orientation of the level 4 branches for Pixie Tree (R2 = 0.52) and Tompa Tree (R2 = 0.39). The length (ℓ) and number (n) of offspring are explained quite well to very well depending on the level or the tree (0.31 ≤ R2 ≤ 0.99). Specifically, the R2 value of ℓ2 is higher than those of ℓ3 and ℓ4 for Pixie Tree and Wonder Tree, whereas the highest R2 value for ℓ is at level 3 for Tompa Tree. We further notice that the selected model for n2 has a higher R2 value than that for n3, and this for all three trees.

In the biological discussion that follows, we focus on the observed effect of the level in our application of the hierarchical models of tree branching, on the effect of the known differences in traits of Pixie Tree, Wonder Tree and Tompa Tree, and on a possible interaction between these two effects. To begin with, these three trees were chosen from the 15 conifers studied by [5] because the structural complexity of their branching pattern ranged from high (Pixie Tree) to low (Tompa Tree) with Wonder Tree in-between. This was measured with fractal dimension estimates (FD) by [5] [Table 2]: FD = 1.35–2.00 for Pixie Tree; 1.16–1.95 for Wonder Tree; and 1.09–1.47 for Tompa Tree. This difference in structural complexity, as measured by FD, is ‘diffused’, or distributed over our results here. The age of the three trees, or the stage of development that they had reached at the time of their CT scanning, is another important factor to bear in mind when interpreting our results. For example, for a Picea glauca Pixie specimen, the growth rate is known to be less than 2.5 cm per year. With a measured height of 22.1 cm, the ‘equivalent age’ of our specimen (Pixie Tree) may be estimated to be about 9 years. A similar reasoning shows that the estimated ages of Wonder Tree and Tompa Tree are also 9 years. The apparent difference in space occupancy for Pixie Tree and Wonder Tree versus Tompa Tree is reflected by the increasing ratios as a function of FD. It is also appears in the negative versus positive minimum value of the directional cosine c at level 3: min c3 = −0.49 (Pixie Tree), -0.13 (Wonder Tree), and 0.12 (Tompa Tree), while max c3 is equal or close to 1.0 for the three trees (Table 3). It means that offspring branches started to ‘lean down’ as soon as at level 3 in Pixie Tree and Wonder Tree, whereas the Tompa Tree branches show such a pattern only at level 4, min c4 = −0.25. This difference in branching pattern may be related to the small number of level 2 branches for Tompa Tree, n2 = 15, compared to n2 = 46 for Pixie Tree and 24 for Wonder Tree.

Level 2 branch variables are important predictors for the level 4 branch response variables. This may be counter-intuitive because the level 4 branches are directly connected to the level 3 branches, rather than the level 2 branches. Basically, prediction consists in the explanation of the variability contained in data, targeting the spatial distribution and space occupancy of the branches. Here, it means that at the stage of development of the trees, level 2 branches (directly connected to the trunk) were better established than level 3 branches which might show more uncontrolled or noisy variability in their position, orientation, length, and number. It is possible that for the same reason, ℓ4 is more difficult to predict for the three trees, and ℓ2 is better predicted for two of the trees, i.e., Pixie Tree and Wonder Tree.

In closing, the satisfactory model fitting for branch orientation at level 4 for Pixie Tree (R2 = 0.52) and Tompa Tree (R2 = 0.39), but not for Wonder Tree, is a result that departs from the other results showing a greater similarity between Pixie Tree and Wonder Tree, along the gradient of structural complexity of tree branching as measured by FD. It suggests that one or two FDs do not capture all aspects of a tree branching pattern, and the statistical models that we present and apply here provide complementary information.

6 Conclusion

We have shown how to leverage CT scanning data for a tree crown by constructing an analytical representation of the tree branches as a hierarchical set of 3D line segments. The determinants of these line segments are the position, the orientation, the length, and the number of offspring, measured in a hierarchy with three levels. These variables, discussed in Section 3, characterizes the geometry of a tree canopy. If measurable in experiments with multiple trees of large size, they could be explanatory variables in predictive models developed with the same objective as [13, 14]. For instance, branch orientation could be an important predictor and branches orthogonal to the dominant wind direction might have higher epiphyte abundance.

Section 4 has presented statistical tools and distributions for a within tree analysis of the variables in the analytical representation. Results of the statistical modeling for three miniature conifers with different structural complexity of the crown were presented. The four geometric characteristics of the analytical representation defined in Section 3 have the same range, regardless the size and species of the tree. The position x is in the interval (0,1), v is a 3D unit vector, and the length ℓ (relative to that of the trunk) is a positive variable, while the number n of offspring is a non-negative integer. Thus, the models proposed in Subsection 4.1 are widely applicable and could be extended to many other trees.

The different levels of structural complexity of the three tree crowns were corroborated by a between-tree analysis based on the mean values and by the models selected to explain the variation at successive levels of the hierarchy. Section 5 shows important between-tree differences. The models with the poorest explanatory power (median R2 of 8%) are those for branch orientations. In future work, it will be interesting to investigate more complex directional models including predictors such as the direction of the mother branch and that of the trunk.

The three trees considered in Fig 1a, 1c and 1e, are from different species. The statistical analyses presented in this work is conclusive in that it has identified geometric characteristics that explain differences among tree crowns. With several CT scanned trees of the same species, the statistical analysis and modeling would allow a quantification of the within-species variation of the geometric characteristics of crowns. A non-destructive, repeated CT scanning of the same tree can also be envisaged, and would provide supplementary data to investigate the strategy followed by a tree in the growth of its crown.

Last but not least, if several trees were CT scanned for a large number of species, other explanatory variables such as age and physico-chemical properties of the growth medium could be added to the Z-matrices of the models of Section 4. This might lead to a universal modeling applicable to a wide range of tree crowns.

The zipped file DMRSuppMat contains the analytical representations for the three trees condidered in this work.

It also contains an R script, basic_example.R, for fitting length models of Pixie Tree, thus reproducing the results presented in Table 12 in S1 Appendix.

(ZIP)

Click here for additional data file.

(PDF)

Click here for additional data file.

28 Feb 2022

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Esmaiel Jabbari, PhD

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. We note that the grant information you provided in the ‘Funding Information’ and ‘Financial Disclosure’ sections do not match.

When you resubmit, please ensure that you provide the correct grant numbers for the awards you received for your study in the ‘Funding Information’ section.

3. In your Data Availability statement, you have not specified where the minimal data set underlying the results described in your manuscript can be found. PLOS defines a study's minimal data set as the underlying data used to reach the conclusions drawn in the manuscript and any additional data required to replicate the reported study findings in their entirety. All PLOS journals require that the minimal data set be made fully available. For more information about our data policy, please see http://journals.plos.org/plosone/s/data-availability.

Upon re-submitting your revised manuscript, please upload your study’s minimal underlying data set as either Supporting Information files or to a stable, public repository and include the relevant URLs, DOIs, or accession numbers within your revised cover letter. For a list of acceptable repositories, please see http://journals.plos.org/plosone/s/data-availability#loc-recommended-repositories. Any potentially identifying patient information must be fully anonymized.

Important: If there are ethical or legal restrictions to sharing your data publicly, please explain these restrictions in detail. Please see our guidelines for more information on what we consider unacceptable restrictions to publicly sharing data: http://journals.plos.org/plosone/s/data-availability#loc-unacceptable-data-access-restrictions. Note that it is not acceptable for the authors to be the sole named individuals responsible for ensuring data access.

We will update your Data Availability statement to reflect the information you provide in your cover letter

4. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide

5. Your abstract cannot contain citations. Please only include citations in the body text of the manuscript, and ensure that they remain in ascending numerical order on first mention.

6. We note you have included a table to which you do not refer in the text of your manuscript. Please ensure that you refer to Tables 6, 7, ,8 , 9, 10, 11, 12, 13, 14, 15,16, and 17 in your text; if accepted, production will need this reference to link the reader to the Table.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: No

Reviewer #2: Yes

Reviewer #3: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: No

Reviewer #2: Yes

Reviewer #3: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: No

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: To understand this paper one needs to have read Dutilleul et al. (2015). A subset of three trees are used here. CT (computerized tomography) seems to be limited to just ‘miniature’ trees (ca. 20 cm tall and wide). Is there any reason to believe that, ontologically, such models will apply to mature trees? It is difficult to see what this work seeks to achieve apart from applying CT scanning to very small trees. Further, the new paper seems to be a rerun of the earlier search for internal patterns to branching with fractal analysis, one that now looks for statistical correlations.

A main concern is what is meant by a ‘sample’ in this paper. From my reading, the whole tree is scanned to obtain a complete skeletal structure (e.g. Fig. 1a). This not a sample in the statistical sense. In other words, the data analysis does not use repeatedly-taken random samples of branches (and sub-branches) but all the information about each tree is taken together. Since statistical distributions are applied to model the error terms (as part of the different regression fitting), using all data per tree means surely that the individual (derived) values are not independent (see Table 1). They are very highly likely to be spatially autocorrelated. Unless, I have overlooked something, or there is a part missing to the paper, the application of regression and statistical inference cannot be valid. I could find no discussion by the authors of ‘sampling’, ‘independence’ or ‘randomization’. On line 225, ‘sample means’ and ‘sample standard deviations’ are referred to. Examining the data frames in the DMRSuppMat.zip files it would indeed appear that all data are involved after the line-segment conversion. On line 181 it is said that there are ‘eleven underlying random variables’, which left me more confused. Furthermore, how can it be justified that a conditional variable becomes an explanatory one? Much more explanation of the rational is needed. Important issues concerning causality are involved and need addressing.

The main outcome of this analysis is given on lines 305-313: that level 2 branch variables are the most important to tree form. But this will not be a surprise to many plant scientists given that major 1-D segments are being distributed (placed) within a 3-D volume (or, when more constrained, onto a 2-D plane). And that level 4 branch variables depend upon the level 2 ones, but not on those of level 3, is not really “intriguing”, at least to me, rather it is counterintuitive to what we know of the physiology and development of growing trees. Surely this aspect needs to be unravelled and fully investigated before publishing? Perhaps this ‘odd’ result occurs because of a mix of positive and negative correlations between variables which have a spatial auto-correlative component?

In conclusion, reading this paper left me dissatisfied and sceptical about both the method and the results. No particular idea or hypothesis was being properly tested. How does such an analysis help advance our understanding of tree architecture and growth? It is more a CT technical applications report.

Minor points:

1. The Abstract has several language errors which lead to misunderstandings. ‘CT’ needs to be defined. The first sentence is not an appropriate way to begin.

2. The Introduction on page 2 has several syntax errors. On line 12, is ‘diverted’ the right word? And, likewise on line 26, is ‘privilege’? In some sentences there are words missing.

3. Figs 1-3 could be easily combined with a common legend.

4. By ‘beta’ distribution I assume (checking the equation in other texts) that ‘beta-binomial’ is meant. It is not explained why this distribution was used to model position x. The beta-binomial is often employed to cater for over-dispersion. Is that the case here? I am unsure that it is the right error distribution and more care is needed to justify it.

5. The conclusion has at least seven spelling errors.

6. Appendices A-E could readily go into the Supplementary Materials files.

Reviewer #2: The authors used modeling and statistical techniques to investigate the branching pattern of trees. The manuscript and data are well organized and presented. This is an interesting approach to analyze trees structure which showed that their branching patterns have various levels of complexity.

Reviewer #3: In this work, branching pattern of trees has been investigated using a statistical modeling technique. The statistical model was fitted to data derived from CT scans of three tree types. I’d recommend accepting the manuscript after the following minor changes.

-Please discuss how the presented method and derived variables/coefficients for three modeled trees could be expanded to other tree types?

-Please discuss how the results will be changed if the modeling is repeated with more trees of the same type. Are the results reported in this work conclusive for three tree types that were investigated, or the model needs to be fitted to more CT data of the same tree type in the future?

-Please discuss the possibility of developing a universal model that accounts for all parameters including tree age and environmental factors.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

15 Apr 2022

Responses appear in the cover letter and in the Reply to reviewers pdf file.

Submitted filename: Reply_Finalv2.pdf

Click here for additional data file.

27 Jun 2022

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Esmaiel Jabbari, PhD

Academic Editor

PLOS ONE

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #4: (No Response)

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #4: Partly

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #4: No

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #4: No

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #4: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #4: In the MS "Learning how a tree branches out: A statistical modeling approach," the authors analyzed the data from three small trees. I do believe that this effort could be beneficial, but:

a.- Trees are the dominant life form in all the forests, and understanding their architecture is not only crucial for modeling their photosynthetic capacity. Trees are the minimal habitat for canopy organisms such as epiphytes and insects. In this regard, several authors have tried to test how trees' complexity affects epiphyte diversity, see:

Ruiz-Cordova, J. P., V. H. Toledo-Hernández and A. Flores-Palacios. 2014. The effect of substrate abundance in the vertical stratification of Bromeliad epiphytes in a tropical dry forest (Mexico). Flora 209: 375–384.

Flores-Palacios, A. and J. G. García-Franco. 2006. The relationship between tree size and epiphyte richness: testing four different hypotheses. Journal of Biogeography 33: 323–330.

Victoriano-Romero, E., S. Valencia-Díaz, V. H. Toledo-Hernández and A. Flores-Palacios. 2017. Dispersal limitation of Tillandsia species correlates with rain and host structure in a central Mexican tropical dry forest. PLoS ONE 12(2): e0171614.

At least one of these papers measures most of the tree structures.

b.- All trees have a hierarchy of structures; for example, the first trunk division is not independent of the trunk, and so on. So the different structures and their properties must be modeled considering that they are not independent.

c.- The model description is complex (I understand this), making the replicability of the modeling will be limited. Can you provide a simple r-library for the use of your model? Can you give a single number that will describe branching patterns or tree complexity? If yes, how can this number be used for comparing tree species? Will your model help test ecological hypotheses like those in previous articles?

Join Figures 1, 2, and 3.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #4: Yes: A. Flores-Palacios

**********

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

25 Jul 2022

See the Response_to_reviewers pdf file in the submission package

Submitted filename: Response_to_reviewers_R2.pdf

Click here for additional data file.

24 Aug 2022

Learning how a tree branches out: A statistical modeling  approach

PONE-D-21-24626R2

Dear Dr. Rivest,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Esmaiel Jabbari, PhD

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #4: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #4: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #4: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #4: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #4: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #4: Thanks for your answers. I think that the impact of the MS will be better if, in the future, you publish an R library with a friendly user’s manual; in this way, you will win many users.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #4: No

**********

30 Aug 2022

PONE-D-21-24626R2

Learning how a tree branches out: A statistical modeling  approach

Dear Dr. Rivest:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Esmaiel Jabbari

Academic Editor

PLOS ONE