- Email: [email protected]
Statistical analysis of second harmonic generation experiments: a phenomenological model
Chemometrics and Intelligent Laboratory Systems 75 (2005) 45 – 54 www.elsevier.com/locate/chemolab
Statistical analysis of second harmonic generation experiments: a phenomenological model A.H. Welsh a,b,*, R.A. Mansson b, J.G. Frey c, L. Danos c a
Centre for Mathematics and its Applications, The Australian National University, Canberra ACT 0200, Australia b School of Mathematics, University of Southampton, Southampton, SO17 1BJ, UK c School of Chemistry, University of Southampton, Southampton, SO17 1BJ, UK Available online 25 June 2004
Abstract We discuss issues arising in fitting theoretically derived nonlinear models with complex coefficients to data from surface Second Harmonic Generation (SHG) experiments conducted at the air/liquid interface. We explore different parametrisations for the complex parameters and show that the Euler (magnitude and phase angle) parametrisation is preferable to the real and imaginary parts parametrisation, both theoretically and empirically. We emphasise the importance and value of diagnostic plots for evaluating the quality of model fit. We derive approximate standard errors for the parameter estimates and discuss issues of making inference about ratios of parameters. We consider approximate confidence intervals (using the approximate standard errors), profile likelihood intervals, Fieller’s method and bootstrap intervals. Fieller’s method (and the bootstrap intervals) provide useful information on the value of the simpler approximate confidence intervals. We also propose and implement a likelihood ratio test to assess whether a common model can be fitted to several independent data sets. Finally, the methods are applied to data sets obtained from SHG experiments on L-phenylalanine at the air/water interface and toluene. D 2004 Elsevier B.V. All rights reserved. Keywords: Complex coefficients; Confidence intervals; Likelihood ratio test; Nonlinear least squares; Ratios of parameters
1. Introduction The investigation of the behaviour of interfaces between two phases, for example the surface of a liquid or the interface between a solid and a liquid is an important part of modern chemistry with implications for example in surfactants, catalysis, membranes and electrochemistry. Interfacial second harmonic generation (SHG) experiments provide a convenient and effective method of investigating the structure and dynamics of these interfaces even in the presence of overlying bulk phases [1– 5]. In SHG experiments second order non-linear interactions convert laser radiation at frequency x to the harmonic at frequency 2x. This process is electric dipole forbidden in media with inversion symmetry (gases and most solids and liquids) but, at the interface where inversion symmetry is absent, SHG is observed. The electric dipole allowed * Corresponding author. Centre for Mathematics and its Applications, The Australian National University, ACT 0200 Canberra, Australia. Tel.: +61-2-6125-2960; fax: +61-2-6125-5549. E-mail address: [email protected] (A.H. Welsh). 0169-7439/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.chemolab.2004.05.002
surface contribution frequently dominates the higher order bulk terms (e.g. magnetic dipole) and this interfacial selectivity is at the heart of many of the applications of SHG. The ultimate objective of many SHG experiments is to obtain the components of v(2), the surface second order susceptibility tensor, and to interpret them in terms of molecular behaviour. The superscript will be suppressed for the remainder of the paper. This process requires several stages of model building. The first step for experiments conducted at the air/liquid interface can be achieved within a phenomenological model with unknown parameters A, B and C (described in Eq. (1) below) which requires only that the processes being observed are due to a second harmonic process from an azimuthally isotropic surface [6,7]. In this paper we consider the problem of fitting this phenomenological model, a process that entails estimating and making inferences about A, B and C. The SHG experiments that provide the data for this paper were carried out in the Chemistry Department at the University of Southampton. The behaviour of phenylalanine at the air/water interface and toluene at the air/toluene
46
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
interface were investigated using the SHG laser system. The details of the experimental apparatus are given elsewhere [6,8 – 11]. Briefly, pulsed laser radiation at 532 nm is focussed to and reflected from the liquid interface and the intensity of the small fraction of harmonic radiation generated at 266 nm is recorded as a function of the input polarisation angle c of the fundamental beam and the output polarisation angle C of the harmonic beam. The angles c and C define the angle the electric field vector makes with the plane of incidence for the input and harmonic beams, respectively. The polarisation dependence of the SHG signal provides information about the surface nonlinear susceptibility terms and these can be used to investigate the orientation of the molecules at the interface. For each pair of polarisation angles cj and Ck selected by the experimenter, the laser is fired 2000 times and the intensity of the reflected light recorded after each shot. The averaged background intensity is subtracted from the 2000 measurements which are then averaged to produce the response Ykj, called the observed intensity. In addition, we compute the sample variance of the 2000 measurements and define the weight wkj to be the number of shots (2000) divided by the sample standard deviation. After repeating this procedure over a grid of gk fundamental and G harmonic polarisation angles, the data available for modelling are (cj, Ck, Ykj, wkj), j = 1,. . ., gk and k = 1,. . .,G. In many experiments, the weights are nearly constant and, in such cases, it is convenient to set them all equal to 1. Let n = SkG= 1gk denote the total number of observations. The observed intensity Ykj is an observation on the theoretical intensity AEkjA2, where, for the molecules considered in this paper, we can write Ekj as Ekj ¼ Es; j sinðCk Þ þ Ep; j cosðCk Þ;
ð1Þ
parts are zero far from resonance but become significant on resonance). It is convenient to express or parametrize the complex parameters A, B and C as three pairs of real coefficients. The two obvious choices of parametrisation are the Euler parametrisation in which we express the complex coefficients in terms of their magnitudes and phase angles as A = ra exp (i/a), B = rb exp(i/b), and C = rc exp(i/c) and the real and imaginary parts parametrisation in which we write A = xa + iya, B = xb + iyb, and C = xc + iyc. The corresponding sets of real coefficients are (ra, rb, rc, /a, /b, /v)T and (xa, xb, xc, ya, yb, yc)T. In either case, the model is overparametrised because the overall phase of the experiment is not determined. This means that only the relative phases of the parameters can be determined. In practical terms, a nonestimable constraint needs to be applied to the model as there are at most five independent parameters not 6: we set ya = /a = 0 so that A is a real coefficient. We also note that there are additional issues of parameter identification arising from the fact that the theoretical intensity is invariant to taking complex conjugates of the parameters, changing the signs of the magnitudes and to adding 2p to the phase angles. There is a well-known one-to-one mapping between the two parametrisations so they are mathematically equivalent and can often be used interchangeably. It is not straightforward to compare their relative merits and we will defer discussion of this issue until after we have discussed estimation and inference for the coefficients. We denote both vectors of the five real coefficients by Q when referring to either set of coefficients. If we make the assumption that the variability in the data can be represented by independent normally distributed random variables, the model can then be written in the form
with 2
2
Ep; j ¼ A cos ðcj Þ þ B sin ðcj Þ Es; j ¼ C sinð2cj Þ: The unknown parameters A, B and C are complex because they are defined in terms of the non-zero susceptibility tensor components of the surface and the Fresnel coefficients, which are functions of the angles and refractive indicies of the media. The presence of a molecule which absorbs light at the fundamental and/or the harmonic wavelength ensures that the refractive index and the nonlinear susceptibilities are complex. SHG experiments often suffer from low signals and to improve the sensitivity of the experiments the laser wavelength (or the harmonic) is chosen to be resonant with an electronic absorption of the interfacial molecules being investigated. The resonant SHG gives significantly larger signals but the surface second order susceptibility tensor becomes complex (the imaginary
Ykj ¼ AEkj ðQÞA2 þ
ekj ; wkj
j ¼ 1; . . . ; gk ; k ¼ 1; . . . ; G;
where ekj f independent N (0,r2). In all the experiments considered in this paper, gk = 14 and G = 3 so that n = 42. These quantities were chosen to compromise between obtaining detailed information about functional form and taking so many shots that the conditions of the experiment change. The output polarisation angles Ck were set to C1 = 0j (P-polarised data), C2 = 90j (S-polarised data), and another arbitrary output angle C3 (often 45j). In all these cases, the model (1) takes on a simple form. The 14 input polarisation angles were taken equally spaced on [0j, 90j] for all three output polarisation angles. This uniform design is reasonable in the absence of information about the parameters in the model. However, once some knowledge of the parameters is available, the design can be improved on; we will pursue this problem in later work. We consider model fitting and the computation of approximate standard errors for the estimated parameters qˆ in
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
Section 2. We compare various methods for making inferences about ratios of the parameters A, B and C in Section 3 and then compare the Euler and the real and imaginary part parametrisations in Section 4. We apply the methodology to phenylalanine and toluene data in Sections 5 and 6 and present our conclusions in Section 7.
A natural way to estimate the parameters Q is by weighted least squares in which we minimise gk G X X
wkj ðYkj AEkj ðQÞA2 Þ2
fitting models than the sequential estimator. This unanticipated finding seems to be due to the reduced flexibility of fitting models to different parts of the data separately. For both of these reasons, we recommend the use of the simultaneous estimation method rather than the sequential method. However we fit the model, the residuals are pffiffiffiffiffiffi ˆ 2 Þ: w ðY AE ðqÞA kj
2. Estimation
ð2Þ
47
kj
kj
Standard diagnostic methods such as residual plots, normal probability plots, etc. can be used to explore the quality of the fitted model; see for example [12]. Additional diagnostics of particular relevance to nonlinear models are discussed in Section 4 below.
k¼1 j¼1
over Q. Given reasonable initial values for Q, we can use a standard Gauss– Newton algorithm to compute the minimum Qˆ . When the response is scaled the observed intensities are between 0 and 1, and the initial magnitudes in the Euler parametrisation can be any values in [0,1]. Also, the theoretical intensities are periodic functions of the phase angles over [0,p], so the initial phase angles can take values in [0,p]. For both parametrisations, we simply set all parameters initially equal to 1. If we let W = diagonal(w11,w12,. . .,wGgG) be an n n diagonal matrix and V(Q) be the n 5 matrix whose rows are the partial derivatives of AEkj (Q)A2 with respect to the components of Q, then the approximate variance matrix of Qˆ is r2{V(Q)TWV(Q)} 1. This matrix can be estimated by rˆ 2 {V(Qˆ )TWV(Qˆ )} 1, where rˆ 2 is the residual sum of squares divided by the residual degrees of freedom n 5. The square roots of the diagonal terms are the approximate standard errors of Qˆ . An alternative approach motivated originally by the need to circumvent software limitations is to proceed sequentially. This approach exploits the fact that when the output polarisation angle Ck is set to C1 = 0j (Ppolarised data) and C2 = 90j (S-polarised data), the theoretical intensity for the two cases has no parameters in common so that these parameters can be estimated separately and the remaining unknown parameters estimated given these estimates. Comparing the simultaneous and sequential fitting methods, it is clear that the simultaneous method is more efficient in the sense of producing estimators with smaller variance than the sequential method because it uses the same or more data to estimate each parameter. For example, with the Euler parametrisation, the simultaneous method uses roughly twice as much data to estimate each parameter (other than /c) than the sequential method. The effect on the variance in an experiment on phenylalanine is shown in Table 2 below. In practice, our experience has also shown that the simultaneous approach can produce much better
3. Ratios of coefficients The main quantities of interest, particularly when we compare several experiments, are in fact ratios of the coefficients A, B and C. The treatment of the ratios of any pair of coefficients is similar so for definiteness we consider only the ratio between A and B. In the Euler parametrisation we have RAB ¼
ra expði/a Þ ra ¼ expðið/a /b ÞÞ ¼ rR expði/R Þ: rb expði/b Þ rb
The variance of the estimated phase angle is Varð/ˆ R Þ ¼ Varð/ˆ a /ˆ b Þ ¼ Varð/ˆ a Þ þ Varð/ˆ b Þ 2Covð/ˆ a ; /ˆ b Þ and, using Taylor series expansions, the variance of the estimated magnitude is ra 1 ra rˆa VarðˆrR Þ ¼ Var þ ðˆra ra Þ 2 ðˆrb rb Þ cVar rb rb rˆb rb ¼
1 r2 2ra Varðˆra Þ þ a4 Varðˆrb Þ 3 Covðˆra ; rˆb Þ 2 rb rb rb
ð3Þ
and the covariance of the estimated magnitude and phase angle is 1 rˆa ˆ CovðˆrR ; /ˆ R Þ ¼ Cov ; /a /ˆ b cCov ðˆra ra Þ rb rˆb ra 1 2 ðˆrb rb Þ; /ˆ a /ˆ b c Covðˆra ; /ˆ a Þ rb rb 1 r a Covðˆra ; /ˆ b Þ 2 Covðˆrb ; /ˆ a Þ rb rb ra þ 2 Covðˆrb ; /ˆ b Þ: ð4Þ rb It is helpful that the magnitude of the ratio involves only the component magnitudes, and similarly that the phase angle involves only the component phase angles.
48
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
Table 1 Values of RMS curvature for the two parametrisations, based on the phenylalanine data of 05-11-2001 pffiffiffiffi pffiffiffiffi cu F Form Variables cL F raexp(i/a) xa + iya
ra, rb, rc, /b, /c xa, xb, xc, yb, yc
0.2935 0.2935
0.5634 1.1153
To compare ratios produced using the two parametrisations, we need to calculate rR cos(/R) and rR sin(/R). Approximate variances for the estimators of these quantities are VarðˆrR cosð/ˆ R ÞÞccos2 ð/ˆ R ÞVarðˆrR Þ þ rˆR2 sin2 ð/ˆ R ÞVarð/ˆ R Þ 2ˆrR sinð/ˆ R Þcosð/ˆ R ÞCovðˆrR ; /ˆ R Þ and VarðˆrR sinð/ˆ R ÞÞcsin2 ð/ˆ R ÞVarðˆrR Þ þ rˆR2 cos2 ð/ˆ R ÞVarð/ˆ R Þ þ 2ˆrR sinð/ˆ Þcosð/ˆ ÞCovðˆrR ; /ˆ Þ: R
R
Fig. 2. Polirisation dependence of the SHG signal generated for S, P and 45j output angles as a function of the output polarisation angle for 5 mM phenylalanine solution. Fitted model included in the graph.
R
An interesting problem is to test the hypothesis that the ratio RAB is real. This is equivalent to testing the hypothesis H0: /R = pk , or alternatively we have H0: sin(/R) = 0ZH0: sin(/1 /2) = 0. Noting that the variance of sin (/ˆ R), is approximately Var(sin(/ˆR)) c cos2(/ˆR)Var(/ˆR), we can estimate the approximate variance of sin(/ˆR) and then construct an approximate test or confidence interval for sin(/R).
We may also be interested in constructing confidence intervals for rR. The delta method uses an estimate of the variance (3) and a normal approximation to construct an approximate confidence interval for rR directly. This is a reasonable interval provided the denominator in the ratio (rb) is not too small. In this case, other approaches may yield better results. We can consider:
Profile intervals (see for example Ref. [13]). We can construct a confidence interval for rR by mapping the
Fig. 1. Profile plots for the phenylalanine data, for the two parametrisations. (a) Euler parametrisation fitted simultaneously. (b) Real and imaginary parametrisation fitted simultaneously.
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54 Table 2 Parameters estimates and their approximate standard errors for the simultaneous and sequential fitting methods applied to the phenylalanine data Parameter
ra rb rc /a /b /c racos(/a) rbcos(/b) rccos(/c) rasin(/a) rbsin(/b) rcsin(/c) xa xb xc ya yb yc
Simultaneous
Sequential
Estimate
S.E.
0.888575 0.317008 0.588357
0.004526 0.013998 0.005289
Estimate 0.889022 0.326231 0.581158
0.005612 0.016661 0.007233
S.E.
1.141552 0.300660 0.888575 0.131934 0.561964
0.072592 0.102028 0.004526 0.017336 0.017102
1.245884 0.309684 0.889022 0.104141 0.553512
0.088642 0.112251 0.005612 0.024425 0.022073
0.288249 0.174242 0.888575 0.131934 0.561964
0.020615 0.057784 0.004526 0.017336 0.017102
0.288249 0.174242
0.020615 0.057784
0.309163 0.177112 0.855245 0.213150 0.581158 0.240559 0.233141
0.022744 0.061812 0.055789 0.052397 0.014066 0.210633 0.072290
coefficients (ra, rb, rc, /b, /c) to (ra/rb, rb, rc/rb, /b, /c /b)=(rR, rb, rc/rb, /R, /c /b), constructing the profile (likelihood) function and then deriving the interval from the profile function. The profile function is a function of rR alone which is constructed by fixing rR and replacing the other coefficients in the weighted sum of squares criterion (2) by their estimates computed with rR fixed. The confidence interval can either be derived directly from the profile function or from the quadratic approximation to the profile function. In the latter case, the interval is the same as
49
the direct interval derived from a Taylor expansion of the estimated ratio. Fieller’s method ([14]; see for example Ref. [15]). A confidence interval for rR can be constructed directly by solving for p in the inequality ðˆr1 pˆr2 Þ2 Vt 2 ; Varðˆr1 Þ þ p2 Varðˆr2 Þ 2pCovðˆr1 ; rˆ2 Þ or equivalently p2 ðt 2 Varðˆr2 Þ rˆ22 Þ 2pðt 2 Covðˆr1 ; rˆ2 Þ rˆ1 rˆ2 Þ þ t 2 Varðˆr1 Þ rˆ12 z0; where t is the critical value from a Student t or Gaussian distribution, and then replacing the variances and covariances by their estimates. The solution to this inequality can be the whole real line, a bounded interval of the real line, or an empty set (corresponding to complex roots). All of these solutions are meaningful and interpretable: substantial differences between this and the direct interval indicate that the direct interval is unreliable. Bootstrap Intervals ([16]). The bootstrap distribution of the weighted least squares estimators can be obtained by selecting 1000 bootstrap samples, and estimating the model parameters for each sample so that we have an empirical distribution of 1000 estimates. Each bootstrap sample is obtained by sampling independently, with replacement from the model residuals, adjusting the sampled residuals to have constant variance, centring them, and then adding these residuals to the fitted values AEkj(Qˆ )A2. The empirical quantiles of the bootstrap distribution are the endpoints of the bootstrap confidence interval.
Fig. 3. Residual and normal probability pots for the model fitted to the phenylalanine data. (a) Euler parametrisation fitted simultaneously. (b) Real and imaginary parametrisation fitted simultaneously.
50
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
Table 3 Ratios of parameters A, B, and C computed for the Euler and real and imaginary parametrisations. Approximate variance are also shown (in parenthesis), as well as 95% confidence intervals (C.I.) of the ratios calculated by Fieller’s method ˆI Ratio rˆR Fieller C.I. /ˆ R rˆR cos(/ˆR) rˆR sin(/ˆ R) Rˆ S.E. (rˆR cos (/ˆ R)) S.E. (rˆR sin (/ˆ R)) S.E. (rˆR) S.E. (/ˆR) S.E. (Rˆ) S.E. (Iˆ) A/B B/A A/C C/A B/C C/B
2.803002 (0.123023) 0.356760 (0.015658) 1.510265 (0.015226) 0.662135 (0.006675) 0.538803 (0.024156) 1.855967 (0.083208)
(2.5810, 3.0668) (0.3261, 0.3875) (1.4809, 1.5406) (0.6491, 0.6753) (0.4916, 0.5863) (1.7056, 2.0342)
1.141552 (0.072592) 1.141552 (0.072592) 0.30066 (0.102028) 0.300660 (0.102028) 0.840893 (0.114442) 0.840893 (0.114442)
1.166563 (0.222039) 0.148478 (0.019883) 1.442517 (0.048480) 0.632433 (0.020739) 0.359273 (0.038842) 1.237555 (0.195874)
Various confidence intervals for the magnitudes of ratios of coefficients are compared in Section 5 below. We can also explore the ratios of parameters in the real and imaginary part parametrisation. We now have RAB ¼
xa þ iya xa xb þ ya yb ya xb xa yb ¼ þi 2 ¼ xab þ iyab 2 2 xb þ iyb xb þ yb xb þ y2b
which is much more complicated (and less convenient) than in the Euler parametrisation. The expression for the ratio does simplify if the denominator coefficient is real ( yb = 0) but not as much as in the Euler parametrisation. Approximate variances of xˆab and yˆab are given by
Axab T Bxab ˆ VarðQab Þ Varðˆxab Þc and BQab BQab T Byab Byab VarðQˆ ab Þ ; Varðˆyab Þc BQab BQab where Qab=(xa, xb, ya, yb) and Qˆ ab=(xˆa, xˆb, yˆa, yˆb).
4. The choice of parametrisation The choice of parametrisation in a nonlinear model like (1) is subtle and important. Both the Euler parametrisation and the real and imaginary part parametrisation are readily interpretable. The Euler parametrisation is simpler and easier to work with when we examine ratios of parameters and has been more commonly used in this type of work. However, the Euler parametrisation represents complex numbers as products of (functions of) coefficients whereas the real and imaginary part parametrisation represents complex numbers as sums of coefficients. This means that the derivatives of the theoretical intensity may be simpler in the real and imaginary part parametrisation and hence that the quadratic approximation to the weighted least squares criterion (2) on which estimation and inference are based
2.548716 (0.085065) 0.324395 (0.022815) 0.447266 (0.147055) 0.196092 (0.064640) 0.401537 (0.053638) 1.383138 (0.116923)
1.166563 (0.114565) 0.148478 (0.021771) 1.442517 (0.048480) 0.632433 (0.20739) 0.359272 (0.038842) 1.237554 (0.195874)
2.548716 (0.128744) 0.324395 (0.002673) 0.447264 (0.147055) 0.196091 (0.064640) 0.401537 (0.053638) 1.383139 (0.116923)
works better. For any particular set of data, we can explore empirically the quality of the approximations by assessing the curvature and examining profiles of the weighted least squares surface. Bates and Watts [17] identify two sources of curvature: intrinsic curvature (cL) which is determined by the data and parameter-effects (cQ) curvature which depends on the parametrisation of the model. The MASS (Modern Applied Statistics with S-PLUS) library which is available for Splus or R ([18]1), computes these measures and adjusts them with a critical value from the F distribution. Smaller values of curvature are better. Table 1 shows the values of these statistics for the phenylalanine data collected on 05-112001. The parameter effects curvature is smaller for the Euler parametrisation than for the real and imaginary part parametrisation. In addition to measures of curvature, we can explore the profiles of the sum of squares functions to compare the two parametrisations. The profile functions are obtained as described in the previous section by fixing one of the model parameters and minimising the weighted sum of squares function over the remaining parameters. These profile functions are summarised in Fig. 1a and b for the phenylalanine data. The ideal plots are V-shaped. These profile plots suggest that the linear approximation for the Euler parametrisation is reasonable. However, there is curvature in the profiles for some of the parameters when the xa + iya parametrisation is considered. Similar results obtain from other data sets we have considered, suggesting that the Euler parametrisation is preferable to the real and imaginary part parametrisation Table 4 Estimates of ratios of parameters from the model using the profile approach
rb/ra rc/ra ra/rb rc/rb /c /b
Estimate
S.E.
Profile C.I.
0.3568 0.6621 2.8031 1.8560 0.8408
0.015658 0.006675 0.123030 0.083210 0.114440
(0.3246, 0.3865) (0.6487, 0.6756) (2.5880, 3.0810) (1.7100, 2.0440) ( 1.0300, 0.6278)
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
51
Table 5 Bootstrap estimates for the parameters and ratios of parameters from the model Parameter
Summary statistics
ra rb rc /b /c rb/ra rc/ra ra/rb rc/rb /c /b
Empirical quantiles
Mean
Bias
S.E.
2.5%
5%
95%
97.5%
0.8885 0.3176 0.5890 1.1354 0.2848 0.3568 0.6629 2.8114 1.8630 0.8399
0.000139 0.001073 0.000851 0.005271 0.015184 0.000458 0.000839 0.004979 0.005017 0.000926
0.004315 0.012835 0.004959 0.071922 0.106077 0.014770 0.006387 0.114420 0.077270 0.114020
0.87987 0.29050 0.57963 0.98688 0.05023 0.32694 0.65084 2.6067 1.7247 1.0700
0.8815 0.2950 0.5808 1.0127 0.1093 0.3310 0.6524 2.6381 1.7451 1.0279
0.8954 0.3377 0.5970 1.2405 0.4345 0.3796 0.6734 3.0071 1.9952 0.6476
0.8967 0.3411 0.5984 1.2593 0.4528 0.3831 0.6750 3.0446 2.0211 0.6107
and we generally recommend using the Euler parametrisation. Of course, if we use the Euler parametrisation, we can still make inferences about the real and imaginary parts using the fact that the real part is RA = racos(/a) and the imaginary part is IA = rasin(/a). The variance of the estimator of the real part of A is approximately h i VarðRˆ A Þc cosð/ˆ a Þ rˆa sinð/ˆ a Þ 2 32 Varðˆra Þ cosð/ˆ a Þ Covðˆra ; /ˆ a Þ 6 76 4 54 ˆ ˆ Covð/ ; rˆa Þ Varð/ Þ ˆra sinð/ˆ a
5. Phenylalanine data Data were collected from an experiment on a 5 10 3 mol/dm3 solution of phenylalanine and nonlinear models were fitted to these data for the two parametrisations and the simultaneous and sequential estimation approaches. In this experiment we have G = 3 and gk = 14 which gives a total of n = 42 observations. The model for the Euler parametrisation fitted simultaneously is shown in Fig. 2 and coefficient estimates for the models fitted to the phenylalanine data with approximate standard errors are given in Table 2. Diagnostic plots for the models are shown in Fig. 3a and b. When a fitted model is considered appropriate for the data, we can investigate ratios of the parameters A, B, and C. See Section 3 for the details of the calculations. Ratios were investigated for the two models where the parameters are fitted simultaneously only, because the diagnostics for the models fitted sequentially showed deficiencies. The components of the ratios for the parameters A, B, and C are given in Table 3 with their approximate variances using the linear approximation. When considering the Euler parametrisation, the magnitude and phase angle component of the ratio ˆ can be considered separately, and estimates of rˆR and /R and their approximate variances are given in Table 3. The Fieller method produces a bounded set for all of the six ratios of magnitudes. Confidence intervals
3 7 5
aÞ
a
¼ cos2 ð/ˆ a ÞVarðˆra Þ þ rˆa2 sin2 ð/ˆ a ÞVarð/ˆ a Þ 2ˆra sinð/ˆ a Þcosð/ˆ a ÞCovðˆra ; /ˆ a Þ and, similarly, the imaginary part has approximate variance given by VarðIˆA Þcsin2 ð/ˆ a ÞVarðˆra Þ þ rˆa2 cos2 ð/ˆ a ÞVarð/ˆ a Þ þ 2ˆra sinð/ˆ Þcosð/ˆ ÞCovðˆra ; /ˆ Þ: a
a
a
Thus we can map the Euler parametrisation to the real and imaginary part parametrisation, enabling us to make inferences on either scale.
Table 6 Parameter estimates and standard errors (in parenthesis) for the models fitted to the ten toluene data sets. There are multiple experiments performed on two of the days Model parameters Data set 10/04/02 19/07/02 19/07/02 19/07/02 20/07/02 20/07/02 20/07/02 20/07/02 20/07/02 25/06/02
No.
ra
1 2 3 1 2 3 4 5
0.645474 0.349865 0.791656 0.261566 0.508602 0.750223 0.461296 0.803142 0.770863 0.682656
(0.022578) (0.005985) (0.005239) (0.003881) (0.002075) (0.003875) (0.005001) (0.005612) (0.006204) (0.007085)
rb
rc
0.137228 (0.056064) 0.087574 (0.014975) 0.227225 (0.011799) 0.069407 (0.016164) 0.084855 (0.013263) 0.154748 (0.021957) 0.116854 (0.021051) 0.139573 (0.010717) 0.187218 (0.027238) 0.191041 (0.014654)
0.417116 (0.032536) 0.274319 (0.007299) 0.705866 (0.00538) 0.211004 (0.004794) 0.177936 (0.00495) 0.418858 (0.005459) 0.397754 (0.004840) 0.437813 (0.004841) 0.443154 (0.009007) 0.417641 (0.011540)
/b
/c 0.0870644 (0.162950) 0.459937 (0.127292)
1.65588 (0.203921) 0.980186 (0.199370) 0.839556 (0.219400) 1.432762 (0.197714)
0.683764 (0.122716)
1.224589 (0.182336) 0.423377 (0.172441)
52
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
constructed from the profile approach are given in Table 4. The bootstrap results given in Table 5 are based on 1000 samples, and the bias for the parameter estimates is not very large for any of the parameters of interest. The intervals are all quite similar in this example of phenylalanine data.
6. Toulene data We examined data from 10 experiments on the toluene/air interface. Each experiment considered two, three, or four output polarisation angles and 14 input polarisation angles for each of the output angles. The presence of fluoresence in
Fig. 4. Residual plots for the models fitted to the toluene data.
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
toluene can interfere with the SHG signal leading to distortion of the SHG results. The SHG technique is sensitive to organic impurities present at the air/liquid interface and these impurities have an impact on observed intensities. Nonlinear models were fitted to these data sets and, where possible, the model was simplified by excluding any redundant parameters. Specifically, the parameter with the largest nonsignificant p-value (using approximate t-tests on the parameter estimates) was dropped from the model, the model refitted and the process repeated until all the parameters remaining in the model were significant. The final model for each experiment is reported in Table 6. The data are on different scales, so the parameter estimates vary between the experiments: we are interested in comparing ratios of the parameters A, B, and C from different experiments. Residual (see Fig. 4) and normal probability plots for the ten experiments showed systematic patterns for the majority of the models. To improve the fit of the model to the data we considered shifting and scaling the intensities. However, this did not improve the diagnostics in general. The residual plots suggested some potential outliers (or influential points) in the data sets, and we considered removing these and refitting the nonlinear models. In some cases the fitted model was worse, because we had removed points from informative regions of the intensity curves. The fit of the model was deemed adequate in four of the experiments so these four models were used to illustrate the procedure for comparing experiments. The coefficient with the largest magnitude was A, so the other two coefficients were compared to A in the ratios. That is, the two ratios of interest for the four experiments are B/A and C/A, and the two components of these ratios rˆR and /ˆR are given in Table 7 with confidence intervals for both components using the delta method and a Fieller confidence interval for the ratio of magnitudes. The Fieller confidence intervals are again very similar to those calculated by the delta method. We are interested in testing whether the ratios of coefficients are the same across the experiments. A likelihood
53
ratio test can be constructed to test this hypothesis. Under the null hypothesis there are four separate scale parameters and four common parameters making eight model parameters, so we are testing H0: f ðQÞ
4 Y
f ðrAs ; rBA ; rCA ; /B ; /C Þ;
s¼1
where rBA is the ratio of B and A common to all four experiments, against the alternative of a separate model for each data set with 4 5 = 20 parameters HA: f ðQÞ
4 Y
f ðrAs ; rBAs ; rCAs ; /Bs ; /Cs Þ;
s¼1
where rBAs, s = 1,. . .,4 are the ratios of B to A for the four experiments. The likelihood ratio statistic for this test where the four experiments have n1, n2, n3, and n4 observations respectively is ! ðrˆ 21 Þn1 =2 ðrˆ 22 Þn2 =2 ðrˆ 23 Þn3 =2 ðrˆ 24 Þn4 =2 2ln ðrˆ 2 ÞN =2 where N = n1 + n2 + n3 + n4 is the total number of observations over the four experiments, rˆ 2 is estimated from the model under the null hypothesis, and rˆ 2k, k = 1,. . .,4 are estimates of rˆ 2 for the models fitted to each experiment separately. The test statistic is compared to the v2 distribution with 20 8 = 12 degrees of freedom. One of the four experiments detailed in Table 7 has a structural zero for one of the phase angles so the full model has 19 rather than 20 parameters. The test statistic is referred to a v2 distribution with 11 degrees of freedom. The likelihood ratio test statistic is 234.56 which suggests that the common model is not appropriate for these four data sets. There is no theoretical reason for the differences between the different sets of toluene data. The production of reliable standard errors and confidence intervals for ratios of parameters has highlighted the existence of these differences and
Table 7 Ratios of parameters for the toluene models. Confidence intervals are given for the ratios of the magnitudes using both the delta and Fieller method. Some of the phase angles (/ˆ R) are zero because the fitted models were simplified when the phase angles were considered negligible Data set
19/07/2002 Delta Fieller 20/07/2002 Delta Fieller 20/07/2002 Delta Fieller 25/06/2002 Delta Fieller
No.
3 C.I. C.I. 2 C.I. C.I. 3 C.I. C.I. C.I. C.I.
B/A
C/A
rˆR
/Rˆ
rˆR
/Rˆ
0.265353 (0.1436, 0.3871) (0.1440, 0.3875) 0.206270 (0.1492, 0.2634) (0.1491, 0.2633) 0.253317 (0.1646, 0.3420) (0.1644, 0.3418) 0.279850 (0.2363, 0.3234) (0.2365, 0.3237)
1.655880 (1.2562, 2.0556)
0.806696 (0.7641, 0.8493) (0.7647, 0.8500) 0.558311 (0.5426, 0.5740) (0.5427, 0.5741) 0.862253 (0.8347, 0.8898) (0.8351, 0.8902) 0.611789 (0.5770, 0.6466) (0.5772, 0.6468)
0.683764 (0.4432, 0.9243)
0.839556 (0.4095, 1.2696) 1.432762 (1.0452, 1.8203) 0
0
–
0.423377 (0.0854, 0.7614)
54
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
enabled experimenters to address them more rigorously. In particular, the findings presented in this section have led to refinements of the experimental technique which have increased the reproducibility of the results. These experiments will be reported elsewhere. If the differences persist after substantial effort at refinement, then we need to consider the appropriateness of the model and the theory from which it is derived. Whatever the outcome, the analysis presented here will have been an important step towards scientific advance and understanding. 7. Conclusions We have explored statistical issues which arise in fitting a phenomenological model to data from second harmonic generation experiments and then using the model to make inferences. The key issues in model fitting involve the choice of parametrisation for the complex parameters in the model, the method used to fit the model and the use of diagnostics to evaluate the quality of model fit. The key issues in inference are the calculation of standard errors, setting confidence intervals for ratios of parameters and comparing models fitted to separate experiments. The Euler parametrisation has been found to be the better representation of complex numbers. The curvature plots to assess the linear approximation are better, and it is easier to form ratios of the parameters A, B, and C and make inferences on these ratios. Fitting the model simultaneously to all the data is better than the sequential approach, because it provides a better fit (less bias) and is more efficient (standard errors are smaller than in the sequential method). Standard residual and normal probability plots supplemented by profile plots to explore curvature in the model provide useful insights into the quality of model fit. We have derived approximate standard errors for the parameter estimates and various useful transformations of the parameter estimates. The approximate standard errors of the ratios of parameters can be used to set approximate confidence intervals provided the denominator estimate is not too small. These intervals can be compared to Fieller intervals which do not use the same kinds of approximation
but may be the whole real line, a proper interval or the empty set. Substantial differences between the two types of intervals suggest that the approximation behind the simple approximate standard errors is poor. We also considered alternative profile likelihood and bootstrap approaches. All these methods performed very similarly in the examples we considered (supporting the use of the simplest method in these examples) but there is useful information when they differ. Finally, we developed and implemented a likelihood ratio test for comparing models fitted to separate experiments.
References [1] Y.R. Shen, Annual Review of Physical Chemistry 40 (1989) 327 – 350. [2] K.B. Eisenthal, Annual Review of Physical Chemistry 43 (1992) 627 – 661. [3] K.B. Eisenthal, Accounts of Chemical Research 26 (12) (1993) 636 – 643. [4] Y.R. Shen, Surface Science 300 (1 – 3) (1994) 551 – 562. [5] Y.R. Shen, Solid State Communications 102 (2 – 3) (1997) 221 – 229. [6] A.J. Fordyce, W.J. Bullock, et al., Molecular Physics 99 (8) (2001) 677 – 687. [7] V. Mizrahi, E. Sipe, Journal of the Optical Society of America. B, Optical Physics 5 (3) (1988) 660 – 667. [8] A.J. Bell, J.G. Frey, et al., Journal of the Chemical Society. Faraday Transactions 88 (14) (1992) 2027 – 2030. [9] M.J. Crawford, S. Haslam, et al., Chemical Physics Letters 230 (3) (1994) 260 – 264. [10] M.J. Crawford, J.G Frey, et al., Journal of the Chemical Society. Faraday Transactions 92 (8) (1996) 1369 – 1373. [11] S. Haslam, S.G. Croucher, et al., Physical Chemistry Chemical Physics 2 (14) (2000) 3235 – 3245. [12] S. Chatterjee, A.S. Hadi, B. Price, Regression Analysis by Example, 3rd edition, Wiley, New York, 2000, Chap. 4. [13] A.H. Welsh, Aspects of Statistical Inference, Wiley, New York, 1996, pp. 97 – 98. [14] E.C. Fieller, Supplement to The Journal of the Royal Statistical Society 7 (1940) 1 – 64. [15] A.H. Welsh, Aspects of Statistical Inference, Wiley, New York, 1996, pp. 152 – 155. [16] A.C. Davison, D.V. Hinkley, Bootstrap Methods and their Application, Cambridge Univ. Press, 1997, pp. 353 – 358. [17] D.M. Bates, D.G. Watts, Nonlinear Regression Analysis and its Applications, Wiley, New York, 1988, Chap. 7. [18] W.N. Venables, B.D. Ripley, Modern Applied Statistics with S, 4th edition, Springer, New York, 2002, p. 461.
Statistical analysis of second harmonic generation experiments: a phenomenological model A.H. Welsh a,b,*, R.A. Mansson b, J.G. Frey c, L. Danos c a
Centre for Mathematics and its Applications, The Australian National University, Canberra ACT 0200, Australia b School of Mathematics, University of Southampton, Southampton, SO17 1BJ, UK c School of Chemistry, University of Southampton, Southampton, SO17 1BJ, UK Available online 25 June 2004
Abstract We discuss issues arising in fitting theoretically derived nonlinear models with complex coefficients to data from surface Second Harmonic Generation (SHG) experiments conducted at the air/liquid interface. We explore different parametrisations for the complex parameters and show that the Euler (magnitude and phase angle) parametrisation is preferable to the real and imaginary parts parametrisation, both theoretically and empirically. We emphasise the importance and value of diagnostic plots for evaluating the quality of model fit. We derive approximate standard errors for the parameter estimates and discuss issues of making inference about ratios of parameters. We consider approximate confidence intervals (using the approximate standard errors), profile likelihood intervals, Fieller’s method and bootstrap intervals. Fieller’s method (and the bootstrap intervals) provide useful information on the value of the simpler approximate confidence intervals. We also propose and implement a likelihood ratio test to assess whether a common model can be fitted to several independent data sets. Finally, the methods are applied to data sets obtained from SHG experiments on L-phenylalanine at the air/water interface and toluene. D 2004 Elsevier B.V. All rights reserved. Keywords: Complex coefficients; Confidence intervals; Likelihood ratio test; Nonlinear least squares; Ratios of parameters
1. Introduction The investigation of the behaviour of interfaces between two phases, for example the surface of a liquid or the interface between a solid and a liquid is an important part of modern chemistry with implications for example in surfactants, catalysis, membranes and electrochemistry. Interfacial second harmonic generation (SHG) experiments provide a convenient and effective method of investigating the structure and dynamics of these interfaces even in the presence of overlying bulk phases [1– 5]. In SHG experiments second order non-linear interactions convert laser radiation at frequency x to the harmonic at frequency 2x. This process is electric dipole forbidden in media with inversion symmetry (gases and most solids and liquids) but, at the interface where inversion symmetry is absent, SHG is observed. The electric dipole allowed * Corresponding author. Centre for Mathematics and its Applications, The Australian National University, ACT 0200 Canberra, Australia. Tel.: +61-2-6125-2960; fax: +61-2-6125-5549. E-mail address: [email protected] (A.H. Welsh). 0169-7439/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.chemolab.2004.05.002
surface contribution frequently dominates the higher order bulk terms (e.g. magnetic dipole) and this interfacial selectivity is at the heart of many of the applications of SHG. The ultimate objective of many SHG experiments is to obtain the components of v(2), the surface second order susceptibility tensor, and to interpret them in terms of molecular behaviour. The superscript will be suppressed for the remainder of the paper. This process requires several stages of model building. The first step for experiments conducted at the air/liquid interface can be achieved within a phenomenological model with unknown parameters A, B and C (described in Eq. (1) below) which requires only that the processes being observed are due to a second harmonic process from an azimuthally isotropic surface [6,7]. In this paper we consider the problem of fitting this phenomenological model, a process that entails estimating and making inferences about A, B and C. The SHG experiments that provide the data for this paper were carried out in the Chemistry Department at the University of Southampton. The behaviour of phenylalanine at the air/water interface and toluene at the air/toluene
46
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
interface were investigated using the SHG laser system. The details of the experimental apparatus are given elsewhere [6,8 – 11]. Briefly, pulsed laser radiation at 532 nm is focussed to and reflected from the liquid interface and the intensity of the small fraction of harmonic radiation generated at 266 nm is recorded as a function of the input polarisation angle c of the fundamental beam and the output polarisation angle C of the harmonic beam. The angles c and C define the angle the electric field vector makes with the plane of incidence for the input and harmonic beams, respectively. The polarisation dependence of the SHG signal provides information about the surface nonlinear susceptibility terms and these can be used to investigate the orientation of the molecules at the interface. For each pair of polarisation angles cj and Ck selected by the experimenter, the laser is fired 2000 times and the intensity of the reflected light recorded after each shot. The averaged background intensity is subtracted from the 2000 measurements which are then averaged to produce the response Ykj, called the observed intensity. In addition, we compute the sample variance of the 2000 measurements and define the weight wkj to be the number of shots (2000) divided by the sample standard deviation. After repeating this procedure over a grid of gk fundamental and G harmonic polarisation angles, the data available for modelling are (cj, Ck, Ykj, wkj), j = 1,. . ., gk and k = 1,. . .,G. In many experiments, the weights are nearly constant and, in such cases, it is convenient to set them all equal to 1. Let n = SkG= 1gk denote the total number of observations. The observed intensity Ykj is an observation on the theoretical intensity AEkjA2, where, for the molecules considered in this paper, we can write Ekj as Ekj ¼ Es; j sinðCk Þ þ Ep; j cosðCk Þ;
ð1Þ
parts are zero far from resonance but become significant on resonance). It is convenient to express or parametrize the complex parameters A, B and C as three pairs of real coefficients. The two obvious choices of parametrisation are the Euler parametrisation in which we express the complex coefficients in terms of their magnitudes and phase angles as A = ra exp (i/a), B = rb exp(i/b), and C = rc exp(i/c) and the real and imaginary parts parametrisation in which we write A = xa + iya, B = xb + iyb, and C = xc + iyc. The corresponding sets of real coefficients are (ra, rb, rc, /a, /b, /v)T and (xa, xb, xc, ya, yb, yc)T. In either case, the model is overparametrised because the overall phase of the experiment is not determined. This means that only the relative phases of the parameters can be determined. In practical terms, a nonestimable constraint needs to be applied to the model as there are at most five independent parameters not 6: we set ya = /a = 0 so that A is a real coefficient. We also note that there are additional issues of parameter identification arising from the fact that the theoretical intensity is invariant to taking complex conjugates of the parameters, changing the signs of the magnitudes and to adding 2p to the phase angles. There is a well-known one-to-one mapping between the two parametrisations so they are mathematically equivalent and can often be used interchangeably. It is not straightforward to compare their relative merits and we will defer discussion of this issue until after we have discussed estimation and inference for the coefficients. We denote both vectors of the five real coefficients by Q when referring to either set of coefficients. If we make the assumption that the variability in the data can be represented by independent normally distributed random variables, the model can then be written in the form
with 2
2
Ep; j ¼ A cos ðcj Þ þ B sin ðcj Þ Es; j ¼ C sinð2cj Þ: The unknown parameters A, B and C are complex because they are defined in terms of the non-zero susceptibility tensor components of the surface and the Fresnel coefficients, which are functions of the angles and refractive indicies of the media. The presence of a molecule which absorbs light at the fundamental and/or the harmonic wavelength ensures that the refractive index and the nonlinear susceptibilities are complex. SHG experiments often suffer from low signals and to improve the sensitivity of the experiments the laser wavelength (or the harmonic) is chosen to be resonant with an electronic absorption of the interfacial molecules being investigated. The resonant SHG gives significantly larger signals but the surface second order susceptibility tensor becomes complex (the imaginary
Ykj ¼ AEkj ðQÞA2 þ
ekj ; wkj
j ¼ 1; . . . ; gk ; k ¼ 1; . . . ; G;
where ekj f independent N (0,r2). In all the experiments considered in this paper, gk = 14 and G = 3 so that n = 42. These quantities were chosen to compromise between obtaining detailed information about functional form and taking so many shots that the conditions of the experiment change. The output polarisation angles Ck were set to C1 = 0j (P-polarised data), C2 = 90j (S-polarised data), and another arbitrary output angle C3 (often 45j). In all these cases, the model (1) takes on a simple form. The 14 input polarisation angles were taken equally spaced on [0j, 90j] for all three output polarisation angles. This uniform design is reasonable in the absence of information about the parameters in the model. However, once some knowledge of the parameters is available, the design can be improved on; we will pursue this problem in later work. We consider model fitting and the computation of approximate standard errors for the estimated parameters qˆ in
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
Section 2. We compare various methods for making inferences about ratios of the parameters A, B and C in Section 3 and then compare the Euler and the real and imaginary part parametrisations in Section 4. We apply the methodology to phenylalanine and toluene data in Sections 5 and 6 and present our conclusions in Section 7.
A natural way to estimate the parameters Q is by weighted least squares in which we minimise gk G X X
wkj ðYkj AEkj ðQÞA2 Þ2
fitting models than the sequential estimator. This unanticipated finding seems to be due to the reduced flexibility of fitting models to different parts of the data separately. For both of these reasons, we recommend the use of the simultaneous estimation method rather than the sequential method. However we fit the model, the residuals are pffiffiffiffiffiffi ˆ 2 Þ: w ðY AE ðqÞA kj
2. Estimation
ð2Þ
47
kj
kj
Standard diagnostic methods such as residual plots, normal probability plots, etc. can be used to explore the quality of the fitted model; see for example [12]. Additional diagnostics of particular relevance to nonlinear models are discussed in Section 4 below.
k¼1 j¼1
over Q. Given reasonable initial values for Q, we can use a standard Gauss– Newton algorithm to compute the minimum Qˆ . When the response is scaled the observed intensities are between 0 and 1, and the initial magnitudes in the Euler parametrisation can be any values in [0,1]. Also, the theoretical intensities are periodic functions of the phase angles over [0,p], so the initial phase angles can take values in [0,p]. For both parametrisations, we simply set all parameters initially equal to 1. If we let W = diagonal(w11,w12,. . .,wGgG) be an n n diagonal matrix and V(Q) be the n 5 matrix whose rows are the partial derivatives of AEkj (Q)A2 with respect to the components of Q, then the approximate variance matrix of Qˆ is r2{V(Q)TWV(Q)} 1. This matrix can be estimated by rˆ 2 {V(Qˆ )TWV(Qˆ )} 1, where rˆ 2 is the residual sum of squares divided by the residual degrees of freedom n 5. The square roots of the diagonal terms are the approximate standard errors of Qˆ . An alternative approach motivated originally by the need to circumvent software limitations is to proceed sequentially. This approach exploits the fact that when the output polarisation angle Ck is set to C1 = 0j (Ppolarised data) and C2 = 90j (S-polarised data), the theoretical intensity for the two cases has no parameters in common so that these parameters can be estimated separately and the remaining unknown parameters estimated given these estimates. Comparing the simultaneous and sequential fitting methods, it is clear that the simultaneous method is more efficient in the sense of producing estimators with smaller variance than the sequential method because it uses the same or more data to estimate each parameter. For example, with the Euler parametrisation, the simultaneous method uses roughly twice as much data to estimate each parameter (other than /c) than the sequential method. The effect on the variance in an experiment on phenylalanine is shown in Table 2 below. In practice, our experience has also shown that the simultaneous approach can produce much better
3. Ratios of coefficients The main quantities of interest, particularly when we compare several experiments, are in fact ratios of the coefficients A, B and C. The treatment of the ratios of any pair of coefficients is similar so for definiteness we consider only the ratio between A and B. In the Euler parametrisation we have RAB ¼
ra expði/a Þ ra ¼ expðið/a /b ÞÞ ¼ rR expði/R Þ: rb expði/b Þ rb
The variance of the estimated phase angle is Varð/ˆ R Þ ¼ Varð/ˆ a /ˆ b Þ ¼ Varð/ˆ a Þ þ Varð/ˆ b Þ 2Covð/ˆ a ; /ˆ b Þ and, using Taylor series expansions, the variance of the estimated magnitude is ra 1 ra rˆa VarðˆrR Þ ¼ Var þ ðˆra ra Þ 2 ðˆrb rb Þ cVar rb rb rˆb rb ¼
1 r2 2ra Varðˆra Þ þ a4 Varðˆrb Þ 3 Covðˆra ; rˆb Þ 2 rb rb rb
ð3Þ
and the covariance of the estimated magnitude and phase angle is 1 rˆa ˆ CovðˆrR ; /ˆ R Þ ¼ Cov ; /a /ˆ b cCov ðˆra ra Þ rb rˆb ra 1 2 ðˆrb rb Þ; /ˆ a /ˆ b c Covðˆra ; /ˆ a Þ rb rb 1 r a Covðˆra ; /ˆ b Þ 2 Covðˆrb ; /ˆ a Þ rb rb ra þ 2 Covðˆrb ; /ˆ b Þ: ð4Þ rb It is helpful that the magnitude of the ratio involves only the component magnitudes, and similarly that the phase angle involves only the component phase angles.
48
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
Table 1 Values of RMS curvature for the two parametrisations, based on the phenylalanine data of 05-11-2001 pffiffiffiffi pffiffiffiffi cu F Form Variables cL F raexp(i/a) xa + iya
ra, rb, rc, /b, /c xa, xb, xc, yb, yc
0.2935 0.2935
0.5634 1.1153
To compare ratios produced using the two parametrisations, we need to calculate rR cos(/R) and rR sin(/R). Approximate variances for the estimators of these quantities are VarðˆrR cosð/ˆ R ÞÞccos2 ð/ˆ R ÞVarðˆrR Þ þ rˆR2 sin2 ð/ˆ R ÞVarð/ˆ R Þ 2ˆrR sinð/ˆ R Þcosð/ˆ R ÞCovðˆrR ; /ˆ R Þ and VarðˆrR sinð/ˆ R ÞÞcsin2 ð/ˆ R ÞVarðˆrR Þ þ rˆR2 cos2 ð/ˆ R ÞVarð/ˆ R Þ þ 2ˆrR sinð/ˆ Þcosð/ˆ ÞCovðˆrR ; /ˆ Þ: R
R
Fig. 2. Polirisation dependence of the SHG signal generated for S, P and 45j output angles as a function of the output polarisation angle for 5 mM phenylalanine solution. Fitted model included in the graph.
R
An interesting problem is to test the hypothesis that the ratio RAB is real. This is equivalent to testing the hypothesis H0: /R = pk , or alternatively we have H0: sin(/R) = 0ZH0: sin(/1 /2) = 0. Noting that the variance of sin (/ˆ R), is approximately Var(sin(/ˆR)) c cos2(/ˆR)Var(/ˆR), we can estimate the approximate variance of sin(/ˆR) and then construct an approximate test or confidence interval for sin(/R).
We may also be interested in constructing confidence intervals for rR. The delta method uses an estimate of the variance (3) and a normal approximation to construct an approximate confidence interval for rR directly. This is a reasonable interval provided the denominator in the ratio (rb) is not too small. In this case, other approaches may yield better results. We can consider:
Profile intervals (see for example Ref. [13]). We can construct a confidence interval for rR by mapping the
Fig. 1. Profile plots for the phenylalanine data, for the two parametrisations. (a) Euler parametrisation fitted simultaneously. (b) Real and imaginary parametrisation fitted simultaneously.
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54 Table 2 Parameters estimates and their approximate standard errors for the simultaneous and sequential fitting methods applied to the phenylalanine data Parameter
ra rb rc /a /b /c racos(/a) rbcos(/b) rccos(/c) rasin(/a) rbsin(/b) rcsin(/c) xa xb xc ya yb yc
Simultaneous
Sequential
Estimate
S.E.
0.888575 0.317008 0.588357
0.004526 0.013998 0.005289
Estimate 0.889022 0.326231 0.581158
0.005612 0.016661 0.007233
S.E.
1.141552 0.300660 0.888575 0.131934 0.561964
0.072592 0.102028 0.004526 0.017336 0.017102
1.245884 0.309684 0.889022 0.104141 0.553512
0.088642 0.112251 0.005612 0.024425 0.022073
0.288249 0.174242 0.888575 0.131934 0.561964
0.020615 0.057784 0.004526 0.017336 0.017102
0.288249 0.174242
0.020615 0.057784
0.309163 0.177112 0.855245 0.213150 0.581158 0.240559 0.233141
0.022744 0.061812 0.055789 0.052397 0.014066 0.210633 0.072290
coefficients (ra, rb, rc, /b, /c) to (ra/rb, rb, rc/rb, /b, /c /b)=(rR, rb, rc/rb, /R, /c /b), constructing the profile (likelihood) function and then deriving the interval from the profile function. The profile function is a function of rR alone which is constructed by fixing rR and replacing the other coefficients in the weighted sum of squares criterion (2) by their estimates computed with rR fixed. The confidence interval can either be derived directly from the profile function or from the quadratic approximation to the profile function. In the latter case, the interval is the same as
49
the direct interval derived from a Taylor expansion of the estimated ratio. Fieller’s method ([14]; see for example Ref. [15]). A confidence interval for rR can be constructed directly by solving for p in the inequality ðˆr1 pˆr2 Þ2 Vt 2 ; Varðˆr1 Þ þ p2 Varðˆr2 Þ 2pCovðˆr1 ; rˆ2 Þ or equivalently p2 ðt 2 Varðˆr2 Þ rˆ22 Þ 2pðt 2 Covðˆr1 ; rˆ2 Þ rˆ1 rˆ2 Þ þ t 2 Varðˆr1 Þ rˆ12 z0; where t is the critical value from a Student t or Gaussian distribution, and then replacing the variances and covariances by their estimates. The solution to this inequality can be the whole real line, a bounded interval of the real line, or an empty set (corresponding to complex roots). All of these solutions are meaningful and interpretable: substantial differences between this and the direct interval indicate that the direct interval is unreliable. Bootstrap Intervals ([16]). The bootstrap distribution of the weighted least squares estimators can be obtained by selecting 1000 bootstrap samples, and estimating the model parameters for each sample so that we have an empirical distribution of 1000 estimates. Each bootstrap sample is obtained by sampling independently, with replacement from the model residuals, adjusting the sampled residuals to have constant variance, centring them, and then adding these residuals to the fitted values AEkj(Qˆ )A2. The empirical quantiles of the bootstrap distribution are the endpoints of the bootstrap confidence interval.
Fig. 3. Residual and normal probability pots for the model fitted to the phenylalanine data. (a) Euler parametrisation fitted simultaneously. (b) Real and imaginary parametrisation fitted simultaneously.
50
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
Table 3 Ratios of parameters A, B, and C computed for the Euler and real and imaginary parametrisations. Approximate variance are also shown (in parenthesis), as well as 95% confidence intervals (C.I.) of the ratios calculated by Fieller’s method ˆI Ratio rˆR Fieller C.I. /ˆ R rˆR cos(/ˆR) rˆR sin(/ˆ R) Rˆ S.E. (rˆR cos (/ˆ R)) S.E. (rˆR sin (/ˆ R)) S.E. (rˆR) S.E. (/ˆR) S.E. (Rˆ) S.E. (Iˆ) A/B B/A A/C C/A B/C C/B
2.803002 (0.123023) 0.356760 (0.015658) 1.510265 (0.015226) 0.662135 (0.006675) 0.538803 (0.024156) 1.855967 (0.083208)
(2.5810, 3.0668) (0.3261, 0.3875) (1.4809, 1.5406) (0.6491, 0.6753) (0.4916, 0.5863) (1.7056, 2.0342)
1.141552 (0.072592) 1.141552 (0.072592) 0.30066 (0.102028) 0.300660 (0.102028) 0.840893 (0.114442) 0.840893 (0.114442)
1.166563 (0.222039) 0.148478 (0.019883) 1.442517 (0.048480) 0.632433 (0.020739) 0.359273 (0.038842) 1.237555 (0.195874)
Various confidence intervals for the magnitudes of ratios of coefficients are compared in Section 5 below. We can also explore the ratios of parameters in the real and imaginary part parametrisation. We now have RAB ¼
xa þ iya xa xb þ ya yb ya xb xa yb ¼ þi 2 ¼ xab þ iyab 2 2 xb þ iyb xb þ yb xb þ y2b
which is much more complicated (and less convenient) than in the Euler parametrisation. The expression for the ratio does simplify if the denominator coefficient is real ( yb = 0) but not as much as in the Euler parametrisation. Approximate variances of xˆab and yˆab are given by
Axab T Bxab ˆ VarðQab Þ Varðˆxab Þc and BQab BQab T Byab Byab VarðQˆ ab Þ ; Varðˆyab Þc BQab BQab where Qab=(xa, xb, ya, yb) and Qˆ ab=(xˆa, xˆb, yˆa, yˆb).
4. The choice of parametrisation The choice of parametrisation in a nonlinear model like (1) is subtle and important. Both the Euler parametrisation and the real and imaginary part parametrisation are readily interpretable. The Euler parametrisation is simpler and easier to work with when we examine ratios of parameters and has been more commonly used in this type of work. However, the Euler parametrisation represents complex numbers as products of (functions of) coefficients whereas the real and imaginary part parametrisation represents complex numbers as sums of coefficients. This means that the derivatives of the theoretical intensity may be simpler in the real and imaginary part parametrisation and hence that the quadratic approximation to the weighted least squares criterion (2) on which estimation and inference are based
2.548716 (0.085065) 0.324395 (0.022815) 0.447266 (0.147055) 0.196092 (0.064640) 0.401537 (0.053638) 1.383138 (0.116923)
1.166563 (0.114565) 0.148478 (0.021771) 1.442517 (0.048480) 0.632433 (0.20739) 0.359272 (0.038842) 1.237554 (0.195874)
2.548716 (0.128744) 0.324395 (0.002673) 0.447264 (0.147055) 0.196091 (0.064640) 0.401537 (0.053638) 1.383139 (0.116923)
works better. For any particular set of data, we can explore empirically the quality of the approximations by assessing the curvature and examining profiles of the weighted least squares surface. Bates and Watts [17] identify two sources of curvature: intrinsic curvature (cL) which is determined by the data and parameter-effects (cQ) curvature which depends on the parametrisation of the model. The MASS (Modern Applied Statistics with S-PLUS) library which is available for Splus or R ([18]1), computes these measures and adjusts them with a critical value from the F distribution. Smaller values of curvature are better. Table 1 shows the values of these statistics for the phenylalanine data collected on 05-112001. The parameter effects curvature is smaller for the Euler parametrisation than for the real and imaginary part parametrisation. In addition to measures of curvature, we can explore the profiles of the sum of squares functions to compare the two parametrisations. The profile functions are obtained as described in the previous section by fixing one of the model parameters and minimising the weighted sum of squares function over the remaining parameters. These profile functions are summarised in Fig. 1a and b for the phenylalanine data. The ideal plots are V-shaped. These profile plots suggest that the linear approximation for the Euler parametrisation is reasonable. However, there is curvature in the profiles for some of the parameters when the xa + iya parametrisation is considered. Similar results obtain from other data sets we have considered, suggesting that the Euler parametrisation is preferable to the real and imaginary part parametrisation Table 4 Estimates of ratios of parameters from the model using the profile approach
rb/ra rc/ra ra/rb rc/rb /c /b
Estimate
S.E.
Profile C.I.
0.3568 0.6621 2.8031 1.8560 0.8408
0.015658 0.006675 0.123030 0.083210 0.114440
(0.3246, 0.3865) (0.6487, 0.6756) (2.5880, 3.0810) (1.7100, 2.0440) ( 1.0300, 0.6278)
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
51
Table 5 Bootstrap estimates for the parameters and ratios of parameters from the model Parameter
Summary statistics
ra rb rc /b /c rb/ra rc/ra ra/rb rc/rb /c /b
Empirical quantiles
Mean
Bias
S.E.
2.5%
5%
95%
97.5%
0.8885 0.3176 0.5890 1.1354 0.2848 0.3568 0.6629 2.8114 1.8630 0.8399
0.000139 0.001073 0.000851 0.005271 0.015184 0.000458 0.000839 0.004979 0.005017 0.000926
0.004315 0.012835 0.004959 0.071922 0.106077 0.014770 0.006387 0.114420 0.077270 0.114020
0.87987 0.29050 0.57963 0.98688 0.05023 0.32694 0.65084 2.6067 1.7247 1.0700
0.8815 0.2950 0.5808 1.0127 0.1093 0.3310 0.6524 2.6381 1.7451 1.0279
0.8954 0.3377 0.5970 1.2405 0.4345 0.3796 0.6734 3.0071 1.9952 0.6476
0.8967 0.3411 0.5984 1.2593 0.4528 0.3831 0.6750 3.0446 2.0211 0.6107
and we generally recommend using the Euler parametrisation. Of course, if we use the Euler parametrisation, we can still make inferences about the real and imaginary parts using the fact that the real part is RA = racos(/a) and the imaginary part is IA = rasin(/a). The variance of the estimator of the real part of A is approximately h i VarðRˆ A Þc cosð/ˆ a Þ rˆa sinð/ˆ a Þ 2 32 Varðˆra Þ cosð/ˆ a Þ Covðˆra ; /ˆ a Þ 6 76 4 54 ˆ ˆ Covð/ ; rˆa Þ Varð/ Þ ˆra sinð/ˆ a
5. Phenylalanine data Data were collected from an experiment on a 5 10 3 mol/dm3 solution of phenylalanine and nonlinear models were fitted to these data for the two parametrisations and the simultaneous and sequential estimation approaches. In this experiment we have G = 3 and gk = 14 which gives a total of n = 42 observations. The model for the Euler parametrisation fitted simultaneously is shown in Fig. 2 and coefficient estimates for the models fitted to the phenylalanine data with approximate standard errors are given in Table 2. Diagnostic plots for the models are shown in Fig. 3a and b. When a fitted model is considered appropriate for the data, we can investigate ratios of the parameters A, B, and C. See Section 3 for the details of the calculations. Ratios were investigated for the two models where the parameters are fitted simultaneously only, because the diagnostics for the models fitted sequentially showed deficiencies. The components of the ratios for the parameters A, B, and C are given in Table 3 with their approximate variances using the linear approximation. When considering the Euler parametrisation, the magnitude and phase angle component of the ratio ˆ can be considered separately, and estimates of rˆR and /R and their approximate variances are given in Table 3. The Fieller method produces a bounded set for all of the six ratios of magnitudes. Confidence intervals
3 7 5
aÞ
a
¼ cos2 ð/ˆ a ÞVarðˆra Þ þ rˆa2 sin2 ð/ˆ a ÞVarð/ˆ a Þ 2ˆra sinð/ˆ a Þcosð/ˆ a ÞCovðˆra ; /ˆ a Þ and, similarly, the imaginary part has approximate variance given by VarðIˆA Þcsin2 ð/ˆ a ÞVarðˆra Þ þ rˆa2 cos2 ð/ˆ a ÞVarð/ˆ a Þ þ 2ˆra sinð/ˆ Þcosð/ˆ ÞCovðˆra ; /ˆ Þ: a
a
a
Thus we can map the Euler parametrisation to the real and imaginary part parametrisation, enabling us to make inferences on either scale.
Table 6 Parameter estimates and standard errors (in parenthesis) for the models fitted to the ten toluene data sets. There are multiple experiments performed on two of the days Model parameters Data set 10/04/02 19/07/02 19/07/02 19/07/02 20/07/02 20/07/02 20/07/02 20/07/02 20/07/02 25/06/02
No.
ra
1 2 3 1 2 3 4 5
0.645474 0.349865 0.791656 0.261566 0.508602 0.750223 0.461296 0.803142 0.770863 0.682656
(0.022578) (0.005985) (0.005239) (0.003881) (0.002075) (0.003875) (0.005001) (0.005612) (0.006204) (0.007085)
rb
rc
0.137228 (0.056064) 0.087574 (0.014975) 0.227225 (0.011799) 0.069407 (0.016164) 0.084855 (0.013263) 0.154748 (0.021957) 0.116854 (0.021051) 0.139573 (0.010717) 0.187218 (0.027238) 0.191041 (0.014654)
0.417116 (0.032536) 0.274319 (0.007299) 0.705866 (0.00538) 0.211004 (0.004794) 0.177936 (0.00495) 0.418858 (0.005459) 0.397754 (0.004840) 0.437813 (0.004841) 0.443154 (0.009007) 0.417641 (0.011540)
/b
/c 0.0870644 (0.162950) 0.459937 (0.127292)
1.65588 (0.203921) 0.980186 (0.199370) 0.839556 (0.219400) 1.432762 (0.197714)
0.683764 (0.122716)
1.224589 (0.182336) 0.423377 (0.172441)
52
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
constructed from the profile approach are given in Table 4. The bootstrap results given in Table 5 are based on 1000 samples, and the bias for the parameter estimates is not very large for any of the parameters of interest. The intervals are all quite similar in this example of phenylalanine data.
6. Toulene data We examined data from 10 experiments on the toluene/air interface. Each experiment considered two, three, or four output polarisation angles and 14 input polarisation angles for each of the output angles. The presence of fluoresence in
Fig. 4. Residual plots for the models fitted to the toluene data.
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
toluene can interfere with the SHG signal leading to distortion of the SHG results. The SHG technique is sensitive to organic impurities present at the air/liquid interface and these impurities have an impact on observed intensities. Nonlinear models were fitted to these data sets and, where possible, the model was simplified by excluding any redundant parameters. Specifically, the parameter with the largest nonsignificant p-value (using approximate t-tests on the parameter estimates) was dropped from the model, the model refitted and the process repeated until all the parameters remaining in the model were significant. The final model for each experiment is reported in Table 6. The data are on different scales, so the parameter estimates vary between the experiments: we are interested in comparing ratios of the parameters A, B, and C from different experiments. Residual (see Fig. 4) and normal probability plots for the ten experiments showed systematic patterns for the majority of the models. To improve the fit of the model to the data we considered shifting and scaling the intensities. However, this did not improve the diagnostics in general. The residual plots suggested some potential outliers (or influential points) in the data sets, and we considered removing these and refitting the nonlinear models. In some cases the fitted model was worse, because we had removed points from informative regions of the intensity curves. The fit of the model was deemed adequate in four of the experiments so these four models were used to illustrate the procedure for comparing experiments. The coefficient with the largest magnitude was A, so the other two coefficients were compared to A in the ratios. That is, the two ratios of interest for the four experiments are B/A and C/A, and the two components of these ratios rˆR and /ˆR are given in Table 7 with confidence intervals for both components using the delta method and a Fieller confidence interval for the ratio of magnitudes. The Fieller confidence intervals are again very similar to those calculated by the delta method. We are interested in testing whether the ratios of coefficients are the same across the experiments. A likelihood
53
ratio test can be constructed to test this hypothesis. Under the null hypothesis there are four separate scale parameters and four common parameters making eight model parameters, so we are testing H0: f ðQÞ
4 Y
f ðrAs ; rBA ; rCA ; /B ; /C Þ;
s¼1
where rBA is the ratio of B and A common to all four experiments, against the alternative of a separate model for each data set with 4 5 = 20 parameters HA: f ðQÞ
4 Y
f ðrAs ; rBAs ; rCAs ; /Bs ; /Cs Þ;
s¼1
where rBAs, s = 1,. . .,4 are the ratios of B to A for the four experiments. The likelihood ratio statistic for this test where the four experiments have n1, n2, n3, and n4 observations respectively is ! ðrˆ 21 Þn1 =2 ðrˆ 22 Þn2 =2 ðrˆ 23 Þn3 =2 ðrˆ 24 Þn4 =2 2ln ðrˆ 2 ÞN =2 where N = n1 + n2 + n3 + n4 is the total number of observations over the four experiments, rˆ 2 is estimated from the model under the null hypothesis, and rˆ 2k, k = 1,. . .,4 are estimates of rˆ 2 for the models fitted to each experiment separately. The test statistic is compared to the v2 distribution with 20 8 = 12 degrees of freedom. One of the four experiments detailed in Table 7 has a structural zero for one of the phase angles so the full model has 19 rather than 20 parameters. The test statistic is referred to a v2 distribution with 11 degrees of freedom. The likelihood ratio test statistic is 234.56 which suggests that the common model is not appropriate for these four data sets. There is no theoretical reason for the differences between the different sets of toluene data. The production of reliable standard errors and confidence intervals for ratios of parameters has highlighted the existence of these differences and
Table 7 Ratios of parameters for the toluene models. Confidence intervals are given for the ratios of the magnitudes using both the delta and Fieller method. Some of the phase angles (/ˆ R) are zero because the fitted models were simplified when the phase angles were considered negligible Data set
19/07/2002 Delta Fieller 20/07/2002 Delta Fieller 20/07/2002 Delta Fieller 25/06/2002 Delta Fieller
No.
3 C.I. C.I. 2 C.I. C.I. 3 C.I. C.I. C.I. C.I.
B/A
C/A
rˆR
/Rˆ
rˆR
/Rˆ
0.265353 (0.1436, 0.3871) (0.1440, 0.3875) 0.206270 (0.1492, 0.2634) (0.1491, 0.2633) 0.253317 (0.1646, 0.3420) (0.1644, 0.3418) 0.279850 (0.2363, 0.3234) (0.2365, 0.3237)
1.655880 (1.2562, 2.0556)
0.806696 (0.7641, 0.8493) (0.7647, 0.8500) 0.558311 (0.5426, 0.5740) (0.5427, 0.5741) 0.862253 (0.8347, 0.8898) (0.8351, 0.8902) 0.611789 (0.5770, 0.6466) (0.5772, 0.6468)
0.683764 (0.4432, 0.9243)
0.839556 (0.4095, 1.2696) 1.432762 (1.0452, 1.8203) 0
0
–
0.423377 (0.0854, 0.7614)
54
A.H. Welsh et al. / Chemometrics and Intelligent Laboratory Systems 75 (2005) 45–54
enabled experimenters to address them more rigorously. In particular, the findings presented in this section have led to refinements of the experimental technique which have increased the reproducibility of the results. These experiments will be reported elsewhere. If the differences persist after substantial effort at refinement, then we need to consider the appropriateness of the model and the theory from which it is derived. Whatever the outcome, the analysis presented here will have been an important step towards scientific advance and understanding. 7. Conclusions We have explored statistical issues which arise in fitting a phenomenological model to data from second harmonic generation experiments and then using the model to make inferences. The key issues in model fitting involve the choice of parametrisation for the complex parameters in the model, the method used to fit the model and the use of diagnostics to evaluate the quality of model fit. The key issues in inference are the calculation of standard errors, setting confidence intervals for ratios of parameters and comparing models fitted to separate experiments. The Euler parametrisation has been found to be the better representation of complex numbers. The curvature plots to assess the linear approximation are better, and it is easier to form ratios of the parameters A, B, and C and make inferences on these ratios. Fitting the model simultaneously to all the data is better than the sequential approach, because it provides a better fit (less bias) and is more efficient (standard errors are smaller than in the sequential method). Standard residual and normal probability plots supplemented by profile plots to explore curvature in the model provide useful insights into the quality of model fit. We have derived approximate standard errors for the parameter estimates and various useful transformations of the parameter estimates. The approximate standard errors of the ratios of parameters can be used to set approximate confidence intervals provided the denominator estimate is not too small. These intervals can be compared to Fieller intervals which do not use the same kinds of approximation
but may be the whole real line, a proper interval or the empty set. Substantial differences between the two types of intervals suggest that the approximation behind the simple approximate standard errors is poor. We also considered alternative profile likelihood and bootstrap approaches. All these methods performed very similarly in the examples we considered (supporting the use of the simplest method in these examples) but there is useful information when they differ. Finally, we developed and implemented a likelihood ratio test for comparing models fitted to separate experiments.
References [1] Y.R. Shen, Annual Review of Physical Chemistry 40 (1989) 327 – 350. [2] K.B. Eisenthal, Annual Review of Physical Chemistry 43 (1992) 627 – 661. [3] K.B. Eisenthal, Accounts of Chemical Research 26 (12) (1993) 636 – 643. [4] Y.R. Shen, Surface Science 300 (1 – 3) (1994) 551 – 562. [5] Y.R. Shen, Solid State Communications 102 (2 – 3) (1997) 221 – 229. [6] A.J. Fordyce, W.J. Bullock, et al., Molecular Physics 99 (8) (2001) 677 – 687. [7] V. Mizrahi, E. Sipe, Journal of the Optical Society of America. B, Optical Physics 5 (3) (1988) 660 – 667. [8] A.J. Bell, J.G. Frey, et al., Journal of the Chemical Society. Faraday Transactions 88 (14) (1992) 2027 – 2030. [9] M.J. Crawford, S. Haslam, et al., Chemical Physics Letters 230 (3) (1994) 260 – 264. [10] M.J. Crawford, J.G Frey, et al., Journal of the Chemical Society. Faraday Transactions 92 (8) (1996) 1369 – 1373. [11] S. Haslam, S.G. Croucher, et al., Physical Chemistry Chemical Physics 2 (14) (2000) 3235 – 3245. [12] S. Chatterjee, A.S. Hadi, B. Price, Regression Analysis by Example, 3rd edition, Wiley, New York, 2000, Chap. 4. [13] A.H. Welsh, Aspects of Statistical Inference, Wiley, New York, 1996, pp. 97 – 98. [14] E.C. Fieller, Supplement to The Journal of the Royal Statistical Society 7 (1940) 1 – 64. [15] A.H. Welsh, Aspects of Statistical Inference, Wiley, New York, 1996, pp. 152 – 155. [16] A.C. Davison, D.V. Hinkley, Bootstrap Methods and their Application, Cambridge Univ. Press, 1997, pp. 353 – 358. [17] D.M. Bates, D.G. Watts, Nonlinear Regression Analysis and its Applications, Wiley, New York, 1988, Chap. 7. [18] W.N. Venables, B.D. Ripley, Modern Applied Statistics with S, 4th edition, Springer, New York, 2002, p. 461.