Summarizing circular data

Data collected from a circular parametric space cannot be summarized by standard statistics. As example, the arithmetic mean must be discarded because of its misleading behavior (i.e. 10° and 350° have and arithmetic mean of 180°, whereas 0° would be a more realistic estimate), and the standard definition of variance cannot be applied either. The circular mean ( $\overline{{y_{i} }} $ ) can be calculated as Fisher (Reference Fisher1993):

$$\overline{{y_{i} }} {\equals}{\rm atan}\,{\rm 2}\left( {{1 \over n}\mathop \sum\limits_{i{\rm {\equals}1}}^n {\rm sin}\left( {y_{i} } \right),{1 \over n}\mathop \sum\limits_{i{\rm {\equals}1}}^n {\rm cos}\left( {y_{i} } \right)} \right)$$

where atan2( , ) is the multi-valued inverse tangent with arguments as defined between parenthesis, and sin() and cos() are the sine and cosine of the angle, respectively. In a similar way, Fisher (Reference Fisher1993) defined a variance-like statistic for circular data ( $\sigma _{{y_{i} }} ^{2} $ ):

$$\sigma _{{{\rm y}_{i} }} ^{2} {\rm {\equals}1{\minus}}\sqrt {\left( {\mathop \sum\limits_{i{\rm {\equals}1}}^n {\rm sin}\left( {y_{i} } \right)} \right)^{2} {\plus}\left( {\mathop \sum\limits_{i{\rm {\equals}1}}^n {\rm cos}\left( {y_{i} } \right)} \right)^{2} } $$

with values between 0 and 1.

Naive circular model

Given a set of n quantitative data (y) taken around a circular parametric space from 0 to 2π, a Gaussian-like pattern can be anticipated by advocating the central limit theorem. Within this context, both location (μ) and dispersion (κ) parameters of y can be estimated within a Bayesian context:

$$p\left( {\mu {\rm ,}\kappa{\rm \,\mid\,}{\bf{y}} } \right)\propto p\left( {{\bf{y}} {\rm \,\mid\,}\mu {\rm ,}\kappa} \right)p\left( \mu \right)p\left(\kappa \right)$$

The von Mises distribution (Figure 1) must be viewed as a reasonable assumption for the a priori distribution of y ( $p\left( {\rm \bf y{\rm \,\mid\,\mu ,\kappa }} \right)$ ) because (1) it wraps around a circle, (2) is a close approximation to the circular normal distribution (i.e. it results from wrapping the normal distribution around a unit circle) with the same number of parameters (Breitenberger, Reference Breitenberger1963) and (3) has a more tractable expression, mathematically (Fisher, Reference Fisher1993). Assuming y _i as the ith circular record, the likelihood of circular data becomes:

$$p\left( {{\bf{y}}\hskip-2pt\mid\hskip-2pt\mu ,\kappa} \right){\equals}\mathop \prod\limits_{i{\rm {\equals}}1}^n {{e^{{\kappa\,{\rm cos}\left( {y_{i} {\minus}\mu } \right)}} } \over {2\pi \,{\rm I}_{0} \left( \kappa \right)}}$$

where the distribution reduces to a wrapped uniform when κ=0 (Figure 1). It must be noted that κ is a measure of concentration analogous to 1/σ ² (σ ² being the variance parameter of the normal distribution), π refers to the ratio of the circumference of a circle to its diameter, and I₀(κ) is the modified Bessel function of order 0, that is,

$$I_{0} \left( \kappa \right){\equals}\mathop \sum\limits_{i{\rm {\equals}}0}^\infty {{\left( {0.{\rm 25}\,\kappa^{2} } \right)^{i} } \over {\left( {{{\it i}\,!\,}} \right)^{2} }}$$

Figure 1

von Mises probability density functions with mean equal to π and dispersion parameter (κ) equal to 0.1 (fair grey line), 1 (dark grey line) and 5 (black line). Plots were presented on a circular (a) and linear (b) parametric space.

Given our current lack of knowledge about circular data, a priori distributions for μ and κ can be stated as uniform between appropriate bounds:

$$p\left( \mu \right){\equals}{1 \over {{\rm 2}\pi }}({\rm if}\,\mu {\repsilon}[0,2\pi {\rm ]}\,{\rm and}\,0\,{\rm otherwise})$$

$$p\left( \kappa \right){\equals}{1 \over Y}({\rm if}\,\kappa{\repsilon}{\rm [0,}Y{\rm ]}\,{\rm and}\,0\,{\rm otherwise})$$

