1. Causal structure in data

Suppose that (x₁, x₂, … , x_J) are random variables with a covariance matrix G. Suppose further that they depend on the variables (t₁, t₂, …) in a following way,

                        x₁ = p₁₁ t₁ + p₁₂ t₂ + …

                        x₂ = p₂₁ t₁ + p₂₂ t₂ + …

(1)                    …

                        x_J  = p_J1 t₁ + p_J2 t₂ + …

In a vector/matrix notation (1) can be written as x=Pt. (When working with data it corresponds to the rows of X=TP¢). If the t-variables are mutually uncorrelated, the equations (1) imply that we can write G=PDPN, where D is a diagonal matrix containing the variances of the t's.  The number of t-variables equals the amount of columns in P. This amounts to stating that the score vectors (t_i) in (1) are mutually orthogonal, t_iNt_j=0 for i¹j. (Note that this D may be different to the D in (1). D in (1) is a scaling matrix that depends on the procedure and the actual choice of scaling). The coefficients in P (K´A) tell us how the variables are related. Following table shows three examples.

Table 1. Schematic illustration of sizes of loading coefficients, P.

Cells in the table marked by x indicate that the corresponding loading coefficient is different from zero. The unfilled part thus consists of zeros. In the first example there are only non-zero loadings for the first latent variable. The variables x₁-x₃ have non-zero loadings on the second latent variable, while x₄-x₆ on the third. We see this type, when we work with psychometric data. The first latent variable explains all the variables, while the variables can be grouped in relation to the other latent variables. In the second example each of the ten variables can be explained by two latent variables. In the third examples the first latent variable explains x₁-x₃, while the second explains x₄-x₆. It indicates that the two sets of variables are mutually independent.

            The non-zero coefficients can be generated as the ones that are not statistically significant. We can also prescribe certain coefficients to be non-zero and estimate the loading in relation to this pattern of non-zero loading values. In this case we study how well the loading pattern describes the present variables.

            The latent variables are derived as linear combinations of the original variables as shown in the following set of equations,

                        t₁ = r₁₁ x₁ + r₂₁ x₂ + … + r_J1 x_J

(7)                    t₂ = r₁₂ x₁ + r₂₂ x₂ + … + r_J2 x_J

                        …

If the t's are uncorrelated, the equations (7) imply that D=RNGR, where D is the diagonal matrix of the variances of the t's. The equations (7) are useful, when you want to study the independence and conditional independence among the x-variables.

            Sometimes we require a special structure in the coefficient matrices P or R. An example is when we require P to be lower triangular. In that case the equations (1) can be written as

                        x₁ = p₁₁ t₁

                        x₂ = p₂₁ t₁ + p₂₂ t₂

x₃ = p₃₁ t₁ + p₃₂ t₂ + p₃₃ t₃

(8)                                       …

                        x_J  = p_J1 t₁ + p_J2 t₂ + …

This makes interpretation of the loadings simpler. E.g., x₁ is the same as t₁, x₂ depends on x₁ and x₂, and so on. Also, e.g., p₃₃ can be interpreted as the part of x₃ that cannot be described by x₁ and x₂. When P is lower triangular, R will also be.

            The presentation in this section has been based on one block of variables (x₁,x₂, …) or one block of data X that is being described by latent variables (t₁,t₂, …) or one block of score data, T. But we will be working with several blocks of X-data. The same considerations (e.g. patterns of non-zeros) or special structure (e.g. lower triangular loading matrix) apply to each block.

2 Measures of fit and changes

The H-principle and PLS are based on least squares approach. In contrast many statistical procedures use maximum likelihood, ML, as the estimation method. We shall briefly review the ML approach, because the way we work with J-divergence (see below) is similar to the ML approach. In the ML approach the log-likelihood function is often used to measure fit and changes in parameter estimates. It is given by

