Q-Mode Factor Analysis 

of Geochemical and 

Petrologic Data Matrices 

with Constant Row-Sums

GEOLOGICAL SURVEY PROFESSIONAL PAPER 574-G
Q-Mode Factor Analysis 

of Geochemical and 

Petrologic Data Matrices 

with Constant Row-Sums

By A. T. MIESCH

STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

GEOLOGICAL SURVEY PROFESSIONAL PAPER 574-G

An extension of the method of 

Q-mode factor (vector) analysis 

to increase its usefulness in 

geochemical and petrologic 

investigations

UNITEDSTATESGOVERNMENTPRINTINGOFFICE, WASHINGTON: 1976
UNITED STATES DEPARTMENT OF THE INTERIOR 

THOMAS S. KLEPPE, Secretary

GEOLOGICAL SURVEY

V. E. McKelvey, Director

Library of Congress Cataloging in Publication Data

Miesch, Alfred T.

Q-mode factor analysis of geochemical and petrologic data matrices with constant row-sums.

(Statistical studies in field geochemistry) (Geological Survey Professional Paper 574-G)

Bibliography: p.

Supt.ofDocs.no.: I 19.16:574-G

1. Geochemistry-Statistical methods. 2. Petrology-Statistical methods. 3. Vector analysis.

I. Title. II. Series. III. Series: United States Geological Survey Professional Paper 574-G.

QE515.M54 551.9; 0l'82 75-619414

For sale by the Superintendent of Documents, U.S. Government Printing Office

Washington, D.C. 20402 

Stock Number 024-001-02909-0
CONTENTS

Abstract ...........

General introduction.

Part I: Mathematical Basis Where the Number of Factors 

Equals the Number of Variables

Introduction..............................................

Q-mode factor analysis.....................................

Scaled scores and composition loadings .......................

Composition scores........................................

Negative loadings and negative scores ........................

Finding the composition scores for any vector in the factor space.. 

Finding the vector representation for a given composition........

Computation of initial and composition loadings with respect to 

oblique reference vectors .................................

Summary ................................................

Part II. Effects of Reducing the Number of Factors

Introduction..............

Estimates of factor loadings.

Gl

1

2

3

5

8

10

11

11

12

13

13

14

16

17

18

21

21

22

23

23

24

Page 

Interpretation of eigenvalues and communalities ............... G15

Scaling the factor loadings and scores.........................

Goodness-of-fit measures for the factor model.................

An example of the effect of the constant row-sum constraint .....

Common and unique factors ................................

Compositions versus composition scores ......................

Identification of the samples of extreme composition ...........

Modification of a varimax model to an oblique model...........

Testing the plausibility of alternative end-members .............

Summary ................................................

Part III. Applications to Some Petrologic-Mixing Problems

Introduction..............................................

Rhyolite-basalt complex on the Gardiner River, Yellowstone

National Park, Wyo. ....................................

Granitoid intrusive in the southern Snake Range, Nev. ..........

Lavas and pumices of the 1959 summit eruption at Kilauea,

Hawaii ................................................

Layered series of the Skaergaard intrusion, Greenland ..........

Summary ................................................

References cited...........................................

25

25

29

34

39

46

46

ILLUSTRATIONS

Page

FIGURE 1. Diagram showing the relationships between the factor model and the various forms of the original and reproduced data........... G3

2. Factor-variance diagram derived from a closed matrix of random numbers having column means and variances similar to those of

the data matrix for the granitoid intrusive in the southern part of the Snake Range, Nev. ................................. 21

3. Factor-variance diagram for the rhyolite basalt complex on the Gardiner River, Wyo. .................................... 26

4. Q-mode vector diagram for the rhyolite basalt complex on the Gardiner River, Wyo........................................ 27

5. Factor-variance diagram for the granitoid intrusive, southern Snake Range, Nev. ........................................... 30

6. Q-mode vector diagram for the granitoid intrusive ..................................................................... 30

7. Factor-variance diagrams derived from data obtained by mathematical mixing of sample 1 of the granitoid intrusive and randomly

selected samples of dark inclusions from the intrusive, Osceola Argillite, and Pioche Shale ............................... 32

8. Fact or-variance diagram for the lavas and pumices of the 1959 summit eruption at Kilauea, Hawaii............................ 34

9. Q-mode vector diagrams for lavas and pumices of the 1959 summit eruption at Kilauea, Hawaii ............................... 36

10. Graph showing communalities of olivines in two- to five-factor varimax solutions for the lavas and pumices of the 1959 summit

eruption at Kilauea, Hawaii .................................................................................... 37

11. Factor-variance diagram for the layered series of the Skaergaard intrusion ................................................. 4^

12. Diagram showing molar compositions in the system albite-anorthite-orthoclase that have communalities in excess of 0.9988 in the

eight-factor varimax solution for the layered series of the Skaergaard intrusion.......................................... 41

13. Graph showing communalities of olivines in the eight-factor varimax solution for the layered series of the Skaergaard intrusion .... 41

14. Q-mode vector diagram for a series of augite specimens from the layered series of the Skaergaard intrusion ..................... 42

15. Graph showing communalities of theoretical mixtures of ferrous iron oxide and magnetite in the eight-factor varimax solution for

the layered series of the Skaergaard intrusion ...................................................................... 43

16. Diagram showing interpretation of the net gains and losses of minerals in the parent magma plotted against structural height of

sample in the layered series of the Skaergaard intrusion ............................................... 45

in
IV CONTENTS

TABLES

TABLE 1. Original hypothetical data, transformed data, and normalized data used to illustrate the methods of computation................ G4

2-4. Initial loadings, unsealed scores, scale factors, and the estimated normalization factors for the three-factor principal components, 

varimax, and oblique factor models derived:

2. Without transforming the data................................................................................ 6

3. After transforming each variable to proportions of the maxiumum value ............................................ 6

4. After transforming each variable to proportions of the range ...................................................... 7

5. Composition loadings and composition scores for the three-factor principal components, varimax, and oblique factor models

derived without transforming the data............................................................................ 8

6. Composition loadings, scaled scores, and composition scores for the three-factor principal components, varimax, and oblique

factor models derived after transforming each variable to proportions of the maximum value ............................. 9

7. Composition loadings, scaled scores, and composition scores for the three-factor principal components, varimax, and oblique

factor models derived after transforming each variable to proportions of the range ...................................... 9

8. Computed composition scores for some arbitrary vectors in the varimax solution of tables 4 and 7 ............................. 11

9. Vector loadings in the varimax solution of tables 4 and 7 for some given compositions ....................................... 12

10 12. Initial loadings, unsealed scores, scale factors, estimated normalization factors, and reproduced data matrix for the two-factor 

principal components, varimax, and oblique factor models derived:

10. Without transforming the data................................................................................ 15

11. After transforming each variable to proportions of the maximum value ............................................. 16

12. After transforming each varable to proportions of the range ....................................................... 17

13. Eigenvalues of the matrices of coefficients of proportional similarity (cosine theta) derived from the data matrices in table 1, and

CPN ........................................................................................................ 17

14. Composition loadings, composition scores, and reproduced data matrix for the two-factor principal components, varimax, and

oblique factor models derived without transforming the data ......................................................... 18

15. Composition loadings, scaled scores, composition scores, and reproduced data matrices for the two-factor principal components,

varimax, and oblique factor models derived after transforming each variable to proportions of the maximum value .......... 19

16. Composition loadings, scaled scores, composition scores, and reproduced data matrices for the two-factor principal components,

varimax, and oblique factor models derived after transforming each variable to proportions of the range ................... 20

17. Measurements of correspondence between the observed data and the data reproduced from the two-factor models ............... 20

18. Means and standard deviations for four sets of data examined by Q-mode factor analysis..................................... 25

19. Analyses of a rhyolite-basalt complex on the Gardiner River, Yellowstone National Park, Wyo., adjusted so that eight major

oxides sum to 100 percent, and corresponding data reproduced from the factor model ................................... 27

20. First eight eigenvalues of the cosine theta matrix for the rhyolite-basalt complex on the Gardiner River and CPN ................. 27

21. Composition scores for some vectors in figure 4........................................................................ 28

22. Compositions of the end-members for the factor model of the rhyolite-basalt complex on the Gardiner River .................... 28

23. Composition loadings and increments of Na2O and H2O+ for the factor model of the rhyolile-basalt complex on the Gardiner

River ....................................................................................................... 28

24. Goodness-of-fit statistics for the factor model of the rhyolite-basalt complex on the Gardiner River ............................ 28

25. Analyses of selected samples of the granitoid intrusive in the southern Snake Range, Nev., adjusted so that nine major oxides sum

to 100 percent, and corresponding data reproduced from the factor model ............................................. 29

26. First nine eigenvalues of the cosine theta matrix for the granitoid intrusive in the southern Snake Range and CPN .............. 30

27. Average compositions of some sedimentary rocks and dark inclusions from the southern Snake Range and computed communalities 31

28. Values of A 2 for the average compositions of dark inclusions and some sedimentary rocks from the southern Snake Range in models

containing two to nine factors................................................................................... 31

29. Composition scores of some vectors in figure 6 ........................................................................ 33

30. Compositions of the end-members for the Snake Range factor model ...................... .. ............. 33

31. Composition loadings and increments of Na2O and H 2O + for the Snake Range factor model ................................ 33

32. Goodness-of-fit statistics for the Snake Range factor model ............................................................. 33

33. Analyses of lavas and pumices from the 1959 summit eruption at Kilauea, adjusted so that nine major oxides sum to 100 percent,

and corresponding data reproduced from the factor model .......................................................... 35

34. First nine eigenvalues of the cosine theta matrix for the lavas and pumices of the 1959 summit eruption at Kilauea and CPN ....... 35
CONTENTS V

Page

TABLE 35. Compositions of two samples of spatter from the 1959 summit eruption at Kilauea, adjusted so that nine major oxides sum to 100

percent, and composition scores of vectors representing these compositions in the three-factor Kilauea model ............... G35

36. Initial loadings and composition scores for selected vectors in the varimax solution for Kilauea, and compositions of some

volcanic materials............................................................................................. 37

37. Chemical compositions and norms of the end-members for the Kilauea factor model ........................................ 38

38. Composition loadings and increments of K2O for the Kilauea factor model ................................................ 39

39. Goodness-of-^it statistics for the Kilauea factor model.................................................................. 39

40. Analyses of samples from the layered series of the Skaergaard intrusion, adjusted so that nine major oxides sum to 100 percent..... 40

41. First nine eigenvalues of the cosine theta matrix for the layered series of the Skaergaard intrusion and CPN ..................... 41

42. Compositions of end-members for the Skaergaard factor model.......................................................... 42

43. Composition loadings for the Skaergaard factor model, average loadings for each zone of the layered series, and weighted

averages for all zones .......................................................................................... 43

44. Interpretation of chemical and mineralogical modifications of the parent magma required to produce sample 5166 from the

Skaergaard intrusion........................................................................................... 44

STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

Q-MODE FACTOR ANALYSIS OF GEOCHEMICAL AND PETROLOGIC 

DATA MATRICES WITH CONSTANT ROW-SUMS

By A. T. MlESCH

ABSTRACT

Matrices of data representing all the major constituents in a suite of 

rock samples tend to have constant row-sums, a property that has caused 

considerable difficulty in attempts to interpret rock genesis from the 

correlations among matrix columns. Q-mode factor analysis may be used 

to interpret rock genesis from relations among matrix rows, and the 

constant row-sum is a definite asset. The constant can be used to compute 

sealers, which, in turn, can be used to adjust the factor loadings and 

scores to conform with the original data. Therefore, principal 

components and varimax models and a variety of oblique models may be 

derived that can be used to recompute estimates of the data in, for 

example, percent or parts per million. Any real or hypothetical 

composition may also be tested as a possible end-member in a petrologic 

mixing problem.

Other advantages of scaling the factor loadings and scores to conform 

with the original data, and thereby deriving composition loadings and 

scores, are that (1) the loadings sum to unity and may be interpreted as 

proportions of end-members; (2) the scores sum to the constant row-sum 

of the original data matrix and, if none are negative, may be interpreted 

as composition values in units of percent or parts per million; and (3) the 

signs of the loadings become fixed and are dependent only on the choice 

of end-members for the model and on the nature of the compositional 

variation in the rocks under examination.

When the number of factors in a Q-mode factor model is less than the 

number of variables in the data matrix, the model accounts for less than 

the total variation in the data and may be used to reproduce the data 

matrix only approximately. The degree of approximation with which the 

normalized form of the data may be reproduced is indicated precisely by 

the eigenvalues of the matrix of coefficients of proportional similarity 

and by the sample communalities. These measures, however, are only 

approximate indicators of the degree to which the derived model may be 

used to reproduce the original data. The proportion of the variance for 

each variable that the model acocunts for is given by the coefficient of 

determination between the original and reproduced data.

The coefficients of determination have been used to construct 

factor-variance diagrams for four suites of igneous-rock samples. The 

diagrams show the proportion of the total variance in each chemical 

constituent that can be explained by petrologic models containing any 

number of factors, or end-members, up to the number of constituents 

present in the samples. Petrologic models have been developed for all 

four suites by choosing reference vectors with composition scores that 

approach the compositions of end-members thought to have been 

involved in the petrologic system. The resultant models are, for the most 

part, in reasonable accord with geologic observations by other workers 

and serve to quantify various aspects of the processes by which the rocks 

may have originated.

The four suites of samples are from (1) a rhyolite-basalt complex in 

Yellowstone National Park, Wyo., (2) a granitoid intrusive in the 

southern part of the Snake Range, Nev., (3) lavas and pumices from the 

1959 summit eruption of Kilauea, Hawaii, and (4) the layered series of the 

Skaergaard intrusion, Greenland.

GENERAL INTRODUCTION

Work of Chayes (1960) has shown that matrices of 

compositional data on rocks and other materials are 

difficult to examine for genetic implications because such 

matrices either have constant row-sums or are parts of 

larger matrices that have constant row-sums. A test 

designed to overcome the difficulty that the constant 

row-sum presents (Chayes and Kruskal, 1966) was found 

to be ineffective (Miesch, 1969). The situation is a serious 

one because petrologists work with this type of data almost 

constantly and the ultimate goal of their work is to 

understand rock genesis. Mathematical geologists have 

known for some time that the constant row-sum property 

of compositional data was no serious hindrance in the 

method of Q-mode factor analysis introduced by Imbrie 

(1963), but it is now apparent that in order to derive a 

Q-mode model in terms of the original data the constant 

row-sum is a definite asset.

Although the method being used is referred to as factor 

analysis throughout this report, this name may be 

unacceptable to many students and workers in multivariate 

statistics. The diagonal values in the similarity matrices 

have been unity in all applications of the method thus far, 

and the method, therefore, might best be referred to by the 

general terms "vector" or "components" analysis, but the 

term "factor" analysis is used because of its prevalence in 

geologic literature.

The purposes of this report are to describe the 

mathematical basis for an extended form of the Q-mode 

method originally developed by Klovan and Imbrie (1971) 

and to demonstrate its application to four problems in 

igneous petrology. The original method leads to a factor 

model that can be used to reproduce the observed data in a

Gl
G2 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

normalized form only. Because of this, the models are 

difficult to interpret in terms of geologic conditions and 

processes. It will be shown here that if the row-sums in the 

original data are constant, these models can be modified to 

reproduce unbiased approximations of the original data in 

units of, for example, weight percent or parts per million. 

This modification further allows the geochemist or 

petrologist to use some other methods to develop a variety 

of models that are mathematically acceptable and to 

examine each of them for geologic plausibility.

(Note. After this report was written, an abbreviated 

description of the extended Q-mode method (Miesch, 

1975) and computer programs that perform most of the 

basic computations (Klovan and Miesch, 1975) were 

prepared for publication in Computers & Geosciences.)

Among the more significant of these other methods is 

one that enables the geochemist or petrologist to propose 

various compositions as those of end-members in a 

geochemical or petrologic system and then to test these 

propositions. If the compositional variation in a basaltic 

magma, for example, is thought to have resulted from the 

crystallization and redistribution of olivine, it is possible 

not only to test this hypothesis but also to determine the 

olivine composition most likely to have been involved. 

Further, for a given parent-magma composition, one can 

determine the amount of olivine that would have had to 

have been separated from or added to the magma to form 

each sample of basalt.

The compositional variation that will be observed in 

most sample suites will have resulted from processes that 

were more complex than the simple separation of one 

mineral from a magma. Other methods made possible by 

the presence of the constant row-sum enable one to 

determine the minimum number of end-members required 

to account for any given proportion of the variance in each 

compositional variable. This can be done before one 

knows what the end-members might be. The methods have 

clear advantages over eigenvalues, which are commonly

used to determine the number of end-members required, 

and are suggested as alternative means for this purpose.

The principal advantage of the extended form of 

Q-mode factor analysis may be that the derived model is 

more easily interpreted; the composition scores 

(end-member compositions) are in the same units as the 

original data, and the composition loadings (mixing or 

unmixing proportions) sum to unity for each sample and 

have signs that are fixed rather than arbitrary. Thus, not 

only can the investigator examine the plausibility of 

proposed geologic processes more easily but also the results 

of his work can be presented in a form that will be more 

understandable and acceptable to nonmathematical 

colleagues.

In part I the mathematical bases of the methods are 

given where the number of factors in the model is equal to 

the number of columns (representing compositional 

variables) in the data matrix. The mathematical effects of 

reducing the number of factors in the model are described 

in part II, and applications of the method to four 

petrologic mixing problems are discussed in part III.

The methods described in this report are intended to 

supplement those described by Klovan and Imbrie (1971). 

Their computer program gives the method for derivation 

of the unsealed factor scores from which the composition 

scores, described here, are computed. Their program also 

provides an efficient method for deriving a Q-mode factor 

solution for large data matrices. Klovan and Imbrie's 

program was adapted to the U.S. Geological Survey's 

system of statistical programs (STATPAC) by George Van 

Trump.

Parts I and II of this report were read by J. J. Connor, 

L. J. Drew, and J. E. Klovan; part III was read by R. W. 

White. I am grateful to each of these reviewers for 

criticisms that led to improvement of the manuscript. In 

addition, I benefited from discussions of the petrologic 

aspects of the work with Fred Barker, D. E. Lee, G. J. 

Neuerburg, and R. E. Wilcox.

Part I. Mathematical Basis Where the Number of Factors 

Equals the Number of Variables

INTRODUCTION

The type of data matrix being considered is one 

wherein the rows represent rock samples and the columns 

represent chemical or mineralogical variables. Where the 

data are expressed as percentages, each row sums to about 

100, depending on whether all the major constituents 

present in the rocks are represented in the data and on the 

nature of the analytical errors. Large departures from 100

may be eliminated by recomputing the row to sum to 100, 

thereby spreading the effect of the analytical errors over all 

the constituents in proportion to their abundances and 

directing the Q-mode analysis at only the part of the rock 

that is composed of the constituents represented in the 

data. If all row-sums are less than 100, the departures may 

also be eliminated by subtracting the row-sum from 100 

and treating the difference as a new compositional 

variable either some major constituent not represented in
Q-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G3

the data or a composite of constituents that may be 

referred to as "others".

As has been customary in most geologic applications of 

Q-mode factor analysis, the matrix of original data can be 

treated in a numblr of different forms. The first type of 

transformation consists of scaling columns (representing 

variables) of the matrix so that the magnitudes of the 

values are approximately the same from one column to 

another. This can be done by either of two options 

provided in the computer program of Klovan and Imbrie 

(1971). The first option is used to express each column as 

a proportion of the maximum value in the column so that 

no values exceed unity; the second option is used to express 

each value as a proportion of the total range for the 

column so that all values within each column are in the 

range from zero to unity. The purpose of this scaling is to 

give each compositional variable an approximately equal 

weight in the factor analysis. Without scaling, in treating 

conventional petrographic variables, for example, 

variables with large variance, such as SiO2» may 

completely dominate the outcome, and variables with 

small variance, such as Na2<D, may exert little or no 

influence. In effect, the use of the scaling options in 

Klovan and Imbrie's program overcomes the problems 

caused by the fact that the means and variances in most 

geochemical and petrographic data are strongly related, 

and it recognizes the fact that minor compositional 

constituents can be just as diagnostic as major constituents 

in arriving at petrogenetic models.

Whether the original data are transformed by scaling or 

not, the data matrix is also treated in a row-normalized 

form. That is, each row of the matrix is adjusted so that 

the sum of the squares of the values within it is unity. This 

adjustment is done automatically when using Imbrie and 

Purdy's (1962) coefficient of proportional similarity 

(cosine theta); in effect, the matrix of cosine theta is 

computed as the major product of the row-normalized 

data matrix and its transpose. The appropriateness of 

cosine theta as a measure of similarity among rock

compositions was shown by Imbrie (1963). Aside from 

providing an appropriate measure of compostional 

similarity, the use of cosine theta allows treatment of each 

row of the data matrix as a vector of unit length, thereby 

facilitating subsequent computation and interpretation.

The conventional forms for a hypothetical data matrix 

are given in table 1. Tables \A and IAA show, 

respectively, the original data with constant row-sums and 

the same data normalized so that the sum of squares for 

each row is unity. The normalization factors, //, which are 

the square roots of the sums of squares for the 

corresponding rows in the original data matrix, are also 

given. Table IB shows the original data transformed so 

that the maximum value in each column is unity, and table 

1C shows the original data transformed so that each 

column ranges from zero to unity. Tables IBS and ICC 

show the normalized forms of the transformed data.

The relationships among the various forms of the data 

are illustrated in figure 1, which shows that the 

row-normalized form of the data, derived from either the 

original or transformed data, is used to derive the factor 

model. Then, the factor model, developed by the 

conventional procedures, can be used to reproduce the 

row-normalized data. The extended method described here 

allows one to modify the model to reproduce the data in 

their original or scaled forms.

Q-MODE FACTOR ANALYSIS

The purpose of Q-mode factor analysis in geochemistry 

and petrology is to resolve a data matrix into a concise 

model wherein each analyzed sample of rock, soil, or other 

material is viewed as a mixture of a small number of 

theoretical or actual end-member samples. As has been 

shown by Imbrie (1963), the method is used to determine 

the number of end-members that are required to account 

for any given proportion of the total variability in the 

normalized data and then to estimate the proportions of 

each end-member present in each sample. Either the

Compositional data in original units of measurementpercent, parts per million, 

and so forth

Data transformed to proportions of the maximum 

values or proportions of 

the ranges for each variable

Row-normalized data Reproduced 

row-normalized data

Reproduced compositional 

data in original units of 

measurement percent, 

parts per million, and 

so forth

Reproduced transformed 

data as proportions of the 

maximum values or proportions of the ranges for 

each variable

FIGURE I. Diagram showing the relationships between the factor model and the various forms of the original and reproduced data.
G4 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

TABLE \.-Original hypothetical data, transformed data, and normalized data used to illustrate the methods of computation

[tj is the normalization factor for the ith row; tj = ( S*/y ) : Xjj = x/j /tj]

A. Original data, x /, = xf

90 7

50 10

* ' =80 15

60 30

85 10

95 3

xmaXj = 95 30

xminj =50 3

3

40

5

10

5

2

40

2

B. Data transformed to proportion of the maximum, j

0.947 0.233

.526 .333

v , .842 .500

«= .632 1.000

.895 .333

1.000 .100

0.075

1.000

.125

.250

.125

.050

t, Row-sum

90.32 100

64.81 100

81.55 100

67.82 100

85.73 100

95.07 100

cy = Xjj/xmaXj

0.978

1.178

.987

1.209

.963

1.006

C. Data transformed to proportion of the range,

x '. = (x^ ; xminj )/(xmaXj xm

0.889 0.148

0 .259

.667 .444

xij = .222 1.000

.778 .259

1.000 0

'";)

0.026

1.000

.079

.211

.079

0

t,

0.902

1.033

.805

1.046

.824

1.000

A A. Normalized original data

0.996 0.078

.772 .154

.981 .184

y ~ .885 .442

.991 .117

.999 .032

0.033

.617

.061

.147

.058

.021

BB. Normalized transformed data from B

0.968 0.238

.447 .283

.853 .506

*V ~ .522 .827

.929 .346

.994 .099

0.077

.849

.127

.207

.130

.050

CC. Normalized transformed data from C

0.986 0.164

0 .251

.828 .552

xij = .212 .956

.944 .315

1.000 0

0.029

.968

.098

.201

.096

0

end-members are theoretical, as in a principal components 

or varimax model, for example, or they are actual samples 

in the data.

Q-mode factor analysis, therefore, can be viewed as an 

attempt to reduce a data matrix into a smaller matrix that 

may facilitate interpretation in terms of petrologic mixing 

(for example, mixing of sediments, magmas, or magma 

and wallrock), or possibly in terms of petrologic mixing 

and chemical processes (such as solution or precipitation of 

material that tends to be of a specific composition). 

Alternatively, Q-mode factor analysis can be used simply 

as a means for classifying samples or for developing a 

concise description of the compositional variation in a 

sample suite. In using the method for classifying samples 

or for developing a concise description of variation, one is 

capitalizing on the fact that chemical and mineralogical 

variables are seldom, if ever, all unrelated to each other. 

The major chemical and mineralogical variations among 

samples can commonly be described more concisely in 

terms of a few actual or theoretical end-members than in 

terms of each one of the compositional variables.

The general Q-mode factor model is

where xij is an approximation of they'th element in the rth 

row of the normalized data matrix, a'ik is an initial factor 

