Multiple circular–circular regression

2.1 Model

Möbius transformation is a closed group transformation in the complex plane ℂ . For the set Ω = { z : | z | = 1 } , the Möbius transformation function f : Ω → Ω is defined by

f ( z ) = e i θ ( z − a ) 1 − a ¯ z ,

Let θ _y and θ _x be circular random variables. Taking y = exp ( i θ y ) and x = exp ( i θ x ) , the regression of y on x is

y = β 0 x + β 1 1 + β 1 ¯ x ε ,

The regression function, which is a form of Möbius transformation, is a mapping of unit circle | z | = 1 onto itself. Here, β₀ is a rotation parameter while β₁ is a fixed point in the complex plane. The predicted point y given x is the intersection of the unit circle with the line joining −x and β₁. This is shown in Figure 1. The case of | β 1 | > 1 is also discussed by Kato et al. (2008). In this case, first the fixed point is taken to be 1 β 1 ¯ . Then this point is joined to β 1 | β 1 | β 1 | β 1 | x ¯ . The intersection of this line with the unit circle is the predicted point y.

In this model, it can be seen that if | β 1 | is closer to 1, the function results in the y-values concentrating around β 1 | β 1 | . When | β 1 | = 0 , we have y = β 0 x ; which is just a rotation. Thus, the regression function defined above can be used for regression when there is only rotation and also in the case when there is reflection about any axis.

OPEN IN VIEWER

Figure 1:

Left: Regression model demonstrated for one regressor taking β₀=1; Right: Regression model demonstrated for two regressors

Identifiability of a model means that no two sets of parameters give the same distribution. This is important because otherwise there will be various choices of the parameter set which describe the same thing.

Kato et al. (2008) did not consider the identifiability issue of this model. However, the regression curve given by the Möbius transformation is identifiable. This is because for all | β 1 | < 1 , there exists only a single point β₁ inside the disc | z | < 1 through which the projection line (the line joining −x and β₁) is taken. The case for | β 1 | > 1 , as discussed in Kato et al. (2008), seems different from the case of | β 1 | < 1 . However, the geometry in this case can also be explained similarly. This can be shown by comparing the individual coordinates of the predicted point from both the geometries. Without loss of generality, we can take β₁ to be lying on the real axis. This is because any orientation of the circle is equivalent to any other orientation by a rotation. After changing the complex plane π to ℝ 2 , this is directly implied. The intersection of the line joining − x = ( − x 1 , − x 2 ) and β 1 = ( b , 0 ) , where b ∈ ℝ , with the unit circle is y = x 1 ( 1 + b 2 ) + 2 b 1 + b 2 + bx 1 , x 2 ( 1 − b 2 ) 1 + b 2 + bx 1 . Changing this in the complex form yields y = x + β 1 1 + β 1 ¯ x . Thus, both the geometries are similar. In both the cases, as the final projection point is unique, the model remains identifiable.

Mathematically, for a given point x and given β₁, taking β₀=1, the predicted point y is the point at the circle which lies on the line joining −x to β₁. Thus, other possible values of β 1 ′ ∈ π which give the same y for this given x can only be the points lying on this line. This is true for all x ∈ Ω . Hence, for all x ∈ Ω , the points x, β₁ and β 1 ′ should be collinear. But, if we take x ′ ( ≠ x + β 1 1 + β 1 ¯ x ) ∈ Ω , then the three points ( x ′ , β 1 , β 1 ′ ) will again be collinear if and only if β 1 = β 1 ′ . Thus, β₁ is unique and hence the model is identifiable. Here, β₀ is just a rotation parameter which can be dealt with by rotation directly as done in Theorem 7.2 (see the Appendix).

3.1 Model I (MCR1)

Let θ y , θ x j , j = 1 , … , k , be circular random variables. Our aim is to regress θ _y on θ x j 's. Let y = exp ( i θ y ) and x j = exp ( i θ x j ) for all j. Then, y can be modelled as

y = β 0 exp i arg Σ j = 1 k p j x j + β j 1 + β j ¯ x j ε ,

(3.1)

arg ( · ) ∈ [ 0 , 2 π ) ; β 0 , ε ∈ π ; | β 0 | , | ε | = 1 ; β j ∈ π ; and Σ j = 1 k p j = 1 ; 0 < p j < 1 for j = 1 , … , k .

(3.2)

The model is not defined for x = ( x 1 , x 2 , … , x k ) when Σ j = 1 k p j x j + β j 1 + β j ¯ x j = 0 . If k = 2 , then, this problem can be overcome if we define p 1 ≠ 0.5 . If k > 2 , then while checking the maximum of log-likelihood during data analysis, for a particular parameter set, such sample points should be ignored and the maximum likelihood should be divided by the number of sample points considered for any parameter set. The proposed model also becomes unidentifiable when | β j | = 1 for all j. This can be rectified by modifying the parameter space such that at least one | β j | ≠ 1 .

Note that Equation (3.1) is a multiple circular extension of circular–circular regression proposed by Kato et al. (2008).

The case of two regressors is shown in the right panel of Figure 1. Here O is the centre, and

y ′ = exp i arg py 1 + ( 1 − p ) y 2 , y = β 0 y ′ ,

The examples of such model can be the direction of bird migration regressed on wind direction and time of year when all the three variables are circular in nature.

The proof of identifiability for this model is given in Appendix A.

Theorem 7.1 can also be used to improve the fit of the circular–circular regression when there is only one circular covariate. The geometry is shown in the left panel of Figure 1, which is as follows:

Step 1: x ⇒ − x .

Step 2: − x ⇒ y ′ passing through β₁.

Step 3: Rotate y ′ by θ 0 = arg ( β 0 ) to obtain y. (In the figure, β 0 = 1 , and hence y ′ = y .)

In the case of bivariate circular regression as explained in this article, the right panel of Figure 1 shows the regression. The geometry can be described as follows:

Step 1: x ⇒ − x .

Step 2: − x ⇒ y 1 passing through β₁.

Step 3: − x ⇒ y 2 passing through β₂.

