Abstract

The knowledge of accumulated evapotranspiration by seasonal vegetation crops throughout their life cycle can be an important tool in decision-making when considering the economic viability of the crop. This knowledge can help understand how much the plants, subject to specific management, can evapotranspirate at the end of their cycle. This information assists in estimating the quantity of a production variable, for example, the mass of shoot fresh matter, besides indicating a more interesting period for its harvest. The objective of this chapter is, from the daily evapotranspiration estimative throughout the cycle, to model the accumulated evapotranspiration over the entire growth period of the crop. In order to do so, we must understand that the behavior of the response variable, i.e., the accumulated evapotranspiration, over time is not linear and keep in mind that the several observations performed in the same experimental unit have correlations and these correlations are more intense the closer temporally the measurements are. This understanding leads us to the analysis of longitudinal data from the nonlinear mixed effect models perspective.

Keywords: longitudinal data, nonlinear mixed effect model, growth curve, correlation structure, irrigation

### 1. Introduction

With available water in the soil, the water flow through the plants depends only on atmospheric demand. Therefore, physical variables as temperature, relative air humidity, and wind and solar radiation affect directly the evapotranspiration (ET) of a vegetated surface [1]. Besides that, the plants development state may also affect the ET.

Seasonal vegetation crops present a small demand for water, while their root system is small, reaching maximum rates in full growths and decreasing in the final stages of development. For these species, the accumulated ET over the plant cycle is directly related to its productivity. In other words, the greater the accumulated amount of evapotranspirated water, the greater the quantity of the production variables as shoot fresh and dry matter masses, number of leaves, and leaf area [2, 3].

The objective of this chapter is to illustrate the use of nonlinear mixed effect models to fit accumulated ET over time in seasonal vegetation crops.

#### 2. Response profile

Figure 1 presents the accumulated ET from a single lettuce plant over 23 consecutive days, being the first day equivalent to the 35th day after seeding. Throughout this period, the plant's daily ET was measured. Hence, the first day in the graph presents the ET value from the last 24 h, the second day refers to the last 48 h, and so on. The daily ET must be understood as the rate at which the accumulated ET occurs over time. In mathematical terms, the daily ET can be understood as an approximation of the derivative (rate) of the accumulated ET with respect to time.

statistical adjustments and be computationally simpler but does not add any theoretical consideration to the physical and/or biological mechanisms that generated the data. On the other hand, a nonlinear model is associated with some theoretical knowledge regarding the studied phenomenon. Besides the interpretability, these models use few parameters when compared to linear models, thus configuring a

Regarding the accumulated ET over time, we have some physical and biological aspects that we can use to choose a proper model. A function that describes this response variable needs parameters that delimit it between a minimum and a maximum value. In other words, the minimum value should be really close to zero and indicate the beginning of the plant's growth, and the maximum should be a value to which the accumulated ET tends asymptotically as the end of the life cycle approaches. Another aspect to highlight is that the model should present an inflection point, which indicates the day that the accumulated ET rate (daily ET) reaches

There are several functions capable of characterizing the accumulated ET over time in the sense we just described. For example, we will use the four-parameter logistic function (4PL) to describe our response variable. This function is widely used to fit growth or decay data. There are some parametrizations for this function

with ET tð Þ being the accumulated ET over time and ϕ<sup>1</sup>�<sup>4</sup> the model parameters (Figure 2). ϕ<sup>1</sup> is the inferior horizontal asymptote which gives the accumulated ET value when t ! �∞. Biologically, this parameter does not have a consistent interpretation, but it is important in the fitting because it ensures the accumulated ET, in times close to zero, to be very small [2, 3]. ϕ<sup>2</sup> is the superior horizontal asymptote

<sup>1</sup> <sup>þ</sup> exp ð Þ <sup>ϕ</sup><sup>3</sup> � <sup>t</sup> <sup>=</sup>ϕ<sup>4</sup> ½ � (1)

ET tðÞ¼ <sup>ϕ</sup><sup>1</sup> <sup>þ</sup> <sup>ϕ</sup><sup>2</sup> � <sup>ϕ</sup><sup>1</sup>

and gives the accumulated ET value when t ! ∞. This parameter can be

Graphical representation of the 4PL parameters (Eq. (1)). Figure adapted from [4].

more parsimonious description of the data [4].

Modeling Accumulated Evapotranspiration Over Time DOI: http://dx.doi.org/10.5772/intechopen.86913

