Forest Research: Open Access
Open Access

Determination of climate predictor variables effecting on annual cone harvest and seed yield of Korean pine (Pinus koraiensis Siebold & Zucc.) seed orchards

Author info »

Introduction

Korean pine (Pinus koraiensis Siebold & Zucc.) is one of the native pines in Korea. It usually grows in high-altitude natural forests over 1,000 meters, but also grows relatively well in plantation forests of 300 to 800 meters. The total forest area of the species is 230,282 ha, which is 5.3% of the total forest area of Korea [1]. Clonal seed orchards of the species had been established in three regions across the country, with a total area of 90 ha. Over the past 10 years (2009-2018), 11.4 tons of seeds were yielded annually from the seed orchards and used to create 1,024 ha of afforested forests per year [1].

For successful afforestation, a sufficient quantity of seeds and seedlings must be prepared in advance. However, it is not easy to collect and supply the required amount of seeds at the right time in line with the actual process of afforestation, because the seed yield is not constant across years due to environmental influences. Hence, predicting the optimal timing and environmental factors for maximizing seed yield has long been a major concern [2]. Many studies have been conducted in conifer species, such as examining the phenology and reproductive cycle [3-6], the effect of tree age, growth, density [7-9], and modeling for temporal variation of the seed yield [10-14].

Most pines are known to have a three-year reproductive cycle. The cycle is a sequential process of initiation of reproductive bud primordia (two years before the seeding year), flowering and pollination (one year before), and fertilization and seed development (the seeding year) [15]. Many researchers have examined which climate factors in the reproductive cycle have a significant effect on flowering and the consequent seed production of western white pine (Pinus monticola) [16,17], Norway spruce (Picea abies), Scots pine (Pinus sylvestris) [2,18], and white spruce (Picea glauca) [10]. The common result of these studies is that pines flower abundantly in weather conditions where the first summer (one year before the flowering year) is cool and the second (the flowering year) is warm; further, the most significant climate predictors are temperature in May (positive correlation) and summer rainfall (negative correlation) [2]. With respect to the Korean pine, only three studies have reported climate factors related to the number and weight of cones per tree [19-21]. The results of the studies differ from each other, so assumably, the significant common climate factors have not yet been identified. Further relevant studies are needed to confirm this.

The aim of this study was to develop a predictive model for the cone harvest and seed yield of Korean pine seed orchards using weather data from 12 years (2008-2019). Data included monthly mean temperature, monthly total precipitation, and monthly hours of sunshine. Subsequently, sets of climate predictor variables that significantly account for annual variation in cone harvest and seed yield were selected based on a stepwise regression analysis. Following this, the selected variables were examined whether they were associated with phenological reproductive processes in the three-year cycle of pines, so to determine whether they could be the predictive climate factor for cone and seed production of Korean pine seed orchards.

Materials and Methods

Seed orchards description

The seed orchards of Korean pine are located in three regions: Chuncheon-city (37° 52.5′ N, 127° 36.5′ E, 439 m in elevation) and Gangneung-city (36° 35′ N, 128° 49′ E, 889 m) in Gangwon Province and Chungju-City (36° 52′ N, 127° 59′ E, 444 m) in Chungbuk province. They are all clonal seed orchards, which comprise grafts of selected genotypes. The total area spans 90 ha; its details are given in Table 1.

Table 1: Details of three seed orchards of Korean pine.

Data collection and preparation

The data for cone harvest and seed yield over the past 12 years (2008-2019) were obtained from the National Forest Seed Variety Center. They are the annual total quantity (kg) of the harvested cones and purified sound seeds in each seed orchard. However, the distributions of the data highly were skewed, and had serious departures from normality in the Shapiro-Wilk test (W=0.76, p<0.0001) [22]. The raw data were transformed to natural logarithms to make them much less skewed and more symmetric [23]. The log-transformed data were used as response variables in a subsequent regression analysis.