where γ is a large enough boundary for κ.

Uncertainly about location parameter

Previous parametrization can generalize to a more flexible hierarchical structure, the CM model. This replaces μ by a set of systematic (b) and random environmental factors (p _j ) linked by appropriate incidence matrices (X and Z _j ) within an additive structure:

$${\bf{y}} {\equals}{\bf{Xb}} {\plus}\mathop \sum\limits_{j{\rm {\equals}1}}^r {\bf{Z}} _{j} {\bf{p}} _{j} {\plus}{\bf{e}} $$

Note that e refers to the circular residual term, and flexibility is gained by assuming r independent and normally distributed sources of environmental variation with mean 0 and variance $\sigma _{{{\bf{p}} _{j} }}^{2} $ . The joint posterior distribution of all unknown parameters in the model can be obtained from

$$p\left( {\bf{b}} ,{\bf{p} ,{\bsigma } {\rm ,}\kappa{\rm \,\mid\,}{\bf{y}} } \right)\propto p\left( {{\bf{y}} \,\mid\,{\bf{b}} ,{\bf{p}} {\rm ,}\kappa} \right)p\left( {\bf{b}} \right)p\left( {\bf p\,\mid\,\bsigma } \right)p\left( {\bf \sigma } \right)p\left( \kappa \right)$$

where ${\bf p}{\equals}\left[ \bf {p_{1}^{,} p_{2}^{,} \,\ldots\, p_{r}^{,} } \right]^{,} $ and ${\bsigma} {\equals}\left[ {\sigma _{{p_{1} }}^{2} \sigma _{{p_{2} }}^{2} \,\ldots\, \sigma _{{p_{r} }}^{2} } \right]$ . As for naive circular (NC) model, the likelihood of data assumes a von Mises process with systematic and permanent environmental sources of variation:

$${\it p}\left( {{\bf y}\,\mid\,{\bf{b}} ,{\bf{p}} {\rm ,\kappa}} \right){\equals}\mathop \prod\limits_{{\it i{\equals}}\rm 1}^{\it n} {{e^{{{\rm \kappa}{\rm .cos}\left( {{\it y}_{{\it i}} {\minus}\left( {{\bf x}_{{\it i}} {\bf b}{\plus}\sum_{{\it j{\equals}\rm 1}}^{\it r} \nf z_{{{\it i,j}}} {\bf p}_{{\it j}} } \right)} \right)}} } \over {2{\rm \pi I}_{0} \left( {\rm \kappa} \right)}}$$

x _i and z _i,j being appropriate row vectors of incidence inherent to the ith circular record. Within this context, p must be viewed as a combination of multivariate normal (MVN) processes:

$${\it p}\left( {{\bf{p}} \,\mid\,{\bf{\bsigma }} } \right){\equals}\mathop \prod\limits_{{\it j{\equals}\rm 1}}^{\it r} {\rm MVN}\left( {\bf 0_{{\it j}} ,{\bf I}_{{\it j}} \sigma _{{{\rm p}_{{\it j}} }}^{2} } \right)$$

where s _j is the number of elements in p _j , 0 _j is a s _j × 1 vector of zeros, and I _j is a s _j × s _j identity matrix. Given our current lack of knowledge about systematic effects on the circular parametric space, a priori distributions of b are conceptually defined flat between 0 and 2π:

$${\it p}\left( {\bf{b}} \right){\equals}{1 \over {{\rm 2\pi }}}({\rm if}\,{\bf{b}} \in\left[ {{\rm 0,2\pi }} \right]{\rm and}\,{\rm 0\hskip3pt otherwise})$$

although an additional restriction applies to guarantee that $0\leq {\bf{x}} _{{\it i}} {\bf{b}} {\plus}\mathop \sum\limits_{{\it j{\equals}\rm 1}}^{\it r} {\bf{z}} _{{{\it i,j}}} {\bf{p}} _{{\it j}} \leq {\rm 2\pi }$ . As in previous model:

$${\it p}\left( {\kappa} \right){\equals}{1 \over {\rm \gamma }}({\rm if}\,{\rm \kappa{\repsilon}}\,{\rm [0,\it Y]\,\rm and}\,0\,{\rm otherwise})$$

whereas ${\rm p}\left( {\bf{\sigma }} \right)$ are scaled inverse χ ² distributions with hyper-parameters ν and S ² (they become flat when ν=−2 and S ²=0).

Uncertainly about dispersion parameter