in the literature, but we are using the one given by [4]

its maximum value.

Figure 2.

71

We wish to describe the behavior of our variable of interest or response variable, the accumulated ET over time for seasonal vegetation crops. As the amount of evapotranspirated water in 1 day is added to the accumulated ET from previous days, its values over time are equal (in case ET from a whole day is null) or greater than the immediately preceding value.

Observations from the response variable in more than one moment in the same experimental unity constitute what we call as response profile. Therefore, Figure 1 presents the response profile of a lettuce plant over time. This profile, apparently, presents an S-shaped format.

At the first days of observations (Figure 1), the accumulated ET is small because the plant is at the beginning of its growth. As the days go by, the daily ET values increase successively approximately until the 20th day. After the 20th day, the daily ET, i.e., the rate of change of the accumulated ET, begins to decrease. This causes the accumulated ET, which has been growing exponentially, to have a less vigorous growth and, therefore, tending slowly to a maximum value.

Figure 1. Response profile of the accumulated ET for a single lettuce plant over 23 consecutive days.

#### 3. A model for growth data

In order to describe the behavior of the accumulated ET over time, we also need an S-shaped function. In addition, it is expected that the chosen model will be interpretable and parsimonious in the parameters. An empirical option is the polynomial model, which has linear parameters. This type of model can promote good

Modeling Accumulated Evapotranspiration Over Time DOI: http://dx.doi.org/10.5772/intechopen.86913

2. Response profile

Water Chemistry

than the immediately preceding value.

growth and, therefore, tending slowly to a maximum value.

presents an S-shaped format.

3. A model for growth data

Figure 1.

70

Figure 1 presents the accumulated ET from a single lettuce plant over 23 consecutive days, being the first day equivalent to the 35th day after seeding. Throughout this period, the plant's daily ET was measured. Hence, the first day in the graph presents the ET value from the last 24 h, the second day refers to the last 48 h, and so on. The daily ET must be understood as the rate at which the accumulated ET occurs over time. In mathematical terms, the daily ET can be understood as an approximation of the derivative (rate) of the accumulated ET with respect to time. We wish to describe the behavior of our variable of interest or response variable,

the accumulated ET over time for seasonal vegetation crops. As the amount of evapotranspirated water in 1 day is added to the accumulated ET from previous days, its values over time are equal (in case ET from a whole day is null) or greater

Observations from the response variable in more than one moment in the same experimental unity constitute what we call as response profile. Therefore, Figure 1 presents the response profile of a lettuce plant over time. This profile, apparently,

At the first days of observations (Figure 1), the accumulated ET is small because the plant is at the beginning of its growth. As the days go by, the daily ET values increase successively approximately until the 20th day. After the 20th day, the daily ET, i.e., the rate of change of the accumulated ET, begins to decrease. This causes the accumulated ET, which has been growing exponentially, to have a less vigorous

In order to describe the behavior of the accumulated ET over time, we also need

an S-shaped function. In addition, it is expected that the chosen model will be interpretable and parsimonious in the parameters. An empirical option is the polynomial model, which has linear parameters. This type of model can promote good

Response profile of the accumulated ET for a single lettuce plant over 23 consecutive days.

statistical adjustments and be computationally simpler but does not add any theoretical consideration to the physical and/or biological mechanisms that generated the data. On the other hand, a nonlinear model is associated with some theoretical knowledge regarding the studied phenomenon. Besides the interpretability, these models use few parameters when compared to linear models, thus configuring a more parsimonious description of the data [4].

Regarding the accumulated ET over time, we have some physical and biological aspects that we can use to choose a proper model. A function that describes this response variable needs parameters that delimit it between a minimum and a maximum value. In other words, the minimum value should be really close to zero and indicate the beginning of the plant's growth, and the maximum should be a value to which the accumulated ET tends asymptotically as the end of the life cycle approaches. Another aspect to highlight is that the model should present an inflection point, which indicates the day that the accumulated ET rate (daily ET) reaches its maximum value.

There are several functions capable of characterizing the accumulated ET over time in the sense we just described. For example, we will use the four-parameter logistic function (4PL) to describe our response variable. This function is widely used to fit growth or decay data. There are some parametrizations for this function in the literature, but we are using the one given by [4]

$$ET(t) = \phi\_1 + \frac{\phi\_2 - \phi\_1}{1 + \exp\left[ (\phi\_3 - t) / \phi\_4 \right]} \tag{1}$$

