Revisiting dose-response: concepts of hormesis, toxicological thresholds and data analysis

Abstract

Several publications and web-based tools are available in weed science literature to help weed scientists to carry out basic analysis of dose-response studies. Given the nature of the complicated relationship between the explanatory variable (dose as x-axis) evaluated against response variables of interest (y-axis), using regression curves should be the preferred method for handling data analysis. The objective of this manuscript is to provide user-friendly instructions for conducting and analyzing several types of dose-response studies that were lacking in current weed science literature. A better understanding of less commonly used concepts of hormesis and toxicological safety thresholds (no-observable-adverse-effect-level [NOAEL] and lowest-observable-adverse-effect-level [LOAEL]) is needed to help address the potential risks and benefits associated with herbicide use while minimizing environmental impacts. Basic codes available in cost-free R software are provided for data analysis and to foster collaboration among weed scientists.

Introduction

Sustainable agriculture aims to align food production with responsible environmental practices. Integrated weed management stands out as a central strategy in the pursuit of sustainability (Zimdahl 2018). As part of an integrated weed management system, herbicides play a crucial role in providing essential solutions to address various weed control challenges. Nevertheless, applying herbicides and pesticides in general always involves potential safety concerns that need to be carefully considered and properly addressed.

One of the common safety concerns is the unintended off-target movement of herbicides, known as drift. Notably, since the introduction of dicamba-tolerant cotton and soybeans in 2017, the incidence of dicamba off-target movement has surpassed any previous records in the history of U.S. agriculture (Bish et al. 2021). For instance, during the initial growing season of dicamba-tolerant cotton and soybean, state departments of agriculture conducted 2,708 nationwide investigations into alleged dicamba-induced crop injuries (Oseland et al. 2020).

Drift is not the sole pathway through which herbicides can reach unintended vegetation. Herbicide residues may persist in application equipment, including hoses and tanks, posing an ongoing risk of harming subsequently treated crops (Batts et al. 2022). Anticipated trends in weed management suggest an increasing reliance on the development of crop varieties engineered to tolerate a range of herbicides. While multiple herbicide–tolerant crops show potential for improving weed control, the widespread use of these crops raises concerns about herbicide drift and tank contamination.

On the other hand, the sublethal doses that are present in herbicide drift may not always have adverse effects on the nontargeted plants. It is well known that many herbicides have beneficial effects at low (sublethal) doses, yet they become detrimental at higher doses through a dose-response phenomenon called hormesis (Cedergreen 2008a). Belz et al. (2018) expressed growing concerns that the presence of variable hormesis responsiveness may lead to the selection of herbicide-resistant weeds. Given that herbicide-resistant weed populations represent a significant issue stemming from the widespread use of specific herbicides, there should be an increased interest in studying the effects of sublethal herbicide doses on weeds. On the other hand, Agathokleous and Calabrese (2019) argue for the use of hormesis in agriculture as a means to enhance crop productivity, thereby increasing global food supplies and fostering socioeconomic development.

Therefore, it is useful to determine whether exposure to sublethal doses of certain herbicides can induce hormesis in various plant species. It is also of interest to determine the toxicological thresholds, such as the no observed adverse effect level (NOAEL) and the lowest observed adverse effect level (LOAEL).

Despite various publications (Knezevic et al. 2007) and web-based tools available on the subject, the existing weed science literature lacks comprehensive guidelines on how to conduct and analyze dose-response studies, particularly those intended to assess hormesis and toxicological thresholds (NOAEL and LOAEL) of herbicides on both crops and weeds of interest. Hence, the objectives of this manuscript are to provide 1) general recommendation for studying the above-listed concepts, and 2) step-by-step guidance and R statistical software codes for data analysis within the drc package (Ritz et al. 2015).

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 2008). 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 2014). 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. 2011), 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 1972), soybean [Glycine max (L.) Merr.], corn (Zea mays L.), and eucalyptus (Eucalyptus grandis Hill ex Maiden) (Velini et al. 2008). 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 2014). 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. 2020).

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 2014). 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, 2007), 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 2018), 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. 2023). 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. 2007). 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. 2007). 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. 2007). 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. 2007; Van der Vilet and Ritz 2013). 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. 2005) (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. 2005).

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. 2007). 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. 2007). 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. 1998). 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. 2016). 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.

Timing of Data Collection

The timing of data collection is a critical aspect of any experiment and typically aligns with the specific hypotheses and objectives of the study. For instance, in herbicide evaluation trials conducted within a single growing season, visual ratings of percent weed control are routinely gathered at intervals of 1, 2, 4, 8, and 12 wk after treatment (WAT) (Knezevic et al. 2007). However, when the occurrence of herbicide hormesis in plants is reported, it is typically observed at single time points (Cedergreen 2008b) or within a short time window (e.g., 2 to 4 wk) after stimulus exposure (herbicide application). Therefore, it is important that plant response data should be collected weekly for at least 4 WAT. This will allow curve fitting and curve comparisons for each weekly response and perhaps calculate growth rates (Equation 6) for that specific growth period.

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 2023). 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. 2015).

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 (2014). 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 2004), 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.