Parameter κ in NC and CM models assumes a homogeneous dispersion pattern across all circular data. This assumption may not be true if residual heterogeneity holds in the data set (SanCristobal-Gaudy et al., Reference SanCristobal-Gaudy, Elsen, Bodin and Chevalet1998). Canalization plays a key role in current animal and plant breeding programs to lead more uniform products (Mulder et al., Reference Mulder, Bijma and Hill2008; Blasco et al., Reference Blasco, Martínez-Álvaro, García, Ibáñez-Escriche and Argente2017). Selection for low seasonality in sheep must be viewed as the opposite situation where heterogeneous (i.e. wider) lambing seasons inevitably impacts on other variables such as trade, consumption and lamb prices (Chemineau et al., Reference Chemineau, Malpaux, Brillard and Fostier2007). Generalizing Ibáñez-Escriche et al. (Reference Ibáñez-Escriche, Sorensen, Waagepetersen and Blasco2008), κ can be rewritten as the following hierarchical structure:

$${\rm \kappa{\equals}}e^{{{\bf{x}} _{i} {{\bbeta }} {\plus}\mathop \sum\limits_{{\it j{\equals}\rm 1}}^r {\bf{z}} _{{{\it i,j}}} {\bf{\brho }} _{{\it j}} }} $$

where β and ρ _j have the structure of effects than b and p _j , respectively. This expands the Bayesian development to accommodate a circular canalization (CC) model:

$${\it p}\left( {{\bf{b}} ,{\bf{\bbeta }} ,{\bf{p}} ,{\bf{\brho }} ,{\bf{P}} \,\mid\,{\bf{y}} } \right)\propto{\it p}\left( {{\bf{y}} \,\mid\,{\bf{b}} ,{\bf{\bbeta }} ,{\bf{p}} ,{\bf{\brho }} } \right){\it p}\left( {{\bf{b}} ,{\bf{\bbeta }} } \right){\it p}\left( {{\bf{p}} ,{\bf{\brho }} \,\mid\,{\bf{P}} } \right){\it p}\left( {\bf{P}} \right)$$

where P stores permanent environmental covariances between p and ρ terms (see below). Similarly, the a priori distribution of data includes a hierarchic dispersion term:

$${\it p}\left( {{\bf{y}} \,\mid\,{\bf{b}} ,{\bf{\bbeta }} ,{\bf{p}} ,{\bf{\brho }} } \right){\equals}\mathop \prod\limits_{{\it i{\equals}}\rm 1}^{\it n} {{{\it e}^{{{\it e}^{{\left( {{\bf x}_{{\it i}} {\bf{\bbeta }} {\plus}\sum_{{\it j{\equals}\rm 1}}^r {\bf{z}} _{{{\it i,j}}} {\bf{\brho }} _{{\it j}} } \right)}} {\rm cos}\left( {{\it y}_{{\it i}} {\minus}\left( {{\bf{x}} _{{\it i}} {\bf{b}} {\plus}\sum_{{\it j{\equals}\rm 1}}^{\it r} {\bf{z}} _{{{\it i,j}}} \rm p_{{\it j}} } \right)} \right)}} } \over {2{\rm \pi I}_{0} \left( {e^{{\left( {{\bf{x}} _{{\it i}} {\bf{\bbeta }} {\plus}\sum_{{\it j{\equals}\rm 1}}^{\it r} {\bf{z}} _{{{\it i,j}}} {\bf{\brho }} _{{\it j}} } \right)}} } \right)}}$$

and the joint distribution of p _j and ρ _j reduces to a MVN process:

$${\it p}\left( {{\bf{p}} {\rm ,}{\bf{\brho }} \,\mid\,{\bf{P}} } \right){\equals}\mathop \prod\limits_{{\it j{\equals}\rm 1}}^{\it r} {\rm p}\left( {{\bf{p}} _{{\it j}} ,{\bf{\brho }} _{{\it j}} \,\mid\,{\bf{I}} _{{\it j}} \,\otimes\,{\bf{P}} _{{\it j}} } \right){\equals}\mathop \prod\limits_{{\it j{\equals}\rm 1}}^{\it r} {\rm MVN}\left( {\matrix{ {{\bf{0}} _{{\it j}} } \cr {{\bf{0}} _{{\it j}} } \cr } ,{\rm I}_{{\it j}} \,\otimes\,{\bf{P}} _{{\it j}} } \right)$$