with ET tð Þ being the accumulated ET over time and ϕ<sup>1</sup>�<sup>4</sup> the model parameters (Figure 2). ϕ<sup>1</sup> is the inferior horizontal asymptote which gives the accumulated ET value when t ! �∞. Biologically, this parameter does not have a consistent interpretation, but it is important in the fitting because it ensures the accumulated ET, in times close to zero, to be very small [2, 3]. ϕ<sup>2</sup> is the superior horizontal asymptote and gives the accumulated ET value when t ! ∞. This parameter can be

Figure 2. Graphical representation of the 4PL parameters (Eq. (1)). Figure adapted from [4].

interpreted as the maximum accumulated ET estimative that the plant can reach in its final life cycle. ϕ<sup>3</sup> is the curve's inflection point and indicates the day (time) in which the daily ET reaches a maximum value. The corresponding time to ϕ<sup>3</sup> results in an accumulated ET between ϕ<sup>1</sup> and ϕ2. More precisely, the accumulated ET until time ϕ<sup>3</sup> is the mean value of both asymptotes, i.e., ϕ<sup>1</sup> þ ϕ<sup>2</sup> ð Þ=2. ϕ<sup>4</sup> is the scale parameter. The day (time) corresponding to ϕ<sup>3</sup> þ ϕ<sup>4</sup> gives � 0:75 ϕ<sup>2</sup> � ϕ<sup>1</sup> ð Þ of the accumulated ET. Therefore, ϕ<sup>4</sup> indicates how quickly the accumulated ET leaves the proximity of ϕ<sup>1</sup> until it reaches values close to ϕ2. The greater the ϕ4, the slower this occurs.

plants' life cycle. That is, the larger the total evapotranspiration water, the higher the values of the production variables [2, 3]. Hence, the inflection point can indicate the proximity of the harvest day. From this point, the ET rates decrease, and therefore, the plants start to accumulate less quantities of shoot fresh matter over time. Thus, this point can be considered when thinking in terms of the economic

Another important aspect, besides the parameter interpretability, is that the nonlinear 4PL model is more parsimonious than the linear model. In general, nonlinear models use a small number of parameters than the linear ones, which grant them more parsimony. Besides that, as can be observed in Figure 3, in regions outside the data interval, the nonlinear model gives a more trustworthy prediction

Clearly, in order to make sense, an experiment must provide data from more than one experimental unit. In our case, more than one plant should be observed over time, i.e., we must obtain more than one response profile. Studies in which the response variable is observed repeatedly throughout time in the experimental units are called longitudinal studies. This kind of work is common in agriculture when analyzing the increase or decrease of the response variable over time [2, 3, 5–8]. Measurements performed in the same experimental unit are most likely to be correlated. Suppose two plants which its ET is registered daily. If all covariables (fertilization, cultivate, planting season, soil water, and so on) were kept constant over time, plants with high rates of ET in a given day will most likely also have high rates of ET in the next day, the same for plants with smaller ET rates. If in a day the rate is small, probably in the next day it will also be small. In other words, measurements performed in the same experimental unit tend to be similar over time. It's

Besides the correlation between the observations within the same experimental unit, we must consider that, most likely, these correlations are greater for observations performed between neighboring times than those performed between more

In a longitudinal study, the monitoring of the experimental units over time generates correlated dataset. As mentioned, these correlations within the same experimental unit are stronger among neighboring observations. The greater the time distance between two measures, the weaker the correlation between them. Besides that, when we observe experimental units which received the same conditions regarding growth over time and are part of the same treatment, we have a variability among them that we attribute to chance. The treatment effects, the correlations, and the variabilities in a longitudinal study indicate that we need a tool that, in addition to being flexible in specifying a mathematical model, also empha-

In mixed effect models, we select an ordinary function to describe the response variable regarding the studied covariables, that is, the responses of the experimental units in a population. Besides that, specific coefficients of this function can be unique for each experimental unit. In a mixed effect model, we assume that the

viability of the crop.

for the response variable [4].

Modeling Accumulated Evapotranspiration Over Time DOI: http://dx.doi.org/10.5772/intechopen.86913

the individual expression of each plant.

distant times.

5. Mixed effect model

sizes each experimental unit.

73

4. Longitudinal data

