CircNNTSR : An R Package for the Statistical Analysis of Circular, Multivariate Circular, and Spherical Data Using Nonnegative Trigonometric Sums

The statistical analysis of circular, multivariate circular, and spherical data is very important in diﬀerent areas, such as paleomagnetism, astronomy and biology. The use of nonnegative trigonometric sums allows for the construction of ﬂexible probability models for these types of data to model datasets with skewness and multiple modes. The R package CircNNTSR includes functions to plot, ﬁt by maximum likelihood, and simulate models based on nonnegative trigonometric sums for circular, multivariate circular, and spherical data. For maximum likelihood estimation of the models for the three diﬀerent types of data an eﬃcient Newton-like algorithm on a hypersphere is used. Examples of applications of the functions provided in the CircNNTSR package to actual and simulated datasets are presented and it is shown how the package can be used to test for uniformity, homogeneity, and independence using likelihood ratio tests. the values of the maximized log-likelihoods for SNNTS models with the combinations of M 1 = 0 , 1 , 2 , 3 and M 2 = 0 , 1 , 2 , 3 and the values of the test statistic for uniformity, Λ Unif , and asymptotic p value. From these results, the null hypothesis of uniformity is rejected at the 5% signiﬁcance level for all cases except M = (1 , 0) and M = (2 , 0), which correspond to radially symmetric distributions; although these models gave a very poor ﬁt relative to others, except M = (0 , 0).


Introduction
There exist random variables that have supports that differ from a subset of the real line. For example, the direction taken by an animal and the wind direction correspond to random variables with support corresponding to the circumference of the unit circle. In these cases, the random variable θ is an angle with possible values on the circumference of the unit circle in the interval (0, 2π]. Therefore, the density function of a circular random variable must be periodic. By considering a vector of circular random variables, θ = (θ 1 , θ 2 , . . . , θ R ), one can analyze multivariate circular data. Examples of multivariate circular data include the pair of dihedral angles in the backbone of a protein and the wind directions at different monitoring stations. Spherical data consist of a bivariate vector θ s = (θ s1 , θ s2 ) with θ s1 ∈ (0, 2π] and θ s2 ∈ (0, π]. In the case of data for the Earth's surface, θ s1 is identified with a longitudinal coordinate and θ s2 is identified with a latitudinal coordinate. Examples of spherical data include the position of the occurrence of different events on the Earth's surface such as the epicenter of an earthquake or the trajectory of a hurricane. The time at which an event occurs, such as the epileptic seizures of a patient during the day, can also be analyzed as circular data. Regarding possible applications, Ridout and Linkie (2009), Linkie and Ridout (2011), and Lynam et al. (2013), for example, used circular models to study the activity patterns of wild cats using camera-trap data. Fernández-Durán (2004) developed a family of univariate densities for circular data based on nonnegative trigonometric sums (Fejér 1916). This family is referred to as NNTS models. In terms of complex numbers, the probability density function of an NNTS model can be expressed as follows: where i = √ −1, e ikθ = cos(kθ) + i sin(kθ), c k = c rk + ic ik is a complex number with c rk and c ik being the real and imaginary components of c k , respectively, andc k = c rk − ic ik represents the conjugate of c k . Then, f (θ) is the squared norm of a sum of complex numbers and is thus nonnegative. In order for f (θ) to integrate to one, the following restriction on the c parameters must be imposed: implying that c i0 = 0, i.e., c 0 must be a real number. Given this restriction on the c parameters, the case of M = 0 corresponds to the uniform circular density, f (θ) = 1 2π . The number of free c parameters is equal to 2M . This idea of using the squared norm of a sum of complex numbers to define a density function can be extended to specify the joint density of a vector of circular random variables in the case of multivariate circular data and the bivariate random vector that defines spherical data. The NNTS models for multivariate circular and spherical data are referred to as MNNTS and SNNTS models, respectively. In the case of multivariate circular data, the squared norm of a multiple complex trigonometric sum is used. For a bivariate circular random vector, θ = (θ 1 , θ 2 ), let In order for f (θ 1 , θ 2 ) to integrate to one, the following restriction on the c parameters must be imposed: The support of a bivariate circular random vector is the surface of a torus. This result can be extended to circular random vectors of dimension R greater than two by defining with the following restriction on the c parameters: Given this restriction, the parameter space of NNTS and MNNTS densities correspond to the surface of a hypersphere. Fernández-Durán and Gregorio-Domínguez (2014b) demonstrated that the marginal and conditional densities of an MNNTS joint density are members of the MNNTS family. The number of free c parameters of an MNNTS model is equal to 2( R k=1 (M k + 1)) − 2. In the case of spherical data, θ s = (θ s1 , θ s2 ) with θ s1 ∈ (0, 2π] and θ s2 ∈ (0, π]. The SNNTS density is defined as follows: The term sin(θ s2 ) corresponds to the uniform measure on the sphere. The restriction in the parameter space corresponds to the following: The restriction in Equation 8 can also be written as by a transformation of the parameter space (for details see Fernández-Durán and Gregorio-Domínguez 2014a). Therefore, the new transformed space corresponds again to a hypersphere, and the NNTS, MNNTS and SNNTS cases can be developed in a similar manner. The number of free c parameters in an SNNTS model is equal to 2(M 1 + 1)(M 2 + 1) − 2.

