Howard Mark

Mark Electronics

21 Terrace Avenue

Sufffern, NY 10901


Prepared for the FDA Process Analytical Advisory Committee

February 20, 2002


This report describes a new method of testing analytical data for linearity. This method overcomes the shortcomings of the current and proposed recommendations for linearity testing.


Let us begin by discussing what we want to test. The FDA/ICH guidelines, starting from a univariate perspective, considers the relationship between the actual analyte concentration and what generically called the “test result”. This terminology therefore holds good for every analytical methodology from manual wet chemistry to the latest high-tech instrument. In the end, even the latest instrumental methods have to produce a number, representing the final answer for that instrument’s quantitative assessment of the concentration, and that is the test result from that instrument. This is a univariate concept to be sure, but the same concept that applies to all other analytical methods. This is currently the way all analytical results are reported and evaluated. So the question to be answered, for any given method of analysis, is: ”Is the relationship between the instrument readings (test results) and the actual concentration linear?”


The FDA/ICH guidelines provide variant descriptions of the meaning of the term linearity. One definition is: “... ability (within a given range) to obtain test results which are directly proportional to the concentration (amount) of analyte in the sample.”(1) This is an extremely strict definition, one which is unattainable in practice when noise and error are taken into account. Figure 1 illustrates the problem. Here we have a plot of a set of hypothetical data that most of us would agree represents a linear relationship between the test result and the analyte concentration. While there is a line that meets the criterion that “test results are directly proportional to the concentration of analyte in the sample”, none of the data points fall on that line, therefore in the strict sense of the phrase, none of the data representing the test results can be said to be proportional to the analyte concentration. In the face of non-linearity of response, there are systematic departures form the line as well as random departures, but in neither case is any data point strictly proportional to the concentration.


Less strict descriptions of linearity are also provided. One recommendation is visual examination of a plot (unspecified, but presumably also of the method response versus the analyte concentration). This method works fairly well, but is subjective and open to different interpretations. It also is not amenable to application of computerized automated screening methods.


Another recommendation is to use “statistical methods”: calculation of a linear regression line is advised. If regression is performed, the correlation coefficient, slope, y-intercept and residual sum of squares are to be reported. These requirements are all in keeping with their background of being applied to univariate methods of analysis. There is no indication given, however, as to how these quantities are to be related to linearity, and as Anscombe shows(2) they are not. Figure 3 presents two of the data sets from the Anscombe paper. One is a data set which is substantially linear, the other is a data set which is obviously very non-linear. The key point about those two sets of data is that when linear regression is performed, as recommended by the guidelines, all the regression statistics are identical for the two sets of data. Therefore it is immediately obvious that the regression results cannot distinguish between the two cases, since the regression results are the same for both of them.


In fact, the recommendations in the official guidelines for assessing linearity, while well-intended, are themselves not suitable for their intended purpose in this regard, not even for univariate methods of analysis. Therefore, let us propose a definition, that can at least serve as a basis for our own discussions. Let us define linearity as the characteristic of data such that a straight line provides as good a fit (using the least-squares criterion) as any other mathematical function, as a description of the relationship between the method response and the concentration of the analyte.


As a test for linearity, the Durbin -Watson statistic(3-5) has been proposed(6, 7) and is a step in the right direction, but it also has several shortcomings, including the fact that it can be fooled by data that had just the right (or wrong!) characteristics. For a relationship between test results and analyte concentration that is linear but contains random, independent, Normally distributed errors, the expected value of DW is 2. But for the data sequence:


            0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, ...


the value of the Durbin-Watson statistic for this set is also two, even though the sequence is neither random nor Normally distributed but is definitely non-linear.


DW also requires a relatively large number of samples, much larger than is needed for the linearity test to be described here. The method we now present is mathematically sound, more subject to statistical validity testing, based on well-known mathematical principals, is sensitive to smaller numbers of samples and can distinguish different types of non-linearity from each other.


This new method of determining non-linearity can be viewed from a number of different perspectives, and can be considered as coming from several sources. One way to view it is as having a pedigree as a method of numerical analysis(8).


The new method of determining non-linearity (or showing linearity) can also be related to a discussion of derivatives, particularly when using the Savitzky-Golay method of convolution functions. This is not very surprising, once you consider that the Savitzky-Golay convolution functions are also (ultimately) derived from similar considerations of numerical analysis.


