Brereton Chemometrics

mean centred spectra is 10.313, hence the first two components account for 100 ×

(7.944 + 1.178)/10.313 = 88.4 % of the overall variance, so the root mean square

√

error after two PLS components have been calculated is 1.191/(27 × 25) = 0.042

(since 1.191 is the residual error) or, expressed as a percentage of the mean centred data,

√

E% = 0.042/ 10.313/(27 × 25) = 40.0 %. This could be expressed as a percentage of the mean of the raw data = 0.042/0.430 = 9.76 %. The latter appears much lower and is a consequence of the fact that the mean of the data is considerably higher than the standard deviation of the mean centred data. It is probably best simply to determine the percentage residual sum of square error (= 100 − 88.4 = 11.6 %) as more components are computed, but it is important to be aware that there are several approaches for the determination of errors.

The error in concentration predictions for pyrene using two PLS components can be computed from Table 5.13:

•the sum of squares of the errors is 0.385;

•dividing this by 22 and taking the square root leads to a root mean square error of 0.128 mg l−1;

Relative to the standard deviation of the centred data it is even higher. Hence the ‘x’ and ‘c’ blocks are modelled in different ways and it is important to recognise that the percentage error of prediction in concentration may diverge considerably from the percentage error of prediction of the spectra. It is sometimes possible to reconstruct spectral blocks fairly well but still not predict concentrations very effectively. It is best practice to look at errors in both blocks simultaneously to gain an understanding of the quality of predictions.

The root mean square errors for modelling both blocks of data as successive numbers of PLS components are calculated for pyrene are illustrated in Figure 5.13, and those for acenaphthene in Figure 5.14. Several observations can be made. First, the shape of the graph of residuals for the two blocks is often very different, see especially acenaphthene. Second, the graph of c residuals tends to change much more dramatically than that for x residuals, according to compound, as might be expected. Third, tests for numbers of significant PLS components might give different answers according to which block is used for the test.

Figure 5.13

Root mean square errors in x and c blocks, PLS1 centred and pyrene

The errors using 10 PLS components are summarised in Table 5.15, and are better than PCR in this case. It is important, however, not to get too excited about the improved quality of predictions. The c or concentration variables may in themselves contain errors, and what has been shown is that PLS forces the solution to model the apparent c block better, but it does not necessarily imply that the other methods are worse at discovering the truth. If, however, we have a lot of confidence in the experimental procedure for determining c (e.g. weighing, dilution, etc.), PLS will result in a more faithful reconstruction.

5.5.2 PLS2

Figure 5.14

Root mean square errors in x and c blocks, PLS1 centred and acenaphthene

a concentration matrix, C, rather than concentration vectors for each individual compound in a mixture, and the algorithm is iterative. The equations above alter slightly in that Q is now a matrix not a vector, so that

X = T .P + E

C = T .Q + F

The number of columns in C and Q are equal to the number of compounds of interest. In PLS1 one compound is modelled at a time, whereas in PLS2 all known compounds can be included in the model simultaneously. This is illustrated in Figure 5.15.

Table 5.15 Concentration estimates of the PAHs using PLS1 and 10 components (centred).

Figure 5.15

Principles of PLS2

Table 5.16 Concentration estimates of the PAHs using PLS2 and 10 components (centred).

It is a simple extension to predict all the concentrations simultaneously, the PLS2 predictions, together with root mean square errors being given in Table 5.16. Note that there is now only one set of scores and loadings for the x (spectroscopic) dataset, and one set of ga common to all 10 compounds. However, the concentration estimates are different when using PLS2 compared with PLS1. In this way PLS differs from PCR where it does not matter if each variable is modelled separately or all together. The reasons are rather complex but relate to the fact that for PCR the principal components are calculated independently of the c variables, whereas the PLS components are also influenced by both blocks of variables.

In some cases PLS2 is helpful, especially since it is easier to perform computationally if there are several c variables compared with PLS1. Instead of obtaining 10 independent models, one for each PAH, in this example, we can analyse all the data in one go. However, in many situations PLS2 concentration estimates are, in fact, worse than PLS1 estimates, so a good strategy might be to perform PLS2 as a first step, which could provide further information such as which wavelengths are significant and which concentrations can be determined with a high degree of confidence, and then perform PLS1 individually for the most appropriate compounds.