To better understand what we have said, let's observe Figure 3. It brings the observed data shown in Figure 1 represented by the black dots and two other curves representing two different fits made from this data. The first is a fifth-degree polynomial model, and the second is given by Eq. (1). The graph in this figure was extended until the 40th day in order to present the behavior of these fittings over time.

Both the fifth-degree polynomial and the nonlinear 4PL models have fitted well with the data. However, the polynomial model goes to zero after the 30th day, which physically is impossible. Besides that, this model presents five parameters without any physical and/or biological explanation for the phenomenon in this study. The 4PL nonlinear model is a strictly increasing function, being compatible with a variable which is accumulated over time. Nevertheless, this fitting does not allow the infinite growth of the accumulated ET as time grows larger and larger. It limits the accumulated ET to a value which can be understood as a maximum amount in which this plant can evapotranspirate throughout its life cycle.

We also see that, from approximately the 20th day, the accumulation of evapotranspiration water is decreasing. This indicates that this day is the inflection point of the nonlinear 4PL model. For a practical example of the model's inflection point regarding seasonal vegetation crops, suppose that the commercial product of a given crop plant is its leaves. It is known that the production variables present a positive correlation with the amount of total evapotranspired water throughout the

#### Figure 3.

Accumulated ET response profile for a single lettuce plant over 23 consecutive days. The graph was extended until the 40th day to present the behavior of two fittings made for this data, a fifth-degree polynomial and a nonlinear model given by Eq. (1).

Modeling Accumulated Evapotranspiration Over Time DOI: http://dx.doi.org/10.5772/intechopen.86913

plants' life cycle. That is, the larger the total evapotranspiration water, the higher the values of the production variables [2, 3]. Hence, the inflection point can indicate the proximity of the harvest day. From this point, the ET rates decrease, and therefore, the plants start to accumulate less quantities of shoot fresh matter over time. Thus, this point can be considered when thinking in terms of the economic viability of the crop.

Another important aspect, besides the parameter interpretability, is that the nonlinear 4PL model is more parsimonious than the linear model. In general, nonlinear models use a small number of parameters than the linear ones, which grant them more parsimony. Besides that, as can be observed in Figure 3, in regions outside the data interval, the nonlinear model gives a more trustworthy prediction for the response variable [4].

## 4. Longitudinal data

interpreted as the maximum accumulated ET estimative that the plant can reach in its final life cycle. ϕ<sup>3</sup> is the curve's inflection point and indicates the day (time) in which the daily ET reaches a maximum value. The corresponding time to ϕ<sup>3</sup> results in an accumulated ET between ϕ<sup>1</sup> and ϕ2. More precisely, the accumulated ET until time ϕ<sup>3</sup> is the mean value of both asymptotes, i.e., ϕ<sup>1</sup> þ ϕ<sup>2</sup> ð Þ=2. ϕ<sup>4</sup> is the scale parameter. The day (time) corresponding to ϕ<sup>3</sup> þ ϕ<sup>4</sup> gives � 0:75 ϕ<sup>2</sup> � ϕ<sup>1</sup> ð Þ of the accumulated ET. Therefore, ϕ<sup>4</sup> indicates how quickly the accumulated ET leaves the proximity of ϕ<sup>1</sup> until it reaches values close to ϕ2. The greater the ϕ4, the slower

To better understand what we have said, let's observe Figure 3. It brings the observed data shown in Figure 1 represented by the black dots and two other curves representing two different fits made from this data. The first is a fifth-degree polynomial model, and the second is given by Eq. (1). The graph in this figure was extended until the 40th day in order to present the behavior of these fittings over

Both the fifth-degree polynomial and the nonlinear 4PL models have fitted well with the data. However, the polynomial model goes to zero after the 30th day, which physically is impossible. Besides that, this model presents five parameters without any physical and/or biological explanation for the phenomenon in this study. The 4PL nonlinear model is a strictly increasing function, being compatible with a variable which is accumulated over time. Nevertheless, this fitting does not allow the infinite growth of the accumulated ET as time grows larger and larger. It limits the accumulated ET to a value which can be understood as a maximum amount in

We also see that, from approximately the 20th day, the accumulation of evapotranspiration water is decreasing. This indicates that this day is the inflection point of the nonlinear 4PL model. For a practical example of the model's inflection point regarding seasonal vegetation crops, suppose that the commercial product of a given crop plant is its leaves. It is known that the production variables present a positive correlation with the amount of total evapotranspired water throughout the

