Phenological Changes of Corn and Soybeans over U.s. by Bayesian Change-point Model

In this paper, a Bayesian change-point model was used to examine the phenological changes in the predominant crop producing states of U.S over a 33-year period (1981–2013). Changes of phenological observation were categorized into a no-change model and two change models. The change point and intensity of shifts were subsequently estimated under the selected change model. The experiments were conducted in the cropping regions using the state-level crop progress reports issued by the U.S. Department of Agriculture. The results demonstrated that the planted, silking and mature stages of corn were significantly advanced under the change models; the vegetative period was shortened, and the reproductive and growing seasons were lengthened. The soybean phenological metrics followed a similar trend as that of corn, even though more states tended to change under a change model. The underlying drivers of such abrupt changes may be the confounding effects of crop breeding, agronomic management and climate change. Specific events, such as the adoption of genetically engineered crops in 1996–1997, can partly explain the changes in phenology. A comparison with the breakpoints function and Pettitt method demonstrated the feasibility and effectiveness of the Bayesian change-point model on crop phenological change detection.


Introduction
An increasing number of studies have observed interannual variations in the timing of phenological events of various crop species, e.g., corn and soybean [1], in many regions worldwide, including U.S. [2], China [3][4][5] and Europe [6], in the last several decades.Previous studies have shown that crop growth is affected by a variety of factors [7], such as breeding, agronomic management [8], and climatic conditions.Specifically, the responses of crops to climate change have been extensively investigated using satellite images [9], historical phenological archives [10] and meteorological data [11].Long-term phenological data indicate that the date of spring advances by 1-2 days per decade, on average, but varies depending on the study period [12].A trend of earlier crop planting and longer growing period has been observed in the U.S. over the past several decades [10].
When evaluating phenological trends, the most commonly used method is regression of dates of a certain phenological event over a year where change is expressed as the magnitude and direction of the slope in days per year [6].However, while allowing the calculation of the rate of change in a phenological event over time, this method is best suited for time periods with relatively gradual trends [11], and often ignores abrupt changes.Sudden temporal changes in a crop phenology time series are important, because they provide accurate and valuable information for examining the mechanisms and processes of crop responses to biotic or abiotic forces, e.g., extreme climatological events [13].
The detection of abrupt phenological changes can be regarded as a change-point problem, and a number of methods have been developed, including the maximum-likelihood estimation, isotonic regression, quasi-likelihood, and non-parametric regression to all types of Bayesian approaches.Chen et al. [14] identified three types of changes in vegetation photosynthetic activities, namely seasonal changes, gradual trend changes and abrupt changes.Breaks and turning change points were assessed using least-square linear and expanded paired-consecutive regression models.Dose and Menzel [15] and Schleip et al. [16] represented the changes in phenological events throughout a year using three types of models: constant model, linear model and one change-point model, and concluded that the one change-point model is preferred with considerable discrimination of the two other alternatives.
The Bayesian approach provides a coherent and rational framework for distilling uncertainties from the subjective historical records [17], and has the capability of making inferences on the posterior distribution with respect to the change-point [18].Another advantage of Bayesian change-point analysis is the minimal requirement for extra analytic solutions for parameter inference [19].The change detection problem is addressed by choosing an optimal model from a set of optional models, and then the parameters of the selected model are estimated.Chu and Zhao [20] utilized the Bayesian change-point analysis to model the shifts in annual tropical cyclone frequency in the central North Pacific region during 1966-2002.Perreault et al. [21] conducted Bayesian analyses of several change-point models using univariate normal data with a hydrologic time series.Kim et al. [17] used Bayesian change-point analysis to detect a single change-point of the annual maximum daily and subdaily precipitation in South Korea.
In this study, we employed the Bayesian change-point model to infer the crop phenological changes, and analyze the potential factors that induced the changes.The time points and intensities of the changes of corn and soybean phenological metrics in the U.S. were estimated for the period of 1981-2013, with an assumption that an abrupt change occurred during the study period.Our hypothesis is that a maximum one change-point exists in the crop phenology time series and the data follow a normal distribution.For a short time period (e.g., 1981-2013), the single change point assumption is plausible [17].Although it is nearly impossible to isolate the impacts of breeding, agronomic management, and climatic conditions on phenological changes, we sought a correlation between specific events and the phenological changes to device explanations of the detected change points.Finally, a comparison study was conducted to examine the feasibility and effectiveness of our proposed method on the estimation of crop phenological changes.

