Mixture Transition Distribution Modelling of Multivariate Time Series of Discrete State Processes: With an Application to Modelling Flowering Synchronisation with Respect to Climate Dynamics

Abstract

A new approach to assess synchronicity developed in this chapter is a novel bivariate extension of the generalised mixture transition distribution (MTDg) model (we coin this B-MTD). The aim of this chapter is to test MTDg an extended MTD with interactions model and its bivariate extension of MTD (B-MTD) to investigate synchrony of flowering of four Eucalypts species—E. leucoxylon, E. microcarpa, E. polyanthemos and E. tricarpa over a 31 year period. The mixture transition distribution (MTDg) is a method to estimate transition probabilities of high order Markov chains. Our B-MTD approach allows us the derive rules of thumb for synchrony and asynchrony between pairs of species, e.g. flowering of the four species. The latter B-MTD rules are based on transition probabilities between all possible on and off flowering states from previous to current time. We also apply MTDg modelling using lagged flowering states and climate covariates as predictors to model current flowering status (on/off) to assess synchronisation using residuals from the resultant models via our adaptation of Moran’s classic synchrony statistic. We compare these MTDg (with covariates)-based synchrony measures with our B-MTD results in addition to those from extended Kalman filter (EKF)-based residuals.

1. Introduction

Separation or lack of overlap of flowering time in eucalypts has been suggested as a mechanism for maintaining overall ‘generic identity’ of a plant species. If, for example flowering times and pollinators overlap in sympatric species, hybridization can occur between closely related eucalypts species. Therefore examination of long-term synchrony establishes a baseline of flowering behaviour which may assist in detecting recent or future changes. Although Eucalyptus as a genus dominates much of the Australian landscape [1, 2], few studies have quantified eucalypt flowering overlap within or between species, due to the shortage of phenological data in Australia [3, 4, 7]. This chapter examines flowering synchrony over a 31 year period, 1940–1970, at the population level among four eucalypts species—Eucalyptus leucoxylon, E. microcarpa, E. polyanthemos and E. tricarpa [3, 4, 5, 6, 7, 8].

A new approach to assess synchronicity developed in this chapter is a novel bivariate extension of the of the MTDg model [6, 9] (we coin this B-MTD). The aim of this chapter is to test mixture transition distribution (MTD) and an extended MTDg with interactions; and a novel bivariate extension of MTD (B-MTD) to investigate synchrony in phenological data. The MTDg model [6, 9] was the first approach developed to study the multivariate relationship between the probability of flowering with two states of rain and mean temperature via a mixture transition distribution (MTD), assuming, however a different transition matrix from each lag to the present time (our MTDg analysis), thus generalising the MTD approach in [13], (see also [10]) which led to the development of the MARCH software to perform MTD without covariates [11, 12]. The MTDg model is different to MARCH not only in terms of incorporating interactions between the covariates but also in its minimization process, namely using the AD Model BuilderTM [14], which uses auto-differentiation as a minimisation tool. This is shown to be computationally less intensive than MARCH. The assumption Berchtold’s MTD model, namely the assumed equality of the transition matrices among different lags, was a strong assumption, so the idea of the mixture transition distribution model was to consider independently the effect of each lag to the present instead of considering the effect of the combination of lags as in pure Markov chain processes. Specifically, an extended model for MTDg analysis which accommodates interactions was developed in [6], and applied to MTDg modelling of the flowering of four eucalyptus species studied in this chapter, as multivariate time series.

This work extends both MARCH and the work in [15, 16] to allow for differing transition matrices among the lags, i.e. our B-MTD method builds on this approach of the MTDg with interaction model [6, 9]. The MTDg model with interactions showed that the flowering of E. leucoxylon and E. tricarpa behave similarly with temperature (both flower at low temperature) and both have a positive relationship with flowering intensity 11 months ago. The flowering of E. microcarpa behaves differently in that E. microcarpa flowers at high temperature.

Our B-MTD approach developed in this chapter allows us the derive rules of thumb for synchrony and asynchrony between pairs of species. The latter B-MTD rules are based on transition probabilities between all possible on and off flowering states from previous to current time. Synchronisation is also tested using residuals from the resultant models via an adaptation of Moran’s [17, 18] classical synchrony statistic, incorporating MTDg residuals [17, 18, 19].

We also apply the earlier MTDg modelling in [6] using climate covariates and lagged flowering states as predictors to model flowering states (on/off) and thus assess synchronisation using an adaptation of the approach of Moran to the resultant MTDg model and fitted residuals. We compare these MTDg (with covariates)-based synchrony measures with our B-MTD results in addition to those using the extended Kalman filter (EKF) [15, 19], based residuals obtained earlier in [21].

2. The mixture transition distribution (MTDg) and B-MTD models: mathematical formulations in brief

2.1 The MTDg model with interactions between the covariates

The MTD model with covariates was discussed in [6] and developed in [15, 19] to incorporate interactions between the covariates (e.g. rainfall, temperature variants in the case study discussed). The high-order MTD transition probabilities are computed as follows:

