Two-Level Factorial Design with Circular Response: Model and Analysis

Since late thirties, factorial analysis of a response measured on the real line has been well established and documented in the literature. No such analysis, however, is available for a response measured on the circle (or sphere in general), despite the fact that many designed experiments in industry, medicine, psychology and biology could result in an angular response. In this paper a full factorial analysis is presented for a circular response using the Spherical Projected Multivariate Linear model. Main and interaction effects are defined, estimated and tested. Analogy to the linear response case, two new effect plots: Circular-Main Effect and CircularInteraction Effect plots are proposed to visualize main and interaction effects on circular responses.


Introduction
Factorial designs are widely used in experiments involving several factors where it is necessary to investigate the joint effects of the factors on a response variable. These joint effects include either the sole effect of each factor (main) or any interaction between two or more factors. The analysis of factorial designs is well established for a response variable that is measured on the real line. Ideas of this analysis were developed in the late thirties by Fisher in 1935and Yates in 1937[Hinkelmann and Kempthorne, 1994]. Since then this topic was a standard chapter in every Design of Experiments textbook; among these are: Fisher (1960), Kuehl (1994), Hinkelmann and Kempthorne (1994), , Montgomery (2009).
Little attention has been paid, however, to the analysis of factorial designs when the response is measured on a circle (or sphere in general). Two approaches to circular data analysis are generally used in the literature: the intrinsic and the embedding approach. In intrinsic approach, circular data are analyzed in its angular form (polar coordinates with a unit radius), in which angles are considered random variables on the unit circle that follow some directional probability distributions. Most commonly, the von Mises distribution is assumed which plays a key role in statistical inference similar to the role of the normal distribution on the real line. Underwood and Chapman (1985) developed a multi-way analysis of a circular response for crossed designs as an extension of Watson-Williams likelihood ratio test that had been initially developed for the one-way analysis. However, a simulation study of Anderson and Wu (1995) showed that this extension does not suit factorial designs because of the difficulty of identifying what the interaction terms measure and the possibility of getting negative sum of squares for these terms. In the embedding approach, the angles are represented by unit vectors (Cartesian representation). Harrison, Kanji and Gadsden (1986) and Harrison and Kanji (1988) used this approach to develop an ANOVA type approach to circular data. Considering this geometric aspect when decomposing the variance yields sums of squared of norms that would never be negative. This decomposition works only for problems where location and dispersion effects are studied simultaneously [Anderson and Wu, 1995]. Anderson and Wu (1995) proposed a likelihood ratio test statistic for testing the effects from a factorial design on the location of a circular response. They used the angular rotations occurring when moving from one level of one factor to its other level at fixed level of the other factor to evaluate the two-factor interaction effect on the location. If these angular rotations are not the same in magnitude or direction, then there is evidence of interaction. For higher order interaction, a two level surrogate factor is defined such that all combinations of odd high levels are considered as one level and all the combinations of even high levels constitute the second level. If the angular rotation resulting when moving from the high level of this surrogate factor to its low level is zero, then there is no interaction. This zero angular rotation concept is also used to evaluate any main effect. However no effect estimation was given in their paper, while merely relative importance of the significant effects was used as an approach to rank effects. Anderson and Wu (1995) justified the lack of model by intractability reasons.
In the context of regression analysis of circular response on real-line valued independent variables, remarkable link functions were offered to map the real line onto the unit circle (see Gould, 1969, Johnson and Wehrly, 1978, and Fisher and Lee, 1992. However, all of these models suffer from computational difficulties because of either the existing of a multimodal likelihood or an irregular likelihood where the maximum might occur on a very narrow peak [Presnell et al., 1998]. In addition, for factorial design, main and interaction effects are hard to be identified and separated from each other [Anderson and Wu, 1995]. Presnell et al. (1998) emphasized that the embedding approach is more useful 2 k FD with Circular Response 417 in modeling directional data on real-line valued variables. They proposed the Spherically Projected Multivariate Linear (SPML) model for circular response. The EM algorithm or Newton-Raphson algorithm could be used in the estimation phase. Presnell et al. noted that Newton-Raphson did converge faster than the EM algorithm but the latter is easier to program. In addition they indicated that SPML model could work well for designed experiments, but did not show any analysis in that context. In brief, up to date no workable model is known for the analysis of factorial designs with circular response. This results in the lack of complete analysis for factorial designs for circular responses. In this paper, the SMPL model is used to analyze factorial designs with circular response. Utilizing the SPML model, enables the author to offer a complete factorial analysis (estimation, testing and interpretation of the effects) for a circular response in a factorial setup. In addition and analogy to the effect-plots for responses measured on the real line, two new effect-plots are introduced to complete the factorial analysis for circular responses. Definitions of main and interaction effects under SPML model are introduced in Section 3. In addition, two proposed effect plots: Circular-Main Effect (C-ME) and Circular-Interaction Effect (C-IE) plots are introduced. Effect estimation and testing are presented in Section 4. An illustrative example is given in Section 5. Finally, Section 6 concludes the paper. The following section gives the notations of circular data and briefly reviews the SPML model.

Notations and the SPML Model for Circular Data
Let θ i , i = 1, · · · , n, be a random sample of circular data. Then u i = ([c i , s i ]) , is a unit random vector pointing in the direction of the angle θ i with a mean direction unit vector of E(u i )/ρ, where c i = cos(θ i ) and s i = sin(θ i ), i = 1, · · · , n, ρ = ||E(u i )|| and || · || represents the length of the vector. The mean resultant vector u = [ i c i /n, i s i /n] is pointing in the direction of the overall mean direction (average direction of the sample) θ while its length, R = ||u||, is used as a concentration measure. The circular variance is (1 − R). When data points are more concentrated around the mean, the circular variance is close to zero. A close to one value is associated with data that are evenly spread around the circle.
The idea of the embedding approach is that distributions on S p−1 can be obtained by radial projection of distributions on the R p . If y is a vector in R 2 , where Pr(y = 0) = 0, then the distribution of y/||y|| is a random point on the unit circle. Under the assumption that y is bivariate normal (N 2 ) random vector with mean vector µ and variance covariance matrix Σ, Mardia (1972) proved that the unit vector y/||y|| is distributed as projected normal (PN 2 ) with mean vector µ and variance covariance matrix Σ [Mardia and Jupp, 1999, p. 46]. Assuming u i = y i /||y i ||, Presnell et al. (1998) developed the SPML model, where u i ∼ PN 2 (µ i , Σ), µ i is the mean direction vector that is pointing in the direction of the angle i and its length serves as the concentration parameter. Further, SPML model assumes that the mean direction vector is linked linearly to k covariates as follows: x ik ) and each β is of order (k + 1). Note that the dispersion depends on the covariates since it is a function of the concentration parameter η i = ||µ i ||. The SPML model parameters are only identifiable under |Σ| = 1. Presnell et al. (1998) used the identity matrix as the variance covariance matrix and showed that PN(µ i , I) is a good approximation of the von Mises distribution. Using the EM algorithm, the SPML parameters could be estimated. Presnell et al. (1998) indicated that the predicted angle is approximately normal distributed with mean i and approximate variance is the observed Fisher information matrix associated with SPML model. The log likelihood function of the SPML model is where ψ(t) = log[1 + (tΦ(t)/φ(t))], φ(t) and Φ(t) are the standard normal density function and the standard normal distribution function, respectively. Presnell et al. (1998) proposed using the asymptotic likelihood ratio test to compare nested models; and the Akaike Information and/or Bayes Information Criterion (AIC and/or BIC) to select the best model among non-nested models.