Study Areas and Datasets
This study was conducted for the period of 1981-2013 in the predominant producing states of corn and soybeans in the U.S. For corn, 15 states with continuous crop progress stages reports (CPRs) [22] were included: Colorado (CO), Illinois (IL), Indiana (IN), Iowa (IA), Kansas (KS), Kentucky (KY), Michigan (MI), Minnesota (MN), Missouri (MO), Nebraska (NE), North Carolina (NC), Ohio (OH), Pennsylvania (PA), South Dakota (SD), and Wisconsin (WI), which occupy approximately 83.9% of the total U.S. corn harvested area.For soybeans, 15 states were selected: Arkansas (AR), IL, IN, IA, KS, KY, Louisiana (LA), MI, MN, Mississippi (MS), MO, NE, NC, OH, and Tennessee (TN), which account for approximately 81.1% of the total U.S. soybean harvested area.In the selected states, corn planting generally occurs in April and May, and soybean planting occurs in May and June.
The phenological data (i.e., weekly assessments of cumulative state-level CPRs) are provided by the National Agricultural Statistics Service (NASS) of the U.S. Department of Agriculture (USDA).According to the definition of corn and soybean progress stages [23], two important crop progress categories are usually surveyed in the entire growing cycle: the progress of farming activities (i.e., planted and harvested stages) and phenological stages.For corn, the progress starts from the planted stage when seeds are sown into the field, and then transforms into phenological stages (i.e., emerged, silking, dough, dent, and mature stages), and ends at the harvested stage.The progress of soybeans also starts from the planted stage and ends at the harvested stage, with phenological stages (i.e., emerged, blooming, setting pods, and dropping leaves stages) in between.Furthermore, the phenological stages are divided into vegetative and reproductive growth stages.Theoretically, the vegetative growth stages begin with the emergence of corn and soybeans [24].However, the emergence stages are not continuously recorded over the analysis period in the CPRs.Therefore, in this study, the vegetative stages of corn and soybeans are assumed to start from the planted stages [1].For corn, the vegetative period starts from the planted stage and ends at the silking stage, and the reproductive period starts from the silking stage and ends at the mature stage.The entire growing season of corn includes both the vegetative period and reproductive period.For soybeans, the vegetative period (from planted to blooming), reproductive period (from blooming to dropping leaves), and entire growing season (from planted to dropping leaves) are defined similarly.We evaluated six phenological metrics that may be used as indicators of biophysical conditions.For corn, the metrics are the planted stage, silking stage, mature stage, vegetative period, reproductive period, and the entire growing season.For soybeans, the metrics are the planted stage, blooming stage, dropping leaves stage, vegetative period, reproductive period, and the entire growing season.

Table 1.
States and years where the data collection ended before the progress was 50% complete or began after the progress had already reached 50%.

Bayesian Change-Point Model
Given a sequence of observations ( = ( , , … , ) ) on independent random variables ( = ( , , … , )), a single change-point can occur at any time point between 1 and − 1, dividing the sequence of random variables into two parts.Suppose the random variables in both the first and second parts of the sequence follow the normal distribution, the sequence can be represented as follows: where and are the mean and variance of the first part; and and are the mean and variance of the second part.
For such a sequence, three models can be devised: When characterizing the time and magnitude of a single change in a stochastic time series, the Bayesian change-point analysis often consists of two steps [17]: (1) Bayesian model selection, which determines the type of change using the Bayes factor or posterior probability; and (2) change-point parameter estimation, which identifies the change-point and estimates its amplitude under the selected change model (i.e., or ).

Bayesian Model Selection
Bayesian model selection evaluates the posterior probability for each model provided with the observations [19].The posterior probability of each model is more useful than the Bayes factor for selecting a model from more than two models [27].In this study, the procedure of Bayesian model selection in the change-point analysis is performed using the posterior probability through the Bayes factor.The posterior probability ( | ) of model ( = 1,2) is determined as follows: where = 1,2 ; and are the change-points of model ; ( ) and ( , ) are the prior probabilities of model and , , respectively.The model probabilities assign equal weights to the "change" and "no change" alternatives [28], i.e., ( ) = 1/2 and ( ) = ( ) = 1/4 ., ( , ) is the Bayes factor of model , to model : where ( | ) and ( | , ) are the marginal density functions of given the models and , , respectively.Let = ( , ) , = ( , , ) and = ( , , , ) .Here, ∈ where is the parameter space for ( = 0,1,2 ).The marginal density functions could be integrated by the likelihood ( | ) and the prior distribution ( ) in the parameter space, e.g., ( | , ) = ( , | ) • ( ) •

