Hormesis

Hormesis is a process in which exposure to a low dose of a chemical or environmental factor that is damaging at higher doses induces an adaptive beneficial effect on the cell or organism (Mattson Reference Mattson2008). The extensive literature documenting herbicide-induced hormesis in plants provides convincing evidence for the reproducibility of hormesis as a phenomenon (see, for example, Belz and Duke Reference Belz and Duke2014). Despite indications that hormesis could be commercially leveraged to improve crop productivity, it has received limited attention in weed science literature. Nevertheless, faced with the looming gap between the demand for food and limited production availability, there is a growing consensus that the future focus in crop production is shifting from crop protection to crop enhancement (Belz et al. Reference Belz, Cedergreen and Duke2011), which may draw more attention toward the concept of hormesis.

The most commonly reported variable with herbicide-induced hormetic stimulation is an increase in plant dry weight across multiple species, including oat (Avena sativa L.) (Wiedman and Appleby Reference Wiedman and Appleby1972), soybean [Glycine max (L.) Merr.], corn (Zea mays L.), and eucalyptus (Eucalyptus grandis Hill ex Maiden) (Velini et al. Reference Velini, Alves, Godoy, Meschede, Souza and Duke2008). Another common observation includes an increase in plant height or root length in various species such as rough dog‘s-tail (Cynosurus echinatus L.), prickly lettuce (Lactuca seriola L.), rough mannagrass (Glyceria maxima (Hartm.) Holmb.), cotton (Gossypium hirsutum L.), and barley (Hordeum vulgare L.), among others (Belz and Duke Reference Belz and Duke2014). Only a few studies have reported the hormetic effects of herbicides on harvestable yield under field conditions, including increased grain yield from barley and sugar yield from sugarcane (Saccharum officinarum L.) (Pincelli-Souza et al. Reference Pincelli-Souza, Bortolheiro, Carbonari, Velini and Silva2020).

Hormesis is usually not recorded in studies with herbicides, primarily due to the infrequent use of doses that induce hormesis, because those doses typically fall significantly below the recommended levels for effective weed management (Belz and Duke Reference Belz and Duke2014). Yet, it is notable that dose-response studies with herbicides are increasingly gaining popularity, especially after the rise of dicamba drift in 2017. A substantial proportion of these studies use the classical ANOVA approach in data analysis. However, this is not the optimal method for analyzing such data because ANOVA pertains to categorical factors (Tabachnick and Fidell, Reference Tabachnik and Fidel2007), whereas regression focuses on quantitative explanatory variables. Therefore, given the quantitative nature of the dose as an explanatory variable, regression is the most appropriate choice for this type of study, which is further detailed in this manuscript.

Concept of NOAEL and LOAEL Toxicological (Safety) Thresholds

NOAEL represents the dose of a herbicide that does not result in any observable adverse effects on the plant. Conversely, LOAEL is the dose at which detrimental effects become measurable (Greim and Snyder Reference Greim, Snyder, Greim and Snyder2018), serving as a potential “safety” threshold for assessing the herbicide’s impact on the plant. From a practical standpoint, a threshold is a dose of a specific herbicide that triggers a measurable response in the plant species of interest. Considering the inherent variation with field experimentation, weather, and field conditions, a certain level of flexibility must be exercised, such as 1%, 2.5%, and 5% thresholds. For example, a 1% threshold represents a herbicide dose that causes a 1% reduction in response (e.g., dry matter, yield, or other observed variable), compared with a nontreated check (Milosevic et al. Reference Milosevic, Osipitan, Scott and Knezevic2023). The same analogy applies to 2.5% and 5% thresholds.

Furthermore, considering natural variability in agriculture research and experimental error, we suggest that the 1% or 2.5% threshold levels should be denoted as the range for reporting NOAEL. Meanwhile, the 5% threshold should be denoted as LOAEL, which should be in line with traditional statistical significance at 95% (α = 0.05), robust detection of treatment effects despite field experimentation variability, and practical acceptance by crop producers and practitioners. The establishment of NOAEL as a 1% to 2.5% threshold range attributes the observed reduction to random error, rather than treatment effect. This aligns with the NOAEL definition outlined by Alexeeff at al. (2002) as the dose at which no statistically or biologically significant adverse effects occur. Conversely, setting LOAEL at 5% accounts for a significant increase in adversities, which would be indicative of treatment-induced effects amid experimental variabilities and errors, as proposed in our methodology. LOAEL values can be calculated for each plant response of interest (e.g., dry matter reduction, yield reduction, height reduction, or visual injury). To obtain results of NOAEL, if an experiment features multiple calculated values (for several responses), the lowest dose should be considered. Determining these thresholds can also be beneficial for safety assessments, risk management, legal and regulatory compliance, and the protection of both the environment and public health.

