# Transformation and regression statistics of the size-effect method for determining fracture energy and process zone size in quasi-brittle materials

Madura Pathirage, Gilles Pijaudier-Cabot, David Grégoire, Gianluca Cusatis

► To cite this version:

Madura Pathirage, Gilles Pijaudier-Cabot, David Grégoire, Gianluca Cusatis. Transformation and regression statistics of the size-effect method for determining fracture energy and process zone size in quasi-brittle materials. Materials and structures, 2025, 58 (39), 10.1617/s11527-024-02565-x. hal-05044779

HAL Id: hal-05044779

https://hal.science/hal-05044779v1

Submitted on 24 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Transformation and regression statistics of the
size-effect method for determining fracture energy
and process zone size in quasi-brittle materials

Madura Pathirage¹*, Gilles Pijaudier-Cabot³, David Grégoire³,
Gianluca Cusatis²

¹*Gerald May Department of Civil, Construction and Environmental Engineering, University of New Mexico, Albuquerque, 87111, NM, USA.

²Department of Civil and Environmental Engineering, Northwestern University, Evanston, 60208, IL, USA.

³Universite de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LFCR, Anglet, France.

*Corresponding author. E-mail: mpathirage@unm.edu;

Abstract

Accepted manuscript in *Materials and Structures*.

DOI: 10.1617/s11527-024-02565-x.

The final publication is available at:

https://doi.org/10.1617/s11527-024-02565-x.

This paper investigates the regression statistics of the size-effect method to obtain fracture parameters of quasi-brittle materials. The correct nonlinear regression model and assumptions are established and verified using a large dataset of size-effect tests extracted from the literature. The effect of model transformation on the change in error structure is then investigated. Three different transformations are considered, including the one leading to the linear regression recommended by RILEM [38]. The behavior of the nonlinear least squares estimators of the fracture parameters corresponding to the untransformed space, i.e., peak load *P* vs. specimen size *D*, and to each of the three transformations are discussed. Monte Carlo simulations on generated data show that the transformations lead to the violation of model assumptions and to highly skewed error distributions prone to artificial outliers. The paper also shows that the estimator corresponding to the RILEM recommendation is asymptotically biased. The estimators corresponding
to the other transformations are found either asymptotically biased or do not possess the minimum variance property. Finally, simulations show that the least squares point estimates of the unknown fracture parameters differ when a model transformation is used, and that the difference is statically significant. The fitting of the fracture parameters through the size-effect method should only be obtained in the space ($P$ vs. $D$) for which the nonlinear least squares estimator is asymptotically unbiased, mean square consistent, and has minimum variance. The linear regression plot suggested by RILEM should be avoided for the statistical inverse problem of the size-effect method.

**Keywords:** Size-effect law, RILEM, Regression statistics, Asymptotic theory, Model transformation, Error structure

# 1 Introduction

Size-effect is one of the most important consequence of fracture mechanics theory. In quasi-brittle materials, size-effect manifests itself through a reduction in strength when the structural size increases in geometrically similar structures. This effect has been extensively investigated in the literature and was recently adopted in specific design standards [1]. There are two main reasons explaining size-effect [7]. The first cause, considered in this study, is attributed to the stress redistribution due to stable crack propagation and release of stored energy into the fracture front. The second reason is the spatially correlated variability of material properties. This paper will solely focus on the size-effect produced due to crack formation and energy release. The size of the fracture process zone in quasi-brittle materials is non-negligible as compared to the structural size but a complete nonlinear theoretical fracture mechanics framework does not exist. Therefore, an exact form of the size-effect law is not known. However, an asymptotic scaling law was discovered by Bažant in the 1980s [4]. This approximate and simple analytical formula known as Bažant's size-effect law or alternatively as the energetical size-effect law is sufficient for most practical purposes. It involves two key material parameters: the fracture energy $G_f$ characterizing a stable crack propagation and the effective length of the fracture process zone $c_f$ which represents the length of the equivalent linear elastic crack that gives the same unloading compliance as the actual crack in a sufficiently large specimen at peak load.

In structures of high sensitivity to fracture and size-effect, an accurate determination of $G_f$ and $c_f$ is necessary. Scatter in test results is inevitable due to the highly heterogeneous nature of many quasi-brittle materials. However any additional bias in the determination of fracture parameters should be avoided, especially when it comes to the statistical method employed in the inverse analysis. In an effort to provide a standardized test procedure, one of the RILEM recommendations on concrete fracture, published in this same journal in 1990 [38], proposes to use an algebraically rearranged formulation of the size-effect law to a linear regression plot to determine fracture characteristics.

As it will be shown in this paper, this seemingly harmless transformation leads to a severe statistical bias in the least squares point estimates of the fracture parameters
due to the distortion of the error structure. A few other transformations resulting in linear and nonlinear regressions accounting or not for weights have been used and compared in the past [6, 7, 36, 43], but none of these studies rely on rigorous nonlinear regression statistics, rather on an inadequate adaptation of linear regression metrics. Eventually, the choice of the statistical method is left to the researcher's judgment who in most cases follows the RILEM recommendation.

Altering the regression model and consequently the statistical method can often lead to very different conclusions. The difference of results obtained using multiple regression methods can well be statistically significant, and therefore cannot be neglected. Moreover, there exist certain fundamental theoretical requirements in regression statistics that must be respected when constructing a regression model.

This paper deals exclusively with statistics and shows the importance of using the correct statistical regression method for the identification of fracture parameters using the size-effect formula. The only appropriate nonlinear regression model is described. The model assumptions are clearly stated and verified based on a large database of size-effect data collected from the available literature for 16 different types of quasi-brittle materials comprising a total of 71 sets of size-effect tests, for a total of 905 individual fracture tests. The effect of algebraic rearrangement of the size-effect law, i.e., model transformation and the induced change in error structure, is demonstrated analytically and on simulated data through Monte-Carlo simulations. Asymptotic inference for the model parameters are then discussed, ruling out all the regression models providing asymptotically biased or less efficient estimators, which therefore should not be employed for the statistical inverse problem of the size-effect method.

# 2 Regression model and assumptions

## 2.1 Bažant's size-effect law

The problem of scaling for quasi-brittle materials is an intricate and challenging task. It simplifies significantly in the case of large sizes as the problem can be tackled approximately by the powerful method of asymptotic matching on the basis of equivalent linear elastic fracture mechanics. A Taylor series expansion truncated to the first order eventually provides Bažant's size-effect law