Θ
. The likelihood ( | ) of each model can be estimated from the given independent observations directly, and the prior distribution ( ) reflects the prior knowledge about the parameters and must be assessed without using the observations.By referring to Jeffrey's prior [29], the prior distribution can be turned into a particular form of non-informative density by letting ( ) = / , ( ) = / and ( ) = /( • ) , where , and are undefined normalizing constants.The marginal density functions ( | ) and ( | , ) are calculated as follows: where (•) denotes the gamma function; ; and is the length of each part of random variables , which are divided by the change point .Therefore, = , and = − .From Equations ( 2)-( 6), we found that the Bayes factors , ( , ) and the posterior probabilities ( | ) were ill-defined, because there are arbitrary constants floating around in the equations due to ( ) being non-informative priors that are improper and are defined only as arbitrary constant multipliers.These factors are defined by arbitrary constants, i.e., , , and .Thus, , ( ) is defined only by the arbitrary multiplicative constant / [30], and the usual Bayes factor is indeterminate.
In Bayesian model selection or hypothesis testing problems, prior specifications are important.Many solutions to this problem have been proposed, including the conventional prior approach, the Bayes information criterion, intrinsic Bayes factor and fractional Bayes factor [31].In this study, a data-splitting method was employed for the intrinsic Bayes factor, which uses a minimal group of training samples to convert an improper posterior density to a proper posterior density [32].The geometric intrinsic Bayes factor (denoted by , ) [33,34] was used to avoid the dependence on the choice of the minimal group of training samples., is computed as follows: where the size of the minimal training sample is 4; is the length of sample; and ( ) is the th sample.

Change-Point Parameter Estimation
After the Bayesian model selection, the inference about the change-point , the model parameters , and intensity of shifts ∆ = − in a certain model is to be derived.To evaluate the posterior probability ( | , ) of the change-point, must be integrated out of ( , | , ) and vice versa.Using conjugate prior distributions for fixed , and assuming prior independence between and parameters , the posterior distributions of the parameters of each selected model can be derived [21].The marginal posterior distributions of is given as follows: where ( | ) and ( | ) are the discrete uniform distributions on the set 2, ⋯ , − 2 .
The marginal distributions of and of each selected model are finite mixtures of the associated conditional distributions with respect to the posterior mass function of , which can be calculated as follows: where denotes the non-standardized Student's t-distribution, which has a density defined by [35] (p.432).
A positive value of the intensity of shift (∆ ) implies early/shortening of phenological metrics, and a negative ∆ implies delay/lengthening of phenological metrics.

Results
Figure 1 shows the spatial distribution of corn phenological changes, estimated using the Bayesian change-point analysis described above.Six key phenological metrics, including the planted, silking, mature, vegetative, reproductive and growing season, are sequentially visualized in Figure 1a-f.We also fit the linear trends of the phenological metrics using regression (Table 2).A significant p-value was set to be <0.05.Trends are given in days per year (slope) over the period of 1981-2013, as well as the root mean square error (RMSE) for measuring the error of regression.Figure 2 and Table 3 illustrate the results for soybean.