5.5.3 Multiway PLS

Two-way data such as HPLC–DAD, LC–MS and LC–NMR are increasingly common in chemistry, especially with the growth of coupled chromatography. Conventionally either a univariate parameter (e.g. a peak area at a given wavelength) (methods in Section 5.2) or a chromatographic elution profile at a single wavelength (methods in Sections 5.3 to 5.5.2) is used for calibration, allowing the use of normal regression techniques described above. However, additional information has been recorded for each sample, often involving both an elution profile and a spectrum. A series of two-way chromatograms are available, and can be organised into a three-way array, often visualised as a box, sometimes denoted by X where the line underneath the array name indicates a third dimension. Each level of the box consists of a single chromatogram. Sometimes these three-way arrays are called ‘tensors’, but tensors often have special properties in physics which are unnecessarily complex and confusing to the chemometrician. We will use the notation of tensors only where it helps in understanding the existing methods.

Enhancements of the standard methods for multivariate calibration are required. Although it is possible to use methods such as three-way MLR, most chemometricians have concentrated on developing approaches based on PLS, to which we will be restricted below. Theoreticians have extended these methods to cases where there are several dimensions in both the ‘x’ and ‘c’ blocks, but the most complex practical case is where there are three dimensions in the ‘x’ block, as happens for a series of coupled chromatograms or in fluorescence excitation–emission spectroscopy, for example. A simple simulated numerical example is presented in Table 5.17, in which the x block consists of four two-way chromatograms, each of dimensions 5 × 6. There are three components in the mixture, the c block consisting of a 4 × 3 matrix. We will restrict the discussion for the case where each column of c is to be estimated independently (analogous to PLS1) rather than all in one go. Note that although PLS is by far the most popular approach for multiway calibration, it is possible to envisage methods analogous to MLR or PCR, but they are rarely used.

5.5.3.1 Unfolding

One of the simplest methods is to create a single, long, data matrix from the original three-way tensor. In the case of Table 5.17, we have four samples, which could be arranged as a 4 × 5 × 6 tensor (or ‘box’). The three dimensions will be denoted I , J and K. It is possible to change the shape so that any binary combination of variables is converted to a new variable, for example, the intensity of the variable at J = 2 and K = 3, and the data can now be represented by 5 × 6 = 30 variables and is the unfolded form of the original data matrix. This operation is illustrated in Figure 5.16.

It is now a simple task to perform PLS (or indeed any other multivariate approach), as discussed above. The 30 variables are centred and the predictions of the concentrations performed when increasing number of components are used (note that three is the maximum permitted for column centred data in this case, so this example is somewhat simple). All the methods described above can be applied.

An important aspect of three-way calibration involves scaling, which can be rather complex. The are four fundamental ways in which the data can be treated:

1.no centring;

2.centre the columns in each J × K plane and then unfold with no further centring, so, for example, x1,1,1 becomes 390–(390 + 635 + 300 + 65 + 835)/5;

and 4 provide radically different answers; for example, sample 2 has the highest value (=170.75) using method 3, but the lowest using method 4 (= −87.85).

Standardisation is also sometimes employed, but must be done before unfolding for meaningful results; an example might be in the GC–MS of a series of samples, each mass being of different absolute intensity. A sensible strategy might be as follows:

1.standardise each mass in each individual chromatogram, to provide I standardised matrices of dimensions J × K;

2.unfold;

3.centre each of the variables.

Standardising at the wrong stage of the analysis can result in meaningless data so it is always essential to think carefully of the physical (and numerical) consequences of any preprocessing which is far more complex and has far more options than for simple two-way data.

After this preprocessing, all the normal multivariate calibration methods can be employed.

5.5.3.2 Trilinear PLS1