Accumulated ET response profile for a single lettuce plant over 23 consecutive days. The graph was extended until the 40th day to present the behavior of two fittings made for this data, a fifth-degree polynomial and a

which this plant can evapotranspirate throughout its life cycle.

this occurs.

Water Chemistry

time.

Figure 3.

72

nonlinear model given by Eq. (1).

Clearly, in order to make sense, an experiment must provide data from more than one experimental unit. In our case, more than one plant should be observed over time, i.e., we must obtain more than one response profile. Studies in which the response variable is observed repeatedly throughout time in the experimental units are called longitudinal studies. This kind of work is common in agriculture when analyzing the increase or decrease of the response variable over time [2, 3, 5–8].

Measurements performed in the same experimental unit are most likely to be correlated. Suppose two plants which its ET is registered daily. If all covariables (fertilization, cultivate, planting season, soil water, and so on) were kept constant over time, plants with high rates of ET in a given day will most likely also have high rates of ET in the next day, the same for plants with smaller ET rates. If in a day the rate is small, probably in the next day it will also be small. In other words, measurements performed in the same experimental unit tend to be similar over time. It's the individual expression of each plant.

Besides the correlation between the observations within the same experimental unit, we must consider that, most likely, these correlations are greater for observations performed between neighboring times than those performed between more distant times.

#### 5. Mixed effect model

In a longitudinal study, the monitoring of the experimental units over time generates correlated dataset. As mentioned, these correlations within the same experimental unit are stronger among neighboring observations. The greater the time distance between two measures, the weaker the correlation between them. Besides that, when we observe experimental units which received the same conditions regarding growth over time and are part of the same treatment, we have a variability among them that we attribute to chance. The treatment effects, the correlations, and the variabilities in a longitudinal study indicate that we need a tool that, in addition to being flexible in specifying a mathematical model, also emphasizes each experimental unit.

In mixed effect models, we select an ordinary function to describe the response variable regarding the studied covariables, that is, the responses of the experimental units in a population. Besides that, specific coefficients of this function can be unique for each experimental unit. In a mixed effect model, we assume that the

experimental units of a population have the same functional form, but the function parameters may vary among the units.

they are distributed the same way among the experimental units. For a given sample

where each vector ε<sup>i</sup> has ni observations throughout time and we assume that they follow a normal multivariate distribution with mean zero and covariance

Most of the times, we consider <sup>Λ</sup><sup>i</sup> <sup>¼</sup> <sup>σ</sup><sup>2</sup>I, being <sup>σ</sup><sup>2</sup> a constant variance for all <sup>j</sup>

By using nonlinear mixed effect models, we must consider the technical difficulty in the parameter estimation. In a mixed effects linear model, the derivative of the logarithmic of the likelihood function allows, in a simple way, the algorithm implementation, like Newton-Raphson, to obtain the estimative of the models parameters. Nonlinear models can, however, present nonlinear random coefficients which make it impossible to directly explain the parameters from the likelihood function. Methods that depend on linear approximations such as the first-order

Nonlinear mixed effect model analysis can be, preferably, made by the R soft-

ware [13] with the package name [14], also, at the SAS software using PROC NLMIXED. An excellent text to learn how to use these skills is given in [4].

Mixed effect model allows the dependence between the observations to be specified in the model parameters through random effects. In other words, the experimental unit responses from a population tend to follow a nonlinear growth path; however, each experimental unit has its own growth path, and the mixed effect model allows the inclusion of specific coefficients to obtain fitted

Thus, mixed models allow relevant flexibility for the specification of the random effects correlation structure. However, the dependence structure of the observations within the experimental units Λ<sup>i</sup> until now has been considered independent, identically distributed with mean zero and constant variance. Depending on the chosen model, the growth responses can be explained just by including specific coefficients for the experimental units. However, this may not be enough, and, in this case, modeling the residual dependence of the data becomes important.

There are cases where dependence on observations not accommodated by the growth function is not well understood or, sometimes, additional covariables that could explain this dependences are absent from the model. Thus, an important resource to model this dependence is to identify the covariance structure that allows correlation between the residuals in different occasions. Then, let us relax on the assumption that the errors are independent and allow them to have heterosce-

growth curves that align better with the individual responses of these

dasticity and/or are correlated within the experimental units.

Taylor approximation can be used to estimate the model.

6. Covariance structure of Λ<sup>i</sup>