where $\,\otimes\,$ is the Kronecker product, and P _j is the matrix storing permanent environmental variances ( ${\rm \sigma }_{{{\rm p}_{{\it j}} }}^{2} $ and $${\rm \sigma }_{{{\rm \rho }_{{_{{\it j}} }} }}^{2} $$ ) and covariances ( $\sigma _{{p_{j} \rho _{j} }} {\equals}\sigma _{{\rho _{j} p_{j} }} $)

$${\bf P}_{j} {\equals}\left[ {\matrix{ {\sigma _{{p_{j} }} ^{2} } & {\sigma _{{p_{\it j} \rrho _{\it j} }} } \cr {\sigma _{{p_{j} \rho _{j} }} } & {\sigma _{{\rrho _{j} }} ^{2} } \cr } } \right]$$

Scaled inverse Wishart distributions are used as prior distributions for P _j because they enable both scale parameters and covariances to be estimated from data, allow flexibility in the hierarchical model and are computationally convenient (Van Tassell and Van Vleck, Reference Van Tassell and Van Vleck1996).

Markov chain Monte Carlo algorithm

Sampling procedures from the joint posterior distribution will focus on the most complex parametrization (i.e. model CC), as they can be easily adapted to models NC and CM. The algorithm used for posterior sampling was an hybrid approach involving both Metropolis-Hastings (Metropolis et al., Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller1953) for b, β, p and ρ, and Gibbs sampling (Gelfand and Smith, Reference Gelfand and Smith1990) for P. Within this context, random samples for systematic (b and β) and permanent environmental effects (p and ρ) where updated from their marginal posterior distribution derived from