Dose-Response and Associated Curves

In the field of weed science, dose usually refers to the amount of herbicide required to achieve a desirable effect on a specific plant species (Knezevic et al. Reference Knezevic, Streibig and Ritz2007). Consequently, the connection between herbicide dose and plant response holds significant importance for comprehending herbicide efficacy. Furthermore, understanding this relationship is critically necessary for the proper design and interpretation of the dose-response studies.

Typically, the shape of a dose-response curve is sigmoidal (Figure 1A), with upper and lower limits, or asymptotes. The upper limit is established based on the response observed in nontreated plants (control) or those exposed to an extremely low herbicide dose. Conversely, the lower limit is determined by the response levels observed when plants are subjected to a high herbicide dosage (Knezevic et al. Reference Knezevic, Streibig and Ritz2007). It is important to note that limits can be altered, depending on the curve direction (ascending or descending). For a sigmoidal (symmetric) curve, the dose corresponding to the midpoint of plant response observed between the upper and lower limits is usually referred to as the effective dose 50, or ED₅₀; that is, the dose required to result in a 50% reduction in observed response (Knezevic et al. Reference Knezevic, Streibig and Ritz2007). An effective dose (ED) may have any value between 1% and 99%, depending on the desired effect (ED _x ). The ED₁, ED_2.5, and ED₅ represent the effective doses associated with achieving 1%, 2.5%, and 5% of change in the desired plant response, establishing a direct connection between the concepts of toxicological thresholds and the parameters of the effective dose.

Figure 1. Typical dose-response curves. A: Sigmoidal ascending (solid line) and sigmoidal descending (dashed line) curve, described by Equations 1 and 2 (LL model family). B: J-shaped (solid line) and inverted J-shaped (dashed line) curve, described by Equations 3 and 4 (BC and CRS model family). C: U-shaped (solid line) and inverted U-shaped (dashed line) curve, described by Equation 5 (UCRS model family). Both J- and U-shaped curves suggest hormesis response, while a sigmoidal curve does not.

Although multiple models exist for describing sigmoidal dose-response curves (Figure 1A), the log-logistic model with three (LL.3) or four (LL.4) parameters is the one most commonly used (Knezevic et al. Reference Knezevic, Streibig and Ritz2007; Van der Vilet and Ritz Reference Van der Vilet and Ritz2013). The LL.4 function is shown in Equation 1:

([1]) $$f\left[ {x\left( {b,c,d,e} \right)} \right] = {\rm{c}} + {{d - c} \over {1 + {\rm{exp}}\left\{ {\left[ {\log \left( {\rm{x}} \right) - \log \left( {\rm{e}} \right)} \right]} \right\}}}$$

where e represents the ED₅₀, d signifies the upper horizontal asymptote (limit), and the lower horizontal asymptote (limit) is denoted with c. Parameter b indicates the relative slope around e. When the lower limit is 0, the resulting model is LL.3, shown in Equation 2:

([2]) $$f\left[ {x\left( {b,d,e} \right)} \right] = {d \over {1 + {\rm{exp}}\left\{ {{\rm{b}}\left[ {{\rm{log}}\left( {\rm{x}} \right) - {\rm{log}}\left( {\rm{e}} \right)} \right]} \right\}}}$$

Both LL.3 and LL.4 are used to characterize the response in which hormesis does not occur. However, if an initial growth is observed, data can be described with several models, most commonly by the Brain-Cousens model with four (BC.4) and five (BC.5) parameters (Equation 3):

([3]) $$f\left[ {x\left( {b,c,d,e,f} \right)} \right] = c + {{d - c + fx} \over {1 + {\rm{exp}}\left\{ {{\rm{b}}\left[ {{\rm{log}}\left( {\rm{x}} \right) - {\rm{log}}\left( {\rm{e}} \right)} \right]} \right\}}}$$

where, similarly to LL models, parameters c and d represent the lower and upper horizontal asymptotes (limits), respectively. Parameters b and e lack specific interpretation. However, parameter f holds significance because it represents the magnitude or size of the hormetic effect. Larger values of the f parameter indicate a more pronounced or substantial hormesis effect. Fixing the lower horizontal asymptote (limit) at 0 yields a BC.4 equation. Another model commonly used to evaluate and quantify the hormesis is the Cedergreen-Ritz-Streibig (CRS) model with five or six parameters, a modification of the log-logistic curve to account for hormesis (CRS.5, CRS.6), given by Equation 4 (Cedergreen et al. Reference Cedergreen, Ritz and Streibig2005) (Figure 1B) for some fixed, positive value of $\alpha $ :