Planted Stage
The corn planting dates exhibit a simultaneous change in both mean and variance (i.e., model ) in NE and NC (Figure 1a).The corresponding change points of both states occurred in 1995, and early sowing dates advanced by approximately 4.3 and 14.1 days, respectively.The change model was found in PA where a change point occurred in 1998 with an advancement of 6.7 days.The remaining 12 states did not have obvious time breaks in the period of 1981-2013, i.e., the corn planting dates in these states were under the hypothesis.Seven states advanced the soybean planting dates under the change models (model and ) (Figure 2a).These states are in the southern U.S. Five states, i.e., AR, KS, KY, LA and TN, had a single change in mean only (i.e., model ).The corresponding change points occurred in 1995, 1996, 1998, 2003 and 1998, and  No change point was found in the planting dates of corn and soybeans across the northeast region of the U.S. during 1981-2013 (Figures 1a and 2a).It is worth noticing that although MI is divided into two peninsulas by the Straits of Mackinac, and the statewide planting date trends are averaged with two different climate zones, the trends for both corn and soybeans changed gradually (i.e., model ).
Table 2 shows that seven of 15 states exhibit significant trends in the annual planting dates of corn.Almost all states, except for CO, experienced an earlier trend.It is worth emphasizing that CO has the minimum linear regression error among all 15 states, and the RMSE is 4.1 days (Table 2).The corn planted stages in CO changed little in the period of 1981-2013, and the dates ranged from 121 to 137 calendar days (Figure 3c).An abnormal trend for the soybean planted stages exits in NC, which is the only state exhibiting delayed dates among all 15 states (Table 3).The RMSE is 4.7 days, which is the smallest.
In addition, Figure 1a shows that the corn planting dates in NC changed under model with the greatest shift (∆ ) among the three states with changes, which may be a result of the incomplete historical phenology records that surveyed either with 50% starts or with 50% ends.Incomplete records exist in 1981, 1982, 1985, and 1995.The high RMSE (10.52 days, see Table 2) partly illustrates the data missing.Similarly, missing records exist in 2001,2004,2005,2006 and 2010 in MS for soybeans, with the highest RMSE (14.4 days, see Table 3).

Silking/Blooming Stage
The corn silking stages in KS, MO, NE and PA changed under model (Figure 1b).The change points of KS and NE occurred in 1984 with the ∆ values of 13.1 and 4.0 days, respectively, and the change points of MO and PA occurred in 1997 with the corresponding stages advanced by 9.1 and 8.9 days, respectively.The abrupt change in the planted stage may have directly contributed to the following silking stage.For example, the abrupt changes of both stages in PA occurred around 1997.
Figure 2b  Table 2 shows eight of 15 states exhibit clear trends in the annual corn silking stage.All four states under model are part of the abovementioned eight states.Nine of 15 states exhibit significant trends in the soybean blooming stage (Table 3).All eight states under the change models (i.e., model and ), except for IL, are part of those nine states.The results demonstrate that a change point can break in time series either with a significant linear trend or not.However, the probability of the former is significantly greater than the latter.

Mature/Dropping Leaves Stage
Early planting increases the likelihood that plant physiologically mature before the killing frosts occur in the fall [10].Figures 1c and 2c show that any state with a change model ( or ) had consistent advances for both corn and soybean.The number of states with an abrupt change for the soybean dropping leaves stage was more than that for corn.For soybean, six of 15 states had abrupt change, which was threefold higher than that for corn.This result may have resulted from soybeans often being planted after corn planting is complete [1], which means more human factors will disrupt the dates of soybean dropping leaves.
For the soybean dropping leaves stage, we found that there was a link between the results from the Bayesian change-point analysis and linear regression (Figure 2c and Table 3).First, six of 15 states have significant linear trend, and they also have abrupt changes identified by the Bayesian change-point analysis.Second, these six states, both in linear regression and Bayesian change-point analysis, exhibit an earlier trend.The link was also found in the corn mature stage (Figure 1c and Table 2).
Furthermore, for corn, only two states (i.e., KY and PA) had a single change in mean (i.e., model ).For soybeans, four of the six states had changes in both mean and variance (i.e., model ).The states under the hypothesis are AR, LA, MS and TN, and the states under the hypothesis are KY and NC.
Figure 2a-c show that the time series of the three key soybean phenological metrics, including the planted, blooming and dropping leaves stages, have a similar pattern in KY, all under the hypothesis.The change points all occurred in or around 1998.The corresponding ∆ ranged from 8.4-12.4days.The shifting of the soybean blooming and dropping leaves stages may have been caused by the disordered seeding dates.Furthermore, the change point of the corn mature stage is relatively close to that of soybeans, which occurred in 1999.The corn physiological maturity dates in advance were also similar to soybeans (see Figure 2c).The ∆ values of corn and soybeans are 12.6 and 12.4 days, respectively, with only 0.2 days of time difference, which may have been caused by the high correlation between the sowing dates of corn and soybeans.The effects of soybean planting dates, plus the interannual weather fluctuations, may have increased the variances of time series before and after the change point.The change patterns transformed from model to model in AR and TN.Also, as mentioned in Section 4.1, the change model ( or ) of the soybean planted stage is more likely to be found in southeastern U.S., and the dropping leaves stages also had a similar change pattern, suggesting that the maturity is largely associated with the sowing time.

