A key step of any statistical multivariate analysis concerns the choice of variables in line with the main objectives of the study. Usually, the available procedures to face this problem are restricted to a-posteriori statistical analysis, using Bayesian approaches or stepwise selection procedures.

The main objective of the present paper is to revisit a framework where the a-priori choice of variables makes sense under specific conditions and to propose a factor analysis model particularly adapted to structured quantitative big data.

We have associated our complete sample of variables to a mixture of two bipolar Watson distributions defined on the n-sphere,$W\left({\mu }_i,{\xi }_i\right)$ ,$i=1,2$ , where ${\mu }_i$ is a direction parameter and ${\xi }_i$ is a concentration parameter. The likelihood estimates of the direction parameter ${\mu }_i$ is just the first principal component associated of a PCA of cluster i. The identification of the mixture of Watson distribution was obtained by cluster analysis, namely a previous hierarchical cluster analysis followed by a k-means partition of the global sample of variables.

These multivariate data were explained by an alternative factor analysis model potentially delivering directly interpretable solutions without the need of rotations procedures.

The loadings of this factorial model were obtained by regression. The final results concerning communalities of the 16 variables showed that for a great part of them unit variance was quite well explained by the factorial model.

 Keywords: big data, cluster analysis, common factor and residual model, principal component analysis, radioactivity effect, sampling variables, Watson distribution

CA, cluster analysis; CFRM, common factor and residual model; ML, maximum likelihood; PCA, principal components analysis, IVPCA, instrumental variables principal component analysis

Over the last decades, advanced technologies in computer and data science have achieved considerable progress in developing statistical or data mining techniques adapted to analyse structured or non-structured big data. Independently of the different domains covered by the concept of “big data”, nowadays it is quite common to have datasets where the number of variables is much larger than the number of observations. Several applications of high dimensional datasets were analysed in astronomy, chemometrics, climate, finance and genomic.^1–3

 In the present paper the focus will cover situations where the number of observations is pre‐defined, attending the concrete nature of the study, and the randomness concern the choice of variables from a universe of variables following a certain probabilistic model. Such challenge was firstly presented by Hotelling in the context of principal component analysis⁴ and later by Escoufier Y⁵ about the sampling of vectorial variables⁵ and by Gomes in the context of directional probabilistic models associated to an universe of standardized variables defined on the n‐dimensional sphere.⁶ Later Vigneau et al.⁷ have proposed to cluster numerical variables about estimated latent variables, using a k-means type algorithm, obtained a classical factor model estimator in each cluster, being the first principal component of the cluster as a centroid. Such approach was extended to latent variables belonging to a space spanned by external variables, in the context of instrumental variables, in the context of instrumental variables principal component analysis (IVPCA) and multivariate partial least square (PLS) regression.

Recently Xavier B. gave a new contribution for clustering numerical and categorical variables under the hypothesis of a mixture of Von Mises Fisher distributions defined on the n-sphere.⁸

Our paper concerns a cluster analysis of a sample of standardized variables based on a similarity measure of variables j and k defined by

Where $r(x^j,x^k)$ represents the Pearson linear correlation coefficient between variables j and k.

An exponential family of axial distribution defined on the n-spere is an obvious distribution to generate a “bundle of variables” with a given level of interrelation.

The selection of variables from identified sub-groups is, per se, an auspicious way to simplify the learning process of a large set of variables and implicitly leads to a natural dimension reduction, where the first principal component may reflect the “privileged direction” that summarizes sub-groups of variables. These procedures have several applications by reducing the redundancy of variables previously considered for a specific statistical multivariate study, for instance to eliminate certain explanatory variables in a multiple regression model where the multicollinearity problem is present.

In our approach, each sub-group of variables identified by cluster analysis, hierarchical cluster analysis followed by-k-means partition methods⁹ or EM algorithm¹⁰ is considered coming from an axial distribution on the sphere$S_{n-1}$ , the Watson distribution which f.d.p. is defined by

$f\left(x\right)={\left\{{}_{1}F{}_1\left(\frac{1}{2},\frac{n}{2},\xi \right)\right\}}^{-1}\exp \left\{\xi {\left({}^{t}ux\right)}^2\right\}$ ,$x\in S_{n-1},\mu \in S_{n-1},\xi >0$

Where ${}_{1}F{}_1\left(\frac{1}{2},\frac{n}{2},\xi \right)$ is the confluent hypergeometric function defined by