In some ways it also bears a resemblance to the current method of assessing linearity that the FDA and ICH guidelines recommend, that of fitting a straight line to the data, and assessing the goodness of the fit. However, based on the work of Anscombe(2), we see that the currently recommended method for assessing linearity is faulty because it cannot distinguish linear from non-linear data, nor can it distinguish between non-linearity and other types of defects in the data.


But an extension of that method can.


Having defined “linearity” as the characteristic of data such that a straight line provides as good a fit (using the least-squares criterion) as any other mathematical function, as a description of the relationship between the method response and the concentration of the analyte, we see that this seems to be almost the same as the FDA/ICH approach, which we just discredited. But there is a difference. The difference is the question of fitting other possible functions to the data; the FDA/ICH guidelines only specify trying to fit a straight line to the data. This is also more in line with out own proposed definition of linearity. We can try to fit functions other than a straight line to the data, and if we cannot obtain an improved fit, we can conclude that the data is linear


It is possible to fit other functions to a set of data, using least-squares mathematics. In fact, this is what the familiar Savitzky-Golay method does. The Savitzky‑Golay algorithm, however, does a whole bunch of things, and lumps all those things together in a single set of convolution coefficients: it includes smoothing, differentiation, curve-fitting of polynomials of various degrees, least-squares calculations and finally combines all those operations into a single set of numbers that you can multiply your measured data by to directly get the desired final answer.


For our purposes, though, we don’t want to lump all those operations together. Rather, we want to separate them and retain only those operations that are useful for our own purposes. For starters, we discard the smoothing, derivatives and performing a successive (running) fit over different portions of the data set, and keep only the curve-fitting. Texts dealing with numerical analysis tell us what to do and how to do it. Many texts exist dealing with this subject, but we will follow the presentation of Arden(8). Arden points out and discusses in detail many applications of numerical analysis. These methods are all based on using a Taylor series to form an approximation to a function describing a set of data. The nature of the data, and the nature of the approximation considered differs from what we are used to thinking about, however. The data is assumed to be univariate (which is why this is of interest to us here) and to follow the form of some mathematical function, although we may not know what the function is. Using a Taylor series implies that the approximating function that we wind up with will be a polynomial, and perhaps one of very high degree (the “degree” of a polynomial being the highest power to which the variable is raised in that polynomial.)


The concepts of interest to us are contained in Arden’s book in a chapter entitled “Approximation”. This chapter takes a slightly different tack than the rest of the discussion in the book, but one that goes in exactly the direction that we want to go. In this chapter, the scenario described above is changed very slightly. There is still the assumption that there is a single (univariate) mathematical system (corresponding to our desired “analyte concentration” and “test reading”), and that there is a functional relationship between the two variables of interest although again, the nature of the relationship may be unknown. The difference, however, is the recognition that data may have error, therefore we no longer impose the condition that the function we arrive at must pass through every data point. We replace that criterion with a different criterion, and the criterion we use is one that will allow us to say that the function we use to describe the data “follows” the data in some sense. The most common criterion used for this purpose is the “least squares” principle: to find parameters for any given mathematical function that minimizes the sum of the squares of the differences between the data and a corresponding point of the function.


Similarly, many different types of functions can be used. Arden discusses, for example, the use of Chebyshev polynomials, which are based on trigonometric functions (sines and cosines). There are other types of polynomials that could also be used, such as Legendre polynomials, Jacobi polynomials, and others. But by far the simplest to deal with, and therefore the most widely used approximating functions are simple polynomials; they are also convenient in that they are the direct result of applying Taylor’s theorem, since Taylor’s theorem produces a description of a polynomial that estimates the function being reproduced:


            Y = a0 + a1X + a2X2  + a3X3 + ...  + anXn     (Equation 1)


For some data a polynomial can provide a better fit to that data than can a straight line. We present an example of  that result as figure 2, for ease of reference. Higher order polynomials may provide an even better fit, if the data requires it.


The mathematics of fitting a polynomial by least squares are relatively straightforward, and we present a derivation here, one that follows Arden, but is rather generic, as we shall see: Starting from equation 1, we want to find coefficients (the ai) that minimize the sum-squared difference between the data and the function’s estimate of that data, given a set of values of X. Therefore we first form the differences:


            D = a0 + a1X + a2X2 + a3X3 + ... + anXn  - Y    (Equation 2)


Then we square those differences and sum those squares over all the sets of data (corresponding to the samples used to generate the data):


            åiD2 = åi (a0 + a1X + a2X2 + a3X3 + ... + anXn  - Y)2    (Equation 3)