Duration of Vegetative Period
For corn, the abrupt changes in the vegetative period mainly occurred in CO, KS and PA (Figure 1d).These three states had changes only in mean (model ) with a shortened duration of the vegetative period.Table 2 confirms that only these three states have significant shorten trend.Figure 1d shows that the change points for CO, KS and PA occurred in 1998, 2006 and 2010, respectively.The average duration was shortened by approximately seven days.
The trends of soybean vegetative period were much more likely to change than that of corn.As shown in Figure 2d, AR, IL, and NC experienced a shortening of the vegetative period.The abrupt change points occurred in or around 2001, and the average duration was shortened by approximately six days.The remaining four states, including LA, MS, NE and OH, had lengthened vegetative periods.However, the increasing trends are not credible because they may have been induced by the incomplete (CPRs data).In LA and MS, the records of the blooming stages are not complete, especially in the years of 2001, 2006, 2009-2011 and 2013.Non-zero values of the blooming stage were nearly always reported for the week before and after 50% progress was reached.Moreover, these two states have larger RMSE values (see Table 3), which enlarge the gaps of variance before and after the change points, and more likely change with the model .Among the other changed states, LA and MS have the highest absolute values of ∆ , −21.6 and −13.7, respectively.The effects also spread to the soybean reproductive stages (Figure 2e).The reproductive stages in LA and MS simultaneously changed in mean and variance (model ), and the absolute values of ∆ for LA and MS (16.1 and 11.3 days, respectively) are the largest among all the changed states.

Duration of Reproductive Period
In contrast to corn vegetative periods, the reproductive periods of corn for all states in the change model were lengthened (Figure 1e).The changed states were IN, MN, NE, OH, SD and WI concentrating in the north-central U.S., all under the model.The average ∆ is 7.2 days, which is similar to the corn vegetative period with only 0.2 days of time difference.A narrow gap resulted in a few states with an abrupt change in the entire corn-growing season (Figure 1f).
Similarly, most of the states exhibit an increasing trend in the amount of days in the soybean reproductive period.As shown in Figure 2e, seven of the 10 states had lengthened reproductive periods, including AR, IL, IN, IA, KS, MI and OH.The average lengthening is by approximately 8.1 days.The remaining three states, including LA, MS, and NC, have shortening trends.Similar to the soybean vegetative stage, the decreasing trends are not credible because of the incomplete data (see Table 1).In LA and MS, the shortening trends are mainly caused by the missing data in the blooming stage, and the shortening trends in NC are mainly caused by the missing data in the dropping leaves stage.

Duration of Growing Season
In most states, the corn growing season change exhibits gradual trends, except for SD.In SD, the growing season change was under the model, and the length was increased by 11 days (Figure 1f).Moreover, both the growing season and reproductive period in SD have the same change point, 1991.The prolonged trend also can be confirmed by referring Table 2, in which the slopes for all five states with significant trends are positive.
For the soybean growing season, most states changed under the same model as the reproductive period, except for AR, IL, IN, OH and NE (Figure 2f).The growing season (from the planted stage to the dropping leaves stage) in AR, IL, IN and OH transformed from a change model to the no-change model.This is mainly due to soybeans in these areas at the planted and dropping leaves stages have the same or similar models.For example, soybeans in AR at these two stages are under the change models, and under the no-change model in IL, IN and OH.In NE, the growing season is under the change model , but the reproductive period has no change.This may be due to the discrepancy of time series in magnitude between the planted and blooming stages.The trends of the soybean reproductive stage in LA, MS and NC also spread into the growing season.Out of the seven changed states, only three states had shortened growing seasons, and the remaining four states had lengthened growing seasons.The average extension is approximately 8.2 days, which is similar to that of the reproductive period.