Step 4: Join y₁ and y₂ and take y ′ = py 1 + ( 1 − p ) y 2 .

Step 5: The regressed point y ′ is the projection of line joining centre O and y ′ on the circle.

Step 6: Rotate y ′ by α = arg ( β 0 ) to get y.

For the general case, we can obtain y 1 , y 2 , … , y k as given in Steps 3 and 4. Then we can take y ′ = Σ j = 1 k p j y j . Now, the regressed point will be the projection of line joining centre O and y ′ on the circle, which is y ′ . Then, as in Step 6, y can be obtained by rotation of y ′ by α anticlockwise.

Note: The bivariate case and the general case are explained only for | β j | ≤ 1 . The case of | β j | > 1 for any j can be dealt similarly only by incorporating the geometry as described in that case earlier for one regressor for getting the projection point corresponding to that regressor.

3.2.1 Multiple circular regression

The model described above is helpful in the case of multiple circular regressors. After specifying the direction of 0 angle for each covariate, the parameters p_js can be used to gauge the magnitude of effect of each covariate on the response. The higher the value of a particular p_j, the higher is the effect of that particular regressor x_j on the response y. This can be said because p_js are nothing but weight parameters. Thus, the model can also be used to check the relative effect of the regressors on the response variable.

In MCR also, new parameters corresponding to the covariates can be added as specified in single regressor case to get a better fit than before. The removal of the dependency of the parameter p_j's on the orientation of the covariates can be done by adding another parameter A_j for each covariate and hence, the new model can be written as the following:

y = β 0 exp i arg Σ j = 1 k p j A j x j + β j 1 + β j ¯ A j x j ε ,

(3.3)

The regression function can also be written in the form of tangent link function. In such case, the regression function can be written as

θ y = θ 0 + arg exp i arg Σ j = 1 k p j A j x j + β j 1 + β j ¯ A j x j + θ ε ,

(3.4)

where

A j ∈ Ω , A 1 = 1 ; arg ( · ) ∈ [ 0 , 2 π ) ; β 0 , ε ∈ π ; | β 0 | , | ε | = 1 ; β j ∈ π ;

p = ( p 1 , … p k ) ∈ R k : Σ j = 1 k p j = 1 ; 0 < p j < 1 for j = 1 , … , k .

The proof of identifiability of this model proceeds similarly as is done in Theorem 7.1 and Theorem 7.2.

3.2.2 Invariance and rotational equivariance of parameters

Circular regression models must be origin independent. Hence, for every rotation of the covariates and responses, a desirable feature is to get the same predicted mean direction by changing the model parameters. Also, p must be invariant of the rotations. We shall denote the parameters of the previous model by (β₀, A_j's, β_j's, p) and the new model parameters by (β₀(n), A_j(n)'s, β_j(n)'s, p(n)).

The predicted mean direction of an MCR1 model can also be written as

y = Σ j = 1 k p j β 0 A j x j + β j 1 + β j ¯ A j x j Σ j = 1 k p j β 0 A j x j + β j 1 + β j ¯ A j x j ,

(3.5)

where A₁=1. A change in the 0 angle is the same as multiplying by a unit complex number W. It can be seen directly from Equation (3.5) that if the response y is rotated and the new response is Wy, where | W | = 1 , there exists a model which gives the same prediction if all the parameters are the same except forβ₀. Thus, (Wβ₀, A_j's, β_j's, p) gives the same prediction in the new model as was given by (β₀, A_j's, β_j's, p) in the previous model.

When x₁ is rotated to Wx₁, then the new model parameters which give the same predicted mean direction are ( W ¯ β 0 , WA j ′ s , W β j ′ s , p ) .

When x j , j ≠ 1 , is rotated by W, then the new equivalent model parameters are ( β 0 , A 1 , W ¯ A 2 , … , W ¯ A k , β j ′ s , p ) . The meaning of p is completely preserved as it is invariant with respect to any rotation.

Fixing A 1 = 1 means fixing this value with respect to one of the covariates. The model does not depend on which covariate is chosen as x 1 because if in another case x r ( r ≠ 1 ) is chosen as the covariate such that A r = 1 is fixed, then the parameters with respect to the new model which give the same predicted mean direction are ( β 0 A r , A r ¯ A j ′ s , A r ¯ β j ′ s , p ).

The problem with the model MCR1 is that adding a new covariate results in adding four new parameters to the model: p j , A j , b 1 , j , b 2 , j , where β j = ( b 1 , j , b 2 , j ) T . Thus, we propose an alternative multiple circular–circular regression model having an increase in two parameters per covariate. We write

y = β 0 x * + β 1 1 + β 1 ¯ x * ε ,

(3.6)

x * = Σ j = 1 k p j A j x j | Σ j = 1 k p j A j x j | .

(3.7)

θ y = θ 0 + arg x * + β 1 1 + β 1 ¯ x * + θ ε .

(3.8)

The parameters p j 's, β 0 , β 1 , A j 's are defined as earlier. The interpretation of the parameters are the same and the rotational equivariance ofβ₀, p j 's, A j 's is exactly similar to what was discussed in the previous section. By fixing A 1 = 1 , we have locked the location of x₁, hence it can easily be seen that β₁ is rotationally equivariant with x₁. The tangent link function in Downs and Mardia (2002) and the Möbius transformation in Kato et al. (2008) are mathematically equivalent for one covariate. It can be seen that the equivalence of the MCR2 model with Downs and Mardia (2002) is maintained if we consider x * as defined in MCR2 as the covariate in Downs and Mardia (2002). The Figure for the regression function in case of two covariates is shown in Figure 2. The identifiability of this model is proved in Appendix B.

OPEN IN VIEWER Figure 2:

Regression model demonstrated for two covariates

3.3.1 Invariance and rotational equivariance of parameters

As explained in Section 3.2.2, there must exist an equivalent model with respect to the rotation of the circular variables in MCR2 also. From Equation (3.5), it can be directly seen that if y is rotated by W, then the new model parameters giving the same predicted mean direction as the previous one will be (Wβ₀, A_j's, β₁, p_j's).

Note that Equation (3.6) can also be written as