experimental units.

75

ε1i ε2i ⋮ εnii

ε<sup>i</sup> � Nð Þ 0;Λ<sup>i</sup> : (6)

unit i, we can describe the error εij as the time j indexed vector, i.e.,

Modeling Accumulated Evapotranspiration Over Time DOI: http://dx.doi.org/10.5772/intechopen.86913

matrix Λi, i.e.,

times.

ε<sup>i</sup> ¼

The name mixed model comes from the fact that this model combines fixed effects and random effects. A mixed effect model is a parametric model which describes the relations between the response variable and the covariables (fixed effects) and takes into account the individual responses of each experimental unit (random effects). In other words, the fixed effects parameters describe the relations of the response variable and the covariables in an entire population, and the random effects specify the contribution of each individual within the population [4, 9–12].

To illustrate how to write a 4PL nonlinear model with mixed effect, let's assume that we are studying the accumulated ET regarding four levels of water in the soil, W1, W2, W3, and W<sup>4</sup> (there could be more levels of water in the soil or less, but to exemplify, let's assume four treatments). We will consider that the response variable, the accumulated ET, has a normal probability distribution. Consider ETij as the accumulated ET at the j situation, for plant i, with j ¼ 1, 2, ⋯, ni e i ¼ 1, 2, ⋯, N, where ni is the number of observations for the i-th plant and N is the total number of plants. The nonlinear 4PL mixed effect model can be expressed by [2, 3]

$$ET\_{\vec{\eta}} = \phi\_{1i} + \frac{\phi\_{2i} - \phi\_{1i}}{1 + \exp\left[\left(\phi\_{3i} - t\_{\vec{\jmath}}\right) / \phi\_{4i}\right]} + \varepsilon\_{\vec{\eta}} \tag{2}$$

being

$$\boldsymbol{\Phi}\_{i} = \begin{bmatrix} \boldsymbol{\Phi}\_{1i} \\ \boldsymbol{\Phi}\_{2i} \\ \boldsymbol{\Phi}\_{3i} \\ \boldsymbol{\Phi}\_{4i} \end{bmatrix} = \begin{bmatrix} \boldsymbol{\beta}\_{1} + \boldsymbol{\gamma}\_{1}\mathbf{x}\_{1i} + \boldsymbol{\delta}\_{1}\mathbf{x}\_{2i} + \boldsymbol{\zeta}\_{1}\mathbf{x}\_{3i} \\ \boldsymbol{\beta}\_{2} + \boldsymbol{\gamma}\_{2}\mathbf{x}\_{1i} + \boldsymbol{\delta}\_{2}\mathbf{x}\_{2i} + \boldsymbol{\zeta}\_{2}\mathbf{x}\_{3i} \\ \boldsymbol{\beta}\_{3} + \boldsymbol{\gamma}\_{3}\mathbf{x}\_{1i} + \boldsymbol{\delta}\_{3}\mathbf{x}\_{2i} + \boldsymbol{\zeta}\_{3}\mathbf{x}\_{3i} \\ \boldsymbol{\beta}\_{4} + \boldsymbol{\gamma}\_{4}\mathbf{x}\_{1i} + \boldsymbol{\delta}\_{4}\mathbf{x}\_{2i} + \boldsymbol{\zeta}\_{4}\mathbf{x}\_{3i} \end{bmatrix} + \begin{bmatrix} \boldsymbol{b}\_{1i} \\ \boldsymbol{b}\_{2i} \\ \boldsymbol{b}\_{3i} \\ \boldsymbol{b}\_{4i} \end{bmatrix} \tag{3}$$

$$= \boldsymbol{\mathfrak{f}} + \boldsymbol{\mathfrak{y}}\mathbf{x}\_{1i} + \boldsymbol{\delta}\mathbf{x}\_{2i} + \boldsymbol{\xi}\mathbf{x}\_{3i} + \mathbf{b}\_{i}, \tag{4}$$