Drivers of Crop Phenological Changes
The shifts in crop phenological metrics resulted from the cumulative changes of multiple factors, including breeding and agronomic management practices, climate change, as well as any changes due to social and economic factors.The goal for crop breeding is to obtain higher yield potential by improving the vigor and tolerance for germination in suboptimal growing conditions (e.g., cold and wet soils) as well as by increasing the resistance to virus and insect [10,36,37].The development of corn hybrids and soybean varieties have contributed to earlier planting dates and allowed for longer growing seasons resulting in crops needing to survive in cool temperate regions [37].
In 1996, genetically engineered (GE) crops were widely adopted by U.S. farmers [38].Corn and soybean contributed 75% of the global growth in transgenics between 1996 and 1997 [39].The sudden renewal of the cultivars may have increased the likelihood of a change model ( or ) for crop phenology.Figures 1a and 2a show the planting dates for corn or soybean fluctuated near 1995, especially corn in NE and NC, and soybeans in AR.The change points in most other states occurred relatively late (e.g., 1996 and 1998), possibly related to the acceptation of GE in the region, or the combination of two or three drivers.
Better agronomic management practices have positively contributed to crop growth and have also affected crop phenology.Seed treatment, e.g., seeds coated with temperature-activated polymers [40], has helped to fend early season soil-borne insects and pathogens.Herbicides have been used to suppress the growth of weeds.Moreover, modern equipment (e.g., planter) has been improved to distribute seeds more uniformly and to function better under less ideal planting conditions [10].Bruns and Abbas [36] claimed that breeding and agronomic management practices have contributed to corn planting being completed earlier today than it was 30 years ago throughout most of the U.S.
Sacks and Kucharik [1] suggested that "AR, KS and LA have seen particularly large changes in soybean planting dates", which agrees with our results that AR, KS and LA changed under the hypothesis.The mean level differences before and after the change points were approximately 13.2-16.1 days.The gaps of the two parts may have been induced by management factors such as the adoption of the early soybean production system (ESPS) in the mid-southern U.S. The ESPS, which uses early planting of the earlier-maturing varieties, has been promoted as a cropping practice that reduces the risk of high temperatures and drought stress to the crop.The ESPS has obviously influenced the soybean maturity in most states under the change models in the southern U.S. Heatherly [41] listed the soybean planted dates and the corresponding acres in MS from 1971 through 1997, which illustrated the abrupt early occurred in 1994.The result mostly coincides with our results that the change points of soybean planted and dropping leaves stages occurred in 1993 and 1994.
Phenology and climate change interact with each other [13].Changes in phenology have long been regarded as sensitive indicators of climatic change [11].Climate change also may disrupt the regularity of crop phenology.Climatic elements, such as temperature, photoperiod, sunshine hours, solar radiation and precipitation, markedly influence crop growth and phenological responses [42].Cultivars respond differently to climatic elements.Generally, temperature, photoperiod, and water availability (e.g., precipitation and irrigation) are important factors in phenology prediction models [13].Crop germination requires a certain temperature, and most corn and soybeans across the U.S. will not germinate if the soil temperature is below 50 °F.Most of the process-oriented crop growth models use the thermal index to define different phenological stages and to quantify the developmental rates [43].For example, the world food studies (WOFOST) model defines the developmental rate as a crop/cultivar specific function of ambient temperature, which may be modified by photoperiod [44].
Changes in spring have been the most commonly reported changes, with the emphasis on the advancement linked to an increase in temperature [45].Spring warming has contributed to the earlier corn planting in KS, MO and NE [1], and this effect may have been carried directly into the following silking stage (Figure 1b).Evidence has been reported that the increasing temperature during the corn vegetative period leads to a decrease in the length of this stage in growth [4].Three changed states, including PA, CO and KS, exhibit shorten trends in the corn vegetative period (Figure 1d).Increased temperature led to lengthen growing seasons for corn (e.g., in SD), and also for soybeans (e.g., in MI, IA, NE, KS, MS and LA).
Wang et al. [46] found that the trend of spring temperature in North America was not consistent from 1982 to 2006.Annual weather fluctuations affect the phenology trends resulting in year-to-year variation [47].North American agriculture has been exposed to many severe weather events during the past decades [48], e.g., the Midwestern states drought in 1983.The corn silking stage in NE and KS, where the change points occurred in 1984, would be partially a consequence of this event (Figure 1b).The soybean blooming stage in IN (changed in 2011) was greatly affected by the great drought in 2012, in part because the drought conditions developed earlier in IN than in other major Corn Belt states [49].The states changed in both mean and variance (model ), which may have been due to the interannual weather variability.To deal with the climate change, adequate cropping systems and crop cultivar have been considered as important adaptation strategies by farmers [50].
Crop phenology is also related to marketing and economics.Early planting of crops allows full-season hybrids to utilize the entire growing season, thereby increasing the profit by reducing the drying costs, i.e., rising yields.Kucharik [2] suggested that earlier planting explains a significant proportion of the corn yield trends in the western and northern Corn Belt.On the contrary, late planting in a year generally results in lower yields than timely planting [36].Although early planting of corn is recommended [51], long-term trends of early planting are unrealistic to continue indefinitely as it is limited by frozen soils [10].