y = β 0 x * + β 0 β 1 1 + β 1 ¯ x * ε .

Thus, the equivalent new model parameters when x₁ is rotated to Wx₁ are ( W ¯ β 0 , A 1 , WA j ( j ≠ 1 ) , W β 1 , p j ′ s ) . When, x j ( j ≠ 1 ) is rotated to Wx j ( j ≠ 1 ) , then the equivalent new model parameters are ( β 0 , A 1 , … , A j ¯ , … , A k , β 1 , p j ′ s ) .

The choice of A 1 = 1 does not change the prediction because the equivalent new parameters will then be ( β 0 A r , A r ¯ A j , A r ¯ β 1 , p j ′ s ) .

4.1 Maximum likelihood estimators

When β j s are all known for j = 1 , 2 , … , k , and the angular error follows the wrapped Cauchy distribution WC ( ψ ) , the estimation of ψ andβ₀ can be done by the method proposed by Kent and Tyler (1988). This estimation can also be done by the method of moments as described in Kato et al. (2008). When the angular error follows asymmetric generalized von Mises distribution, then the MCR1 model is of the form

θ y = θ 0 + θ E + arg exp i arg Σ j = 1 k p j A j x j + β j 1 + β j ¯ A j x j + θ ε ′

θ y = θ 0 + θ E + arg x * + β 1 1 + β 1 ¯ x * + θ ε ′ .

If the angular error follows VM ( 0 , κ ) , then the MLEs ofβ₀ and κ reduce to

θ ^ 0 = arg ( C + iS ) ,

S = Σ r = 1 n sin θ y r − arg Σ j = 1 s p j x r , j + β j 1 + β j ¯ x r , j , C = Σ r = 1 n cos θ y r − arg Σ j = 1 s p j x r , j + β j 1 + β j ¯ x r , j , and R = Σ r = 1 n cos θ y r − θ ^ 0 − arg Σ j = 1 s p j x r , j + β j 1 + β j ¯ x r , j .

S = Σ r = 1 n sin θ y r − arg x r * + β 1 1 + β 1 ¯ x r * , C = Σ r = 1 n cos θ y r − arg x r * + β 1 1 + β 1 ¯ x r * , and R = Σ r = 1 n cos θ y r − θ ^ 0 − arg x r * + β 1 1 + β 1 ¯ x r * .

L = − n log 2 π − n log ( I 0 ( κ ) ) + κ Σ r = 1 n cos ( θ y r − θ ^ y r ) ,

(4.1)

Now, for the first model, Σ r = 1 n cos ( θ y r − θ ^ y r ) is maximized when θ 0 is the circular mean of θ y r − arg Σ j = 1 s p j x r , j + β r , j 1 + β j ¯ x r , j s, and hence this is equivalent to maximizing

Q = Σ r = 1 n y r y c , r

y c , r = Σ j = 1 s p j x r , j + β r , j 1 + β j ¯ x r , j Σ j = 1 s p j x r , j + β r , j 1 + β j ¯ x r , j .

In case of the second model, Σ r = 1 n cos ( θ y r − θ ^ y r ) is maximized when θ 0 is the circular mean of θ y r − arg x r * + β 1 1 + β 1 ¯ x s, and hence this is equivalent to maximizing

Q = Σ r = 1 n y r y c , r

y c , r = x r * + β 1 1 + β 1 ¯ x r * .

OPEN IN VIEWER Figure 3:

Plots for Likelihood Function. [Top row:] Graphs of likelihood (L) against (Left:) A θ 0 and (Right:) b₁. [Middle row:] Graphs of likelihood (L) against (Left:) b₂ and (Right:) b₂ near the local maximum. [Bottom row:] Graph of likelihood (L) against p

Another very useful alternative for finding the global maximum is to use the bootstrap restarting method discussed in Wood (2001). The steps required for approximating the MLE by bootstrap restarting for our set-up are the following:

Step 1: Take an initial parameter set q₀.

Step 2: Starting from q₀, find a local maximum q_i and calculate L ( q i , x , y ) , the log-likelihood function.

Step 3: Take a bootstrap sample ( y * , x * ) from the dataset.

Step 4: Find q i * taking q i as the starting parameter.

Step 5: Check if L ( q i * , x , y ) > L ( q i , x , y ) . If it is so, then take q 0 = q i * = q i + 1 , else take q 0 = q i and repeat the steps from 1 to 5 until convergence.

As already mentioned, for the parameter set, taking the parameters related to are enough because we can get the MLE of β0 and κ directly from these parameters.

If angular error follows VM ( 0 , κ ) , then − E d 2 L d κ 2 = n d d κ A ( κ ) , − E d 2 L d θ 0 2 = nkI 1 ( κ ) , and

− E d 2 L d θ 0 d κ = − E d 2 L d θ j d κ = − E d 2 L dr j d κ = − E d 2 L dp j d κ = 0 .

Similarly, if the angular error follows C * ( ψ ) , then

− E d 2 L d ψ d θ 0 = − E d 2 L d ψ d θ k = − E d 2 L d ψ dr k = − E d 2 L d ψ dp k = 0 .

4.3 Improving the fit in case of single regressor

In case of circular–circular regression of Kato et al. (2008), one circular variable is regressed on only one circular regressor. The fit of the regression can be improved by using the model proposed in the article by adding one new parameter β j and correspondingly p j at each step and calculating the accuracy in terms of the sum of cosines of the angles between the fitted and observed value. A new pair of parameters ( β j , p j ) can be added till the improvement in fitting is more than a predefined value. This is similar to adding the powers of x in the linear regression case so that the fitting is improved. This case can also be equated to the case when in multiple linear regression, orthogonal polynomials of the regressor x are added and penalty is calculated until the improvement in fitting is less than a predefined value. Here, also a method similar to adjusted R 2 may be used to calculate the optimum number of β j s.

Here we intend to define the circular version of R 2 statistic, and denote it by R C 2 . For illustration, we assume that arg ( ε j ) follows VM ( 0 , κ ) distribution, the von Mises distribution with location 0 and concentration parameter κ , with probability density function (pdf) given by