Maximum likelihood estimation and likelihood ratio tests
Maximum likelihood estimation of the c parameters of the NNTS, MNNTS, and SNNTS models is implemented numerically by using a Newton-like algorithm on manifolds. In the NNTS, MNNTS, and transformed SNNTS cases, the c parameter space (smooth Riemann manifold) corresponds to a hypersphere as defined by the restrictions in Equations 2, 4, 6, and 9. This algorithm is described in detail in Fernández-Durán and Gregorio-Domínguez (2010b) for circular data. The modifications for multivariate and spherical data are directly obtainable and can be consulted in Fernández-Durán and Gregorio-Domínguez (2014a,b).
By using the maximum likelihood method of estimation, it is possible to apply different likelihood ratio tests. In practice, tests of uniformity, homogeneity, and independence are among the most important tests for circular, multivariate circular, and spherical data. The uniformity null hypothesis states that the data follow a uniform distribution. Given the definition of NNTS, MNNTS, and SNNTS models, it is equivalent to test that the M parameters are equal to zero. For example, in the case of circular data, the uniformity null hypothesis is equivalent to H 0 : M = 0 and a likelihood ratio test can be implemented by comparing the maximized log-likelihood of the uniform model, which is an NNTS model with M = 0, against the maximized log-likelihood of an NNTS model with M > 0. The likelihood-ratio test statistic Λ Unif (see Fernández-Durán and Gregorio-Domínguez 2014c) asymptotically follows a chi-squared distribution with twice as many degrees of freedom as the number of free c parameters in the non-uniform model. The same definition can be used for multivariate circular and spherical uniform tests. One disadvantage of the NNTS test of uniformity is the need to select a value for M which is related to the maximum number of peaks considered to exist in the population. Fernández-Durán and Gregorio-Domínguez (2014c) compare the NNTS uniformity test with other more conventional tests and recommend to select the larger value when considering two consecutive values for M . In the cases of multivariate circular and spherical data, it is relevant to test for independence among the components or subsets of components of the random vector. Again, a likelihood ratio test can be implemented by comparing the maximized log-likelihood of the considered joint model against the maximized log-likelihood of the corresponding model assuming a certain structure of stochastic independence among the components of the random vector. For example, in the case of testing for independence among the components of a random vector, i.e., under the null hypothesis that the joint density is the product of the univariate marginal densities, the likelihood ratio test statistic, Λ Ind , is defined as follows: wherel M is the maximized log-likelihood of the joint model with M = (M 1 , . . . , M R ) and l M k is the maximized log-likelihood for the k-th component of the random vector. The test statistic, Λ Ind , asymptotically follows a chi-squared distribution with degrees of freedom In the case of spherical data, one can only test for independence among the latitude and longitude components of the spherical random vector.
When comparing K different populations, a homogeneity test is typically relevant. Again, a homogeneity likelihood ratio test can be constructed by comparing the maximized loglikelihood of the homogeneous model when fitting an NNTS (MNNTS or SNNTS) model to the pooled data from the K different populations to the sum of the maximized log-likelihoods of the NNTS (MNNTS or SNNTS) models fitted to the data of each of the K different populations being compared. Interested readers can consult Fernández-Durán and Gregorio-Domínguez (2010a) for the construction of a likelihood ratio homogeneity test in the case of circular data. The test statistic, Λ Hom , asymptotically follows a chi-squared distribution. If the homogeneous model uses an