with the parameters β, γ, δ, ζ representing the fixed effects and b<sup>i</sup> the random effects in the model. xki with k ¼ 1; 2; 3 are indicative covariables of treatments or groups and may have values zero or one. The parameter β is the reference level in the study. When x1<sup>i</sup> ¼ x2<sup>i</sup> ¼ x3<sup>i</sup> ¼ 0, thus, ϕ<sup>i</sup> ¼ β þ b<sup>i</sup> and the i-th plant belongs to the treatment W1. When x1<sup>i</sup> ¼ 1 and x2<sup>i</sup> ¼ x3i ¼ 0, ϕ<sup>i</sup> ¼ β þ γ þ b<sup>i</sup> and the i-th plant belongs to the treatment W2. If x2<sup>i</sup> ¼ 1 and x1<sup>i</sup> ¼ x3i ¼ 0, ϕ<sup>i</sup> ¼ β þ δ þ b<sup>i</sup> and the i-th plant belongs to the treatment W3. And lastly, if x3<sup>i</sup> ¼ 1 and x1<sup>i</sup> ¼ x2<sup>i</sup> ¼ 0, ϕ<sup>i</sup> ¼ β þ ζ þ b<sup>i</sup> and the i-th plant belong to the treatment W4. The random effects b<sup>i</sup> are considered independent among the plants and are normally distributed with mean zero and covariance matrix Ψ (b<sup>i</sup> � Nð Þ 0; Ψ ). In this case, the covariance matrix is given by

$$
\Psi = \begin{bmatrix}
\sigma\_1^2 & \sigma\_{12} & \sigma\_{13} & \sigma\_{14} \\
\sigma\_{12} & \sigma\_2^2 & \sigma\_{23} & \sigma\_{24} \\
\sigma\_{13} & \sigma\_{23} & \sigma\_3^2 & \sigma\_{34} \\
\sigma\_{14} & \sigma\_{24} & \sigma\_{34} & \sigma\_4^2
\end{bmatrix} \tag{5}
$$

with σ<sup>2</sup> 1, σ<sup>2</sup> 2, σ<sup>2</sup> 3, and σ<sup>2</sup> <sup>4</sup> being the variances of random effects b1<sup>i</sup>, b2<sup>i</sup>, b3<sup>i</sup>, and b4<sup>i</sup>, respectively, and σ12, σ13, σ14, σ23, σ24, and σ<sup>34</sup> are the covariances between them. The variances of the random effects indicate how a model parameter varies between the experimental units. Frequently, we suppose that the errors within the groups εij are independent between the observations of the same experimental unit and that

experimental units of a population have the same functional form, but the function

the contribution of each individual within the population [4, 9–12].

ETij <sup>¼</sup> <sup>ϕ</sup>1<sup>i</sup> <sup>þ</sup> <sup>ϕ</sup>2<sup>i</sup> � <sup>ϕ</sup>1<sup>i</sup>

The name mixed model comes from the fact that this model combines fixed effects and random effects. A mixed effect model is a parametric model which describes the relations between the response variable and the covariables (fixed effects) and takes into account the individual responses of each experimental unit (random effects). In other words, the fixed effects parameters describe the relations of the response variable and the covariables in an entire population, and the random effects specify

To illustrate how to write a 4PL nonlinear model with mixed effect, let's assume that we are studying the accumulated ET regarding four levels of water in the soil, W1, W2, W3, and W<sup>4</sup> (there could be more levels of water in the soil or less, but to exemplify, let's assume four treatments). We will consider that the response variable, the accumulated ET, has a normal probability distribution. Consider ETij as the accumulated ET at the j situation, for plant i, with j ¼ 1, 2, ⋯, ni e i ¼ 1, 2, ⋯, N, where ni is the number of observations for the i-th plant and N is the total number of plants. The nonlinear 4PL mixed effect model can be expressed by [2, 3]

1 þ exp ϕ3<sup>i</sup> � tj

β<sup>1</sup> þ γ1x1<sup>i</sup> þ δ1x2<sup>i</sup> þ ζ1x3<sup>i</sup> β<sup>2</sup> þ γ2x1<sup>i</sup> þ δ2x2<sup>i</sup> þ ζ2x3<sup>i</sup> β<sup>3</sup> þ γ3x1<sup>i</sup> þ δ3x2<sup>i</sup> þ ζ3x3<sup>i</sup> β<sup>4</sup> þ γ4x1<sup>i</sup> þ δ4x2<sup>i</sup> þ ζ4x3<sup>i</sup>