f 0 , κ ( z ) = 1 2 π I 0 ( κ ) exp κ cos ( z ) ,

I j ( κ ) = 1 π ∫ 0 π exp ( κ cos z ) cos ( j z ) dz ,

κ ^ = A − 1 1 n Σ r = 1 n cos θ ^ ε r , if Σ i = 1 n cos θ ^ ε r > 0 , 0 otherwise ,

R C 2 = 1 2 1 + 1 n Σ r = 1 n cos θ ^ ε r .

However, a more objective method to check the significance of change can be through the log-likelihood. Given the model and assuming angular error following VM ( 0 , κ ) distribution, the likelihood based testing for the addition of parameters can be done in the following way. Under the model given in equations (3.1) and (3.2), the log-likelihood ratio test statistic for H 0 : p s = 0 against H 1 : p s ≠ 0 can be used to check for any significant difference by adding new parameter β s .

Let ℓ 0 denote the maximum of log-likelihood of y given x under H 0 , and ℓ 1 denote the same under H 1 . For using Likelihood Ratio Test (LRT), we need to find the distribution of Λ = − 2 ( ℓ 0 − ℓ 1 ) . However, if p s = 0 , then the parameter β s becomes unidentifiable. This is because β s can take any value when p s = 0 . Moreover, even if p s ≠ 0 and β s is equal to any one of β 1 , … , β s − 1 , then the parameters p j 's become non-identifiable. Hence, the regularity condition of model identifiability for the standard asymptotic chi-squared distribution of the LRT statistic is not satisfied by the model under null hypothesis. Thus, we need an alternative approach.

McLachlan (1987) discussed a problem where the LRT statistic was bootstrapped to estimate the number of components in a normal mixture. Here, we consider a similar methodology in our article. First, based on the sample, MLEs of the model parameters were found using ( s − 1 ) number of β 's, and for the same sample, MLEs were found by using s number of β 's. Based on the estimates from the first model, K bootstrap samples of the same sample size as of the original sample were generated. Our case is a bit different from McLachlan (1987) as the bootstrap sample is bivariate ( θ x , θ y ) . Let the distribution of θ x be known. Then, an independent bootstrap sample for θ x is generated from the distribution of θ x . This generation of θ x from its distribution is considered because then the fitted model will not only be good for the sample data but it will be good over the whole sample space of θ x . Based on each of these θ x , θ y is generated using the estimates Φ ^ based on the model under null hypothesis. This bootstrap sample hence obtained is used to estimate the parameters under both the models. The same process is repeated independently K times, and Λ is calculated each time. The value of the jth order statistic of the K replications can be taken as an estimator of the quantile of order j / ( K + 1 ) , and the P -value can be assessed with respect to the ordered bootstrap replications of Λ .

Under the assumption of von Mises distribution for angular error,