Where $\tau (\cdot )$ is the Tau function,⁶ $\mu$ is a directional parameter, and $\xi$ is the concentration parameter. So, each sample of variables bounded by a “double cone” is supposed a sample of a Watson distribution on the n-spere$w\left({\mu }_i,{\xi }_i\right)$ , i=1,…K.

Consequently, the global sample is a realization of a random vectorial variable having as distribution a mixture of Watson distributions. Additionally, a focused statistical analysis must be made concerning other isolated variables detecting if they are potential discordant variables under the hypothesis of such distributional mixture previously identified, evaluating the effective role of such variables in view of their nature and the study’s objectives.¹¹

Our proposal is no longer related to the well-known problem of “choice of variables à posteriori” but with the choice of variables a priori supposing the goodness of fit to the proposed probabilistic model under the particular context of quantitative variables where the statistical standardization of data makes sense.¹²

From likelihood estimators of parameters $({\mu }_i, {\xi }_i)$ , i=1,…,k we have proposed the formulation of an alternative factorial model.

Where the columns of matrix F are the directional parameters ${\mu }_i$ , $A$ is the loading matrix and $$U$$ the residual matrix which include noise and variables generally weakly correlated with factors and so demanding further statistical analysis facing the previous objectives of a specific study.

The proposed model was applied in multiple contexts^13,14 being an alternative factorial model called factorial model in common factors and residuals (FCFR). The performance of this model is illustrated in the present paper using the so-called Amiard Fish data under radioactivity.

The weapons testing program of several nations regardless of the type of blast, has increased the radioactivity of the seas.

The present study concerns the metabolism of radio strontium of fishes. Strontium, if radioactive, may influence the blood cell formation of many fishes. Over the last decades it has been clearly demonstrated that several fission products are potential hazards from a public health point of analysis. The classical data described here was delivered by Amiard laboratory, related to the Aquatic Ecotoxicology research developed during the last decades.¹⁵

The sample was divided into three aquariums under the same conditions of radioactivity. However, the three aquariums were subject to increase durations of contact with the radioactive pollutant:

 is the aquarium with fishes numbered 1 to 8.

 contains fishes numbered 9 to 17.

And  contains fishes numbered 18 to 24. Fish 17 died during the experiment.

Each fish was referenced by 16 characteristics divided into two groups, the first nine measured at the end of the experiment

Group 1 – Radioactivity characteristics:

Variable 1 – eye radioactivity

Variable 2 – gill radioactivity

Variable 3 – radioactivity of capping

Variable 4 – fin radioactivity

Variable 5 – liver radioactivity

Variable 6 – radioactivity of digestive tract

Variable 7 – kidneys radioactivity

Variable 8 – scale radioactivity

Variable 9 – muscle radioactivity and

Group 2 – Size features:

Variable 10 – weight

Variable 11 – length

Variable 12 – standard length

Variable 13 – head width

Variable 14 – width

Variable 15 – muzzle width

Variable 16 – eye diameter

Statistical analysis

The strong heterogeneity of empirical standard deviations of variables under study, from a minimum of 0.96 (variable 16) to a maximum of 259.09 (variable 6), would justify a previous identification of potential multivariate outliers using the minimum covariance determinant criteria. However, considering the main objectives of our study we have just standardized our data, giving the same weight to each variable.

Hence the variables will be represented on a sphere  (23 active observations). A hierarchical cluster analysis of the sixteen variables followed by a k-means partition, identified two clusters of variables 6: nine radioactivity variables (Group 1) and seven size variables (Group 2).

This means that we have associated our complete sample of variables to a mixture of two Watson distributions.

In the context of sampling variables, the goodness of fit methods for the bipolar Watson distribution was applied to check if the clusters of variables obtained by the previous algorithms come from a Watson distribution.

If x comes from a bipolar Watson distribution,$W_n\left(\mu ,\xi \right)$ , then for large $\xi$ it was shown¹⁶ that $2\xi (1-({}^{\tau }\mu x)x){)}^2\sim x^2{}_{\left(n-1\right)}$ . Simulation statistical research have shown that approximation to $x^2$ distribution it works for moderate values of $\xi$ .⁶

The parameters $\left({\mu }_i,{\xi }_i\right)$ , i=1,2 were estimated by ML method: the estimate ${\widehat{\mu }}_i$ is just the first principal component of group i, (i=1,2) and ${\widehat{\xi }}_i$ is obtained from the equation $Y\left({\xi }_i\right)=\frac{w_i}{p_i}$ ¹⁷ were:

$w_i$ is the highest eigenvalue of principal components of group i, i=1,2, $p_i$ is the number of variables of group i and $Y\left(\xi \right)$ is defined by

$Y\left(\xi \right)=\frac{d}{d\xi }In{F}_1\left(\frac{1}{2},\frac{n}{2},\xi \right)$ ¹⁷

So the factor matrix $\widehat{F}$ .can be written by $\widehat{F}=\left[{\widehat{u}}_1⋮{\widehat{u}}_2\right]$ where ${\widehat{u}}_1$ and ${\widehat{u}}_2$ are not, in general, orthogonal vectors.

Our main objective is to construct an alternative factor analysis model $X=F{}^{t}A+U$ where $A$ is the matrix of loading $\left(px2\right)$ and $U$ is the residual matrix $\left(nxp\right)$ which contains variables supposed not correlated with the factors.

The estimation of loadings was obtained by regression giving the coordinates of variables along the privileged direction generated by vectors ${\widehat{u}}_{i\left(s\right)}$ .

Let be $X$ the standardized data set and considering the theoretical correlation matrix $R$ defined by

So $E(\widehat{{}^{t}FF})={}^t\widehat{F}\widehat{F}=\left[\begin{array}{cc}{}^t\widehat{u}{}_1{\widehat{u}}_1& {}^t\widehat{u}{}_1{\widehat{u}}_2\\ {}^t\widehat{u}{}_2{\widehat{u}}_1& {}^t\widehat{u}{}_2{\widehat{u}}_2\end{array}\right]=\left[\begin{array}{cc}1& r\left(F_1,F_2\right)\\ r\left(F_2,F_1\right)& 1\end{array}\right]$

And the term $AE({}^{t}FF){}^{t}A$ will be estimated by $\widehat{A}E({}^{t}FF){}^t\widehat{A}$

The diagonal elements of this estimated term will give the communalities of our factor analysis model, and so the part of the unit variances of the original variables that were explained by the model.

The elements out of the diagonal of this matrix will give the linear correlation between variables reproduced by our model.