The functions in CircNNTSR
The R (R Core Team 2016) package CircNNTSR (Fernández-Durán and Gregorio-Domínguez 2016) is available from the Comprehensive R Archive Network (CRAN) at http://CRAN. R-project.org/package=CircNNTSR and includes functions to perform maximum likelihood estimations, calculate and plot the density and distribution functions, and perform simulations of NNTS, MNNTS, and SNNTS models.
The names of functions for circular data begin with nnts. For example, the function nntsdensity calculates the density of an NNTS model with given c and M parameters. Similarly, the names of functions for spherical data begin with snnts, and those for multivariate circular data begin with mnnts. The package CircNNTSR also includes functions to consider grouped or incidence circular data in which only the total number of observations in different intervals of the random variable is provided. The names of functions for incidence data in the interval (0, 2π] end with interval0to2pi, and those in the interval (0, 1], which is commonly used when considering data on the occurrence of events in time, end with interval0to1. The function that implements the Newton-like algorithm on a hypersphere to obtain the maximum likelihood estimates of the c parameters is named nntsmanifoldnewtonestimation, with an s and an m adding to the beginning of the name for the spherical and multivariate circular cases, respectively. For example, the following commands fit an NNTS circular model with M = 3 to data on the directions taken by 76 turtles after a treatment. The data were analyzed by Stephens (1969) and Fisher (1993, dataset B.3).
R> library("CircNNTSR") R> data("Turtles_radians", package = "CircNNTSR") R> cnnts <-nntsmanifoldnewtonestimation(data = Turtles_radians, M = 3, + iter = 1000) R> cnnts $cestimates k cestimates 1 0 0.28645175+0.00000000i 2 1 0.11655438-0.12669303i 3 2 -0.14659000-0.16080633i 4 3 0.01079598+0.00065866i The value of M represents the maximum number of modes in the data; the exploratory analysis of the data could give an indication of this value. Also, the degree of concentration of data around the modes should be considered. For highly concentrated data, large values of M are needed in order to get a good fit. The value of iter depends upon the complexity of the model, it should be large enough for the norm of the gradient to be near zero, i.e., within the desired error margin. The function nntsmanifoldnewtonestimation returns a list with a matrix containing the c parameter estimates (cestimates). The first column of the matrix contains the index of the c parameters. The other elements in the list are the values of the maximized log-likelihood (loglik), the Akaike's information criterion (AIC), the Bayesian information criterion (BIC), and the norm of the gradient (gradnormerror). The function nntsmanifoldnewtonestimation has two additional arguments, initialpoint and cinitial, that allow the user to begin iterations of the Newton algorithm from a different point. Generally, this option is used to check the convergence of the Newton algorithm from different initial points. The function nntsrandominitial can be used to randomly select a value for cinitial. For example, one could type R> cnnts <-nntsmanifoldnewtonestimation(data = Turtles_radians, M = 3, + iter = 2000, initialpoint = TRUE, cinitial = nntsrandominitial(M = 3)) The following commands are used to plot the fitted NNTS model and the histogram of the data that are presented in Figure 1. Note that the values of the fitted NNTS density are equal at zero and 2π satisfying the periodicity condition for the density of a circular random variable.