([4]) $$f\left[ {x\left( {b,c,d,e,f} \right)} \right] = c + {{d - c + f{\rm{exp}}\left( { - {1 \over {{x^\alpha }}}} \right)} \over {1 + {\rm{exp}}\left\{ {{\rm{b}}\left[ {{\rm{log}}\left( {\rm{x}} \right) - {\rm{log}}\left( {\rm{e}} \right)} \right]} \right\}}}$$

Parameters c and d are still interpreted as lower and upper limits, the steepness of the curve following the maximal hormetic effect is denoted by the size of b. Parameter e provides a lower bound of the ED₅₀ level. The upper limit of the hormesis effect is determined by d + f. Within drc package are three variations of this model: CRS.5a, CRS.5b, and CRS.5c, for $\alpha $ being equal to 1, 0.5, and 0.25, respectively. For CRS.6, the last parameter is a freely varying $\alpha $ . Note that not all features, such as EDx calculation, are available for the model with the freely varying $\alpha $ (Cedergreen et al. Reference Cedergreen, Ritz and Streibig2005).

These CRS models describe the J-shaped (Figure 1B) hormesis data, while the U-shaped (Figure 1C) model is given by Equation 5:

([5]) $$f\left[ {x\left( {b,c,d,e,f} \right)} \right] = cd - {{d - c + f{\rm{exp}}\left( { - {1 \over {{x^\alpha }}}} \right)} \over {1 + {\rm{exp}}\left\{ {{\rm{b}}\left[ {{\rm{log}}\left( {\rm{x}} \right) - {\rm{log}}\left( {\rm{e}} \right)} \right]} \right\}}}$$

The statistical test for hormesis in each of the models presented above is given by the test of f ≠ 0, with the additional requirement of f being positive to indicate the presence of hormesis. While each of the models has its own set of advantages and disadvantages, these specific pros and cons are beyond the scope of this manuscript.

Optimal Dose Selection

Determining an appropriate number of doses in dose-response studies is a critical part of the planning process. From a statistical standpoint, the common rule of thumb is that the number of doses must be at least one or two doses (e.g., data points) higher than the number of parameters of the intended equation (model) used to describe (graph) the relationships between the dose and plant’s response (Knezevic et al. Reference Knezevic, Streibig and Ritz2007). For example, at least seven doses should be used when fitting an equation with six parameters (e.g., Equation 5). If time, space, and funding allow, the ideal scenario should have 10 doses (including a nontreated check), which would be handy for robust analysis. While selecting doses, it is beneficial to have about three data points for the lower and upper ends of the sigmoidal curve and perhaps four data points for the slope region of the curve. However, this might be hard to achieve from a practical standpoint, especially when testing multiple application timings with various crop or weed growth stages, which can limit (or reduce) the overall number of treatments. Nevertheless, it is essential to maintain a balance, ensuring that the number of doses is not fewer than seven. This would allow a description of the response using the regression model with up to six parameters, without the risk of overfitting the curves and inflating the estimated parameters, thus maintaining the integrity of the analysis. It is important to note that a higher number of data points typically yields a more precise estimation of parameters, and reduces standard errors of model parameters and standard errors of estimated ED values.

Considering the relative nature of plant response to a specific dose, the set of tested doses should be carefully selected and perhaps targeted toward the specific zone of the curve. To achieve this, a thorough review of the existing literature must be conducted, or perhaps one could rely on the previous experience of the researcher. When no published information exists, nor prior knowledge is available, conducting one or more preliminary screening trials can be a valuable strategy to determine how the organism responds to the initial set of doses.

The distribution of doses should also align with the specific objectives of the study. For example, when investigating a potential hormetic effect, it is advisable to include several doses that are lower or around the NOAEL and coupled with a few doses on the higher end of the spectrum. This would likely cover a suitable range for detecting any hormetic responses that may occur at lower doses, while the few higher doses would induce the inhibitory phase of hormesis.

On the other hand, if the primary objective is to describe the NOAEL and LOAEL, at least two doses lower than the expected NOAEL should be coupled with a few mid-range doses and one dose at the higher end, for a total of seven doses. This distribution enables a more precise determination of the dose-response curve, slope (b), and the ED₅₀, but it also enhances the accuracy in predicting the NOAEL and LOAEL values.

