Paleoceanography - Fundamentals of Factor Analysis

1 Fundamentals of Factor Analysis: Satellite Image I. Correlation A. Variance and Covariance 1. Variance (of x is denoted Sx2; variance of y is denoted Sy2) is a measure of the scatter of values of a variable about its mean: 2 n S = ∑ ( x − x) / n 2 x i =1 i 2. Covariance (of x and y) expresses the relationship between two variables (a measure of the scatter of values of points in a plane relative to the centroid of the data set): n S = ∑ (x − x )( y − y) / n 2 xy i =1 i i B. Correlation coefficient n 1 x −x y −y 1 ∑ (x − x )(y − y) r= )= )( ∑( n −1 n −1 S S SS n i i i=1 i i i =1 x _ y x y where _ Sx = std. dev. of x = (Σ(xi-x )2/n-1)1/2 and _ Sy = std. dev. of y = (Σ(yi-y )2/n-1)1/2 i.e. Sx and Sy relate the deviations of points from the average relative to the "range" (actually std. dev.) of the observations. In other words, the Covariance divided by the Variance. 3. r2 is "the variance of Y accounted for by its covariance with x" (usually expressed in % units). C. Relation to linear regressions (x on y; y on x; others) 1. Common linear regression of y on x: y = A + Bx Let S = Σ (yi - A - B xi)2 Set ∂S/∂A =0; ∂S/∂B =0; solve for A and B. 2. Matrix math solution of Linear Regression: for eq'n Ax = b (m eq'ns, n unknowns), if columns of A are linearly independent, then: x = ( A A) A b T −1 T a. For example, for the simple linear regression 2 y = C + Dx, where we want to fit pairs of data xi, yi we want to find C, D that mimimize Σ [yi - (C + Dxi)]2 In matrix form, we write the equation y = C + Dx as: _ _ _ _ _ _ | y1 | | 1 x1 | | 1 x2 | | C | | y2 | | . . | | D | = | . | _ _ | . . | | . | |_ 1 xn_| |_ yn _| _ _ _ _ i.e. A x = b Numerical example: _ | | | | |_ _ 1 1 1 1 1 1 2 3 4 5 _ | | | | __| _ _ _ | C | | D | _ _ _ | 4.1 | | 5.9 | | 7.8 | | 10.3 | |_ 12.1 | _ _ = 14 12 10 8 6 4 2 0 0 gives C = 1.92 and D = 2.04 1 2 3 4 5 6 Similarly, to solve the equation y = A + Bx +Cx2: _ _ 2 | 1 x1 x1 | _ _ | 1 x2 x22 | | A | | . . | |_B_| | . . | |_C_| _ _ 2 |_ 1 xn xn _| _ _ _ _ | y1 | | y2 | | . | | . | |_ yn _| _ _ = 3 b. Simple matrix formulas also allow you to compute the estimated uncertainties of the regression coefficients and the correlation coefficients. II. Correlation in n dimensions A. Multiple linear regression (e.g. y = a + bx1 + cx2, the equation for a plane in 3D space) B. r-matrix (later we will refer to this as the matrix Σ) Property 1 2 3 4 5 1 1.00 0.86 0.45 0.83 0.45 2 0.86 1.00 0.74 0.23 0.64 C. Ellipsoids and eigenvectors 3 0.45 0.74 1.00 0.78 0.57 4 0.83 0.23 0.78 1.00 0.39 5 0.45 0.64 0.57 0.39 1.00 (ρij) X2 X1 Bivariate scatter diagram with equal-density contour lines. The inner ellipse encloses 66% of the data points, The outer ellipse 95% of the data points. Figure by MIT OCW. Adapted from source: Joreskog et al. Geological Factor Analysis (1976). A square matrix = x eigenvector λ eigenvalue x One way to find eigenvectors: 1st eigenvector = major axis of ellipsoid 2nd eigenvector = largest minor axis of ellipsoid etc. This works more or less as if we did a regression to get a line that "explains" 4 most of the variance (the dominant linear trend of the data in n-dimensional space), subtracted that regression from the data, then perform another regression to find the next most important contribution to the variance, and so forth. 1. This procedure works because we can rewrite the equation (II.C.) as: (A - λI) x = 0 (Equation II.C.1) In other words, the unknown vector x is orthogonal to all row vectors of (A-λI). The expression for the determinant of A-λI is a polynomial of degree equal to the number of rows and columns of the square matrix A. The roots of the polynomial are the eigenvalues. If A is a real, square, and symmetric matrix, the roots are always real. However, these p eigenvalues may not always be different and some may equal zero. If two or more eigenvalues are equal, we say that they are multiple eigenvalues; otherwise we say the eigenvalue is distinct. Example: consider a 2x2 symmetrical singular (i.e., determinant = 0) matrix; equation II.C.1 is then: | r11-λ | r21 r12 | r22-λ | = 0 (r11-λ)(r22-λ) - r12r21 = 0 λ2 - λ(r11+r22) + (r11r22 - r12r21) = 0 i.e., a simple quadratic equation (also: remember, by assumption of symmetry, r12=r21) Things get out of hand quickly as the matrix gets bigger (that's what computers are for!). 2. Once the eigenvalues are known, the eigenvectors can be calculated from (a.) a. A unique solution cannot be obtained for an eigenvector. If x is a solution, so is cx, where c is a scalar. By convention, eigenvectors are therefore always normalized (unit length). b. Eigenvectors associated with different eigenvalues are orthogonal. Multiple solutions are possible for multiple eigenvalues, but it is always possible to choose an eigenvector which is orthogonal to all other eigenvectors. 3. If the eigenvalues λ are placed as the elements of a diagonal matrix Λ, and the eigenvectors are collected as columns into the matrix U, then eq. (II.C.1) becomes: A U = U Λ _1 2 3.. 1| | 2| | 3| | | n|_ n x n Square Matrix n_ | | | | | | | _| _1 2 n_ | EV1 EV2 … | | ↓ ↓ | | . . | | . . | = | . . | | . . | | . . | |_ . . _| 5 _1 2 3... n_ _1 2 3... n_ | EV1 EV2 … | |λ1 0 0… | | ↓ ↓ | |0 λ2 | | . . | | . | | . . | | | | . . | | . | | . . | | | | . . | | | |_ . . _| |_ λ n| 4. The matrix U (the eigenvector matrix) is square orthonormal (the matrix is nxn and the eigenvectors are of unit length), so U UT = UT U = I, therefore: A = U Λ UT Therefore any symmetric matrix such as the correlation coefficient table can be expressed in this form: Σ = U Φ UT m x m m x N N x N N x m r-table eigenvectors eigenvalues eigenvectors(transpose) III. Principle Components Analysis: basically, the principle components are the eigenvectors as outlined above: A. PCA is inherently variance-oriented; it accounts for maximum variance of all the observed variables. In other words, PCA accepts that large part of the total variance of a variable is important and common to other observed variables. B. Factor analysis (below) is correlation-oriented; it accounts for the maximum intercorrelation of variables. In other words, factor analysis allows for a considerable "amount of uniqueness" to be present in the data and utilizes only that part of a variable that takes part in correlation with other variables; i.e., to account for covariance rather than variance. Example: imagine a lake with three stream systems draining into it. Suppose that each of the drainage basins has a distinctive (average) sediment composition, and that the sediments within the lake can be described as linear combinations of the sediments derived from each stream system. What is interesting in this case is the distribution of the three stream system component sediments, not the mineralogy or chemical composition per se (i.e., we don't care if Al correlates with Si). So we adopt a framework that lets us describe the raw data (chemical or mineralogical composition) in terms of sums of source components, and map the distribution of the components. 6 Source of example: Joreskog et al. 7 IV. Species space and sample space A. Samples in species space: PROPERTY 2 sample 1 PROPERTY 1 sample 2 PROPERTY 3 B. R Mode (species in sample space): sample 2 PROPERTY 1 sample 1 PROPERTY 2 sample 3 8 C. What is sample space? 1. Consider the species percentage data: PaL Sample 1 17.4 Sample 2 36.4 PaR 17.4 22.7 Gbu 39.1 27.3 Ndu 26.1 13.6 Average 20.1 33.2 19.9 26.9 2. Express each species percentage in terms of its deviation from the average for the samples: PaL Sample 1 -9.5 Sample 2 +9.5 PaR -2.7 +2.7 Gbu +5.9 -5.9 Ndu +6.2 -6.2 3. Transpose: PaL PaR Gbu Ndu Sample 1 -9.5 -2.7 +5.9 +6.2 Sample 2 +9.5 +2.7 -5.9 -6.2 4. Row-normalize (to unit vector length): PaL PaR Gbu Ndu Sample 1 -0.71 -0.71 +0.71 +0.71 Sample 2 +0.71 +0.71 -0.71 -0.71 length 13.42 3.77 8.38 8.80 9 5. Plotting the species as vectors in sample space: Note that the PaL and PaR vectors coincide, as do the Gbu and Ndu vectors. 6. In sample space, species that are highly correlated fall in the same vector region, i.e., there is a very small angle between highly correlated species. a. The correlation coefficient is the cosine of the angle between the two PROPERTY vectors in sample space: Let a and b any two PROPERTY vectors (with a common origin). Then aTb ρ = cos θ = __________ |a| |b| 10 So in this case, the correlation coefficient table is: PaL PaR Gbu Ndu PaL +1 +1 -1 -1 PaR +1 +1 -1 -1 Gbu -1 -1 +1 +1 Ndu -1 -1 +1 +1 So (PaL and PaR) and (Gbu and Ndu) are positively correlated whereas (PaL and GBu), (PaL and Ndu), (PaR and Gbu) and (PaR and Ndu) are negatively correlated. This method of calculating the correlation coefficient is mathematically equivalent to the "normal" definition of the correlation coefficient given at the beginning; i.e., this sample space diagram is the basis for the definition of the correlation coefficient! When there are only two samples as in this example, this definition seems as there are only three possibilities: r = -1, 0, +1. This is because any two data pairs define a straight line (unless they coincide, in which case no line is defined) and this line has either a positive slope or a negative slope. But when there are more samples (i.e. multidimensional sample space), the data pairs no longer have to define a line, the vectors can occur in any region of sample space, and r can have any value –1 < r < 1. If A is the matrix with the row-normalized (i.e., unit length so that |a| = |b| =1) vectors of species data (relative to the mean for each species) in sample space, then the correlation coefficient table is the inner product of A: Σ=ΑΑ Τ 7. Likewise, one can define a similarity coefficient as the cosine of the angle between samples in PROPERTY space; it is a measure of how similar two samples are to one another in their PROPERTY composition. We will call the similarity coefficient matrix H. 8. What we are going to do next is reduce the number of "species" by combining them into "factors". In other words, we are going to reduce our system from having (say) 31 species counts for each sample to having (say) 6 "factor loadings" for each sample. Continuing on from equation II.C.4., a math theorem says that a symmetric matrix like the correlation matrix Σ can be expressed as: Σ = A Φ AT Ψ + where "Σ is the p by p population covariance matrix of the observed variables" (the symmetrical correlation coefficient matrix), "A is the p by k matrix of factor loadings, Φ is the k by k covariance matrix for the factors (if the factors are in standardized form, this is a correlation matrix with ones in the diagonal), and Ψ is the p by p residual covariance matrix." Source: Joreskog et al. 9. If we force the solution to be orthogonal (factors are uncorrelated), then Φ = I and hence: Σ = A AT + Ψ 11 Ψ is now a diagonal matrix (i.e. there is covariance between the factors) 10. A math theorem tells us how to calculate the Singular Value Decomposition (SVD) of a matrix: Let us assume that we have an N x p (N>p) data matrix X. The product moments are defined as: major product moment (MPM) V = X XT minor product moment (mPM) U = XT X Then we can decompose the data matrix as: X=VΓU where V is an N by r matrix with orthonormal columns U is a p by r matrix with orthonormal columns, and Γ is a diagonal matrix of order r by r with positive diagonal elements γ1, γ2, ...γr called singular values of X. The major product moment XXT, which is square, symmetric, and of order N by N has r positive eigenvalues and (N-r) zero eigenvalues. The positive eigenvalues are γ12, γ22, ..., γr2 and the corresponding eigenvectors are u1, u2, ..., ur. The minor product moment XTX, which is square, symmetric, and of order p by p has r positive eigenvalues and (p-r) zero eigenvalues. The positive eigenvalues are γ12, γ22, ..., γr2 and the corresponding eigenvectors are v1, v2, ..., vr. The positive eigenvalues of X XT and XT X are the same, namely γ12, γ22, ..., γr2 . Furthermore, if vm is an eigenvector of X XT and um and eigenvector of XT X corresponding to one of these eigenvalues γm2, then the following relationships hold between um and vm : vm = (1/γm) X um and um = (1/γm) XT vm These relationships make it possible to compute vm from um and viceversa, i.e. V = X U Γ-1 and U = XT V Γ-1. The analysis of the minor product moment XTX is referred to as R-mode analysis and that of the major product moment XXT as Q-mode analysis. 11. Major problem with this analysis: there are an infinite number of solutions! If we find a solution to the problem (a set of orthogonal vectors that describe the data), then any rotation of that set of vectors is also a solution! How do we choose any one of these? a. Rotation of a set of row vectors in the matrix X can simply be done through the operation: Y = X R where R is a rotation (transformation) matrix conforming to the requirement that R RT = I for a rigid rotation. 12 b. In the end, the choice is arbitrary (this is one of the problems with factor analysis). However, certain more-or-less reasonable choices can be made: i. As much as possible, try to make the factors simple, i.e., have a few high loadings and many zero or near-zero loadings. ii. Rotate the axes so as to put as much variance as possible into the factors (VARIMAX criterion). i.e., we try to "explain" as much of the variance as possible with the fewest possible factors. This solution is often favored because it is an objective solution (i.e., it arises untouched by human hands). That doesn't necessarily make it a better solution, however. 12. Brief summary of factor analysis steps: a. Start with a data table with rows of samples and columns of species b. Row-normalize the data table so that each row has unit length. (it is now referred to as W). c. Calculate the minor product moment (mpm) = WT W d. Calculate the eigenvalues (Λ) and eigenvectors (U) of the mpM. e. We can then express each sample in terms of the eigenvectors by multiplying: W U (the result of this calculation is a list of the eigenvector composition of each sample; e.g., each sample can be expressed as a linear combination of the eigenvectors f. We now choose to simplify the properties of our data set by throwing out all but the first few (dominant) eigenvectors, and then describe our data set in terms of these. This is now the factor matrix V, which describes the loadings of each species to each factor. g. Because we have thrown out some information, the factor loading matrix A = W V does not fit the original data exactly. h. In order to improve the fit somewhat, we will now rotate the factor matrix V so that it accounts for as much of the variance as possible. This is done through an iterative process: we rotate around axis one to find the maximum, then rotate around axis two to find the maximum, etc. for all of the axes (= # of factors). This process is repeated until the result converges (no further improvement in fit is found). i. Our final rotated factor matrix (the Varimax solution) is F, a list of the composition of each species described as a linear combination of the individual factors. The samples are now described as 13. One problem with this type of analysis is that the number of factors is arbitrary! The geologist's approach to choosing number of factors: make maps. Keep adding factors until additional factors become "unmappable". 14. Negative factor problem: the mathematics don't specify that you only can add things; the math is just as happy subtracting a factor. In chemical and physical reality, a sample can only be the sum of positive components (philosophically, it may be possible to subtract components, e.g., incongruently dissolve something away to create a "negative" component of a different composition than the preexisting components). Various approaches to this problem: (a) ignore the problem (most common), (b) rotate solution so that factors are (mainly) positive (but which rotation should you choose?) 15. Another problem: summation to 100% and the "...and everything else factor"