$$\begin{array}{l}P\left(X_t=i_0|X_{t-1}=i_1,\dots X_{t-f}=i_f,C_1=c_1,\dots C_e=c_e,M_1=m_1,\dots M_l=m_l\right)\\ =\sum _{g=1}^f{\lambda }_{{\mathit{gqi}}_g{i}_0}+\sum _{h=1}^e{\lambda }_{f+h}{d}_{{\mathit{hj}}_h{i}_0}+\sum _{u=1}^l{\lambda }_{f+e+u}{s}_{{\mathit{uv}}_u{i}_0}\end{array}$$ E1

where ${\lambda }_{f+e+u}$ is the weight for the interaction term, $q_{i_g{i}_0}$ is the transition probability from modality $i_g$ observed at time t−g and modality $i_0$ observed at time t in the transition matrix Q, $s_{{\mathit{uv}}_u{i}_0}$ is transition probability between covariate $h_1$ and covariate $h_2$ interaction term ( $v_u=d_{h_1{j}_{h_1}}\times d_{h_2{j}_{h_2}}$ ) and $X_t$ , and where ${\sum }_{g=1}^{f+e+l}{\lambda }_g=1$ and where ${\lambda }_g\ge 0$ .

We refer the reader to [6], and further works in the seminal book by Hudson and Keatley [7] for further mathematical details.

2.2 The bivariate mixture transition distribution (B-MTD)

Let $\left\{X_t\right\}$ and $\left\{Y_t\right\}$ be sequences of random variables (say two (flowering intensity) time series) taking values in the finite set N = {1, …, k}. In a f ^th-order Markov chain, the probability that $\left\{X_t,Y_t\right\}=\left\{i_0,i_0^{\prime }\right\}$ , ( $i_0,i_0^{\prime }\in N$ ) depends on the combination of values taken by $X_{t-f},\dots ,X_{t-1},Y_{t-f},\dots ,Y_{t-1}$ . In the MTD model, the contributions of the different lags are combined additively. Then a bivariate MTD model, which we denote by B-MTD, has the following formulation:

$$\begin{array}{l}P\left(\left\{X_t,Y_t\right\}=\left\{i_0,i_0^{\prime }\right\}|X_{t-1}=i_1,\dots ,X_{t-f}=i_f,Y_{t-1}=i_1^{\prime },\dots ,Y_{t-f}=i_f^{\prime }\right)\\ =\sum _{g=1}^f{\lambda }_g{q}_{i_g,i_g^{\prime },i_0,i_0^{\prime }}\end{array}$$ E2

where $i_f,\dots ,i_0,i_f^{\prime },\dots ,i_0^{\prime }\in N$ , the probabilities $q_{i_g,i_g^{\prime },i_0,i_0^{\prime }}$ are elements of an $m^2\times m^2$ transition matrix $Q=q_{i_g,i_g^{\prime },i_0,i_0^{\prime }}$ , each row of which is a probability distribution (i.e. each row sums to 1 and the elements are nonnegative) and $\lambda ={\left({\lambda }_f,\dots ,{\lambda }_1\right)}^{\prime }$ is a vector of lag parameters. To ensure that the results of the model are probabilities, that is, $0\le {\sum }_{g=1}^f{\lambda }_g{q}_{i_g{i}_g^{\prime }{i}_0{i}_0^{\prime }}\le 1$ the vector $\lambda$ is subject to the constraints ${\sum }_{g=1}^f{\lambda }_g=1$ and ${\lambda }_g\ge 0$ .

Covariates and interaction terms can be added to the bivariate MTD (B-MTD) as follows:

$$\begin{array}{l}P(\left\{X_{1,t},\dots ,X_{n,t}\right\}=\left\{i_{1,0},\dots ,i_{n,0}\right\}\mid X_{t-1}=i_1,\dots ,X_{t-f}=i_f,Y_{t-1}=i_1^{\prime },\dots ,Y_{t-f}=i_f^{\prime },\\ C_1=c_1,\dots ,C_e=c_e,M_1=m_1,\dots ,M_l=m_l)\\ =\sum _{g=1}^f{\lambda }_g{q}_{i_{1,g},\dots ,i_{n,g},i_{1,0},\dots ,i_{n,0}}+\sum _{h=1}^e{\lambda }_{f+h}{d}_{{\mathit{hj}}_h,i_{1,0},\dots ,i_{n,0}}+\sum _{u=1}^l{\lambda }_{f+e+u}{s}_{{\mathit{uv}}_u{i}_{1,0},\dots ,i_{n,0}}\end{array}$$ E3