The problem now is to find a set of values for the ai that minimizes SD2 with respect to each ai. We do this by using the usual procedure of taking the derivative of SD2 with respect to each ai and setting each of those derivatives equal to zero. Note that since there are n + 1 different ai. we wind up with n + 1 equations, although we only show the first three of the set:


(åiD2 ) / a0 = (å (a0 + a1X + a2X2 + a3X3 + ... + anXn  - Y)2 ) / a0 = 0  (Equation 4a)


(åiD2 ) / a1 = (å (a0 + a1X + a2X2 + a3X3 + ... + anXn  - Y)2 ) / a1 = 0  (Equation 4b)


(åiD2 ) / a2 = (å (a0 + a1X + a2X2 + a3X3 + ... + anXn  - Y)2 ) / a2 = 0  (Equation 4c)





Now we actually take the indicated derivative of each term and separate the summations. Noting that (åi F2) = 2 åi F F (where F is the inner summation of the ai X):


2 a0 å(1) + 2a1 åiX + 2a2 åiX2 + 2a3 åiX3 + ... + 2an å Xn - 2åiY = 0  (Equation 5a)


2a0 åiX +2a1 åiX2 + 2a2 åiX3 + 2a3 åiX4 + ... + 2an åi Xn+1 - 2 åiXY = 0  (Equation 5b)


2a0 åiX2 + 2a1 åiX3 + 2a2 åiX4+ 2a3 åiX5 + ... + 2an åiXn+2 - 2 åiX2Y = 0  (Equation 5c)




Dividing both sides of equations 5 by two eliminates the constant term and subtracting the term involving Y from each side of the resulting equations puts the equations in their final form:


            a0 å(1) + a1 åiX + a2 åiX2 + a3 åiX3 + ... + an åi Xn = åiY  (Equation 6a)


            a0 åiX +a1 åiX2 + a2 åiX3 + a3 åiX4 + ... + an åi Xn+1 = åiXY  (Equation 6b)


            a0 åiX2 + a1 åiX3 + a2 åiX4+ a3 åiX5 + ... + an åiXn+2 = åiX2Y  (Equation 6c)




The values of X and Y are known, since they constitute the data. Therefore equations 6 comprise a set of n + 1 equations in n + 1 unknowns, the unknowns being the various values of the ai since the summations, once evaluated, are constants. Therefore, solving equations 6 for the ai as simultaneous equations results in the calculation of the coefficients that describe the polynomial (of degree n) that best fits the data in the least squares sense.


In principle, the relationships described by equations 6 could be used directly to construct a function that relates test results to sample concentrations. In practice there is an important consideration that must be taken into account. This consideration is the possibility of correlation between the various powers of X. We find, for example, that the correlation coefficient of the integers from 1 to10 with their squares is 0.974, a rather high value.


Correlation effects are of concern to us. Our goal, recall, is to formulate a method of testing linearity in such a way that the results can be justified statistically. Ultimately we will want to perform statistical testing on the coefficients of the fitting function that we use. In fact, we will want to use a t‑test to see whether any given coefficient is statistically significant, compared to the standard error of that coefficient. We do not need to solve the general problem, however, just as we do not need to create the general solution implied by equation 1. In the broadest sense, equation 1 is the basis for computing the best-fitting function to a given set of data, but that is not our goal. Our goal is only to determine whether the data represent a linear function or not. To this end it suffices only to ascertain whether the data can be fitted better by any polynomial of degree greater than 1, than it can by a straight line (which is itself a polynomial of degree 1). To this end we need  to test a polynomial of any higher degree. While in some cases, the use of more terms may be warranted, in the limit we need test only the ability to fit the data using only one term of degree greater than one. Hence, while in general we may wish to try fitting equations of degrees 2, 3, ... m (where m is some upper limit less than n), we can begin by using polynomials of degree 2, i.e., quadratic fits.


A complication arises. We learn from considerations of multiple regression analysis, that when two (or more) variables are correlated, the standard error of both variables is increased over what would be obtained if equivalent but uncorrelated variables are used. This is discussed by Daniel and Wood (see page 55 in(9)), who show that the variance of the estimates of coefficients (their standard errors) is increased by a factor of:


            VIF = 1 / (1 - R2)    (Equation 7)