loading for the /th sample on the Ath factor, and/*} is the 

Arth unsealed factor score for the yth variable. The initial 

factor loadings, the a-'k 's, are the conventional loadings 

for the principal components, varimax, or oblique models 

and are those that are provided by widely available 

computer programs. The unsealed factor scores, the/*} 's, 

are provided by the computer program of Klovan and 

Imbrie (1971) for the principal components model and 

varimax model.

The matrix of normalized data, x\j , is N x M in size 

M). The matrix of initial loadings, a'^ , is 

, and the matrix of unsealed factor scores, 

f£ , is m x M. Where m, the number of factors in the 

model (the number of end-members), is set equal to M, the 

number of chemical or mineralogical variables, the loading
g-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G5

matrix is equal to the data matrix in size, and the model 

can be used to reproduce the normalized data matrix 

exactly. In most actual applications, however, m is less 

than M, and the reproduced data, x-j , are only approximations of the original normalized data, x -j .

For purposes of illustration and verification of the 

extended form of the Q-mode method, the data given in 

tables IA, IB, and 1C will be used to generate a variety of 

factor models with m set equal to M. Verification of the 

models, therefore, will consist of using them to reproduce 

the original data exactly, except for round-off errors.

The first step in generating a factor model is to compute 

a matrix of sums of squares and cross products among the 

columns of the normalized original or transformed data 

matrices. Klovan and Imbrie (1971) have shown that (1) the 

eigenvalues of this matrix are equal to the nonzero 

eigenvalues of a matrix of coefficients of proportional 

similarity (cosine theta), and (2) the eigenvectors of the 

same matrix are the unsealed principal component factor 

scores. The coefficient of proportional similarity (cosine 

theta) is a measure of similarity between the compositions 

of two samples, / and p (Imbrie and Purdy, 1962):

Ex'x

IJ PJ

x ,2 

PJ

(2)

The unsealed varimax factor scores are computed as the 

product of the matrix of principal component scores and 

the varimax transformation matrix (Klovan and Imbrie, 

1971).

Other details of the procedures for deriving the initial 

factor loadings are given by Imbrie (1963) and Klovan and 

Imbrie (1971). It will suffice here to point out that the only 

important difference among the principal components, 

varimax, and various oblique models is that the vectors 

representing the samples are described in terms of different 

sets of reference axes. (See Imbrie, 1963.) The initial factor 

loadings, a£, and the unsealed scores, fkj, for the 

principal components and varimax models, and for the 

type of oblique model developed by Imbrie (1963), which 

uses extreme sample vectors as reference axes, are given in 

tables 2, 3, and 4. Those in table 2 were derived using the 

original data without transformation (table I A). Those in 

table 3 were derived from the data transformed to 

proportions of the maximum value in each column (table 

IE). Those in table 4 were derived from the data 

transformed to proportions of the total range for each 

column (table 1C). The products of the 0/£ and fk'j 

matrices for each of the nine models represented are 

exactly equal to the corresponding normalized data 

matrices in tables 1AA, IBB, and ICC. The principal 

component and varimax loadings and scores in tables 2, 3,

and 4 were derived using the computer program given by 

Klovan and Imbrie (1971). The loadings for the oblique 

models were derived using an addition to the program that 

follows the technique described by Imbrie (1963) and 

Manson and Imbrie (1964). Because m equals M in these 

examples, the scores for the oblique models are simply the 

normalized forms of the original and transformed data and 

are taken directly from tables 1AA, IBB, and ICC.

SCALED SCORES AND 

COMPOSITION LOADINGS

Models of the type represented in tables 2, 3, and 4 are 

difficult to interpret because the factor scores, the f£ 's, 

are normalized and, therefore, dimensionless. Also, the 

initial loadings, the a"k 's, do not sum to unity across k. It 

is preferable to express the scores in units conformable to 

those of the transformed data or those of the original data, 

commonly weight percent concentration. Klovan and 

Imbrie (1971) scale the factor scores by multiplying them 

by the square root of M, but this procedure is only for 

purposes of standardization. Alternative, and realistic, 

scaling is possible if the row-sums of the original data 

matrix are constant across all rows, so that

?*// = K> 

j

(3)

where K = 100 if the original data are percentages. 

Similarly, factor scores in units of the original data should 

sum to the same constant, inasmuch as the scores should 

represent the compositions of real or theoretical endmember samples:

= K, ' <3a>

The transformed data may be regarded as having been 

derived from a general expression found useful for 

computer programming:

- xmnin} ) / (xmax- irij ), (4)

where xmaxj and xmirij are the maximum and minimum 

values in they'th column of the original data matrix if the 

data were transformed to proportions of the total range. If 

transformed to proportions of the maximum value, all 

values of xmirij are set equal to zero, and, if not 

transformed, all values of xmirij are zero and all values of 

xmaXj are unity. Similarly, the transformed scores can be 

defined as

(5)

and this can be rearranged to

fk = fkj + (6)
G6 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

TABLE 2. Initial loadings, a £ , unsealed scores, f k'- , scale factors, sk , and the estimated normalization factors, i\, for the three-factor principal 

components, varimax, and oblique factor models derived without transforming the data

(The products of the matrices of a,-£ and/ty are equal to the normalized data in table IAA]

Model

Principal 

components ........

Varimax .............

Oblique ..............

I nitial loadings, a('^

0.988

.871

.996 a = 

<k .960

.994

.981

0.851

.416

a " - ' 797 

* .608

.826

.871

0.901

o

.800

1.000

-0.138

.481 

-.087 

.066 

-.105

-.159

0.420

.873 

.443 

.493 

.442

.409

-0.004

1.000

-.012 

0 

.020

0

-0.068

-.095 

.031 

.273 

-.035

-.110

0.316

.255 

.410 

.623 

.349

.274

0.112

0

.373 

1.000 

.200

0

Unsealed scores, fkj

0.973

f* = -.186 kj -.136

0.888 

/" = .388 

kj .247

0.999

V" = -772 

k} .885

0.172 

.199 

.965

-0.280 

.030 

.960

0.032

.154 

.442

0.152 

.962 

-.226

-0.365 

.921 

-.136

0.021

.617 

.147

Scale factors, sk

77.053 

s k = 102.499 

165.778

412.169 

sk = 74.647 

93.393

95.068

s = 64.807 

67.823

Estimated normalization 

factors' JJ

90.4

64.8 f. = 81.5

' 67.8 

85.7

95.1

90.3

64.8 

- 81.6 

' ~ 67.8 

85.7

95.0

90.4

64.8

'' = 67.8 

85.7

95.1