where ${\lambda }_{f+e+u}$ is the weight for the interaction term, $s_{{\mathit{uv}}_u{i}_{1,0},\dots ,i_{n,0}}$ is the transition probability between covariate $h_1$ and covariate $h_2$ interaction term ( $v_u=d_{h_1{j}_{h_1}}\times d_{h_2{j}_{h_2}}$ ) and $\left(X_t,Y_t\right)$ , and where ${\sum }_{g=1}^{f+e+l}{\lambda }_g=1$ . For example, if both $X_t$ and $Y_t$ are time series that constitute random realizations of two states {0, 1} and the covariates $C_1,\dots ,C_e$ are also defined by bivariate states {0, 1}, then the set of all possible states for $\left(X_t,Y_t\right)$ is {(0, 0), (0, 1), (1, 0), (1, 1)}. Hence the transition matrix $Q=\left[q_{i_g{i}_g^{\prime }{i}_0{i}_0^{\prime }}\right]$ is a $4\times 4$ matrix as specified below.

Previous state $\left(X_{t-1},Y_{t-1}\right)$	Current state (Xt , Yt )
	0, 0 (1)	0, 1 (2)	1, 0 (3)	1, 1 (4)
0, 0 (1)	(1, 1)	(1, 2)	(1, 3)	(1, 4)
0, 1 (2)	(2, 1)	(2, 2)	(2, 3)	(2, 4)
1, 0 (3)	(3, 1)	(3, 2)	(3, 3)	(3, 4)
1, 1 (4)	(4, 1)	(4, 2)	(4, 3)	(4, 4)

The transition matrices $D_h=\left[d_{{\mathit{hj}}_h{i}_0{i}_0^{\prime }}\right]$ , h = 1,…, e, are $2\times 4$ matrices as below.

	Covariate state
$\left(X_t,Y_t\right)$	0 (1)	1 (2)
0, 0 (1)	(1, 1)	(1, 2)
0, 1 (2)	(2, 1)	(2, 2)
1, 0 (3)	(3, 1)	(3, 2)
1, 1 (4)	(4, 1)	(4, 2)

2.3 Synchrony analysis using Moran’s approach

Moran in [17, 18] suggested that if two series x_t and y_t are synchronous, and if x_t can be estimated by a model f(x), the residuals from series x_t fitted to f(x), and the residuals from series y_t , fitted with the same model, but with observations, y_t , then, f(y) will be positively correlated. The synchrony of two series can then be examined by testing the significance of the correlation of these two series of residuals (using the same model). Moran used an autoregressive integrated moving average (ARIMA) model to test synchrony. Moran’s theorem suggests that if two (or more) populations sharing a common linear density-dependence (in a so-called renewal process) are disturbed with correlated noise, they will become synchronised with a correlation matching the noise correlation (see details in [4], and also [6, 15, 21]).

In this chapter we adopt the kth order linear stochastic difference to assess synchrony. Goodness of fit of the second order AR (k = 2) model is obtained. The series of residuals can then be found by subtracting the predicted (fitted species) value from the observed series. In summary, synchrony (or otherwise) of two series can be established by performing a test of significance on the correlation coefficient calculated from the two series of residuals as follows:

• Calculate the residuals for say, E. leucoxylon (Leu) using its AR (2) model. We denote this residual series by R1.

• Calculate residuals for say, E. tricarpa (Tri) using this same model. We denote this residual series by R2.

• Calculate the Pearson correlation coefficient between the residual series R1 and R2 and test for its significance at p < 0.05.

Further details on how Moran’s method is used and adapted in the case of the MTDg-based models are given in [15] (see also Section 4.8). We use the functionals and parameterisations from the mixture transition distribution (MTD) analysis as the basis of our EKF modelling approach. EKF is likewise a method to estimate the past, present and future status of non-linear time series data by minimising the mean square error. We will also test whether EKF better detects asynchronous species pairs, given EKF estimates the Kalman gain and covariance matrix at each time point [15, 19].

3. Data

Flowering data were sourced from the Box-Ironbark Forest near Maryborough, Victoria, in particular the flowering records of E. leucoxylon, E. microcarpa, E. polyanthemos and E. tricarpa (1940 and 1971). Flowering intensity was calculated by using a rank score (from 0 to 5) based on the quantity and distribution of flowering [4, 20, 23].

Flowering intensity scores were dichotomised into two discrete states, namely on and off (1/0) flowering (Figure 1) as in [6]. One temperature variant, mean monthly diurnal temperature (MeanT), in addition to the monthly rainfall (Rain) were included as climate covariates in the MTDg models; along with the temperature by rain interaction effect. We used discrete state low/high (lower than median temperature vs higher than median temperature) for the temperature variable dichotomies and less/more (less than the median rainfall vs more than the median rainfall) for the rainfall variable. The cut-points for the states or low/high categories of each climate covariate are shown in Table 1.

Figure 1.

Flowering of the four eucalypts species.

Climate variables	Low (less)	High (more)
Mean diurnal temp (°C)	≤13.84	>13.84
Rain (mm)	≤40.45	>40.45

Table 1.

Cut-points for climate variables based on medians.

4. Results

4.1 Bivariate MTD (B-MTD) discrete states results