− 2 ( ℓ 0 − ℓ 1 ) = 2 n log ( I 0 ( κ ^ s − 1 ) − κ ^ s − 1 A ( κ ^ s − 1 ) − log ( I 0 ( κ ^ s ) + κ ^ s A ( κ ^ s ) .

d d κ log ( I 0 ( κ ) ) − κ A ( κ ) = − κ dA ( κ ) d κ

Other circular distributions for ε will work as well. One such example can be the wrapped Cauchy distribution as described earlier by Kato et al. (2008). When angular error follows the wrapped Cauchy distribution with 0 mean, then arg( ε ) ∼ C * ( ψ ) , where 0 ≤ ψ ≤ 1 . In this case, the model is a good fit if ψ is closer to 1.

Standard model selection criteria like the bias-corrected Akaike information criterion (AICC) can also be used to select the appropriate candidate for the model; see Lund (1999) in this context. For von Mises error distribution,

AICC = 2 n log I 0 ( κ ^ ) − 2 n κ ^ + n ( n + l ) n − l − 2 ,

4.4 Donut-plot

A new type of plot is also introduced in this present article. This plot is called donut-plot. The problem with normal x versus y graph in circular cases is that it depends on the choice of ‘0’. The spokeplot, introduced by Zubairi et al. (2008) and as shown in Kato et al. (2008), is a very good pictorial representation of observed versus estimated values. But, when the sample size is large, usually a lot of lines criss-cross each other and in that case the plot is not very clear. See the spokeplot at the left in the bottom row of Figure 5.

The donut-plot is used to plot the estimated value and the corresponding observed value. Let the estimated angle be θ y e and the observed angle be θ y . Then, this plot shows the point ( 1 + cos ( θ y e − θ y ) ) ( cos θ y e , sin θ y e ) . Thus, if the fit is good, then more points will be near the circumference of the circle with centre at 0 and radius =2, while a bad fit results in more points inside the unit radius disc. Here, as in spokeplot, we do not need to delete the outliers. The donut-plot can be further improved by using different coloured dots for clockwise and anticlockwise deviations. This is because it also shows the direction in which the estimated points are deviating from the observed values. In the donut-plots given in Section 5, we have used solid circles for clockwise deviations and empty circles for anticlockwise deviations. Thus, the donut-plot provides an alternative to the x versus y graph and the spokeplot.

4.5 Simulation

Simulation studies are done for the MCR case when the angular variable θ y is regressed on two independent angular covariates θ x 1 and θ x 2 . The simulation is performed both for the MCR1 and MCR2 models. The method employed for finding the global maximum in this article is by giving different initial values to the optim function in R. The simulation results, for the MCR1 model, are shown in Table 1 and the simulation results for the MCR2 model are shown in Tables 2 and 3. Table 1 shows the estimates when the data is generated from the MCR1 model and estimates are the MLE based on the MCR1 model. Table 2 shows the estimates when the data is generated from the MCR2 model and estimates are the MLE based on the MCR2 model. In the tables, the sample standard errors (se's) for the estimates of linear variables and sample circular variances for the estimates of circular variables are given in parentheses. The definition of sample circular variance is the one used in equation (2.3.3) of Mardia and Jupp (2000). The sample circular variance is 1 − R ¯ , where R ¯ is the mean resultant length of the unit complex numbers corresponding to the MLEs of the circular parameters. The mean of the estimates are the mean of MLE found from 1 000 runs and the standard deviations (SDs) are the SD of the MLEs found from the 1 000 runs. Consequently the se's are calculated. We carried out extensive computations, but reported the results for only three parameter combinations and for sample sizes 20 (small sample case, similar to Example 1 in Section 5.1), 80 (moderate sample size, close to the situation of Example 2 of Section 5.2), and 200 (large sample case). The estimates are quite close, even for a small sample case, and the estimates improve as the sample size increases. However, the estimate of κ is not quite good when the sample size is small (like 20), the estimate improves with the increase in sample size. For the MCR2 model, we have also changed p and it can be seen from Table 2 and Table 3 that for low p, the se's of the MLEs are high.

In Table 1, β₁ and β 2 are represented by ( b 11 , b 12 ) and ( b 21 , b 22 ) , respectively, A θ 0 represent the argument of A 2 , θ 0 gives the argument ofβ₀ and p represents the weight of the first covariate. In Table 2, ( b 1 , b 2 ) represents β₁, θ 0 is the argument ofβ₀, p is the weight p 1 of the first covariate and A θ 0 is the argument of A 2 .

OPEN IN VIEWER

Table 1:

Estimates (and se in parentheses) of the parameters for MCR1

	Estimates
True value	n=200	n=80	n=20
κ=3	3.095 (0.009)	3.318 (0.015)	5.621 (0.080)
b11=0.7	0.707 (0.002)	0.704 (0.005)	0.572 (0.012)
b12=0	0.000 (0.002)	0.002 (0.005)	−0.003 (0.012)
b21=0.5	0.459 (0.005)	0.389 (0.005)	0.225 (0.011)
b22=0	0.005 (0.006)	0.003 (0.009)	0.001 (0.013)
p=0.4	0.411 (0.001)	0.427 (0.001)	0.449 (0.002)
Aθ0=0	0.008 (0.036)	0.003 (0.104)	0.009 (0.302)
θ0=0	−0.005 (0.018)	0.001 (0.049)	−0.002 (0.182)
κ=3	3.101 (0.008)	3.337 (0.015)	5.634 (0.090)
b11=0.7	0.707 (0.002)	0.710 (0.005)	0.585 (0.012)
b12=0	0.006 (0.002)	0.004 (0.005)	0.003 (0.012)
b21=0.5	0.461 (0.002)	0.386 (0.005)	0.217 (0.011)
b22=0	−0.015 (0.005)	0.015 (0.009)	0.016 (0.014)
p=0.4	0.411 (0.001)	0.430 (0.001)	0.449 (0.002)
Aθ0=π/3	1.024 (0.034)	1.067 (0.109)	1.085 (0.336)
θ0=0	0.015 (0.017)	−0.020 (0.053)	−0.025 (0.206)
κ=1.5	1.564 (0.005)	1.729 (0.008)	2.837 (0.038)
b11=0.7	0.706 (0.005)	0.670 (0.008)	0.517 (0.015)
b12=0	0.002 (0.005)	−0.006 (0.008)	0.013 (0.014)
b21=0.5	0.388 (0.005)	0.282 (0.008)	0.187 (0.014)
b22=0	0.001 (0.009)	0.010 (0.012)	−0.009 (0.015)
p=0.4	0.432 (0.001)	0.451 (0.001)	0.470 (0.001)
Aθ0=0	0.004 (0.102)	0.032 (0.222)	0.010 (0.338)
θ0=0	−0.008 (0.050)	−0.022 (0.116)	−0.006 (0.239)

OPEN IN VIEWER

Table 2:

Estimates (and se in parentheses) of the parameters for MCR2

	Estimates
True value	n=200	n=80	n=20
κ=3	3.08 (0.009)	3.21 (0.015)	4.32 (0.051)
b1=0.3	0.28 (0.001)	0.25 (0.003)	0.14 (0.007)
b2=0.0	0.00 (0.003)	−0.01 (0.005)	−0.01 (0.008)
p=0.15	0.16 (0.001)	0.14 (0.002)	0.24 (0.003)
θ0=0	−0.01 (0.07)	0.03 (0.2)	0.02 (0.55)
Aθ0=0	0.01 (0.07)	−0.02 (0.2)	−0.03 (0.55)
κ=3	3.06 (0.008)	3.22 (0.015)	4.32 (0.053)
b1=0.6	0.60 (0.001)	0.59 (0.002)	0.52 (0.004)
b2=0.0	0.00 (0.002)	0.00 (0.003)	−0.03 (0.006)
p=0.25	0.25 (0.001)	0.25 (0.002)	0.28 (0.003)
θ0=0	0.00 (0.01)	−0.00 (0.04)	−0.04 (0.15)
Aθ0=0	0.00 (0.01)	−0.00 (0.04)	0.02 (0.19)
κ=3	3.08 (0.008)	3.22 (0.015)	4.23 (0.051)
b1=0.3	0.30 (0.001)	0.30 (0.002)	0.28 (0.009)
b2=0.0	0.00 (0.002)	0.00 (0.003)	0.02 (0.005)
p=0.3	0.30 (0.001)	0.30 (0.001)	0.32 (0.003)
θ0=0	0.01 (0.01)	−0.02 (0.02)	0.02 (0.13)
Aθ0=π/3	1.04 (0.01)	1.06 (0.02)	1.04 (0.13)

OPEN IN VIEWER

Table 3:

Estimates (and se in parentheses) of the parameters for MCR2

	Estimates
True value	n=200	n=80	n=20
κ=1.5	1.54 (0.005)	1.63 (0.008)	2.13 (0.022)
b1=0.6	0.59 (0.002)	0.57 (0.003)	0.43 (0.012)
b2=0.0	0.00 (0.004)	0.00 (0.006)	0.00 (0.013)
p=0.3	0.30 (0.001)	0.31 (0.002)	0.36 (0.003)
θ0=0	0.00 (0.02)	−0.01 (0.07)	0.02 (0.36)
Aθ0=0	0.00 (0.02)	0.00 (0.07)	0.03 (0.36)
κ=1.5	1.54 (0.005)	1.61 (0.008)	2.19 (0.023)
b1=0.3	0.29 (0.002)	0.28 (0.003)	0.25 (0.008)
b2=0.0	0.00 (0.003)	−0.00 (0.004)	−0.04 (0.010)
p=0.3	0.30 (0.002)	0.31 (0.002)	0.36 (0.003)
θ0=0	−0.02 (0.03)	0.00 (0.09)	−0.07 (0.23)
Aθ0=π/3	1.06 (0.03)	1.04 (0.09)	1.05 (0.22)
κ=1.5	1.53 (0.005)	1.62 (0.008)	2.13 (0.022)
b1=0.3	0.24 (0.003)	0.18 (0.006)	0.10 (0.012)
b2=0.0	0.00 (0.005)	0.00 (0.007)	0.00 (0.011)
p=0.15	0.17 (0.002)	0.20 (0.003)	0.32 (0.004)
θ0=0	0.01 (0.20)	0.01 (0.43)	0.07 (0.69)
Aθ0=0	−0.01 (0.19)	−0.00 (0.42)	−0.05 (0.69)

4.5.1 Simulation for robustness

In this article, we have considered von Mises distribution for angular error. Kato et al. (2008) used the wrapped Cauchy distribution for angular error. Cauchy distribution is a heavy-tailed distribution. Thus, we have compared the case when the angular error follows von Mises distribution with zero mean to the case when the angular error follows the wrapped Cauchy distribution with real parameter. We have simulated the data from the MCR2 model for this comparison. We have used an accuracy measure mentioned below for checking which among the wrapped Cauchy and von Mises distribution is more robust distribution for the error. For this, first we have taken the true distribution of angular error to be von Mises and then found the accuracy at different sample sizes by assuming both the wrapped Cauchy and von Mises distribution for the angular error. Then, we have taken the wrapped Cauchy distribution as the true distribution for the angular error and then again found the accuracy assuming both the wrapped Cauchy distribution as well as von Mises distribution for angular error.

OPEN IN VIEWER

Table 4:

Accuracy values when angular error follows von Mises distribution

,	n=20		n=80		n=200
κ	VM	WC	VM	WC	VM	WC
1.5	121.17	119.87	49.77	48.71	14.56	12.62
3	163.15	162.19	66.05	64.97	17.44	16.34
4.5	176.71	176.26	79.97	70.47	18.22	17.65
6	182.71	182.49	73.35	73.06	18.67	18.21

For all the simulations, true parameter values are fixed at β 1 = 0.3 , p = 0.3 , β 0 = 1 and A 2 = 1 . In Table 4, the error follows von Mises distribution with zero mean direction and different values of the concentration parameter ( κ = 1.5 , 3 , 4.5 , 6 ) . For each of these cases, 1 000 simulations are done at three sample sizes ( n = 20 , 80 , 200 ) . The robustness is measured with respect to the total accuracy in prediction. The total accuracy is defined as the sum of cosines of the angles between the predicted mean direction and the observed mean direction, that is,

Accuracy = Σ i = 1 n cos ( θ y r − θ ^ y r ) .

The mean value of the Accuracy over 1 000 simulations is used to check the robustness. This accuracy measure is defined using circular distance, see Sengupta et al. (2013) for details. In Table 5, the distribution for angular error is considered to be the wrapped Cauchy with different values of parameter ( ρ = 0.2 , 0.4 , 0.6 , 0.8 ) at the above mentioned sample sizes. For each of the cases, 1 000 simulations are done and the total accuracy is found by assuming that the angular error follows the wrapped Cauchy and then by assuming that the angular error follows the von Mises distribution. Then, we have taken the case when angular error follows the wrapped t distribution with the scale parameter 1 and mean as 0 at different degrees of freedom ( ν = 1 , 3 , 5 , 7 ) at different sample sizes. Then, we have obtained the mean total accuracy based on 1 000 simulations by assuming first that the angular error follows wrapped Cauchy distribution and again assuming that the angular error follows von Mises distribution. The results are given in Table 6. Note that ν = 1 for the wrapped t reduces to wrapped Cauchy with ρ = e − 1 ≃ 0.368 , and the performance of the wrapped t for ν = 1 in Table 6 is expected to be in between those corresponding to ρ = 0.2 and ρ = 0.4 in Table 5, closer to the case corresponding to ρ = 0.4 . Our results are obtained accordingly.

From the results in Tables 4, 5 and 6, it can be seen that the von Mises distribution not only gives better results than the wrapped Cauchy distribution when the true density of the angular error is von Mises but also when the true density is the wrapped Cauchy or the wrapped t . This might be for our specific choice of Accuracy, which is a term in the expression of the log-likelihood for von Mises distribution we maximize to find the MLE of von Mises parameters; see Equation (4.3).

OPEN IN VIEWER

Table 5:

Accuracy when angular error follows the wrapped Cauchy distribution

,	n=20		n=80		n=200
ρ	VM	WC	VM	WC	VM	WC
0.2	46.91	45.70	23.76	20.96	10.09	7.12
0.4	82.53	81.33	34.98	34.28	11.51	9.65
0.6	121.29	121.23	49.49	48.52	14.02	12.73
0.8	160.41	159.73	64.52	64.24	16.72	16.08

OPEN IN VIEWER

Table 6:

Accuracy when angular error follows the wrapped t distribution

,	n=20		n=80		n=200
ν	VM	WC	VM	WC	VM	WC
1	68.44	54.45	30.50	23.44	10.93	7.26
3	89.22	88.19	38.25	36.59	12.03	9.70
5	99.42	98.29	41.75	40.16	12.63	10.76
7	104.39	103.59	43.50	41.85	13.00	11.29

5.1 Example: Single regressor

The data analysis for the case of one circular regressor and one circular response variable is used here to show the advantage of the model in improving the fit in the case of a single regressor. The data were also used in Kato et al. (2008). The wind direction was measured each day at a weather station in Milwaukee for 21 consecutive days at 6 a.m. and 12 noon (Johnson and Wherly (1977), Table 2). The response variable is the wind direction at 12 noon and the covariate is the wind direction at 6 a.m. Watson's test (1983) indicates no significant departure of the angular error from VM ( 0 , κ ) distributions in both one parameter and two parameter cases at 1% significance level (p-values are 0.248 and 0.523, respectively).

The MLEs for the parameters when only one β j is used, come out to be: β ^ 1 = 0.055 + 0.259 i , β ^ 0 = 0.753 + 0.658 i , κ ^ = 1.113 and R c , 1 2 = 0.743 . Here, the estimate of β₁ is obtained using the estimates of r and θ 1 .

When both β₁ and β 2 are considered in the model, the MLE for the parameters come out to be β ^ 1 = − 0.637 +0.597i, β ^ 2 = 0.648 + 0.113 i , β ^ 0 = 0.878 + 0.478 i , p ^ = 0.501 , κ ^ = 1.710 and R c , 2 2 = 0.823 . Here, β₁ and β 2 are estimated by taking β 1 = r 1 e i θ 1 and β 2 = r 2 e i θ 2 and then estimating the values of r 1 , r 2 , θ 1 , θ 2 . The estimates of all the parameters are given in Table 7.

The estimates of SD of the parameters reported in the tables are found by taking the inverse of the hessian matrix at the MLEs.

OPEN IN VIEWER

Table 7:

Estimates of parameters (SD in parentheses) for the Milwaukee wind data when (a) only β₁ is used; (b) β₁ and β₂ are used

only β1 is used	β1 and β2 are used
r→ 0.265 (0.192)	r1→ 0.658 (0.004)
θ1→ 1.362 (0.539)	θ1→ 0.171 (0.026)
θ0→ 0.719 (0.297)	r2→ 0.874 (0.003)
κ→ 1.113 (0.380)	θ2→ 2.389 (0.026)
	p→ 0.499 (0.002)
	θ0→ 0.501 (0.249)
	κ→ 1.710 (0.478)

By looking at R c , 1 2 and R c , 2 2 , it can be said that the second model fits the data better than the first one. The graph between circular distance ( 1 − cos ( θ y e − θ y )), and predictors for both the cases are shown in the top row of Figure 4. Comparing the two plots, it can be said that the fit is better in the second case. Then, the bootstrapping procedure for Λ was used to approximate the true distribution of Λ under the null hypothesis as described in Section 4.3 with K = 25 . First, Watson's test was performed on the data and it could not reject uniformity for θ x . Then, θ y 's were generated based on θ x 's generated from uniform distribution over [ 0 , 2 π ) . The estimate of 0.95-sample quantile based on 25 such bootstrapping comes to be 12.629 and the 0.9-sample quantile estimate comes to be 9.137. The observed value of Λ from the sample is 9.265. Hence, the choice of the first model is not rejected at 5% significance level while the null hypothesis is rejected and the second model is chosen at the 10% significance level. The AICC statistic as defined in Lund (1999) is −3.126 and −1.775 for the first and the second cases, respectively, showing the preference for the first model, but the difference is not large.

The spokeplots are shown in the middle row of Figure 4, while the proposed donut-plots are shown in the bottom row of Figure 4. The left figure in each row corresponds to the model with only one β, and the right hand figure is for two βs. The R codes of this data analysis can be obtained from the authors upon request.

OPEN IN VIEWER Figure 4:

Plots for Milwaukee data. [Top row:] Graph of circular distance versus predictors when (Left:) only β₁ is used and when (Right:) β 1 , β 2 are used. [Middle row:] Spokeplot when (Left:) only β₁ is used and when (Right:) β₁ and β 2 are used. [Bottom row:] Donut-plot when (Left:) only β₁ is used and when (Right:) β₁ and β 2 are used

We consider a similar weather data with one circular response and two circular covariates to illustrate the MCR application of the model. The data used is the wind direction data in Ambikapur district of Chhattisgarh from September 2011 to June 2012. The data is taken from the website of Agricultural Meteorology Division, India Meteorological Department, Ministry of Earth Sciences (http://www.imdagrimet.gov.in/winddirection). The dataset consists of 81 points, each corresponding to a different date and a particular hourly time corresponding to the date. The circular covariates considered here are the the date of the year and the time of the day. For the time of the day, 00:00 hrs is considered to be the zero angle; for the time of the year, 8th June (the date corresponding to the first data point) is considered to be the zero angle. Data analysis is done for both the MCR1 and MCR2 models.

OPEN IN VIEWER

Figure 5:

Plots for Chhattisgarh data. [Top row:] Graph of circular distance versus predictors. [Bottom row:] (Left:) Spokeplot; (Right:) Donut-plot

OPEN IN VIEWER Figure 6:

Plots for Chhattisgarh data.(MCR2) [Top row:] Graph of circular distance versus predictors. [Bottom row:] (Left:) Spokeplot; (Right:) Donut-plot

In case of MCR1 model, the goodness of fit test using Watson's test (1983) indicated no significant departure from the von Mises distribution for the angular error at 1% significance level (p-value = 0.42 ). Consequently, the MLEs of the parameters ( β 0 , β 1 , β 2 , p , κ ) are obtained. Here, β₁ is the parameter corresponding to time of the day, and β 2 is the parameter corresponding to the date of the year. The MLEs of the parameters are: β ^ 0 = − 0.963 + 0.271 i , β ^ 1 = 1.921 − 1.773 i , A ^ 2 = − 0.218 + 0.976 i , β ^ 2 = − 0.1 + 1.136 i , p ^ = 0.469 κ ^ = 1.293 . The graph between circular distances and response values is shown in the top row of Figure 5. The spokeplot excluding the outliers is shown in the left of the bottom row, while the donut-plot is shown in the bottom right of Figure 5. Here, β₁ and β 2 are estimated by first taking β 1 = r 1 e i θ 1 and β 2 = r 2 e i θ 2 and then estimating the values of r 1 , r 2 , θ 1 , θ 2 . The estimates of all the parameters are given in the left panel of Table 8.

OPEN IN VIEWER

Table 8:

Estimates of parameters (SD in parentheses) of the Chhattisgarh data when there are two covariates

MCR1	MCR2
b1→ 1.921 (0.257)	b1→ 0.076 (0.237)
b12→ −1.773 (0.192)	b2→ −1.222 (0.096)
b21→ −0.100 (0.202)	p→ 0.247 (0.073)
b22→ 1.136 (0.032)	θ0→ 5.627 (0.160)
p→ 0.469 (0.004)	Aθ0=→ 4.878 (0.138)
Aθ0=→ 1.791 (0.173)	κ→ 1.05 (0.191)
θ0→ 2.867 (0.134)
κ→ 1.293 (0.208)

When data analysis is done by using teh MCR2 model, Watson's test (1983) indicates that there is no significant departure from the assumption of the angular error following von Mises distribution 1% level (p-value = 0.23 ). Then the MLEs of the parameters are obtained. The MLEs for this case are tabulated in the right panel of Table 8. The graph between circular distances and response values is shown in the top row of Figure 6. The spokeplot excluding the outliers is shown in the left of bottom row, while the donut-plot is shown in the bottom right of Figure 5. The R codes of this data analysis can be obtained from the authors upon request.

7.1 Appendix A

The identifiability of the MCR1 can be shown through the following theorems.

Theorem 7.1. Let β 1 , β 2 , … , β k , β 1 ′ , β 2 ′ , … , β k ′ ∈ π and β 1 ≠ β 2 ≠ … ≠ β k . Then,

arg Σ j = 1 k p j x + β j 1 + β j ¯ x = arg Σ j = 1 k p j ′ x + β j ′ 1 + β j ′ ¯ x for all x ∈ Ω ,

( β 1 , β 2 , … , β k ; p 1 , p 2 , … , p k ) = ( β 1 ′ , β 2 ′ , … , β k ′ ; p 1 ′ , p 2 ′ , … , p k ′ ) ,

Proof: Note that

arg Σ j = 1 k p j x + β j 1 + β j ¯ x = arg Σ j = 1 k p j ′ x + β j ′ 1 + β j ′ ¯ x for all x ∈ Ω

Σ j = 1 k p j x + β j 1 + β j ¯ x Σ j = 1 k p j ′ x + β j ′ 1 + β j ′ ¯ x = t for all x ∈ Ω ,

Σ j = 1 k p j cos 2 θ j Σ j = 1 k p j ′ cos 2 θ j ′ = Σ j = 1 k p j sin 2 θ j Σ j = 1 k p j ′ sin 2 θ j ′ ,

β j = r j exp ( i Φ j ) , β j ′ = r j ′ exp ( i Φ j ′ ) , x = exp ( i α ) , θ j = tan − 1 r j sin ( Φ j − α ) 1 + r j cos ( Φ j − α ) , θ j ′ = tan − 1 r j ′ sin ( Φ j ′ − α ) 1 + r j ′ cos ( Φ j ′ − α ) .

Σ j = 1 k Σ l = 1 k p j p l ′ sin 2 ( θ l ′ − θ j ) = 0 ,

(7.1)

z − z ¯ = 0 ,

z = Σ j = 1 k Σ l = 1 k p j p l ′ exp 2 i ( θ l ′ − θ j ) = Σ j = 1 k Σ l = 1 k p j p l ′ exp ( 2 i θ l ′ ) exp ( 2 i θ j ) .

exp ( 2 i θ j ) = 1 + β j x ¯ 1 + β j ¯ x ,

Hence, the given statement is true only when

( β 1 , β 2 , … , β k ; p 1 , p 2 , … , p k ) = ( β 1 ′ , β 2 ′ , … , β k ′ ; p 1 ′ , p 2 ′ , … , p k ′ ) .

Theorem 7.2. Let β 1 , β 2 , … , β k , β 1 ′ , β 2 ′ , … , β k ′ ∈ π , and β 0 , β 0 ′ ∈ Ω . Then,

β 0 exp i arg Σ j = 1 k p j x + β j 1 + β j ¯ x = β 0 ′ exp i arg Σ j = 1 k p j ′ x + β j ′ 1 + β j ¯ x ′ for all x ∈ Ω

Proof: Let, β 0 / β 0 ′ = exp ( i γ ) . Then, rotating the axes to an angle γ anticlockwise, it can be seen from the previous theorem that the result will hold uniquely. This is because then the new β j 's will be nothing but β j exp ( − i γ ) .

Note: Theorem 7.1 is enough to say that the model is identifiable in case of MCR with more than one regressors because if there are no two models corresponding to ( x , x , … , x ) ∈ Ω k then, there cannot be two different models for all ( x 1 , … , x k ) ∈ Ω k .

The identifiability of the MCR2 model when | β | ≠ 1 can be proved through the following theorem.

Theorem 7.3. Let, A 1 , A 2 , … , A k , A 1 ′ , A 2 ′ , … , A k ′ ∈ Ω , β 1 ( 1 ) , β 1 ( 2 ) ∈ π . Then

x * ( 1 ) + β 1 ( 1 ) 1 + β 1 ( 1 ) ¯ x * ( 1 ) = x * ( 2 ) + β 1 ( 2 ) 1 + β 1 ( 2 ) ¯ x * ( 2 ) for all x ∈ Ω ,

(7.2)

implies ( A 1 , A 2 , … , A k , β 1 ( 1 ) , p 1 , … , p k ) = ( A 1 ′ , A 2 ′ , … , A k ′ , β 1 ( 2 ) , p 1 ′ , … , p k ′ ) where x * ( m ) = Σ j = 1 k p j ( m ) A j ( m ) x j | Σ j = 1 k p j ( m ) A j ( m ) x j | for m = 1 , 2 and A 1 = A 1 ′ = 1 . Then, if we take x j = A j − 1 x for each j, we have x * ( 1 ) = x and x * ( 2 ) = x Σ j = 1 k p j ′ A j ′ A j − 1 | Σ j = 1 k p j ′ A j ′ A j − 1 | . Denoting A = Σ j = 1 k p j ′ A j ′ A j − 1 | Σ j = 1 k p j ′ A j ′ A j − 1 | , Equation (7.3) implies