when there is correlation between the variables, where R represents the correlation coefficient between the variables and we use the term VIF, as is sometimes done, to mean Variance Inflation Factor. Thus we would like to use uncorrelated variables. Arden describes a general method for removing the correlation between the various powers of X in a polynomial, based on the use of Chebyshev or other types of orthogonal polynomials, as we briefly mentioned above. But this method is unnecessarily complicated for our current purposes. In addition, Chebyshev polynomials (along with the Legendre and other types of polynomials) have a limitation: they require that the data be uniformly (or at least symmetrically) spaced along the X-axis, a requirement that real data may not always be able to meet. Since, as we shall see, we do not need to deal with the general case, we can use a simpler method to orthogonalize the variables, based on Daniel and Wood, who showed how a variable can be transformed so that the square of that variable is uncorrelated with the variable. This is a matter of creating a new variable by simply calculating a quantity Z and subtracting that from each of the original values of X. Z is calculated using the expression (see page 121 in(9)):


                   (Equation 8)


where the summations are taken over all the samples. Then the set of values (X ‑ Z)2 will be uncorrelated with X, and estimates of the coefficients will have the minimum possible variance, making them suitable for statistical testing. This calculation also has the inherent characteristic of not imposing special requirements on the sample distribution.


In his discussion of using these approximating polynomials, Arden presents a computationally efficient method of setting up and solving the pertinent equations. But we are less concerned with abstract concepts of efficiency than we are with achieving our goal of determining linearity. To this end, we point out that the equations 6, and indeed the whole derivation of them starting from equation 1, is familiar to us, although in a different context. We are all familiar with using a relationship similar to equation 1; in using spectroscopy to do quantitative analysis, one of the representations of the equation involved is:


     C = b0 + b1X1 + b2X2 + ...       (Equation 9)


which is the form we commonly use to represent the equations needed for doing quantitative spectroscopic analysis using the MLR (Multiple Linear Regression) algorithm. The various Xi in equation 9 represent entirely different variables. Nevertheless, starting from equation 9, we can derive the set of equations for calculating the MLR calibration coefficients, in exactly the same way we derived equations 6 from equation 1. This derivation is presented in (5) and in (10).


Because of this parallelism between the situations we can set up the equivalencies:


            a0 = b0

            a1 = b1             X1 = X

            a2 = b2             X2 = X2

            a3 = b3             X3 = X3



and we see that by replacing our usual MLR-oriented variables X1, X2, X3, etc. with X, X2, X3, etc., respectively, we can use our common and well-understood mathematical methods (and computer programs) to perform the necessary calculations. Furthermore, along with the values of the coefficients, we can obtain all the usual statistical estimates of variances, standard errors, goodness of fit, etc. that MLR programs produce for us. Of special interest is the fact that MLR programs compute estimates of the standard errors of the coefficients, as described by Draper and Smith (see, for example, page 129 in (5)). This allows testing the statistical significance of each of the coefficients, which, as we recall, are now the coefficients of the various powers of X that comprise the polynomial we are fitting to the data.


This is the basis of our tests for non-linearity. We need not use polynomials of high degree since our goal is not necessarily to fit the data as well as possible. Especially since we expect that well-behaved methods of chemical analysis will produce results that are already close to linearly related to the analyte concentrations, we expect non-linear terms to decrease as the degree of the fitting equation used increases. Thus we need only fit a quadratic, or at most a cubic equation to our data to test for linearity, although there is nothing to stop us from using equations of higher degree if we choose. Data well-described by a linear equation will produce a set of coefficients with a statistically-significant value for the X1 term and non-significant values for the coefficients of X2 and (if used) X3 or higher.


One recipe for performing the calculations can be expressed as follows:


1) We start with a set of test results (X) and corresponding analyte concentrations (Y)


2) Compute Z according to equation 8


3) From the set of test results (X) create a set of numbers from the values of X by calculating (X ‑ Z)2. Consider this set of values as a new variable, X2.


4) Perform Multiple Linear Regression calculations of X and X2 against Y, the analyte concentrations. The program should calculate the coefficients of X and X2, regression statistics to evaluate the overall fit of the model to the data, calculate the standard error of the coefficients and from those the t-values for the coefficients.


As an example, these concepts are applied to the Anscombe data(2). Figure 3 reproduces two of the plots from Anscombe’s paper, the ones showing a “normal”, linear relationship and the one showing non-linearity. We recall that the linear regression statistics for both sets of data were essentially identical. Table 1 shows the results of applying this to both the “normal” data (Anscombe’s X1, Y1 set) and the data showing non-linearity. We computed the nature of the fit using both a straight-line (linear) fit only as was done originally by Anscombe, and we also fitted a polynomial using the quadratic term a well. It is interesting to compare results for the four cases.