Four eucalypts species, E. leucoxylon, E. microcarpa, E. polyanthemos and E. tricarpa were modelled using the order 1 B-MTD model discussed in Section 2.2—without the inclusion of covariates (such as temperature (variants) and rainfall). These species were paired as follows: E. leucoxylon and E. microcarpa (LeuMic); E. leucoxylon and E. polyanthemos (LeuPol); E. leucoxylon and E. tricarpa (LeuTri) and so on; hence 6 pairs were modelled via B-MTD (see Table 2) for the corresponding bivariate transition probabilities (see also Figure 2).

Species	Previous state	Current state
		(0, 0)	(0, 1)	(1, 0)	(1, 1)
LeuMic	(0, 0)	0.6667	0.2280	0.1053	0.0000
	(0, 1)	0.0000	0.6000	0.1333	0.2667
	(1, 0)	0.0845	0.0376	0.8357	0.0423
	(1, 1)	0.0000	0.0612	0.4490	0.4898
LeuPol	(0, 0)	0.6970	0.0303	0.2626	0.0101
	(0, 1)	0.4444	0.3889	0.0000	0.1667
	(1, 0)	0.0562	0.0000	0.7921	0.1517
	(1, 1)	0.1309	0.0952	0.1429	0.6310
LeuTri	(0, 0)	0.6947	0.1263	0.1053	0.0737
	(0, 1)	0.0455	0.3636	0.0000	0.5909
	(1, 0)	0.2203	0.0085	0.7034	0.0678
	(1, 1)	0.0069	0.0069	0.1736	0.8125
MicPol	(0, 0)	0.7637	0.1429	0.0879	0.0055
	(0, 1)	0.1818	0.6705	0.1023	0.0455
	(1, 0)	0.2737	0.0000	0.6842	0.0421
	(1, 1)	0.0714	0.2141	0.3572	0.3573
MicTri	(0, 0)	0.7975	0.0316	0.1329	0.0380
	(0, 1)	0.2232	0.7500	0.0179	0.0089
	(1, 0)	0.1090	0.0182	0.5819	0.2909
	(1, 1)	0.0000	0.4259	0.0000	0.5741
PolTri	(0, 0)	0.7464	0.1739	0.0797	0.0000
	(0, 1)	0.0719	0.7842	0.0360	0.1079
	(1, 0)	0.3067	0.0400	0.6400	0.0133
	(1, 1)	0.0370	0.1482	0.4074	0.4074

Table 2.

Transition matrices for the 6 B-MTD models.

Figure 2.

Subcomponents of possible transitions.

The possible states for any pair of species is the set {(0, 0), (0, 1), (1, 0), (1, 1)}, where no flowering is represented as 0 (state = 0 = no flowering) and flowering is represented by a 1 (state = 1 = flowering). Since lag order 1 B-MTD models were used, the mixing probability $\lambda$ is equal to 1.0.

The corresponding transition matrices for the 6 B-MTD models are given in Table 2. These transition profiles are also shown schematically as flow diagrams in Figures 3–4, and also as transition signatures in Figures 5–6. These shall be discussed in more detail later. The transitions to differing states (from Table 2) are shown as arrows (transitions A to F) in the schematic diagram of Figure 2. The exact probabilities of such transitions are given by the off diagonal elements of Table 2 and also shown above or below the arrows in Figures 3 and 4.

Figure 3.

Diagram of transition probabilities for synchronous pairs: LeuTri and LeuPol.

Figure 4.

Diagram of transition probabilities for asynchronous pairs: PolTri, LeuMic and MicPol.

Figure 5.

Transition probabilities from (0, 0) and (1, 1) states for 6 species pairs.

Figure 6.

Transition probabilities from (0, 1) and (1, 0) states for 6 species pairs.

The transitions have the following intuitive interpretation and associated probability (sum), which are derived from the subcomponents of the transition matrices Q (see Table 2).

• A: transition of both species off to one species on: q(0, 0;0, 1) + q(0, 0, 1, 0)

• B: transition of both species on to one species off: q(1, 1;0, 1) + q(1, 1, 1, 0)

• C: species switching states: q(0, 1;1, 0) + q(1, 0; 0, 1)

• D: transition of one species off to both species off: q(0, 1;0, 0) + q(1, 0;0, 0)

• E: transition of one species on to both species on: q(0, 1;1, 1) + q(1, 0;1, 1)

• F: transition of one species on/off to both species off/on: q(0, 0;1, 1) + q(1, 1;0, 0)

In this chapter we shall demonstrate that transitions that lead towards both species being off or both species being on (states D, E or F), are considered to be synchronising. However, transitions that lead towards only one species being on or off (flowering) (A and B) and where within a species pair flowering switches (transitions C) are considered to be asynchronous.

Note that the probabilities of staying in the same state; e.g. both species continuing to be in a non-flowering state (a (0, 0) to (0, 0) transition); one species flowering off and the other species in the pair with flowering on, (a (0, 1) to (0, 1) transition); one species on the other in the pair off (a (1, 0) to (1, 0) transition); and both species continuing to flower (a (1, 1) to (1, 1) transition) are not shown on Figure 2. These to same states transitions, are given for each species, by the diagonal elements in the transition matrices (from previous to current states) in Table 2; and are also shown in Figures 3 and 4 as numbers (positioned next to the 4 states as boxes).