Factorial Designs, SPML Model and SPML-Effect Model
The radial projection of the circular data facilitates the use of multivariate regression, where the response variables are the Cartesian coordinates of the angles. Using SPML model in the context of factorial design enables us to define, estimate and test the main and interaction effects in a unified manner.
Consider a 2 2 factorial design with interaction, SPML model is where u ijk is the unit vector associated with angle θ ijk , which is the k th replicate that is measured at the i th level of factor A, and the j th level of factor B, i = 1, 2, j = 1, 2, and k = 1, · · · , r (r the number of replicates). The x's are coded variables that are set at −1 if the corresponding factor is at its low (1) level and 1 if it is at its high (2) level. ε ijk ∼ PN(0, I) is an error row vector of order 2 associated with θ ijk while β's are the suitable intercept and slope parameters associated with factor A and B in the sine and cosine equations (Cartesian representation). The SPML model decomposes the main (or interaction) effect into two sub-effects associated with the Cartesian representation. The SPML slope parameters can be multiplied by two to be analogous to an effect of a two-unit change in a factor changing from −1 to 1. The main/interaction effect is then the combined effect of its sub-effects that will affect the direction of the mean response vector and/or its length.
Figure 1(a) shows high and low fitted mean response vectors of factor A, u 2.. and u 1·· and a possible effect of factor A. The main effect is then the change that occurs in the direction and/or the length of the fitted mean response vector when changing factor A from high to low averaged over all levels of factor B and all its replicates. This main effect is a joint location-dispersion effect, where the location effect is the angular difference between the two angels θ A 2 and θ A 1 associated with u 2.. and u 1.. , while the dispersion effect is the difference in the length of these two vectors. Differences that occur in the cosine and sine coordinates when changing from high to low construct two important sine and cosine coordinates which define the resultant main-effect vector [A] = u 2.. − u 1.. . Figure 1(b) shows the proposed Circular-Main Effect (C-ME) plot. The graph depends on graphing the two mean vectors with the low-mean vector translated to begin from the endpoint of the high-mean vector and mirrored/reflected (graphed in the opposite direction of its normal direction) to reflect the fact that it is multiplied by a negative sign. There will be no main effect if [A] = 0, i.e., both coordinates of the resultant main effect vector are zeros. In Figure 1(b) a nonzero A main effect is present. To have a zero A main effect both its slopes in SPML has to be zero (see Appendix 1). In Appendix 2, possible main effect scenarios are shown.
Regarding the interaction effect, it is hard to visualize the SPML slopes of sub-interaction effects in a graph like Figure 1(a), but it could be visualized in a plot similar to C-ME plot. Figure 2 depicts the Circular-Interaction Effect (C-IE) plot that shows the four fitted mean response vectors that are involved in AB interaction effect [for more details see Appendix 3]. In the proposed C-IE plot no interaction occurs if the resultant vector of the four vectors is the zero-vector. Figure 2 shows a nonzero interaction.
Analogy to the factorial design in the linear response, the simple effect of a factor is the mean response change that occurs when moving from the high level   Hinkelmann (2004) in the linear response case, the sub-interaction is the "average difference" between the simple effects of one factor at levels of other factor (or combinations of other factors) in the experiment. Higher order interactions are definable also in terms of simple effects of the lower-factor interactions (see example in Appendix 3). Analogy to the linear response case, the SPML-effect model could be defined as where τ A , τ B and τ AB represent the effects of A, B and their interaction, respec-tively, at the different levels i = 1, 2, j = 1, 2, whereas µ · is the grand mean. τ A (τ B ) is the deviation of µ i· (µ ·j ) from the grand mean, whereas τ AB is the difference of the deviation of µ ij mean from the grand mean and the sum of τ Ai and τ Bj at suitable i and j levels.