We find that in all four cases, the coefficient of the linear term is 0.5. In Anscombe’s original paper, this is all he did, and obtained the same result, but this was by design. Using the polynomial, the fact that we obtained the same coefficient demonstrates that the quadratic term was indeed uncorrelated to the linear term.


The improvement in the fit from the quadratic polynomial indicated that the square term was indeed an important factor in fitting the data. The coefficient obtained for the quadratic term is comparable in magnitude to the linear term, as we might expect from the amount of curvature of the line we see in Anscombe’s plot(2). The coefficient of the quadratic term for the “normal” data is much smaller than for the linear term.


As we expected, for the “normal”, linear relationship, the t-value for the quadratic term for the linear data is not statistically significant. This demonstrates our contention that this method of testing linearity is indeed capable of distinguishing the two cases, in a manner that is statistically justifiable.


The performance statistics, the SEE and the correlation coefficient, show that including the square term in the fitting function for Anscombe’s non-linear data set, gives essentially an almost perfect fit. Indeed, the fit is so good that it is fair to say that we have probably reproduced his generating function: the values of the coefficients obtained are probably the ones he used to create the data in the first place. The small SEE and very large t-values of the coefficients are indicative of the fact that we are near to having only computer round-off error as operative in the difference between the data he provided and the values calculated from the polynomial that included the second-degree term.


This is the basis for our new test of linearity. It has all the advantages we described: it gives an unambiguous determination of whether any non-linearity is affecting the relationship between the test results and analyte concentration. It provides a means of distinguishing between different types of non-linearity, if they are present, since only those that have statistically-significant coefficients are active.



This new linearity test provides all the statistical tests that the current FDA/ICH test procedure recommends. and it also provides information as to whether, and how well, the analytical method gives a good fit of the test results to the actual concentration values. It can distinguish between different types of non-linearities, if necessary, while simultaneously evaluating the overall goodness of the fitting function. As the results from applying it to the Anscombe data show, it is eminently suited to evaluating the linearity characteristics of small data set as well al large ones.


This new test also provides all the statistical tests that the current FDA/ICH test procedure recommends. and therefore also provides information as to whether, and how well, the analytical method gives a good fit of the test results to the actual concentration values.



Table 1 - the results of applying the new method of detecting non-linearity to Anscombe’s data sets, both the linear and non-linear, as described in the text.



Coefficient when using only linear term

t-values for the coefficients using only linear term


Coefficient including square term

t-values for the coefficients  including the square term







Results for non-linear data (figure 3B)












Linear term






Square term






























Results for linear data (figure 3A)


















Linear term






Square term




























Figure 1 - A representation of linear data


Figure 2 - A quadratic polynomial can provide a better fit to a nonlinear function over a given region than a straight line can; in this case the second derivative of a Normal absorbance band


Figure 3 - Anscombe’s plots reproduced, as described in the text. Figure 1A: the data showing a linear relationship. Figure 1B: the data showing a non-linear relationship.






















1. ICH-Q2A; "Guideline for Industry: Text on Validation of Analytical Procedures" (1995)


2. Anscombe, F. J.; Amer. Stat.; 27 p.17-21 (1973)


3. Durbin, J., Watson, G. S.; Biometrika; 37 p.409-428 (1950)


4. Durbin, J., Watson, G. S.; Biometrika; 38 p.159-178 (1951)


5. Draper, N., Smith, H.; "Applied Regression Analysis"; 3 ed; John Wiley & Sons; New York (1998)


6. Ritchie, G. E., Roller, R. W., Ciurczak, E. W., Mark, H., Tso, C., MacDonald, S. A.; Journal of Pharmaceutical and Biomedical Analysis; (Part A: in press) (2002)


7. Ritchie, G. E., Roller, R. W., Ciurczak, E. W., Mark, H., Tso, C., MacDonald, S. A.; Journal of Pharmaceutical and Biomedical Analysis; (Part B: in press) (2002)


8. Arden, B. W.; "An Introduction to Digital Computing"; 1st ed; Addison-Wesley Publishing Co., Inc.; Reading, Mass. (1963)


9. Daniel, C., Wood, F.; "Fitting Equations to Data - Computer Analysis of Multifactor Data for Scientists and Engineers"; 1 ed; John Wiley & Sons; (1971)


10. Mark, H.; "Principles and Practice of Spectroscopic Calibration"; John Wiley & Sons; New York (1991)