An examination of the transition probabilities for the species pairs in Table 2 shows that there is a significantly high propensity (probability) to remain in the same (bivariate) state as the previous state (see highlighted transition probabilities on the diagonals). For synchronous species pairs, such as LeuPol, and LeuTri the likelihood of species switching flowering state (states C), i.e. transition from one species flowering in a pair previous state = (0, 1) to the other species flowering, current state = (1, 0) never occurs (transition probability = 0.0000); or the likelihood of the transition from one species flowering to the other species flowering (i.e. a (1, 0) to (0, 1) transition) is rare (0.0000 ≤ transition probability ≤ 0.0085). For asynchronous species pairs such as LeuMic, MicPol, and PolTri, their switching probabilities are significantly higher in that at least one of the transition probabilities from (0, 1) to (1, 0); or from (1, 0) to (0, 1) is greater than 0.036, with associated probability ≥0.076.

Overall for synchronous pairs the probabilities of one species flowering to both or no species flowering, i.e. one off to both off, or one on to both on are high (>0.30). The latter are delineated by D and E transitions in Figure 2 and Table 4. Overall for asynchronous pairs there are high probabilities of both off (or on) to one off (or on). The latter transitions are delineated by A and B in Figure 2, with probabilities given in Tables 3 and 4.

Transition names	Description	Probability: sum of subcomponents	Threshold for synchrony	Threshold for asynchrony	Rules
A	Both off to one off	q(0, 0;0, 1) + q(0, 0;1, 0)	<0.30▼	≥0.30▲	P(A or B) > 0.8 for asynchrony
B	Both on to one on	q(1, 1;0, 1) + q(1, 1;1, 0)	<0.50▼	≥0.50▲
C ϕ	Switching	q(0, 1;1, 0) + q(1,0;0, 1)	<0.05▼	≥0.05▲
D ϕ	One off to both off	q(0, 1;0, 0) + q(1, 0;0, 0)	≥0.40▲	<0.40▼	P(D or E) > 0.65 for synchrony
E	One on to both on	q(0, 1;1, 1) + q(1, 0;1, 1)	≥0.40▲	<0.40▼
F	Both on (off) to both off (on)	q(0, 0;1, 1) + q(1, 1;0, 0)	≥0.08▲	<0.08▼

Table 3.

Descriptions and rules of (a) synchrony based on the transitions A-F.

^ϕ

Events or transitions C and F do not occur often.

Transition probability sums	P(A)	P(B)	P(A or B)	P(C)	P(D)	P(E)	P(D or E)	P(F)
Synchronous (S) pairs
LeuTri	0.232	0.181		0.008	0.266	0.659		0.081
LeuPol	0.293	0.238		0.000	0.501	0.318		0.141
Asynchronous (A) pairs
PolTri	0.254	0.556		0.076	0.379	0.121		0.037
LeuMic	0.333	0.510		0.171	0.085	0.309		0.000
MicPol	0.231	0.571		0.102	0.455	0.088		0.077
Neither S nor A
MicTri	0.165	0.426		0.036	0.332	0.300		0.038

Table 4.

Transition probabilities of events A to F for each species pair categorised into synchronous and asynchronous (or neither) species pairs.

In summary the transitions that lead to both species being off (no flowering) or both species being on (flowering) (transitions D, E or F), are considered to be synchronizing. However, transitions that lead to only one species being on or off (no flowering) (transitions A and B) and where a species pairs’ flowering status switches (transitions C) are considered to be asynchronous.

We now provide a rule for synchrony (or asynchrony) based on subcomponent (sums) of the transition probabilities derived from the B-MTD model:

• Two species are synchronous if P(D or E) > 0.65, i.e. P(one on to both on) + P(one off to both off) > 0.65,

• Two species are asynchronous if P(A or B) > 0.80, i.e. P(both off to one off) + P(both on to one on) > 0.8.

The transitions have the following interpretation and probabilities (Tables 3 and 4):

• A: transition of both species off (in the past state) to one species flowering (on) in the current state;

• B: transition of both species on to one species off;

• C: species switching states;

• D: transition of one species off to both species off;

• E: transition of one species on to both species on;

• F: transition of one species on/off to both species off/on.

According to the rules given in Table 3, the synchronous pairs are LeuTri and LeuPol (with P(D or E) > 0.65); asynchronous pairs are: PolTri, LeuMic and MicPol (with P(A or B) > 0.80) and a species pair that is neither synchronous nor asynchronous is MicTri.

In summary we have a simple rule for (a) synchrony, which in agreement with the work of [6] (see also [25]), using the synchronisation theory of Moran that:

• E. leucoxylon flowering is synchronous with both E. polyanthemos and E. tricarpa, but asynchronous with E. microcarpa.

• E. microcarpa is synchronous with none of three species; specifically it is asynchronous with both E. leucoxylon and E. polyanthemos (and has no relationship with E. tricarpa).