Comparison of Different Methods
Aside from the Bayesian change-point analysis, Pettitt's test [52] and frequentist inference have been widely used for detecting change points in time series.The Pettitt method is a rank-based and distribution-free test for detecting a significant change in time series.The null hypothesis (no-change) against the alternative (an upward or downward change point) in a given year is tested at the 5% significance level.The frequentist inference is the main alternative approach to Bayesian inference.The breakpoints function [53], a frequentist inference based method, is available in the R package strucchange [54].The method first determines the number of change points; and if the number is greater than or equal to one, then the F statistics is used to detect the optimal change point (only one).Table 4 lists the change points (if existing) detected by these three methods, covering six key corn phenological metrics in the 15 states mentioned above.The comparison results for soybeans are presented in Table 5. Figure 3 shows special cases using these three methods.1995 1993 1995 1997+ 1997 1997 1997 1997 1997 2001+ 1997/2001 2001 1990+ 1991 1990    As shown in Tables 4 and 5, despite most of the results of the breakpoints function (BP), Pettitt method (PE) and Bayesian change point model (BC) are highly similar, there still are some differences.Figure 3A shows the time series of the corn planted stage in MN and the corresponding change points detected by these three methods.The change points detected by BP and PE occurred in 1996, while BC estimated no change point.In the Bayesian model selection step of BC method, the posterior probabilities of (i.e., the no-change model), (i.e., the mean change model) and (i.e., the mean and variance change) are 0.51, 0.42 and 0.06, respectively.( | ) (no-change model) is greater than the sum of ( | ) and ( | ) (change models) with a minimum value, which demonstrates that the Bayesian approach to model selection allows uncertainty in both model-specific parameters and the models themselves.Such cases can also be observed in IA, KY, MN, SD and WI in the corn planted stage, CO, MI and MN in the corn silking stage, CO, MN and WI in the corn growth season, IA in the soybean planted stage, MO in the soybean vegetative period, and KY in the soybean reproductive period.In Pettitt's test, the change point is determined by the sequence arrangement of elements in the time series, not the value itself.Therefore, the Pettitt method is not only less sensitive to outliers, but also insensitive to normal fluctuations in values.Figure 3B shows the normal fluctuations of time series of the corn vegetative stage in CO, Both BP and BC detected change points, while PE found none.The change points detected by BP and BC occurred in 1998.The difference among the posterior probabilities of Bayesian model selection is obvious.The posterior probability of the no-change model ( ( | )) is only 0.19.Furthermore, the Pettitt method is more sensitive to breaks in the middle of the time series [55].Thus, the Pettitt method shows excessive detection on the change point of crop phenology time series.Figure 3D shows the time series of the soybean blooming stage in NC, in which all three methods detected change points.The change point detected by BP occurred in 1997 (1997 is in the middle of the time period of 1981-2013), while the ones detected by BP and BC occurred in 2003.The posterior probability of the no-change model tends to be zero, meaning the time sequence certainly changed under a change model.However, due to the differences among the internal mechanism of these three methods, the change points are most likely different.
Figure 3C shows the time series of the corn mature stage in OH.The change point detected by BP occurred in 1987, whereas PE and BC detected no abrupt change.The posterior probability of the no-change model is 0.78.Such cases also occurred in the planted stage in KS, the reproductive stage in MI and the growth season in IN and OH for corn, and in the blooming stage in NE, the dropping leaves stage in NE and the growth season in AR for soybeans.
Therefore, these three methods have their own advantages and disadvantages.Also, errors caused by subjective surveying, missing data (see Table 1), and data interpolation, are inevitable.In the absence of prior knowledge of data, we suggest a comprehensive analysis with the integration of Bayesian change-point model and other methods (e.g., BP and PE) to improve the effectiveness of phenology change detection.
Figure 3E shows the time series of the soybean reproductive period in KS.More than one abrupt break were detected by BP method.If only one change point was allowed in the analysis, the optimal change point would have occurred in 2000, which is consistent with the results of PE and BC.Such cases were also found in other states and key phenological metrics, e.g., the corn planted stage in PA, the corn reproductive period in IN, and the soybean blooming stage in MI.It also suggests that multiple change points are possible, although this research mainly considers a single change point.