Finally, representing by $R^{\ast }$ the empirical correlation matrix,$(E(\widehat{{}^{t}UU})=R^{*}-\widehat{A}E(\widehat{{}^{t}FF}){}^t\widehat{A}$

Previous hierarchical cluster analysis showed that it was quite realistic to consider just two groups of variables, so that from an arbitrary initial partition into two clusters of equal size, we have achieved a local optima solution using a variant k-means method “la méthode des nuées dynamiques” where the distance function is defined by$D(x,\mu ,\xi )=Const+Lo{g}_1{F}_1\left(\frac{1}{2},\frac{n}{2},\xi \right)-\xi ({}^t\mu x{}^{t}x\mu )$

A stable solution was obtained at the fourth interaction

Group1: radioactivity variables$\left(p_1=9\right)$

Group2: size variables$\left(p_2=7\right)$

The first principal component of standardized variables of group 1 is the ML estimate of directional parameter ${\mu }_1$ and the respective concentration parameter ${\xi }_1$ was estimated by $Y\left(\xi \right)=\frac{w_1}{9}$ where $w_1$ is highest eigenvalue of group 1’s PCA,$w_1=5.03$ , giving${\widehat{\xi }}_1=25.99$

Hence the inertia explained by first principal component is 55.86%.

Similarly, ${\widehat{\mu }}_2$ is the first principal component of PCA of group 2 and ${\widehat{\xi }}_2$ is the solution of equation
$Y\left(\xi \right)=\frac{w_2}{7}$ were $w_2$ is highest eigenvalue of such PCA giving ${\widehat{\xi }}_2=69.98$ . The inertia explained by the first principal component is now equal to 84.16%

Representation of observations on the first principal plan

The factor matrix $F=\left[{\widehat{u}}_1⋮{\widehat{u}}_2\right]$ allowed the representation of fishes on the first factorial plan (Table 1 & Figure1). The linear correlation coefficient between ${\widehat{u}}_1$ and ${\widehat{u}}_2$ ,$r=-0.356$ , explain the non-orthogonality of factors.

Figure 1 Representation of fishes on first factorial plan.

Estimation of loadings and representation of variables on the first factorial plan

The loadings of our model (Table 2) were estimated by regression obtaining the coordinates of variables along the two privileged directions (Figure2) and the communalities, so the part of unit variance of each variable explained by the model.

Figure 2 Representation of variables on first principal plan.

From Table 3 we may conclude that, for group 1, all the variables except radioactivity of digestive tract, kidneys and muscle radioactivity, contribute to the first factor. The exam of the communalities show that these variables are not quite well explained by the model.

We emphasize the fact that the loading matrix reveals the “simple structure” underlying the Amiard data set, where two correlated factors really explain the behaviours of this data set. It means that this model is potentially competitive in such situations, providing directly interpretable solutions, avoiding rotation procedures.

It is quite interesting to check the performance of this model to explain the correlation between variables for each cluster

From$\widehat{R}=\widehat{A}{}^t[{\widehat{u}}_1⋮{\widehat{u}}_2][{\widehat{u}}_1⋮{\widehat{u}}_2]{}^t\widehat{A}$

Table 4 & Table 5 compare the original intra linear correlations and the correlations reproduced by the model.

Interpretation of factor analysis outputs

The Variables {1,2,3,4,5,8} contribute 92.6% to the inertia associated to first factor. All these five variables are strongly negatively correlated with factor 1. In general terms, the factor explains the “radioactivity effect” on the fishes in direct relationship with the duration of such contamination. So, the fishes with smaller factor score, are the most contaminated and the fishes with larger score are the less radio contaminated (Figure1). Complementary, the “size variables” gave a similar contribution to factor 2 (Table 3) and all of them are strongly negatively correlated with such factor. So, the second factor discriminate the smallest fishes of aquarium 1 {1,2,3,4} from the larger fishes {5,6,7,8}. The factorial representation (Figure 1) doesn’t suggest different levels of contamination in these two groups of fishes, except the effect on Variable 4 (fin radioactivity) where the smallest fishes compared to the largest ones, registered, on average, 20% more contamination. And except Variable 8 (scale radioactivity) where the registered variation, among these two sub-groups was about 55%.

The fishes belonging to aquarium 2 (intermediate duration of radio contamination) presented a relatively homogenous behaviour in relation to first factor (Figure1). However, in this aquarium, fish 14 present a clear isolated position, being the smallest fish of the global sample and become particularity affected at the live and digestive tract level (Table 6).

As it was expected, most of the fishes of aquarium 3 {18,19,20,22, 23} presented the highest degree of contamination (Figure 1). In contrast fish 21 and fish 24 suffered a quit smaller contamination considering that these fishes belong to group of the bigger fishes (Table 6, data set). In fact, our previous statistical analysis pointed out the fact that the two factors are negatively correlated.

The starting focus of a classical factor analysis is the correlation matrix which describes the interrelation between the variables under study. The classical factor analysis is defined by:

Where F is the matrix of the latent variables supposed non-correlated and potentially enable to explain a great part of the correlation between the observed variables. In this model, the residual matrix V contains unknown non-correlated variables and also not correlated with the latent variables, representing the specific component of each one of the original variables. The loading matrix informs about the importance of latent factors in their relationship with the variables of the model.

Such hypothesis lead to the results $\Sigma =A{}^{t}A+{\gamma }^2$ where ${\gamma }^2=E({}^{t}VV)$ being $\Sigma$ the variance and covariance matrix.

In practical terms if R is the empirical correlation matrix, the main objective of classical factor analysis model is to find the loadings A and the estimate of variance and covariance of residual terms in order to minimize the difference between R and $\Sigma$ .

Joreskog¹⁸ studied the relative performance of alternative approaches to such optimization problems, namely by the comparison of least square method where we have:

Min $\Delta$ with$\Delta =\frac{1}{2}Tr{(R-\Sigma )}^2$

and generalized least squares where we have:

Min $G$ with$G=\frac{1}{2}Tr{(I_p-R^{-1}\Sigma )}^2$

Or by the maximum likelihood method where we have:

Min $M$ with$M=Tr({\Sigma }^{-1}R)-logdet({\Sigma }^{-1}R)-p$

In the context of large sample sizes, the classical properties of ML estimators have shown the relative preference of this last choice, under the hypothesis of an approximately underlying multinormal distribution of observations. Implicitly, such factor analysis model uses a rational to justify a specific choice of the variables, in general quite connected with the objectives and the priorities of the statistical study.

In the present paper we have supposed that the statistical study doesn’t justify any sampling from individuals because they are previously quite well defined. So, the real problem concerns how to sample variables, namely where we have a very big number of variables. First of all, the treatment of quantitative structured multivariate data suggests a variable cluster analysis procedure just to try to downgrade the problem’s complexity. It means to identify sub-clusters of variables especially intercorrelated. The described framework leads us to the formulation of a probabilistic distribution model associated to each sub-cluster previously identified. In the context of a large number of variables from quantitative structured data, the usual heterogeneity of statistical descriptive indicators or units of measure, justified the previous standardization of data. So it was quite obvious to think about a probabilistic model defined on the sphere n-dimensional, $S_{n-1}$ : The bipolar Watson distribution $W\left(u_1{\xi }_1\right)$ with a direction parameter $u_1$ and a concentration parameter ${\xi }_1$ seemed to be a natural model to explain the stochastic behaviour of each sub-cluster of variables. The fact that the ML estimate $u_1$ of was the first principal component of PCA of sub-cluster i, gave the basis to propose a factorial model derived from the identification of a mixture of Watson distribution on the n-sphere.

The potential presence of some variables considered as statistical discordant in relation to any Watson component of the mixture, justifies a specific statistical analysis to investigate their incorporation in the model facing the main objective of the study. Or eventually asking for an adequate transformation if we are detecting a non-linear correlation with some of the remaining variables. Sometimes it could be more adequate to consider such transformed variable as supplementary variables, so not included on the construction of factors, but afterword’s projected on the factorial plan if the linear correlation coefficients with each factor have some statistically significant meaning. The proposed model

$X=F{}^{t}A+U$ where $F=\left[{\widehat{u}}_1;{\widehat{u}}_2\dots \right]$ is such that the factors are not necessarily orthogonal vectors, meaning that our approach is enables to deliver outputs directly interpretable without a specific rotation procedure.

We have no analytic results to calculate the convenient size of each sub-sample in order to stabilize the associate factor, depending on the value of the concentration parameter and on the dimension of n-spere. However, the simulation work already developed constitutes an interesting platform to face such problem in practical terms.⁶

Recent developments extending our research to qualitative variables⁷ clearly shows the renewed interest of this topic in the age of big data.

There is yet an interesting open problem regarding the joint sampling of individuals and variables in Hilbert spaces.

None.

The author declares that there are no conflicts of interest.

None.

1. Clarke R, Resson H, Wang A, et al. The properties of high-dimensional data spaces: implications for exploring gene and protein expression data. Nature Rev cancer. 2008;8(1):37–49.
2. Johson RA, Wichern D W. Applied multivariate data analysis. 6th ed. New Jersey, USA: Prentice Hall; 2007.
3. Johsonstone I M, Titterington DM. Statistical challenge of high dimensional data. Phil Trans R Soc A. 2009;367:4237–4253.
4. Hotelling H. Analysis of a complex of statistical variables into principal components. J Educ Psychol. 1933;24:417–441.
5. Escoufier Y. These de Doctorat d’État Sciences. France; 1970.
6. Gomes P. Distribution de Bingham sur la n-spere: une nouvelle approache de l’analyse fcatorielle. France; 1987.
7. Vigneau E, Qannari EM. Clustering of variables around latent components. Communication in Statistics–simulation and computations. 2003;32:1131–1150.
8. Bry X, Cucala L. A von Mises-fisher mixture model for clustering numerical and categorical variables. Advances in Data Analysis. 2022;16(2):429–455.
9. Diday E. A new method in automatic classification and pattern recognition the method of dynamic clouds. Journal of Applied Statistics. 1971;19(2):19–33.
10. Figueiredo A, Gomes P. Performance of the EM algorithm on the identification of a mixture of distributions defined on the hypersphere. 2006.
11. Figueiredo A, Gomes P. Discordancy test for the bipolar Watson distribution defined on the hypersphere. Communication in Statistics – simulation and computations. 2006;145–153.
12. Bert DJ. Goodness-of-fit and discordancy tests for samples from the Watson distribution on the sphere. Austral Statistics. 1986;28(1):13–31.
13. Figueiredo A, Gomes P. Clustering of variables based on Watson distribution on hypersphere: a comparation of algorithms. Communication in Statistics – simulation and computation. 2015;2622–2635.
14. Figueiredo A, Gomes P. Classificação de variáveis definidas na hiperesfera através de um modelo de mistura. Proceeding SPE congress. Porto. 2012.
15. Triquet C, Amiard J Mouneyrac C. Aquatic ecotoxicology – Advancing tools for dealing with emerging risks. 1st ed. Elsevir Science. 2015.
16. Mardiak V. Statistics of directional data. London: Academic Press; 1972.
17. Gomes P. Contribution to the problem of the choice of variables in PCA. Montpellier, France: Technical Report nº 8505; 1985.
18. Joreskog KG. Factor analysis by least squares and maximum likelihood methods in statistical methods for digital computers. Wiley. 1977.