• E. polyanthemos flowering is synchronous only with E. leucoxylon; and asynchronous with both E. microcarpa and E. tricarpa.

• E. tricarpa flowering is synchronous only with that of E. leucoxylon; and is asynchronous with E. polyanthemos (and has no relationship with E. microcarpa).

We can view Figure 5 as the transition signatures from past states, where both species flowering is off or both species flowering is on, for synchronous pairings (LeuTri or LeuPol) and the asynchronous species pairs (PolTri, LeuMic and MicPol). Figure 6 likewise delineates transition signatures from past states, where only one species of the pair is flowering. These signatures (Figures 5 and 6) distinctly differ according to whether a species pair is synchronous or asynchronous.

For MicTri the associated sum of the probabilities for transitions A and B (both off/on to one off/on) is 0.591 (see Table 4), which is close to the threshold for synchrony of 0.65. Note that the more sophisticated MTDg modelling approach in Section 4.2 which incorporates covariates (mean temperature and rainfall) with interactions, shows that indeed E. microcarpa and E. tricarpa are synchronous (Tables 6 and 7), wherein the MTDg model allows for prior lag 1 to lag 12 month flowering effects and climate covariates (see also Table 7 and Figure 7).

Figure 7.

Synchrony relationships among the four eucalypts species.

4.2 Moran tests on residuals of the MTDg models incorporating climatic covariates

In this section synchronisation among species pairs is tested using Moran’s correlation method on the cross-residuals, based on MTDg models which incorporate both climate covariates and lagged effects of previous flowering. This work is based on [16], where MTDg models allowing interactions were fitted to the same four species. We present here only MTDg models with two covariates, namely, mean temperature and rainfall.

Parameters of the MTDg models are shown in Table 5. Significant lag effects of previous flowering states (lag j, where j = 1, ..., 12 months), and of the climatic covariates (meanT and rain) and their interaction (meanT*rain) are also given in Table 5. The estimated parameters for the MTDg models generally show a (positive) 1 month lag effect and 9, 11 and 12 months lag effects of previous flowering status (Table 5).

Species	lag 1	lag 9	lag 10	lag 11	lag 12	Temp variable	Rain	Temp × rain
E. mic	0.534	-	-	0.032 ϕ	0.275	0.136	-	-
E. poly	0.530	0.060	-	0.160	0.105	0.091	0.009	0.045
E. leu	0.611	-	-	0.124	0.042	0.202	-	-
E. tri	0.617	0.059	0.009	0.096	-	0.157	0.062	-

Table 5.

MTDg mixing probabilities of MeanT and rain models.

^ϕ

Covariate effects above 0.03 are considered significant.

- indicates cells with zero probabilities.