Furthermore, the set of doses should be arranged in even step-wise increments (Knezevic et al. Reference Knezevic, Streibig and Ritz2007). For instance, if seven herbicide doses are used, they should be structured as 0, 50, 100, 150, 200, 250 and 300 g ae ha⁻¹, while in preliminary types of studies, a set of doses may consist of doubled increments (e.g., 0, 5, 10, 20, 40, 80, and 100 g ae ha⁻¹). Ideally, doses should be administered in the following pattern: 0, 1×, 2×, 4×, 8×, 16×, 32×, 64×, and 128×, with x being the lowest tested dose. This uniform progression simplifies the interpretation of the dose-response curve and ensures that each dose level contributes evenly to the overall understanding of the response pattern, making the analysis more robust and precise.

Data Variables and Collection Timing

In dose-response studies aimed at investigating plant tolerance or hormesis effects, it is essential to report several key variables (responses) and the timing of data collection.

Data Variables

Dry matter (DM) and relative biomass serve as robust biological indicators for assessing plant growth (Knezevic et al. Reference Knezevic, Sikkema, Tardif, Hamill, Chandler and Swanton1998). DM can be reported as an absolute value or expressed relative to a nontreated control (a percent of the nontreated check). Plant height is also one of the common variables through which hormesis occurrence has been documented. Additional parameters can include growth and leaf stages, leaf area index (LAI), and harvest index, which may help explain how the stressors affect plant development. Moreover, LAI and DM data can be further used to calculate the net assimilation rate (NAR), as a useful measure of a plant’s photosynthetic efficiency (Sudhakar et al. Reference Sudhakar, Latha and Reddy2016). The NAR can be calculated with Equation 6:

([6]) $$NAR\left( {g/{m^2}/day} \right) = {{{W_2} - {W_1}} \over {{t_2} - {t_1}}} \times {{{\rm{ln}}\left( {{A_2}} \right) - {\rm{ln}}\left( {{A_1}} \right)} \over {{A_2} - {A_1}}}$$

where W ₁ and W ₂ are DM at times t ₁ and t ₂, while the ln(A ₁) and ln(A ₂) are the natural logarithms of leaf area at times t ₁ and t ₂. As indicated, to calculate NAR, two destructive harvests are required.

In studies that focus on a crop species, data related to yield and its components hold significant importance, particularly when studying crops grown for grain (corn, wheat, rice, soybean). For crops grown for biomass production, such as silage corn and leguminous plants including alfalfa, the DM measure assumes the role of yield. When the subject of investigation is a weed, DM and seed production can serve as an equivalent measure of yield.

Finally, percent visual injury is one of the most common data variable collected by weed scientists. Typicality, visual injury ratings are assigned using a scale from 0 to 100, where 0 signifies no injury and 100 indicates plant death, providing a basis to quickly quantify the extent of damage herbicide inflicted on the plant. However, visual injury ratings may be a biased assessment because they can vary greatly among researchers. On the other hand, visual injuries are usually the most visible plant response to herbicides, which is very important when estimating NOAEL values.

R Software and Data Analysis

R is a free and open-source software widely used for statistical computing, available for download from the CRAN website (R Core Team 2023). It provides a versatile platform for tasks such as data analysis, statistical modeling, and visualization. To work with R effectively, many analysts use RStudio, an integrated editor (development environment) that streamlines the entire data analysis process, available on the POSIT website (Posit Team Reference Team2023). RStudio serves as a user-friendly space where data manipulation, visualization, and statistical operations come together, making it an essential tool for data analysis professionals. Many useful packages can be used within RStudio, however, in this manuscript the focus will be on the drc (dose-response curve) package (Ritz et al. Reference Ritz, Baty, Streibig and Gerhard2015).

Data Organization and Input

Before loading data into RStudio, it is essential to ensure that the dataset is structured in a tidy manner, following the principles outlined by Wickham (Reference Wickham2014). Tidy data entails organizing the data so that each variable is represented as a separate column, each observation corresponds to a distinct row, and each cell is a single value, which typically occurs within an Excel spreadsheet (Table 1).

Table 1. Example of data organization in Excel.csv^a.

^a Abbreviation: .csv, comma-delimited file format.