Histogram of Turtles_radians
The following command produces the plot of the fitted bivariate MNNTS density and the corresponding univariate marginal densities in the left plot in Figure 2. The x axis corresponds to the φ dihedral angle, and the y axis corresponds to the ψ dihedral angle. The Ramachandran plot, the two dimensional plot in Cartesian coordinates with angle φ along the horizontal axis and angle ψ along the vertical axis (Ramachandran, Ramakrishnan, and Sasisekharan 1963), is also included on the right side of Figure 2.
R> mnntsplotwithmarginals(cestimates = cmnnts$cestimates, M = c(2, 2), + theta = 45) R> plot(ProteinsAAA, xlab = "x", ylab = "y", xlim = c(0, 2 * pi), + ylim = c(0, 2 * pi), cex = 0.5) Finally, the following commands fit an SNNTS model with M = (M 1 , M 2 ) = (3, 1) for a dataset on the measurements of magnetic remanence from specimens of red beds from the Bowen Basin, Queensland (Fisher, Lewis, and Embleton 1987, dataset B.5). The dataset contains 52 pairs of angles in declination-inclination coordinates that are transformed into longitude-latitude coordinates.  When using an algorithm to implement the Newton-like algorithm on the hypersphere, it is important to test different initial points to avoid local maxima. These initial points can be selected randomly using the function nntsrandominitial. If several random initial points result in the same maximum, this increases the certainty that the algorithm did not get trapped in local maxima. In our experience, for many datasets that we have analyzed, one gets consistently the same maximum when using different random initial points. Difficult datasets in which the algorithm found different maxima were detected when running the algorithm from five to ten different random initial points. For these cases, the algorithm must be executed more times from different random initial points until one gets many replicates resulting in the same maximum. In the next section, we present examples with actual datasets using the package CircNNTSR.

Circular univariate data: NNTS
We consider data on the time of occurrence of earthquakes of intensity greater than 6.0 on the Richter scale with the epicenter occurring on the coast of the Pacific Ocean in Mexico from 1920 to 2002 (Fernández-Durán and Gregorio-Domínguez 2016, EarthquakesPacificMexico-gt6 dataset). There are a total of 241 observations. We apply a uniformity test by fitting NNTS models with M = 0, 1, . . . , 6. Table 1 includes the values of the maximized loglikelihood,l, the likelihood ratio uniformity test statistic, Λ Unif , the degrees of freedom, d.f ., and the p value for the uniformity test. From the results in Table 1, the null hypothesis of uniformity is not rejected at the significance level of 5% although for the tests with M = 2, 4, 6 the p values are close to 5%. By applying other tests of uniformity included in the circular package (Agostinelli and Lund 2013), the null hypothesis of uniformity is rejected by the tests of Kuiper and Watson but not rejected by the Rayleigh's test when applied to the EarthquakesPacificMexicogt6 dataset at a 5% significance level. Rejection of the null hypothesis of uniformity could be related to the existence of aftershocks of intensity greater than 6.0 on the Richter scale.
By considering only the 76 earthquakes of intensity greater than 7 on the Richter scale (Fernández-Durán and Gregorio-Domínguez 2016, dataset EarthquakesPacificMexicogt7), Table 2 includes the results for tests of uniformity; the p values in Table 2 are greater than those in Table 1, thus clearly not rejecting the null hypothesis of uniformity. For the case of earthquakes of intensity greater than 7.0 on the Richter scale where the existence of large aftershocks is very low, the uniformity tests of Kuiper, Rayleigh and Watson also do not reject the null hypothesis of uniformity. These results support the idea of using a homogeneous Poisson process to model the times of occurrence of high-intensity earthquakes.