'*} = 1/2 <'.«$/**>.

' k

2 From table IAA.

TABLE 3. Initial loadings, a jk , unsealed scores, fkj, scale factors, s k , and the estimated normalization factors, t,-, for the three-factor principal 

components, varimax, and oblique factor models derived after transforming each variable to proportions of the maximum value

[The products of the a£ and/k'j matrices are equal to the normalized transformed data in table 1 BB\

Model

Principal 

components ........

Varimax .............

Oblique ..............

Initial loadings, afk

0.969

.706 

.990 

a!k = -851 

.991

.927

0.904

.261 

 .749 

a ik = .371 

.843

.951

0.880

0

a"= '578 

* 0

.754

1.000

-0.234

.627 

-.037 

.294 

-.131

-.324

0.256

.927 

.300 

.339 

.307

.226

-0.006

1.000

-.019 

0

.032

0

-0.084

-.328 

.135 

.435 

-.015

-.189

0.343

.268 

.591 

.864 

.441

.211

0.185

0

.549 

1.000

.317

0

Unsealed scores, /jy

0.880 

/" = -.407 

* -.245

0.978 

/". = -173 

W .117

0.994

2f" = -447 

kj .552

0.416 

.408 

.813

-0.128 

.055 

.990

0.099

.283 

.827

0.231 

.817 

-.528

-0.164 

.983 

-.076

0.050

.849 

.207

Scale factors, s^

0.950 

s k = 16.001 

-4.983

1.212 

s = 1.742 

* 2.650

1.006

sk = 1.178 

1.209

Estimated normalization 

factors ', Tf

0.978

1.179 

A _ .987 

' 1.209 

.963

1.007

0.978

1.178 

7 _ -987 

' 1.210 

.963

1.006

0.978

1.178

- .988 

' 1.209

.963

1.006

2 From table Iflfl.

Substitution of the right side of equation 6 in equation 3a 

gives:

2 (fkj (xmaxj-xmirij) + xmin^ = K. (7)

Although equations 5 and 6 define the relation between the 

values of fkj and fkj , their absolute values have not been 

established. This is done by defining fkj as equal to skf'k'j , 

where sk is a scale factor and fkj is an unsealed score as
Q-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G7

TABLE 4. Initial loadings, a-^, unsealed scores, f(], scale factors, sk , and estimated normalization factors,^;, for the three-factor principal 

components, varimax, and oblique factor models derived after transforming each variable to proportions of the range

[The products of the a// and/^'matrices are equal to the normalized transformed data in table ICC]

Model

Principal

components ........

Varimax .............

Oblique ..............

Initial loadings, a$

0.963

.234

.985 

°ik .603

.995

.906

0.984

.022

.822 

°ik .204

.942

.999

0.949

0

n" - -704 a it 

1 0

.876

1.000

-0.249

.830

.094 

.643

-.095

-.375

0.012

.979

.105 

.244

.088

-.026

-0.006

1.000

-.020 

0

.032

0

0.101

.506

-.147 

-.473

.037

.1%

0.178

.202

.559 

.948

.324

.016

0.173

0

.583 

1.000

.321

0

Unsealed scores, fkj

0.906

f'k'j =-374 

.196

1.000

/". = -.025 

.016

1.000

.212

0.400

.608 

-.686

-0.014

.051 

.999

0

.251 

.956

0.138

.700 

.701

0.026

.998 

-.050

0

.968 

.201

Scale factors, sk

0.792

sk = 1.720 

2.656

0.987

sk = 1.179 

1.747

1.000

sk = 1.033 

1 1 .WTU 046

Estimated normalization 

factors ',"//

0.902

1.032

.805 

tf 1.045

.823

1.000

0.902

1.033

.805 

'' ~ 1.046

.823

1.001

0.902

1.033

.805 

'' ~ 1.046

.824

1.000

2 From table ICC.

previously defined. Equation 7 may now be written as

kj (xmax j-xmin j) + = K. (8)

The scale factors for converting the unsealed scores to 

scores that are conformable to the transformed data may 

be derived by rearrangement of equation 8 to

K-

/'.' (xmax- xmin. )\

kj j j /

(9)

where, again, xmaXj and xmirij are the maximum and 

minimum values in the y'th column of the original data 

matrix if the data were transformed to proportions of the 

total range. If the transformation was to proportions of 

the maximum value, then all values of xmirij used in 

equation 9 are set equal to zero, and if no transformation 

was made, all values ofxmirij are set equal to zero, and all 

values of xmaxj are set equal to unity.

The Q-mode factor analysis model in equation 1 can be 

rewritten as

aik

or

(10)

(11)

where a-k = o-k /sk and f-k is a scaled factor score. 

Equation 11, like equation 1, may be used to derive an 

estimate of the normalized form of the original or 

transformed data matrix or may be used to derive the 

normalized form of the data exactly, for models where m 

= M. However, the model may be changed to one that 

reproduces the transformed data (or the original data if the 

data were not transformed) if both sides are divided by the 

quantity So A. giving K '*

Although not proven here, it has been observed and 

abundantly verified in the present study that when m = M 

the quantity \/^a-k is, except for round-off errors, 

precisely equal to the normalization factor, t( , which is the 

square root of the sum of squares of the rth row in the 

original or transformed data matrix, depending on whether 

or not a transformation of the data was made. This can be 

seen by comparing the values of ?,- in tables 2, 3, and 4 with 

the corresponding values of t\ in table 1. Therefore, the 

left side of equation 12 may be written simply as x\j , 

which is the approximation of an element in the 

transformed data matrix or in the original data matrix if no 

transformation was made. The quantity a'ik /'Eai'l( in 

equation 12 may be written as aik : a matrix analogous to 

that of aik is referred to as the composition matrix by 

Imbrie and VanAndel (1964, p. 1141) where it is derived
G8 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

TABLE 5. Composition loadings, a ik , and composition scores, fkj = fkj, for the three-factor principal components, varimax, and oblique factor

models derived without transforming the data

[The products of the a^ and/^.. matrices are equal to the original data in table \A\

Model

Principal

components ........

Varimax .............

Oblique ..............

Composition loadings, atf.

1.159

.732 

a 1.054 u iir 

lk .854

1.106

1.211

0.186

.065

aik = -158 

.100

.172

.201

0.856

0

0.567

Oik =

'* 0

.722

1.000

-0.122

.304 

-.069

.044

-.088

-.148

0.508

.758

.484 

.448

.508

.521

-0.005

1.000

-.016

0

.026

0

-0.037

-.037 

.015

.112

-.018

-.063

0.306

.177

.358 

.452

.320

:279

0.150

0

.449

1.000

.253

0

Composition scores,/^- = /.'.

74.99

fkj = fkj = -19.04

-22.49

365.93 

fkj = / ' = 29.00 

' kj 23.04

95

'4 = fkj = 50

60

13.29

20.41

159.92

-115.35 

2.25 

89.62

3

10

30

11.72

98.63

-37.43

-150.58 

68.75 

-12.66

2

40

10

1 From table 1/1.

for an oblique model. The quantity aik will be referred to

here as a composition loading whether it pertains to an 

oblique model, a varimax model, or a principal 

components model. Equation 12 then becomes

X;; = (13)

where x\. refers to the transformed data if a 

transformation was made, or to the original data if there 

was no transformation. The m values of a^k sum to unity 

across k and may be interpreted as proportions. Table 5 

gives values of aik and /fa for the principal components, 

varimax, and oblique models derived without transformation of the data. The products of the corresponding aik 

and/^ matrices are equal to the original data matrix in 

table \A. The matrices of a ik and fkj in table 6 were 

derived after transforming the data to proportions of the 

maximum values, and their products are equal to the 

transformed data matrix in table \B. The matrices of ajk 

and/ky in table 7 were derived after transforming the data 

to proportions of the total ranges for each variable, and 

their products are equal to the transformed data matrix in 

table 1C.

COMPOSITION SCORES

If the data were transformed prior to the derivation 

of the OIK and fkj matrices, the scaled scores must be 

converted to composition scores in order to reproduce the 

original data, or an approximation of the original data,

from the factor model. Each side of equation 13 is 

multiplied by the quantity (xmaXj- xmirij) and increased 

by xmin to obtain

Xjj(xmaXj -xmirij) + xmitij =

^(aikfkj (xmaXj -xmin.)) + xmin . (14)

By virtue of the fact that aik = a-k /2tf^ , it follows that

Therefore, equation 14 may be written as

x'ij(xmaXj xminj) + xmin- =

^i(aikfkj (xmaxj-xminj)) + J^a

and

(15)

+ xmn y = 

fli* (fkj (xmaXj-xminj) + xmin^. (16)

The effect of the values ofxmaXj and xminj in equation 16 

is to reverse the transformation that consisted of scaling 

each variable to a proportion of the total range for the 

variable. Or, if all values of xmin. are set equal to zero, the 

effect is to reverse the transformation that consisted of 

scaling each variable to a proportion of the maximum
2-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G9

TABLE 6. Composition loadings, a jk , scaled scores, fkj, and composition scores, fy , for the three-factor principal components, varimax, 

and oblique factor models derived after transforming each variable to proportions of the maximum value

[The products of the a,-^ and ffa matrices are equal to the transformed data in table Ifl. The products of the a,-^ and /£.- matrices are equal to the original

data in table LA]

Model

Principal 

components ........

Varimax .............

Oblique ..............

Composition loadings

0.998

.876

1.029 

fl*~ 1.083 

1.005

.982

0.730

.254

.610 

fl*~ .370

.670

.790

0.856

0

0.567 

a* 0

0.722

1.000

-0.014

.046

-.002 

.022 

-.008

-.020

0.144

.627

.170 

.235

.170

.131

-0.005

1.000

-.016 

0

.026

0

"ik

0.016

.078

-.027 

-.106 

.003

.038

0.127

.119

.220 

.394

.160

.080

0.150

0

.449 

1.000

.253

0

Scaled scores,/^

0.836

f'kj =-6.515 

1.223

1.185

4'. = -301 

1 .309

1.000

'/' = .526 J kj 

.632

0.395

6.533 

-4.050

-0.155

.096 

2.624

0.100

.333 

1.000

0.219

13.074 

2.633

-0.199

1.714 

-.202

0.050

1.000 

.250

Composition scores, /^y

79.39

A/ =-618.92 

116.18

112.62

fkj = 28.57 

29.34

95

2fkj= 50 

60

11.84

195.98 

-121.51

-4.65

2.89 

78.73

3

10 

30

8.77

522.94 

105.33

-7.97

68.54 

-8.07

2

40 

10

1 From table Ifl.

2 From table \A.

TABLE 7 . Composition loadings, a jk , scaled scores, fk'., and composition scores, fkj , for the three-factor principal components, varimax, and 

oblique factor models derived after transforming each variable to proportions of the range

[The products of the a,-^ and ffi matrices are equal to the transformed data in table 1C. The products of the a,-£ and fy matrices are equal to the original

data in table \A]

Model

Principal

components ........

Varimax .............

Oblique ..............

Composition loadings

1.096

.304 

1.001

aik = .795

1.034

1.144

0.899

.023

a - - 671 

ik .216

.786

1.012

0.856

0

a - °'567 

ik 0

.722

1.000

-0.131

.499 

.044

.391

-.045

-.217

0.009

.858

.072 

.216

.061

-.021

-0.005

1.000

-.016 

0

.026

0

"ik

0.034

.197 

-.045

-.186

.011

.074

0.092

.119

.258 

.567

.153

.009

0.150

0

.449 

1.000

.253

0

Scaled scores, f/fi

0.718

fkj = -.644

.522

0.987 

fkj = -.030

.027

1.000 

l fk'j= 0

.222

0.317

1.045

-1.822

-0.014 

.060

1.745

0 

.259

1.000

0.109

1.204

1.861

0.026 

1.177

-.088

0 

1.000

.211

Composition scores,/^;

82.31

fkj= 21.02

73.47

94.40 

A/= 48.65

51.24

95

A/ = 50

60

11.55

31.23

-46.20

2.62 

4.62

50.11

3 

10

30

6.14

47.75

72.72

2.98 

46.73

-1.35

2 

40

10

1 From table 1 C.

2 From table 1/4.

value for the variable. In either case, because of equations 

4 and 6, we can now write:

*» = ?«»/«. (17>

where x^ refers to the original data in table L4, the a/# *s 

are composition loadings, and the fy 's are composition 

scores for the principal-components, varimax, or oblique 

factor models.
G10 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

The composition scores are in units of the original data 

and sum to K across j (eq 3a). For the oblique model 

developed by Imbrie (1963), the composition scores, when 

m = M, are simply the compositions of the samples 

selected to serve as end-members in the factor solution. For 

the principle component and varimax models, the derived 

composition scores for the kih factor may be regarded as 

the composition of a theoretical sample that may serve as 

an end-member for the corresponding model. However, if 

any value of fy is negative or in excess of K, the derived 

composition of the kih factor is, of course, unrealistic and 

should be viewed as only an approximate indication of the 

nature of some theoretical end-member sample.

The composition scores, fk - , derived from the data in 

table 1 are given in tables 5, 6, and 7 for the principal 

components, varimax, and oblique models. Those in table 

5 were derived without transformation of the data. Those 

in table 6 were derived after the data were transformed to 

proportions of the maximum values, and those in table 7 

were derived after the data were transformed to 

proportions of the total range for each variable. The 

matrices of composition loadings, a/A. , in tables 5, 6, and 7 

may be multiplied by the corresponding matrices of 

composition scores, fk . t to produce the matrix of original 

data in table IA.

NEGATIVE LOADINGS AND 

NEGATIVE SCORES

In most factor analysis applications the signs of the 

factor loadings are regarded as arbitrary, and it is entirely 

valid to change all the signs in any column of the initial 

loading matrix. However, if factor scores are computed 

and the signs are changed in the kih column of the loading 

matrix, it is necessary to also change the signs in the kth 

row of the score matrix. When these sign changes are made 

for the composition loadings and scores, the scores will 

sum to minus K, rather than K, and the requirement in 

equation 3a is violated. This situation is not corrected by 

changing the signs in the initial loadings and unsealed 

scores from which the composition loadings and scores 

were derived. When the sign changes are made for the 

initial loadings and unsealed scores, the scale factor, sk , 

also changes in sign, as may be seen from equation 9. 

Because of this, the sign change is cancelled out (eq 10) 

before the composition loadings and scores are derived. 

Therefore, unlike the conventional (initial) loadings, the 

composition loadings and scores have fixed signs that are 

determined by the choice of reference axes and the nature 

of the compositional variation in the rocks under 

examination.

In the use of Q-mode factor analysis to develop a genetic 

model, the composition scores for the selected reference 

vectors should be entirely nonnegative because these scores 

represent chemical or mineralogic compositions of

end-members. However, if the factor solution is to be used 

as a device for summarizing geochemical or petrologic data 

or for purposes of sample classification, negative 

composition score values can be perfectly acceptable. I 

have found that the composition scores for the 

conventional varimax reference axes are commonly 

negative in part but that the varimax axes can still be useful 

for reference purposes and the set of scores for each axis is 

still indicative of the general compositional nature of the 

theoretical end-member. For example, a series of scores 

consisting of some value larger than the constant K for 

SiC>2, and of negative values for all other constituents 

indicates that a more realistic end-member might have the 

composition of a siliceous magma or a pure quartz sand.

In situations where negative composition scores are 

unacceptable, they can be avoided by the choice of a 

different set of reference axes. The method for finding the 

composition scores for any vector that may serve as a 

potential reference axis within the factor space of the 

Q-mode solution is demonstrated in the following section. 

After m reference axes with all nonnegative composition 

scores have been selected, there will be no difficulty in 

determining the composition loadings for each sample with 

respect to them as long as the selected axes are sufficiently 

independent. (That is, they are not colinear or 

approximately colinear in a two-factor space, coplanar in a 

three-factor space, and so forth.)

The composition loadings, aik , always sum to unity 

across k and their signs indicate the required addition or 

subtraction of the kih end-member to form the fth sample 

according to the model. Where either the addition of an 

end-member, as indicated by a positive loading, or the 

subtraction of an end-member, as indicated by a negative 

loading, is unreasonable, a new choice of reference axes 

may be desirable.

Composition scores for the first axis of a principal 

component model tend to be all nonnegative, as do the 

composition loadings on the first end-member. The scores 

and loadings for subsequent axes tend to be positive and 

negative in about equal number.

Although some of the composition scores for varimax 

models are commonly negative, the composition loadings 

are generally almost entirely positive. As the varimax axes 

are changed to oblique axes by moving them toward the 

vectors representing the actual samples, their composition 

scores become increasingly positive and the composition 

loadings become increasingly negative. The first of these 

properties can be useful in the modification of varimax 

models to form models with oblique reference axes that 

have all nonnegative composition scores, as described in 

part II.

Oblique models that employ the extreme sample vectors 

as reference axes, as used by Imbrie (1963), have all 

nonnegative composition scores if m = M, or if for some
6-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS Gil

other reason the model accounts for all the variability in 

the data. Where this is not true, the composition scores 

differ from the actual compositions of the reference 

samples and some may be negative. The composition 

loadings for this type of oblique model can be entirely 

nonnegative, but commonly some of them are slightly 

negative.

The signs of the composition scores and loadings for 

other types of oblique models depend entirely on the 

reference axes chosen to describe a particular set of sample 

data. These other types of oblique models may be 

developed by choosing alternative oblique reference axes. 

The choice of reference axes may be based either on 

examinations of the composition scores for selected vectors 

or on the determination of the vector representations for 

given compositions.

FINDING THE COMPOSITION

SCORES FOR ANY VECTOR IN

THE FACTOR SPACE

The oblique model developed by Imbrie (1963) is one 

wherein actual samples of extreme composition serve as 

end-members, or reference samples, for describing the 

compositional variations within the entire sample series. In 

many studies, however, there is no reason to expect that 

the actual end-members which contributed to the 

formation of a rock body are represented in the sample 

data. Thus, the extreme samples might be rejected as 

plausible end-members, and a search begun for alternative 

real or hypothetical end-members that have acceptable 

composition scores as well as acceptable composition 

loadings of all samples with respect to them. It is an easy 

matter to determine the composition scores for any vector 

that may potentially serve as a reference axis.

As an example, we may consider an arbitrary vector 

within the three-factor space of the varimax solution 

represented in tables 4 and 7. The initial loadings of the 

vector on the varimax axes are taken as

a/,' = 0.03162, a/2 = 0.03162, a"3 = 0.9990. (18)

These loadings (which define a vector close to but not 

coincident with the third varimax axis) and the unsealed 

varimax factor scores from table 4 are used to derive the 

normalized composition of the vector according to 

equation 1. Because the normalized composition is to be 

regarded as a set of unsealed scores, the notation is 

appropriately changed from x -j toffi . The derived values

are

; = 0.047, //I =0.999, f' =-0.018. (19) 'k2

Equation 9, then, gives sk = 1.586, and the scaled factor 

scores for the vector are

//, = 0.074, /' =1.584, /'=-0.029, (20) *3

and, from equation 6, the composition scores for the 

vector are

fkl = 53.33, fa = 45.77, 4,= 0.90. (21)

Because none of the composition scores is negative, the 

vector may serve as a more acceptable reference axis than 

the third varimax axis which has a negative value for /33 

(table 7).

Composition scores from some other arbitrary vectors 

within the same factor space are given in table 8. Note that 

the scores for colinear vectors having loadings that are 

opposite in sign are exactly the same. Note also that some 

vectors, such as vectors 3 and 4, may have composition 

score values that are negative or in excess of K and 

impossible to interpret in terms of rock composition.

TABLE 8. Computed composition scores, f.., for some arbitrary vectors 

in the varimax solution of tables 4 and 7

Loadings, Composition scores,/;;

*=1 * = 2 * = 3 j=\ j=2 j=3

1

2

3

4

5

6

7

0.577

-0.577

.271

-.271

.924

.271

.382

0.577

-.577

-.924

.924

.271

-.924

.654

0.577

-.577

.271

-.271

.271

-.271

.654

68.31

68.31

12.06

12.06

81.38

30.38

62.85

14.49

14.49

-13.77

-13.77

8.55

16.07

16.97

17.20

17.20

101.70

101.70

10.08

53.55

20.18

FINDING THE VECTOR

REPRESENTATION 

FOR A GIVEN COMPOSITION

Although the search for possible end-members in a 

compositional system may consist of the determination of 

composition scores for various vectors in the factor space, 

as described in the previous section, it is commonly more 

satisfactory to begin with the composition, Zj, of some 

material that may have been involved in the system and 

determine whether or not the composition can be 

satisfactorily represented as a vector in the m-dimensional 

factor space. Any composition that can be represented in 

the factor space can be arrived at by the mixing or 

unmixing of compositions equal to the composition scores 

of any other m vectors in the same space. Any composition 

that cannot be represented in the factor space could not 

have been an end-member in the w-component 

compositional system.

The procedure cannot be fully illustrated here because 

the example factor solutions have m = M, and any 

composition of the same M components can be completely 

represented in the m-dimensional space. However, the 

procedure begins with the transformation and normalization of the given composition in the same manner used 

to transform the original data before derivation of the
G12 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

factor solution, using the same values of xmaxj and xmin} 

(table 1). Representing the transformed and normalized 

composition by a row of zj, the relation given in equation 

1 can be written as

Equation 22 in matrix notation is 

Z = AF.

(22)

(23)

Following Klovan and Imbrie (1971, p. 64), both sides of 

equation 23 are postmultiplied by the transpose of the 

unsealed score matrix to give

ZF' = AFF'. (24)

However, the unsealed scores are in normalized form, and, 

if they are for orthogonal reference axes, such as the 

principal component or varimax axes, the product of F and 

F' is an identity matrix. Consequently, equation 24 

becomes

A = ZF'.

This translates to

(25)

(26)

which shows that the initial loadings on the principal 

components or varimax axes are derived simply by 

multiplying the array of normalized composition values by 

the transposed matrix of unsealed scores.

The sum of squares of the m values of a'k is the 

communality h 2 of the vector representing the given 

composition and indicates how well the composition fits 

into the factor space of the sample data and the degree to 

which the given composition in original units (generally 

percent or parts per million) may be reproduced from the 

factor model. To reproduce the original composition, the 

derived values of ak are divided by the scale factors, sk , for 

the principal components or varimax axes, and adjusted to 

sum to unity across k, giving ak . The original composition, 

Z: , is then estimated by

£ ^ ~ f fyi\ Z- =2jOkfkj , (£r)

where the values of fkj are the composition scores for the 

principal components or varimax axes. Wherever h 2 is 

unity, the quantities (z-z)j will be zero for each variable, 

and, where h 2 is less than unity, the absolute values of (z-z)j 

will be correspondingly larger.

The vector representations of some arbitrarily selected 

compositions in the varimax model of tables 4 and 7 are 

given in table 9. Because m = Min this model, the value of 

h 2 is unity for each composition.

In real applications of the method, where m is less than 

M, the acceptability of the given composition z , as that of 

a possible end-member in a petrologic system will depend 

on its derived communality, h 2.

TABLE 9. Vector loadings, a^, in the varimax solution of tables 4 

and 7for some given compositions, z,j

Phase, /'

Composition, Zj; Vector loadings,

7=1 7=2 7 = 3 k=\

1

2

3

4

5

6

7

63

75

99

0

12

14

6

24

3

0

82

5

76

10

13

22

1

18

83

10

84

0.325

.744

.995

-.361

-.344

-.294

-.388

0.366

.668

-.054

.189

.939

.130

.920

0.873

-.023

-.084

.913

-.020

.947

.057

COMPUTATION OF INITIAL

AND COMPOSITION LOADINGS

WITH RESPECT TO OBLIQUE

REFERENCE VECTORS

Determination of the initial factor loadings with respect 

to any oblique reference axes follows the technique 

described by Imbrie (1963). It is illustrated here using the 

varimax model represented in tables 4 and 7. For the 

purpose of illustration only, the oblique reference axes will 

be taken as vectors 5, 6, and 7 in table 8; the loading matrix 

of the new reference axis system is

0.924 0.271 0.271

A=a-'k = .271 -.924 -.271

.382 .654 .654. (28)

The inverse of A is

1.306 0.0 -0.541

A' 1 = .859 -1.531 -.990

-1.621 1.531 2.836. (28a)

The product of the matrix of initial varimax loadings in 

table 4 and the inverse matrix in 28a gives the matrix of 

initial loadings on the new reference axes

a., = ik

1.007

.542

.258

-1.061

.781

1.256

0.254

-1.190

.695

1.078 

.361 

.064

-0.039

-.408

1.037

2.337

.322

-.469. (29)

The matrix of unsealed factor scores, fk'., for the new 

reference axes is then determined as the product of the 

matrix of initial loadings, a'ik , in 28 and the matrix of 

unsealed varimax scores, fkj , in table 4. The scale factors, 

sk , for the new reference axes are then computed, using the 

new scores and the values of xmaxj andxmirjj from table 1 

in equation 9. The scale factors are

s = 0.7562, s2 =-1.5040, *3 =0.7595. (30)

The composition loadings for the new model are derived by 

dividing the columns of 29 by the corresponding values of
0-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G13

Sk in 30 and then adjusting each row to sum to unity. The 

composition loadings are

1.199 -0.152 -0.047

.739 .816 -.556

aik = .274 -.372 1.098

-1.465 -.749 3.214

.849 -.198 .349

1.661 -.043 -.619. (31)

The original data in table \A can be reproduced exactly, 

except for differences due to round off, by multiplying the 

composition loading matrix in 31 by the matrix of 

composition scores for vectors 5, 6, and 7 given in table 8.

SUMMARY

Q-mode factor analysis appears to be a method that is 

well suited for exploring the nature of the compositional 

variation in a suite of rock samples. It can be used in at 

least three ways: (1) to develop a geochemical or petrologic 

model, (2) to develop a scheme of classification, or (3) to 

summarize complex multivariate data. Applications of the 

second and third types can employ any set of reference 

axes, regardless of their composition scores. The 

development of tenable geochemical or petrologic models, 

however, requires that all composition scores be in the

range from zero to K, and the acceptability of the model 

may depend on the consistency of the composition 

loadings, especially their signs, with other geologic 

observations. These conditions will usually require 

experimentation using various sets of reference vectors. 

Potential reference vectors can be examined by either 

finding the composition scores for selected vectors or 

finding whether or not certain compositions of interest can 

be satisfactorily represented as vectors in the factor space.

The ultimate test of the mathematical validity of any 

factor model is the closeness with which it can reproduce 

the original data. The only test of its geologic validity is 

whether or not it conforms to other geologic observations 

regarding the rock unit, or units, from which the series of 

rock samples was taken.

All the models developed here can be used to reproduce 

the normalized, transformed, or original data exactly 

because the number of factors they contain was made 

equal to the number of variables represented in the data 

matrix. In most real applications, m is less than M because 

the purpose of the analysis is to reduce the problem, and it 

cannot be expected that the data will be reproduced 

exactly. However, the closeness with which the data can be 

reproduced may or may not be satisfactory to the 

investigator, and the test is recommended.

Part II. Effects of Reducing the Number of Factors

INTRODUCTION

Part I has shown that the constant row-sum of 

matrices of compositional data on rocks and other 

materials can be used to advantage in Q-mode factor 

analysis. Where the matrix row-sums are constant, it is 

possible to scale the factor loadings and scores, which 

together comprise the factor model, so that they are 

conformable with the original data, which are commonly 

in units of weight percent or parts per million. Without 

such scaling, the loadings are difficult to interpret in terms 

of proportions of end-members, and the scores cannot 

easily be interpreted in terms of rock compositions. 

Moreover, the products of the matrices of loadings and 

scores are not approximations of the original data in 

percent or parts per million but of the data in their 

normalized form. If the loadings are scaled and adjusted to 

sum to unity for each sample, and, if the scores are scaled 

and detransformed (if a data transformation was used in 

their derivation), the products of the matrices of loadings 

and scores approach or equal the matrix of original data. 

Examples of factor models were given in part I for which 

the number of factors was equal to the number of variables

(columns) in the data matrix, and the products of the 

matrices of composition loadings and composition scores 

equaled the original data exactly, except for differences 

due to round off. The purpose here in part II is to show 

how the method and the mathematical basis are affected 

when the number of factors is reduced to something less 

than the number of variables. Thus, the computations and 

the models to be developed are more representative of most 

real factor analysis applications wherein the objective is to 

reduce, or simplify, the problem by developing a less 

complex model that explains some major part of the 

compositional variability in a suite of rock samples. The 

minor part of the variability left unexplained is ascribed to 

measurement errors and to the effects of geologic processes 

that have had only minor effects on variation in rock 

composition. The same hypothetical data used in part I 

(table 1) will be used here, but the models will be developed 

with two factors rather than three. The methods will be 

applied to real data in part III.

The initial factor loadings and unsealed scores for the 

principal components and varimax models developed here, 

as those in part I, were derived by the method of Klovan 

and Imbrie (1971).
G14 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

ESTIMATES OF FACTOR LOADINGS

The geometric representation of a Q-mode factor 

model is a system of vectors with a common origin. Where 

the number of factors, m, is equal to the number of 

variables, M, the vectors are of unit length, and each 

represents a row of either the normalized original data 

matrix or the normalized matrix of transformed data 

where a transformation was used. The cosine of the angle 

between any two of the vectors is exactly equal to the 

coefficient of proportional similarity (cosine theta) for the 

two samples that the vectors represent. Because the 

minimum value of cosine theta for nonnegative data is 

zero, none of the sample vectors is separated by more than 

90 degrees from the others. The positions of the vectors 

may be described in terms of any reference axis system, but 

the conventional systems are the principal components 

axes, the varimax axes, and the oblique axes as used by 

Imbrie (1963).

The principal components axes are positioned in the 

vector system so that the sum of squares of the projections 

of the sample vectors on one of the axes is as high as 

possible. In other words, one of the principal axes is placed 

near the center of the vector cluster, and all the sample 

vectors have positive projections on it. The second 

principal axis is positioned orthogonal to the first but is 

oriented so that the sample vectors have as large a 

projection on it as possible considering the orthogonality 

constraint. Subsequent axes are positioned in the 

M-dimensional space so that they are orthogonal to the 

preceding axes and so that the projections of the sample 

vectors on each of them are as large as possible. 

Consequently, the first column of the initial factor loading 

matrix for the principal-components model, which gives 

the loadings on the first axis, contains relatively large 

positive values, and subsequent columns contain both 

positive and negative values of decreasing absolute 

magnitude (tables 2, 3, and 4).

The sum of the squares of the fcth column of the matrix 

of initial loadings on the principal components axes is 

equal to the kth eigenvalue of the cosine theta matrix 

(Harman, 1967, p. 167). Consequently, the eigenvalues 

indicate the number of dimensions occupied by the cluster 

of sample vectors within the M-dimensional space. For 

example, if the first two eigenvalues of the cosine theta 

matrix are large and subsequent eigenvalues are small, the 

sample vectors cluster about a plane. This has been taken 

to indicate that a major part of the compositional variation 

can be described in terms of two end-members and as 

justification for reducing the number of factors, m, to two 

so that m is less than M. Interpretation of the eigenvalues is 

discussed further in the following section.

In the derivation of the varimax model, where m = M, 

the principal components axes are rotated simultaneously, 

while maintaining their mutual orthogonality, so that the 

total sum of the squares of the projections of the sample 

vectors on all the axes is as large as possible. This is the 

varimax criterion. In the derivation of the oblique model, 

where m - M, using the method of Imbrie (1963), the 

reference axes are taken as the M sample vectors that occur 

at the extremes of the sample vector cluster and, if 

possible, enclose all the others.

If the eigenvalues of the cosine theta matrix indicate that 

the cluster of sample vectors occurs largely in m 

dimensions, where m is less than M, only the first 

m-principal-components axes are rotated in deriving the 

varimax model, and the sample vectors are then projected 

into m-dimensional space. Following this projection, the 

sample vectors tend to be of less than unit length. The 

squares of the vector lengths, computed as the sums of the 

squares of the initial loadings, are referred to as the sample 

communalities and indicate the degree to which each 

sample conforms with the final factor model.

These procedures are illustrated here using the same 

hypothetical data (table 1) used in part I, wherein all 

models were such that the number of factors equalled the 

number of variables. As before, the data matrix is treated 

in three ways: (1) the original data (table L4), (2) the data 

transformed so that each variable is expressed as a 

proportion of the maximum value for that variable (table 

IB), and (3) the data transformed so that each variable is 

expressed as a proportion of the total range for that 

variable (table 1C). The row-normalized forms of the data 

are given, respectively, in tables IAA, IBB, and ICC. The 

sample vector cluster previously defined for the principal 

components model is projected into two-dimensional space 

merely by eliminating the third column of the initial factor 

loading matrix. (Compare tables 10,11, and 12 with tables 

2,3, and 4 of part I.) The third row of each unsealed factor 

score matrix is also eliminated. Only two of the principal 

components axes are rotated in deriving the varimax 

model. The projections of the sample vectors on the two 

varimax axes are given in tables 10, 11, and 12. 

Computation of the initial loadings and unsealed scores for 

both the principal components and the varimax models 

follows the methods given by Klovan and Imbrie (1971).

According to the method of deriving the oblique factor 

model described by Imbrie (1963), the m end-member 

samples are first chosen as those that have the largest 

projections on each of the varimax reference axes. As 

shown in part I, the m by m matrix of initial varimax 

loadings, A, is then inverted to give/4' 1 , and the entireNby 

m matrix of initial loadings postmultiplied by A1 gives the
TABLE

Q-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G15

10. Initial loadings, a^ , unsealed scores, fy, scale factors, s k , estimated normalization factors, Tj, and reproduced data matrix, Xjj, for 

the two-factor principal components, varimax, and oblique factor models derived without transforming the data

[The products of the matrices of a-f. and/j.y are approximations of the normalized data in table \AA]

Model

Principal 

components 2

Varimax .....

Oblique ......

Initial loadings,

0.988

.871

.996

 'i = - 960 

.994

.981

0.891

.443

.869

a :if = TC-J

111 .IjL

.877

.898

0.975

0

.908

a ik = .663

.933

1.000

"Ik

-0.138

.481

-.087 

.066 

-.105

-.159

0.448

.891

.495

.600

.479

.427

0.036

1.000

.120

.355

.091

0

Unsealed scores, f£j

0.973 0.172 0.152 

3V ~ -.186 .199 .962

f" _ 0.906 0.029 -0.422

V .400 .262 .878

s 0.98 0.14 0.00

J kj ~ .76 .25 .60

Scale factors, fy

s - k 102.5

194.8

Sk 64.9

89.3 S I, = ,. .

* 62.1

Estimated normalization 

factor ' , /}

87.2

62.5

7 _ 82.9 

76.4 

84.3

89.5

87.1

62.5

r 82.7 

' if. i

7o.J 

84.2

89.4

87.0

62.1

? _ 82.6

' 7A 1

/O.I

83.9

89.3

Reproduced data matrix (product of matrices of a£ and/jy)

0.99

.76

.99

*ii ~ .92

.99

.98

0.14

.25

.15

.18

.15

.14

0.02

.60

.07

.21

.05

.00

2The communalities (A,2 = 2 a^2) are 1.00, 0.99, 1.00, 0.93, 1.00, and 0.99. 

 'From reproduced data matrix.

initial factor loading matrix for the oblique model. If any 

of these oblique loadings exceeds unity, an iterative 

procedure is begun to select different end-member samples 

for which all initial oblique loadings are unity or less.

A different iterative procedure has been used to derive 

the initial oblique loadings in tables 10, 11, and 12. The 

criterion in the iteration is that none of the oblique 

composition loadings, rather than the initial loadings, 

exceeds unity. The unsealed scores for the oblique models 

are taken directly from the reproduced data matrices 

(tables 10,11,12). The alternative procedure is discussed in 

a subsequent section dealing with the identification of 

samples of extreme composition.

INTERPRETATION OF 

EIGENVALUES AND COMMUNALITIES

The first three eigenvalues, A , of the 6 by 6 matrices of 

cosine theta derived from the original and transformed 

data matrices (tables L4, IB, and 1C) and their cumulative 

proportions of TV, here after abbreviated as CPN, are given

in table 13. The values of CPN indicate the proportions of 

the total sums of squares in the normalized data matrices 

that will be explained by principal components, varimax, 

or oblique models having the corresponding numbers of 

factors. It was shown in part I that three-factor models 

account for all the variability in the data matrices. That is, 

the final factor models could be used to reproduce the 

original, transformed, and normalized data matrices 

exactly, except for differences due to round off. The values 

of CPN in table 13 show that two-factor models will 

explain 98 percent, 94 percent, or 91 percent of the 

variability (as a sum of squares) in the normalized 

data matrices, depending on whether or not a data 

transformation was used and on which type of 

transformation it was. For example, when the initial factor 

loadings for any of the three models represented in table 10 

are multiplied by the corresponding matrices of unsealed 

scores, the product is an estimate of the normalized 

original data in table IAA, referred to as the reproduced 

data matrix. The total sum of squares in the reproduced 

data matrix of table 10 is 98 percent of N = 6, the total
G16 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

TABLE II. Initial loadings, a ik', unsealed scores, fk'j', scale factors, sk, estimated normalization factors, Tj, and reproduced data matrix, "x;j', for 

the two-factor principal components, varimax, and oblique factor models derived after transforming each variable to proportions of the 

maximum value

[The products of the matrices of a/^ and/^! are approximations of the normalized transformed data in table \BB]

Principal 

components 2

Varimax 2 .....

Oblique ......

Initial load in

0.969

.706

.990

°ik ~ .851 

.991

.927

0.936

.241

.845 a'.' = CA t

ik .546 

.898

.951

0.954

0

.799

a* = .403 

.882

1.000

5 s- alk

-0.234

.627

-.037 

.294 

-.131

-.324

0.341

.913

.517

.716 

.440

.244

0.119

1.000

.353

.677 

.246

0

Unsealed scores, /£. 

f" _ 0.880 0.416 0.231 

kj -.407 .408 .817

f = 0.958 0.120 -0.260

kj .148 .570 .808

3y» = 0.95 0.25 -0.05

^ .37 .55 .68

Scale factors, s^

0.950 

* ~ 16.001

1.188

Sk ~ 1.575

1.044

Sk ~ 1.268

Estimated normalization 

factors ',7J

0.99

1.28

.96 

'' ~ 1.09 

.97

1.05

1.00

1.28

.96

*i ~ 1.09 

.97

1.05

0.99

1.27

.96

'' ~ 1.09 

.96

1.04

Reproduced data matrix (products of matrices of a'.'k aud f^'.)

0.95

.37

7" -89

IJ .63

.93

.95

0.31

.55

.40

.47

.36

.25

0.03

.68

.20

.44

.12

-.05

*> "y

2The communalities (A,.=2fl^2) are 0.99, 0.89, 0.98, 0.81, 1.00, and 0.96.

 'From reproduced data matrix.

sum of squares in the normalized data of table \AA. 

Similarly, the total sum of squares in the reproduced data 

matrix of table 11 is 94 percent of that in the normalized 

data of table \BB, and the total sum of squares in the 

reproduced data matrix of table 12 is 91 percent of that in 

the normalized data of table ICC. For those who check 

these assertions by computation, the minor differences that 

will be observed are due solely to round off. Thus, the 

eigenvalues of the cosine theta matrix indicate the overall 

degree to which principal components, varimax, and 

oblique Q-mode models will sueceed in reproducing the 

normalized forms of the original data.

The communalities of the samples, as previously 

mentioned, are computed as the sums of squares of the 

initial loadings on the m principal components or varimax 

axes (tables 10, 11,12). The sample communalities are also 

equal to the sums of squares of the corresponding rows in 

the data matrices reproduced from the initial loadings and 

unsealed scores (tables 10, 11, 12). Thus, the communalities, like the eigenvalues of the cosine theta

matrices, indicate the degree to which the models will 

succeed in reproducing the original normalized data. In 

fact, the average sample communality is equal to the 

corresponding value of CPN.

It will be shown in a following section on goodness-of-fit 

measures that the eigenvalues and sample communalities 

can be misleading where they are used as indicators of the 

degree to which the factor model can be used to reproduce 

the original data.

SCALING THE FACTOR 

LOADINGS AND SCORES

The appropriate factors, sk, for scaling the initial factor 

loadings and scores to conform with the original data were 

computed from equation 9 of part I and are given in tables 

10, 11, and 12.The composition loadings, a/£, derived by 

dividing the initial loadings, a/£, by the corresponding 

scale factors to give a'^ and adjusting the values of a ik ' so 

that they sum to unity across k, are given in tables 14, 15,
FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G17

TABLE 12. Initial loadings, a^, unsealed scores, f^'j, scale factors, s k, estimated normalization factors, ?j, and reproduced data matrix, x fj, for 

the two-factor principal components, varimax, and oblique factor models derived after transforming each variable to proportions of the range

[The products of the matrices of a^ and f/^ are approximations of the normalized transformed data in table ICC]

Model

Principal

components 2

Varimax .....

Oblique ......

Initial loadings,

0.963

.234

.985

a 'i'k = .603

.995

.906

0.991

-.054

.898 

a'i'k = -357

.970

.979

1.021

0

.947

a ut ~ .418

1.010

1.000

"Ik

-0.249

.830

.094

.643

-.095

-.375

0.083

.861

.414

.8UO

.239

-.054

0.161

1.000

.541

.963

.341

0

Unsealed scores, /^.'

0.906 0.400 0.138

f*J ~-.374 .608 .700

0.979 0.177 -0.101 J kj -_.Q54 .706 .706

3 f » 0.96 0.13 -0.14

' -.10 .60 .61

Scale factors, sk

s = °- 792

k 1.720

5 _ i.ooo

* 1.035

1.087

Sk "1.290

Estimated normalization 

factors, 7/

0.93

1.28

.77

'' = .88

.83

1.08

0.93

1.29

.75

'/' = 87 .o/

.83

1.08

0.94

1.29

.77

'' = .88

.84

1.09

Reproduced data matrix (product of matrices of a£ and/")

Xjj =

0.97

-.10

.86

.31

.94

.96

0.23

.60

.45

.63

.34

.13

-0.04

.61

.20

.53

.07

-.14

2The communalities (h2. = 2tf^2) are 0.99, 0.74, 0.98, 0.78, 1.00, and 0.96. 

3From reproduced data matrix.

and 16. The scaled scores, f^., derived simply by

multiplying the unsealed scores, j^y, by the scale factors, 

are also given in tables 14, 15, and 16. The products of the 

matrices of a ik and/^ are approximations of the original 

data where no data transformation was made,1 or of the 

transformed data where transformations were used? Where 

a transformation was not used (table 14), the scaled scores 

are equal to the composition scores. Where the data were 

transformed to proportions of the maximum values for 

each variable (table 15), the composition scores are derived 

by multiplying the scaled scores by the corresponding 

values of xmaxj from table 1. Where the data were 

transformed to proportions of the ranges for each variable 

(table 16), the composition scores are derived by 

multiplying the scaled scores by the quantity xmaXj 

and adding xmin.

T ABLE 13. Eigen values of the matrices of coefficien ts of proportional 

similarity (cosine theta) derived from the data matrices in table 1, 

and CPN i

Data not transformed 

table \A

A

5.5986 

.2987 

.1027

CPN

0.93 

.98 

1.00

Data transformed 

to proportion of 

maximum values, 

table IB

A

4.9840

.6577 

.3583

CPN

0.83 

.94 

1.00

Data transformed 

to proportion of 

ranges, table 1C

A

4.1256 

1.3225 

.5518

CPN

0.69 

.91 

1.00

1 Compare tables \A and 14.

2 Compare tables IB and 15, and 1C and 16.

Cumulative proportions of N.

GOODNESS-OF-FIT MEASURES 

FOR THE FACTOR MODEL

Because the number of factors, m, in the models 

represented in tables 10 12 and 14 16 is less than the 

number of variables, M, less then 100 percent of the 

variability in the data will be explained by the models, as
G18 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

TABLE 14. Composition loadings^ a ik , composition scores, fkj = f'^, 

and reproduced data matrix, x^ , for the two-factor principal 

components, varimax, and oblique factor models derived without 

transforming the data

[The products of the matrices of a^ and/^y are approximations of the original data in table I A]

Model

Principal 

components ....

Varimax . .

Oblique . . .

Composition loadings, a^

1.117 

.707 

1.070 

aik = .951 

1.086

1.139

0.399

.142

.369 

* ~ .295 

.379

.412

0.950

0

.840

aik = -565

.877

1.000

-0.117 

.293 

-.070 

.049 

-.086

-.139

0.601

.858

.631 

.705 

.621

.588

0.050

1.000

.160

.435

.123

0

Composition scores, /£ = /£. 

fkj =f'kj = 75 -° 13.3 H.7 

-19.1 20.4 98.6

fkj =A'= 176.5 5.7-82.2 

26.0 17.0 57.0

fkj =fkj = 88 12 0

47 15 37

Reproduced data matrix (product of matrices of a^ and /jy)

86

47

x..= 82

70

83

88

12

15

13

14

13

12

2

37

6

16

4

0

Row-sum

100

99

101

100

100

100

the eigenvalues in table 13 indicate. Consequently, the data 

matrices reproduced from the composition loadings and 

composition scores (tables 14, 15, 16) are not equal to the 

original data in table \A. Appropriate comparisons are 

made as shown in table 17.

The values of d* and dj* in table 17 are, respectively, the 

mean and standard deviation of the factor-model 

residuals, and the values of r and rj are, respectively, the 

correlations and coefficients of determination for the 

reproduced and original data. If the model explains all the 

data perfectly, all values of d * and d*j* are zero, and all 

values of ry and rj are unity. In real applications of the 

method, the values of dj are found to be small in 

comparison to the original data values, indicating that the 

models are unbiased. The coefficients of determination 

have been found to be particularly informative in that they 

show the proportions of the total variance in each variable 

explained by the model. That is

(32)

2. where s(x)j is the variance in theyth column of the original 

data matrix.

_ ^

The values of ry- in table 17 show that none of the factor 

models explains much of the variance in variable 2. This 

shortcoming of the two-factor models is not evident from 

the eigenvalues given in table 13, and the indication is clear 

that the eigenvalues alone can be misleading measures of 

the mathematical adequacy of a factor model. It will be 

shown in the following section and in part III that 

diagrams factor-variance diagrams showing the proportion of the total variance in each variable, r2, 

accounted for by factor models with up to M reference 

axes, are useful devices for selecting the number of 

end-members that the model should include.

The reason that the original data cannot be reproduced 

exactly by factor models that have m < M is related to the 

equality f, = 1/2afk ' , where t/ is the quantity that was 

used to convert the fth row of the original or transformed 

data to the normalized form. Where less than the total 

variability in the data is explained by the model (that is, 

where the number of factors, m, is less than the rank of the 

matrix of cos 0), this equality does not hold true. However, estimates of the normalization factors will approach 

the correct values as the proportion of the variability 

explained by the model approaches one, just as the 

reproduced data matrices approach the matrix of original 

data. The estimates of the normalization factors, f}, are 

given in tables 10, 11, and 12 for comparison with the 

corresponding values of tj in table 1.

AN EXAMPLE OF THE EFFECT

OF THE CONSTANT 

ROW-SUM CONSTRAINT

Inspection of factor-variance diagrams in real 

applications shows that as the number of factors, m, in the 

model is increased, the model explains somewhat greater 

portions of the variance in all, or nearly all, of the 

compositional variables. Generally, however, as will be 

demonstrated in part III, the addition of a factor to the 

model will increase the variance accounted for in a few of 

the variables to a greater degree than it will increase the 

accountable variance in others. The particular variables 

that respond to the same factor depend, at least in part, on 

the correlations among the variables, and this raises the 

question of the effect of the constant row-sum constraint 

that was described by Chayes (1960). Chayes showed that 

if the row-sums of the data matrix are constant across all 

rows, each variable must show some negative correlation 

with at least one of the others, even if the variables have no 

genetic relationships whatsoever with each other. The 

question here is the degree to which the constant row-sum 

constraint affects the configuration of the factor-variance 

diagrams. The question cannot be fully answered at this 

time, but it is possible to show the nature of a diagram, 

derived from simulated data, wherein nearly all the 

intercolumn correlation is a result of the constraint.
Q-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G19

TABLE 15. Composition loadings, aik , scaled scores, fkj , composition scores, fkj, and reproduced data matrices for the two-factor principal 

components, varimax, and oblique factor models derived after transforming each variable to proportions of the maximum value

[The products of the matrices of a,^. and fA are approximations of the transformed data in table IB. The products of the matrices of a/£ and fy are approximations of the

original data in table \A\

Model

Principal 

components .

Varimax ......

Oblique ......

Composition loadin

1.015

.950

1.002 

a* ~ .980

1.008

1.021

0.784

.259

.684

Q it ~ cni '* .503

.730

.838

0.907

0

.733

Oft A"\t\ K .420

.813

1.000

gs, aik

-0.015

.050

-.002 

.020

-.008

-.021

0.216

.741

.316

.497

.270

.162

0.093

1.000

.267

.580

.187

0

Scaled scores f^;

f > 0.835 0.395 0.219 

kj -6.51 6.52 13.06

f _ 1.138 0.143 -0.309

kj .233 .898 1.273

f 0.992 0.261 -0.052

kj .469 .697 .862

Composition scores, /^y

79.3 11.9 8.8 

'V ~ -618 196 522

108.1 4.3 -12.4

'*>~ 22.1 26.9 50.9

94 8-2

JkJ 45 21 34

Reproduced data matrices

(Product of matrices of a,-£ and f^j )

X. 

.94

.47

.85

.69

.89

.99

0.30

.70

.38

.52

.35

.26

0.03

.86

.19

.48

.12

-.05

(Product of matrices of a,-£ and/^-)

90

44

65

85

94

9 1

21 34

12 8

16 19

10 5

8 -2

Row-sum

100

99

101

100

100

100

The exercise was carried out using Chayes and Kruskal's 

(1966) equation 7 to search for the column variances in an 

open matrix with no intercolumn correlation that would 

yield, on adjustment of all rows to sum to 100, a closed 

matrix with column means and variances equal to those 

computed from data given by Lee and Van Loenen (1971, 

table 5) for 81 samples of a granitoid intrusive in eastern 

Nevada. (The computed means and standard deviations 

are given in table 18 of part III.) Chayes and Kruskal's 

equation 7 yielded the following variances for the open 

matrix:

Column

1

2

3

4

5

6

7

8

9

Representing

SiO2

A1203

Fe2 03

FeO

MgO

CaO

Na2O

K2O

H2° +

122.417

-2.2181 

.2394 

.4104 

.3104 

.9272

-.0658 

.4306 

.0950. (33)

A 500 by 9 matrix was constructed from a population 

of normally distributed psuedorandom numbers with

column means equal to the computed means for the 

granitoid intrusive and column variances equal to these 

open variances where they are nonnegative or equal to zero 

where they are negative. After setting all negative values 

that appeared in the 500 by 9 matrix to zero, each matrix 

row was adjusted to sum to 100 by dividing through by the 

row-sum and multiplying by 100. The matrix then became 

closed in the terminology of Chayes (1960).

The column means and standard deviations of the closed 

matrix were

Column Mean Standard deviation

71.29 3.59

15.20 1.78

1.15 .48

1.19 .68

.94 .53

2.37 1.00

3.75 .44

8 3.43 .78

9 .69 .31. (34)

The matrix of closed data was then used to derive 

varimax factor models containing two to eight factors 

(only the first eight eigenvalues of the matrix of cosine
G20 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

TABLE 16. Composition loadings, aik , scaled scores, fk - , composition scores, f kj, and reproduced data matrices for the two-factor principal 

components, varimax, and oblique factor models derived after transforming each variable to proportions of the range

[The products of the matrices of a^ and fy are approximations of the transformed data in table 1 C. The products of the matrices of a^ and /^y are approximations of the

original data in table \A\

Model

Principal 

components .

Varimax ......

Oblique ......

Composition loadings, a^

1.135

.379 

.958 

fl * ~ .670

1.046

1.236

0.925

-.069

.692 

"*- .314

.808

1.056

0.883

0

.675 

a'k ~ .340

.779

1.000

-0.135

.621 

.042 

.330

-.046

-.236

0.075

1.069

.308 

.686

.192

-.056

0.117

1.000

.325 

.660

.221

0

Scaled scores, /^y

,., 0.718 0.317 0.109 

/w ~ -.643 1.045 1.203

f , 0.979 0.177 -0.101 

/w ~-.056 .731 .731

,, 1.043 0.141 -0.152 

/^~-.129 .774 .787

Composition scores,/^:

82.3 11.6 6.1 

J*i~ 21.1 31.2 47.7

, 94.1 7.8 -1.8 

Jkj 47.5 22.7 29.8

.97 7 -4 

kj 44 24 32

Reproduced duta matrices

(Product of matrices of a^ and /^y)

0.90

-.13

Xf'j = .66

.27

.78

1.04

0.22 -d.04

.77

.35

.56

.28

.14

.79

.15

.47

.06

-.15

(Product of matrices of a^ and/^y)

91

44

80

*V ~ 62

85

97

9

24

12

18

11

7

0

32

8

20

4

-4

Row-sum

100

100

100

100

100

100

TABLE 17 . Measurements of correspondence between the observed data x^ , and the data reproduced from the two-factor models, x fj

[dj is the mean of the six values of rf,y = (*,y - Xjj); dj* is the standard deviation of the six values of rf,y. ry is the correlation coefficient between ?fy and

Xjj for the six pairs; rj is the coefficient of determination]

Transformation

None (table 14) ......

To proportions of the

maximum for each

variable (table 15) ..

To proportions of the

range for each

variable (table 16) ..

di*

Variable (/)

1 2 3

-0.6 0.7 -0.1

-.1 .1 .0

-.2 .9 -.7

dr

Variable (/)

1 2

5.5 8.4

3.2 7.9

2.6 7.7

3

2.9

4.8

5.9

rj

Variable (/)

1 2 3

0.94 0.28 0.98

.98 .41 .93

.99 .49 .90

i

0.88

.96

.98

O 2

Variable (/)

2

0.08

.17

.24

3

0.96

.86

.81

theta were greater than zero), and each of the models was 

used to reproduce an estimate of the original matrix. 

Although the models employed the varimax axes as 

reference vectors, use of the principal components axes or 

any oblique axes would have led to the same reproduced 

data and the same goodness-of-fit measures. The 

goodness-of-fit measures included the coefficients of 

determination, eight sets of rj 's, which were used to

construct the factor-variance diagram of figure 2. In 

constructing figure 2, all values of rj- for the single factor 

model were set equal to zero, because a model with only 

one factor would lead to a reproduced data matrix with no 

intracolumn variance.

The factor-variance diagram shows that, for the most 

part, the addition of a factor to the model results in an 

important increase in the accountable variance for only
g-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G21

1 2345678 

NUMBER OF FACTORS

FIGURE 2. Factor-variance diagram derived from a closed matrix of 

random numbers having column means and variances similar to those 

of the data matrix for the granitoid intrusive in the southern part of 

the Snake Range, Nev.

one variable. As will become apparent in part III, this is far 

different from most of the situations encountered in 

treating real geochemical and petrographic data and 

suggests that the effect of the constant row-sum on the 

nature of the factor-variance diagrams is exceedingly 

minor. However, even if this were not the case, the factorvariance diagrams would still provide a correct appraisal of 

the degree to which each factor model will account for the 

total variance in the original data.

COMMON AND UNIQUE FACTORS

In general, the coefficients of determination, r?, increase 

for each variable with each increase in m, the number of 

factors in the model. This is in accord with the 

conventional concept of common factors. (See, for 

example, Harman, 1961, p. 15.) Common factors, 

representing processes or effects, are those that cause 

variability in more than one variable. However, the 

situation is not uncommon where the addition of one 

factor to the model will cause only small increases in rj for 

all but one variable but will cause a large increase in r,2 for 

that one variable. In conventional factor analysis 

terminology, this factor is said to be unique (Harman,

1967, p. 15). In some situations, more than one unique 

factor may appear.

It has been conventional in factor analysis applications 

to use common factors only, excluding those that are 

unique. If a major portion of the variation in a 

compositional variable (/' = u) can be explained only with a 

unique factor, and the unique factor is not included in the 

model, the model will, of course, fail to account for some 

increment of this variable. The increment, which can be 

ascribed to a unique process or condition, is that amount, 

// , which must be added to the reproduced value for the 

unique variable, xiu , to give the observed value for the 

sample, xiu , as in

X. =

iu

K + I. (35)

IU

The increment that must have been added to or subtracted 

from the fth sample by some process attributed to the 

unique factor can be estimated by rearrangement of 

equation 35 and solving for Iiu . Similarly, if two of the 

compositional variables (/' = u and j = v) respond to 

unique factors not included in the model, the increments, 

Iju and 7/v , can be estimated by rearrangement and 

simultaneous solution of equations 36a and 36b:

 "* in 

K+ /,. +/....

K + /,.. + /,

(36a)

(36b)

IU

COMPOSITIONS VERSUS COMPOSITION 

SCORES

Generally in this report, the composition scores have 

been referred to as properties of the reference axes and 

denoted by fy . However, the composition scores are 

essentially analogous to reproduced data designated by Xy 

in equation 17 of part I. That is, each sample vector 

representing a composition, x-- (l<y <M), has a set of 

composition scores, Xy , which represents the original 

composition modified to correspond with the factor 

solution. Thus, the only difference is in notation. The 

symbol, fa .- , is used for the scores of the k\h reference 

vector, and x^ is used for the scores (reproduced data) of 

the Ah sample vector.

When the number of factors in the model is less than the 

number of variables in the data matrix, most of the sample 

vectors are in different positions and are less than unit 

length after projection into the m-dimensional factor 

space. Consequently, the vectors no longer represent the 

exact compositions of the corresponding samples. That is, 

the composition scores, xfj , for the vectors are not exactly
G22 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

equal to the sample compositions, x^ . The same is true for 

the unsealed scores, x-j , and the normalized data, #;", 

and for the scaled scores, #J , and the transformed data, 

x-j . The degree of departure between the unsealed scores 

and the normalized data is inversely related to the sample 

communality.

The composition scores for the sample vectors or for 

reference axes that represent selected end-member 

compositions may be viewed as compositions that have 

been modified to conform with a compositional system. 

The degree of modification required will tend to be large 

when the vector communality is small and may be cause for 

declaring that a particular composition is inconsistent with 

the compositional series.

IDENTIFICATION OF THE SAMPLES 

OF EXTREME COMPOSITION

After it has been determined that the sample vectors 

occur mostly within a factor space of m dimensions and 

that the factor model, therefore, should contain m endmembers, it may be helpful to identify the m samples of 

extreme composition. By samples of extreme composition, 

we mean the m samples having compositions that can be 

combined in various proportions, none exceeding unity, to 

give approximations of the compositions of all the other 

samples.

In Q-mode analyses of data matrices with variable 

row-sums, the investigator has no means of scaling the 

initial loadings and scores and computing the composition 

values; in a search for extreme samples, his only recourse is 

to find the sample vectors on which all the initial loadings 

for the other sample vectors are less than one that is, the 

procedure suggested by Imbrie (1963). This will frequently 

serve to identify the same samples as being of extreme 

composition as would have been identified using the 

composition loadings. Examination of the initial varimax 

loadings in table 12, however, shows that this will not 

always be true. According to Imbrie's procedure, these 

loadings indicate that samples 1 and 2 are of extreme 

composition; the initial loadings of all samples with respect 

to the oblique vectors represented by them are all one or 

less:

1.000

0

.928

.409

.989

.980

0

1.000 

.391 

.896 

.182

-.157. (37)

Although samples 1 and 2 are, indeed, of extreme 

normalized composition, the fact that they are not of

extreme actual composition becomes evident when the 

composition loadings are derived. The scale factors to be 

used for deriving the composition loadings are computed 

from the unsealed scores for samples 1 and 2, taken from 

the reproduced matrix of normalized data at the bottom of 

table 12:

/,", = 0-97 = °-23 =-o.o4

=-o.io /" =o.60 /:; = o.6i. (38) 22

The scale factors, from equation 9, are 0.931 and 1.290. 

Division of the respective columns of matrix 37 through by 

these values and adjustment of each row to sum to unity 

give the composition loadings as follows:

aik =

1.000

0

.767

.387

.883

1.131

0

1.000

.233

.613

.117 

-.131, (39)

These values indicate that sample 6 is beyond the range of 

compositions represented by samples 1 and 2, a fact that 

also seems apparent from the matrix of original data in 

table \A. Alternatively, returning to the initial varimax 

loading matrix in table 12 and selecting samples 6 and 2, 

rather than 1 and 2, as the end-members for the oblique 

model, the initial oblique loadings are as given in table 12, 

and the oblique composition loadings, which are all in the 

range from zero to plus one, are as given in table 16. The 

revised scale factor, s^ for the first end-member (sample 6) 

is found, using equation 9, to be 1.087 (table 12).

The products of the matrices of the composition 

loadings and composition scores for the oblique model are 

equal to the corresponding products for the principal 

components and varimax models to two significant figures 

(table 16). Thus, the oblique model is consistent with the 

other types.

The fundamental differences between the procedure 

suggested here and the method for deriving the oblique 

model given by Imbrie (1963) are that here (1) the reference 

vectors are taken as the vectors representing the samples of 

extreme composition rather than those of extreme 

normalized composition, and (2) the scale factors for 

adjusting the factor loadings are computed from the 

unsealed composition scores for the extreme samples 

rather than from their original or transformed 

compositions. The advantage of this second modification 

is that the reproduced data from the oblique model will be 

exactly the same as the data reproduced from the principal 

components or varimax models. In fact, derivation of the 

scale factors from the unsealed composition scores rather
Q-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G23

than from the original or transformed compositions will 

lead to the same reproduced data regardless of any choice 

of reference axes within the /w-dimensional factor space.

MODIFICATION OF A VARIMAX 

MODEL TO AN OBLIQUE MODEL

Commonly, the computed set of composition scores for 

one or more of the varimax factor axes contains one or 

more negative values. If this is viewed as undesirable for 

the purpose of the investigation and if the composition 

scores computed for the sample vectors tend to be all 

positive, it may be helpful to move the reference axes, one 

at a time, from the varimax positions by small increments 

toward the sample vectors. Computer programs can be 

easily prepared to do this, wherein, after each increment, 

the composition scores can be determined by the procedure 

give in part I, and the process terminated when all scores 

are nonnegative. The procedure leads to a model that 

approximates the conventional varimax solution, but the 

reference axes are oblique. The advantage of the oblique 

model is that the end-members for the model are more 

likely to be of realistic composition.

If the /w-dimensional model does not account for large 

proportions of the variances for each variable, however, 

the sample vectors may have composition scores that are 

not all positive, and movement of the reference axes 

toward them will not have the desired effect. The only 

recourse then is to increase m.

In general, as the reference axes are moved from the 

varimax positions towards the sample vectors, the 

composition scores of the axes move towards the range of 

zero to K, and the loadings of the sample vectors on each 

of them become increasingly variable. As the axes are 

moved past sample vectors, the composition loadings move 

outside the range of zero to one.

TESTING THE PLAUSIBILITY 

OF ALTERNATIVE END-MEMBERS

As pointed out in part I, the unsealed varimax scores 

may be used to find the vector representation of any given 

composition in the varimax space and to determine 

whether or not the composition could be that of an 

end-member in the compositional system. For purposes of 

illustration, we may suppose that the varimax model 

derived after transforming the variables to proportions of 

their maximum values (table 15) is unsatisfactory because 

one of the composition scores, f{3 , is negative and because 

the composition scores for the varimax axes, in general, are 

not sufficiently close to those believed, on the basis of 

other evidence, to have been involved in the genesis of the 

petrologic system under study. Then it may be desirable to

test the following five observed, theoretical, or 

hypothetical phases as possible end-member compositions:

60

10

5

90

25

30

80

90

10

25

10

10

5

0

50. (40)

The transformed values in this case are obtained by 

dividing the value of zfj by the corresponding values of 

from table \A, giving

Phase

1

2

3

4

5

Z U = 0.632

.105

.053

.947

.263

1.000

2.667

3.000

.333

.833

0.250

.250

.125

0

1.250. (41)

The transformed compositions are normalized to give

Phase 

I

2

3

4

5

Z 2j =

0.523

.039

.018

.943

.172

0.827

.995

.999

.332

.546

0.207

.093

.042

0

.820. (42)

Applying equation 25 of part I, the initial loadings on the 

varimax axes are computed for each phase. These are 

found to be

Phase 

1

2

3

4

5

0.546

.133

.126

.944

.018

0.716

.648

.605

.329

.999

h2 =

0.811

.438

.382

.998

.999 (43)

The values of a'£ may be converted to composition loadings 

by dividing by the corresponding scale factors for the 

varimax model (table 11) and adjusting to sum to unity. 

The composition loadings can then be postmultiplied by 

the varimax composition scores,/^ , from table 15 to give 

composition scores of the 5 phases. These are

Phase

1

2

3

4

5

x 2j = 

x 3j = 

x4j = 

x5j =

65.37

40.50

40.71

90.23

24.13

15.54

22.10

22.04

8.99

26.41

19.10 

37.40 

37.24 

0.81 

49.45. (44)

Any of the phase composition scores could be the 

composition of an end-member in the petrologic system if
G24 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

the composition loadings in the final model are subject to 

plausible geologic interpretation. However, the composition scores differ from the phase compositions tested to a 

degree indicated by h2 in 43. Phases 1, 2, and 3 might be 

rejected because the revised compositions are appreciably 

different from the original compositions tested. (Compare 

40 and 44.) Phases 4 and 5 may be accepted as possible 

end-members if the composition loadings of the other 

samples on the new reference axes are geologically 

plausible.

If the decision were made to describe the samples in 

terms of phases 4 and 5, one would form the following 

partition of the matrix in 43:

 _ 0.944 0.329 

°ik ~ 0.018 0.999 '

(45)

and then proceed as described in equations 28-31 of part I.

SUMMARY

The recommended steps to be followed in the 

development of a geochemical or petrologic model using 

the extended form of Q-mode factor analysis can be 

summarized as follows. It is assumed that the rows of the 

data matrix have a constant sum.

1. Use the Q-mode factor analysis computer program 

(CABFAC) of Klovan and Imbrie (1971), rotating 

2 to M (or 10 if M> 10) varimax axes. Using the 

unsealed varimax scores, output from CABFAC, 

compute all the varimax scale factors and varimax 

composition scores for each rotation.

2. Using the appropriate scale factors, adjust each matrix 

of initial varimax loadings to a matrix of varimax 

composition loadings.

3. Multiply each matrix of varimax composition loadings 

by the corresponding matrix of varimax composition 

scores to produce M-\ approximations of the 

original data, generally in units of percent or parts 

per million.

4. For each of the M 1 approximations, generate the 

measures of correspondence to the original data that 

are given in table 17.

5. Use the measures of rj to construct a factor-variance 

diagram and select the number of end-members, m, 

that the model should contain in order to account for 

as much of the variance in as many of the variables 

as possible. The principle of parsimony (Imbrie, 

1963) should be applied.

6. The next step is to select the reference axes or endmember compositions, and here geologic criteria

play an important part. For an /n-dimensional 

model, it is necessary to form an m by m matrix of 

the initial loadings of the reference vectors on the 

varimax axes. One or more of the reference vectors 

may be those representing samples in the original 

data matrix if those samples are thought to represent 

end-members or extremes in the petrologic system. 

They also may be the varimax axes or any other 

vectors within the varimax space as long as the composition scores for the axes or vectors approach the 

compositions of materials thought possibly to have 

been end-members in the petrologic system. The reference vectors may also be vectors with composition 

scores approaching the compositions of materials 

that had been tested for fit in the factor space. In 

brief, if the reference vectors are not the varimax 

axes or vectors representing the original samples, 

they are selected either by examining the composition scores of selected vectors or by finding vectors 

with composition scores approaching the compositions of materials tested and found to be representable in the factor space. The techniques are given in 

parts I and II.

7. When an m by m matrix of initial loadings of the reference vectors on the varimax axes has been selected, 

the initial loadings on the reference axes (generally 

oblique axes) are derived using the method given 

here, which is originally from Imbrie (1963).

8. The unsealed scores of the new reference axes are computed by multiplying the matrix of initial loadings on 

the varimax axes by the matrix of unsealed varimax 

scores. The scale factors are then computed for the 

reference axes and their composition scores are 

derived.

9. The scale factors and the initial loadings on the new 

reference axes (derived in step 7) are then used to 

compute the composition loadings on the new 

reference axes. The composition loadings, especially 

their signs, are examined for geologic plausibility. If 

they are not acceptable, return to step 6 for the selection of alternative reference axes.

10. As a partial check for computational errors, postmultiply the matrix of composition loadings on the 

selected reference axes by the matrix of their composition scores. Errors are present if the departure 

of the product matrix from the corresponding 

matrix derived in step 3 is beyond that which can be 

attributed to round off.

The steps outlined here were followed in application of the 

Q-mode technique to some petrologic mixing problems, as 

described in part III.
Q-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS

Part III. Applications to some Petrologic-Mixing Problems

G25

INTRODUCTION

The extended form of Q-mode factor analysis has been 

applied to four problems in igneous petrology to illustrate 

how it may be used and the kinds of information it may 

yield. Problems in igneous petrology were chosen because 

igneous systems generally exhibit compositional variations 

that are simpler than those in sedimentary and 

metamorphic systems and because relatively precise 

analyses of igneous suites are available in abundance. The 

simpler compositional variations in igneous rocks are 

attributed to the relative simplicity of the systems in which 

they formed, increasing the likelihood that a small number 

or realistic end-members can be found that will account for 

a large part of the compositional variation.

The four sample suites chosen for examination are from

(1) a Pleistocene rhyolite-basalt complex on the Gardiner 

River, Yellowstone National Park, Wyo. (Fenner, 1938),

(2) a Jurassic granitoid intrusive in the southern part of the 

Snake Range, Nev. (Lee and Van Loenen, 1971), (3) lavas 

and pumices of the 1959 summit eruption at Kilauea, 

Hawaii (Murata and Richter, 1966), and (4) the layered 

series of the Skaergaard intrusion, Greenland (Wager and 

Brown, 1968). Analyses of the samples are from the 

references cited, but only the determinations on 

constituents that are essential to the principal rock-forming 

minerals were used in this exercise; these constituents are 

SiO2 , A1 2O3 , Fe 2O 3 , FeO, MgO, CaO, Na2O, K2O and 

H2 O+. The same nine constituents were used in all four 

studies so that the results would be more easily comparable 

from one to another. Subsequent studies should, perhaps, 

include such constituents as MnO, TiO2 , P2O5 , and 

Cr 2O 3 , especially if they might be diagnostic of processes 

that may have led to important amounts of variation in the 

rocks under study.

Means and standard deviations for the compositional 

data of all four sample suites are given in table 18.

Prior to the factor analyses, each data matrix was 

transformed so that each column (chemical variable) 

ranged from zero to one (that is, to proportions of the 

range).

The first step in the examination of each sample suite 

was the construction of a factor-variance diagram. The 

diagram, as shown in part II, is merely a plot of the values 

of rf, as given in table 17, versus the number of factors in

7

the model. The values of ry were derived from varimax 

models but would be identical if a principal components 

model or some oblique model had been used. Factorvariance diagrams show immediately and concisely the

complexity of the petrologic system and indicate which 

compositional variables can be explained, and to what 

degree, by models with any given number of end-members.

TABLE 18. Means and standard deviations (in weight percent) for four 

sets of data examined by Q-mode factor analysis

[All analyses were adjusted before computation of means and standard deviations so that the nine 

constituents summed to 100]

Data sets

Constituents Rhyolite-basalt Granitoid Lavas and pumice Layered series of

complex, Yellow- intrusive in the from the 1959 Skaergaard

stone National Park, southern Snake summit eruption intrusion,

Wyo. Range, Nev. at Kilauea, Hawaii Greenland

(15 samples) 1 (81 samples) 2 (22samples) 3 (19samples) 4

Means

SiO2 ......

A1 2 03 ....

Fe 2O3 ....

FeO ......

MgO ......

CaO ......

Na 2O . . . . .

K 2O ......

H2 O +

... 62.62

... 14.60

. . . 2.07

. . . 4.65

. . . 4.07

. . . 5.98

... 3.11

. . . 2.40

0.50

71.72

14.97

1.18

1.14

.94

2.32

3.69

3.38

.67

49.89

12.12

1.88

10.14

13.09

10.34

2.00

.49

.03

45.34

12.52

4.29

19.70

5.65

9.00

2.69

.31

.50

Staadard deviations

SiO 2 . . 

A1 203 

Fe 20 3 

FeO .. 

MgO.. 

CaO .. 

Na 2O . 

K 2O .. 

H2 O +

7.69

1.15

.51

2.58

2.34

3.00

.31

1.34

.10

3.11

1.08

.50

.64

.56

.97

.32

.73

.31

1.26

1.39

.70

.70

4.12

1.23

.25

.06

.03

6.%

6.00

3.16

10.08

6.15

2.38

1.25

.15

.38

' Data from Fenner(1938, table 6), excluding samples Y.P. 904, Y.P. 1024, and Y.P. 1025. 

j? Data from Lee and Van Loenen (1971,table 5), excluding samples 71, 12, 77, 81, 85, and 87.

3 Data from Murata and Richter (1966, table 1) and Wright (1973, table 2). Same selection 

of samples as used by Wright.

4 Data from Wager and Brown (1968, table 5).

RHYOLITE-BASALT COMPLEX 

ON THE GARDINER RIVER, 

YELLOWSTONE NATIONAL PARK, WYO.

Exposures of the rhyolite-basalt complex occur in an 

area about 76 m wide and 274 m long within the canyon of 

the Gardiner River in the northern part of Yellowstone 

National Park (Fenner, 1938). The complex forms part of 

the Obsidian Creek Member of the Plateau Rhyolite of 

Pleistocene age in the third volcanic cycle of Christiansen 

and Blank (1972, p. B13). The thickness of the complex is 

difficult to determine, but it has been measured as less than
G26 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

30 m at one point. The complex consists of basalt 

impregnated with rhyolite and rhyolite containing 

abundant inclusions of basalt. The rhyolite penetrating the 

basalt forms veinlets and pseudodikes, and an 

intermingling of rhyolitic and basaltic constituents is 

evident from the groundmass and phenocrysts throughout 

most of the complex. Fenner showed that the compositions 

of most analyzed samples can be explained as mixtures of 

rhyolitic and basaltic end-members; the principal 

mathematical evidence is in the straight-line relationships 

as seen on silica-variation diagrams. According to the 

interpretation of Fenner, the complex formed by the 

extrusion of a rhyolitic lava that attacked an older basalt 

which constituted the floor and walls of the Gardiner River 

canyon. More than 68 percent of some basalt has been 

replaced by rhyolite, which entered the basalt as gaseous 

emanations from the rhyolite lava. In Fenner's 

interpretation, the gaseous emanations concept was used to 

overcome the difficulties of penetrating the basalt with 

rhyolitic liquid.

Fenner's interpretation of the origin of the complex was 

debated by Wilcox (1944), who explained the complex as 

the product of the mixing of rhyolitic and basaltic lavas 

that erupted simultaneously. Wilcox's principal evidence, 

rejected by Fenner (1944), consisted of observations of a 

chill phase in the rhyolite in contact with country rock but 

not in contact with basalt and the presence of quartz and 

orthoclase xenocrysts from the rhyolite in the basalt.

The factor analysis to be described here does not offer 

any new evidence with regard to the problem of the 

complex's origin, but it does emphasize and show some 

exceptions to the linearity of the compositional gradation 

between rhyolite and basalt. It also provides estimates of 

the relative proportions of basalt and rhyolite in each of 

the analyzed samples. The principal conclusion reached 

through the analysis is that the rhyolitic end-member, 

whether liquid rhyolite or gaseous emanations, had to have 

a composition near that of the least contaminated rhyolites 

now present in the complex. This tends to support the 

interpretation of Wilcox (1944). The exceptions to the 

linear compositional gradient from rhyolite to basalt 

involve Na2 O, H2O +, and Fe 2O3 . The first two of these 

are late-stage volatile constituents that appear to have 

behaved independently of other constituents, and the 

Fe2O 3 probably formed with oxidation after the mixing 

had occurred.

Fifteen of the 18 analyses given by Fenner (1938, p. 

1466) were used to derive the factor-variance diagram. 

Analyses of samples 1025, 904, and 1024 were not used 

because the source of the samples has been a matter of 

debate (Wilcox, 1944, p. 1067) and Fenner referred to 

them as varied from the standard composition. The 

factor-variance diagram is given in figure 3 and shows the 

proportion of the total variance in each constituent that

1.0 SiO2 , A12 03 , FeO. MgO. CaO, K2O

1 2345678 

NUMBER OF FACTORS

FIGURE 3. Factor-variance diagram for the rhyolite-basalt complex on 

the Gardiner River, Yellowstone National Park, Wyo.

can be explained by models containing one to nine factors, 

or end-members. The curves for all constituents originate 

at zero variance for one factor, inasmuch as single-factor 

models would lead to reproduced data matrices with no 

intracolumn variation. All but two of the curves rise 

sharply, indicating that two-factor models can account for 

most of the chemical variation in the complex. That is, a 

linear combination of two end-members in varying 

proportions can account for nearly all of the variation in 

six of the major oxides and 85 percent of the variation in 

Na2O. The variations in H2 O+ and Fe2O3 , however, can 

be explained equally well only if the model contains at least 

four factors. These constituents appear to represent factors 

that are largely unique; that is, each of these constituents 

appears to have behaved independently of the others, and 

each may have been controlled by some process that 

affected it alone. The process that accounted for the 15 

percent unexplained variance in Na2O and most of the 

variance in H2 O+ is interpreted as migration of late-stage 

volatile constituents in the lavas. The process controlling 

Fe2 O3 was probably later oxidation.

Because of the apparent importance of postvolcanic 

oxidation, the Fe2O3 in Fenner's analyses was recomputed 

and combined with FeO to give total iron as FeO, usually 

designated by 'FeO.' The complete set of analyses,
g-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G27

TABLE 19. Analyses of a rhyolite-basalt complex on the Gardiner River, 

Yellowstone National Park, Wyo., adjusted so that eight major oxides 

sum to 100 percent, and corresponding data reproduced from the factor 

model

[Analyses from Fenner (1938)]

Fenner (1938, p. 1476), it does not vary much more than 

several of the samples accepted as taken from the rhyolitebasalt complex. It seems probable that some kind of 

genetic relation exists between the samples that vary from

Sample^No. Si02 M^ , FeO, MgO CaQ ^Q ^Q ^^

A. Adjusted analyses, xy (in percent)

969 51.64 16.25 10.41 7.44 10.53 2.77 0.52 0.44 

918 54.33 16.06 9.49 6.70 8.98 2.87 1.04 .53 

948 54.49 15.74 9.49 6.75 9.30 2.76 .98 .49 

985 55.07 15.72 9.40 6.27 9.25 2.77 1.13 .40 

512 55.33 15.74 9.40 6.34 8.94 2.61 1.13 .52 

925 58.66 15.31 7.96 5.35 7.28 3.13 1.58 .72 

514 59.81 14.97 7.76 5.09 7.02 2.94 1.97 .45 

909A 62.24 14.82 6.79 4.27 6.09 3.27 2.02 .51 

909B 64.94 14.11 5.78 3.45 5.15 3.36 2.66 .56 

543 65.92 14.00 5.38 3.19 4.78 3.13 2.98 .61 

467 67.30 13.94 4.99 2.55 4.22 3.22 3.26 .53 

980 68.06 14.20 4.30 1.95 4.16 3.58 3.22 .53 

505 72.23 13.13 3.26 1.02 2.22 3.37 4.16 .61 

554 75.48 12.71 1.85 .37 1.10 3.58 4.59 .31 

914 75.72 12.70 1.72 .40 .83 3.44 4.80 .37

B. Reproduced dnti (in percent)

969 51.49 16.30 10.66 7.51 10.39 2.77 0.43 0.44 

918 54.81 15.78 9.42 6.50 9.07 2.87 1.02 .53 

948 53.85 15.% 9.81 6.82 9.48 2.76 .84 .49 

985 53.94 15.% 9.80 6.81 9.47 2.77 .85 .40 

512 53.63 16.01 9.92 6.91 9.60 2.61 .79 .52 

925 59.65 14.99 7.56 4.99 7.08 3.13 1.89 .72 

514 59.30 15.15 7.82 5.20 7.36 2.94 1.79 .45 

909A 62.36 14.63 6.62 4.21 6.07 3.27 2.34 .51 

909B 65.24 14.19 5.55 3.34 4.92 3.36 2.85 .56 

543 65.51 14.19 5.50 3.30 4.87 3.13 2.88 .61 

467 67.06 13.97 4.94 2.85 4.27 3.22 3.15 .53 

980 68.15 13.73 4.45 2.45 3.75 3.58 3.37 .53 

505 71.66 13.26 3.24 1.46 2.45 3.37 3.% .61 

554 76.72 12.57 1.45 .00 .54 3.58 4.82 .31 

914 76.63 12.60 1.51 .05 .60 3.44 4.80 .37

adjusted so that the eight oxides sum to 100 for each 

sample, is given in table 19A. The recomputed 

factor-variance diagram is essentially the same as the one 

given in figure 3, except for the absence of the curve for 

Fe 2O3 . 

The eigenvalues of the cosine theta matrix for the data in 

table 19A are given in table 20 and, like the factor-variance 

diagram, show that somewhere between two and four 

factors are required in the factor model. The eigenvalues, 

however, do not indicate the constituents whose behaviors 

depart from that of the others. 

A plot of the sample vectors in the two-dimensional 

varimax space is given in figure 4 and shows that the 15 

samples used in the analysis as well as the 3 samples (1025, 

904, and 1024) that vary from the standard have 

communalities near 1. The largest departures from 1 

among the 15 samples used in the analysis are for samples 

925 and 554, both with communalities of 0.94. The largest 

departure among all 18 samples is for sample 1025 with a 

communality of 0.90. Although sample 1025 is one of the 

samples that varies from the standard, as pointed out by

they may differ in age and occurrence.

TABLE 20. First eight eigenvalues of the cosine 

theta matrix for the rhyolite-basalt complex on the 

Gardiner River and CPN l

No. Eigenvalue CPN

11.353 0.7569 

2 3.298 .9767 

3 .258 .9939 

4 .080 .9993 

5 .007 .9997 

6 .003 .9999 

7 .001 1.0000 

8 .001 1.0000

Total......... 15.001

'Cumulative proportions of N.

Unmetamorphosed 

basalt 

Basalts metamorphosed at 

a distance from rhyolite

/ ^

o 8 *&<» ^\. 

"~ "TOOT \v Basalts modified by impregnation

INITIAL VARIMAX LOADING, .j' a(

o 

in o<

l ....| 1025

r r loT^^-^^ "^ \< with riivoiitic materia|

IP/TV \

\l // // //yf / .^ >^ 

^

Basalt greatly modified 

II // / ^ /'X \ <& bv C0ntact with rhvolite

11 // / ^/^.^^^ \ Rhyolite containing 

// // V\/SxX'^ «\ ^ disintegrated remnants 

If // .X^'^^ ^-~~-~~~~~~~~~~~'^\ of basalt

554 rhyolites i i i i i i i i i

-0.2 0 0.5 1.0 

INITIAL VARIMAX LOADING, «,.£

FIGURE 4. Q-mode vector diagram for the rhyolite-basalt complex on 

the Gardiner River, Yellowstone National Park, Wyo. Numbers at ends 

of vectors are the sample numbers with prefix Y.P. in table 19. Vectors 

/ and // are the varimax axes. Notes on vectors give interpretation of 

Fenner (1938). Wilcox (1944) interpreted each sample (except 904, 

1024, and 1025) as resulting from the mixture of basaltic and rhyolitic 

lavas.

Composition scores for some of the vectors represented 

in figure 4 are given in table 21. It may be seen that vectors 

outside of the positive quadrant formed by the two 

varimax axes have at least one negative score value. This is 

true for two of the samples that vary from the standard 

(1025 and 904) as well as for the two most rhyolitic samples
G28 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

TABLE 21. Composition scores, xjj, for some vectors in figure 4 

(in percent)

Vectors SiO2 A1 2O3 'FeO' MgO CaO Na2O K2O H2O+

904 46.50 17.03 12.48 9.03 12.33 2.53 -0.44 0.54 

1025 46.54 17.02 12.46 9.02 12.31 2.53 -.43 .54 

/ 49.41 16.60 11.42 8.14 11.19 2.63 .07 .54 

1024 49.60 16.57 11.34 8.08 11.12 2.64 .10 .54 

969 51.47 16.29 10.66 7.51 10.39 2.71 .43 .54 

77 76.00 12.66 1.70 .02 .81 3.62 4.70 .49 

554 76.69 12.56 1.45 -.19 .54 3.64 4.82 .49

from the complex (914 and 554). Because the negative 

score values are small, however, they can be set to zero, 

and the scores can then be recomputed to sum to 100 

without seriously affecting the validity of subsequent 

arithmetic operations.

The linearity of the compositional gradation between 

rhyolite and basalt, as recognized by Fenner and verified 

by the factor-variance diagram of figure 3 and as 

represented in Fenner' s silica variation diagrams and in the 

vector diagram of figure 4, makes it difficult to accept any 

other origin for the complex than the mixing of two 

end-members of nearly constant composition. One of these 

end-members can be taken as basalt, represented by sample 

969 (fig. 4). This is the most basaltic sample that is 

unquestionably part of the complex. The other 

end-member could be the gaseous emanations referred to 

by Fenner, but it is not likely that such emanations would 

have a constant composition that could be represented by a 

vector with positive scores at the other end of the vector 

system, and there is little doubt that the other end-member 

could be so represented. Consequently, Wilcox's 

alternative explanation that the complex formed by the 

mixing of rhyolitic and basaltic lavas seems probable, and 

the second end-member is taken as rhyolitic lava similar in 

composition to the composition scores of the vector 

representing sample 554 (fig. 4). The composition scores 

for vector 969 and the adjusted scores for vector 554 are 

given as end-member compositions and designated by 

'969' and '554' in table 22.

TABLE 22. Compositions of the end-members, f^j, for the factor model 

of the rhyolite- basalt complex on the Gardiner River (in percent)

End- SiO2 A1 2O3 FeO MgO CaO Na2O K2O H2O+ Total

'969' 51.47 16.29 10.66 7.51 10.39 2.71 0.43 0.54 100.00 

'554' 76.54 12.54 1.45 .00 .54 3.63 4.81 .49 100.00

The composition loadings for each of the 15 samples 

accepted as taken from the rhyolite-basalt complex are 

given in table 23, along with the increments (eq. 36a, 36b) 

of Na2O and H2O+ that must be added by independent 

processes to account for these constituents. For example,

mixing of 0.666 parts basaltic lava and 0.334 parts rhyolitic 

lava, with the addition of 0.12 percent Na2 O and 0.20 

percent H2 O+. When the end-members of table 22 are 

mixed in these proportions and the increments of Na2 O 

and H2 O+ are added, the total mixture recomputed to 

sum to 100 is as given in table 19B. The entire matrix of 

recomputed data in table 19B is formed by postmultiplying 

the matrix of aik in table 23 by the matrix in table 22, 

adding the increments of Na2O and H2 O+ in table 23, 

and adjusting each row of the product to sum to 100.

TABLE 23. Composition loadings, a ik , and increments of Na 2O and 

H2 O+ for the factor model of the rhyolite-basalt complex on the 

Gardiner River

Composition loadings, Increments (in percent)

Sample aik 

No Y P

End-member End-member .. _ 

,%9 , ,554 , Na20 H20 +

969 1.000 0.000 0.06 -0.10 

918 .866 .134 .04 .00 

948 .907 .094 -.04 -.04 

985 .905 .095 -.03 -.14 

512 .918 .082 -.18 -.01 

925 .666 .334 .12 .20 

514 .691 .309 -.06 -.08 

909A .562 .438 .16 -.01 

909B .446 .554 .14 .05 

543 .440 .560 -.09 .09 

467 .379 .621 -.06 .02 

980 327 .673 .26 .03 

505 .194 .806 -.08 .11 

554 .000 1.000 -.06 -.18 

914 .006 .994 -.19 -.12

The correspondence of the matrix of reproduced data in 

table 19B to the matrix of original data in table 19>1 is 

measured by the goodness-of-fit statistics given in table 24. 

Although the correspondence is good, departures from 

perfect correspondence are attributed to (1) minor errors 

present in the analytical data, (2) some small amount of 

compositional variation in the actual basaltic and rhyolitic 

end-members, and (3) other processes that may have 

caused compositional variation in the complex, such as 

contamination of the lavas by other materials and later 

alterations brought about, perhaps, by weathering.

basalt complex on the Gardiner River

[dij=xij-x..; dj*= mean of rf,y; dj* = standard deviation of djf, TJ= correlation coefficient for 

'xjj and xj;; rj= coefficient of determination]

Goodness- Variable (/)

of-fit 

statistic SiO2 A1 2O3 FeO MgO CaO Na2O K2O H2O +

dj*. ...... -0.08 -0.01 0.02 0.08 0.01 0.00 -0.02 0.00

d".* ...... .80 .23 .28 .33 .31 .00 .21 .00

rj . ....... 1.00 .98 1-00 .99 1.00 1.00 .99 1.00 

r} ..... 99 97 .99 .98 .99 1.00 .98 1.00
Q-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G29

GRANITOID INTRUSIVE IN THE 

SOUTHERN SNAKE RANGE, NEV.

An exposure of granitoid intrusive rocks in the Snake 

Creek-Williams Canyon area of the southern Snake Range 

of eastern Nevada was interpreted by Lee and Van Loenen 

(1971) to have resulted from the crystallization of a magma 

after it intruded and assimilated one or more kinds of 

sedimentary country rock. The granitoid intrusive is 

exposed over an area of 34 km 2 and varies from 

granodiorite to quartz monzonite in composition. The 

predominant minerals are quartz, microcline, plagoiclase 

(An 5 -An35 ), and biotite; hornblende and a number of 

other accessory minerals are also present. The intrusive is 

undeformed and domelike and probably has not been 

eroded to a depth of much more than 300 m. Its contacts 

with the overlying sedimentary rocks are roughly 

concordant with the bedding in the sediments at most 

localities, but exceptions that point clearly to sloping 

do occur. The eastern part of the intrusive contains 

abundant dark inclusions that Lee and Van Loenen 

referred to as xenoliths of Pioche Shale (Cambrian). They 

range in size from several centimetres in width to about 

8,000 m2 in area and are especially abundant where the 

intrusive comes into contact with the Pioche Shale. They 

tend to be ellipsoidal in shape and oriented with the long 

axis parallel to the intrusive margins. The preferred 

orientation was interpreted by Lee and Van Loenen to 

indicate that the magma flowed before solidification.

Sedimentary rocks in contact with the intrusive at the 

surface are the Osceola Argillite of Misch and Hazzard 

(1962), which is of Precambrian age and about 230 m 

thick, the Prospect Mountain Quartzite of Precambrian Z 

and Early Cambrian age, and the Pioche Shale and Pole 

Canyon Limestone of Cambrian age. The Prospect 

Mountain Quartzite is about 1,000 m thick, the Pioche 

Shale is 90 -120 m thick, and the Pole Canyon Linestone is 

about 600 m thick. The Pioche Shale contains a limestone 

unit near its base, which ranges from 1.5 to 7.5 m in 

thickness and which is an informally named unit referred 

to locally as the Wheeler linestone. None of the 

sedimentary rocks shows evidence of appreciable 

metasomatism from the intrusive, whose principal effects 

on the sediments appear to have been some deformation 

and thermal metamorphism within about 1 m of the 

contacts. The Prospect Mountain Quartzite has the 

composition of a mature shelf sediment. The Pioche Shale 

is highly variable in composition both laterally and 

stratigraphically. The Pole Canyon Limestone is a rather 

pure limestone, whereas the Wheeler limestone, a bed in 

the Pioche Shale, tends to be more argillaceous and quite 

variable from one place of occurrence to another.

Chemical data on 81 samples of intrusive rocks and 6 

samples of the dark inclusions are given in table 5 of Lee

and Van Loenen (1971), along with normative mineral 

compositions and partial modal analyses for all samples. 

However, as for the other applications, only the data on 

the nine major oxides were used for the factor analysis. Lee 

and Van Loenen have shown that the concentrations of 

most of the chemical constituents vary with the CaO 

contents of the samples. The data on the nine major 

oxides, adjusted to sum to 100 for each sample, are given 

in table 25A for the five samples with the lowest and five 

samples with the highest CaO contents.

TABLE 25. Analyses of selected samples of the granitoid intrusive in the 

Snake Range, Nev, adjusted so that nine major oxides sum to 100 

percent, and corresponding data reproduced from the factor model

[Analyses from Lee and Van Loenen (1971, table 5)]

Sample 

No. SiO> 2 A12O3 Fe2O3 FeO MgO CaO Na2O K2O H2O

A. Adjusted analyses, x« (in percent)

1

2

3

4

5

80

82

83

84

86

76.33

75.91

76.66

75.44

76.22

66.20

67.89

64.12

62.84

63.84

13.48

13.15

13.01

14.05

13.79

16.65

16.19

17.21

18.17

17.49

0.33

.60

.42

.27

.38

2.03

2.23

2.52

2.54

1.93

0.23

.24

.21

.44

.44

2.03

1.32

2.32

2.54

3.05

0.17

.12

.21

.32

.10

2.23

1.42

2.21

2.13

2.13

0.47

.52

.55

.61

.64

4.16

4.25

4.33

4.37

4.57

3.72

5.18

3.73

4.11

3.93

3.76

4.05

3.42

3.76

3.86

4.42

4.08

4.34

4.21

4.33

2.23

2.02

2.72

2.34

2.34

0.84

.19

.86

.54

.16

.70

.64

1.16

1.32

.88

B. Reproduced data (in percent)

1

2

3

4

5

80

82

83

84

86

75.98

74.98

75.95

75.27

76.26

66.75

67.35

66.61

65.35

65.69

13.29

13.32

13.30

13.52

13.40

16.78

16.46

16.78

17.03

17.05

0.50

.55

.50

.60

.51

1.96

1.84

1.97

2.10

2.09

0.25

.32

.25

.37

.27

2.10

1.94

2.11

2.28

2.26

0.11

.17

.11

.23

.13

1.84

1.70

1.85

2.01

2.00

0.82

93

.82

1.03

.85

4.04

3.77

4.06

4.35

4.33

3.72

5.18

3.73

4.11

3.93

3.76

4.05

3.24

3.76

3.86

4.49

4.36

4.49

4.32

4.49

2.06

2.25

2.04

1.80

1.83

0.84

.19

.86

.54

.16

.70

.64

1.16

1.32

.88

The factor-variance diagram derived from the oxide 

data on all 81 intrusive samples is given in figure 5, and 

resembles, to some extent, that for the rhyolite-basalt 

complex on the Gardiner River. That is, a large part of the 

variance for most constituents can be accounted for by 

models with two factors, or end-members, and the 

variances in Na2O and H2O+ call for unique factors or 

for processes that controlled their variances alone. An 

important difference from the diagram for the 

rhyolite-basalt complex, however, is the generally lower 

proportions of the total variances accounted for by 

two-factor models. The reason for this may be partly due 

to the fact that the data of Lee and Van Loenen are based 

on an analytical method acknowledged to be less precise 

than conventional methods (Shapiro and Brannock, 1956, 

p. 19), but also the reason may be partly due to variation in 

the compositions of the end-members and to other geologic 

processes that caused minor compositional variations in 

the intrusive. Nevertheless, two-factor models that include 

unique processes to account for the variabilities in
G30 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

23456789 

NUMBER OF FACTORS

FIGURE 5. Factor-variance diagram for the granitoid intrusive, 

southern Snake Range, Nev.

and H 2O+ can explain more than 70 percent of the 

variance in each constituent.

The eigenvalues of the cosine theta matrix are given in 

table 26.

TABLE 26. First nine eigenvalues of the cosine 

theta matrix for the granitoid intrusive in the 

southern Snake Range and CPN1

No.

1

2

3

4

5

6

7

8

9

Eigenvalue

67.445

10.230

1.013

.708

.572

.473

.318

.177

.064

CPN

0.8327

.9590

.9715

.9802

.9873

.9931

.9970

.9992

1.0000

Total.........81.000

1 Cumulative proportions of N.

If the interpretation of Lee and Van Loenen (1971) that 

the intrusive formed by solidification of a magma after 

assimilation of sedimentary country rocks is correct, the 

major questions raised pertain to which particular 

sedimentary units were assimilated and to what degree. 

The average compositions of the five major sedimentary

units in contact with the intrusive at the surface and of the 

dark inclusions within the intrusive are given in table 27; 

the communalities of vectors representing these 

compositions in the two-factor varimax space are also 

listed. The same vectors are shown diagrammatically in 

figure 6. It is apparent that none of the sedimentary units is 

highly compatible with the compositional series formed by 

the intrusive samples in two-factor space; the composition 

scores of the vectors representing sedimentary units in 

figure 6 depart widely from the actual compositions of the 

units.

In contrast to the sedimentary units, the average 

composition of the dark inclusions within the intrusive has 

a large communality within the two-factor space (table 27) 

and is represented by a vector close to the margin of the 

range of vectors representing samples of the intrusive (fig. 

6).

The consequences of increasing the number of factors in 

the model are examined in table 28, which gives the 

commumalities of vectors representing the dark inclusions 

and of each of the sedimentary units in factor space of two 

to nine dimensions. As indicated in the table, the 

communalities of the limestone units and the Prospect 

Mountain Quartzite remain small until the model contains 

almost as many factors as there are chemical variables. 

When m = M, any combination of the nine variables can 

be perfectly represented in the factor space.

A

Average Wheeler 

limestone of 

local usage

Average dark inclusion 

'/Intrusive sample 84 

/Intrusive sample 86

0.5

-0.3

Average Osceola Argillite of 

Misch and Hazzard (1962)

Average Pioche Shale

All vectors representing 

samples of the intrusive 

occur in this region Pole Canyon 

Limestone

Average Prospect Mountain Quartzite

-0.3 0 0.5 

INITIAL VARIMAX LOADING,a"2

+ 1.0

FIGURE 6. Q-mode vector diagram for the granitoid intrusive. Vectors / 

and // are the varimax axes. See text for explanation of vectors A and 

B.
0-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS 031

TABLE 27. Average compositions, in percent, of some sedimentary rocks and dark inclusions from the southern Snake Range and computed

communalities, h 2

Rock

Limestone ' . . . . 

Limestone '.....

Quartzite ' . . . . 

Shale 1 ........

Stratigraphic unit

Pole Canyon 

Limestone Wheeler limestone (of local 

usage).

of Misch and 

Hazzard

(1962). 

Prospect 

Mountain 

Quartzite. 

Pioche Shale ......

Si02

4.5 

31.4

62.9

91.6 

. 63.6

. 62.9

AI 203

1.1

2.9 

16.2

4.1 

15.9

16.5

Fe203

0.33 

1.6

4.3

.79

3

0

2.4

FeO

0.25

5.2 

1.9

.25 

3.3

3.0

MgO

2.0 

16.5

2.2

.23

1

0

2 3

CaO

91.0 

39.6

4.7

.01 

1.9

4.2

Na20

0.12 

.3

1.5

.13

4.0

K20

0.37 

.42

3.0

2.1 

5.3

2.1

H20 +

0.07 

.19

2.2

.38

2.7

1.0

Total

99.74 

98.11

no on

99.59

OQ CO

98 40

*2

0.218 

.436

.656

.227 

.420

QA7

1 Average composition is from Lee and Van Loenen (1971, table 9). Average compositions for the limestones were computed on a CO2-free basis.

2 Average of six analyses from Lee and Van Loenen (1971, table 5).

As shown in table 26, the average communality of the 81 

vectors representing samples of the intrusive in the 

two-factor space is 0.9590. Table 28 shows that, in the 

two-factor space, only the vector representing the average 

dark inclusion has a communality this large. Of the 

sedimentary units, only the communality for the average 

Osceola Argillite is this great, but it requires a factor space 

of 5 dimensions. The communalities for the average Pioche 

Shale are second largest.

Although the evidence presented so far points strongly 

toward a model containing the dark inclusions as one of 

the end-members, a further comparison of the suitabilities 

of the dark inclusions, the Osceola Argillite, and the 

Pioche Shale was made by use of computer simulation. 

The assumption was made that the actual mafic 

end-member in the two-dimensional compositional system 

can be represented by a vector in the position of the vector 

representing the average dark inclusion in figure 6. The 

fact that the average Osceola Argillite and the average 

Pioche Shale are not so represented is attributed, for the 

present purpose, to errors in sampling that is, the 

samples from which the averages were obtained are not 

perfectly representative of the sedimentary units. The

composition loadings of all 81 intrusive samples were then 

determined with respect to the two reference vectors that 

represent the dark inclusions and intrusive sample 1 (fig. 

6). The first matrix of simulated data, 81 by 9 in size, was 

then derived by mixing the composition of intrusive sample 

1 with randomly selected analyses of individual dark 

inclusions, given by Lee and Van Loenen (1971, table 5), in 

the proportions indicated by the composition loadings. Six 

analyses of individual inclusions were available. Two other 

simulated data matrices were formed using six analyses of 

the Osceola Argillite and nine analyses of Pioche Shale 

(Lee and Van Loenen, 1971, table 1). The factor-variance 

diagrams derived from the three simulated data matrices 

are given in figure 7. There is no doubt that the diagram 

simulating the assimilation of dark inclusions bears the 

closest resemblance to the diagram derived from the 81 

intrusive samples (fig. 5), and the results of the experiment 

are accepted as further evidence that the major part of the 

compositional variation in the intrusive can be correctly 

ascribed to the dark inclusions. Although Lee and Van 

Loenen (1971) interpreted the inclusions, apparently on the 

basis of field evidence, to be xenoliths of Pioche Shale, it is 

obvious from figure 6 that they are much closer to the

TABLE 28. Values o/h2 for the average compositions of dark inclusions and some sedimentary rocks from the southern Snake Range in models

containing two to nine factors

[Line in body of table separates values of h? above and below the value of 0.9590 (table 26)]

Number of factors in model, m

Rock Stratigraphic unit 5

Shale ...............

Quartzite ............

Limestone ...........

Limestone ...........

and Hazzard (1962) ........

...... Pioche Shale ................

...... Prospect Mountain Quartzite . .

...... Pole Canyon Limestone. ......

0.967

.656 

.420

.227

.436 

.218

0.968

.895 

.792 

.379

.469

.237

0.969

.903 

.816 

.422

.471 

.237

0.977 

.964

.926 

.453

.471 

.292

0.985 

.989

.939 

.459

.586 

.292

0.987 

.989

.941 

.890

.669 

.293

0.992 

.996

.942 

.943

.852 

.686

.000

.000 

.000 

.000

.000 

.000
G32 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

234 

NUMBER OF FACTORS

2 3 

NUMBER OF FACTORS

2 3 

NUMBER OF FACTORS

FIGURE 7. Factor-variance diagrams derived from data obtained by mathematical mixing of sample 1 of the granitoid intrusive and randomly 

selected samples of (A) dark inclusions from the intrusive, (B) Osceola Argillite, and (Q Pioche Shale.

Osceola Argillite in bulk chemical composition. The source 

of the inclusions, therefore, appears to be unknown; 

judging from Lee and Van Loenen's (1971) petrographic 

descriptions, even an igneous origin is possible.

If the inclusions are igneous, they may represent 

segregations that were separated from the magma that 

formed the intrusive. The composition loadings (mixing 

proportions) for at least some of the samples, in this 

situation, would be negative with respect to the 

end-member representing the inclusions. All that can be 

said of the composition of the magma is that it would be 

represented by a vector somewhere within or close to the 

plane of sample vectors represented in figure 6. This 

vector, and the magma composition, could be fixed more 

closely only if some assumption were made as to which, if 

any, samples formed by incorporation of inclusions rather 

than separation of inclusions from the magma. However, 

we shall proceed with the assumption that all mixing

proportions are positive, as is necessarily the situation if 

Lee and Van Loenen (1971) are correct in their 

interpretation that the inclusions are of sedimentary origin. 

The composition scores for the vectors represented in 

figure 6 are given in table 29. It may be seen that vectors 

lying outside the positive region formed by vectors A and B 

have composition scores that are partly negative and that 

all the vectors representing samples of the intrusive occur 

in the positive quadrant formed by the varimax axes. If it is 

assumed that the compositional variation in the intrusive 

has been due mostly to the mixing of two end-members in 

positive proportions rather than to the separation of some 

constituent, one of the end-members must have been of a 

composition that can be represented by a vector between 

that of intrusive sample 1 and vector B (fig. 6). The other 

end-member must have been of a composition that could 

be represented by a vector somewhere between vector A 

and the vector representing sample 84. The composition
g-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G33

TABLE 29. Composition scores, xy, of some vectors in figure 6 (in 

percent)

Vector SiO2 A1 2O3 Fe2O3 FeO MgO CaO Na2O K2O H2O +

Average Pole 

Canyon 

Limestone. 51.19 22.02 4.24 4.99 4.56 9.07 3.94 -1.81 1.80 

Average 

Wheeler 

Limestone. 53.60 21.20 3.88 4.53 4.13 8.28 3.91 -1.20 1.67 

A...... 58.39 19.55 3.17 3.63 3.29 6.71 3.85 .00 1.41

Average dark 

inclusion.. 65.17 17.22 2.17 2.36 2.09 4.50 3.76 1.70 1.03 

Intrusive 

sample84. 65.57 17.09 2.11 2.29 2.02 4.36 3.75 1.80 1.01 

Intrusive 

sample86. 65.69 17.05 2.09 2.26 2.00 4.33 3.75 1.83 1.00 

Average 

Osceola 

Argillite .. 65.90 16.98 2.06 2.22 1.97 4.26 3.75 1.88 .99 

Average 

Pioche 

Shale..... 69.39 15.78 1.54 1.57 1.35 3.11 3.70 2.76 .80 

Intrusive 

sample?.. 76.37 13.38 .51 .26 .12 .83 3.61 4.51 .41 

Intrusive 

samplel.. 76.42 13.37 .50 .25 .11 .81 3.60 4.52 .41 

B...... 77.05 13.15 .41 .13 .00 .61 3.60 4.68 .38 

II...... 77.56 12.97 .33 .04 -.09 .44 3.59 4.81 .35 

Average 

Prospect 

Mountain 

Quartzite 86.67 9.85 -1.02 -1.67 -1.70 -2.54 3.47 7.10 -.15

loadings for the 81 intrusive samples will not vary greatly 

with any choice of possible end-members within these 

limits. However, inasmuch as the dark inclusions and 

intrusive sample 1 are known to exist, the vectors 

representing these materials were used as reference vectors. 

That is, the intrusive is interpreted as having formed by 

assimilation of material similar in composition to the 

average dark inclusion into a magma similar in 

composition to intrusive sample 1 . The composition scores 

of the reference vectors are brought together in table 30, 

where they are interpreted as end-member compositions.

TABLE 30. Compositions of the end-members for the Snake 

Range factor model (in percent)

End-member SiO2 A1 2C>3 Fe2C>3 FeO MgO

Magma........ 76.42 13.37 0.50 0.25 0.11

Inclusions 65 - 17 17 -22 2 - 17 2 - 36 2 -09

End-member CaO Na2O K2O H2O+ Total

Magma ....... 0.81 3.60 4.52 0.41 99.99

Inclusions..... 4.50 3.76 1.70 1.03 100.00

The composition loadings for the 81 intrusive samples with 

respect to these end-members are given in table 31, along 

with the increments (eq 36a, 36b) of Na2O and H2O + 

that must have been added to each sample by independent 

geochemical processes in order to account for the 

variabilities in these constituents. Postmultiplication of the 

matrix of loadings in table 31 by the matrix of scores in 

table 30. addition of the increments of Na->O and H-.O +

TABLE 3 1 . Composition loadings and increments of Na 2O and H2O + 

for the Snake Range factor model

a . Increments 

'* (in percent)

Magma Inclu- Na2O H2O + 

sions

1 1.000 0.000 0.13 0.44 

2 .966 .034 1.64 -.24 

3 .999 .001 .14 .45 

4 .940 .060 -.07 -.09 

5 .991 .009 .32 -.25

6 .944 .056 .12 .44 

7 .9% .004 .11 -.05 

8 .998 .002 -.09 .04 

9 .928 .072 .10 .03 

10 .952 .048 -.01 -.05

11 .970 .030 -.12 -.04 

12 .909 .091 -34 -.33 

13 .881 .119 1.25 .18 

14 .892 .108 -.41 -.30 

15 .764 .237 -.03 .17

16 .806 .194 -.23 -.01 

17 .793 .207 -.23 -.11 

18 .793 .207 -.35 -.28 

19 .786 .215 -.09 .06 

20 .739 .261 .18 .01

21 .763 .237 -.22 .00 

22 .686 .314 -.05 .19 

23 .794 .206 -.22 -.08 

24 .647 .353 -.03 .22 

25 .813- .187 -.02 -.28

26 .780 .221 -.23 -.20 

27 .651 .349 .16 -.23 

28 .658 .342 -.14 .04 

29 .670 .330 -.45 .07 

30 .659 .341 -.04 -.08

31 .630 .370 -.01 .30 

32 .635 .365 .07 .18 

33 .613 .387 -.16 _.06 

34 .568 .432 -.06 _ oi 

35 .580 .420 -.03 .19

36 .637 .363 -.35 .26 

37 .558 .442 .38 .28 

38 .596 .404 .60 -.15 

39 .573 .427 -.02 .44 

40 .448 .552 -.07 .12 

41 .563 .437 -.05 _.15

a .. Increments 

' (in percent)

Sample 

Magma Inclu- NajO H2O + 

sions

42 0.593 0.407 -0.58 0.17 

43 .497 .503 -.26 .22 

44 .593 .407 .57 -.27 

45 .540 .460 .04 -.15 

46 .672 .328 -.24 -.38

47 .5% .404 -.48 .15 

48 .532 .468 .06 .15 

49 .564 .436 -.15 .20 

50 .516 .484 .07 -.07 

51 .583 .417 -.16 -.22

52 .602 .398 -.15 .01 

53 .634 .367 -.36 .02 

54 .589 .411 -.05 -.15 

55 .654 .346 -.24 -.31 

56 .545 .455 -.03 .19

57 .602 .399 .37 .13 

58 .371 .629 -.06 -.43 

59 .604 .3% -.31 -.11 

60 .423 .577 .25 .18 

61 .256 .744 .11 -.24

62 .520 .480 -.37 .04 

63 .443 .557 .05 .15 

64 .364 .636 .02 -.37 

65 .368 .632 .64 -.20 

66 .413 .587 .15 .15

67 .505 .495 .16 .26 

68 .378 .623 .15 -.25 

69 .254 .746 .46 1.43 

70 .327 .673 -.05 -.15 

73 .251 .750 -.10 -.22

74 .237 .763 .01 -.13 

75 .243 .757 -.11 .00 

76 .170 .830 .12 -.30 

78 .129 .871 .33 .47 

79 .125 .875 .12 -.17

80 .127 .873 .01 -.25 

82 .197 .803 .32 -.27 

83 .121 .879 -32 .20 

84 .036 .964 .01 .32 

86 .046 .954 ,\\ -.12

to the product matrix and, then, adjustment of each row of 

the product matrix to sum to 100 yield the matrix of 

reproduced data partially given in table 25B. 

The correspondence between the original data and the 

reproduced data, given partially in table 25, is measured by 

the goodness-of-fit statistics given in table 32. As before, 

the lack of perfect correspondence is attributed to 

analytical errors, variation in the compositions of the 

actual end-members, and other geologic and chemical 

processes that caused minor compositional variations in 

the intrusive. It seems likely that these processes included

TABLE 32. Goodness-of-fit statistics for the Snake Range factor model 

[dy = Xjj-Xjj; dr- mean of dy, d~* = standard deviation of rf,y; TJ = correlation coefficient 

for ^ . and Xjj; r? = coefficient of determination]

Good- Variable (/)

statistic Si°2 A1 2°3 Fe2°3 Fe° M«° CaO Na2O K2O H2O +

dj ...... 0.06 -0.02 0.00 -0.02

r Qf\ Q1 ft^ ftft

/ /....... .93 .87 .72 .78

-0.01 0.01 0.00 -0.01 0.00 

.22 .23 .00 .28 .00 

.92 .97 1.00 .93 1.00 

.85 .94 1.00 .86 1.00
G34 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

the sloping and assimilation of sedimentary country rock 

as reported by Lee and Van Loenen (1971). The major 

variation in the intrusive, however, seems certain to be 

related to the dark inclusions of unknown source and 

origin.

LAVAS AND PUMICES

OF THE 1959 SUMMIT ERUPTION

AT KILAUEA, HAWAII

On November 14, 1959, Kilauea Volcano on the Island 

of Hawaii began a series of eruptions from Kilauea Iki near 

its summit, that consisted of 17 phases and lasted until the 

20th of December. About 3 weeks later another series of 

eruptions began on the flank of the volcano and continued 

almost without interruption for a period of more than a 

month. The volcanic activity and the detailed chemistry, 

mineralogy, and petrology of the volcanic products have 

been described by Murata and Richter (1966), Richter and 

Moore (1966), Murata (1966), Richter and Murata (1966), 

and Richter, Eaton, Murata, Ault, and Krivoy (1970). The 

volcanic products consisted principally of lavas and 

pumices of tholeiitic basalt. The less siliceous samples of 

the lavas contain phenocrysts of olivine (Fa 13) only, and 

the more siliceous samples contain clots of clinopyroxene 

and plagioclase, along with a more ferrous olivine (Fa 18 ) 

(Richter and Murata, 1966, p. D3 D4). The groundmass 

portions of the lavas are either glassy or consist of 

fine-grained plagioclase, clinopyroxene, and opaque 

minerals.

Murata and Richter (1966, p. Al) summarized their 

interpretation of the origin of the summit lavas as follows:

The hottest and most mafic lavas were produced in the 1959 summit 

phase of the eruption. The compositional variation (7 to 19 percent MgO) 

among most summit lavas is ascribable solely to fractional crystallization 

of olivine (Fa 13 ), the first silicate mineral to separate from the primitive 

magma. A direct relationship between the olivine content of summit lavas 

and the rate of lava discharge suggests that strong currents of magma 

erode beds of previously sedimented olivine crystals lying on the bottom 

of the magma chamber. The coolest, olivine-poor summit lavas contained 

some small phenocrysts of clinopyroxene and plagioclase as well as 

olivine. The compositions of these lavas indicated a slight overall 

accumulation of clinopyroxene, the second mineral to start separating 

from the cooling magma.

A mathematical study of the compositional variations in 

the lavas and pumices from the summit eruption was made 

by Wright (1973), who derived a model that accounts for 

nearly all of the variation in 11 oxide constituents by the 

mixture of two magmas and addition of chromite and 

olivine of varying composition (Fa9 4 to Fa 13.9 , or Fa 13.1 

to Fa 18.9 in weight percents as given by Wright). Wright's 

model, therefore, contained five end-members: (1) magma 

A, (2) magma B, (3) an impure chromite, (4) a magnesiumrich olivine (Fajo), and (5) a less magnesium-rich olivine 

(Fa30). One possible difficulty with the model is that it

calls for the addition of olivine and chromite to the parent 

magmas for all samples, rather than the separation of these 

minerals by fractional crystallization and settling as might 

be expected to occur in the upper parts of a magma 

chamber. As noted in the preceding paragraph, Murata 

and Richter (1966) regard the separation of olivine as the 

dominant process that caused compositional variation 

among the lavas. Another possible difficulty with Wright's 

model is that, although it includes a chromite end-member, 

it does not include the more prominent minerals, 

clinopyroxene and plagioclase.

The factor-variance diagram for the lavas and pumices 

of the summit eruption (fig. 8) was derived from the 

analytical data on the same samples as those used by 

Wright (1973, table 1) which, except for four samples, are 

from Murata and Richter (1966, table 1). Because of the 

apparent chemical distinctiveness of samples S -1 and 

S - 2, as noted by Wright (1973, p. 850), the data on these 

samples were not used in the derivation of the 

factor-variance diagram, but they were later tested for 

their conformity with the compositional series formed by 

the other 22 samples. Wright (1973, table 5) also withheld 

the data on these samples from the data to be explained by 

the model and used them to represent the compositions of 

his two end-member magmas. The analyses used in the 

derivation of figure 8 are given in table 33A.

34567 

NUMBER OF FACTORS

FIGURE 8. Factor-variance diagram for the lavas and pumices of the 

1959 summit eruption at Kilauea, Hawaii.
g-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G35

TABLE 33. Analyses of lavas and pumices from the 1959 summit 

eruption at Kilauea, adjusted so that nine major oxides sum to 100 percent, and corresponding data reproduced from the factor model

[Analyses from Murata and Richter (1966) and Wright (1973))

Sample SiO2 A12O3 FeaO;J FeO MgO CaO NaaO K2O H2O +

A. Adjusted analyses, x - (in percent)

S-3

S-3A

S-4

S-5

S-7

S-8

S-9

Iki 51

S-10

S-ll

S-12

S-13

S-14

S-14A

S-15

S-19

S-20

S-21

S-22

S-24

S-25

S-25A

S-3

S-3A

S-4

S-5

S-7

S-8

S-9

Iki 51

S-10

S-ll

S-12

S-13

S-14

S-14A

S-15

S-19

S-20

S-21

S-22

S-24

S-25

S-25 A

51.16

51.18

50.42

47.77

49.58

50.77

49.78

51.08

50.03

49.43

51.29

51.04

51.07

51.00

47.95

48.71

49.00

47.88

48.27

48.17

51.08

51.00

51.76

50.84

50.26

47.73

49.74

50.69

49.68

50.68

50.02

49.35

51.42

50.67

51.82

51.31

47.91

48.64

49.05

47.95

48.24

48.17

51.35

51.28

13.22

13.18

12.30

9.74

11.89

12.84

11.95

12.97

12.46

11.73

13.79

13.65

13.74

13.64

10.06

10.77

11.27

9.94

10.12

10.31

13.57

13.58

14.08

13.20

12.46

9.74

11.80

12.85

11.84

12.98

12.21

11.54

13.78

13.08

14.48

13.69

9.94

10.79

11.26

10.00

10.24

10.26

13.66

13.52

2.78

1.88

1.51

1.61

1.40

1.34

1.40

1.66

1.44

1.47

2.38

1.81

4.39

2.74

1.40

1.74

1.44

.39

1.27

1.85

2.42

2.00

B.

2.08

2.38

.72

.33

.23

.43

.61

2.15

1.74

1.83

2.24

2.54

3.63

2.36

.36

.72

.85

.40

.11

.57

2.06

1.72

9.09

9.96

10.37

10.54

10.60

10.51

10.58

10.26

10.76

10.57

9.55

10.04

7.77

9.39

10.78

10.37

10.62

10.70

10.81

10.29

9.58

9.96

9.12

9.18

11.84

19.98

14.06

10.80

13.72

9.62

11.96

14.51

8.51

9.21

8.83

8.98

19.32

17.07

15.81

19.72

18.75

18.71

9.14

9.24

11.85

11.87

10.91

8.45

9.95

10.93

10.05

11.70

10.72

9.87

11.55

11.41

11.39

11.38

8.51

9.18

9.48

8.36

8.68

8.64

11.34

11.34

2.18

2.21

2.11

1.56

1.99

2.16

2.00

2.18

2.09

1.95

2.30

2.27

2.26

2.27

1.61

1.74

1.87

1.59

1.65

1.65

2.25

2.22

0.52

.52

.48

.36

.46

.56

.48

.51

.49

.46

.57

.56

.55

.54

.37

.42

.51

.42

.41

.39

.55

.56

0.08

.02

.05

.01

.07

.09

.04

.03

.04

.01

.06

.00

.00

.05

.00

.00

.00

.00

.03

.00

.07

.09

Reproduced data (in percent)

9.84

9.64

10.28

10.82

10.76

10.50

10.42

9.86

10.27

10.24

9.71

9.49

8.43

9.62

10.79

10.40

10.24

10.74

10.98

10.58

9.88

10.19

7.20

9'.95

12.02

20.20

13.86

10.68

13.85

10.55

12.76

14.85

8.11

10.34

6.44

8.45

19.61

17.13

15.69

19.40

18.60

18.69

8.44

8.75

12.10

11.27

10.67

8.26

10.13

11.06

10.12

11.10

10.45

9.83

11.81

11.15

12.30

11.72

8.43

9.16

9.57

8.49

8.73

8.70

11.72

11.63

2.37

2.19

2.07

1.57

1.96

2.16

1.%

2.16

2.02

1.89

2.31

2.17

2.40

2.29

1.60

1.75

1.84

1.61

1.67

1.66

2.29

2.27

0.52

.52

.48

.36

.46

.56

.48

.51

.49

.46

.57

.56

.55

.54

.37

.42

.51

.42

.41

.39

.55

.56

0.06

.02

.04

-.01

.06

.08

.03

.03

.03

.01

.04

.00

-.04

.03

-.01

-.01

.00

 .01

.02

 .01

.05

.08

The factor-variance diagram in figure 8 shows that the 

compositional variation in the lavas and pumices of the 

1959 summit eruption at Kilauea can be almost completely 

accounted for by the mixing, or unmixing, of five factors, 

or end-members. This is in accord with the fact that 

Wright's model, containing five end-members, accounts 

for the observed data almost perfectly. However, the 

diagram shows further that the final two factors are 

required mainly to account for about 15 percent of the 

variance in K2O, 23 percent of the variance in FeO, and 29 

percent of the variance in Fe2O3 . The FeO and Fe2O3 

contents of the lavas and pumices have undoubtedly been 

modified by oxidation after they reached the atmosphere, 

and no attempt will be made to account for this process in 

the derived model. This leaves K2O as the only other 

variable that a three-factor model would not explain

almost entirely. Therefore, a three-factor (end-member) 

model will be derived and increments of K^O will be added 

to, or subtracted from, each sample, as required (eq 35), to 

account for the additional 15 percent of the variation in 

this constituent.

The eigenvalues of the cosine theta matrix for the data in 

table 33^4 are given in table 34 and show that the 

three-factor models will account for 98.56 percent of the 

total variability, as a sum of squares, in the normalized 

data matrix. They also indicate that the average sample 

communality in the three-factor space is 0.9856.

TABLE 34. First nine eigenvalues of the cosine 

theta matrix for the lavas and pumices of the 1959 

summit eruption at Kilauea and CPN^

No.

1

2

3

4

5

6

7

8

9

Total....

Eigenvalue

16.868

4.347

.469

.209

.086

.016

.003

.002

.001

......22.001

CPN

0.7667

.9643

.9856

.9951

.9990

.9997

.9999

1.0000

1.0000

'Cumulative proportions of N.

The compositions of samples S - 1 and S - 2 (table 35^4), 

used as end-member compositions by Wright (1973), were 

tested by the method illustrated in parts I and II and were 

found to have communalities in the three-factor space of 

0.9710 and 0.9460, respectively. These values indicate that 

sample S - 1 conforms with the compositional series of the 

summit lavas almost as well as the average sample in the 

series; sample S-2 does not conform as well. The 

composition scores for the vectors representing these 

samples in the three-factor space are given in table 355 

and, except for iron and MgO, are close to the actual 

sample compositions.

TABLE ^. Compositions of two samples of spatter from the 1959 

summit eruption at Kilauea, adjusted so that nine major oxides sum to 

100 percent, and composition scores of vectors representing these 

compositions in the three-factor Kilauea model

[Analyses from Murata and Richter (1966)]

Sample No. Si02 A12O3 Fe2O3 FeO MgO CaO Na2O K2O H2O +

A. Adjusted compositions, Jt(y, iu percent

S-l

S-2

51.57 

51.70

13.11 2.64 9.15 8.34 12.32 

14.15 1.44 10.33 7.47 11.93

2.20 

2.38

0.57 

.62

0.09 

.00

B. Composition scores,?;;, in percent

S-l 

S-2

51.93 

50.94

14.30 2.06 9.82 6.52 12.28 

13.37 2.56 9.43 9.51 11.39

2.42 

2.22

.58 

.54

.09 

.02

A vector diagram for the three-factor solution is given in 

figure 9A. The diagram was constructed after restoring all 

represented vectors to unit length. Interpretation of the
G36 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

A. VARIMAX AXIS/ IS VERTICAL ENTIRE VECTOR SYSTEM HAS BEEN ROTATED 

TO THE LEFT ABOUT VARIMAX AXIS///

FIGURE 9. Q-mode vector diagrams for lavas and pumices of the 1959 summit eruption at Kilauea, Hawaii. Solid dots with Roman numerals 

represent the three varimax axes. Solid dots with sample numbers represent sample vectors. Open circles represent vectors discussed in text. Shaded 

areas indicate the portions of the vector space where vectors have all nonnegative composition scores. All points have been projected from an 

upper hemisphere in a direction perpendicular to the plane of the diagram.

diagram is similar to interpretation of the stereograms used 

in structural geology. Points in the diagram represent the 

intersections of the vectors with the upper surface of a 

sphere formed by the varimax axes (/, //, and /// in fig. 9). 

The shaded areas identify the parts of the sphere wherein 

all contained vectors have composition scores that are 

entirely nonnegative. All vectors outside of these areas 

have one or more composition score values that are 

negative. Thus, any composition represented by a vector 

far outside of these areas must be regarded as incompatible 

with the compositional series formed by the 22 lava and 

pumice samples and not possibly that of one of the three 

end-members for the factor model.

The small shaded area near the negative end of varimax 

axis // in figure 9A reflects a continuation of the larger 

shaded area at the positive end of varimax axis // on the 

lower half of the sphere. As pointed out in part I, colinear 

vectors have identical composition scores. In order to show 

the area of nonnegative composition scores in a contiguous 

form, the entire vector system was mathematically rotated 

60 degrees to the left about varimax axis ///. The rotated 

system is shown in figure 9B.

The vector configuration in figure 9B shows the 

compositional series of the lavas and pumices as it varies in 

three dimensions. The major dimension roughly parallels

the plane of varimax axes / and //, and the lesser dimension 

is across this plane. The general nature of the 

compositional variation in each direction can be seen by 

inspection of the composition scores for the three varimax 

axes given in table 36.

The more siliceous samples are represented by vectors 

near varimax axis / and the more mafic samples by vectors 

near axis //. The composition of sample S-5 (table 33), 

for example, is nearly identical to the composition scores 

for axis //(table 36). Variation across the plane of varimax 

axes /and //is less pronounced than variation along it and 

appears to be caused by relatively minor variation in 

certain constituents, especially Fe2O3 , FeO, and H2O+. 

(See varimax axis /// in table 26.)

Even though the varimax axes serve as an adequate 

reference system for describing the sample vector 

configuration, there is no reason to believe that their 

composition scores even approximate the compositions of 

any materials actually involved in the magmatic 

differentiation. Alternative reference vectors having more 

probable composition scores must be sought.

In view of the observation of Murata and Richter, 

(p. G34), regarding the importance of fractional 

crystallization of olivine (Fa 13) in determining the 

compositional variation in the 1959 summit lavas, olivine
Q-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G37

TABLE 36. Initial loadings, a^, and composition scores, x^, for selected vectors in the varimax solution for Kilauea, and compositions of some

volcanic materials

[Leaders (...) indicate no data]

Vector

j

j[

HI

A ........

B ........

C ........

Material

Initial loadings, a'/.

1.0

.0

.0

.0

. -.4146

.0

0.0 

1.0 

.0 

.9682 

.8750 

.8660

Nature

0.0 

.0 

1.0 

-.2500 

-.2500 

-.5000

Composition scores, Xji, in percent

SiO2

A.

52.80 

47.66 

49.26

47.52 

45.94 

47.31

Si02

A1203 Fe203 FeO MgO

Vectors in the varimax solution

15.30 

9.68 

12.59 

9.42 

7.65 

9.03

A1203

2.46 

1.38 

7.21 

.82 

.25 

.03

Fe203

9.39 

10.81 

5.32 

11.30 

11.94 

12.03

FeO

3.61 

20.39 

13.43 

21.01 

26.25 

21.94

Compositions,

MgO

CaO Na2O

2

13.14 ;

8.20 

10.18 

8.02 

6.48

7.75

5.60

.55 

.90 

.52 

.20

.47

in percent

CaO Na2O

K2O

0.63 

.38 

.48 

.37 

.29 

.36

K20

H2O +

0.09 

-.01 

-.36 

.02 

.01 

.07

H2O +

B. Composition of volcnnic materials

a. ........ . Average oceaniteof \lacdonald 46.4 8.5 2.5 9.8 20.8 7.4 1.6 0.3

(1968, p. 502) 

b .......... Crystal-rich fraction separated from lava.

Contaminated with glass (Murata and 

Richter, 1966, table 6, column 1)..... 45.76 '7.59 12.20 24.50 6.59 1.22 .28

'includes Cr2O3 .

must be considered as one possible end-member in the 

three-factor model. The complete series from fayalite to 

forsterite was tested at increments of one molar percent for 

compatibility with two- to five-factor solutions; the 

resultant communalities are shown in figure 10 and 

indicate that Fa 15 is the most compatible with the 

compositional system formed by the lava and pumice 

samples as represented in a three-factor solution. The 

computed communality, however, is only 0.9834, and four 

of the nine composition score values for the representative 

vector are negative. The positions of vectors representing 

olivines that range from Fa 0 to Faioo are shown in figure 

9; all the magnesium-rich varieties lie outside of the areas 

of nonnegative composition scores. This and the fact that 

the iron-rich varieties all have low communalities (fig. 10) 

are interpreted as indications that the addition or 

separation of any pure olivine alone cannot help to account 

for compositional variations in the lavas and pumices of 

the 1959 summit eruption. Rather, magnesium-rich olivine 

may have separated from the parent magma along with 

some other phases. If the amounts of olivine and of other 

phases that were separated were highly correlated, their 

combined composition might represent, in effect, that of a 

single end-member.

The composition scores for vector B in figure 9 were 

computed because this vector lies near one extreme of the 

elongate area of nonnegative composition scores and 

because it occurs relatively close to the vector representing 

the most compatible olivine composition (Fa, 5 ). The 

composition scores for vector B are given in table 36, and 

the normative mineral composition is given in table 37. The

'0 10 20 30 40 50 60 70 80 90 100 

FAYALITE CONTENT OF OLIVINE, IN MOLAR PERCENT

FIGURE 10 Communalities of olivines in two- to five-factor varimax 

solutions for the lavas and pumices of the 1959 summit eruption at 

Kilauea, Hawaii.

composition scores are closely similar to the composition 

of crystalline material separated from one of the flank 

lavas by Murata and Richter (1966, p. A20) which
G38 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

TABLE 37. Chemical compositions, fk ., and norms of the end-members for the Kilauea factor model

(in percent)

End-member SiO2 A12O3 Fe2O3 FeO MgO CaO Na2O K2O H2O+ Total

A. Compositions, fkj

A ......

B

C ......

Endmember

47.52 

45.94 

47.31

Or

9,

7. 

9.

Ab

.42 0.82 11.30 

.65 .25 11.94 

.03 .03 12.03

Wo

Di

En

21.01 

26.25 

21.94

Fs

8.02 

6.48

7.75

Hy

En

1.52

1.20 

1.47

Fs

0.37 0.02 

.29 .01 

.36 .07

01 

 - Mt

Fo Fa

100.00 

100.01 

99.99

Total

B. Norms

A ......

B ......

C ......

2.19 

1.71 

2.13

12.86 

10.15 

12.44

17.79 9.20 

14.63 7.31 

16.98 8.97

6.15

5.05 

5.93

2.36 

1.68 

2.39

5.75 

3.59 

4.39

2.21 

1.19

1.77

28.33 11.98 1.19 

39.76 14.56 .36 

31.06 13.83 .04

100.01 

99.99 

99.93

contained plagioclase, clinopyroxene, olivine (Fa J4 8 ), and 

contaminations of glass. The composition of the crystalline 

material is given in table 36B. The composition of the 

normative olivine derived from the composition scores for 

vector B is Fa 2o (table 37), about the same as the most 

magnesium-rich normative olivines computed for the 

summit lavas. (See Murata and Richter, 1966, p. A5.) 

Vector B will be accepted as representing one of the three 

end-members for the factor model. The end-member is 

interpreted to have been a mixture that tended to consist of 

about one-half olivine plus pyroxene and plagioclase in 

about equal amounts.

The composition scores were computed and were 

examined for a series of vectors located along the 

lowermost margin of the shaded area of figure 9, in search 

of a vector that might represent the glomeroporphyritic 

clots of clinopyroxene and plagioclase, which are 

associated with a relatively iron-rich olivine (Fajg), as 

described by Richter and Murata (1966, p. D4). Such clots 

were also described by White (1966, p. 270, 308). It was 

found that vector C, as represented in figure 9, has 

composition scores (table 36) that yield a normative olivine 

(table 37) of Fa 24. The normative composition also shows 

plagioclase, pyroxene, and about 2 percent orthoclase. 

This vector was taken as representing the second of the 

three end-members for the factor model. The principal 

reason for selecting it was that the normative olivine 

computed from the composition scores is a few molar 

percent richer in fayalite than the normative olivine of 

vector B; thus, the first two end-members, although 

hypothetical and of unknown actual mineral composition, 

have normative compositions that are at least roughly in 

accord with the principal mineral assemblages described by 

Richter and Murata (1966). These are (1) a less fayalitic 

olivine (Fa 13) that is not in equilibrium with the enclosing 

glass (vector B), and (2) a slightly more-fayalitic olivine 

(Faj 8 ) associated with clots of plagioclase and 

clinopyroxene (vector C). The only major discrepancy 

apparent is that the first of these end-members

(represented by vector B) contains normative plagioclase 

and pyroxene as well as olivine. It is possible, however, 

that any plagioclase and pyroxene that may have actually 

separated from the magma along with the olivine either 

could be present in the fine-grained groundmasses of the 

lavas or may have settled into a hotter magma and been 

partially or wholly resorbed.

The third and final end-member for the factor model 

will represent the parent magma, and its selection will 

determine the composition loadings for each of the 24 lava 

and pumice samples that is, the proportions of each 

end-member that must be combined with or extracted from 

the others to approximate the composition of each sample. 

A positive loading will indicate the net addition of the 

end-member to the parent magma, even though it may also 

have been precipitated and separated from the parent to 

some degree. A negative loading will indicate a net loss of 

the end-member from the parent magma, even though it 

may also have been added to some lesser degree. This kind 

of interpretation is necessary in order to reconcile the 

model with the concept of a magma chamber that is both 

precipitating mineral phases at least throughout its upper 

part, losing these phases by gravitational settling, and 

receiving precipitated phases from above.

According to the petrographic descriptions of the lava 

and pumice samples that were given by Richter and Murata 

(1966, p. D3 - D4), the more fayalitic olivine (Fa 18 ) and 

clots of plagioclase and pyroxene were found only in 

samples S-l, S-3, and S-25; they were not found in 

samples S - 9 or S - 20. Therefore, a parent magma must 

be selected such that the composition loadings on vector C 

(representing the more fayalitic olivine, plagioclase, and 

pyroxene) will be positive for samples S-l, S-3, and 

S-25 and zero or negative for samples S - 9 and S - 20. 

This requires that the composition of the parent magma be 

such that it can be represented by a vector that lies in the 

parent magma plane indicated in figure 9. The plane was 

drawn through vector B so that it separates the vectors 

representing samples S-25 and S -9.
Q-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS 

TABLE 38. Composition loadings, a ik , and increments o/K2O/or the Kilauea factor model

G39

Sample No. -

S 1

S 2

S 3

S 3A

S 4

S 5

S 7

S 8

S-9

Iki 51

S-IO

S 11

A

3.386

3.806

3.344

3.532

2.622

1.660

1.929

2.355

2.380

3.228

2.606

2.585

«ik

B

-2.850

-2.064

-2.701

-2.020

-1.738

-.045

-1.461

-2.107

-1.363

-1.951

-1.569

-1.087

C

0.464

-.741

.357

-.512

.116

-.616

.532

.753

-.017

-.276

-.037

-.499

Increment of

K

f\

(in percent)

-0.02

.08

-.07

-.02

-.02

-.02

-.02

.03

.00

-.02

.00

.00

Sample No.

S 12

S-13

S 14

S 14A

S 15

S 19

S 20

S 21

S 22

S 24

S 25

S 25A

A

3.481

3.708

5.246

3.603

1.726

2.316

2.553

1.787

1.499

2.032

3.250

2.822

"ik

B

-2.461

-1.890

-2.465

-2.353

-.170

-.617

-.898

-.201

-.452

-.315

-2.432

-2.452

C

-0.019

-.818

-1.781

-.250

-.556

-.699

-.655

-.587

-.046

-.717

.182

.630

Increment of

K2O

(in percent)

0.00

.03

-.04

-.02

-.02

.00

.07

.03

.00

-.01

-.02

.00

Selection of the vector in the parent magma plane that 

may represent the composition of the parent magma is 

difficult, but I assumed that the lavas and pumices were 

erupted from the upper parts of the magma chamber and 

that each sample formed from parent magma that had 

undergone at least some net loss in the less fayalitic olivine. 

Thus, the parent magma vector must be to the right 

(toward the positive end of varimax axis If) of a plane 

through vector C and the vector representing sample S - 5. 

Vector A (fig. 9) was selected for this reason, but its 

composition scores (table 36) are, for the most part, 

strikingly similar to the average composition of oceanites 

as given by Macdonald (1968, p. 502). The only notable 

exceptions are the score values for Fe2O3 and FeO 

(compare vector A and material "a" in table 36), but the 

total iron values are in almost perfect accord. Macdonald 

(1968, p. 511) judged oceanite to be the most likely single 

parent for the Hawaiian lavas.

The composition loadings for all samples with respect to 

reference vectors A, B, and C in figure 9 are given in table 

38; also included are the increments of K^O that must be 

added to each sample to account for the additional 15 

percent of the variability in this constituent. Sample S - 1, 

for example, is viewed as having originated from 3.386 

parts original magma (vector A) from which 2.850 parts of 

the less fayalitic olivine as well as plagioclase and pyroxene 

(vector B) have been subtracted and 0.464 parts of the 

more fayalitic olivine plus plagioclase and pyroxene (vector 

Q have been added. The sum of these three loadings is 

unity; the source of the 3.386 parts original magma, 

presumably, was a lower part of the magma chamber 

where magma was displaced upwards by settling crystalline 

phases from above. The magnitudes of the loadings on the 

vector representing the parent magma (vector A in table 

38) indicate that the model requires a magma chamber at 

least several times the volume of the erupted materials.

Combination of the three end-members represented in 

table 37 in the proportions indicated in table 38, addition 

of the corresponding increment of K2O, and minor 

adjustment of each analysis to sum to 100 yield the

Goodnessof-fit 

statistic

d}.. ......

d" .......

r. ....... .

J 

r* ........

J

Variable (/)

SiO2

0.04 

.28 

.98 

.96

A1 203

0.03 

.29 

.98 

.96

Fe203

-0.01 

.39

.84 

.71

FeO

0.02 

.35 

.88

.77

MgO

-0.11 

.82 

.98

.97

CaO

0.04 

.33 

.97 

.94

Na2O

0.00 

.06 

.97 

.95

K2O

0.00 

.00 

1.000 

1.000

H2O +

0.01 

.01 

.96 

.93

reproduced data given in table 33B. Correspondence of the 

reproduced data to the original data given in table 33^4 is 

measured by the goodness-of-fit statistics given in table 39. 

The values of rj show that the model explains more than 90 

percent of the variance in all constituents except Fe 2O3 

and FeO.

TABLE 39. Goodness-of-fit statistics for the Kilauea factor model

[djj = Xjj-Xjj-, d? = mean of d,y; rfj* = standard deviation of dtf, TJ = correlation coefficient for 

% and *(y; r? = coefficient of determination]

LAYERED SERIES OF THE

SKAERGAARD INTRUSION,

GREENLAND

The Skaergaard intrusion in eastern Greenland has been 

the object of some classical studies in igneous petrology by 

L. R. Wager and his colleagues, and it is generally 

acknowledged as having attained its compositional 

variation almost entirely by fractional crystallization of a 

basaltic magma and by redistribution of crystals through 

settling and convection. The geology, petrology, and 

chemistry of the intrusion have been described in detail by 

Wager and Brown (1968), who also tabulated selected 

analytical data. The intrusion is exposed over an area of 50 

km 2 and is interpreted to have the form of an inverted 

cone. Most of it consists of a layered series that formed by 

gravity stratification of plagioclase, pyroxene, and olivine. 

The layered series has been subdivided into zones and 

subzones on the basis of cumulus minerals, those that 

precipitated and collected at the bottom of the magma. 

The three zones of the exposed part of the layered series are 

a lower zone (LZ) containing cumulus olivine, a middle 

zone (MZ) in which cumulus magnesian olivine is absent,
G40 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

and an upper zone (UZ) containing cumulus olivine that is 

richer in iron than that in the lower zone (Wager and 

Brown, 1968, p. 30). Wager and Brown estimated that 70 

percent of the layered series is hidden beneath the part that 

is exposed (1968, p. 204).

The factor-variance diagram derived from 19 analyses 

of specimens from the layered series (table 40) is given in 

figure 11. It shows that a petrologic model with five 

end-members could account for more than 85 percent of 

the variance in each compositional variable; each addition 

of an end-member beyond five will result in only small 

increases in the accountable variance for each constituent. 

However, an extensive search for five end-members failed 

to reveal a set that could be accepted as geologically 

plausible. Similarly, a search was made for sets of six to 

eight end-members. Again, no models that seemed 

-geologically acceptable could be found. The most common 

reason for the rejection of models containing five to eight 

factors was the appearance of negative composition 

loadings on the reference vector representing the 

composition of the parent magma, which would indicate 

that some samples formed by subtraction of the parent 

magma from crystallizing minerals. The composition of 

the parent magma was taken as that of specimen 4507 of 

Wager and Brown (1968, p. 150-151 and table 4) with the 

value for HsCH- set to 0.25 percent (Wager and Brown, 

1968, p. 192), and the composition was adjusted so that the 

nine constituents summed to 100 percent.

TABLE 40. Analyses, xjj, of samples from the layered series of the 

Skaergaard intrusion, adjusted so that nine major oxides sum to 100 

percent

[Analyses from Wager and Brown (1968, table 5 and fig. 15)1

Zone of layere 

series and 

sample Nos.

UZc:

5166

1881

4139

4142

1974

UZb:

4145

4272

UZa:

5181

1907

5322

5321

MZ:

3661

3662

LZa:

2308

2307

LZb:

4077

LZc:

5109

5108

4087

d

Si02

50.83

50.00

46.55

46.02

45.63

46.68

45.58

47.62

46.07

29.39

54.65

46.82

49.43

24.43

47.33

46.84

50.98

40.44

46.18

A1 203

9.61

8.89

9.65

8.22

8.04

12.39

11.20

15.29

14.35

1.71

22.44

15.47

18.50

5.63

13.75

16.99

25.71

3.36

16.67

Fe203

3.80

4.21

5.95

4.22

4.15

2.17

4.75

2.97

3.86

8.68

2.47

3.50

2.59

15.60

5.16

1.54

1.12

2.73

2.12

*

FeO

21.57

23.71

24.48

27.77

26.12

24.01

19.30

16.48

17.13

44.60

5.38

15.24

9.75

37.60

13.70

10.55

4.10

23.39

9.43

j

MgO

0.09

1.25

.44

.26

.39

1.81

4.15

4.35

5.70

10.88

.51

6.51

5.39

9.23

6.59

9.71

1.91

26.35

11.83

CaO

9.29

7.69

9.32

10.46

10.92

9.22

10.98

8.73

8.77

4.16

8.26

9.41

10.44

6.49

10.67

11.41

11.82

2.34

10.62

Na2O

3.19

2.74

2.50

2.24

2.88

3.12

3.59

3.86

3.44

.21

5.69

2.54

3.55

.61

2.34

2.48

3.56

.49

2.09

K2O

0.63

.35

.50

.49

.33

.37

.32

.28

.34

.06

.50

.29

.14

.07

.23

.20

.36

.08

.27

H20 +

0.99

1.17

.59

.31

1.55

.23

.12

.42

.35

.32

.09

.23

.21

.33

.22

.29

.44

.83

.78

1234567 

NUMBER OF FACTORS

FIGURE 11. Factor-variance diagram for the layered series of the 

Skaergaard intrusion.

The only geologically acceptable factor model that could 

be derived contained nine end-members, equal to the 

number of compositional variables used in the analysis. 

The disadvantage of a nine-factor model is that none of the 

variance in each constituent can be attributed to either 

analytical error, to compositional variation within each 

end-member, or to other geologic and petrologic processes 

not included in the model. Although the resultant model 

accounts for the observed data exactly, the three other 

general sources of variation are almost certain to have 

exerted at least a small amount of influence on the 

variation in the data. Consequently, the exactness with 

which the observed data are reproduced is, to some degree, 

a false indication of the correctness of the model. 

Nevertheless, the model can be expected to at least 

approximate a general process by which the Skaergaard 

intrusion may have formed.

The eigenvalues of the cosine theta matrix for the data in 

table 40 are given in table 41; the corresponding values of 

CPN when plotted against the number of factors display a 

smooth curve, having no prominent breaks that might 

indicate the proper number of end-members to include in 

the model.

The concept of the general process by which the layered 

series of the Skaergaard intrusion formed is that of a 

parent magma which was modified in composition during
p-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G41

TABLE 41. First nine eigenvalues of the cosine 

theta matrix for the layered series of the 

Skaergaard intrusion and CPN '

No.

1

2

3

4

5

6

7

8

9

Eigenvalue

15.240

1.935

.688

.613

.241

.131

.091

.038

.022

CPN

0.8021

.9040

.9402

.9724

.9851

.9920

.9968

.9988

1.0000

Total............18.999

'Cumulative proportions of N.

the process of solidification from the bottom of the magma 

chamber upward. Modification resulted from the 

crystallization of mineral phases and redistribution of 

these phases throughout the chamber by gravity settling 

and convection (Wager and Brown, 1968, p. 101). The 

problem of deriving a factor model, therefore, is one of 

finding mineral phases that (1) are conformable with the 

compositional series formed by the samples from the 

layered series and (2) can be represented by reference 

vectors on which the sample vectors have composition 

loadings that are subject to plausible interpretation. Above 

all, as stated previously, the loadings on the vector 

representing the composition of the parent magma must be 

all nonnegative. As pointed out in part I, however, when 

the number of factors in the model equals the number of 

compositional variables (m = M), any combination of the 

variables fits the model perfectly (all communalities are 

exactly one), and testing of mineral phases for conformity 

with the samples is impossible. Because of this, the mineral 

phases were tested within eight-factor space, even though 

the phases found to be satisfactory were then used in the 

nine-factor model. The validity of this procedure seemed 

apparent from the examination of the feldspar system.

The feldspar system of albite-anorthite-orthoclase was 

examined throughout at equal increments of 2 molar 

percent. That is, 1,326 compositions within this system 

were tested using a computer program based on the 

procedures given in parts I and II. Of these, only 54 

showed communalities equal to or greater than the average 

sample communality (0.9988) for the eight-factor solution 

(table 41). Nearly all of these 54 feldspar compositions 

included 2 molar percent orthoclase (fig. 12). This finding 

is in excellent accord with the data and the observations of 

Wager and Mitchell (1951, p. 146-147), who found, by 

analysis of three plagioclase specimens, 2 molar percent 

orthoclase in two specimens and 3 molar percent 

orthoclase in the third. They observed further that during 

fractionation of the magma the change in the K2O content 

of the plagioclase was small. Two of the end-members for

Anorthite

Albite Orthoclase

FIGURE 12. Molar compositions in the system albite-anorthiteorthoclase that have communalities in excess of 0.9988 in the eightfactor varimax solution for the layered series of the Skaergaard 

intrusion.

the nine-factor model are taken as Ab98Oro2 

An 98 Or 02 .

Cumulus olivines from the layered series vary in 

composition from Fa33 to almost pure fayalite (Wager and 

Brown, 1968, p. 28). Testing of the complete range of 

olivine compositions throughout the series, at increments 

of one molar percent, showed that all compositions have 

high communalities in the eight-factor space (fig. 13). 

Therefore, two additional end-members for the nine-factor 

model were taken as pure forsterite and pure fayalite.

The most abundant pyroxenes in the layered series are

1.000 r

0.999 -

0.998

0.997 -

0.996

Range of layered series olivines

10 20 30 40 50 60 70 80 90 100

Fa CONTENT OF OLIVINE, IN MOLAR PERCENT

FIGURE 13. Communalities of olivines in the eight-factor varimax 

solution for the layered series of the Skaergaard intrusion.
G42 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

augites that vary in composition between diopside 

(CaMg(SiO3 ) 2 ) and hedenbergite (CaFe(SiO3 ) 2 ). An 

Mg-Ca-Fe diagram given by Wager and Brown (1968, p. 

39), however, shows that these ideal compositions are not 

those of the actual end-members in the augite series. The 

compositional series formed by five augite specimens for 

which analyses are available (Wager and Brown, 1968, p. 

42) was examined by means of a special Q-mode analysis 

using the same procedures described throughout this 

report. The two-factor varimax solution is represented by 

the vector diagram given in figure 14 and shows that the 

augites can be well represented as a two-component 

system that is, all vectors representing pyroxene 

specimens are of near unit length. The theoretical extreme 

end-members in the system are represented by the two 

extreme vectors with all nonnegative composition 

scores vectors A and B in figure 14. Because vector A is 

close to vectors 13 and 15 (representing augite specimens 13 

and 15 of Wager and Brown), no particular advantage is 

gained by using the theoretical composition that vector A 

represents as an end-member composition in the factor 

model. Wager and Brown (1968, p. 39) regarded specimen 

15 as the extreme in the series and suggested that the 

composition of specimen 13 had been affected by processes 

of exsolution. Consequently, the compositions of the 

end-members for the augite series, used in the nine-factor 

model for the layered series, are taken as the composition

Parent

magma' . .

Forsterije^ . .

Fayalite^ . . .

Ca-poor

pyroxene' . Mg-augite4 . .

Fe-augite4 . . 

Iron ore'. . . .

Anorthite" . .

Albite6 .....

flcj (in percent)

SiO2

49.10

42.69

29.48

52.70

51.68

47.23 

0

43.62

68.65

A1 203

17.58

0

0

2.10

3.33

.94 

0

36.28

19.42

Fe203

1.35

0

0

1.00

1.52

0.60 

31.74

0

0

FeO

8.62

0

70.52

18.10

.48

31.83 

68.26

0

0

MgO

8.80

57.31

0

22.70

22.16

.14 

0

0

0

CaO

11.62

0

0

3.20

20.21

18.% 

0

19.76

0

Na2O

2.42

0

0

.20

.59

.26 

0

0

11.57

K20

0.26

0

0

0

.03

.03 

0

.34

.36

H2O +

0.25

0

0

0

0

0 

0

0

0

INITIAL VARIMAX LOADING, an

FIGURE 14. Q-mode vector diagram for a series of augite specimens 

from the layered series of the Skaergaard intrusion. Vectors A and B 

represent, respectively, the theoretical Fe- and Mg-augite extremes. 

Numbers on other vectors refer to specimen numbers of Wager and 

Brown (1968, p. 42).

of augite specimen 15 adjusted to sum to 100 and as the 

composition scores for vector B in figure 14. The 

end-members are referred to, respectively, as Fe-augite and 

Mg-augite, and their compositions are given in table 42.

Calcium-poor pyroxenes are far less abundant than 

augites in the layered series (Wager and Brown, 1968, table 

5) but are present there as well as in the border phases of 

the intrusion (Wager and Brown, 1968, p. 42). Analyses of 

two Ca-poor pyroxenes, a bronzite from the border phase 

and an inverted pigeonite from the middle zone of the 

layered series, are given by Wager and Brown (1968, p. 42), 

and the average, adjusted to sum to 100, is taken as the 

composition of another end-member for the nine-factor 

model (table 42).

TABLE 42. Compositions, fjq , of end-members for the Skaergaard 

factor model

1 Sample 4507 of Wagner and Brown (1968, table 7) with H2O+ set to 0.25 percent, adjusted 

sum to 100.

2 Theoretical.

3 Average of two analyses from Wager and Brown (1968, table 1), adjusted to sum to 100.

4 The composition of the Mg-augite is the theoretical extreme in the series of five analyses given 

by Wager and Brown (1968, table 1). The composition of the Fe-augite is taken as the analysis of 

augite specimen 15 from Wager and Brown (1968, table 1), adjusted to sum to 100.

' Theoretical mixture of iron ores. (See text.)

6 Anorthite is theoretical AnggOr^- Albite is theoretical AttygO^-

The reference vectors chosen thus far for the nine-factor 

model represent the parent magma, two extremes in the 

plagioclase series, two extremes in the olivine series, two 

extremes in the augite system, and a Ca-poor pyroxene. It 

seems that the remaining reference vector should represent 

the composition of the iron-ore minerals, principally 

ilmenite and magnetite, which constitute the only other 

cumulus mineral group of quantitative importance 

throughout the layered series. Although ilmenite is thought 

to have formed a cumulus mineral phase slightly before 

magnetite in the differentiation process, ilmenite and 

titaniferous magnetite tend to be found together (Wager 

and Brown, 1968, p. 49). It is appropriate, therefore, to 

represent these minerals by a single reference axis in the 

factor model. The composition of magnetite was taken as 

69 percent Fe2O3 and 31 percent FeO; because the data for 

TiO2 were not used in this exercise, the composition of 

ilmenite in terms of the constituents that were used is 100 

percent FeO. Neither of these compositions has a high 

communality in the eight-factor space for the layered series 

rocks, but, as shown in figure 15, a mixture of 46 percent 

magnetite and 54 percent FeO has a communality of
FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G43

1.000

0.950

0.900

46-percent Fe 3O4

0 10 20 30 40 50 60 70 

FeO, IN PERCENT

80 90 100

100 90 80 70 60 50 40 30 

Fe 3 O4 , IN PERCENT

20 10

FIGURE 15. Communalities of theoretical mixtures of ferrous iron oxide 

and magnetite in the eight-factor varimax solution for the layered 

series of the Skaergaard intrusion.

almost precisely 1 (h2 > 0.9999). Using the ideal compositions of magnetite and ilmenite, this mixture is that present 

in a cluster consisting of 29 percent magnetite and 71 

percent ilmenite, having a chemical composition of about 

32 percent Fe2O3 and 68 percent FeO. This chemical 

composition is accepted as that of the final end-member 

for the Skaergaard model (table 42).

Composition loadings for the 19 samples from the 

layered series with respect to the nine selected reference 

vectors are given in table 43. The loadings on the vector 

representing the parent magma range from about 0.4 to 

6.2, indicating that the model calls for some samples to 

have formed with a net addition of the constituents of the 

cumulus minerals to the parent magma (an < 1.0) and 

others to have formed with a net loss of these constituents 

(an > 1.0). The greatest losses of the constituents of the 

cumulus minerals, according to the model, were from the 

uppermost parts of the upper zone (UZc) and from the 

lowermost parts of the lower zone (LZa). As the cumulus 

minerals precipitated and either settled or were swept away 

by convection currents, they were replaced by influxes of 

magma from elsewhere in the magma chamber. Some of 

this additional magma, according to the model, must have 

been derived from outside of the volume occupied by the 

exposed layered series, because the exposed layered series 

as a whole requires about 1.454 parts parent magma (table 

43C). This value is the weighted average of the 

composition loadings on the parent material reference 

vector, using the volume weights for the zones and 

subzones of the layered series given by Chayes (1970, p. 

4). The most likely source of the additional parent material 

is the part of the magma chamber now occupied by the 

hidden part of the layered series.

The composition loadings of the 19 layered series 

samples on the vectors representing cumulus minerals 

(table 43) should not be interpreted as gains and losses due 

to precipitation of the end-members in the following 

two-component compositional series: forsterite-fayalite,

TABLE 43. Composition loadings, 

model, average loadings for each 

weighted averages for all zones

aik. for the 

zone of the

Skaergaard factor 

layered series, and

Zone of 

layered series1 Forsterite Ca-poor Fe-augite Anorthite 

and sample Parent pyroxene 

Nos magma Fayalite Mg-augite Iron ore

A. Composition loadings.

UZc:

5166

1881

4139

4142

1974

UZb:

4145

4272

UZa:

5181

1907

5322

5321

MZ:

3661

3662

LZc:

2308

2307

LZb:

4077

LZa:

5109

5108

4087

UZc ....

UZb. . . .

UZa....

MZ ....

LZ .....

3.975

4.686

2.344

1.257

6.204

.933

.469

1.669

1.400

1.280

.378

.907

.825

1.322

.886

1.171

1.772

3.301

3.123

3.693

.701

1.182

.866

1.929

3.171

-1.132

3.026

4.020

-2.814

1.992

1.741

-.385

.898

-.480

2.507

.683

-1.871

 .931

.172

-.878

-.653

-1.946

-1.457

1.255

1.867

.635

-.594

-.949

-3.922

.577

-3.630

-4.442

2.459

-2.035

-1.794

.404

 .946

.845

-2.787

-.771

1.957

1.017

-.400

1.043

.520

2.341

1.459

B.

 1.792

-1.915

-.621

.593

-.997

-1.029

.198

 .775

-1.048

.017

-.626

-.640

-.100

-.392

.084

-.624

-.096

.485

.119

.132

.088

.041

.466

.092

-8.776

.931

-7.995

-9.876

4.851

-4.830

 3.873

.637

-2.234

1.146

-6.009

-J.746

4.286

2.193

-.639

2.168

1.036

4.467

2.990

aik

8.484

-1.747

7.999

10.276

-5.840

5.043

4.332

-.803

2.212

-1.280

6.144

1.753

 4.396

-2.325

.714

-2.238

-1.402

-5.283

-3.567

0.244

-.085

.345

.392

-.256

.185

.254

.010

.140

.185

.253

.125

-.092

.370

.138

-.066

-.065

-.184

-.145

-0.864

-1.674

 .280

.201

-2.457

.122

.190

 .405

-.153

-.500

.548

.063

-.201

-.480

 .034

-.196

-.166

-1.398

-.950

Albite

-0.283

-.754

-.035

.219

-1.164

.218

.323

-.028

.075

-.281

.590

.081

.007

-.284

.031

-.092

-.084

 .764

 .545

Zone averages

-0.527

-.633

 .258

.195

.156

-4.173

-4.352

 1.615

1.270

2.036

3.834

4.687

1.568

-1.322

-2.350

0.128

.220

.147

.017

.008

-1.015

.156

 .127

 .069

 .537

-0.404

.270

J389

.044

 .290

C. Weighted averages for all zones 2

1.454 -0.266 0.248 0.014 0.425 -0.552 0.051 0.288 -0.085

1 From Wager and Brown (1968, p. 15 and table 5).

2 Weights used for UZc, UZb, UZa, MZ, and LZ were, respectively, 6.0, 11.4, 6.2, 37.1, and 

39.3 (Chayes, 1970, p 4).

Mg-augite - Fe-augite, and anorthite-albite. Rather, the 

loadings indicate the net effects of the precipitation of 

species of varying composition within each of these series. 

For example, the first sample (No. 5166) in table 43 formed 

from the parent magma partly by the precipitation of 

olivine of varying composition and partly by the 

concentration of the constituents of olivine in the residual 

magma; and so the net effect was the gain of 3.171 parts 

forsterite and the loss of 3.922 parts fayalite. This amounts 

to a net loss of 0.751 parts olivine, even though the 

constituents of olivine were both gained and lost during the 

course of differentiation. A detailed interpretation of the 

chemical and mineralogical changes in the parent magma 

that produced sample 5166 is given in table 44. It shows 

that the precipitation of olivine led to a net loss of FeO and 

net gains in MgO and Si(>2, even though there are no 

means of specifying the compositions of individual olivine 

crystals that were added to and lost from the parent 

magma. Similarly, the precipitation of compositions in the
G44 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

TABLE 44. Interpretation of chemical and mineralogical modifications of the parent magma required 

to produce sample 5166 from the Skaergaard intrusion

Minerals and

mineral SiO2 

groups

FeO MgO CaO Na2O K2O H2O Total

A. 3.9752 parts of parent magma

1.9519 0.6988 0.0537 0.3427 0.3498 0.4619 0.0962 0.0103 0.0099 3.9752

B. Gains aad losses by precipitation of minerals and mineral groaps

Olivine . .

Ca-poor

pyroxene . . .

Augite .

Iron

ore . . .

Plagioclase . .

Total . . .

.. 0

. .

 .

. 0

_

. -1.

.1975

,5420

5277

5716

4438

0

-.0216

-.2125

0

-.1686

-.6027

0

-.0103

-.0828

.0775

0

-.0156

-2.7659

-.1862

2.6586

.1666

0

-.1269

1.8174

-.2335

-1.9330

0

0

-.3491

0

-.0329

-.1652

0

-.1708

-.3689

0

-.0021

-.0295

0

-.0328

-.0644

0

0

-.0001

0

-.0040

-.0041

0

0

0

0

0

0

-0.7510

-1 .0286

-.2922

.2441

-1.1478

-2.9755

C. Net results as percentages

50.81 9.61 3.81 21.58 0.07 9.30 3.18 0.62 0.99 99.97

1 Compare with analyses of sample 5166 given in table 40. Subtract total in B from corresponding value in A and multiply by 100.

augite series caused a net loss of all chemical constituents 

except FeO. The net gain in FeO (table 44) may have 

resulted from the introduction of Fe-rich augite or its 

dissolved constituents into that part of the magma 

chamber at a late stage in the differentiation process, 

whereas augites that formed earlier and settled out were 

Mg-rich. The composition loadings for anorthite and albite 

show that the parent magma underwent a net loss of 

plagioclase in the formation of sample 5166, even though it 

is possible that some plagioclase, of undetermined 

composition and in either crystalline or dissolved form, 

was introduced into that part of the magma chamber 

sometime during the differentiation process. In fact, 

Wager and Brown (1968, table 5) reported the presence of 

cumulus plagioclase in this sample. The net result of the 

introduction and loss of crystalline or dissolved plagioclase 

from that part of the magma, however, is interpreted, on 

the basis of table 44, to have been a loss in all the 

constituents of plagioclase.

Interpretation of the positive composition loadings of 

table 43 as proportions of mineral phases in dissolved as 

well as crystalline form is important and necessary. For 

example, the composition loadings for sample 5166 on the 

iron-ore end-member is + 0.244, indicating that more than 

24 percent iron ore (comprised of 32 percent Fe2O3and 68 

percent FeO, table 42) has been added to this sample by the 

differentiation process. However, Wager and Brown 

(1968, table 5) reported only 6 percent iron ore in sample 

5166. If the iron ore had been added entirely in crystalline 

form, it is unlikely that it could have been resorbed into the 

magma, and it should be present in the rock. If added at 

least partly in dissolved form, however, it could provide 

the iron needed for the formation and precipitation of

olivine and pyroxene (tables 43, 44). Thus, some of the 

positive composition loadings on mineral end-members are 

interpreted to reflect residual concentration of the 

constituents of these minerals in the magma. Residual 

concentration of iron in the Skaergaard magma has been 

firmly established (Wager and Brown, 1968, p. 242).

Several features of the model presented in tables 42 and 

43 support its creditability, and others must be 

rationalized. Notable among the latter are that (1) the 

model does not account for the effects of other cumulus 

minerals present in minor amounts, such as apatite, and (2) 

it leaves none of the variance in each constituent to be 

attributed to the effects of analytical error, to minor 

petrologic processes, or to compositional variation within 

each of the end-member mineral phases. Regarding 

variation in the end-member mineral phases, the 

compositions of the Ca-poor pyroxenes in the Skaergaard 

intrusion, for example, are not perfectly represented by the 

end-member composition given in table 42, and the olivine 

compositions are not perfectly explained as mixtures of 

ideal forsterite and fayalite. At best, the model can be 

accepted as only an approximation of the effects of 

magmatic differentiation on each sample, in spite of the 

fact that it accounts for the observed data exactly.

On the other hand, the eight end-member compositions 

representing the average compositions of Ca-poor 

pyroxene, the iron-ore minerals, and the end-members in 

the olivine, augite, and plagioclase series are at least good 

approximations of the compositions certain to have been 

involved in the differentiation of the layered series. The 

least certain end-member composition is that of the parent 

magma, which was taken as interpreted by Wager and 

Brown (1968). The model, therefore, may be viewed as a
g-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G45

test of their interpretation, and the test shows no evidence 

that they were wrong.

If the model is accepted as approximately correct, the 

net gains and losses of the constituents of each mineral in 

the parent magma can be shown for each sample. In figure 

16 the gains and losses are plotted against the respective 

structural elevations of each sample in the layered series. 

The plots show that the constituents of olivine and 

Ca-poor pyroxene were largely lost from the parent magma 

in the upper parts of the magma chamber and that they 

were both lost from and added to the parent magma in the 

lower parts. Addition is interpreted to have been by crystal 

settling from upper parts of the chamber and by the 

transfer of crystalline material and magma with 

convection. Augite and plagioclase were mostly lost from 

the parent magma in the lower parts of the magma 

chamber, but they were both lost from and added to the 

upper parts. Inspection of the relative loadings on the 

Mg-augite and Fe-augite end-members (table 43) indicates 

that the added augites tended to be rich in iron. Similarly, 

the added plagioclases tended to be richer in sodium than 

in calcium. The addition of augite and plagiolcase, 

therefore, is interpreted to have been by concentration in

residual magma and transfer of the magma either by 

convection or by upward displacement due to settling 

crystals. Iron-ore minerals tended to be lost from the lower 

parts of the lower zone (LZa and LZb), which accords with 

Wager and Brown's (1968, p. 71) observation that 

magnetite and ilmenite are not present there as cumulus 

minerals. Above these zones, however, the constituents of 

these minerals were concentrated in the residual magma or 

were added as cumulus minerals perhaps both.

The concentration of iron-ore minerals in the residual 

magma in the upper parts of the magma chamber is 

significant with respect to the volume of the hidden part of 

the layered series. The requirement that the exposed part of 

the series formed from 1.454 parts parent magma implies 

that the hidden part forms at least 31 percent of the series 

as a whole. However if the model is correct with respect to 

iron-ore minerals, the part of the magma chamber now 

occupied by the exposed layered series received an addition 

of about 0.051 parts (table 43Q, or 5 percent, iron-ore 

minerals, either as crystalline magnetite and ilmenite or as 

dissolved phases. Judging from the communalities of the 

total range of mixtures (fig. 15), the iron-ore minerals 

would contain 46 percent magnetite, indicating the

UZcUZb

. UZaIUJ

^

z

HI

Hi

I

< 

DC

D

MZ,

LZcX

l

LZb

LZa

2500-

2000-

1500-

1000-

500-^

Olivine

-0.5 0 0.5

1 l

<

Ca-poor 

pyroxene

-0.5 0 0.5

1 i

-

Augite

-0.5 0 0.5

1 i

M

Iron ore

-0.5 0 0.5

l i

*

" 

Plagioclase

-3.0 -2.0 -1.0 0 1.0 2.0 3.0

I l l l l l l l l l l l i l

FIGURE 16. Interpretation of the net gains and losses of minerals in the parent magma plotted against structural height of sample in the layered 

series of the Skaergaard intrusion. The horizontal scale, representing gains and losses, is in units of proportions. The amounts of parent magma 

that contributed to the formation of each sample are given in table 43.
G46 STATISTICAL STUDIES IN FIELD GEOCHEMISTRY

addition of 2.3 percent magnetite to the exposed layered 

series, presumably from the part of the layered series that 

is hidden. If the original magma of the hidden layered 

series had the composition assumed for the parent magma 

(table 42), it contained a maximum of 2 percent magnetite 

and would have had to more than equal the exposed 

layered series in volume in order to have supplied the 

constituents of magnetite in the amounts required. If only 

one-half of the dissolved magnetite in the hidden layered 

series were transferred to the exposed part by residual 

concentration or by any other mechanism, the hidden 

layered series would necessarily have composed 70 percent 

of the total intrusion in order to supply the required 

amount. The estimate of 70 percent, however, is given only 

to illustrate how such estimates will vary with the assumed 

amounts of iron depletion in the hidden zone; the only 

absolute requirement of the model is that the hidden 

layered series forms at least one-half of the layered series as 

a whole. As pointed out before, Wager and Brown gave an 

estimate of 70 percent; Chayes (1970) suggested that 81 

percent is more correct.

SUMMARY

The four series of igneous rock analyses to which the 

extended method of Q-mode factor analysis has been 

applied vary a great deal in complexity. The series of 

analyses from the layered series of the Skaergaard 

intrusion is complex, not because the analyses for each 

constituent exhibit large variation among the samples 

(table 18) but because the compositional variation within 

the intrusion was caused by a complex interaction of 

many processes specifically, the precipitation of minerals 

of changing composition from a melt that also changed in 

composition as the processes continued. In contrast, 

although similar processes must have occurred to some 

extent in the formation of the rhyolite-basalt complex on 

the Gardiner River, the major part of the compositional 

variation in the complex was caused by only two 

processes the eruption of a rhyolite lava and the 

incorporation of basalt or basaltic lava. The 

factor-variance diagrams given in figures 3 and 11 clearly 

and immediately reflect this difference in origin. 

Moreover, the diagrams for the granitoid intrusive in 

eastern Nevada (fig. 5) and for the lavas and pumices of 

the 1959 summit eruption at Kilauea (fig. 8) reflect origins 

of intermediate complexity that are in at least fair accord 

with those arrived at after exhaustive field and laboratory 

study by other workers. After the samples and the data 

have been selected, the factor-variance diagrams are 

derived by a straightforward procedure with complete 

objectivity, and if the samples are representative and the 

analyses reasonably correct, they can be regarded as 

fundamental descriptors of the rock bodies or series from 

which the samples came.

The Q-mode factor analysis procedures that follow 

derivation of the factor-variance diagram may or may not 

contain some element of subjectivity. If the purpose of the 

factor analysis is only to condense geochemical or 

petrologic data into fewer variables without significant loss 

of information, then the choice of end-members is 

unimportant, and the varimax axes will serve as well as any 

other possible reference vectors. All that is required is to 

present the composition loadings for each sample with 

respect to each varimax axis in tables or on maps or to 

use these loadings to classify the samples into 

compositional groups. As a classification technique, the 

Q-mode method has an important advantage over common 

cluster methods and many petrographic diagrams that have 

been derived in that the classification scheme is not forced 

into a two-dimensional framework. Where the varimax 

axes are used as reference vectors for the compositional 

system, subjectivity enters the procedure only if and when 

an interpretation is made of what the varimax axes may 

represent. If no such interpretation is made, the varimax 

axes serve only as useful and conventional reference axes. 

Of course, the choice of the number of varimax reference 

axes could be partly subjective if, as can be expected to 

happen on some occasions, the factor-variance diagram 

does not clearly indicate how many there should be.

In the derivation of petrologic or geochemical models, it 

will usually be necessary to reject the principal components 

and varimax axes as representative of the compositions of 

the actual end-members in the compositional system. If the 

proper number of end-members in the system is clear from 

the factor-variance diagram, the search for more realistic 

end-members can be made either by the examination of the 

composition scores of selected vectors or by finding the 

vector representations of compositions of interest in the 

problem. If the proper number of end-members is not 

clear, the same procedures may have to be followed using 

factor spaces of several different dimensions. It will always 

be possible to derive a number of models that are 

mathematically satisfactory. The ultimate test of the model 

is its geologic plausibility, which must be inferred from the 

compositions of the selected end-members and the 

composition loadings of the samples with respect to them. 

The best reason for developing a factor model may be to 

subject a model developed by other means to rigorous 

quantitative examination.

REFERENCES CITED

Chayes, Felix, 1960, On correlation between variables of constant sum:

Jour. Geophys. Research, v. 65, no. 12, p. 4185-4193 

____1970, On estimating the magnitude of the hidden zone and the

compositions of the residual liquids of the Skaergaard layered series:

Jour. Petrology, v. 11, no. 1, p. 1 - 14. 

Chayes, Felix, and Kruskal, William, 1966, An approximate statistical

test for correlations between proportions: Jour. Geology, v. 74, no.

5, pt. 2, p. 692-702.
Q-MODE FACTOR ANALYSIS WITH CONSTANT ROW-SUMS G47

Christiansen, R. L., and Blank, H. R., Jr., 1972, Volcanic stratigraphy

of the Quaternary rhyolite plateau in Yellowstone National Park:

U.S. Geol. Survey Prof. Paper 729- B, 18 p. 

Fenner, C. N., 1938, Contact relations between rhyolite and basalt on

Gardiner River, Yellowstone Park: Geol. Soc. America Bull., V. 49,

no. 9, p. 1441-1484. 

 1944, Rhyolite-basalt complex on Gardiner River, Yellowstone

Park, Wyoming A discussion: Geol. Soc. America Bull., v. 55, no.

9, p. 1081-1096. 

Harman, H. H., 1967, Modern factor analysis [2d ed. rev.]: Chicago,

University of Chicago Press, 474 p. 

Imbrie, John, 1963, Factor and vector analysis programs for analyzing

geologic data: Office Naval Research, Geography Branch, Tech.

Rept. 6 [ONR Task No. 389-135], 83 p. 

Imbrie, John, and Purdy, E. G., 1962, Classification of modern

Bahamian carbonate sediments, in Classification of carbonate

rocks A symposium: Am Assoc. Petroleum Geologists Mem. 1,

p. 253 - 272. 

Imbrie, John, and VanAndel, T. H., 1964, Vector analysis of

heavy-mineral data: Geol. Soc. America Bull., v. 75, no. 11, p.

1131-1155. 

Klovan, J. E., and Imbrie, John, 1971, An algorithm and FORTRAN-IV

program for large-scale Q-mode factor analysis and calculation of

factor scores: Jour. Internal. Assoc. for Mathematical Geology, v. 3,

no. 1, p. 61-77. 

Klovan, J. E., and Miesch, A. T., 1975, Extended CABFAC and

QMODEL computer programs for Q-mode factor analysis of

compositional data: Computers & Geosciences, v. 1, no. 3 [1976]. 

Lee, D. E., and Van Loenen, R. E., 1971, Hybrid granitoid rocks of the

southern Snake Range, Nevada: U.S. Geol. Survey Prof. Paper 668,

48 p. 

Macdonald, G. A., 1968, Composition and origin of Hawaiian lavas, in

Coats, R. R., Hay, R. L., and Anderson, C. A., eds., Studies in

volcanology: Geol. Soc. America Mem. 116, p. 477-522. 

Manson, Vincent, and Imbrie, John, 1964, FORTRAN program for

factor and vector analysis of geologic data using an IBM 7090 or

7094/1401 computer system: Kansas Geol. Survey Spec. Distrib.

Pub. 13, 46 p. 

Miesch, A. T., 1969, The constant sum problem in geochemistry, in

Merriam, D. F., ed., Computer applications in the earth sciences 

An international symposium: New York, Plenum Press, p.

161-176. 

Miesch, A. T., 1975, Q-mode factor analysis of compositional data:

Computers & Geosciences, v. 1, no. 3 [1976]. 

Misch, Peter, and Hazzard, J. C., 1962, Stratigraphy and metamorphism

of Late Precambrian rocks in central northeastern Nevada and

adjacent Utah: Am. Assoc. Petroleum Geologists Bull., v. 46, no. 3,

p. 289-343. 

Murata, K. J., 1966, An acid fumarolic gas from Kilauea Iki, Hawaii:

U.S. Geol. Survey Prof. Paper 537-C, 6 p. 

Murata, K. J., and Richter, D. H., 1966, Chemistry of the lavas of the

1959 60 eruption of Kilauea Volcano, Hawaii: U.S. Geol. Survey

Prof. Paper 537-A, 26 p. 

Richter, D. H., Eaton, J. P., Murata, K. J., Ault, W. U., and Krivoy,

H. L., 1970, Chronological narrative of the 1959-60 eruption of

Kilauea Volcano, Hawaii: U.S. Geol. Survey Prof. Paper 537 - E, 73

PRichter, D. H., and Moore, J. G., 1966, Petrology of the Kilavea Iki lava

lake, Hawaii: U.S. Geol. Survey Prof. Paper 537-B, 26 p. 

Richter, D. H., and Murata, K. J., 1966, Petrography of the lavas of the

1959-60 eruption of Kilauea Volcano, Hawaii: U.S. Geol. Survey

Prof. Paper 537-D, 12 p. 

Shapiro, Leonard, and Brannock, W. W., 1956, Rapid analysis of silicate

rocks: U.S. Geol. Survey Bull. 1036-C, p. 19-56. 

Wager, L. R., and Brown, G. M., 1968, Layered igneous rocks:

Edinburgh and London, Oliver & Boyd, Ltd., 588 p. 

Wager, L. R., and Mitchell, R. L., 1951, The distribution of trace

elements during strong fractionation of basic magma A further

study of the Skaergaard intrusion, East Greenland; Geochim. et

Cosmochim Acta, v. 1, p. 129-208. 

White, R. W., 1966, Ultramafic inclusions in basaltic rocks from Hawaii:

Contr. Mineralogy and Petrology, v. 12, no. 3, p. 245 314. 

Wilcox, R. E., 1944, Rhyolite-basalt complex on Gardiner River,

Yellowstone Park, Wyoming: Geol. Soc. America Bull., v. 55, no. 9,

p. 1047- 1080. 

Wright, T. L., 1973, Magma mixing as illustrated by the 1959 eruption,

Kilauea Volcano, Hawaii: Geol. Soc. America Bull., v. 84, no. 3,

p. 849- 858.

* U.S. GOVERNMENT PRINTING OFFICE 1976-777-034/ 1 7