In Table 1, the variables “Herbicide,” “Stage,” “Dose,” and “Yield” are clearly separated into distinct columns. Each row represents a specific observation, providing information about a single treatment as a unique combination of herbicide, growth stage, and dose, and observed yield as a response. Additionally, each cell within the table contains a single, discrete value, ensuring that the data are structured to facilitate efficient data analysis and interpretation. In an ideal experimental setup, it is recommended to have enough replicates to enhance the robustness and reliability of the analysis. Aiming for six to eight replicates is advantageous for achieving optimal statistical power and accuracy. However, recognizing practical constraints, a minimum of four replicates is recommended to ensure reasonable confidence in the observed trends and outcomes. Replicates, represented by the “Replicate” column in the dataset, contribute to the overall reliability of the study by accounting for variability and providing a basis for more comprehensive statistical analyses. The inclusion of “Plot” information further aids in distinguishing and tracking individual experimental units, contributing to the overall clarity and organization of the dataset.

Although data can be stored in various file formats, the preferred format is comma-delimited (.csv). There are several ways to load the data in RStudio. In this manuscript, we will focus on one specific method.

Case 1. Hormesis

When investigating potential herbicide hormesis effects, it is essential to obtain both visual (graphic) and statistical evidence (f parameter estimation) of the phenomenon. If the fit lacks one of the two pieces of evidence, hormetic effect should not be claimed.

The dataset used in the following three examples was gathered in field experiments concerning the effects of increasing doses of glyphosate (0, 224.7, 449.4, 898.8, 1,347.5, 1,796.5, 2,310, and 2,887.5 g ae ha⁻¹) on several weed species, including Amaranthus tuberculatus (Moq.) (AMATU) (Case Study 1a with data set AMATU), Convolvulus arvensis (L.) (CONAR) (Case Study 1b with data set CONAR), and Melilotus officinalis (L.) (MEUOF) (Case Study 1c with data set MEUOF) (Knezevic, unpublished data). Each experiment used a randomized complete block design. The observed response in each experiment was plant DM relative to an untreated check (100%) determined at 28 d after treatment. Specific details regarding the experimental site and procedures are omitted because the focus of this paper is not to discuss the biological interpretation of the results. Therefore, the following are the three examples of likely resulting scenarios involving determination of hormesis (case studies 1a, 1b, and 1c).

Case Study 1a. No Hormesis and No Increase in Response

We used the AMATU dataset to illustrate a scenario in which a no-hormesis and no visual increase in response was observed.

The initial step in the procedure involves loading the drc package by using the library function (Table 2, Step 1, Line 01). Line 02 assigns the object name of our choice to the data file. In this step, it is essential that the .csv file is located in the R project directory, otherwise an error message will occur. Line 03, function head, shows the first six rows of the imported dataset. Conversely, function tail (Table 2, Line 04) is used to display the last six rows of the dataset, which helps confirm that the data were correctly read into the working environment.

Table 2. R codes, outputs, and comments for Case Study 1a.^a

Once the data are loaded successfully, the next step involves curve fitting using the drm function (Table 2, Line 05). The fct argument within the drm function specifies the model of our choice to be fitted to the data. When hormesis is anticipated, either through preanalysis or examination of a scatterplot, it is advisable to initiate the fitting process with the BC.4 model. If in the latter steps BC.4 turns out to be a poor choice, it can be changed in the above-mentioned fct argument. Executing Line 05 will not produce a specific output. Line 06 invokes the lack-of-fit test. In the particular example of fitting the BC.4 model to our AMATU data, the lack-of-fit test resulted in a highly insignificant P-value (0.7860; i.e., no significant lack of fit), indicating that the data have been appropriately described by the selected model. Therefore, proceed by obtaining parameter estimates using the summary function (Table 2, Line 07). After obtaining parameter estimates and corresponding P-values, hormesis parameter f is estimated to be 2.2, with a highly insignificant P-value (0.8976) (Table 2). Despite the lack-of-fit test of the BC.4 hormesis model being highly insignificant, that is not enough proof of a hormesis effect, because the test and an indicator is the estimation of the f parameter (2.2, P = 0.8976) (Table 2).

Consequently, the next step is describing the data by the LL.4 model, again using the drm function (Table 2, Line 08). Note that the fct argument is now changed from BC.4 to LL.4. Once again, the lack-of-fit test is being performed (Table 2, Line 09), yielding a highly insignificant P-value (0.6911; Table 2), suggesting that there is no significant lack of fit, thus the chosen model fits the data well. Proceed by obtaining the parameter estimates by running the summary function (Table 2, Line 10). Finally, the visualization of both models is produced by the plot function (Table 2, Line 11; Figure 2). Additionally, a relative effective dose for a chosen level (1 to 99) can be obtained by invoking the ED function (Table 2, Line 12).

