Geographically weighted models are widely used for analyzing the spatial data. There is a problem with spatial data processing of extracting a potentially lower number of unobserved variables while a set of correlated variables is observed. The factor analysis is commonly used to overcome this problem. A bandwidth selection is a main difficulty during the identification a spatial heterogeneity of factor loadings. In the paper an original bandwidth selection criterion is proposed. It is based on the testing the difference between factor loadings of global and geographically weighted model. Using the simulated data it is shown that the criterion proposed allows to define accurately the appropriate number of nearest neighbors. The proposed approach is used to analyze real data on performance metrics of Russian universities.
R T
Image Processing, Geoinformation Technology and Information Security / A. Timofeeva, K. Tesselkina where an orthogonal matrix contains the eigenvectors of R , and is a diagonal matrix which entries are the eigenvalues of R . The eigenvalues of diagonal imply the variances of the corresponding principal components. The eigenvectors in are column vectors and represent the loadings of each variable on the corresponding principal component.
If the number of principal components equal to the number of variables, the decomposition (1) perfectly reproduces the correlation matrix R . By reduction m variables in q -dimensional sub-space ( q m ) the correlation matrix is represented as
Rˆq q qqT where q denotes the matrix with the first q columns, i.e. the loadings on the first q principal component, and q is a diagonal matrix which entries are the first q eigenvalues of R . The principal components are sorted in decreasing order of eigenvalues so the first principal components keep the most important information from the data set.
Component scores in q -dimensional sub-space are found by multiplying the original data matrix X by loading matrix . q The best rank q approximation to X is Xˆ q X qqT . A standard result in linear algebra states that
T qqT (q)(Tq) where (q) denotes the matrix with the first q columns removed.
The error matrix is
To assess the quality of the reconstitution of X with q components, the dissimilarity between X and Xˆ q is usually evaluated. The most popular coefficient used for evaluating the quality of PCA model is the residual sum of squares [3] E X Xˆ q X T X qqT X (q)(Tq) .
RESSq X(q)(Tq) where M
is the square root of the sum of all the squared elements of the matrix M .
Thus, mathematically, PCA depends on the eigen-decomposition of positive semidefinite matrices. Its main goal is to extract the important information from the data using the correlation between the variables and to represent it as a set of orthogonal principal components in the sub-space of lower dimension.
EFA model assumes that the relationship between the measured variables is due to the effect of some unobservable (latent) factors. The input information is a correlation matrix R for all variables. It can be represented as [4]
R T (3) where is a factor loading matrix reflecting the relationship between the variables and factors, is a correlation matrix of q factors, is a covariance matrix of the unique factors.
The presence of unique factors in EFA model (3) is the main difference from the model of PCA (1). It is due to the fact that extracted latent factors do not fully (with some errors) describe the correlation between the observed variables. Uniqueness is the variance that is ‘unique’ to the variable and not shared with other variables. The independence of unique factors is assumed, so the matrix is a diagonal with uniqueness on the diagonal.
The matrices , , are estimated. In contrast to the PCA model the matrix is of a particular interest, but loadings are not uniquely determined, so the rotation procedure is used, so that the resulting factor structure has a meaningful interpretation. With orthogonal rotation the independence between the latent factors is assumed. So matrix is given as identity. There are a number of factor extraction methods for estimating loadings and uniqueness, for example, principal factor solution, minimum residual, maximum-likelihood method.
Minimum residual method is based on ordinary least squares (OLS). The loss function is
FOLS , tr R T .
2
Here tr M is the trace of a square matrix M . The OLS-estimates ˆ , ˆ are arguments at which the minimum of loss function (4) is achieved.
The usage of local weighting as part of regression estimation initially was proposed by [5]. This approach is widely used in spatial data analysis [6] and known as geographically weighted model. (2) (4)
Image Processing, Geoinformation Technology and Information Security / A. Timofeeva, K. Tesselkina
Geographically weighted model identifies spatial differences in the relationship between factors by constructing a regression model at each control point for geographically closed observations. The proximity regulated by assigning larger weights to closest points and reducing weights for observations as they move away from the control point. Thus, the weight is determined as a function of distance from the control point to the objects. The regression is estimated over the local subregion which volume is determined by the weight function parameter (a bandwidth).
Regardless of the form of weight function specified, the local correlation matrix is
R(i) (i) (i) (i)T
R(i) (i) (i) (i)T (i) with respect to the local subregion of the i -th control point. The scores for the i -th control point on the m variables are x(i) (i) where x(i) is a vector of variable values at i -th control point.
Similarly geographically weighted EFA model is defined as and values of matrices (i) , (i) are estimated with respect to the local subregion of the i -th control point.
The choice of bandwidth value has a decisive influence on the estimation quality [6]. If someone takes bandwidth too large, then almost all observations will be included in the model, so it will be coincidental to the global model without geographical weighting. Thus it will not be possible to describe the change of the explanatory factors impact depending on the spatial location of the objects. Otherwise, too small bandwidth leads to the overfitting problem: the model perfectly predicts the training data, but drastically fails on some new datasets.
The choice of a bandwidth parameter cannot be based on the common fitting indicators (like R-squared, the mean square error, etc.). So for optimal bandwidth selection the cross-validation technique is often used when the training and quality evaluation both are produced on the distinct sample data [7]. There are some other approaches to solve this problem, for instance, Akaike information criterion [7], and the Lagrange multiplier test [8].
In recent studies on the bandwidth selection [8, 9] for geographically weighted models two essentially different methods to the weighting function construction are considered: • • with a fixed local area radius; with a given number of nearest neighbors.
The second one is considered to be adaptive because it allows adjusting to varying density of the spatial location of the objects. Thus in the neighborhood of one control point the objects may be more concentrated than for the other p oints where they are more distant from each other. In such cases the latching of the local area radius leads to the fact that for some co ntrol points the regression will be estimated on a very large number of observations, for the others – on very small.
We selected the adaptive approach. Therefore, the task is to determine the optimal number of nearest neighbors, taking into account features of EFA.
The criterion ‘goodness of fit’ is proposed for geographically weighted PCA model construction (see [2]). The criterion is based on the minimization of the residual sum of squares. It is calculated by the formula (2) for a global model. For a local PCA model (5) the residual sum of squares at the i -th control point is (5) (6) where a superscript (i) denotes the values that are calculated in a local subregion of the i -th control point. The values of the residual sum of squares are summed for all the control points to calculate the goodness of fit statistics:
RESSq(i) X (i)(i)
(q)((i)qT) n GOF RESSq(i) .
i1
A set of control points can be selected in different ways. Leave-One-Out Cross-Validation (LOOCV) is the simplest procedure for cross-validation. It suggests that only one observation is selected as a control point from the data set, while other observations are considered as a training sample. The procedure is repeated until all the objects will be alternately selected as control points. The advantage of this approach is a computational speed. Often the model structure based on LOOCV leads to the overfitting problem and the forecast error underestimation [10]. In our case the wrong selection of bandwidth parameter may cause such problems.
Image Processing, Geoinformation Technology and Information Security / A. Timofeeva, K. Tesselkina
The more complicated procedure is a Monte-Carlo cross-validation (MCCV) [10]. It assumes that the whole sample is separated randomly into training and check samples. Nevertheless, this choice could increase the computational time. Therefore, we have chosen LOOCV procedure.
There are some problems with usage the cross-validation technique for evaluating the quality of EFA models. It is mainly used in settings where the goal is prediction, and one wants to estimate how accurately a predictive model will perform in practice. The task of variable reduction is not quite a prediction. Of course, we can use the loss function RESSq for PCA, and FOLS , for EFA. But goodness of fit is not so important for EFA, the loadings are more interesting. For this reason a new criterion for bandwidth selection is proposed. It is based on testing the difference between global and local factor loadings.
A statistical inference for comparing global and local factor loadings is based on the information about mean values of loadings and their standard deviation. We need replications of sample data to get this information. One way to get it is to take a sample of the same size n from the rows of data matrix X with replacement. Let we have L replications of sample data. For each l th replication we estimate loading matrix l of global EFA model (3) and loading matrix l(i) of the geographically weighted EFA model (6). We can calculate matrices containing the average values of loadings for all replications The k, j th elements of matrices l and l(i) are denoted by m kj and m k(ji) .
Similarly, we can calculate the variance of loadings for all replications. Let us denote them as v kj and v k(ji) . So the test statistic comparing the means is well known. It is given by
SS 1 n q m
SSk(ji) . n q m i1 j1 k 1 (7) (8) We propose the significance test statistics for optimal bandwidth selection
The maximum value of the significance test statistics (8) corresponds to the optimal number of nearest neighbors. On the one hand it is evident that for the global model (maximum number of nearest neighbors) the numerator of (7) will be minimal, and vice versa. So we would expect that the geographically weighted EFA model with the smallest number of nearest neighbors will be the best by test statistics (7). But on the other hand a small number of nearest neighbors results in very large loadings variation. Thus, the denominator of (7) will increase with a decrease of the number of nearest neighbors. Essentially the significance test statistics (8) is a trade-off between differences in the average factor loadings and their variation.
To calculate statistics (8) we need to compare multiple EFA models. These comparisons require columns of factor loading matrices to be properly aligned. However, the most rotation criteria do not uniquely define a factor loading matrix. This is referred to an alignment problem [4]. The most popular method for aligning a factor loading matrix against another is to minimize the sum of squared differences of factor loadings in the two matrices. Further, in simulation study, this approach was used. We compared the dissimilarity of loading column of global EFA model and one of local models, initial and with the opposite sign. Column reflection (an operation when the signs of values in column are replaced with the opposite) of original column was carried out in cases when reflected loading column corresponded to less value of sum of squared deviations.
There are some problems with the calculation of the statistics (8). Firstly, it is necessary to conduct the Ln-fold factor analysis. With a large number of control points and repeated replications, this procedure takes a very long time. A smaller number of control points can be taken to reduce the calculation time. Another way is to replicate the sample using the jackknife method. However, the decrease in the number of replications may result in loss of the quality of an optimal bandwidth selection.
Secondly, factors can be extracted using various methods. There is a problem with the starting values during the optimization of the log likelihood. The uniqueness is technically constrained to lie in [0, 1], but there are some problems with near-zero values, and the optimization is typically done with a lower bound of 0.005. Sometimes it is unable to optimize the likelihood from certain starting value, because the algorithm does not converge. If we try to increase or decrease the lower bound for uniqueness during the optimization, it allows a solution to be converged. However, such lower bound selection is practically not efficient in the case of Ln-fold factor analysis. So more simple factoring methods should be used.
The method of principal axes may be used in the cases when maximum likelihood solutions fail to converge. However, it is based on the iterative algorithm, so it works rather slowly. If the procedure of factor analysis is repeated many times, the speed
Image Processing, Geoinformation Technology and Information Security / A. Timofeeva, K. Tesselkina of implementation of the factoring procedure is very important. Therefore, optimization procedures are more preferable. In addition, they produce even better solutions for some examples. Minimum residual method based on OLS usually has no problems with convergence and tends to produce better solutions.
Further, two approaches to the bandwidth selection are compared in a simulation study. For identification of PCA and EFA models one can use statistical packages, for instance, the free software for statistical analysis R [11]. The function princomp {stats} performs a principal components analysis on the given numeric data matrix, function efa{EFAutilities} performs exploratory factor analysis. The algorithms of optimal bandwidth selection are implemented using R.
The main purpose of the simulation study is concluded in comparison of approaches mentioned above in precision of the bandwidth selection. A simple one-factor model was chosen. The factor F is standard normally distributed. It affects three variables x1 , x2 , x3 . So the EFA model has the form x1 b1F 1, x2 b2 F 2 , x3 b3 F 3 (9) where b1,b2 ,b3 are factor loadings, 1 , 2 , 3 are random errors. The simulated random error i was chosen as a normally distributed variable with variance that equals to 0.82 bi2 .
The case with certain local centers of object’s concentration was considered for modeling spatial heterogeneity. The whole number of those centers was equal to six, and each of them was represented by a circle with the same number of observations. 2 All center’s locations were chosen randomly within the unit square 0,1 . The sample size was n 300 . The true value of number of nearest neighbors was 50. Different spatial location of objects towards the centers was set in two ways.
In Model 1, the radius of a circle with homogeneous observations was fixed. It was set to 0.05.
Model 2 has a distance from the objects to the center of the local area multiplied by a coefficient 1 0.2K , where K is a serial number of the area, K 1, ,6 . The number of objects belonging to the K-th region was set to 29 6K . This ensures that there are areas with varying density of spatial location of objects.
Equal loadings were set for all objects in the same local area. Their values are shown in Table 1.
The total number of experiments was equal to 50 for each model. The values of both statistics GOF and SS were calculated on the simulated data for fixed number of nearest neighbors. The number of nearest neighbors was set from 20 to 150 with the step of 10. The optimal number of nearest neighbors based on GOF statistics was selected as argument of the goodness of fit statistics minimum constructed as a result of the LOOCV procedure. The optimal number of nearest neighbors based on SS statistics was selected as argument of the significance test statistics maximum constructed as a result of the MCCV procedure. The number L of replications of sample data was set to 100. The size of random sample with replacement was n 300 . The final results are presented in Table 2.
Image Processing, Geoinformation Technology and Information Security / A. Timofeeva, K. Tesselkina 50. Thus, the use of the goodness of fit criterion leads to an explicit overestimation of the bandwidth. In addition, SS statistics also provide a more accurate determination of the number of nearest neighbors. The interquartile range of the optimal bandwidth for it is 30.5. While the optimal number of nearest neighbors, according to the GOF criterion, has an interquartile range of 95 and 67.5, that is 2-3 times more than the results of use the second criterion. Fig. 1.
The dependence on the goodness of fit statistics and the significance test statistics on the number of nearest neighbors
for model 1 (a) and model 2 (b).
It should be noted that values of the goodness of fit statistics vary greatly. We see a lot of local minima and maxima. This fact complicates the use of optimization routines to find the best values of the bandwidth. At the same time, the dependence of the significance test statistics on the number of nearest neighbors appears smoother. This fact allows us to develop more effective optimization algorithms than direct-search method on the grid.
The Ministry of Education and Science of the Russian Federation initiated monitoring of the effectiveness of universities in 2012. Since then, all Russian universities are obliged to provide information on their activities on a number of indicators. The decision on the effectiveness of the university is made depending on whether the university is able to reach thresholds for most indicators. The leadership of different universities can differently determine the priority indicators. It is interesting to determine the structure of universities’ efficiency taking into account regional differences in their activities.
We will rely on the model (9) for determining the structure of performance indicators. We are interested in the factor F that is the overall efficiency of the university. So the observed indicators of activities can be interpreted as x1 is a financial and economic activity: income of the educational organization from all sources per one NDP; x2 is a level of wages of the teaching staff: the ratio of the salary of PPP to the average wage for the economy of the region; x3 is an employment of graduates: the proportion of graduates who have found employment during the calendar year following the year of release, in the total number of graduates of the educational organization who have studied the main educational programs of higher education.
Coefficients b1,b2 ,b3 show the extent to which an particular performance indicator determines effectiveness.
The data from monitoring the effectiveness of educational institutions of higher education for 2015, downloaded from the pages of each individual university, were used as an information base [12]. 571 universities are represented in the sample. They provided information on performance indicators, branches of universities are not included in the analysis. The optimal number of nearest neighbors was chosen based on the SS statistics. The MCCV procedure was used. Control points were located in administrative centers within six federal districts: Central Federal District (CFD), Northwestern Federal District (NWFD), Volga Federal District (VFD), Ural federal district (UFD), Siberian Federal District (SFD), Far Eastern Federal District (FEFD). A total of 67 control points are set. The number L of replications of sample data was set to 300. The size of random sample with replacement was 571. The optimal number of nearest neighbors was 119. The average values of loadings for all replications denoted by m (b1), m (b2), m (b3) are presented graphically in Figure 2 using pie charts. Image Processing, Geoinformation Technology and Information Security / A. Timofeeva, K. Tesselkina
Geographical variations of the factor loadings of universities’ efficiency.
As can be seen from Fig. 2, the structure of indicators of the effectiveness of HEIs varies greatly depending on the territory. Thus, for the regions of Siberia and the Far East, the income of the educational organization has the greatest weight. For many European regions of the country, financial and economic activity does not have such a significant contribution. In most cases, the level of wages of scientific and pedagogical workers is most important in the formation of performance of universities. Only for one northern region of Russia (Arkhangelsk region), employment of graduates predominates in the structure of performance indicators Consequently we can conclude that for most universities the key performance indicators are financial, while the employment of graduates as a result of educational activities does not play such a significant role.
In the paper we propose an original criterion to determine the bandwidth for geographically weighted factor analysis estimation. For this purpose the authors developed a software implementation using the statistical framework R. The investigation of the accuracy of the criterion that determines the optimal number of neighbors is based on the results of the experiments. They show that proposed significance test statistics determines the optimal number of nearest neighbors more accurately. Furthermore the dependence of significance test statistics on the number of nearest neighbors appears smoother. T his makes it possible to develop more effective optimization algorithms for automatic bandwidth selection. The proposed approach is used for education monitoring problems.