$$ \sigma_N = c_N \frac{P}{bD} = \frac{\sqrt{EG_f}}{\sqrt{g'_0 c_f + g_0 D}} \quad (1) $$

where $\sigma_N$ is the nominal stress at failure, $c_N$ is a dimensionless coefficient introduced for convenience, and $P$ is the ultimate load. $b$ and $D$ are the width and characteristic dimension of the considered geometry, respectively. $E$ is the elastic modulus, $g(a_0/D) = g_0$ where $g$ is the dimensionless energy release rate function evaluated at the initial crack length $a_0$, and $g'(a_0/D) = g'_0$ is its derivative evaluated at the same point.

In the size-effect method, the objective is to quantify the effect of the structural size $D$ on the nominal strength $\sigma_N$ in geometrically similar structures. For such structures, the coefficients $b, c_N, g_0$ and $g'_0$ have fixed values for a given geometry. In addition, the
elastic modulus is size-independent and constant. D is the only fixed regressor, in the sense that it can take a variety of values but it is not a random variable. One can safely assume that any random variation in the measurement of the regressor variable D is much smaller than the range of that variable when applying the size-effect method. Let's also assume that there is no significant errors in the measurements of b, c_N, and a_0. Note that if errors in measurements of either the regressor or the other geometrical dimensions are considered, one would need to develop an errors-in-variables model [2]. This would however only add complexity to the regression analysis and would probably not improve it by much. At last, the physical unknown parameters to be identified are G_f and c_f.

## 2.2 Regression model

One can first rearrange Equation 1 and write

$$P = \frac{bD}{c_N} \frac{\sqrt{EG_f}}{\sqrt{g'_0 c_f + g_0 D}}. \quad (2)$$

In regression statistics notation, the nonlinear regression model writes as follows

$$y = f(x, \theta) + \varepsilon \quad (3)$$

where y = P is the response variable and x = D is the regressor variable. θ = [θ_1 θ_2]^t = [G_f c_f]^t denotes the 2 × 1 unknown parameter vector. An additive error term ε is considered and represents by construction the discrepancy between the response variable and the expectation function f which writes as

$$f(x, \theta) = \frac{c_1 \sqrt{\theta_1 x}}{\sqrt{c_2 \theta_2 + x}} \quad (4)$$

where c_1 = b√E/(c_N√g_0) and c_2 = g'_0/g_0. Note that the regression model is nonlinear since the expectation function is nonlinear in the unknown parameters.

## 2.3 Model assumptions

There are seven assumptions related to this regression model.

(i) The expectation function is correct. This assumption is not unreasonable as the size-effect law, although approximate by nature, has been theoretically derived in several different ways and has been extensively validated against experimental data.

(ii) The response is the expectation function to which a disturbance term is added. This is a common assumption which allows the probability density function of the random variable y to be simply calculated from the probability density function of the random variable ε describing the observed errors.

(iii) The error term is independent of the expectation function. Again, theoretical and experimental evidences support the fact that there are no major omitted variable
bias in the size-effect law. Therefore, there is no apparent reason to account for
correlation between the error term and the regression function.

(iv) Each error has zero mean. This assumption is mainly a simplifying one, for reduc-
ing the number of model parameters. Indeed, a disturbance with non-zero mean
common to all observations can be accounted for by introducing an additional term
in the expectation function which might depend or not on the regressor variable.
As an example, a systematic bias in the errors is caused by the effect of gravity, viz.
the self-weight of the tested specimens, if not properly accounted for [33, 38]. Quan-
titatively, this assumption implies that each error is a centered random variable,
i.e., with null expectation $E[\varepsilon] = 0$ and therefore $E[y] = f(x, \theta)$.

(v) The errors are independently distributed. This reflects the fact that the errors obtained by varying the size of specimens are independent of one another.

There are two additional assumptions, which are potentially subject to discussions.
They will be tentatively entertained and their validity will be checked by examining
residuals using a large data set collected from the literature (see Section 3). These
assumptions are as follows.

(vi) The errors have equal variances. Let us assume that the variance of the error exist and let's denote it $\text{Var}[\varepsilon] = \sigma^2$ where $\sigma$ is the standard deviation and is identical for all the errors.

(vii) The errors have a normal distribution. The normal distribution occurs frequently in the natural world when collecting data without further manipulation or transformation. The reason is that experimental scatter results in general from a combination of various errors from several sources with no dominating source. In that case, the celebrated central limit theorem ensures that the error distribution will tend more and more to the normal distribution as the number of sources of errors increases, no matter what the probability distribution of the separate errors may be. This is the reason why Equation 1 was rearranged to Equation 2 in the first place since it is the ultimate load that is measured experimentally, not the nominal stress at failure.

Note that the normality assumption is in reality not required because the least squares estimation method introduced next does not necessitate errors to follow a normal distribution. The same comment is valid for the asymptotic analysis presented in Section 5.2. Nevertheless, this assumption is convenient because it allows one to generate simulated data from a known parametric distribution. In that case, the least squares estimator becomes identical to the maximum likelihood estimator.

## 2.4 Least squares estimation

The most widely used statistical method to obtain the estimate of the error term is the unweighted least squares method, which is employed here. Let us denote by $n_i$ the number of replications of the response $y$ for each value of $x_i$. For a given $x_i, y_{ij}$ are the observations for $j = 1, ..., n_i$. The total number of observations is $N = \sum_{i=1}^{m} n_i$ where $m$ is the number of different values that the regressor takes. Equation 3 can therefore
be rewritten to account for replications as follows

$$y_{ij} = f(x_i, \boldsymbol{\theta}) + \varepsilon_{ij} \quad (5)$$

Without recalling the likelihood function and assuming a normal distribution of errors for a data set of N points including replicates, one can write the least squares function S(θ) as

$$S(\boldsymbol{\theta}) = \sum_{i=1}^{m} \sum_{j=1}^{n_i} [y_{ij} - f(x_i, \boldsymbol{\theta})]^2 \quad (6)$$

Any $\hat{\boldsymbol{\theta}} = [\hat{\theta}_1 \ \hat{\theta}_2]^t$ that minimizes the least squares function is called a least squares estimator of $\boldsymbol{\theta}$. $\hat{y} = f(x, \hat{\boldsymbol{\theta}})$ denotes the predicted or fitted value of y for a given x. The residuals are therefore defined as $\hat{\varepsilon}_{ij} = y_{ij} - \hat{y}_i$ where $\hat{y}_i = f(x_i, \hat{\boldsymbol{\theta}})$ for $i = 1, ..., m$ and $j = 1, ..., n_i$.

# 3 Residual diagnostics

## 3.1 Alternative regression model for residual diagnostics

A residual diagnostic is typically performed to verify model assumptions (vi) and (vii). To prove the generality of the results shown in this paper, a large set of experimental data on various quasi-brittle materials was gathered from the literature. This set comprises 71 size-effect tests on different geometries and sizes of specimens made of regular concrete [6, 8, 11, 18, 19, 21, 23, 37], dam concrete [22], self-compacting and consolidating concrete [16, 26], steel fiber reinforced self-compacting concrete [20], mortar [6, 29], hardened cement paste [27], limestone [3, 35], shale [30, 35], sandstone [13, 45], marble and granite [17], composite [31, 39, 40], bone [28], acrylonitrile butadiene styrene [34], aluminium [5], and closed-cell polymeric foam [9, 32]. Fortunately, the tedious task of extracting all the information required to use the size-effect method (Equation 2) can be avoided. Indeed, the expectation function can be rewritten in the following form

$$f'(x, \boldsymbol{\theta}') = \frac{\theta_1' x}{\sqrt{\theta_2' + x}} \quad (7)$$

for which only peak loads and specimen sizes are necessary. The new unknown parameter vector $\boldsymbol{\theta}' = [\theta_1' \ \theta_2']^t$ is chosen to incorporate the different geometrical features. The least squares estimator of $\boldsymbol{\theta}'$ is not of interest since $\hat{y} = f(x, \hat{\boldsymbol{\theta}}) = f'(x, \hat{\boldsymbol{\theta}})$; only the residuals are sought. For each size-effect test, the least squares function defined by Equations 6 and 7 was minimized numerically using a direct search algorithm. Careful attention was given to the numerical convergence of the least squares estimator, and therefore of the residuals, for each set of tests.

The diagnostic of residuals in terms of variance and distribution can be performed at two levels: (i) by using residual plots and normal quantile-quantile plots, and (ii) by using test statistics. In general, a statistical analysis of data is typically performed at
the first level only through a visual assessment of these plots [2]. If there is evidence of
a trend, or if the visual assessment is not conclusive due to various reasons including
a limited number of data points, one can employ hypothesis testing. In the next
sections, test statistics are also used to ensure objectivity and have an additional level
of confidence to support the conclusions of this study.

## 3.2 Residual plots and homogeneity of variance

The diagnostic for residual variance constancy can be verified at different levels. One approach could consist in computing the standard deviation for each size since size-effect tests almost always contain replicates for a given regressor value. While the number of replications strictly allows a calculation of an estimate of the variance, it generally leads to inaccurate variance estimates. A more efficient method consists in plotting the data with respect to the regressor variable and visually assess if the spread of the data tends to systematically increase or decrease with respect to the regressor. Figures 1(a) and (d) show examples of two different sets of size-effect data corresponding to end notch flexure specimens made of twill composite [40] and eccentric-compression specimens made of regular concrete [6], respectively. The dashed line indicate the fitted mean response. A better way to visualize the residual heterogeneity are residual plots. In these plots, the residuals are plotted against the regressor variable (Figures 1(b) and (e)) or against the predicted values (Figures 1(c) and (f)) to account for possibility that the predicted values are not a linear combination of the regressor. The residuals can be considered homogeneous if the plots roughly show a horizontal band of points, which is for example the case for the tests on composite (Figures 1(b) and (c)). Unsatisfactory residual plots to support the constant error variance assumption are the ones exhibiting funnel-shaped bands often widening for increasing regressor or predicted values. Figures 1(e) and (f) clearly display the heterogeneous nature of the error estimates. Unfortunately, the trends of the bands of residual points were not always clear for the 71 datasets analyzed, due for instance to the limited number of specimens tested or the presence of outliers. To ensure objectivity, test statistics were performed to reveal any serious violations of assumption (vi). There are several possible tests that could be used to assess variance homogeneity [42]. Here, a Brown-Forsythe test statistic was performed on each set of data to reject or not the constant variance assumption (see Appendix 7.1 for more details). This test has the advantage to be non-parametric as it does not require the error terms to originate from a normal distribution, which is assumption (vii) not yet verified. The objective is to perform the following hypothesis test

$H_0$: the variance is constant.

$H_1$: the variance is not constant.

Assumption (vi) is supported if one fails to reject the null hypothesis $H_0$. The power of this test statistic was first assessed for the size-effect regression model. While a minimum sample size for a test statistic is not required, its power, i.e., the probability that the test detects some effect when there is one, could be too low if the sample size is too small. To obtain a rough idea of this minimum sample size, Monte Carlo
<figure><img src="image_2.png" /><figcaption>**Fig. 1:** Fitting and residual plots for two sets of size-effect data: twill composite specimens under end notch flexure [40] (a) experimental data and fitting, (b) residuals vs. regressor, (c) residuals vs. predicted values; regular concrete under eccentric-compression [6] (d) experimental data and fitting, (e) residuals vs. regressor, (f) residuals vs. predicted values.</figcaption></figure>

simulations were performed on pseudo-data generated using the size-effect law assuming a normal distribution for the error term, and with predefined fracture parameters (see Section 3.4 for more details on the generation of pseudo-data). Different numbers of replicates for a given regressor value, from 3 to 250, with 3 and 4 different regressor values, and different levels of variance are considered. The Brown-Forsythe test was found powerful enough for sets of data of size 9 and above. The alpha level, i.e., the probability of rejecting the null hypothesis when the null hypothesis is true, was taken as 0.1 rather than the traditional 0.05 to better ensure that the test has sufficient power, as suggested by others [44]. In the database, there are 6 datasets with only 8 samples for which the test statistic was still performed together with a visual inspection of the residual plots.

The results of the test statistic on the 71 datasets reveal three tendencies: (1) there are 61 sets with constant variance as the Brown-Forsythe test failed to reject $H_0$; (2) there are 8 sets with increasing variance, and (3) there are 2 sets with unclear trends. For the last two tendencies, $H_0$ was rejected and the residual plots were used to visually assess the behavior of variance.

Based on these results, it appears that assumption (vi) is appropriate and reflects the real behavior of variance in size-effect testing results. In other words, improving the
regression model (Equation 3) by systematically incorporating a model for variance does not seem justified as the vast majority of the size-effect datasets was characterized by a homogeneous variance regardless of material type, specimen geometry, or experimental test configuration. If there is evidence that the variance for instance increases as the regressor value increases, one could use a model for variance as a function of the regressor, with however additional model parameters to be fitted together with the unknown parameters in the expectation function.

## 3.3 Normal quantile-quantile plots and normality

Assumption (vii) related to the normality of errors can be verified through quantile-quantile plots, i.e., plots of the quantiles of the estimates of the errors versus the theoretical quantile values drawn from a normal distribution. These plots are of course only relevant to the datasets with constant variance. One visually assesses departures from a straight line suggesting that the normality assumption is inappropriate. Figure 2(a) shows an example of error estimates following a normal distribution, corresponding to the size-effect tests conducted on mortar under three-point bending [29] for a total of 20 samples. Quantile-quantile plots are typically sufficient to diagnose deviations from normality. As a matter of fact, statistical tests such as the Anderson-Darling

<figure><img src="image_3.png" /><figcaption>**Fig. 2:** Quantile-quantile plots: (a) normal distribution, mortar under three-point bending [29], (b) non-normal distribution, concrete with notch-to-depth ratio of 0.3 under three-point bending [11], (c) non-normal distribution with very small sample size, acrylonitrile butadiene styrene under three-point bending, orientation C [34].</figcaption></figure>

test statistic [42] are generally not powerful enough to detect small or moderate deviations from normality when dealing with small sample sizes characteristic of size-effect datasets. In order to confirm this, the power of the Anderson-Darling test statistic was evaluated through Monte Carlo simulations, similarly to the approach described in Section 3.2, with an alpha level of 0.1. The formal hypothesis test is written as

$H_0$: the errors follow a normal distribution.

$H_1$: the errors do not follow a normal distribution.
The simulations confirmed that the Anderson-Darling test is not powerful enough for
any of the sample sizes in the 61 datasets with constant variance considered. In other
words, the failure to reject the null hypothesis could well be a false negative result.
Nevertheless, the statistic is still meaningful if the null hypothesis $H_0$ is rejected
regardless of the sample size since it would be an indication of a highly non-normal
error distribution. Therefore, in order to complement the visual inspection of the
quantile-quantile plots, the Anderson-Darling test statistic was applied to all the 61
datasets (see Appendix 7.2 for more details on the test statistic).

$H_0$ was rejected for 5 datasets. Figure 2(b) shows an example of deviation from the normal distribution for a set of size-effect tests performed on concrete with notch-to-depth ratio of 0.3 under three-point bending with 16 samples [11]. Another example of violation of assumption (vii) is illustrated in Figure 2(c) with a very small sample size, only 8 data points, of a size-effect dataset on acrylonitrile butadiene styrene under three-point bending with a certain “C” orientation (see reference [34] for more details). The visual inspection of the quantile-quantile plots for the remaining 56 datasets did not reveal any serious deviation from normality.

It seems therefore reasonable to assume that the error term in the regression model
is normally distributed, i.e., assumption (vii) is retained. Once again, this assump-
tion is not needed to obtain point estimates of the parameters through least squares
minimization but it is a convenient assumption to generate pseudo-data and perform
numerical simulations as it described in the next sections of the paper.

3.4 Generation of simulated data

In order to draw more general conclusions on the effects of transformation of the size-
effect law on the parameter estimates, one cannot use a single set of data. In this
work, generating pseudo-data are used for two main reasons: (i) the true parameters
are known, and therefore one can compare the least-squares estimates with these
true values, and (ii) Monte Carlo simulations can be performed when a parametric
statistical distribution of the residuals is known.

In a typical simulation study, it is assumed that both the form of the deterministic
component of the model $E[y]$ and the nature of the stochastic component $\varepsilon$ are known.
Data are generated using a set of predetermined values of the parameters $\theta$, allowing
only the values of $\varepsilon$ to change randomly from set to set. By this means, many sets of
simulated data are produced, each of which provides a set of least squares estimates of
the parameters which may then be examined for their bias, variance, and distributional
properties. The simulated data were here constructed using the specimen geometry
and test configuration listed in Grégoire et al. [21]. More specifically, a three-point
bending configuration with geometrically similar specimens was considered with the
following parameters: span-to-depth ratio $s/D = 2.5$, out-of-plane thickness $b = 50$
mm, dimensionless energy release rate and derivative $g_0 = 2.963$ and $g'_0 = 18.948$,
respectively, elastic modulus $E = 37$ GPa, and coefficient $c_N = 3s/2D = 3.75$. The
values of the constants $c_1$ and $c_2$ defined in the expectation function in Equation 4
were then computed. The number of replications $n_i = 6$ for $i = 1, ..., 4$ was assumed for
each specimen size. Four values of the specimen size or regressor $D = x_i$ were chosen
($m = 4$), where $x_1 = 50$ mm, $x_2 = 100$ mm, $x_3 = 200$ mm, $x_4 = 400$ mm, for a total
of $N = \sum_{i=1}^{m} n_i = 24$ data points. Note here that a geometric progression (1:2:4:8) was respected without referring to a proper design of experiments. It is important to emphasize that the use of optimal sampling theory would not change the conclusions listed in the paper as the structural size sampling does not affect the variance of the data. From the variance assumption perspective only, it was shown that at least 3 replicates and 3 different structural sizes are needed for a Brown-Forsythe test to be done with an alpha level of 0.1. This requirement is here respected. The simulated data were generated through the following formula

$$y_{ij} = f(x_i, \boldsymbol{\theta}^0) + \varepsilon_{ij} \quad (8)$$

for $i = 1, \dots, 4$ and $j = 1, \dots, 6$ where the true fracture parameters were set to be $\theta_1^0 = 40 \text{ N m}^{-1}$ and $\theta_2^0 = 15 \text{ mm}$. In accordance with Sections 3.2 and 3.3, the errors $\varepsilon_{ij}$ were assumed to be independent and identically distributed (i.i.d.) following a normal distribution with mean 0 and variance $\text{Var}[\varepsilon] = 0.0625$, realistic with respect to the scatter typically encountered in concrete testing for instance. The pseudo-data generated in this subsection will be utilized throughout the rest of the paper. Let us keep in mind that pseudo-data are used mainly because the true parameters are known and Monte Carlo simulations can be performed. The conclusions listed throughout the paper remain of course valid for real datasets.

# 4 Transformation and change in error structure

The objective of this section is to demonstrate that the error structure is changed significantly when the regression model is manipulated algebraically, which eventually violates certain assumptions established previously.

## 4.1 Untransformed model P vs. D

The untransformed model (P vs. D) with a normally distributed error $\varepsilon$, with $E[\varepsilon] = 0$ and $\text{Var}[\varepsilon] = \sigma^2$ writes as

$$y = E[y] + \varepsilon \quad (9)$$

where $y$ represents the peak load P, $x$ designates the specimen size D. The expectation function is defined as

$$E[y] = \frac{c_1 \sqrt{\theta_1} x}{\sqrt{c_2 \theta_2 + x}} \quad (10)$$

The least squares estimators $\hat{\theta}_1$ and $\hat{\theta}_2$ can be obtained by minimizing $S(\boldsymbol{\theta})$ defined as

$$S(\boldsymbol{\theta}) = \sum_{i=1}^{m} \sum_{j=1}^{n_i} \left( y_{ij} - \frac{c_1 \sqrt{\theta_1} x_i}{\sqrt{c_2 \theta_2 + x_i}} \right)^2 \quad (11)$$

Let us see how various transformations routinely used in the literature affect the error structure of the untransformed regression model.
## 4.2 Transformed model $\sigma_N$ vs. $D$

The fracture parameters are sometimes identified by fitting the size-effect law in the space ($\sigma_N$ vs. $D$). This necessitates a transformation of load $y = P$ to nominal stress $y^{(1)} = \sigma_N = (c_3/x)y$ where $c_3 = c_N/b$. The transformed regression model in the space ($\sigma_N$ vs. $D$) writes as

$$y^{(1)} = E[y^{(1)}] + \varepsilon^{(1)} \quad (12)$$

where

$$E[y^{(1)}] = \frac{c_3 c_1 \sqrt{\theta_1}}{\sqrt{c_2 \theta_2 + x}} \quad (13)$$

and

$$\varepsilon^{(1)} = \frac{c_3}{x} \varepsilon \quad (14)$$

One can observe that the error structure is changed by the transformation. Indeed, while the correct mean is conserved $E[\varepsilon^{(1)}] = 0$, the variance $Var[\varepsilon^{(1)}]$ now varies with $x$. Therefore, a well-defined least squares function must account for such a dependency of the error term on the regressor. In other words, the correct least squares estimator of the parameters in the transformed space ($\sigma_N$ vs. $D$) is defined by minimizing the following weighted least squares function identical to Equation 11

$$S(\theta) = \sum_{i=1}^{m} \sum_{j=1}^{n_i} \left(\frac{x_i}{c_3}\right)^2 \left(y_{ij}^{(1)} - \frac{c_3 c_1 \sqrt{\theta_1}}{\sqrt{c_2 \theta_2 + x_i}}\right)^2 \quad (15)$$

where the weights are proportional to $x_i^2$, i.e., are increasing as the tested specimen sizes increase.

## 4.3 Transformed model $\ln(\sigma_N)$ vs. $\ln(D)$

Another space in which the fracture parameters are sometimes identified is the $(\ln(\sigma_N), \ln(D))$ space. By taking the Napierian logarithm of $\sigma_N$, the transformed model writes as

$$\begin{aligned} y^{(2)} &= \ln(y^{(1)}) \\ &= \ln(E[y^{(1)}] + \varepsilon^{(1)}) \\ &= \ln(E[y^{(1)}]) + \ln\left(1 + \frac{\varepsilon^{(1)}}{E[y^{(1)}]}\right) \\ &= E[y^{(2)}] + \varepsilon^{(2)} \end{aligned} \quad (16)$$

where

$$E[y^{(2)}] = \ln(E[y^{(1)}]) \quad (17)$$

and

$$\varepsilon^{(2)} \approx \frac{\varepsilon^{(1)}}{E[y^{(1)}]} \quad (18)$$
for small enough errors $\varepsilon$ and using $\ln(1+u) \approx u$ for small enough $u$. One can note that $E[\varepsilon^{(2)}] \approx 0$ but $\mathrm{Var}[\varepsilon^{(2)}]$ varies with $E[y^{(1)}]$ and $x$. In this case, the correct weighted least squares function to consider, assuming small errors, write as

$$S(\theta) \approx \sum_{i=1}^{m} \sum_{j=1}^{n_i} \left( \frac{c_1 \sqrt{\theta_1} x_i}{\sqrt{c_2 \theta_2 + x_i}} \right)^2 \left[ \ln(y_{ij}^{(1)}) - \ln\left(\frac{c_3 c_1 \sqrt{\theta_1}}{\sqrt{c_2 \theta_2 + x_i}}\right) \right]^2 \quad (19)$$

This time, the weights depend not only on the regressor but on the expectation function, viz. the unknown fracture parameters.

## 4.4 Transformed model $1/\sigma_N^2$ vs. D

This transformation is the one recommended by RILEM [38] and by far the most used transformation to identify fracture parameters using the size-effect law. The transformed model in the space $(1/\sigma_N^2 \text{ vs. } D)$ writes as

$$\begin{aligned} y^{(3)} &= \frac{1}{(y^{(1)})^2} \\ &= \frac{1}{(E[y^{(1)}] + \varepsilon^{(1)})^2} \quad (20) \\ &= \frac{1}{(E[y^{(1)}])^2} \left( 1 + \frac{\varepsilon^{(1)}}{E[y^{(1)}]} \right)^{-2} \end{aligned}$$

in which case the error term became multiplicative. For small enough errors $\varepsilon$ and using $(1+u)^{-2} \approx 1-2u$ for small enough $u$, one can write

$$y^{(3)} = E[y^{(3)}] + \varepsilon^{(3)} \quad (21)$$

where

$$E[y^{(3)}] = \frac{1}{(E[y^{(1)}])^2} \quad (22)$$

and

$$\varepsilon^{(3)} \approx -\frac{2\varepsilon^{(1)}}{(E[y^{(1)}])^3} \quad (23)$$

For this transformation, $E[\varepsilon^{(3)}] \approx 0$ but $\mathrm{Var}[\varepsilon^{(3)}]$ depends on the unknown parameters and on the regressor. The correct weighted least squares function in the space $(1/\sigma_N^2 \text{ vs. } D)$ for small enough errors therefore read as

$$S(\theta) \approx \sum_{i=1}^{m} \sum_{j=1}^{n_i} \left( \frac{c_3 c_1 \sqrt{\theta_1}}{\sqrt{c_2 \theta_2 + x_i}} \right)^6 \left( \frac{x_i}{2c_3} \right)^2 \left[ \frac{1}{(y_{ij}^{(1)})^2} - \left( \frac{c_2 \theta_2}{(c_3 c_1)^2 \theta_1} + \frac{1}{(c_3 c_1)^2 \theta_1} x_i \right) \right]^2 \quad (24)$$

where the weights depend on the regressor and the parameters.
4.5 Distributional properties of errors

To illustrate the effects of the three transformations, the simulated data in Section 3.4 were used. Figure 3(a) shows one realization of the simulated data fitted by the model in the untransformed space. Figures 3(b-d) show how the three transformations described earlier affect the simulated data and the model response. For all the figures, the solid dots are the generated pseudo-data, whereas the dashed lines represent the regression model response using the fitted parameters $\hat{\theta}$ in the space ($P$ vs. $D$), i.e., by minimizing the least squares function listed in Equation 11. Note that using the true parameters $\theta^0$ to obtain the regression model response is adequate only if the total sample size $N$ is large enough. In that case, $\hat{\theta}$ converges in a probabilistic sense to $\theta^0$ (the asymptotic behavior will be analyzed in-depth in the next section). One can observe that the error structure is severely affected for the three transformations. The magnitude of the errors associated with transformations ($\sigma_N$ vs. $D$) and ($\ln(\sigma_N)$ vs. $\ln(D)$) decreases for increasing regressor values. This is due to the division by $x$ of the untransformed model error term (see Equation 14). In other words, the weights that should be included to correctly obtain the estimators of the parameters by minimizing Equation 15 have a visual justification: one must assign more weights for less scattered data. On the other hand, the RILEM transformation (Figure 3(d)) seems to preserve the homogeneity of variance property but artificial outliers are generated.

Quantitatively, one can characterize the distributional properties of the errors through the third and fourth moments about the mean, i.e., skewness and kurtosis, respectively. Skewness is a measure of the asymmetry of a distribution around the mean, which is zero for the errors in spaces ($P$ vs. $D$) and ($\sigma_N$ vs. $D$), and approximately zero in spaces ($\ln(\sigma_N)$ vs. $\ln(D)$) and ($1/\sigma_N^2$ vs. $D$) for small enough errors. Positive skewness means that the error spread out more to the right of the mean than to the left. On the other hand, kurtosis characterizes the presence of outliers in a distribution. Since the comparison between distributional properties is being made with respect to the error distribution in the untransformed space, which is a normal distribution with zero mean and constant variance, it is more relevant to measure the excess kurtosis (= kurtosis - 3). The skewness and excess kurtosis of the estimates of $\varepsilon$, $\varepsilon^{(1)}$, $\varepsilon^{(2)}$, and $\varepsilon^{(3)}$ was evaluated through Monte Carlo simulations with 5,000 different realizations, with the same $\theta^0$ but with $\varepsilon$ that were i.i.d. following a normal distribution with mean 0 and constant variance.

Results show that about in the untransformed space, only 2% of the error distributions were found highly skewed, i.e., when the skewness values are greater than 1 or less than -1 [10]. In space ($\sigma_N$ vs. $D$), about 32% of the error distributions were highly skewed, with 16% positive and 16% negative skewness values. In space ($\ln(\sigma_N)$ vs. $\ln(D)$), about 31% of the error distributions were found highly skewed, with about 27% of them having negative skewness values. Last but not least, in space ($1/\sigma_N^2$ vs. $D$), about 41% of the 5,000 simulations were characterized by highly skewed error distributions, all with positive skewness values. About 1.5%, 45%, 33%, and 28% of the error distributions were characterized by a high kurtosis, i.e. when the excess kurtosis is larger than 2, in the untransformed, ($\sigma_N$ vs. $D$), ($\ln(\sigma_N)$ vs. $\ln(D)$), and ($1/\sigma_N^2$ vs. $D$) spaces, respectively.
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3:** Transformation of the size-effect law: (a) *P* vs. *D*, (b) *σ*<sub>*N*</sub> vs. *D*, (c) ln(*σ*<sub>*N*</sub>) vs. ln(*D*), and (d) 1/*σ*<sub>*N*</sub><sup>2</sup> vs. *D*. The solid dots are the generated pseudo-data, whereas the dashed lines represents the regression model response using the fitted parameters *θ* on the space (*P* vs. *D*).</figcaption></figure>

In addition, a Brown-Forsythe test statistic for homogeneity of variance was per-
formed on each of the 5,000 simulations with a significance level of 0.1. There were
about 8.5%, 96%, 80%, and 12% of the simulations for which the null hypothesis of
constant variance was rejected. The percentage for the untransformed case simply
reflects the amount of Type I error, i.e. when a true null hypothesis is rejected, defined
by the alpha level of 0.1 chosen. A greater number of Monte Carlo simulations would
exhibit a converging trend towards 10%. Simulations not reported here showed a slow
convergence in the studied problem. The vast majority of the error distributions were
characterized by a heterogeneous variance in spaces $(\sigma_N \text{ vs. } D)$ and $(\ln(\sigma_N) \text{ vs. } \ln(D))$,
illustrated on one example in Figures 3(b) and (c). The percentage related to space
$(1/\sigma_N^2 \text{ vs. } D)$ confirms the fact that variance is not affected by the linear transformation.
Similarly, an Anderson-Darling test statistic was performed for each simulation
with an alpha level of 0.05 to assess the normality of errors. Results showed that 4.5%,
65%, 42%, and 40% of the tests rejected the null hypothesis of normality.

To summarize, each of the transformations $(\sigma_N \text{ vs. } D)$, $(\ln(\sigma_N) \text{ vs. } \ln(D))$, and $(1/\sigma_N^2 \text{ vs. } D)$ violates assumption (iii) that states that the error term should be independent of the expectation function, and assumption (v) which states that the errors should be independently distributed. More importantly, all the three transformations lead to highly skewed error distributions with a probability of about 30% to obtain a significant number of outliers. The constant variance property is lost except for the case $(1/\sigma_N^2 \text{ vs. } D)$. The normality property is lost with a probability of about 50% for all the transformations. It is important to mention that the variance and normality conditions might in some cases be statistically significant for the transformed forms of the experimental data. However, a correct use of a regression model requires all the 7 assumptions to be respected.

In regression statistics, there are in general three main reasons that motivate a
model transformation [12]: (1) to achieve linearity, (2) to obtain errors approximately
normally distributed, and (3) to achieve a constant error variance. While the RILEM
transformation (1/σ<sub>N</sub><sup>2</sup> vs. D) is motivated by the first reason and does not violate the
third reason, the error distribution does not follow a normal distribution, and is prone
to artificial outliers, which in turns affect the values of the fitted parameters if the
correct weights are not considered.

5 Asymptotic results for large sample size

As seen previously, the error structure changes and several assumptions of the regression model are violated if one uses the unweighted least squares method in transformed spaces. This demonstration is enough of a proof to prefer a nonlinear regression on size-effect data in the space (P vs. D). From a practical perspective, it is however interesting to understand how the least squares point estimates of the parameters are changed. The least squares estimators are different if the data is fitted in the untransformed or the transformed spaces without using the right expression of least squares function S(θ). The objective of this section is to quantify this difference for large sample size and assess if the difference is statistically significant.

5.1 Law of large numbers and point estimates

Using the law of large numbers, the approximation error in the fitted parame-
ters can be made arbitrarily small when the number of tested specimens is large
enough. A new set of simulated data was generated using this time 250,000 repli-
cates for each of the four specimen sizes x₁, ..., x₄, for a total of N = 1,000,000 data
points. To simplify the notation, let us drop the double summation accounting for
replicates and use the subscript t to designate a data point. In the untransformed
space, the nonlinear least squares estimator is obtained by minimizing Equation 11.
$\hat{\theta}_{\infty} = [\hat{\theta}_{1,\infty} \hat{\theta}_{2,\infty}]^t = [40.0006 \ 14.9977]^t$ was found. One can define a different nonlinear least squares estimator by minimizing

$$S^{(1)}(\theta) = \sum_{t=1}^{N} \left( y_t^{(1)} - \frac{c_3 c_1 \sqrt{\theta_1}}{\sqrt{c_2 \theta_2 + x_t}} \right)^2 \quad (25)$$

which corresponds to fitting the size-effect law directly in the space $(\sigma_N \text{ vs. } D)$. $\hat{\theta}_{\infty}^{(1)} = [\hat{\theta}_{1,\infty}^{(1)} \hat{\theta}_{2,\infty}^{(1)}]^t = [40.003 \ 15.0283]^t$ was obtained. A third nonlinear least squares estimator can be defined by minimizing

$$S^{(2)}(\theta) = \sum_{t=1}^{N} \left[ \ln(y_t^{(1)}) - \ln\left(\frac{c_3 c_1 \sqrt{\theta_1}}{\sqrt{c_2 \theta_2 + x_t}}\right) \right]^2 \quad (26)$$

which is equivalent to finding the fracture parameters directly in the space $(\ln(\sigma_N) \text{ vs. } \ln(D))$. $\hat{\theta}_{\infty}^{(2)} = [\hat{\theta}_{1,\infty}^{(2)} \hat{\theta}_{2,\infty}^{(2)}]^t = [40.8701 \ 16.4084]^t$ was obtained from the simulated data. Finally, if one were to follow the RILEM recommendation, a fourth nonlinear least squares estimator could be defined by minimizing

$$S^{(3)}(\theta) = \sum_{t=1}^{N} \left[ \frac{1}{(y_t^{(1)})^2} - \left( \frac{c_2 \theta_2}{(c_3 c_1)^2 \theta_1} + \frac{1}{(c_3 c_1)^2 \theta_1} x_t \right) \right]^2 \quad (27)$$

in the space $(1/\sigma_N^2 \text{ vs. } D)$. $\hat{\theta}_{\infty}^{(3)} = [\hat{\theta}_{1,\infty}^{(3)} \hat{\theta}_{2,\infty}^{(3)}]^t = [42.0466 \ 19.0027]^t$ was obtained. Aside from the asymptotic estimates in spaces $(P \text{ vs. } D)$ and $(\sigma_N \text{ vs. } D)$ which converged to the true parameters $[40 \ 15]^t$, the estimators in spaces $(\ln(\sigma_N) \text{ vs. } \ln(D))$ and $(1/\sigma_N^2 \text{ vs. } D)$ are asymptotically biased, i.e., the estimates do not tend to the true parameters as the sample size tends to infinity. Let us determine whether this asymptotic bias is statistically significant. In other words, let us assess whether $\hat{\theta}_{\infty}^{(2)}$ and $\hat{\theta}_{\infty}^{(3)}$ are within the marginal confidence intervals and confidence contours of $\hat{\theta}$ through classical asymptotic theory and Monte Carlo simulations. The properties of each nonlinear least squares estimators will be discussed later in Section 6.

## 5.2 Asymptotic theory applied to the size-effect model

Asymptotic results on nonlinear regression models with unweighted least squares functions are well-established. A full derivation of the asymptotic analysis can be found for instance in Ivanov, 1997 [25] or to some extent in Seber and Wild, 1989 [41]. Here, a brief and simplified summary of it is described to explain which equations were implemented numerically for Monte Carlo simulations shown next.

The regression model writes as

$$y_t = f(x_t, \theta^0) + \varepsilon_t \quad (28)$$
where $\theta^0 = [\theta_1^0 \theta_2^0]^t$ are the true values of $\theta$ and $f$ is defined in Equation 10. The errors $\epsilon_t$ are i.i.d. with $E[\epsilon_t] = 0$ and $\text{Var}[\epsilon_t] = \sigma^2$. The assumption of normality of errors is not necessary and will not be used here. Let us define $\mathbf{y} = [y_1 \dots y_N]^t$ and $\mathbf{f}(\theta) = [f_1(\theta) \dots f_N(\theta)]^t$ where $f_t(\theta) = f(x_t, \theta)$ and the Jacobian of $\mathbf{f}$ as

$$ \mathbf{F}(\theta) = \frac{\partial \mathbf{f}(\theta)}{\partial \theta^t} \qquad (29) $$

where for $1 \le t \le N$

$$ \frac{\partial f_t(\theta)}{\partial \theta_1} = \frac{c_1 x_t}{2\sqrt{\theta_1} \sqrt{c_2 \theta_2 + x_t}} \qquad (30) $$

$$ \frac{\partial f_t(\theta)}{\partial \theta_2} = - \frac{c_1 c_2 \sqrt{\theta_1 x_t}}{2(c_2 \theta_2 + x_t)^{3/2}} \qquad (31) $$

Let us denote $\mathbf{F} = \mathbf{F}(\theta)$ and $\mathbf{F}^0 = \mathbf{F}(\theta^0)$. $\hat{\theta}$ is the least squares estimate of $\theta^0$ which can be obtained by minimizing the error sum of squares

$$ \begin{aligned} S(\theta) &= [\mathbf{y} - \mathbf{f}(\theta)]^t [\mathbf{y} - \mathbf{f}(\theta)] \\ &= \|\mathbf{y} - \mathbf{f}(\theta)\|^2 \end{aligned} \qquad (32) $$

In a small neighborhood of $\theta^0$, one can perform a Taylor expansion of $f$ to the first order and write the least squares function as

$$ \begin{aligned} S(\theta) &\approx \|\mathbf{y} - \mathbf{f}(\theta^0) - \mathbf{F}^0(\theta - \theta^0)\|^2 \\ &= \|\varepsilon - \mathbf{F}^0\beta\|^2 \end{aligned} \qquad (33) $$

where $\mathbf{y} - \mathbf{f}(\theta^0) = \varepsilon$ and $\beta = \theta - \theta^0$. Using the results for linear model regression, the least squares estimate writes as

$$ \hat{\beta} = [(\mathbf{F}^0)^t \mathbf{F}^0]^{-1} (\mathbf{F}^0)^t \varepsilon \qquad (34) $$

When $N$ becomes large, $\hat{\theta}$ becomes almost certain to lie within a small neighborhood of $\theta^0$. It can be shown [25] that

$$ \begin{aligned} \hat{\theta} - \theta^0 &= \hat{\beta} + o_p \left( \frac{1}{\sqrt{N}} \right) \\ &= [(\mathbf{F}^0)^t \mathbf{F}^0]^{-1} (\mathbf{F}^0)^t \varepsilon + o_p \left( \frac{1}{\sqrt{N}} \right) \end{aligned} \qquad (35) $$

The notation $o_p(1/\sqrt{N})$ denotes a matrix valued random variable $\mathbf{R}_N$ with the property that each element $R_{ijN}$ satisfies $\lim_{N \to \infty} P(|\sqrt{N}R_{ijN}| > \epsilon) = 0$ for any $\epsilon > 0$. where the operator $P()$ should be read as “probability of”. Equation 35 implies that $\hat{\theta} - \theta^0$ is asymptotically a linear combination of the i.i.d. $\varepsilon$. An appropriate version of the central limit theorem can
be used [25] to show that $\hat{\theta} - \theta^0$ is normally distributed. For large $N$, $\hat{\theta}$ converges almost surely to $\theta^0$, i.e., $\hat{\theta}$ is asymptotically unbiased. One can write the following approximate expression

$$ \hat{\theta} - \theta^0 \sim N_2(0, V) \qquad (36) $$

where $N_2$ denotes a double-variate normal distribution and

$$ V = \sigma^2 \left[ (F^0)^t F^0 \right]^{-1} \qquad (37) $$

is the asymptotic variance-covariance matrix.

<figure><img src="image_7.png" /><figcaption>Fig. 4: Distribution of $\hat{\theta}$ for (a) and (d) $N = 6 \times 4 = 24$, rejected $H_0$; (b) and (e) $N = 50 \times 4 = 200$, rejected $H_0$; (c) and (f) $N = 100 \times 4 = 400$. failed to reject $H_0$. Bin sizes: 0.5 N m⁻¹ for $\hat{\theta}_1$ and 0.5 mm for $\hat{\theta}_2$.</figcaption></figure>

## 5.3 Sample size for asymptotics

In general, the estimation of standard error metrics to obtain marginal confidence intervals on simulated data should not be done using a very large sample size $N$. The reason is that the standard error becomes arbitrarily small such that the smallest difference between the estimate, i.e., $\hat{\theta}^{(2)}$ and $\hat{\theta}^{(3)}$, and the null hypothesis, i.e., $\theta^0 \approx \hat{\theta} \approx \hat{\theta}^{(1)}$, becomes statistically significant. A remedy would be to decrease the alpha
level as the sample size increases, which then brings subjectivity to the analysis. A different strategy was chosen here. The sample size N was increased incrementally with $n_i = 6, 50, 100$ replicates for each of the four regressor values ($i = 1, ..., 4$). For each sample size, the least squares estimate $\hat{\theta}$ was obtained by minimizing $S(\theta)$ in the untransformed space. This step was repeated 5,000 times to obtain a distribution of $\hat{\theta}$, with $\varepsilon_t$ i.i.d. following a normal distribution with $E[\varepsilon_t] = 0$ and $Var[\varepsilon_t] = \sigma^2$ defined in Section 3.4. For large enough N, the distribution of $\hat{\theta}$ converges to a double-variate normal distribution based on the asymptotic theory. To quantify the threshold N at which the asymptotic results start to be valid, an Anderson-Darling test statistic was performed for each value of N with an alpha level of 0.05, with the following hypothesis

$H_0$: $\hat{\theta}$ follows a double-variate normal distribution.

$H_1$: $\hat{\theta}$ does not follow a double-variate normal distribution.

Figure 4 shows the distributions of $\hat{\theta}$ for $N = 24, 200, 400$. The threshold for total sample size was found to be $N = 4 \times 100 = 400$, for which the Anderson-Darling test statistic failed to reject the null hypothesis.

## 5.4 Marginal confidence intervals

Using this threshold value, N no longer tends to infinity, one can now estimate the corresponding asymptotic variance-covariance matrix as

$$ 
\begin{aligned}
\hat{V} &= \mathrm{Var}[\hat{\theta}] \\
&= s^2 (\hat{F}^t \hat{F})^{-1}
\end{aligned}
\qquad (38) $$

where $\hat{F} = F(\hat{\theta})$, and $s^2$ is the estimate of the variance $\sigma^2$ calculated as

$$ s^2 = \frac{S(\hat{\theta})}{N} \qquad (39) $$

Note that in nonlinear regression statistics, using $N-p$ degrees of freedom with $p=2$ being the number of parameters as opposed to N does not lead to an unbiased estimate of $\sigma^2$. Since the estimates of the parameters follow a double-variate normal distribution, the estimate of $\sigma^2$ becomes asymptotically unbiased [24]. The diagonal elements of the asymptotic variance-covariance matrix are estimated by $\hat{V}_{11}$ and $\hat{V}_{22}$. They correspond to the estimates of the asymptotic variances of $\hat{\theta}_1$ and $\hat{\theta}_2$, respectively. Let us denote by $\mathrm{se}(\hat{\theta}_1) = \sqrt{\hat{V}_{11}}$ and $\mathrm{se}(\hat{\theta}_2) = \sqrt{\hat{V}_{22}}$ the estimates of the asymptotic standard errors on the model parameters, respectively. By analogy to the normal linear regression case, the marginal confidence intervals for $\hat{\theta}_p$ with an alpha level of 0.05 can
be computed [24] and write as

$$ \hat{\theta}_p \pm \text{se} \left( \hat{\theta}_p \right) \sqrt{\frac{N}{N-2}} t_{0.975} \quad (40) $$

where $p = 1, 2$ and $t_u$ is the $u$ percentile of a Student's t-distribution with $N - 2$ degrees of freedom.

Monte Carlo simulations were performed to generate 5,000 distributions of $\hat{\theta}$. For each distribution, the marginal confidence intervals were calculated and the presence of $\hat{\theta}_{\infty}^{(2)}$ and $\hat{\theta}_{\infty}^{(3)}$ within the interval limits was checked.

The probabilities for $\hat{\theta}_1^{(2)}$ and $\hat{\theta}_2^{(2)}$ to be within the confidence limits were 77.16% and 66.02%, respectively. The probabilities for $\hat{\theta}_1^{(3)}$ and $\hat{\theta}_2^{(3)}$ to be within the marginal confidence intervals were low or very low, 23% and 2.56%, respectively. Figure 5(a) shows the marginal confidence limits represented by dashed lines for one $\hat{\theta}$ distribution.

<figure><img src="image_8.png" /><figcaption>Fig. 5: Asymptotic inference with an alpha level of 0.05: (a) a simulation showing the marginal confidence intervals and (b) another simulation showing the ellipsoidal confidence contour.</figcaption></figure>

## 5.5 Confidence contours

The two dimensional rectangular region defined by the marginal confidence limits does not however characterize a joint confidence. Indeed, the individual parameters might be within the limits but the parameter vector might not be contained collectively in the joint region. One solution is to build a confidence region based on the least squares function $S(\theta)$, which measures the distance between the data to the model response
for any $\theta$. By setting $S(\theta) = \text{constant}$, one can construct an exact confidence contour. Although the contour is correct, an exact alpha level cannot be specified because of the non-linearity of the regression model. For large sample size $N$, however, a very close approximation of the alpha level can be obtained as $\hat{\theta}$ becomes sufficiently close to $\theta^0$, because the Taylor expansion in Equation 33 becomes more accurate. In that case, the confidence contour with an alpha level of 0.05 can be derived [14] and is defined by

$$S(\theta) \le S(\hat{\theta}) \left(1 + \frac{2}{N-2} F_{2,N-2}^{0.05}\right) \quad (41)$$

where $F_{2,N-2}^\alpha$ is the upper $\alpha = 0.05$ quantile for Fisher's $F$ distribution with 2 and $N-2$ degrees of freedom. For $N=400$, the right hand side of Equation 41 simplifies as $\hat{S}_{0.05} = 1.0152S(\hat{\theta})$. Finally, the confidence region can be calculated by finding $\theta$ such that

$$\sum_{t=1}^{N=400} \left( y_t - \frac{c_1 \sqrt{\theta_1} x_t}{\sqrt{c_2 \theta_2 + x_t}} \right)^2 \le \hat{S}_{0.05} \quad (42)$$

Once again, Monte Carlo simulations were performed to generate 5,000 distributions of $\hat{\theta}$. For each distribution, the confidence region was this time calculated and the presence of $\hat{\theta}_\infty^{(2)}$ and $\hat{\theta}_\infty^{(3)}$ within the region was checked.

The probability for $\hat{\theta}_\infty^{(2)}$ to be within the confidence contour was 75.18% whereas the probability for $\hat{\theta}_\infty^{(3)}$ to be within the confidence region was almost null (1.08%). Figure 5(b) shows the ellipsoidal confidence contour observed for one of the $\hat{\theta}$ distributions of the Monte Carlo simulations, where the RILEM estimator is clearly out of the 95% confidence region. On a side note, one can notice the orientation and symmetrical shape of the contour which indicate that $\hat{\theta}_1$ is as well determined as $\hat{\theta}_2$.

# 6 Properties of the nonlinear least squares estimators

## 6.1 Asymptotic bias

Four different nonlinear least squares estimators were defined in Section 5.1 corresponding to four different least squares functions $S(\theta)$, $S^{(1)}(\theta)$, $S^{(2)}(\theta)$, and $S^{(3)}(\theta)$. The estimators defined in spaces $(P \text{ vs. } D)$ and $(\sigma_N \text{ vs. } D)$ are asymptotically unbiased, i.e., the estimates tend to the true values when the sample size tends to infinity. On the other hand, the estimators in spaces $(\ln(\sigma_N) \text{ vs. } \ln(D))$ and $(1/\sigma_N^2 \text{ vs. } D)$ are asymptotically biased. Biased estimators should not be employed in regression statistics. In addition, the bias of these two estimators is statistically significant as the probabilities to obtain the asymptotic estimates within the confidence region with an alpha level of 0.05 are only 75% and 1% for the estimators in $(\ln(\sigma_N) \text{ vs. } \ln(D))$ and $(1/\sigma_N^2 \text{ vs. } D)$, respectively. In other words, even from a practical point of view, the size-effect law should not be fitted in spaces $(\ln(\sigma_N) \text{ vs. } \ln(D))$ or $(1/\sigma_N^2 \text{ vs. } D)$, the latter being suggested by RILEM.
## 6.2 Minimum variance

The asymptotic theory ensures that the nonlinear least squares estimator corresponding to the space (P vs. D) is mean square consistent, i.e., its variance converges to zero when the sample size tends to infinity. This can for instance be seen visually in Figure 4.

What about the estimator corresponding to space ($\sigma_N$ vs. D)? It was shown on one example that the estimator is asymptotically unbiased (see Section 5.1). However the error structure is changed (see Equation 14), which means that the two estimators cannot be identical. To understand the difference between these two estimators, let us reconsider the asymptotic theory described earlier. Let us add a few simple modifications to the theory to account for regressor-dependent weights.

The least squares function $S^{(1)}(\theta)$ defined in Equation 25 can be rearranged as

$$ 
\begin{aligned}
S^{(1)}(\theta) &= \sum_{t=1}^{N} \left(\frac{c_3}{x_t}\right)^2 \left(y_t - \frac{c_1 \sqrt{\theta_1} x_t}{\sqrt{c_2 \theta_2 + x_t}}\right)^2 \\
&= [\mathbf{y} - \mathbf{f}(\theta)]^t \mathbf{A}^2 [\mathbf{y} - \mathbf{f}(\theta)]
\end{aligned} 
\quad (43) $$

where $\mathbf{A}$ is a square matrix with general term $A_{ii} = c_3/x_i$ and $A_{ij} = 0$ for $i \neq j$, $1 \le i, j \le N$. One can then write

$$ 
\begin{aligned}
S^{(1)}(\theta) &= [\mathbf{y} - \mathbf{f}(\theta)]^t \mathbf{A}^t \mathbf{A} [\mathbf{y} - \mathbf{f}(\theta)] \\
&= [\mathbf{A}\mathbf{y} - \mathbf{A}\mathbf{f}(\theta)]^t [\mathbf{A}\mathbf{y} - \mathbf{A}\mathbf{f}(\theta)]
\end{aligned} 
\quad (44) $$

By defining $\mathbf{z} = \mathbf{A}\mathbf{y}$ and $\mathbf{g}(\theta) = \mathbf{A}\mathbf{f}(\theta)$, one can further write

$$ 
\begin{aligned}
S^{(1)}(\theta) &= [\mathbf{z} - \mathbf{g}(\theta)]^t [\mathbf{z} - \mathbf{g}(\theta)] \\
&= \|\mathbf{z} - \mathbf{g}(\theta)\|^2
\end{aligned} 
\quad (45) $$

which corresponds to the least squares function of the following nonlinear regression model

$$ \mathbf{z} = \mathbf{g}(\theta^0) + \eta \quad (46) $$

where $\eta = \mathbf{A}\mathbf{\epsilon}$ is the error term. Let us recall that $\mathbf{y} - \mathbf{f}(\theta^0) = \mathbf{\epsilon}$ where $\mathbf{E}[\mathbf{\epsilon}] = 0$ and $\mathrm{Var}[\mathbf{\epsilon}] = \sigma^2$. One can therefore deduce that

$$ \mathbf{E}[\eta] = 0 \quad (47) $$

whereas the variance reads as

$$ 
\begin{aligned}
\mathrm{Var}[\eta] &= \mathrm{Var}[\mathbf{A}\mathbf{\epsilon}] \\
&= \mathbf{A}\mathrm{Var}[\mathbf{\epsilon}]\mathbf{A} \\
&= \mathrm{Var}[\mathbf{\epsilon}]\mathbf{A}^2 \\
&= \sigma^2\mathbf{A}^2
\end{aligned} 
\quad (48) $$
Applying the asymptotic theory to the nonlinear regression model (Equation 46), one can deduce that

$$ \hat{\theta}^{(1)} - \theta^{0} \sim N_{2} \left( 0, \sigma^{2} \left[ (G^{0})^{t} G^{0} \right]^{-1} \right) \qquad (49) $$

where $G^0 = G(\theta^0)$ and

$$ \begin{aligned} G(\theta) &= \frac{\partial g(\theta)}{\partial \theta^t} \\ &= A \frac{\partial f(\theta)}{\partial \theta^t} \\ &= AF \end{aligned} \qquad (50) $$

The nonlinear least squares estimator in space $(\sigma_N \text{ vs. } D)$ is indeed asymptotically unbiased. The asymptotic variance-covariance matrix simplifies as

$$ \begin{aligned} V^{(1)} &= \sigma^2 \left[ (AF^0)^t AF^0 \right]^{-1} \\ &= \sigma^2 \left[ (F^0)^t A^2 F^0 \right]^{-1} \end{aligned} \qquad (51) $$

The estimate of this asymptotic variance-covariance matrix is simply defined as

$$ \begin{aligned} \hat{V}^{(1)} &= \mathrm{Var} \left[ \hat{\theta}^{(1)} \right] \\ &= \left( s^{(1)} \right)^2 \left( \hat{F}^t A^2 \hat{F} \right)^{-1} \end{aligned} \qquad (52) $$

where

$$ \left( s^{(1)} \right)^2 = \frac{S^{(1)} \left( \hat{\theta} \right)}{N} \qquad (53) $$

Therefore, the estimator is also mean square consistent.

There is however a difference between the two least squares estimators: the expressions of the asymptotic variance-covariance matrices in Equation 37 and in Equation 51 differ from a factor $A^2$. Since the spectral norm of this matrix is always smaller than 1, the variance of the estimator in $(\sigma_N \text{ vs. } D)$ is always larger than the variance of the estimator in $(P \text{ vs. } D)$. In other words, the nonlinear least squares estimator in $(\sigma_N \text{ vs. } D)$ does not possess the minimum variance property and is thus less efficient than the estimator in $(P \text{ vs. } D)$.

## 6.3 Small sample size

In general, most of the size-effect tests performed in laboratory involve a small or a very small number of fracture tests including replicates. As a matter of fact, 53 out of the 71 datasets considered in this study only involved 12 specimens or less. For such small sample sizes, there is no theoretical result known about the superiority of one least
squares estimator over another. However there is no motivation in using a less efficient estimator or an asymptotically biased estimator like the one suggested by RILEM when the current computational capabilities to perform nonlinear optimization of a specified objective function are nowadays widespread.

As an illustration of the possible differences in fracture parameters one might obtain for small sample sizes, a set of simulated data was generated using this time 4 replicates for each of the four specimen sizes $x_1, ..., x_4$, for a total of $N = 16$ data points. The errors were assumed to be i.i.d. and drawn from a normal distribution with mean 0 and variance 0.25. The point estimate in space (P vs. D) was obtained by minimizing Equation 11, and $\hat{\theta} = [\hat{\theta}_1 \ \hat{\theta}_2]^t = [40.039 \ 16.1753]^t$ was found. The same data was fitted in the space $(1/\sigma_N^2 \text{ vs. } D)$ following the RILEM recommendation by minimizing Equation 27. $\hat{\theta}^{(3)} = [\hat{\theta}_1^{(3)} \ \hat{\theta}_2^{(3)}]^t = [51.8167 \ 37.3466]^t$ was obtained. The differences between the fitted values using the RILEM recommendation and the true values of $\theta$, i.e., $\theta^0 = [\theta_1^0 \ \theta_2^0]^t = [40 \ 15]^t$ are consequent. Indeed, the fracture energy is larger than its true value by about 30% while the effective length is larger by 150%. This can be easily explained by observing the fitted curves shown in Figure 6. The RILEM transformation generates artificial outliers in space $(1/\sigma_N^2 \text{ vs. } D)$ that were nonexistent in the untransformed space. The reader can refer to Section 4.5 for more details on the changes in error structure. As a consequence, the outliers make the unweighted fitting in the space $(1/\sigma_N^2 \text{ vs. } D)$ incorrect. The approximated expression of the weights that should have been considered for a correct fitting is given in Equation 24. This example serves only for illustration purposes and one might find

<figure><img src="image_9.png" /><figcaption>**Fig. 6:** Generated data fitted in spaces (P vs. D) and $(1/\sigma_N^2 \text{ vs. } D)$: (a) representation in space (P vs. D) and (b) representation in space $(1/\sigma_N^2 \text{ vs. } D)$.</figcaption></figure>
other examples in which the fitting in the untransformed or the transformed spaces provide similar fracture parameters for very small size samples. However, in a case like the one described above, the simple presence of artificially generated outliers and the spurious effects they produce should convince the reader on the importance of using the only theoretically well-justified regression framework developed in this paper.

Finally, let us mention that the asymptotic theory can still be applied as a first approximation to real small sample size regression analysis, for which the practical significance is preferred to the statistical significance. One could lift the rigorous Anderson-Darling test statistic and determine whether the departure of the distributions of the parameters from normality is serious or not for practical purposes. This could be done for instance by simply calculating the skewness and excess kurtosis of the distributions. An alternative to the asymptotic theory applied to small sample sizes is to use the bootstrap distributions of the parameters [15] and define the confidence region from these distributions. There is however no theoretical result about the superiority of the bootstrapping method or the asymptotic theory for small size samples [24].

# 7 Conclusion

In the vast majority of contributions using the size-effect law to determine fracture parameters, the linear regression plot $(1/\sigma_N^2$ vs. $D)$ recommended by RILEM [38] is employed. Other transformations, i.e., $(\sigma_N$ vs. $D)$ and $(\ln(\sigma_N)$ vs. $\ln(D))$, are sometimes also used [7]. In this paper, the effect of such model transformations was investigated. The following conclusions can be drawn.

* Each of the transformations violates the assumptions that the error term should be independent of the expectation function, and that the errors should be independently distributed.

* All the three transformations lead to highly skewed error distributions with a non-negligible probability to obtain a significant number of outliers.

* The nonlinear least squares estimators in spaces $(\ln(\sigma_N)$ vs. $\ln(D))$ and $(1/\sigma_N^2$ vs. $D)$ are asymptotically biased. The difference with the correct point estimates is statistically significant with an alpha level of 0.05.

* The nonlinear least squares estimator in space $(\sigma_N$ vs. $D)$ is asymptotically unbiased and mean square consistent, but is less efficient than the estimator in the untransformed space ($P$ vs. $D$) because it does not possess the minimum variance property.

In summary, the fracture parameters should only be obtained in the space ($P$ vs. $D$). The linear regression plot suggested by RILEM [38] should be avoided for the statistical inverse problem of the size-effect method.

# Declarations

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
# Acknowledgments

This research was performed within the University of New Mexico supported jointly by Université de Pau et des Pays de l'Adour and Northwestern University. Financial support from the Investissement d'Avenir French program (ANR-16-IDEX-0002) is gratefully acknowledged.

# Appendix

## 7.1 Brown-Forsythe test statistic

The Brown-Forsythe test is a nonparametric test for nonconstant variance. It consists in grouping the residuals into two groups. The first group contains the $N_1$ lowest values of the regressor, while the second group is composed of the $N_2$ highest values, such that $N_1 + N_2 = N$. The null hypothesis is $H_0$: the variance is constant, whereas the alternative hypothesis is $H_1$: the variance is not constant. The test statistic is calculated as

$$BF = \frac{\bar{d}_1 - \bar{d}_2}{s_{BF} \sqrt{\frac{1}{N_1} + \frac{1}{N_2}}} \qquad (54)$$

where $\bar{d}_1$

$$s_{BF} = \sqrt{\frac{\sum_{t=1}^{N_1} (d_{1,t} - \bar{d}_1)^2 + \sum_{t=1}^{N_2} (d_{2,t} - \bar{d}_2)^2}{N_1 + N_2 - 2}} \qquad (55)$$

and

$$d_{1,t} = |\varepsilon_{1,t} - \tilde{\varepsilon}_1| \text{ and } d_{2,t} = |\varepsilon_{2,t} - \tilde{\varepsilon}_2| \qquad (56)$$

$\varepsilon_{1,t}$ and $\varepsilon_{2,t}$ are the residuals in groups 1 and 2, respectively. $\tilde{\varepsilon}_1$ and $\tilde{\varepsilon}_2$ are the medians of the first and second groups, respectively. Finally, $\bar{d}_1$ and $\bar{d}_2$ are the mean values of the quantities $d_{1,t}$ and $d_{2,t}$, respectively, for $t = 1, ..., N_1$ and $t = 1, ..., N_2$. The quantity $BF^2$ is approximately distributed as an $F_{1,N_1+N_2-2}$ distribution.

## 7.2 Anderson-Darling test statistic

The Anderson-Darling test calculates the difference in area between the residuals and a fitted line of the normal distribution. The null hypothesis is $H_0$: the residuals follow a normal distribution, whereas the alternative hypothesis is $H_1$: the residuals do not follow a normal distribution. The statistic can be computed as

$$A^2 = -N - \sum_{t=1}^{N} \frac{2t-1}{N} [\ln(F_Z(\varepsilon_t)) + \ln(1 - F_Z(\varepsilon_{N+1-t}))] \qquad (57)$$

where $F_Z$ is the cumulative distribution function of the normal distribution. One can note that the statistic $A^2$ can be seen as a squared distance with larger weights in the tails of the distribution, and therefore smaller $A^2$ values indicate a better fitting of the distribution to the residuals.
References

[1] ACI-Committee (2019). Building code requirements for structural concrete (aci 318-19) and commentary. American Concrete Institute.

[2] Bates, D. M. and Watts, D. G. (1988). *Nonlinear regression analysis and its applications*. John Wiley & Sons.

[3] Bažant, Z., Gettu, R., and Kazemi, M. (1991). Identification of nonlinear fracture properties from size effect tests and structural analysis based on geometry-dependent r-curves. In *International journal of rock mechanics and mining sciences & geomechanics abstracts*, volume 28, pages 43–51. Elsevier.

[4] Bažant, Z. P. (1984). Size effect in blunt fracture: concrete, rock, metal. *Journal of engineering mechanics*, 110(4):518–535.

[5] Bažant, Z. P., Lee, S.-G., and Pfeiffer, P. A. (1987). Size effect tests and fracture characteristics of aluminum. *Engineering Fracture Mechanics*, 26(1):45–57.

[6] Bažant, Z. P. and Pfeiffer, P. A. (1987). Determination of fracture energy from size effect and brittleness number. *ACI Materials Journal*, 84(6):463–480.

[7] Bažant, Z. P. and Planas, J. (1997). *Fracture and Size Effect in Concrete and Other Quasibrittle Materials*. CRC press.

[8] Bazant, Z. P. and Xu, K. (1991). Size effect in fatigue fracture of concrete. *ACI Materials Journal*, 88(4):390–399.

[9] Bažant, Z. P., Zhou, Y., Zi, G., and Daniel, I. M. (2003). Size effect and asymptotic matching analysis of fracture of closed-cell polymeric foam. *International journal of solids and structures*, 40(25):7197–7217.

[10] Bulmer, M. G. (1979). *Principles of statistics*. Courier Corporation.

[11] Çağlar, Y. and Şener, S. (2016). Size effect tests of different notch depth specimens with support rotation measurements. *Engineering Fracture Mechanics*, 157:43–55.

[12] Carroll, R. J. and Ruppert, D. (2017). *Transformation and weighting in regression*. Chapman and Hall/CRC.

[13] Chen, Y., Han, X., Hu, X., and Zhu, W. (2020). Statistics-assisted fracture modelling of small un-notched and large notched sandstone specimens with specimen-size/grain-size ratio from 30 to 900. *Engineering Fracture Mechanics*, 235:107134.

[14] Draper, N. R. and Smith, H. (1998). *Applied regression analysis*. John Wiley & Sons.
[15] Efron, B. and Tibshirani, R. J. (1994). *An introduction to the bootstrap*. Chapman and Hall/CRC.

[16] Eskandari, H., Muralidhara, S., Raghuprasad, B., and Venkatarama Reddy, B. (2010). Size effect in self consolidating concrete beams with and without notches. *Sadhana*, 35:303–317.

[17] Fathy, A. (1992). *Application of fracture mechanics to rocks and rocky materials*. PhD thesis, PhD Thesis, Dep. de Ciencia de Materiales, ETS de Ingenieros de Caminos ....

[18] Gettu, R., Bazant, Z. P., and Karr, M. E. (1990). Fracture properties and brittleness of high-strength concrete. *ACI Materials Journal*, 87(6):608–618.

[19] Gettu, R., Saldívar, H., and Kazemi, M. T. (1998). Implications of the size effect method for analyzing the fracture of concrete. *International journal of solids and structures*, 35(31-32):4121–4132.

[20] Ghasemi, M., Ghasemi, M. R., and Mousavi, S. R. (2019). Studying the fracture parameters and size effect of steel fiber-reinforced self-compacting concrete. *Construction and Building Materials*, 201:447–460.

[21] Grégoire, D., Rojas-Solano, L. B., and Pijaudier-Cabot, G. (2013). Failure and size effect for notched and unnotched concrete beams. *International Journal for Numerical and Analytical Methods in Geomechanics*, 37(10):1434–1452.

[22] He, S., Plesha, M., Rowlands, R., and Bažant, Z. (1992). Fracture energy tests of dam concrete with rate and size effects. *Dam Engineering*, 3(2):139–159.

[23] Hoover, C. G., Bažant, Z. P., Vorel, J., Wendner, R., and Hubler, M. H. (2013). Comprehensive concrete fracture tests: description and results. *Engineering Fracture Mechanics*, 114:92–103.

[24] Huet, S., Bouvier, A., Poursat, M.-A., Jolivet, E., and Bouvier, A. (2004). *Statistical tools for nonlinear regression: a practical guide with S-PLUS and R examples*. Springer.

[25] Ivanov, A. (2013). *Asymptotic theory of nonlinear regression*, volume 389. Springer Science & Business Media.

[26] Karamloo, M., Mazloom, M., and Payganeh, G. (2017). Effect of size on nominal strength of self-compacting lightweight concrete and self-compacting normal weight concrete: A stress-based approach. *Materials Today Communications*, 13:36–45.

[27] Karihaloo, B. L., Abdalla, H., and Xiao, Q. (2003). Size effect in concrete beams. *Engineering fracture mechanics*, 70(7-8):979–993.
[28] Kim, K.-T., Bažant, Z. P., and Yu, Q. (2013). Non-uniqueness of cohesive-crack stress-separation law of human and bovine bones and remedy by size effect tests. *International Journal of Fracture*, 181(1):67–81.

[29] Le Bellégo, C., Dubé, J. F., Pijaudier-Cabot, G., and Gérard, B. (2003). Calibration of nonlocal damage model from size effect tests. *European Journal of Mechanics-A/Solids*, 22(1):33–46.

[30] Li, W., Jin, Z., and Cusatis, G. (2019). Size effect analysis for the characterization of marcellus shale quasi-brittle fracture properties. *Rock Mechanics and Rock Engineering*, 52:1–18.

[31] Li, W., Qiao, Y., Fenner, J., Warren, K., Salviato, M., Bažant, Z. P., and Cusatis, G. (2021). Elastic and fracture behavior of three-dimensional ply-to-ply angle interlock woven composites: Through-thickness, size effect, and multiaxial tests. *Composites Part C: Open Access*, 4:100098.

[32] Linul, E., Marşavina, L., Sadowski, T., and Kneć, M. (2012). Size effect on fracture toughness of rigid polyurethane foams. *Solid State Phenomena*, 188:205–210.

[33] Mercuri, M., Pathirage, M., Gregori, A., and Cusatis, G. (2023). Influence of self-weight on size effect of quasi-brittle materials: generalized analytical formulation and application to the failure of irregular masonry arches. *International Journal of Fracture*, pages 1–28.

[34] Nurizada, A. and Kirane, K. (2020). Induced anisotropy in the fracturing behavior of 3d printed parts analyzed by the size effect method. *Engineering Fracture Mechanics*, 239:107304.

[35] Pijaudier-Cabot, G., Hajimohammadi, A., Nouailletas, O., La Borderie, C., Padin, A., and Mathieu, J.-P. (2022). Determination of the fracture energy of rocks from size effect tests: Application to shales and carbonate rocks. *Engineering Fracture Mechanics*, 271:108630.

[36] Planas, J. and Elices, M. (1990). Fracture criteria for concrete: mathematical approximations and experimental validation. *Engineering fracture mechanics*, 35(1-3):87–94.

[37] Rao, T. M. and Rao, T. G. (2013). Size effect of plain concrete beams–an experimental study. *International Journal of Research in Engineering and Technology*, 2(06):1047–1051.

[38] RILEM, D. R. (1990). Size-effect method for determining fracture energy and process zone size of concrete. *Materials and Structures*, 23:461–465.

[39] Salviato, M., Kirane, K., Ashari, S. E., Bažant, Z. P., and Cusatis, G. (2016). Experimental and numerical investigation of intra-laminar energy dissipation and
size effect in two-dimensional textile composites. *Composites Science and Technology*, 135:67–75.

[40] Salviato, M., Kirane, K., Bažant, Z. P., and Cusatis, G. (2019). Mode i and ii interlaminar fracture in laminated composites: a size effect study. *Journal of Applied Mechanics*, 86(9):091008.

[41] Seber, G. and Wild, C. (1989). *Nonlinear regression*. John Wiley & Sons.

[42] Sheskin, D. J. (2003). *Handbook of parametric and nonparametric statistical procedures*. Chapman and hall/CRC.

[43] Tang, T., Bažant, Z. P., Yang, S., and Zollinger, D. (1996). Variable-notch one-size test method for fracture energy and process zone length. *Engineering Fracture Mechanics*, 55(3):383–404.

[44] Wickens, T. D. and Keppel, G. (2004). *Design and analysis: A researcher's handbook*. Pearson Prentice-Hall Upper Saddle River, NJ.

[45] Zhang, S., An, D., Zhang, X., Yu, B., and Wang, H. (2021). Research on size effect of fracture toughness of sandstone using the center-cracked circular disc samples. *Engineering Fracture Mechanics*, 251:107777.