Figure 2. Dose-response curves of log-logistic (solid line) and Brain-Cousens (dashed line) models displayed together with the same AMATU [Amaranthus tuberculatus (Moq.)] dataset. The Brain-Cousens curve shows initial increase in response with no data points in region to support it. Log-logistic model displays an adequate fit to the data. Commands and equation parameters can be found in Table 2.

Based on statistical evidence, including the f parameter estimate, corresponding P-value, and a visual inspection of the plotted models, one can conclude that none of the tested glyphosate doses resulted in an increase in AMATU response. Therefore, hormesis did not occur.

Case Study 1b. No Hormesis, No Statistical Significance Despite Visual Evidence of Increase in Response

To illustrate a scenario in which a visual increase in response is observed without statistical hormesis conformation, the CONAR was used. Step 1 is the same as it was in Case Study 1a. In Step 2 (curve fitting), BC.4 is used similarly with the drm function (Table 3, Line 01). Like in the previous example, the lack-of-fit test is performed (P = 0.9782; Table 3, Line 02) followed by obtaining parameter estimates (Table 3, Line 03). As shown in the output Table 3, the f estimate is equal to 0.49334, with its P-value (0.4163334) indicating the lack of statistical evidence of f being significantly different than zero, suggesting no hormesis.

Table 3. R codes, outputs, and comments for Case Study 1b.^a

The next step involves further testing of some alternative hormesis models such as BC.5 or CRS.4 (or others). For example, the CRS.4b model was fitted (Table 3, Line 04), whose lack-of-fit test (Table 3, Line 05) showed an insignificant P-value (0.9726; Table 3), whereas the f parameter estimate (Table 3, Line 06) yielded a P-value of 0.207168 (Table 3), again suggesting the lack of statistical evidence for hormesis.

While both fitted hormesis models demonstrate highly insignificant P-values for the lack-of-fit test, confirming hormesis cannot be conclusively established. Despite the data being well described by one of the hormesis models and an apparent increase in response (Figure 3), this alone does not offer sufficient evidence to support the claim because the statistical test for hormesis requires f > 0. This example emphasizes the inherent challenges in interpreting hormesis, where a visual indication of increased response may not align with statistical confirmation. Consequently, it underscores the significance of considering both visual and statistical evidence in the analytical process of hormesis evaluation. It also serves as a reminder that modeling tools are only instruments in the hands of researchers, while the ultimate decision should be made through a comprehensive evaluation by weighing statistical metrics, graphical insights, and overarching research objectives.

Figure 3. Dose-response curves of Brain-Cousens (BC, solid line) and Cedergreen-Ritz-Streibig (CRS, dashed line) models displayed together on a CONAR [Convolvulus arvensis (L.)] dataset. The CRS curve overestimates the upper limit (untreated check response). The BC curve displays an adequate fit to the data. Commands and equation parameters can be found in Table 3.

Case Study 1c. Hormesis Confirmed with Both Statistical and Visual Evidence

To demonstrate a scenario in which both a visual increase in response is observed along with statistical evidence, we will use the MEUOF dataset. Step 1 is the same as in the previous two examples. Assuming no prior knowledge of the data, the LL.4 model is initially fitted (Table 4, Line 01), followed by the modelFit test (Table 4, Line 02). In this instance, a significant lack of fit is identified (P = 0.0465), indicating that the LL.4 model does not adequately describe the data, which is plotted for reference and visuals (Table 4, Line 10). The combination of statistical and visual evidence (Figure 4) makes it evident that the LL.4 model is not the suitable choice for this specific data.

Table 4. R codes, outputs, and comments for Case Study 1c.^a

Figure 4. Dose-response curves with four parameter log-logistic (solid line) and five parameter Brain-Cousens (dashed line) models displayed together on MEUOF [Melilotus officinalis (L.)] data set. The log-logistic curve overestimates the upper limit (untreated check response) and exhibits high deviation from the actual data points. The Brain-Cousens curve adequately fit the data. Commands and equation parameters can be found in Table 4.

The hormesis model of choice (BC.5) is fitted next (Table 4, Line 03). The modelFit test (Table 4, Line 04) suggests that there is no significant lack of fit (P = 0.8190; Table 4), indicating that the data is well described by the BC.5 model. Parameter estimates obtained (Table 4, Line 05) using the summary function, showed an estimate of the f parameter (0.48473) with a significant P-value (0.010647).