Estimation and Testing Factor Effects in SPML Model
Estimation and testing of the SPML parameters is documented in Presnell et al. (1998). In the factorial design case and as a result of the orthogonality of the X columns, we have X X = 2 f rI p , where f is the number of factors in the experiment and p = k + 1 is the number of parameters in one equation. This simplifies the SPML likelihood estimation equations as following:

Φ(t) is the standard normal distribution form and φ(t) is the standard normal density form.
To test the significance of a main or interaction effect, the asymptotic likelihood ratio test is used. Generally in factorial designs, interaction effects are tested first, if it was significant the main effects of the factors involved in this interaction are not to be tested.

Illustrative Example
The effect of four two-level factors (A: location of a butt weld to the flywheelfixed or random, B: flywheel radius grade-low or high, C: flywheel thickness grade-high or low, D: size of the counter weight attached at 0-low or high) on the angle where a corrective weight should be attached to balance the engine part is under study. Each combination of the factors is replicated ten times. The data set is published in Anderson and Wu (1996). A MINITAB14 macro is written to fit the SPML model and to calculate the likelihood ratio test. The macro programs the EM algorithm and graphs the C-ME and C-IE plots. Figure 3 shows C-ME plot of the main effects considered in fitted full 2 4 FD model. It is clear how factor B and D affect mainly the dispersion of the 2 k FD with Circular Response 423 response, while factor C affects both location and dispersion. On the other hand, main effect of factor A is relatively much smaller than the other main effects and hence of less importance.  Table 1 contains different statistics of different models that are of interest in any initial step of a factorial analysis. The best model in terms of AIC is the model that contains up to all three factor interactions, while in terms of BIC, the best model is the model that contains only up to all two factor interactions. To compare different models in terms of the likelihood ratio principle, one can use the concept of deviance. Deviance allows one to test whether the fitted reduced model is significantly worse than the saturated model (lack of fit test). By definition the deviance is twice the negative difference between the log likelihood functions of the reduced and the saturated model, which is asymptotically chi-squared distributed with degrees of freedom equal to the difference of the number of parameters in the two models. Differences in deviances would be used to compare two nested models [Myers et al. 2002, p 114]. In terms of the deviance analysis, the model without any higher than two-factor interactions is reasonable and does not suffer from lack  Figure 4 for the fitted full model. The relatively most important interaction is the BC while the least important one is BD. The remaining four two-way interactions seem to be of same importance. Table 1 contains different statistics of different models that are of interest in any initial step of a factorial analysis. The best model in terms of AIC is the model that contains up to all three factor interactions, while in terms of BIC, the best model is the model that contains only up to all two factor interactions. To compare different models in terms of the likelihood ratio principle, one can use the concept of deviance. Deviance allows one to test whether the fitted reduced model is significantly worse than the saturated model (lack of fit test). By definition the deviance is twice the negative difference between the log likelihood functions of the reduced and the saturated model, which is asymptotically chi-squared distributed with degrees of freedom equal to the difference of the number of parameters in the two models. Differences in deviances would be used to compare two nested 424 Alyaa Roshdy Zahran 13 set should contain B, C, D, CD and BC. This model has AIC=460.337, BIC=480.012 and 2log(L) =-218.168 (p-value=0.733592). Note that this model has the least BIC statistics among all the models considered in Table 1.  [Myers et al., 2002, p. 114]. In terms of the deviance analysis, the model without any higher than two-factor interactions is reasonable and does not suffer from lack of fit (p-value = 0.66113). This model is the best model in terms of BIC, too. Although the model without all the two-factor interaction does not suffer from lack of fit (p-value = 0.06480), models without CD or BC do suffer from lack of fit (p-value = 0.005424 and 0.031635, respectively). Hence, keeping CD and BC is recommended, which leads to keep also main effect B, C, and D. The deviance analysis shows also that the model without main effect A fits reasonably (p-value = 0.436049). Models that do suffer from lack of fit (given an asterisk in Table 1) indicate that dropping the chosen effect is not recommended. Accordingly, the best SPML model for this data set should contain B, C, D, CD and BC. This model has AIC = 460.337, BIC = 480.012 and 2 log(L) = −218.168 (p-value = 0.733592). Note that this model has the least BIC statistics among all the models considered in Table 1. Wu (1995, 1996) considered the same data set to study a location-only and a dispersion-only effect, respectively. Their final locationinfluential factors are A, C, D, CD and AD, while their dispersion model contained B, C, D, and CD. Our SPML model contains their dispersion-factors,