(3)                    N(S) = c - N (log (|S| + tr(S^-1S) )/2,

where c is a constant and S the sample covariance matrix, S=X^TX/N. From a theoretical point of view it is comfortable to work with the function (3). If S₁ is the maximum likelihood estimate for S under null-hypothesis and S₂ the estimate under a (reduced) hypothesis, the theory tells us that the difference

(4)                    2(N(S₂) - N(S₁)) = - N (log(|S₂|/|S₁|) + tr ((S₁^-1 - S₂^-1)S))

approximately will follow a x²-distribution with degrees of freedom equal to the difference between the number of parameters at the two hypothesis. (4) is often used to test the effect of adding the number of parameters by 1 and comparing (4) to the critical value of x²-distribution with one degree of freedom. The disadvantage of using this procedure is that the function (3) is unstable, when data has reduced rank. It is easy to show examples, where arbitrary small score vectors become significant, when using (3) or (4), although they have no predictive ability. It is instructive to look closer at the differential of (3),

(11)      dN = N tr( S^-1 - S^-1S S^-1)d S

The equation shows that optimal S is one, where S is a generalised inverse of S^-1, S^-1 = S^-1S S^-1.

The equation (5) shows that S^-1S should be close to the identity matrix I. When we work with reduced rank, we also need that SS^-1 is close to I. This can be obtained by working with J-divergence,

(12)      J = tr( S S^-1 +  S^-1S - 2 I)/2

The theory associated with (6) is the same or similar to (3) and (4). It is derived from testing equality of two covariance matrices, see ref 10. It is instructive to look at (6) in terms of the eigen values of S S^-1. If (l_i) are the eigen values, it is simple to show that

(7)                    J = ½ S (1 - l_i)²/l_i

The equation (7) shows that the ratio of variances from S and S^-1 should be both close to one and different from zero, if J in (7) is to be small. The differential of J is given by

(8)      dJ = ½ tr( S^-1 - S^-1S S^-1)d S

If we compare (5) and (8) we see that (8) will be close to (5), when the inverse S^-1 is close to S^-1. But the advantage of working with (6) and (8) is that the expressions can be used also in the case of reduced rank in data. But it should be noted, that the theory is different in the case of reduced rank. Furthermore, we work with the expressions in the case of reduced rank in the same way as in case of full rank. Thus, the criteria we use are well motivated, when data are not of full rank.

3 Sequential estimation of loading coefficients

We shall briefly describe the situation, when revising estimates of the loading matrix. The algorithm has given a loading matrix P=P₀, such that the sample covariance matrix S is approximately S@P₀P₀^T. Here we suppose that the score vectors have been scaled to unit length. The algorithms associated with the H-principle always provides with a matrix R=R₀ such that R₀^TP₀= I, where I is the identity matrix with number of diagonal elements equal to the number of score vectors. The task we consider here is to find a simple or a significant form of the loading matrix P. The matrix P that we find may contain only significant non-zeros. It may also contain the best estimates according to some pattern of non-zeros, cf. Table 1. Or it may contain a combination of these two, i.e., significant loading values in a given pattern of non-zero loading values.

            We shall use the following notation:

            S=P₀P₀^T,        S⁺=R₀R₀^T,     S=PP^T,           S⁺=RR^T.

The matrix R is computed as R=P(P^TP)^-1. Note that the number of columns in P can be smaller than the number of columns in P₀. There are three steps in filling out P:

1.      Find the next element in P that should be considered as non-zero.

2.      Estimate the values of all non-zero elements in P.

3.      Judge the improvement and closeness to P₀.

There are many ways to handle each of the three steps. We shall here only consider one approach of each.

Step 1. We are to find a new loading matrix P₁, P₁=P + c E_ij, where E_ij has 1 as element (i,j) and zero for others. It is natural to choose the next non-zero element that gives the maximal increase in [tr(P₁P₁^TS⁺)-tr(PP^TS⁺)]. If we solve the equation ¶[tr(P₁P₁^TS⁺)-tr(PP^TS⁺)]/¶c=0, it is shown in Appendix 1 that a solution is given by

(9)      c = - p_j^Tu_i/(u_ii)

Here p_j is the j-th column of P and u_i the i-th column of U=S⁺=(u_ij). It gives an optimal increase as,

(10)      [tr(P₁P₁^TS⁺)-tr(PP^TS⁺)] = (p_j^Tu_i)²/(u_ii)

The task is thus to find the element in P that gives the maximal value of

(11)                  max  (p_j^Tu_i)²/(u_ii),        for (i,j) within the pattern of actual or required non-zeros.

Step 2. Previous step gave a candidate for the next non-zero element. We need to revise the estimates of present non-zero elements in the light of the new one. If we differentiate J with respect to an element in P, we get from (8)

                        ¶J/¶p_ij = ½ tr( S^-1 - S^-1S S^-1)(e_i p_j^T + p_j e_i^T)

Here we need an expression for ¶S/¶p_ij. From S =PP^T=p₁ p₁^T + p₂ p₂^T + …, we derive ¶S/¶p_ij =e_i p_j^T + p_j e_i^T. Here e_i=(0,.,0,1,0,.), where the 1 is at the i-th element. If we require ¶J/¶p_ij=0, we get

tr( S^-1 p_j e_i^T) = tr(S^-1S S^-1 p_j e_i^T) =  tr(S^-1S RR^T p_j e_i^T) = tr(S^-1S R e_j e_i^T)

This gives

                        u_i^T p_j = (S^-1S R)_ij                   ( the (i,j)-th element of S^-1S R).

If we write out the linear equations, we get

(12)                 u_1i p_1j + u_2i p_2j + ... + u_Ji p_Jj = (S^-1S R)_ij

The equation (12) is used to find p_j. Some of the p-values in (12) can be zero and therefore are not in the equations. There is one equation for each non-zero element in p_j. Therefore (12) will give us a quadratic coefficient matrix that is positive definite. The part of the coefficient matrix U=S⁺ that is used in (12) will typically be of reduced rank. Therefore, some care must be given in finding the solution p_j. Note that both R=P(P^TP)^-1 and E^-1=RR^T depend on P. Thus, some few iterations are necessary to get convergence. The starting values are the previous values of non-zero elements in P and the value c in (9) for the estimate of the last non-zero one. The task of (12) is to provide with a small adjustment of the values of the non-zero elements.

Step 3. The starting point is the expression for the J-divergence, (8), that we need to truncate to fit data with reduced rank. The products S S^-1 and S^-1S  are both positive definite. They can therefore be written as CC^T and C^T-1C^-1, resp. It gives

                        J = tr( S S^-1 +  S^-1S - 2 I)/2 = tr( (C - C^T-1)(C - C^T-1)^T)/2

When we truncate, we use C as C@R^TP₀ og C^T-1@(R₀^TP)^T. It gives

(13)                 J @ tr( (R^TP₀ - (R₀^TP)^T) (R^TP₀ - (R₀^TP)^T)^T)/2 = S f_ij²/2,

where F=R^TP₀ - (R₀^TP)^T. Both R^TP₀ and (R₀^TP)^T= P^TR₀ approach the unity matrix of appropriate size, when P is filled out (and R computed as R=P(P^TP)^-1). The filling out of P stops, when J is sufficiently small or F can be judged as a random matrix.