Furthermore, it is also useful to conduct another measure of comparison between the two models (LL.4 and BC.5) discussed above, which should be based on the Akaike information criterion (AIC) test (Table 4, Lines 06 and 07). The AIC values for the LL.4 model and BC.5 model were 259.279 and 248.9205 (Table 4), respectively. A lower AIC value suggests a better balance between model complexity and fit to the data. In this case, the BC.5 model exhibits a substantially lower AIC by more than 10, providing robust evidence to discard the LL.4 model in favor of the BC.5 model (Burnham and Anderson Reference Burnham and Anderson2004), signifying its capacity to describe the observed data patterns better.

Finally, to further confirm these findings, the anova function for final model comparisons (within the drc package) (Table 4, Line 08) can be used. This function evaluates whether a more complex model (BC.5) would provide a significantly better fit than the simpler model (LL.4). In essence, it tests the null hypothesis that the larger (more complex) model is not significantly better in describing data than the simpler one. The results strongly support the rejection of this null hypothesis, indicating that the inclusion of the hormesis parameter (f) contributes significantly to improving the model fit (P = 0.002; Table 4). Additionally, the mselect function (Table 4, Line 09) helps in model selection by evaluating several model fit criteria: maximum log likelihood value; AIC; estimated residual variance; and the P-value derived from a lack-of-fit test. As shown in the output (Table 4, Line 09) the BC.5 model exhibits the best performance across these metrics, boasting the highest log likelihood, lowest AIC value, lowest residual variance, and a nonsignificant lack-of-fit test P-value among the considered models. Therefore, based on these criteria, the BC.5 model emerges as the most suitable choice for the data. This further confirms our choice of the BC.5 model as the most suitable option for the dataset. However, although the mselect function is valuable, we do not advise relying solely on its output when modeling potential hormesis data. Modeling necessitates a comprehensive approach that goes beyond numerical metrics alone. It should incorporate the evaluation of different metrics, statistical tests, and visual inspection of plots. This holistic approach ensures a thorough understanding of the underlying patterns in the data and helps in making informed decisions regarding model selection and interpretation. The fitted model is visualized by employing the plot function (Table 4, Line 10). This confirms the occurrence of hormesis.

Therefore, based on the multiple levels of statistical analysis described above (Table 4) and visual examination of the fitted curve (Figure 4), glyphosate-induced hormesis is confirmed in MEUOF, which was expressed as an increase in relative DM. Finally, a similar procedure can be used for determining hormesis in other plant species, including crops of interest, with DM or crop yield as the response variable.

Case Study 2. Estimating Toxicological (Safety) Thresholds

The estimation of safety thresholds (NOAEL and LOAEL) in herbicide applications is a critical aspect of agricultural research, ensuring responsible and effective weed management. As indicated earlier, the NOAEL was arbitrarily assigned within the range of ED₁ to ED_2.5, accounting for general variability and potential errors. On the other hand, LOAEL is denoted as ED₅, considering traditional statistical significance (α = 0.05) and a practical acceptance level (e.g., for yield loss).

Case Study 2a. Estimating NOAEL

A NOAEL estimate is typically derived from the most sensitive response in plants, usually a visual estimate of injury. Choosing the most sensitive response (e.g., visual injury) ensures detection of adverse effects at the lowest possible dose, maintaining a conservative and protective approach in safety threshold determination.

Sample data were used to illustrate the safety threshold doses of dicamba drift on Roundup Ready soybean (Bayer Crop Science, St. Louis, MO). Ten doses were tested, including 0, 0.0112, 0.014, 0.019, 0.028, 0.056, 0.112, 0.56, 5.6 and 56 g ae ha⁻¹. Crop visual injury, as the most sensitive response, was assessed at 21 d after treatment (Knezevic unpublished data). For readability, doses are expressed in milligrams.

The procedure of loading in the data (Step 1) is the same as in other case studies. Subsequently, the LL.4 model is fitted (Table 5, Line 01), followed by the modelFit test (Table 5, Line 02). Parameter estimates are obtained by employing the summary function (Table 5, Line 03). The fitted model is plotted for inspection (Figure 5; Table 5, Line 04). Finally, estimates of the threshold doses, with corresponding 95% confidence intervals (CIs), are obtained using the ED function (Table 5, Line 05). Notice that the argument type is set to “relative” (Table 5, Line 05), which forces the ED to be calculated as a percent change in response (1%, 2.5%, and 5% in this case) between the estimated lower and upper limit (0 and 87, obtained from by applying the summary function). Estimated ED₁ is 0.32 (95% CI: 0.19 to 0.44) mg ae ha⁻¹, whereas ED_2.5 is 1.48 (95% CI: 0.90 to 2.05) mg ae ha⁻¹. Therefore, it can be concluded with 95% confidence that the dicamba dose ranging from 0.19 and 2.05 mg ae ha⁻¹ will cause no observable adverse effects to soybeans (at 21 d after treatment). It is common for the estimated upper limit of a curve not to reach 100%, which represents the maximal possible response. Therefore, deriving the ED values with the argument type set to “absolute” is often beneficial. This approach estimates the ED values as a percent change in response between 0 and 100, regardless of the estimated upper and lower limits, which may provide a better representation of the true NOAEL value. Now, estimated ED₁ is 0.35 (95% CI: 0.21 to 0.48) mg ae ha⁻¹, and estimated ED_2.5 is 1.62 (95% CI: 0.99 to 2.25) mg ae ha⁻¹ (Table 5, Line 06). In this case, the relative calculations of ED values are more conservative, as the estimated upper limit is 87%, which is lower than 100%.