Conclusion
The SPML model is used in the factorial design setup to help in modeling circular data with predictors measured on the real line in a relatively easy way than the models proposed in the circular regression setup. SPML solved the dilemma of not having a suitable model to complete a factorial analysis for factorial data with a response measured on the circle. It enables us to determine and visualize the factor effects on the response. The proposed C-ME and C-IE plots with the circular response plays the same role the effect plots with a linear response is playing in the factorial analysis. The deviance analysis is used to identify adequate model with significant effects.
The use of SPML model in the factorial design setup along with the two effectplots and the effect interpretation, proposed in this paper, enables the researcher to perform a complete analysis of a factorial experiment for a circular response. For the cosine equation, we have This holds if one of the following three solutions satisfied: (1) β cos 0 = 0 and β cos A = 0 (2) ||u 2 || = ||u 1 || and β cos A = 0 (3) ||u 2 || = ||u 1 || = 0.
The first and the third solutions are not acceptable (trivial solution), since it means SPML model is not significant, while the second solution means that the dispersion is fixed and sub-effects are zeros. Note that fixed dispersion is satisfied if Appendix 2 Figure A2.1 represents two scenarios where main effect affects both dispersion and location while only one of [A] sine or cosine coordinates equals zero. Figure A2.2 shows possible scenarios in which the main effect increases the dispersion without any effect on the location while both slopes are nonzero. Possible Scenarios in which the main effect increases the location without any effect on dispersion are depicted in Figure A2.3. Note that in all of these three figures if we switched the high mean response vector (blue vector) with the low response vector (red vector), the effect on location or/and dispersion is still present but with the opposite sign. Scenario d in Figure A2.3 shows that when the high and low vectors of the factor are 2π apart from each other, [A] will exactly lie on the high vector and with same direction and dispersion.

Appendix 3
The SPML model of the full 2 4 FD for the example data Sub-effects for cosine equation: [A] cos = {c 2.... − c 1.... },  Other main effects, 2-way interactions and 3-way interactions are defined like the above ones with its suitable subscript. In addition, the definitions of the sine sub-effects are similar to those of the cosine counterparts where the sine function is used instead of the cosine one.