$$\eqalign{{\it p}\left( {{\bf{b}} ,{\bf{\bbeta }} ,{\bf{p}} ,{\bf{\brho }} \,\mid\,{\bf{y}} ,{\bf{P}} } \right) \cr &\hskip -80pt\propto\mathop \prod\limits_{{\it i{\equals}}\rm 1}^{\it n} {{{\it e}^{{{\it e}^{{\left( {{\bf{x}} _{{\it i}} {\bf{\bbeta }} {\plus}\sum_{{\it j{\equals}\rm 1}}^{\it r} {\bf{z}} _{{{\it i,j}}} {\bf{\brho }} _{{\it j}} } \right)}} {\rm cos}\left( {{\rm y}_{{\it i}} {\minus}\left( {{\bf{x}} _{{\it i}} {\bf{b}} {\plus}\sum_{{\it j{\equals}\rm 1}}^{\it r} {\bf{z}} _{{{\it i,j}}} {\bf{p}} _{{\it j}} } \right)} \right)}} } \over {2{\rm \pi I}_{0} \left( {{\it e}^{{\left( {{\bf{x}} _{{\it i}} {\bf{\bbeta }} {\plus}\sum_{{\it j{\equals}\rm 1}}^{\it r} {\bf{z}} _{{{\it i,j}}} {\bf{\brho }} _{j} } \right)}} } \right)}}\cr &\hskip -80pt\mathop \prod\limits_{{\it j{\equals}\rm 1}}^{\it r} {\it e}^{{{\minus}{1 \over 2}\left[{{\bf{p}} _{{\it j}} {\bf{\brho }} _{{\it j}} } \right] ’\lsquo {\bf{P}} _{{\it j}}^{{{\rm {\minus}1}}} \left[ {{\bf{p}} _{{\it j}} {\bf{\brho }} _{{\it j}} } \right]}} $$

For each parameter to be sampled (θ _k ), a uniform proposal distribution between between θ _{k,−1 −} π/4 and θ _k,−1+π/4 was assumed, where θ _k,−1 was the value obtained in the previous iteration. These boundaries were defined in preliminary analyses to guarantee an acceptance rate greater than 10% in all θ _k . Once updated for p and ρ, random samples from the marginal posterior distribution of the (co)variance matrix P can be obtained as previously reported by Van Tassell and Van Vleck (Reference Van Tassell and Van Vleck1996) in multivariate linear models.

Example on field data

Lambing distribution throughout the year can be viewed as a representative example of circular quantitative data from livestock industry. This research focused on the re-scaled day of lambing from 1 January (0) to 31 December (regular years, 2π–2π/365; leap years, 2π–2π/366). A total of 110 116 records from 21 630 Ripollesa ewes were collected from 25 purebred flocks between years 1976 and 2017. The Ripollesa is a medium-sized sheep (rams, 75 to 90 kg; ewes, 50 to 65 kg of BW) and the most abundant native ovine breed in Catalonia, the north-east region of Spain. Purebred animals are white with dark brown or black marks on both head and legs, with white semi-fine wool and closed fleece. The Ripollesa breed is a representative example of meat-type sheep management under semi-extensive Mediterranean conditions in northeast Spain, the ‘pascual’-type lamb being its typical commercial product (22 to 24 kg BW at slaughter; Esquivelzeta et al., Reference Esquivelzeta, Fina, Bach, Madruga, Caja, Casellas and Piedrafita2011).

Except for NC parametrization, the analytical model included two systematic and two permanent environmental effects. More specifically, b accounted for ewe age at lambing (<3, 3, 4, 5, 6 and >6 years) and lambing interval (first lambing, <210 days, between 210 and 239 days, between 240 and 329 days, and ⩾330 days); those thresholds balanced the number of records per level, and they had been previously used in this breed for both research and genetic evaluation purposes (Casellas and Caja, Reference Casellas and Caja2014). Permanent environmental effects accounted for flock-year-season with 763 levels, and the environmental contribution of each ewe with 21 630 levels. Although the estimation of the additive genetic merit is the ultimate goal (see companion paper), genetic effects have not been considered in this part. We focus on the fit of the different models and genetic effects should be captured by to the permanent environmental effect inherent to each ewe, if any. Note that flocks were kept under continuous mating systems with multiple sires and paternity testing was not implemented in this breed until 2011; that originated a high percentage of unknown sires in the pedigree file (~95%). Incompleteness of genealogical information discouraged the implementation of animal models under complex genetic covariance structures (Id-Lahoucine and Casellas, Reference Id-Lahoucine and Casellas2017) such as the CC model. If genetic effects are included, fewer ewes with complete pedigree information and close relatives have to be used in order to constrain the estimation problem to a reasonable paradigm. In such a reduced data set, inevitably, the adequacy of each model would be difficult to assess.

All the unknown parameters in models were sampled from their marginal posterior distributions by appropriate Gibbs sampling (Gelfand and Smith, Reference Gelfand and Smith1990) or Metropolis-Hastings (Metropolis et al., Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller1953). Three independent Markov chain Monte Carlo (MCMC) with 510 000 iterations were run for each analysis, and the first 10000 iterations were discarded as burn-in. The suitability of the length of the burn-in period was evaluated by both visual inspection and the method of Raftery and Lewis (Reference Raftery and Lewis1992) on the sampling path of κ (NC model) and $\sigma _{{{\rm p}_{{\it j}} }}^{2} $ (CM and CC models). A total of 10 000 samples of model parameters were saved from each chain with a lag interval of 50 iterations. All analyses were performed in a server with two Intel^® Xeon^® E5-2400 (Foxen SA, Madrid, Spain) processors with a CentOS 6.6 Linux operating system. In terms of central processing unit time, the most complex model (CC) spent almost 35 days for each independent 510 000-iteration MCMC.

Model comparison

Model performance was evaluated on two different ways, goodness of fit and predictive ability. Moreover, goodness of fit was assessed twice, by the deviance information criterion (DIC; Spiegelhalter et al., Reference Spiegelhalter, Best, Carlin and van der Linder2002) and the prediction of angular data in terms of mean square error (MSE) between phenotypic and predicted records. The expectation of the predictive distribution of a given ith record ( $\tilde{y}_{i} $ ) under NC model reduced to μ, whereas this expanded to ${\bf x}_{{\it i}} {\tilde{\bf b}}_{{\it i}} {\plus}\mathop \sum\limits_{{\it j{\equals}\rm 1}}^{\it r} {\bf z}_{{{\it I,j}}} {\tilde{\bf p}}_{{\it I}} $ under CM and CC models. Note that ${\tilde{\bf b}}_{{\it I}} $ and ${\tilde{\bf p}}_{{\it I}} $ were posterior modal estimates of b and p _j . The MSE was calculated in each iteration of the MCMC process as

$${\rm MSE}{\equals}{{\sum_{{\it i{\equals}\rm 1}}^{\it n} \left( {{\it y}_{{\it i}} {\minus}{\tilde{y}}_{{\it i}} } \right)^{2} } \over {\it n}}$$

On the other hand, predictive ability relied on the prediction of future observations given a trained data sample (Matos et al., Reference Matos, Thomas, Gianola, Pérez-Enciso and Young1997). Two new data sets were created by randomly keeping 50% of records from the initial data set, and both were analyzed under all three models. The MSE was computed for each MCMC iteration between expectations from the predicted distribution and the removed records.