Table 5. R codes, outputs, and comments for Case Study 2a.^a

Figure 5. Dose-response curve with four parameter log-logistic model (LL.4). The curve adequately fits the data. Commands and equation parameters can be found in Table 5.

Case Study 2b. Estimating LOAEL

From a practical standpoint, we propose estimating LOAEL using crop yield, at least with crops that produce grains as the final yield. When crops grown for biomass are the focus of research, DM serves as an equivalent of yield. Determining the lowest dose of a herbicide with an adverse effect on crop yield is crucial for farmers, extension educators, and weed science practitioners. Moreover, knowing the LOAEL values can be critically important when addressing drift complaints and for environmental and pesticide regulation agencies when establishing legislation, and for other stakeholders involved in ensuring the safety of chemical use in agricultural practices.

Data for this example were obtained in a field experiment aimed at studying the effects of clethodim tank contamination on subsequently treated corn (Knezevic unpublished data). The study used a randomized complete block design with four replications and tested nine different clethodim doses, including 0, 0.133, 0.268, 0.531, 1.062, 2.124, 8.494, 16.988, and 39.976 g ae ha⁻¹. The observed response was grain yield, expressed relative to an untreated check (100%).

Step 1 is the same (reading in the data) as in other case studies, and the curve fitting process is similar to that of previous examples. With a prior inspection of the data, we assume there is no initial increase in response and that the lower limit is not zero. Therefore, the drm model of choice is LL.4 (Table 5, Line 01), followed by the modelFit (Table 5, Line 02), which yielded an insignificant P-value (0.7136), suggesting an adequate description of data.

Proceed with parameter estimate (Table 5, Line 03). The fitted model is plotted for visual inspection (Table 5, Line 04; Figure 6). Finally, threshold effective dose values are estimated using the ED function (Table 5, Line 05). Note that the additional argument “interval = delta” specifies the 95% confidence interval estimation (delta-method, default 95%). Considering the very subtle changes in response (1% to 5%) is detected, it is useful to report the confidence intervals along with estimated ED values. In this case, the estimated ED₅ is calculated as 7.61 g ae ha⁻¹, with a corresponding 95% CI ranging from 7.28 to 7.94. Therefore, the results suggested with 95% confidence that the true lowest dose of clethodim that will cause a significant yield reduction (yield LOAEL) is between 7.28 and 7.94 g ae ha⁻¹. In other words, a dose lower than 7.28 g ae ha⁻¹ is not expected to cause a 5% yield reduction.

Figure 6. Dose-response curve with four parameter log-logistic model (LL.4). The curve adequately fits the data. Commands and equation parameters can be found in Table 6.

Table 6. R codes, outputs, and comments from Case Study 2b

It is important to note that the same procedures can be used to estimate the NOAEL and LOAEL values for any other responses, including visual injury, dry matter, height reduction, or any other growth or yield parameter of interest. This is achieved by specifying the respective variable within the drm function (Table 5, Line 02).

Practical Implications

The R software, coupled with the drc package, offers substantial value to users in the weed science community. In particular, the ability to estimate ED values, compare multiple curves and models, make the drc package a versatile tool for researchers. The R codes shared in this manuscript could foster collaboration, encourage dialogue and knowledge exchange among researchers. Users can leverage the codes as a basis for the analysis of their data, and to promote discussions, refinements, and improvements. Collaborative efforts driven by shared codes will contribute to the collective advancement of methodologies in the field of weed science.

Furthermore, the basic approach provided in this manuscript, which emphasized the enhancement of research in several dose-response concepts (e.g., hormesis) should help weed scientists move away from using traditional ANOVA approaches. The simple calculation of NOAEL and LOAEL, as demonstrated with the R codes provided here, will contribute to advancements in methodology and can be beneficial in environmental conservation and addressing ecotoxicological concerns.