From Tables 5 and 6 we observe that mean diurnal temperature (meanT) has a significant effect on flowering for all species; rain impacts significantly only on E. tricarpa (Tri) and an interaction effect between rain and meanT exists for E. polyanthemos (Pol). Overall, flowering increases as temperature (MeanT) increases for E. microcarpa; and flowering decreases as temperature increases for both E. leucoxylon and E. tricarpa. Rainfall positively impacts the flowering of E. tricarpa (i.e. flowering increases with more rainfall). Interestingly E. polyanthemos exhibits increased flowering at low meanT when there is contemporaneous below average rainfall and at high meanT with above average rainfall (see the transition probabilities to flowering for the interaction effect of E. polyanthemos (i.e. (0.88, 0.12, 0.20, 0.96)) in Table 6.

Species	Climate effects	Previous flowering		Temperature		Rain		Temperature by rain interaction
	(temp/rain)	Off	On	Low 1	High 2	Less 3	More 4	Low/less	Low/more	High/less	High/more
E. mic	(+/−)	0.00	1.00	0.00	1.00	0.39	0.28	-	-	-	-
E. poly	Inter-action	0.01	1.00	0.00	0.34	0.94	0.03	0.88	0.12	0.20	0.96
E. leu	(−/+)	0.05	1.00	1.00	0.00	0.88	0.94	-	-	-	-
E. tri	(−/+)	0.00	1.00	1.00	0.00	0.00	1.00	-	-	-	-

Table 6.

Transition probabilities of flowering for the meanT and rain MTDg models.

¹

Cut point for low temperature states: MeanT ≤13.83°C.

²

Cut point for high temperature states: MeanT >13.84°C.

³

Cut point for less rain: rain ≤40.44 mm.

⁴

Cut point for more rain: rain >40.45 mm.

Note that ‘-’ indicates cells with zero probabilities.

In what follows we denote the species used to estimate the parameters for the MTD-based equation as the ‘Model species’ and the species fitted with these estimated parameters as the ‘Fitted species’. Table 7 gives the resultant significant Moran correlations based on the residual series from the MTDg-based model and fitted species equations. Significant Moran correlations from both the MTDg (and the EKF models show that (a)synchronous pairings found via the MTD and EKF models in [15, 16, 17, 18, 19] generally agree (Tables 7 and 8); refer also to Figure 7, where a solid line indicates synchronous pairs and a dashed line indicates asynchronous pairs of species.

Model species	mic	pol	leu	tri
Synchronous fitted species	tri (0.14)	leu (0.14)	pol (0.16)	mic (0.15)
			tri (0.11)
Asynchronous fitted species			mic (−0.14 ϕ )

Table 7.

Significant Moran correlations (in brackets) from the MTDg models.

^ϕ

A negative and significant correlation indicates an asynchronous species pair.

Model species	mic	pol	leu	tri
Synchronous fitted species	tri (0.12)	leu (0.19)	pol (0.18)	leu (0.26)
			tri (0.33)
Asynchronous fitted species	leu (−0.17 ϕ )	mic (−0.10 ϕ )

Table 8.

Significant Moran correlations (in brackets) from the EKF models.

^ϕ

A negative and significant correlation indicates an asynchronous species pair.

Table 7 shows significant positive MTDg-based correlations (P < 0.006) for the following (model species: fitted species) pairs—(LeuPol), (PolLeu), (LeuTri), (MicTri) and (TriMic), indicating that E. leucoxylon is synchronous with E. polyanthemos, in agreement with the rules of synchrony described earlier (Tables 3 and 4). E. leucoxylon is synchronous with E. tricarpa; and that E. microcarpa and E. tricarpa are synchronous. The synchrony of the latter species pair (MicTri) however, contrasts the results of Moran-based results on raw intensity profiles which indicate that E. microcarpa and E. tricarpa were neither synchronous or asynchronous (Table 4). It is noteworthy however, that for this species pairing, E. microcarpa and E. tricarpa (i.e. MicTri or TriMic), the associated sum of the probabilities for transitions D and E (one species off/on to both species off/on) is 0.591 (Table 4), which is close to the threshold for synchrony of 0.65 (Tables 3 and 4).

Tables 7 and 8 shows significant negative-based correlations (P < 0.001) for the following (model species: fitted species) pairs; (LeuMic), (PolMic) and (MicLeu) indicating that that E. leucoxylon is asynchronous with E. microcarpa and E. microcarpa is asynchronous with E. polyanthemos (only via the EKF-based residuals) (Figure 7 RHS); in agreement with the rule for asynchrony (Table 4) and Moran-based AR analysis of the flowering intensities.

Both the MTDg- and EKF-based models show that E. tricarpa is not asynchronous with E. polyanthemos (Tables 7 and 8). Note that for this species pairing E. tricarpa and E. polyanthemos (i.e. TriPol and PolTri) the associated sum of the probabilities for transitions P(A) and P(B) (both species off/on to one species off/on) is equal to 0.802 (Table 4), which is just above the to the threshold for asynchrony of 0.80.

4.3 Principal component analysis on the (λ, Q, d, s) parameters

In this section a novel approach which invokes a principal component analysis (PCA) of the resultant (λ, Q, d, s) parameters (Section 2.2) which details the weight λ, q, d and s parameters from the MTD (n = 4) models) is performed. The resultant two dimensional PCA axis plots (Figure 8) of the rotated (λ, Q, d, s)-based PCs provides an informative visualisation of the synchronous and asynchronous species groupings (of n > 2 species) allowing for interpretation of the main climate drivers and climatic profiles (e.g.+/− or ( −/+)) detailed in Table 6.

Figure 8.

Distances in the (λ, Q, d, s) parameters among the 4 species—without interaction terms (left) or with interaction terms (right).

The resulting parameters estimated from the MTDg models with and without interaction terms can be compared among all four species using Figure 8, which shows that the separation of E. tricarpa (−/+) and E. microcarpa (+/−) from other species along the horizontal axis 1, is due to the effect of mean temperature. Although E. leucoxylon is affected by the similar lag 1 and 11 month flowering terms as E. polyanthemos, E. leucoxylon (−/+) commences flowering at low temperature and shuts down at high temperatures. E. microcarpa begins flowering at high temperature (+/−). Figure 8 also displays the similarity (synchronicity) of E. leucoxylon and E. polyanthemos.

5. Discussion and conclusion

The highest degree of synchrony (via the B-MTD rules of synchrony, the MTD models and Moran AR method) occurs between E. leucoxylon and E. tricarpa; then followed by E. polyanthemos and E. leucoxylon which indicates the potential for intense competition for potential pollinators, and therefore the prospect for a high level of hybridization. Both these species pairs were shown to be synchronous by Keatley et al., [20]; with E. leucoxylon and E. tricarpa having 6 years of no overlap (and a long term mean synchrony value of 0.62); and E. polyanthemos and E. leucoxylon having 5 years of the 31 years (between 1940 and 1970) with no overlap (long term mean synchrony value of 0.51); as quantified in [20]. The degree of synchrony or overlap of flowering was however determined using the method outlined in [22] which measures the extent of overlapping in the flowering periods among pairs of individuals in a population.

E. leucoxylon is the only species to synchronise flowering with E. tricarpa, as shown by all three methods, namely the B-MTD rules of synchrony, MTD models and Moran’s AR method. Synchrony between E. leucoxylon and E. tricarpa, may be explained in terms of niche/competition and also facilitation may be a factor, due to their different modes of flower production. This agrees with the findings of [20]. Interestingly the MTD models discussed here (see also [6, 16, 24]) show that the climatic drivers or signature of E. leucoxylon and E. tricarpa is similar with respect to temperature, in that both exhibit decreased flowering with increased temperature.

Likewise E. leucoxylon is the only species to synchronise flowering with E. polyanthemos. E. leucoxylon and E. polyanthemos sometimes occur in the same geographical area; and earlier studies have shown they overlap significantly [20]. From the flowering behaviour indices of Keatley and Hudson in [23], E. leucoxylon and E. polyanthemos were shown to have temporally separated months of peak flowering, September and November, respectively; likewise their flowering commencement months May and October, respectively. These two species can occur in the same geographical area and their flowering period. Differentiation of these two species is based on their differing months of peak flowering as well as their separated months of most probable flowering; October and November, respectively. Likewise their flowering commencement months differ, May and October, respectively [23].

The least degree of synchrony (via the B-MTD rules of synchrony, the MTD models and Moran method) is shown in this chapter to occur between E. leucoxylon and E. microcarpa; then followed by E. polyanthemos and E. microcarpa. Our results agree with the findings in [20], which established that a cross between E. leucoxylon and E. microcarpa is impossible. In terms of climatic signatures: the flowering of E. microcarpa behaves differently from E. leucoxylon and E. tricarpa. E. microcarpa flowers at higher temperature and its flowering has a significant and positive relationship with flowering a year ago, refer also to the results reported in [23].

Eucalyptus tricarpa and E. polyanthemos were shown in this chapter also to be asynchronous (discordant or out of phase). This is in agreement with conclusions reported in [2]. The MTDg model found a significant interaction between two climate variables, mean temperature and rainfall on the flowering of E. polyanthemos. As flowering is viewed as either ‘off’ or ‘on’ this interaction appears to be delineating E. polyanthemos’ flowering period. It usually commences flowering in late spring—as mean temperature is increasing and rainfall is decreasing and ceases in early summer; just prior to the warmest mean temperature and lowest rainfall.

Specific temperature thresholds for commencement and for the cessation of flowering for the four species studied here, have been established, see [5, 7, 8]. For example, E. microcarpa was shown to flower at high temperatures, and E. leucoxylon and E. tricarpa both at lower temperatures. The flowering of E. polyanthemos was shown to be impacted by both rainfall and temperature, with increased flowering when conditions were either cool and dry, or hot and wet—indicative of a rainfall by temperature interaction.

Moran residual analysis and the B-MTD analysis described in this chapter showed that E. tricarpa and E. microcarpa did not exhibit a significant synchronous nor an asynchronous relationship. However, for this species pairing, the associated sum of the probabilities for transitions A and B (both off/on to one off/on) is 0.591, which is close to the threshold for synchrony of 0.65. Indeed the more sophisticated MTDg modelling approach which incorporates covariates (mean temperature and rainfall) with interactions, showed that E. microcarpa and E. tricarpa are synchronous, wherein the MTDg model allows for prior lag 1 to lag 12 month flowering effects and climate covariates.

SOM-based clustering [4] and Moran AR (2) tests also found that E. polyanthemos was asynchronous to E. microcarpa and E. tricarpa, in agreement with the extended Kalman filter (EKF)-based synchrony measures in [15, 21]. Note also it was demonstrated in [20] that E. polyanthemos and E. microcarpa have 25 years with no overlap (with a long term mean synchrony value of 0.29). Note that the more sophisticated MTDg modelling approach which incorporates covariates (mean temperature and rainfall) with interactions, showed that indeed E. microcarpa and E. tricarpa are synchronous, wherein the MTDg model allows for prior lag 1 to lag 12 month flowering effects and climate covariates.

Recently synchronisation of eucalypt flowering is shown to be a complex mechanism that incorporates all the flowering elements—flowering duration, timing of peak flowering, and the timing of start and finishing of flowering, as well as possibly specific climate drivers for flowering [4]. The four species studied were shown to be influenced by temperature and rainfall and as a consequence their flowering phenology will change in response to climate change. This in turn will have an impact on species interactions and community [4].

Extensions of the B-MTD models to allow for climate covariates and for the comparison of more than 2 species at a time (a so-called multivariate M-MTD) is the topic of future work. Other forthcoming research is to examine the timing and a/synchronisation of the within species phenostages of both budding and flowering. Refer to earlier work using wavelets [26] and Generalized Additive Model for Location, Scale and Shape (GAMLSS) [27] to model the relationship between climate (mean monthly minimum, maximum temperatures and rainfall) during bud development and the flowering cycles of Eucalyptus leucoxylon and E. tricarpa from the Maryborough region of Victoria between 1940 and 1962. Monthly behaviour (start, peak, finish, monthly intensity, duration and success) in budding and flowering was assessed using, as in this current chapter, the indices of Keatley in [23].