Mathematical Geology, Vol. 24, No. 6, 1992
Spatial R e l a t i o n s h i p s o f M u l t i v a r i a t e D a t a 1 E. C. Grunsky 2 and F. P. Agterberg 3
Spatial factor analysis (SFA) is a multivariate method that determines linear combinations of variables with maximum autocorrelation at a given lag. This is achieved by deriving estimates of auto-~ cross-correlations of the variables and calculating the corresponding eigenvectors of the covariance quotient matrix. A two-point spatial factor analysis model derives factors by the formation of transition matrix U comparing auto-~cross-correlations at lag "'0," R o, with those at a specified lag " d , " Rd, expressed as Ud = P~-IRd. The matrix Ua can be decomposed into its spectral components which represent the spatial factors. The technique has been extended to include three points of reference. Spatial factors can be derived from the relationship:
I
Udl~:[RoRd31
Ua2]
k Rd3Ro]
lIRdl ] L Ra~]
where the factors of Ual and Uaz predict the relationships of the variables over three lag distances and orientations from matrices Ra~, Ra2, and Rd3. Estimates of auto-~cross-correlation can be achieved using irregularly spaced or gridded data from experimental correlograms or function estimates determined from the experimental correlograms. The technique has applications in studies where several variables have different spatial ranges or zones of influence.
KEY WORDS: spatial factors, autocorrelation, crosscorrelation, correlogram, gridded data, irreg-
ularly spaced data. INTRODUCTION Multivariate relationships of data have been traditionally described by techniques such as principal components analysis which are based on linear combinations of variables derived from maximum variance. Such techniques do not consider the spatial relationships of the data. Multivariate relationships of spatially-based data have been studied by Agterberg (1966, 1974), Myers (1982, 1988), Royer 1988), Switzer and Green (1984), Wackernagel (1988), Grunsky and Agterberg Received August 4, 1991; accepted February 26, 1992. 2CSIRO, Perth, Western Australia, 6014, Australia; Present address: Geological Survey Branch, Ministry of Energy, Mines and Petroleum Resources, 553 Superior St., Victoria, B.C., Canada, V8V 1X4. 3Geological Survey of Canada, 601 Booth St., Ottawa, Canada, KIA OE8. 731 0882-8121/92/0800-0731$06.50/1 © 1992InternationalAssociationfor MathematicalGeology
732
Grunsky and Agterberg
(1988), Grunsky (1989), and Chung and Fabbri (in press). Gmnsky and Agterberg (1988) have shown that the application of a technique termed spatial factor analysis, produces linear combinations of variables that have maximum correlation based on auto-/cross-correlation estimates at specific lag intervals. Suppose the column vector of m variables, X, has lag 0 auto-/cross-correlations of R0, and lag d auto-/cross-correlations of R d where m
all
a12
...
aim m
a21
a22
...
a2m
am2
...
amm__
Ro =
_aml
--Fdl 1
Fdl 2 ...
Fdl m --
Fd21
Fd22
•..
Fd2 m
-Faml
Fdm2
...
Fdmm__
e d =
Fdo represents estimates of auto-/cross-correlations of variables Xs and Xj for lag d using a distance function F. Suppose that ~ is a column vector with m (unknown) elements. Then c~'X has lag 0 auto-/cross-correlations of c~'R0ce and lag d auto-/cross-correlations of c~'Rdc~. The problem is to maximize c~'Rdc~ subject to the constraint that c~'Roc~ --- 1. The solution satisfies: RdC~ = XRoc~
or
RolRdC~ = Xc~
where X is a Lagrange multiplier and c~ is an eigenvector of the transition matrix Ud such that Ua = Ro 1Rd" For m variables, Ua can be decomposed into m separate spectral components. Estimates of the overall noise of the system, and a measure of the goodness-of-fit (squared multiple correlation coefficient, R 2) to estimate the predictive power of each component can be computed. For the approach to be valid, Rd and Ro must be positive definite. Grunsky and Agterberg (1988) discuss some of the problems involved in obtaining positive definite estimates of R o and Rd. A correction procedure can be applied to modify Rd, or R 0 such that they are positive definite. These adjustments are usually small. A detailed explanation of the terms and the decomposition of Ud can be found in Grunsky and Agterberg (1988). The spectral components of Ud are equivalent to forming linear corn-
Spatial Relationships of Multivariate Data
733
binations based on maximum auto-/cross-correlation in the form of a quotient matrix as defined by Chung and Fabbri (in press). The preceding model is determined by the relationships between the variables Xoi (i = 1, 2 . . . . . m) and X u (j = 1, 2 . . . . . m) at two points Po and P~ that are a distance d apart. This two-point model can be written in the form:
Z6
= Zl
Ua + E ;
(1)
where Z 6 and Z I are row vectors consisting of standardized variables corresponding to Xoi and Xu; E D is a row vector consisting of residuals, If $1 represents a column vector of signal values corresponding to Z1, then premultiplication of both sides o f Eq. (1) by S l yields:
S1Z6 = S1Zl Ud + SIE6 This gives Ud = Rot Ra if S 1 and E 0 are uncorrelated or St E6 is a null matrix. Each column of Ua represents a set of regression coefficients by which the value of a variable at location Po is predicted from the values of all variables at location P1. If the residual variance for variable Z~ is written as 0 2, then the kth column Of UdhaS okR 2 o-1 as its variance-covariance matrix. The variances of the variables are proportional to the elements along the main diagonal of this matrix. Grunsky and Agterberg (1991a,b) have provided computer programs for the implementation of the spatial factor technique. The method can be extended to a three-point model whereby the simultaneous relationships between two points (PI, P2) and a third reference point (Po) can be predicted. This is illustrated in Fig. 1 which shows two possible configurations of three points in relation to each other. Of particular interest, are the relationships of different values of d and/or corresponding anisotropic functions (Fig. ta,b). This paper will only consider the case shown in Fig. la, with identical values of d, but with anisotropic relationships between the three points.
Po0)
d
i
P~(J/
/
P2(k) Fig. 1. Two possible graphical arrangements of auto-/crosscorrelation functions for the three locations; i, j and k represent the m variables at the three points. (a) anisotropic case: three points (P0, P1, P2), non-collinear, equidistant, (b) anisotropic case: three points (Po, P~, P2), non-collinear, non-equidistant.
Po(i)
P2(k)
d
PI(j)
734
Grunsky and Agterberg
Consider two locations, PI and P2, which are at distances d 1 and d2 from some reference location Po (see Fig. 1). In analogy with Eq. (1), the statistical model that expresses the relationships between Z0, the vector of variables at location Po, and Z I and Z2 representing two other vectors of variables (at locations P1 and P2 respectively) can be stated as,
Z~ = ZlUa, + Z~Ud2 + E~
(2)
The multiple regression model gives as a solution:
Ud2J
Rd3 Ro J
[-Rd2J
(3)
where Ral, Rd2, and Rd3 are matrices of auto-/cross-correlations at lags d 1, d2, 2~1/2, respectively. These three sets of auto-/cross-correlations andd3 = ( d ~ + d 2J are for the directions pointing from P0 to PI, Po to P2, and P1 to P2, respectively. The two- and three-point methods are similar to co-kriging in which a variable Xi at a point (Po, Fig. la) is related to the variables Xs or Xj and Xk at one or two other points (P1, P2; Fig. la) by estimating coefficients from the auto-/cross-correlations for the distances between the points. In the case of only one variable, this procedure is analogous to punctual kriging. Theoretically, the advantages of the three-point procedure are as follows. The two-point method is based on auto-/cross-correlations that represent either a multivariate average for all directions (isotropic function) or an average based on a search angle with a specific direction at its center. In general, it is more difficult to estimate the directional auto-/cross-correlation function because it is based on fewer points. If the third point in the three-point procedure is on the same line as the other two points, the method provides an extension of the twopoint procedure which is similar to the use of a second order Markov scheme in multiple time series analysis (Quenouille, 1957). The three-point method then can be applied in either the isotropic situation assuming that the same auto-/ cross-correlation function is valid in all directions, or in the anisotropic situation using single directions only. The third point does not have to be situated on the same line as the other two points. If there exists a clearly defined anisotropy in a region, the second and third point may be placed on lines pointing in the principal directions of this anisotropy (cf. Fig. la,b). A shorter distance between two points can be selected for the direction in which the auto-/cross-correlation function decreases rapidly (Fig. lb). The three-point model has twice as many coefficients as the two-point model. Half of these are for each of the two principal directions of the anisotropy. Theoretically, these two sets of coefficients can be treated as the single set of coefficients in the two-point model. The three-point model then is equiv-
Spatial Relationships of Multivariate Data
735
alent to separate two-point models for the two main directions of the anisotropy. The advantage of the three-point approach would be that the two sets of coefficients are estimated simultaneously, so that the effect of possible interaction between the sets of variables for the two directions is minimized. This procedure is equivalent to performing multiple regression on the combination of two sets of explanatory variables instead of analyzing them separately. Drawbacks of the three-point model are that it has to be based on directional auto-/cross-correlation functions which are difficult to estimate. Also, doubling the number of coefficients increases the possibility of the occurrence of multicollinearity and violations of the required condition of positive definiteness. M O D E L S OF A U T O - / C R O S S - C O R R E L A T I O N The spatial factor procedure requires that estimates of auto-/cross-correlation be obtained for the variables over the spatial range of the data. Auto-/cross-correlation relationships can be expressed such that a random variable Xi(i = 1 . . . . . m), with zero mean, at location k, in geographical space is related to Xj(j = 1 . . . . . m) at another location 1 by the function Fo(dkl ) with: Xik = Fij(dkl)Xjl
-t- Yik
(4)
Both k and l go from 1 to n where n denotes total number of observations (Agterberg, 1970). Fo.(dkl) is a function of distance dkl between the two points. This function is considered to decay to zero at the maximum distance of the "neighborhood" (range) for which it is defined. The residual is a random variable Y/k at point k. It satisfies E(Y/~) = 0, and Y/~ is assumed to be independent of X:-t. For simplicity, the subscripts k and I can be deleted. The choice of functions Fo(d) used to estimate auto-/cross-correlation is dependent on the nature of the spatial relationships. A number of models can be used such as the spherical, Gaussian, or exponential models. These functions have specific characteristics that are useful for describing the behavior of regionalized variables (Joumel and Huijbregts, 1978). A specific requirement of auto-/cross-correlation functions is that they satisfy the following conditions. For the autocorrelation of a random variable, Fii o r crosscorrelation of a pair of random variables, Fq, and lag, d: (1) ]Fij(d) I <_ 1, and (2) Fii(d) = F i i ( - d ) . As well, deterministic functions that obey the properties required of regionalized variables within defined limits can also be applied. Quadratic functions were used as estimates of spatial auto- and cross-correlation by Agterberg (1970) and Gmnsky and Agterberg (1988, 1989). Although quadratic functions are not positive definite, their use has shown them to be suitable estimates of the autocovariance function when used within the limits in which they exhibit positive definite characteristics (Haining, 1987).
736
Grunsky and Agterberg
In this study, three functions of auto-/cross-correlation were used to define the spatial variation of a multi-element geochemical dataset. For a pair of random variables, Xi, Xj, and distance d between any two locations, the functions are defined as: exponential model: Fij(d) = aij exp [ - d • bij], Gaussian model: Fij(d) = aij exp [ - ( d / b i j ) 2 ] , and quadratic function: Fij(d ) = aij + bijd + cijd 2. aij, bij, and cij are estimates of the model parameters, for the exponential, Gaussian, and quadratic models. The function estimates used here differ from the method used in a previous study (Grunsky and Agterberg, 1988) in which estimates were made using the method of ordinary least squares applied to irregularly spaced data. In this study, the model parameters were determined using the estimates of auto-/cross-correlation that were computed over specified lag intervals for gridded data. The function parameters were then calculated using the downhill simplex method of multidimensional minimization (Press et at., 1986, p. 289) whereby successive estimates of the function parameters are chosen based on the sum of squares between the observed autocorrelations and the function estimates. The parameter estimates converge as the sum of squares converges to a minimum. As the parameters are calculated, a sum of squares measure (SSQ) is computed such that for a given pair of variables (Xi, Xj), N
SSQ = ~1 [rij(d) - F~j(d)]2
(5)
i=
where N is the number of lags used to compute the function estimates, rij(d) is the observed auto-/cross-correlation estimate for Xi and Xj at lag d, and Fo.(d ) is the function estimate for the current set of parameters for Xi and Xj. The function with the smallest SSQ value is chosen as the function used to estimate the auto-/cross-correlations. In the examples used here, the choice of auto-/cross-correlation functions has been based on this minimum sum of squares criterion. Spatial data can be regularly or irregularly distributed over an area. When data are regularly distributed, better measures of auto-/cross-correlation can be obtained since the distance between points is exact and the lag intervals can be precisely defined. When the data points are irregularly distributed, the determination of auto-/cross-correlation at specific lag intervals is based on averaging the values of samples that fall within the lag interval. This procedure tends to smooth the estimates of the correlations and is thus not as precise as estimates derived from the irregularly spaced data. Gridding of the data also smooths the noise associated with each independent random variable with the result that the model of Eq. (4) is no longer valid. However, examination of the smoothed data is valuable for interpreting the general patterns of the data that are known to exist from empirical studies such as knowledge of the geology and the inferred geochemical relationships from it. As well, in cases where estimates of auto-/
Spatial Relationships of Multivariate Data
737
cross-correlation are difficult to obtain from irregularly spaced data, it may be necessary to regularize the data into a grid form. In this study, spatial factors have been calculated for both the gridded, and ungridded, irregularly spaced data, and the differences between them are examined. APPLICATION OF THE SPATIAL FACTOR TECHNIQUE TO THE BEN NEVIS DATA Archean lode gold deposits are typically associated with alteration in the form of carbonatization, sulphur and potassium enrichment, and proximity to shears and faults (Colvine et al., 1988). A suite of Archean age metavolcanics in the Ben Nevis area, Ontario, Canada contains three significant spatial patterns: a primary magmatic variation of the volcanics, a large regional zone of carbonatization centered along a major north-south fault, and smaller zones of K~O and S enrichment associated with Au-Cu mineralization. The geology of the area is shown in Fig. 2. Figures 3a-d display wire mesh plots of SiO~, CO~, S, and K~O, respectively. Topographically high regions reflect relative geochemical increases, and low regions represent relative geochemical decreases. The topographically high areas of Fig. 3a represent lithologies that are relatively enriched in SiO~ and outline the rhyolitic lithologies. Areas of carbonatized
mineral occurrences " ~ ~ maflc & Interrnedlate intrusive rocks fauit ~ [-2r~ maflc & intermediate volcanic rocks granitic rocks [ ' ~ - g ] ~ fslslc volcanic rocks
1 Canagsu Mine Deposit 2 Croxall Property
Fig. 2. Geology of the Ben Nevis Area. Note the two areas of mineralization, the Canagau Mines Deposit (1), and the Croxall Property (2).
738
Grunsky and Agterberg
a
\ ¢t,
\
b
Fig. 3. (a) SiO2 map. SiO 2 outlines the part of the lithological variation of the volcanics, distinguishing the SiO2 enriched felsic volcanics from the SiO2 depleted mafic volcanics and intrusions. (b) CO2 map. CO2 enriched areas occur in a north-south trend in the eastern part of the map area and in the western part of the map area around the Croxall property. (c) Sulphur map. Sulphur occurs along a north-south trend, partly enveloped by CO 2 in the eastern part of the map area. Several S-rich spikes occur throughout the area. A large S halo occurs around the Croxall property. (d) Potassium map. K20 shows a pattern associated with the felsic volcanics as well as alteration associated with the Croxall property and the Canagau Mines deposit.
Spatial Relationships of Multivariate Data
G
d \
Fig. 3. Continued
739
740
Grunsky and Agterberg
(C02) rock show up as topographically high features in Fig. 3b. Sulphur enrichment (S) associated with zones of pyrite and other sulphide mineralization shows up as peaks in an otherwise flat geochemical terrain in Fig. 3c, and Fig. 3d shows potassium enrichment (K20) representing lithological variation similar tO SiO 2 and alteration associated with CO2, (cf. Figs. 3a,b). Previous multivariate studies on the data are discussed in Grunsky (1986, 1988, 1989). Grunsky and Agterberg (1988) have shown that the application of the spatial factor technique to this geochemical dataset can be used to determine linear combinations of geochemical variables with similar spatial characteristics that represent primary magmatic trends, and zones of alteration potentially associated with mineralization. This approach can assist, along with other exploration methods, in the selection of potential exploration targets for ore mineralization. Many areas of mineralization often display a broad multi-element signature with significant addition or depletion of elements. This multi-element signature is often spatially larger than the deposit itself and is thus an important target in mineral exploration. The chemical patterns of the variables exhibit may reflect the activity of several geological processes. SiO2, A1203, FeO, MgO, CaO, Na20, and K20 lithogeochemical patterns generally reflect the compositional variation due to fractionation of magmas in the volcanic environment. It is the dominant geochemical pattern in the area. Oxides such as A1203, Fe203, CaO, Na20, and K20 usually reflect the primary compositional variation as well but their patterns also display the effects of other secondary processes (i.e., metamorphism, alteration). Elements such as CaO, Na20, and K20 reflect both the primary compositional variation of the volcanics as well as other secondary effects related to alteration. The process of hydrothermal alteration and mineralization has resulted in larger patterns that are manifested by CO2 which surround the mineralized S rich zones. CaO also reflects alteration that is associated with the mineralization near the Canagau Property. A zone of CaO enrichment occurs around the zone of carbonatization. Several sulphide occurrences exist throughout the area with minor amounts of base metal mineralization. The more significant mineralized sulphide occurrences have associated alteration within the surrounding host rocks. These occurrences are the primary targets for exploration and a multivariate approach might better assist in distinguishing these areas. SPATIAL DISTRIBUTION OF T H E BEN NEVIS DATA The geochemical dataset of the Ben Nevis area is composed of 825 samples that are irregularly spaced throughout the area (see Grunsky and Agterberg, 1988). For visual presentation, the data were gridded by dividing the area into 51 by 31 grid cells with 250-m spacing. The estimate of each variable on each
Spatial Relationships of Multivariate Data
741
grid node was based on the inverse squared distance weighting (Shepard, 1968) from the eight nearest neighbors at a maximum of four lag distances (1000 m). In this study, six elements were used (SiO2, FeO, CaO, K20, CO2, S) and auto-/cross-correlation estimates were obtained for these elements. Aitchison (1986) has discussed the problems of the interpretation of closed geochemical data (constant sum) and Pawlowsky (1989) has discussed the problem of closed data in the context of regionalized variables (also see article by Pawlowsky and Burger in this issue). This will not be dealt with here but will be the subject of future investigations. The following estimates of auto-/cross-correlation have been calculated for these elements: experimental correlograms and cross-correlograms derived from irregularly spaced data, experimental correlograms and cross-correlograms derived from gridded data, function estimates of the auto-/cross-correlations derived from irregularly spaced data, and function estimates of the auto-/crosscorrelation derived from data that were gridded using an inverse square smoothing technique (Shepard, 1968).
E S T I M A T E S OF E X P E R I M E N T A L C O R R E L O G R A M S Correlograms and cross-correlograms were computed for both the irregularly spaced and gridded data. The irregularly spaced data were binned into 120-m lag intervals with a tolerance (range of influence) of 60 m to a maximum lag of 6000 m. Binning of the data is a smoothing procedure, with a range of influence limited to the lag interval tolerance. Estimates of auto-/cross-correlation can be made for specific lag intervals and also within zones of specific orientations. Thus, correlogram and crosscorretogram estimation for 0 + 90 ° orientation results in an isotropic correlogram estimate. Anisotropic estimates that have been used in this study include the orientations, 0 _+ 15 ° (east-west), 90 _+ 15 ° (north-south), 45 ___ 15 ° (northeast-southwest), 135 ± 15 ° (southeast-northwest). Estimates of the correlations are poor for the smallest lag of 120 m since there are not enough samples available to estimate the auto-/cross-correlation. Subsequently, correlations at this lag were not considered. If a model for the regionalized variables is not assumed, then there can be no projected intercept at lag zero and hence the signal to noise ratio cannot be estimated. This also means that Ro can not be constructed. In this study, no autocorrelation model was considered and the intercept at lag zero was assigned the value of the ordinary correlation coefficient. Alternatively, visual estimates of the intercept at lag zero could have been used or the noise could be taken as the value derived from the auto-/crosscorrelogram estimate of the first lag.
742
Grunsky and Agterberg
E S T I M A T E S OF A U T O - / C R O S S - C O R R E L A T I O N F U N C T I O N S Estimates of function parameters were made for the exponential, Gaussian, and quadratic functions using both the irregularly spaced and gridded data. Grunsky and Agterberg (1988, 1989) used ordinary least squares for estimating auto-/cross-correlation quadratic function parameters based on a maximum lag limit or "neighborhood." Estimates of the function parameters were made over a specified range of neighborhoods (500 m, 1000 m, 1500 m, 2000 m) and orientations, (0 + 90 °, isotropic), (0 + 15 °, 90 + 15 °, 45 + 15 °, and 135 _+ 15 °, anisotropic). Function estimates for a pair of variables, Xi and Xj, are based on the average of the cross-correlation of Xi with Xs at distance d, and of Xj with X/at distance d. Two properties of these function estimates are that they estimate the auto-/cross-correlation most precisely in the midpoint range of the neighborhood and that the lag 0 estimate decreases as the neighborhood range increases. These properties are due to the averaging process of the least squares calculations. Estimates of the auto-/cross-correlation function parameters were computed differently for this study. Auto-/cross-correlation estimates were first obtained by creating correlograms and cross-correlograms for the data as described above, using 120-m lag intervals with a tolerance of 60 m to a maximum lag of 6000 m. As discussed above, the function parameters were then calculated using the downhill simplex method of multidimensional minimization (Press et al., 1986, p. 289) whereby successive estimates of the function parameters are chosen based on the sum of squares between the observed autocorrelations and the function estimates. For functions computed from the irregularly spaced data, for a given lag, d, the estimates of cross-correlation between variable Xi(0) and Xj(d) are not necessarily the same as the cross-correlation between Xj(0) and Xi(d). In this notation, (0) refers to a point Po and d to another point PI which is at distance d from Po. Thus, prior to computing the parameters of the functions, the crosscorrelations between two variables are averaged. The difference in correlations is due to the selection of samples based on distance and the binning of the samples into discrete lag intervals. The first auto-/cross-correlation lag estimate was generally poor and resulted in poor function estimates when included in the estimation of function parameters. Therefore, function estimates for the irregularly spaced data ignored the first lag estimate of the auto-/cross-correlation. For functions determined from the regularly spaced, gridded data auto-/ cross-correlation estimates were determined from a regularized set of points, spaced at 250 m, derived from the initial set of 825 irregularly distributed points as described previously. The function estimates were obtained from the lagged estimates of the auto-/cross-correlation in the same manner as for the irregularly spaced data. Similarly, for a given variable or pair of variables, the estimate of
Spatial Relationships of Multivariate Data
743
the auto-/cross-correlation for a particular lag was determined from the least squares minimum of the three functions. As an example, for SiOz, if the least squares minimum was determined from the exponential model, then that model was used to compute the autocorrelation estimate at lag 0 and lags dl, d2, and d3 for the construction of Ro, Rt, R2, and R3. In many cases, the sum of squares minima for the functions were very close in values and suggests that any one of the models describes the auto-/ cross-correlations equally well. In other cases however, the distinctions between functions are significant. Figure 4 shows the experimental correlograms and the plotted function estimates at orientations 0 +_ 15 ° and 90 _+ 15 ° for the ungridded Ben Nevis data for SiO2 and CO2. Figure 5 shows the experimental correlogram estimates and function estimates of the same variables for the gridded dataset. Auto-/ cross-correlogram estimates are displayed as small crosses joined together. Function estimates based on the three autocorrelation models, quadratic, exponential, and Gaussian are shown as squares, circles, and triangles, respectively. Figures 4 and 5 are clearly distinct from each other. First, the estimates of auto-/cross-correlation are less at lag 0 for the irregularly spaced data (Fig. 4), indicating that there is a component of noise in the unsmoothed regionalized variables. Figure 5 displays smoother and more continuous autocorrelation coefficients for the gridded data. In this case, the noise has been incorporated into the signal through the averaging process of the gridding procedure. The correlogram of Fig. 4a shows that SiO2 exhibits erratic features and indicates that the data displays a lack of stationarity. The function estimates suggest that the signal decays to zero at about 2000 m. A hole effect is apparent at a lag of 4000 m. Such features can add complications to the spatial factor analysis, however, for the purposes of this paper only lags less than 4000 m were considered. Because of the geology of the area, the corresponding geochemistry is anisotropic, and the function parameters vary with respect to orientation. Examination of experimental autocorrelation plots indicated that for most of the variables, the range of the variable was less than 4000 m. Figures 4a and b show the autocorrelation functions for SiO 2 at 0 ___ 15 ° and 90 + 15 °. Figure 4a shows considerable fluctuation of the autocorrelation over the range of 6000 m and indicates that SiO 2 varies considerably in the east-west direction. This implies a lack of stationarity in the data. The function estimates decay to zero at about 2000 m. Figure 4b shows a different pattern. The observed autocorrelation estimates are less erratic but the signal decays to zero at about 1200 m. The autocorrelation coefficient also increases significantly toward the origin. This is also reflected in the estimates for the exponential and Gaussian models. The difference in the rates of decay between the two often-
744 o
Grunsky and Agterberg Ben Nevis E x p e r i m e n t a l
Correlogram
Lag=120m
Si
SI
7"IQuadralic a= 0.0905 b= -0.0418 c= 0.0052
o 6
OExponenfial a= 0.1549 b= -1.1878
o
AGaussian a: 0.1103 b = 1.1525
"90. ~6
A n g l e : 0.0 Toler.= 15.0 Coy.= 54.76 R = 1.00
V
o
v - '~ - ' ~ S
"--
~ ' ' ' ' v -
+Observed Correlation
o
%oo
o;o
i
~8o
i
t
~,o
~.ao ,8o
3~o
Distance
o
9-
Ben
Nevis E x p e r i m e n t a l
~;o
(d)
Correlogram
Lag:120m
~;o
~io
Si
8.00
Si
mQuadrafic a = O. 1697 b = - 0 , 1 285 c= 0.0191
o i 00.
©Exponential a= 0.8616 b = -4.7,397
o
~Gaussian a= 0.4652 b= 0.5485
"9o.
o ol o o
_T_ _ ~ T p ~
~ . ~
~L~.
~
A n g l e = 90.0 Toler.= 15.0 Coy,= 54.76 R = 1.00
~ o
+Observed Correlation
o ~.00
0.~*0
1.80
2.40
3.20
4.00 Distance
4.so
5.so
6.40
7.20
8.00
(d)
Fig. 4. Function estimates and experimental correlograms of irregularly spaced data (lag = 120 m, Ben Nevis area): (a) SiO2 orientation 0 -I- 15 ° (anisotropic), (b) SiO2 orientation 90 + 15 ° (anisotropic), (c) COz orientation 0 _+ 15 ° (anisotropic), and (d) CO2 orientation 90 + 15 ° (anisotropic). Observed autocorrelations are shown as small crosses and the function estimates are denoted by symbols with the legend located on the right side of the diagram.
Spatial Relationships of Multivariate Data
3en Nevis E x p e r i m e n t a l C o r r e l o g r a m
745 Lag=120m
C02 C02
DQuadratic a= 0.4388 b= -0.2647 c= 0.0351
c~
©Exponential a= 0.5582 b= -1.1662 AGaussian a = 0.3931 b= 1.3094
~o k,- o o
Angle= 0.0 Taler. = 15.0 Cov.= 3.30 R = 1.00
,5 o
c~
+Observed
Correlation
o
,5 10.0o
0.80
1.so
2.4o
o Ben Nevis E x p e r i m e n t a l
&2o
4.oo
Distance (d)
Correlogram
4.8o
5.60
6.4.0
Lag=t20m
7.zo
8.00
C02 CO2
EX)uadrofic a= 0.2935 b= -0.0933 c= 0 . 0 0 8 2 ©Exponential a= 0.5224 b= -0.4468
o c~
AGaussion a= 0.2412 b= 2.9622 o
Angle= 90.0 Taler.= 15°0 Coy.= 3 . 3 0 R = 1.00 +Observed Correlation
0.00
0.80
1.60
2.40
3.20
4.1)0 4.80 Distance (d)
5.60
6.:,0
7.20
8.0D
Fig. 4. Continued
tations reflects the geochemical (and geological) anisotropy of the area. Both figures show low autocorrelation estimates at lag zero indicating considerable noise for this variable, possibly indicating a lack of stationarity. The hole effect is not apparent in the north-south direction.
746
Grunsky and Agterberg
~. Ben Nevis Dataset G L a g = 2 5 0 m
SI
Si
I~Quadratic a = 0.5535 b= - 0 . 1 8 6 9 c= 0.01 69
o (5
OExponenfial a= 0.7753 b= -0.6592
o (5
~Gaussian a= 0.5922 b= 1.7492
~o
o 6
A n g l e = 0.0 Toler.= 15.0 Coy.= 2 5 . 5 8 R = 1.00 +Observed Correlofion
o 1 IO.O0
o.8o
,.~o
2.~o
3.~o
,.bo
Disfance (d)
~.~o
5.6o
~.. Ben Nevis Dafoset G L o g = 2 5 0 m
o.~,o 7.~o
Si
8.00
Si
E]Quad ratic a= 0.4806 b= -0.2265 c= 0.0227
o= 6
OExponential a = 1.0001 b= -1.4522
o
~Goussion o = 0.8951 b= 0.8174
~o
o. (3
Angle= 90.0 Toler. = 15.0 Coy.= 2 5 . 5 8 R = 1.00
[!3 []
8.
+Observed Correlation °0.00
o8o
~6o
2~o
3~o
4~o
48o
56o
o.~,o 72o
8.00
Disfance (d)
Fig. 5. Function estimates and experimental correlograms of gridded data (lag = 250 m, Ben Nevis area): (a) SiO2 orientation 0 _+ 15 ° (anisotropic), (b) SiO2 orientation 90 _ 15 ° (anisotropic), (c) CO2 orientation 0 _+ 15 ° (anisotropic), and (d) CO2 orientation 90 _+ 15 ° (anisotropic). Observed autocorrelation are shown as small crosses and the function estimates are denoted by symbols with the legend located on the fight side of the diagram.
Spatial Relationships of Multivariate Data o ,2
747
Ben Nevis Datoset G L a g = 2 5 0 m
C02
C02
[]Quadratic a = 0.5454 b= - 0 . 2 0 5 4 c = 0.0166
i o
©Exponential a= 1.0000 b= - 1 . 0 9 6 9 AGaussian a : 0.7811 b= 1.2112
~o
°i
Angle= 0.0 Taler. = 15.0 Cov.= 1.05 R = 1.00
[] o
o
[!] i
,5
o o 10.00
+Observed Correlation
I
0;0
,;o
2;0
~i0
,;0 Distance
(d)
,;o
5.80
oo Ben Nevis Dotoset G L a g = 2 5 0 m
~0
~io
8.00
C02 C02
E]Ouodretic o = 0.6360 b= - 0 . 1 4 5 6 c= 0.0105
o o
OExponenfiai a = 0.6572 b= - 0 . 2 6 6 3
o o
AGaussian a= 0.4764 b = 5.1317
u_o
o
A n g l e = 90.0 T a l e r . = 15.0 Coy.= 1.05 R = 1.00
c5
o
o
+Observed Correlation
o c5
I0.0,
0.80
1.60
2.40
3.zo
4.~)0
4.BO
5.60
8.40
7.ZO
8.00
Distance (d) d
Fig. 5. Continued
Figures 4c and d show the autocorrelation functions for C O 2 at 0 + 15 ° and 90 + 15 °, respectively. Figure 4c shows a distinctive Gaussian curve and suggests that the autocorrelation o f CO2 is more continuous at the shorter lag intervals. The function decays to zero at about 2400 m. Figure 4d also displays a curve with Gaussian characteristics. However, at this orientation (90 +_ 15°),
748
Grunsky and Agterberg
the curve decays to zero at about 5000 m and reflects the continuous northsouth zone of carbonatization that occurs in the southern part of the map area (cf. Fig. 3b). Figures 5a-d show the estimates of autocorrelation for the gridded data of both SiO2 and CO2. As expected, the curves are smoother with fewer fluctuations which is directly the result of the averaging process of gridding. Figures 5a and b show similar characteristics to Figs. 4a and b. For SiO2 at 0 __+ 15 °, the function decays slowly to about 6000 m. The hole effect at 4000 m is also present but the increase of the autocorrelation coefficient is relatively less than the surrounding estimates in comparison to Fig. 4a. Estimates of the autocorrelation at lag 0 are also considerably higher than for the irregularly spaced data. Figure 5b shows the estimates for the orientation 90 _ 15 °. As in Fig. 4b, the autocorrelation decays rapidly to zero at about 1600 m. Estimates of the autocorrelation at lag 0 are also significantly greater than in Fig. 4b. Figures 5c and d show the autocorrelation estimates for the gridded data for CO2. Figure 5c shows the estimates for the orientation 0 + 15 ° which decays to zero at about 2000 m and is similar to Fig. 4c. Figure 5d shows the estimates for the orientation 90 + 15 ° in which the autocorrelation decays very slowly and maintains a significant value past the range considered here ( > 6000 m) which confirms the continuity of the main zone of carbonatization in the north-south direction. The lagged estimates of the irregularly spaced data display patterns that follow the observed geological features. However, these estimates have considerable noise at lag zero and can cause problems in applying the spatial factor method. As an alternative, the gridded data can be used for the estimates of the auto-/cross-correlation, although the noise can no longer be assumed to be independent and becomes part of the signal through the gridding process. T H R E E - P O I N T S P A T I A L F A C T O R ANALYSIS A three-point spatial factor analysis was applied to the Ben Nevis data for three lags and two orientations using the configuration of points in Fig. la. The analysis was applied to both the gridded and irregularly spaced data for both the experimental and function estimates of auto-/cross-correlation. Because estimates of auto-/cross-correlation cannot be estimated for lag 0, in the case were experimental correlograms were used, R o was constructed using two methods. The first method consists of using the ordinary correlation coefficient. This implies that there is no noise at lag zero as the diagonal elements are equal to 1.0. The second method uses the auto-/cross-correlation coefficient estimate taken from the first lag of the correlogram. In this case, the amount of noise is overestimated. As well, at smaller lags, the estimate of the auto-/cross-correlation is less reliable because fewer points were used in the calculation. Eval-
Spatial Relationships of Multivariate Data
749
uation of spatial factors using both methods can assist in seeing how the noise component influences the linear relationships of the variables. EXPERIMENTAL
I R R E G U L A R L Y SPACED 750 m, R 2 : 9 0 + 15 ° , d2 = 750 m,
CORRELOGRAMS:
DATA (RI: 0 + 15 ° , d l = R3:45
+ 15 °, d3 = 1060 m)
The results of the analysis are shown in Table 1 where the R0, R~, R2, R3, UI, and U2 matrices are listed along with the eigenvalues, X. U1 represents the spatial associations of the three points for the orientation 0 _+ I5 ° and corresponding lag distance dl (points Po-Pl of Fig. la). U2 represents the spatial associations of the three points for the orientation 90 + 15 ° and corresponding lag distance d2 (points Po-P2of Fig. la). In this example, the lag distances for the two pairs of points are the same. Table 1 shows that the solution for U~ has five positive and one negative eigenvalues. All auto-/cross-correlation matrices R o, R~, R~, R3, were checked to ensure positive definiteness. Small adjustments were made to RI, R 2, and R 3. However, R03 and R~2 are not necessarily positive definite and there is no simple computational procedure for ensuring such a condition. Because U~ has only five eigenvalues greater than 0, the results cannot be interpreted in a meaningful way. Similarly, there are only four positive eigenvalues for U2, and the corresponding eigenvectors and R 2 coefficients cannot be interpreted. E X P E R I M E N T A L C O R R E L O G R A M S : G R I D D E D DATA (RI: 0 + 15 ° , d l = 750 m, R 2 : 9 0 + 15 ° , d2 = 750 m, R 3 : 4 5 + 15 ° , d3 = 1060 m)
The results of the analysis from correlograms derived from gridded data did not have negative eigenvalues. Some of the results are shown in Table 2. In this example, R o was derived from the ordinary correlation coeffÉcients. R 0 could have also been derived using the auto-/cross-correlation estimates of the first lag. None of the matrices required adjustments for positive definite solutions. Examination of the components for UI (points Po-PI of Fig. 1) show that the first three components are the most significant. The value Q (see Gmnsky and Agterberg, 1988) determines the relative significance of the U matrix and the components, and shows the order of significance to be components 1, 3, and 2, respectively. The sum of the values of Q for the components is less than or equal to the value of Q for the'U-matrix. The row of R 2 coefficients of the first component is dominated by CaO, S i Q , and CO2. Figure 6a shows a wire mesh map of the first component of U~. The positive relief of the map reflects zones of relative CaO enrichment that surround the main zone of carbonatization (c.f. Fig. 3b). Examination of the amplitude vector T in Table 2 shows that the
750
Grunsky and Agterberg
I
I I I
I
II
I
I
II
<
I
I
II
If
I
I l l
I
I
I
I
6
I I I
II
II
I I I I I
I
i
+l
li
I I I I
II
I I I
f
II
I
I
. ~ ~ II
I .o r~
~ ~ : I
~1
I
I
I I I
~. I
I
~.~ I
.-.~.
I
I
I
8 +i ~..~ I
~: ~
~. ~. I
I
. . . . . . I
~
I
I
II
II
Z "~
I"
I'
"
I"
"
'~
'
I'
I. . . .
i'1"i''1'1""
~
"
I
Spatial Relationships of Multivariate Data
751
Table 2. Three-Point Spatial Factor Analysis"
UI .3546 .0440 - .0798 - .0875 .0162 - .0129 Ut eigenvalues .4123 Component I 2 3 Component
-.1308 .1665 .0548 .0809 -.0047 .0136
-.1403 -.0403 .3027 .1045 -.0282 -.0165
.1231 .0867 .0099 .1802 -.0308 .0157
.0734 .0472 -.0352 -.0804 .3176 .0060
-.0923 -.0609 -.0280 .0198 -.0193 .3326
.3359
.3104
.2330
.2119
.1509
Amplitude eigenvectors T -.6324 .4865 .6593 .1481 -.1682 .0834 .5165 -.4559 -.3375
-.0529 -.0281 .2704
-.4344 -.1718 -.6349
.2207 -.9019 -.3248
Squared multiple cowelation coefficients (R2) SiO2 FeO CaO K20 .0550 .0264 .0779 .0013 .0602 .0320 .0740 .0000 -.0051 -.0056 .0039 .0013 .0268 .0235 .0145 .0115
CO2 .0366 .0430 -.0063 .0527
Q .0494 .0363 .0131 .0197
U 1 2 3
S .0990 .0087 .0903 -.0105
U2 .0759 -.0066 -.0428 .1607 -.0265 .0050 U2eigenvalues .4318 Components 1 2 Component U 1 2
.2231 .3017 .0729 -.1405 .0669 -.0192
.0342 -.0127 .2131 -.1254 -.0121 .0231
.0557 .0161 .0280 .3584 .0037 .0543
.1733 .2045 .0129 .0281 .3660 -.0414
.0786 .0272 .0893 .0821 -.0362 .1746
.3852
.2590
.2097
.1222
.0817
Amplitude eigenvectors T -.2436 .7538 -.0907 .6357 -.5752 -.5611
-.2238 1.0537
t.8069 .0787
Q .0455 .0208 .0248
SiO2 .0253 .0001 .0252
-.3519 .3223
Squared multiple co,elation coefficients (R2) FeO CaO KzO CO 2 .0303 .0199 .0745 .1081 .0115 .0008 .0000 .1073 .0188 .0191 .0745 .0008
S .0151 .0049 .0102
aBen Nevis correlogram estimates based on gridded data (SiO2, FeO, CaO, K20, CO2, S). Underlined coefficients indicate the significant variables.
752
Grunsky and Agterberg
first component has positive score contributions from CaO and FeO, which also contribute to the positive scores in the map pattern of Fig. 6a. The first two components of U2 (points Po-P2 of Fig. 1) are the most significant. The second component is more significant than the first and is dominated by K20. The first component is nearly as significant and is dominated by CO2. Figure 6b displays a wire mesh map of the first component. The positive relief of the map closely resembles that of Fig. 3b which outlines the CO2 abundance of the area. Examination of the amplitude eigenvectors of Table 2 shows that CaO and CO2 contribute to the positive scores of the first component. Figure 6c shows a wire mesh map of the third component which is dominated by K20 as shown by the R 2 coefficients. The amplitude eigenvector for the third component indicates a strong positive loading for K20. The pattern closely resembles the map of K20 in Fig. 3d.
FUNCTION ESTIMATES OF A U T O - / C R O S S - C O R R E L A T I O N F R O M IRREGULARLY SPACED DATA (RI: 0 + 15 °, dl = 750 m, R 2 : 9 0 +_ 15 ° , d2 = 750 m, R 3 : 4 5 + 15 ° , d3 = 1060 m) The results of the spatial factor analysis applied to function estimates based on irregularly spaced data are shown in Table 3. The first eigenvalue for U1 is very large (58.03) and indicates that the result is unrealistic and may represent a condition of multicollinearity. Additionally, the last eigenvalue is less than zero and this violates the requirement of real positive eigenvalues for a second order Markov scheme (Quenouille, 1957) and the results may not be interpretable. For U2 there are two complex eigenvalues which also violates the requirement of a solution of positive real eigenvalues for a meaningful solution. The failure to derive linear combinations of variables using these functions possibly suggests that they are poor estimates of the auto-/cross-correlations of the variables. As stated above, the directional estimates of auto-/cross-correlation are not as good as isotropic estimates since they are based on fewer points. Combined with an increased number of coefficients in the solution of Uj and U2, the possibility of multicollinearity can have an adverse effect on obtaining a successful solution.
FUNCTION ESTIMATES OF A U T O - / C R O S S - C O R R E L A T I O N F R O M GRIDDED DATA (RI: 0 _+ 15 ° , dl = 750 m, R 2 : 9 0 _+ 15 ° , d2 = 750 m, R 3 : 4 5 +_ 15 ° , d3 = 1060 m) The results of the spatial factor analysis applied to the functions generated from the gridded data are shown in Table 4. The solution for U1 (points Po-Pl of Fig. 1) required slight adjustments for R o and the eigenvalues are all positive. The second component is the most significant and is dominated by a positive
Fig. 6. (a) Map of the first component of U~ for lag = 750 m based on correlogram estimates from gridded data. The pattern reflects zones of CaO and FeO enrichment. (b) Map of the second component of U2 for lag = 750 m based on correlogram estimates from gridded data. The pattern outlines zones of K20 enrichment. (c) Map of the first component of U2 for lag = 750 m based on correlogram estimates from gridded data. The pattern outlines zones of CO2 enrichment. (d) Map of the second component of U~ for lag = 750 m based on function estimates from gridded data. The pattern outlines the zones of S, CaO, and CO2 enrichment.
754
Grunsky and Agterberg
Fig. 6. Continued
Table 3. Three-Point Spatial Factor Analysis~
U1 -1.1724 -.6987 -.9421 .4988 .0022 -2.8183 U1 eigenvalues 58.0320
2.1993 1.4374 1.1816 -1.2762 -.0409 5.5794
-2.2334 - 1.6180 -.2891 1.9557 .1480 -7.2780
-7.6519 -3.7651 -4.9336 3.9563 .1785 -20.0012
-2.2984 - 1.0944 -1.5254 1.0670 .5940 -5.6886
20.6729 10.6874 11.8493 -11.6892 -.7503 55.3459
1.0097
.5121
.2878
.0768
-.0464
1.1474 .5201 .8223 -.4541 .1146 2.8477
-2.1650 -.8564 - 1.6330 .9577 .0752 -5.6267
3.0267 1.5974 1.9649 -1.6562 -.2817 7.2278
7.1162 3.3960 4.5396 -3.5806 -.0795 20.0557
2.1190 1.1317 1.1102 -1.1539 .31t2 5.4701
-20.9262 -10.7316 - 12.3832 11.5133 .7413 -54.7290
.0476
-57.3462
Complex roots (r, i) (.56202397E + 00, .82824874E - 01) Complex roots (r, i) (.56202397E + 00, -.82824874E - 01) U2 eigenvalues .5620
.5620
.3020
.1301
"Ben Nevis function estimates based on irregularly spaced data (SiO2, FeO, CaO, K20, C O 2, S),
Spatial Relationships of Multivariate Data
755
Table 4. Three-Point Spatial Factor Analysisa
Ut eigenvNues
1.1266
.7023
Component 1 2 3 4
-.4433 -.5607 .4207 -.5749
.2590 .4517 -.3200 .4614
.4956
.4377
Amplitude eigenvec~ T -.8389 .0704 .8931 -.2608 -.7349 -.0166 .1958 -.4858
.2223
.2085
-.9932 .8607 .3045 .8689 .5996 .1514 .2758 -.3471
Noise component (1 - t/n)*lO0 = 17.0632 Component U 1 2 3 4 U2 eigenvalues
Q .3084 .0577 .0665 .0374 .0351
SiO2 .2802 .0100 .0563 .0344 .0614
.5958
.4645
Squamdmultiple correlation coefficients (R2) FeO CaO K20 CO 2 .1646 .3918 .1481 .5128 -.0020 .0489 .0016 .2075 .0220 .1276 .0032 .0374 .0165 .0901 -.0004 .1035 .0345 .0152 .0435 .0177 .3144
.2627
S ,5003 .1483 .1983 .0081 .0303
.1693 -.1992
~Ben Nevis function estimates based on gridded data (SiO2, FeO, CaO, K20), CO2, S). Underlined coefficients indicate the significant variables.
association between S and CaO. The amplitude eigenvector o f U 1 also indicates a positive loading contribution from CO2. This is graphically shown in Fig. 6d which shows a wire mesh diagram reflecting a combination of S (Fig. 3c), CO2 (Fig. 3b), and CaO (not shown) patterns. The first component is dominated by CO2 and S. The third component is characterized by an inverse relationship between CO2 and CaO. The fourth component shows a positive association between SiO2 and K20. The solution for U2 (points Po-P2 of Fig. 1) contains one negative eigenvalue and, subsequently, the linear combinations o f the variables for U2 are not meaningful.
CONCLUDING
REMARKS
The spatial factor technique can determine linear combinations o f variables that have m a x i m u m autocorrelation for a given lag interval. Thus, the technique can be used to explore the relationships o f regionalized variables over various spatial domains with specific orientations and lag intervals. Further to the initial application o f quadratic functions o f auto-/cross-correlation estimates outlined in G m n s k y and Agterberg (1988), we have demonstrated that other measures o f auto-/cross-correlation can be incorporated into the spatial factor technique. As in traditional principal components analysis many multi-element relationships
756
Grunsky and Agterberg
(linear combinations of variables) can be explained through known models of geological processes whereas other linear combinations cannot be explained and may be the result of random processes. The advantage of using the estimates of auto-/cross-correlation from the irregularly spaced data is that an estimate of noise can be attained and the measures of autocorrelation follow the functions of the autoregressive models. The spatial factor technique could not be applied effectively to auto-/crosscorrelations derived from correlograms since the lag 0 estimate could not be obtained. If no noise is assumed, then the ordinary correlation coefficient can be substituted for the auto-/cross-correlation coefficients at lag 0. The use of gridded data has the advantage that the estimates of auto-/crosscorrelation are smoother and the results of the two and three points spatial factor methods yield more stable results, particularly where lack of stationarity is a problem. The main disadvantage of the use of the gridded data is that the associated noise of each data point can no longer be estimated and thus the noise becomes correlated. Alternatively, the auto-/cross-correlation estimates of the first lag could be used, but the noise component would be greater than it actually is. The use of experimental estimates of auto-/cross-correlation from the correlograms provides a measure of comparison for testing the "reality" of the results compared with function estimates. In many cases, the use of the modeled function estimates for evaluating U1 and U2 failed to provide meaningful results. This may be due to poor estimates of the functions from the ungridded, irregularly spaced data, poor auto-/ cross-correlation estimates of the correlograms, a lack of stationarity of the data or multicollinearity of the system of equations, or spurious spatial correlations due to the fact that compositional data are dealt with. Further research is needed to see if actual results are validated when considering these aspects. Despite these difficulties, the method can be a useful adjunct at studying spatially-based phenomena in which the multivariable data can be described in linear combinations based on maximum auto-/cross-correlation.
ACKNOWLEDGMENTS Assistance from the Division of Exploration Geoscience of the Commonwealth Scientific Industrial Organization (CSIRO) of Australia and the Geological Survey Branch, British Columbia Ministry of Energy, Mines and Petroleum Resources is gratefully acknowledged. This paper is Geological Survey of Canada contribution No. 50591. Many helpful suggestions and comments were given by Harri Kiiveri, Division of Mathematics and Statistics, CSIRO, Australia, Vera Pawlowsky, Department of Applied Mathematics, University of
Spatial Relationships of Multivariate Data
757
C a t a l u n y a , B a r c e l o n a , Spain, and C h a n g - J o C h u n g , G e o l o g i c a l Survey o f C a n ada, O t t a w a .
REFERENCES Agterberg, F. P., 1966, The Use of Multivariate Markov Schemes in Petrology: J. Geol., v. 79, p. 764-785. Agterberg, F. P., 1970, Autocorretation Functions in Geology, in D. F. Memam (Ed.), Geostatistics: Plenum, New York, p. 113-142. Agterberg, F. P., 1974, Geomathematics: Elsevier, Amsterdam, 596 p. Aitchison, J., 1986, The Statistical Analysis of Compositional Data: Chapman and Hail, New York, 416 p. Chung, C. F., and Fabbri, A. G., in press, Two Approaches to Parsimony of Multi-channel Remote Sensed Data, in G. Hill (Ed.), Advanced Data Integration in Mineral and Energy Resources Studies: USGS Bulletin. Colvine, A. C., et al., 1988, Archean Lode Gold Deposits in Ontario: Ontario Geological Survey, Miscellaneous Paper 139, 136 p. Grunsky, E. C., 1986, Recognition of Alteration in Volcanic Rocks Using Statistical Analysis of Lithogeochemical Data: J. Geochem. Expl., v. 25, p. 157-183. Grunsky, E. C., 1988, Multivariate and Spatial Analysis of Lithogeochemical Data from Metavoicanics with Zones of Alteration and Mineralization in Ben Nevis Township, Ontario: Unpublished Ph.D. thesis, University of Ottawa. Grunsky, E. C., t989, Spatial Factor Analysis: A Technique to Assess the Spatial Relationships of Multivariate Data, in F. P. Agterberg and G. F. Bonham-Carter (Eds.), Statistical Applications in the Earth Sciences: Geological Survey of Canada Paper 89-9, p. 329-347. Gmnsky, E. C., and Agterberg, F. P., 1988, Spatial and Multivariate Analysis of Geochemical Data from Metavolcanic Rocks in the Ben Nevis Area, Ontario: Math. Geol., v. 20, n. 7, p. 825-861. Gmnsky, E. C., and Agterberg, F. P., 1989, The Application of Spatial Factor Analysis to Unconditional Simulations with Implications for Mineral Exploration, in Applications of Computers and Operations Research in the Mineral Industry (APCOM), Chap. 20:1989 Proceedings, p. 194-208. Grunsky, E. C., and Agterberg, F. P., 1991a, FUNCORR: A FORTRAN 77 Program for Computing Multivariate Spatial Autocorrelation: Comput. Geosci., v. 17, n. 1, p. 115-131. Gmnsky, E. C., and Agterberg, F. P., 1991b, SPFAC: A Fortran 77 Program for Spatial Multivariate Analysis: Comput. Geosci., v. 17, n. 1, p. 133-160. Haining, R., 1987, Trend Surface Models with Regional and Local Scales of Variation with an Application to Aerial Survey Data: Technometrics, v. 29, n. 4, p. 461-469. Journel, A. G., and Huijbregts, C. J., 1978, Mining Geostatistics: Academic Press, London, 6 p. LeMaitre, R. W., 1982, Numerical Petrology, Statistical Interpretation of Geochemical Data: Elsevier, New York, 281 p. Myers, D. E., 1982, Matrix Formulation of Co-kriging: Math. Geol., v. 14, p. 249-257. Myers, D. E., 1988, Some New Aspects of Multivariate Analysis, in A. G. Fabbri, C. F. Chung, and R. Sinding-Larsen (Eds.), Statistical Treatments for Estimation of Mineral and Energy Resources; Proc. NATO Conference Held at I1 Ciocco, Italy, July 1986, Reide!, Dordrecht, p. 669-687. Pawlowsky, V., 1989, Co-kriging of Regionalized Compositions: Math. Geol., v. 21, p. 512-521. Press, W. H., Flannery, B. P., Teukolsky, S. A., and Vetterling, W. T., 1986, Numerical Recipes, The Art of Scientific Computing: Cambridge University Press, 818 p.
758
Grunsky and Agterberg
Quenouille, M. H., 1957, The Analysis of Multiple Time-Series: Hafner, New York, 105 p. Royer, J.-J., 1988, Geochemical Data Analysis, in A. G. Fabbri, C. F. Chung, and R. SindingLarsen (Eds.), Statistical Treatments for Estimation of Mineral and Energy Resources: Proc. NATO Conference Held at I1 Ciocco, Italy, July 1986, Reidel, Dordrecht, p. 89-112. Shepard, D., 1968, A Two-Dimensional Interpolation Function for Irregularly Spaced Data: Proceedings of the 23rd National Conference, Association for Computing Machinery, Brandon, Princeton, N.J., p. 517-524. Switzer, P., and Green, A. A., 1984, Min/Max Autocorrelation Factors for Multivariate Spatial Imagery: Tech. Report No. 6, April 1984, Department of Statistics, Standard University, 14 p. Wackernagel, H., 1988, Geostatistical Techniques for Interpreting Multivariate Spatial Information, in A. G. Fabbri, C. F. Chung, and R. Sinding-Larsen (eds.), Statistical Treatment for Estimation of Mineral and Energy Resources: Proc. NATO Conference Held at I1 Ciocco, Italy, July 1986, Reidel, Dordrecht, p. 393-349.