Multivariate circular data: MNNTS
The pollution monitoring network of Mexico's Central Valley consists of 36 monitoring stations located in Mexico City and neighboring states that measure pollution variables, such as We fitted trivariate MNNTS models with M parameters up to M 1 = 3, M 2 = 3, and M 3 = 3 (see Table 6). The best model according to BIC is the one with M 1 = 3, M 2 = 2, and M 3 = 3.
A model that has the same orders M 1 = 3, M 2 = 2, and M 3 = 3, but assumes independence among the three components has a much larger BIC than the dependent joint MNNTS model. Table 7 provide the maximized log-likelihood values of the univariate marginal models. All of the likelihood ratio tests for different combinations of M 1 , M 2 and M 3 clearly reject the null hypothesis of independence, with p values less than 10 −6 .
As another example of multivariate circular data, we simulated 200 realizations from a bivariate circular uniform distribution with independent components and fitted MNNTS to these data with M 1 = 0, 1, . . . , 3 and M 2 = 0, 1, . . . , 3. This dataset is included in CircNNTSR with the name DataUniformBivariate200obs. As expected, the best model in terms of BIC is the one with M 1 = 0 and M 2 = 0, which corresponds to a bivariate circular uniform distribution.     greater than 0.05.
The following commands can be used to generate 500 realizations from an MNNTS model with M = (2, 5, 1) and randomly selected c parameters (cparsim): R> cparsim <-mnntsrandominitial(M = c(2, 5, 1), R = 3) R> simulations <-mnntssimulation(nsim = 500, cpar = cparsim, M = c(2, 5, 1), The simulation algorithm is based on the acceptance-rejection method of simulation as suggested by Fernández-Durán and Gregorio-Domínguez (2014b). Since the Marsenne-Twister pseudorandom number generator is the default in R, when a massive number of realizations from an (MS)NNTS model is required, the proposed acceptance-rejection algorithm can be implemented in parallel to create independent random streams in different machines by applying the dynamic creation of pseudorandom number generators proposed by Matsumoto and Nishimura (2000). : Scatterplot of measurements of magnetic remanence from 52 specimens of red beds from the Bowen Basin, Queensland (see Fisher et al. 1987, dataset B.5 Figure 4 includes the scatterplot for the 52 measurements of magnetic remanence from specimens of red beds from the Bowen Basin, Queensland (Fisher et al. 1987, dataset B.5.). Based on Figure 4, it is relevant to test for spherical uniformity using the CircNNTSR dataset DataB5FisherSpherical. Table 10 includes the values of the maximized log-likelihoods for SNNTS models with the combinations of M 1 = 0, 1, 2, 3 and M 2 = 0, 1, 2, 3 and the values of the test statistic for uniformity, Λ Unif , and asymptotic p value. From these results, the null hypothesis of uniformity is rejected at the 5% significance level for all cases except M = (1, 0) and M = (2, 0), which correspond to radially symmetric distributions; although these models gave a very poor fit relative to others, except M = (0, 0).

Summary
The R package CircNNTSR includes functions for the statistical analysis of circular, multivariate circular and spherical data using models based on nonnegative trigonometric sums.
Since the parameter space of NNTS, MNNTS and SNNTS models corresponds to the surface of a hypersphere, maximum likelihood estimation is implemented by a Newton-like algorithm on manifolds. (MS)NNTS models are useful to analyze datasets that present skewness and multiple modes but that are not highly concentrated around the modes. For extremely highly concentrated datasets one can consider more parsimonious models such as those based on, for example, the von Mises distribution and its generalizations. Comparisons in terms of AIC and BIC criteria among (MS)NNTS models and alternative models, such as those included in the circular (Agostinelli and Lund 2013) and movMF (Hornik and Grün 2014) R packages, can be implemented to assess the most adequate model for a particular dataset. A future extension to the CircNNTSR package under current development is the inclusion of R functions for the statistical analysis of (MS)NNTS models from a Bayesian perspective considering Markov chain Monte Carlo (MCMC) algorithms.