Fourteen years (2006-2019) of meteorological data for three climate factors, namely, monthly mean temperature (MMTp), monthly total precipitation (MTPr), and monthly hours of sunshine (MHSs), were obtained from website of the Korea Meteorological Administration (https://data.kma.go.kr/cmmn/main.do). They are the recorded data from the weather stations closest to each of the three seed orchards. A total of 102 predictor variables were created and used for a subsequent regression analysis. They included 30 climate variables for 10 months (from January to October) of a given seeding year (t0), 36 for 12 months (from January to December) of one year before the given seeding year (t1), and 36 for 12 months (from January to December) of two years before the given seeding year (t2). The names of the 102 variables were written in the following manner: the year (t0, t1, t2), the climate factor (MMTp, MTPr, MHSs), and the number of months (1 to 12) were arrayed sequentially in the name. For example, the name “t0_MMTp_5” denotes the climate predictor variable for the “monthly mean temperature” of “May” in the seeding year “t0”.

A new data set was prepared and used for a subsequent regression analysis, in which the response variables (the log-transformed cone and seed production) in a given year in a given seed orchard corresponded to the 102 predictor variables (climate variables) in the same year in the same seed orchard.

Data analysis

All statistical analyses were performed using SAS JMP® Pro (https://www.jmp.com). Distributions for the time series of the cone harvest, seed yield, log-transformed cone harvest and logtransformed seed yield were analyzed. Additionally, time series plots were prepared and analyzed to determine whether represent a trend and seasonality. Pearson’s simple and partial correlation coefficients between the log-transformed data and 102 climate variables were calculated [22,24].

The following model was used for the stepwise regressing analysis:

where Yi, i=1,…,36, denotes the i-th observation of the cone harvest or seed yield, b0 is an intercept, Xij is the respective value of the j-th climate predictor variable Xj, and bj, j=1,…,102, denotes the linear coefficient. The error term Ɛi reflects factors that are not accounted for in the model.

A stepwise regression method was used to select candidate regressions of the best climate variables for the log-transformed cone harvest and seed yield, in which a threshold value of 0.05 was used to specify whether a given variable can be entered into the model or whether the variable can remain in the model [24]. The plot for the residuals using predicted values from the selected regression was examined to test whether the residuals were random. The Durbin-Watson test was performed to examine autocorrelation between consecutive residuals. The variance inflation factor (VIF) for each predictor variable was checked to assess multicollinearity of the regression [24,25].

Examination of associations between predictor variables and phenological reproductive processes

Associations between the predictor variables of the selected regression models and phenological processes during the threeyear reproductive cycle of the pines were examined as follows; a time-series table of phenological processes during the three-year reproductive cycle of western white pine was prepared by reference to the table presented in the paper of Eis [17]. Then, the predictor variables were assigned to the corresponding periods in the threeyear reproductive cycle to examine which variables were associated with the phenological processes of which periods, and also the associations were compared with those of relevant studies in western white pine [16,17] and Korean pine [19].

Results

Annual cone harvest and seed yield

The annual mean of cone harvest, seed yield, log-transformed cone harvest (LN (cone)) and seed yield (LN (seed)) between 2008-2019 was 25,803 kg, 3.999 kg, 9.7, and 8.0, respectively (Table 2). It appeared that Chungju and Gangneung seed orchards annually harvested more cone than Chuncheon, but had a lower seed yield. However, there was no significant difference among seed orchards in the annual mean of cone harvest (F=0.6387, p=0.5348), seed yield (1.3315, 0.2779), LN (cone) (0.0423, 0.9586), and LN (seed) (2.3553, 0.1106).

Table 2: Descriptive statistics for cone harvest and seed yield, and their logarithms.

The 12-year time series plot for the mean value of the cone harvest, seed yield, LN (cone), and LN (seed) are shown in Figures 1 and 2. The values did not show a distinct increasing or decreasing trend. Moreover, the results from the time-series regression analysis showed that all slope of the fitting regression line were close to zero. This indicates that the time series had no trend component. Conversely, at first glance, the time series plot seemed likely to have a seasonal cycle with a three-year interval between peaks. However, the value of the autocorrelation function of the time series was less than 0.33, indicating that there is no seasonality. Therefore, no further analysis was conducted in this study, because 12 years of time series data on seed production were found to have no significant trend or seasonal components.

Figure 1: Time series plot for mean of cone harvest, seed yield (left), LN (cone) and LN (seed) (right). Red line: Cone harvest and LN (cone); Blue line: Seed yield and LN (seed); Black dotted line: Mean value.

Figure 2: Plot for observed values of LN (cone) (left) and for residuals (right) by predicted values. Blue line denotes mean value.

Climate variable correlations

The pair-wise correlation matrix between 17 variables (2 response variables and 15 climate predictor variables) is present in Table 3. In the simple correlation, 10 climate variables were significantly related to the LN (cone) and LN (seed), in which five variables were significantly correlated with the LN (cone), and seven with the LN (seed), and two significantly correlated with both. In the partial correlation, however, five variables were significantly related to LN (cone) and LN (seed), in which two were significantly related to LN (cone), and three to LN (seed). The two variables, t2_MTPr_2 and t2_MTPr_10, were significantly related to the LN (cone) and LN (seed) in both correlations. Each of the 15 variables had the same sign of correlation coefficient to both the LN (cone) and LN (seed); that is the same direction of the relationship, even though the coefficient value was not the same. Among the 15 variables, one was MMTp, 10 were MTPr, and the remaining four were MHSs.

Table 3: Pearson’s simple and partial correlations between 17 variables.

Regression model for LN (Cone) and LN (Seed)

Two models for LN (cone) and LN (seed) were selected from the stepwise regression analysis as follows:

LN (cone)=10.0375 + (-0.0212) × (t2_MTPr_2) + (-0.0086) × (t2_MHSs_6)+ (0.0077) × (t2_MTPr_10) + (0.0089) × (t2_MHSs_10), LN (seed)=7.3969 + (-0.0278) × (t2_MTPr_2) + (-0.0009) × (t2_MTPr_7) + (-0.0059) × (t1_MTPr_5) + (0.0083) × (t1_MHSs_6) + (-0.0036) × (t1_MTPr_10) + (0.0091) × (t1_MTPr_11) + (-0.0049) × (t0_MTPr_10).

The model for LN (cone) with four climate predictor variables (t2_MTPr_2, t2_MHSs_6, t2_MTPr_10, and t2_MHSs_10) was highly significant at a level of p=0.05 for the F-ratio (F=12.2, p ≤ 0.0001), accounting for 62.6% of the total variation in the LN (cone) (R²=0.6263, RMSE=0.6165) (Table 4). The plot for the observed values by the predicted values from the model parameter estimates showed that the model fit the observed values relatively well, even though there were some deviations between the values (Figure 2). The plot for the residuals by the predicted values showed that the residuals were randomly scattered around zero with no obvious pattern. The Shapiro-Wilk W test for normality of the model was not significant at a level of p=0.05 (W=0.9756, p=0.6310). All the VIF values for four variables were much smaller than 5 (Table 5), indicating that the model for LN (cone) has no significant multicollinearity with four predictor variables. The two variables (t2_MTPr_2 and t2_MHSs_6) were negative for the LN (cone), while the other two (t2_MTPr_10, and t2_MHSs_10) were positive (Table 5). Based on the partial regression coefficient and the absolute value of the standardized coefficient, t2_MTPr_2 was the most important predictor variable (0.2211, 0.5397), and the relative importance of the remaining three variables was similar (0.0877-0.1964, 0.2979-0.3862).

Table 4: Goodness-of-fit of regression models for logarithms of cone harvest and seed yield.

Table 5: Parameter estimates of model for cone harvest (LN (cone)).

In the model for the LN (seed), seven climate predictor variables (t2_MTPr_2, t2_MTPr_7, t1_MTPr_5, t1_MHSs_6, t1_MTPr_10, t1_MTPr_11 and t0_MTPr_10) were included, and the model was highly significant at a level of p=0.05 (F=21.4, p ≤ 0.0001), accounting for 80.3% of the total variation in LN (seed) (R²=0.8425, RMSE=0.3673) (Table 4). The plot for the observed values by the predicted values from the model parameter estimates showed that the model fitted well with the observed values, even though there were a few deviations between the observed and the predicted values (Figure 3). The plot for the residual by the predicted value showed that the residuals of the model were randomly scattered around zero with no obvious pattern (Figure 3). The Shapiro-Wilk W-test for normality of the model was not significant at a level of p=0.05 (W=0.9641, p=0.2865). The four climate predictor variables (t2_MTPr_2, t1_MTPr_5, t1_MTPr_10, and t0_MTPr_10) were negatively related to LN (seed), while the other three variables (t2_MTPr_7, t1_MHSs_6 and t1_MHPr_11) were positive (Table 6). All VIF values for the seven variables were much smaller than 5, indicating that the model for LN (seed) had no significant multicollinearity with seven variables (Table 6). Based on the partial regression coefficient and the absolute value of the standardized coefficient, t2_MTPr_2 was the most important predictor variable in the model (0.3637, 0.7884), and the relative importance of the remaining six variables was similar (0.0306-0.0583, 0.1860-0.4326).

Figure 3: Plot for observed values of LN (seed) (left) and for residuals (right) by predicted values. Blue line denotes mean value.

Table 6: Parameter estimates of model for seed yield (LN (seed)).

The regression results can be summarized as follows: first, approximately 63% of the total variation in LN (cone) can be explained by the fitting model consisting of four climate predictor variables of two years before the seeding year. A higher cone harvest would be expected if there was less precipitation in February, longer sunshine in June and October, and more precipitation in October. Second, approximately 84% of the total variation in LN (seed) can be accounted for by the fitting model composed of seven climate predictor variables (two of two years before the seeding year, four of one year before the seeding year, one of the seeding year). A higher seed yield would be expected if there was less precipitation in February of the two years before, May and October of one year before, and October of the seeding year. Additionally, a higher seed yield would be expected if more precipitation occurred in July of two years before and November of one year before, and if longer sunshine occurred in June of one year before. Third, February precipitation for two years before the seeding year is the climate predictor variable that most affects both cone harvest and seed yield. It would be expected that the lower the February precipitation, the higher the cone harvest and seed yield.

Associations of predictor variables with phenological reproductive processes

Ten climate predictor variables of the regression models for the LN (cone) and LN (seed) were assigned to five phenological phases in the three-year reproductive cycle (Table 7). All four variables associated with LN (cone) were assigned only to the first phenological phase. Conversely, the seven variables related to LN (seed) were assigned to three phenological phases and one interphase; two of the seven variables were assigned to the first phenological phase, two variables to both the second phase and interphase, and one variable to the third phase. The assignment pattern of the climate variables for LN (cone) is similar to that of Rehfeldt et al. [16], but the pattern for LN (seed) is similar to that of Eis [17] and Joo et al. [19].

Table 7: Associations of climate predictor variables with phenological processes of three-year reproductive cycle.

It is noteworthy in this study that many climate variables (five out of 10 variables) were assigned to the period of the first phenological phase where the processes of the LSB development and the initiation of pollen and cone bud take place. Further, in particular, the assigned variables explained a large proportion of total variation in each of the LN (cone) and LN (seed) (62.6% and 39.7%, respectively). Rehefeldt et al. [16] and Eis [17] also showed similar results, where four of six variables and eight of 16 variables were assigned to the same phase in each study. According to Owens [15], the LSB is a different type of vegetative bud of pine species compared with that of other conifers. It contains all types of buds, such as pollen and seed cone bud, short-shoot bud, and lateral branch bud, which is found on new elongated shoots in the following year. These facts indicate that the normal formation and development of LSB two years before the seeding year may be a prerequisite for determining the amount of flowering, cone harvest and seed yield. Consequently, it is thought that the six climate variables assigned to the first phenological phase of the LSB development and initiation of pollen and seed cone bud in this study may be good predictors for cone harvest and seed yield of Korean pine.

Particularly noteworthy is the February precipitation among the five climate variables allocated to the first phenological phase, which is the most influential variable that accounted for 22.1% and 36.3% of the total variation of LN (cone) and LN (seed), respectively, and is the limiting factor that negatively affected both of them the most. February in Korea is the period between late winter and early spring, so the temperature is low. Moreover, the average temperature in February in the three seed orchards over the past 12 years was -2.12°C (± 2.51). Therefore, it can be inferred that the higher the precipitation at this time, the lower the temperature, which would inhibit or delay the LSB bursting and development, and eventually lead to a reduction in flowering and consequently cone and seed production.

Five of seven variables related to LN (seed) were assigned to the following phenological periods; namely, two variables were assigned to the period of flowering and pollination in one year before the seeding year, two to the period of seed cone dormancy in the same year, and one to the period of the cone and seed maturation in the seeding year. These three groups of variables assigned to each period can be ranked in the following order based on the proportion of total variation in LN (seed) explained by them: May precipitation and Jun sunshine hours (9.1%) in one year before the seeding year, October and November precipitation (7.0%) in the same year, and October precipitation (5.8%) in the seeding year. The assignment of the five variables to the phenological periods in this study seems to be slightly similar to that of Eis [17] rather than that of Joo et al. [19], except for October precipitation in the seeding year. In addition, this study and Eis [17] showed that the variables related to the LN (seed) and cone number were not assigned to the period of fertilization in Jun and July, the seeding year. Conversely studies by Rehfeldt et al. [16] and Joo et al. [19] showed that dry weather conditions and July temperature were correlated with the cone number and the weight of cone and seed. The comparison between the studies indicates that there were no common climate variables assigned to the phenology period of one year before the seeding year and the seeding year. Further studies are needed to ascertain whether the difference between the studies is caused by different statistical methods used in each study or by different local weather conditions in each study site, to determine the more generalized climate variables related to seed yield of Korean pine [15-26].

Discussion

Given that the value of R² and RMSE of the regression model and the association of predictor variables of the model with the phenological phases, the model for LN (seed) seems to be statistically more reliable in explaining the total variation of LN (seed) and also practically more applicable to determine the optimal timing of seed production rather than the model for LN (cone). However, the model for LN (seed) in this study has the weakness of fewer observations (36 data) compared to the number of predictor variables, resulting in the possibility of over-fitting of the model. In general, it has been suggested that there should be at least 10 observations per variable in a regression analysis [25], which means that at least 70 observations of seed yield would be needed in the model for LN (seed) with seven climate variables in this study. Therefore, a Least Absolute Shrinkage and Selection Operator analysis [26] was additionally conducted to reduce the number of predictor variables of the model, where the number of variables included in the shrunk new model was set to four appropriate variables for the 36 observations in this study [26]. The result of the additional analysis showed that the new model comprised four variables (t2_MTPr_2, t1_MTPr_5, t1_MHSs_6, t0_MTPr_10) with a R² of 0.5217 and a square root of the mean squared prediction error of 0.5643, became lower in predictive precision than that of the previous one (0.8425, 0.3239). All considered, it is difficult to determine which of the two models would be better in explaining the variation in seed yield in the seed orchards of Korean pine, and it will only be possible after additional data are obtained and analyzed.

It seems that the best way would be to obtain more observations within approximately five years and then build a new model. The number of observations should be more than ten times the number of variables included in the new model. However, within five years, it is impossible to obtain more than the minimum number of observations only from the three seed orchards of Korean pine. Instead, it is desirable to obtain data from many forests with density and age similar to that of the seed orchards in multiple regions with different weather conditions in Korea and adjacent countries, if possible. Subsequently, a new model should be re-built using recently developed modeling methods. In such case, the results of this study are expected to be used as preliminary data, namely as training or comparative dataset for future research.

Conclusion

In this study, we built regression models that could significantly explain large proportion of total variation in cone harvest and seed yield in three seed orchards of Korean pine over the past 12 years (2008-2019) using climate predictor variables. The results showed that the model for seed yield would be better than that of cone harvest in predicting the quantity of seed production and making mid-term plans for seed collection. The model for seed yield indicated that seed yield was significantly related to seven climate variables (temperature and sunshine hours) assigned to seven phenological periods for the LSB bursting, LSB development, pollen and cone bud dormancy (in two years before the seeding year), flowering, pollination, seed cone dormancy (in one year before the seeding year), and cone and seed maturation (in the seeding year). Particularly, precipitation during the period of the LSB bursting two years before the seeding year was a major limiting climate variable that had the greatest effect on seed yield. The remaining variables assigned to other phenological periods were the minor variables that significantly and independently affected the seed yield.

The model built in the study seems to be meaningful in that it clearly present which climate variables are associated with the seed production of Korean pine and also to what extent they affect it, despite the possibility of over-fitting due to relatively many variables compared to the number of observations. It is hoped that further studies on a more advanced predictive model will be carried out based on these results.

Acknowledgement

We greatly appreciate SAS JMP Korea for allowing us to use the JMP Pro software for statistical analysis of this study. We also thank the anonymous reviewers for their useful comments on this manuscript.

References

1. Pukkala T, Hokkanen T, Nikkanen T. Prediction models for the annual seed crop of Norway spruce and Scots pine in Finland. Silva Fenn. 2010;44(4):629-642.
2. Owston PW. The shoot apex in eastern white pine: its structure, seasonal development, and variation within the crown. Can J Bot. 1968;47(7):1181-1188.
3. Owens JN, Molder M. Development of long-shoot terminal buds of western white pine (Pinus monticola). Can J Bot. 1977;55:1308-1321.
4. Owens JN, Molder M. Seed-cone differentiation and sexual reproduction in western white pine (Pinus monticola). Can J Bot. 1977;55(20):2574-2590.
5. Philipson JJ. Predicting cone crop potential in conifers by assessment of developing cone buds and cones. Forestry. 1997;70(1):87-96.
6.  Ordonez JL, Retana J, Espelta JM. Effects of tree size, crown damage and tree location on post-fire survival and cone production of Pinus nigra trees. For Ecol Manag. 2005;206(1-3):109-117.
7. Ayari A, Khouja ML. Ecophysiological variables influencing Aleppo pine seed and cone production: a review. Tree Physiol. 2014;34(4):426-437.
8.  Nguyen TT, Tai DT, Zhang P, Razaq M, Shen HL. Effect of thinning intensity of tree growth and temporal variation of seed and cone production in a Pinus koraiensis plantation. J For Res. 2019;30(3):835-845.
9. Krebs CJ, LaMontagne JM, Kenney AJ, Boutin S. Climatic determinants of white spruce cone crops in the boreal forest of southwestern Yukon. Bot. 2012;90(2):113-119.
10. Bisi F, Von Hardenberg J, Bertolino S, Wauters LA, Imperio S, Preatoni DG, et al. Current and future conifer seed production in the Alps: testing weather factors as cues behind masting. Eur J For Res. 2016;135(4):743-754.
11.  Mencuccini M, Piussi P, Sulli AZ. Thirty years of seed production in a subalpine Norway spruce forests: Patterns of temporal and spatial variation. For Ecol Manag. 1995;76:109-125.
12. Mutke S, Gordo J, Gil L. Variability of Mediterranean stone pine cone production: yield loss as response to climate change. Agric For Meteorol. 2005;132:263-272.
13. Nussbaumer A, Waldner P, Apuhtin V, Aytar F. Impact of weather cues and resource dynamics on masting occurrence in the main forest tree species in Europe. For Ecol Manag. 2018;429:336-350.
14. Owens JN. The reproductive biology of western white pine. Forest Genetic Councils of British Columbia. Extension Note No. 2004: 04.
15. Rehfeldt GE, Stage AR, Bingham RT. Strobili development in western white pine: periodicity, prediction, and association with weather. For Sci. 1971;17:454-461.
16. Eis S. Association of western pine cone crops with weather variables. Can J For Res. 1975;6(1):6-12.
17. Nikkanen T, Ruotsalainen S. Variation in flowering abundance and its impact on the genetic diversity of the seed crop in a Norway spruce seed orchard. Silva Fenn. 2000;34(3):205-222.
18. Joo YT, Chon SK, Chung DJ. The effect of some meteorological factors on the seed characteristics in Korean white pine (Pinus koraiensis S. et Z.): The weight of cone and seed per cone. Korean J Agric Forest Meteorol. 1999;1:20-29.
19. Kim IH, Shin MY, Kim YC, Chon SK. Evaluation of reproductive growth in a mature stand of Korean pine under simulated climatic condition. Korean J Agric Forest Meteorol. 2001;3:185-98.
20. Shin MY, Jang YS, Han SU, Kim YC. Effects of local climatic conditions on the yearly cone production in progeny test stands of Korean white pine. Korean J Agric Forest Meteorol. 2002;4(3):141-150.
21. Rushing H, Karl A, Wisnowski J. Design and analysis of experiments by Douglas Montgomery: A Supplement for Using JMP . Cary, NC: SAS Institute Inc. 2013.
22. Feng C, Wang H, Lu N, Tu XM. Log transformation: Application and interpretation in biomedical research. Statistics in Medicine. 2002;32:230-239.
23. Field A, Miles J, Field Z. Discovering Statistics Using R. Newbury Park, CA: SAGE Publications Ltd. 2012; p: 992.
24. Troutt MD. Regression, 10k rule of thumb for. Encyclopedia of Statistical Sciences. New York: Wiley, USA. 2006.
25. Klimberg R, McCullough P. Fundamentals of Predictive Analytics with JMP. (2nd edn). Cary, NC: SAS Institute Inc., USA, 2016.