with the parameters β, γ, δ, ζ representing the fixed effects and b<sup>i</sup> the random effects in the model. xki with k ¼ 1; 2; 3 are indicative covariables of treatments or groups and may have values zero or one. The parameter β is the reference level in the study. When x1<sup>i</sup> ¼ x2<sup>i</sup> ¼ x3<sup>i</sup> ¼ 0, thus, ϕ<sup>i</sup> ¼ β þ b<sup>i</sup> and the i-th plant belongs to the treatment W1. When x1<sup>i</sup> ¼ 1 and x2<sup>i</sup> ¼ x3i ¼ 0, ϕ<sup>i</sup> ¼ β þ γ þ b<sup>i</sup> and the i-th plant belongs to the treatment W2. If x2<sup>i</sup> ¼ 1 and x1<sup>i</sup> ¼ x3i ¼ 0, ϕ<sup>i</sup> ¼ β þ δ þ b<sup>i</sup> and the i-th plant belongs to the treatment W3. And lastly, if x3<sup>i</sup> ¼ 1 and x1<sup>i</sup> ¼ x2<sup>i</sup> ¼ 0, ϕ<sup>i</sup> ¼ β þ ζ þ b<sup>i</sup> and the i-th plant belong to the treatment W4. The random effects b<sup>i</sup> are considered independent among the plants and are normally distributed with mean zero and covariance matrix Ψ (b<sup>i</sup> � Nð Þ 0; Ψ ). In this case, the covariance

� �=ϕ4<sup>i</sup>

� � <sup>þ</sup> <sup>ε</sup>ij (2)

b1i b2i b3i b4<sup>i</sup>

(3)

(5)

¼ β þ γx1<sup>i</sup> þ δx2<sup>i</sup> þ ζx3<sup>i</sup> þ bi, (4)

parameters may vary among the units.

Water Chemistry

being

matrix is given by

with σ<sup>2</sup>

74

1, σ<sup>2</sup> 2, σ<sup>2</sup>

3, and σ<sup>2</sup>

ϕ<sup>i</sup> ¼

ϕ1i ϕ2i ϕ3i ϕ4<sup>i</sup>

Ψ ¼

σ2

σ<sup>12</sup> σ<sup>2</sup>

σ<sup>13</sup> σ<sup>23</sup> σ<sup>2</sup>

respectively, and σ12, σ13, σ14, σ23, σ24, and σ<sup>34</sup> are the covariances between them. The variances of the random effects indicate how a model parameter varies between the experimental units. Frequently, we suppose that the errors within the groups εij are independent between the observations of the same experimental unit and that

<sup>1</sup> σ<sup>12</sup> σ<sup>13</sup> σ<sup>14</sup>

σ<sup>14</sup> σ<sup>24</sup> σ<sup>34</sup> σ<sup>2</sup>

<sup>2</sup> σ<sup>23</sup> σ<sup>24</sup>

<sup>3</sup> σ<sup>34</sup>

4

<sup>4</sup> being the variances of random effects b1<sup>i</sup>, b2<sup>i</sup>, b3<sup>i</sup>, and b4<sup>i</sup>,

they are distributed the same way among the experimental units. For a given sample unit i, we can describe the error εij as the time j indexed vector, i.e.,

$$
\varepsilon\_{\mathbf{i}} = \begin{bmatrix} \varepsilon\_{1i} \\ \varepsilon\_{2i} \\ \vdots \\ \varepsilon\_{n\_i i} \end{bmatrix},
$$

where each vector ε<sup>i</sup> has ni observations throughout time and we assume that they follow a normal multivariate distribution with mean zero and covariance matrix Λi, i.e.,

$$
\mathfrak{e}\_{\mathbf{i}} \sim N(\mathbf{0}, \Lambda\_{\mathbf{i}}).\tag{6}
$$

Most of the times, we consider <sup>Λ</sup><sup>i</sup> <sup>¼</sup> <sup>σ</sup><sup>2</sup>I, being <sup>σ</sup><sup>2</sup> a constant variance for all <sup>j</sup> times.

By using nonlinear mixed effect models, we must consider the technical difficulty in the parameter estimation. In a mixed effects linear model, the derivative of the logarithmic of the likelihood function allows, in a simple way, the algorithm implementation, like Newton-Raphson, to obtain the estimative of the models parameters. Nonlinear models can, however, present nonlinear random coefficients which make it impossible to directly explain the parameters from the likelihood function. Methods that depend on linear approximations such as the first-order Taylor approximation can be used to estimate the model.

Nonlinear mixed effect model analysis can be, preferably, made by the R software [13] with the package name [14], also, at the SAS software using PROC NLMIXED. An excellent text to learn how to use these skills is given in [4].