Some of the most interesting theoretical developments in chemometrics over the past few years have been in so-called ‘multiway’ or ‘multimode’ data analysis. Many such methods have been available for some years, especially in the area of psychometrics, and a few do have relevance to chemistry. It is important, though, not to get too carried away with the excitement of these novel theoretical approaches. We will restrict the discussion here to trilinear PLS1, involving a three-way x block and a single c variable. If there are several known calibrants, the simplest approach is to perform trilinear PLS1 individually on each variable.

Since centring can be fairly complex for three-way data, and there is no inherent reason to do this, for simplicity it is assumed that data are not centred, so raw concentrations and chromatographic/spectroscopic measurements are employed. The data in Table 5.17 can be considered to be arranged in the form of a cube, with three dimensions, I for the number of samples and J and K for the measurements.

Trilinear PLS1 attempts to model both the ‘x’ and ‘c’ blocks simultaneously. Here we will illustrate the use with the algorithm of Appendix A.2.4, based on methods proposed by de Jong and Bro.

Superficially, trilinear PLS1 has many of the same objectives as normal PLS1, and the method as applied to the x block is often represented diagrammatically as in Figure 5.17, replacing ‘squares’ or matrices by ‘boxes’ or tensors, and replacing, where necessary, the dot product (‘.’) by something called a tensor product (‘ ’). The ‘c’

Figure 5.17

Representation of trilinear PLS1

block decomposition can be represented as per PLS1 and is omitted from the diagram for brevity. In fact, as we shall see, this is an oversimplification, and is not an entirely accurate description of the method.

In trilinear PLS1, for each component it is possible to determine

• a scores vector (t), of length I or 4 in this example;

• a weight vector, which has analogy to a loadings vector (j p) of length J or 5 in this example, referring to one of the dimensions (e.g. time), whose sum of squares equals 1;

•another weight vector, which has analogy to a loadings vector (k p) of length K or 6 in this example, referring to the other one of the dimensions (e.g. wavelength) whose sum of squares also equals 1.

Superficially these vectors are related to scores and loadings in normal PLS, but in practice they are completely different, a key reason being that these vectors are not orthogonal in trilinear PLS1 influencing the additivity of successive components. Here, we keep the notation scores and loadings, simply for the purpose of retaining familiarity with terminology usually used in two-way data analysis.

In addition, a vector q is determined after each new component, by

q=(T .T )−1.T .c

so that

cˆ = T .q

where T is the scores matrix, the columns of which consist of the individual scores vectors for each component and has dimensions I × A or 4 × 3 in this example if three PLS components are computed, and q is a column vector of dimensions A × 1 or 3 × 1 in our example.

A key difference from bilinear PLS1 as described in Section 5.5.1 is that the elements of q have to be recalculated afresh as new components are computed, whereas for two-way PLS, the first element of q, for example, is the same no matter how many components are calculated. This limitation is a consequence of nonorthogonality of individual columns of matrix T.

The x block residuals after each component are often computed conventionally by

where resid,a xij k is the residual after a components are calculated, which would lead to a model

Sometimes these equations are written as tensor products, but there are numerous ways of multiplying tensors together, so this notation can be confusing and it is often conceptually more convenient to deal directly with vectors and matrices, just as in Section 5.5.3.1 by unfolding the data. This procedure can be called matricisation.

In mathematical terms, we can state that

A

where unfolded pa is simply a row vector of length J.K. Where trilinear PLS1 differs from unfolded PLS described in Section 5.5.3.1 is that a matrix Pa of dimensions J × K can be obtained for each PLS component given by

Pa =jpa .kpa

and Pa is unfolded to give unfolded pa .

Figure 5.18 represents this procedure, avoiding tensor multiplication, using conventional matrices and vectors together with unfolding. A key problem with the common implementation of trilinear PLS1 is that, since the scores and loadings of successive components are not orthogonal, the methods for determining residual errors are simply an approximation. Hence the x block residual is not modelled very well, and the error matrices (or tensors) do not have an easily understood physical meaning. It also means that there are no obvious analogies to eigenvalues. This means that it is not easy to determine the size of the components or the modelling power using the x scores and loadings, but, nevertheless, the main aim is to predict the concentration (or c block),