Conclusions
In this study, the Bayesian change-point model was employed to analyze the time series of key crop phenological metrics in the U.S. We found that the corn planted dates, silking stages, and mature stages were significantly advanced under the change models.The corn vegetative period in most changed states was shortened, while the reproductive period and growing season were lengthened.The soybean phenological metrics followed the same trends as those of corn.The sudden renewal of cultivars during 1996-1997 partly increased the likelihood of abrupt changes in the time series of crop planting dates.The phenological metrics of soybean are more likely to change than those of corn.This may be because soybeans are planted after corn, thus allowing their phenology to be impacted by corn.The states with changed soybean planted stages are mainly located in the southern U.S., due to the establishment of the ESPS in this area.This influence continued into the maturity stage.The comparisons between the Bayesian change-point analysis, breakpoints function and Pettitt method demonstrated the feasibility and effectiveness of the Bayesian change-point model on crop phenological change detection.

( 1 )
a no-change model ( ) where = and = ; (2) a mean change model ( ) where ≠ and = ; and (3) a mean and variance change model ( ) where ≠ and ≠ .Obviously, the first model assumes no change occurs, and the latter two models assumes otherwise.
the early sowing dates advanced by 14.3, 13.2, 11.6, 16.1 and 11.1 days, respectively.The soybean planting date in NE changed similarly as corn, under the hypothesis.The change point occurred in 1996 with an advancement of 8.4 days.Springtime warming may have greatly contributed to the changes on both the corn and soybean planting dates in NE[10].In addition, the planting date in MS also exhibited a change in both mean and variance.The date has advanced by approximately 21.6 days since 1993.Overall, the change points of corn occurred around 1996, and the change points of soybean occurred near 1997, except for LA in 2003 and MS in 1993.
shows that the soybean blooming stages in IL exhibit a change in both mean and variance (model ), with the change point in 2009 and the ∆ value of 2.9 days.The same change model was also found in IN, where the change point occurred in 2011 and was advanced by 3.8 days.The other six states, namely AR, KS, KY, MI, NC and TN, changed under model with the corresponding change points in 1997, 1995, 1998, 1997, 2003 and 2000, respectively, and were advanced by 21.1, 10.3, 8.4, 7.4, 8, 14.7 days, respectively.The effect of abrupt change from the soybean planted stage has been extended to the blooming stage, which is obvious in KY, KS, AR and TN.

Figure 3 .
Figure 3. Change points of crop phenology detected by different methods.linear = linear regression; BP = breakpoints function; PE = Pettitt method; BC = Bayesian change-point analysis.(A) Changes points detected by BP and PE; (B) Change points detected by BP and BC; (C) Change points detected by BP; (D) Change points detected by all three methods; (E) Change points by all three methods, while more than one change point for BP.

Table 2 .
Slope (s) and RMSE (r) of linear regression for corn in the predominant producing states.The abbreviation of states is defined in Section 2.

Table 3 .
Slope (s) and RMSE (r) of linear regression for soybean in the predominant producing states.The abbreviation of states is defined in Section 2.

Table 4 .
Comparison of the breakpoints function (BP), Pettitt method (PE) and Bayesian change-point analysis (BC) for detecting corn phenological changes in the predominant producing states.
-indicates there is no change point; + indicates there is more than one change points, and the change point shows in the table is the optimal change point.

Table 5 .
Comparison of the breakpoints function (BP), Pettitt method (PE) and Bayesian change-point analysis (BC) for detecting soybean phenological changes in the predominant producing states.
-indicates there is no change point; + indicates there is more than one change points, and the change point shows in the table is the optimal change point.