Fighting Collinearity in QSPR Equations for Solution Kinetics with the Monte Carlo Method and Total Weighting

A Monte Carlo method is used in addition to functional and individual weighting to overcome multicollinearity problems in multiple linear regression equations applied as quantitative-structureproperty-relationships, allowing the estimation of correct coefficient confidence intervals. The method was applied to rate constants for the Menschutkin reaction between Et3N and EtI in monoand di-alcohols, at 25.00 oC. Results show that the use of our methodology produces a significant improvement upon confidence interval estimates regardless of the level of collinearity present. Addition of weighting shows additional advantages, increasing the overall consistency of the regression process.


Introduction
The search for quantitative structure-property relationships (QSPR) to interpret and predict solvent effects upon reaction rates often resorts to linear model equations, of the general form, ln k = f(x 1 , x 2 ,…, x n ) (1)   where k stands for rate constant and x n is a descriptor of a given type of substrate-solvent or solvent-solvent interaction.
Early models included just one or two descriptors but it is now well established that, in general, four descriptors are needed to account for the main physicochemical features underlying the reaction process.These descriptors include Lewis acidity and basicity, dipolarity/polarizability and also a cavity term related to the work needed to form a cavity to accommodate the substrate. 1,2][5] Subsequently, several multiparametric model equations corresponding to different combinations of descriptors have been used, leading to a greater understanding of the type of interactions present and their effect upon the rate of a given reaction, thus providing useful mechanistic clues and, most of all, generating sound equations to accurately predict rate constants for similar reactions. 5he most useful method to obtain this information has been found to be multiple linear regression (MLR).Ideally, from a mathematical point of view, any descriptors in a given MLR model equation should be orthogonal to each one of the others and, desirably, also to their linear combinations threshold values being, respectively, r 2 > 0.5 and R 2 > 0.8. 6However, more often than not, they present significant degrees of collinearity, especially when only a closely related set of solvents is used, e.g., a set of monoalcohols, and/or small variabilities are observed for the chosen descriptors, within the data set.
In general, this fact does not affect significantly regression coefficients, except in cases where descriptors present high correlation coefficients.Even when collinearity does not affect regression coefficients, consequences of collinearity become apparent through the often inflated magnitude of coefficient uncertainties, s i , whose values are calculated from the diagonal values of the variancecovariance matrix.From these, confidence intervals for the coefficients (CI i ), can subsequently be calculated through the use of t-values for the chosen confidence level (e.g., 68, 95, or 99%), according to, CIs are particularly important, for interpretative and predictive purposes, especially if one aims to interpolate reliable estimates of predicted k values. 7 number of methods proposed over time (principle component analysis, factor analysis, etc.) although successful in overcoming this problem, are more difficult to apply and change the descriptors' structure in their effort to orthogonalize descriptors' axes, leading to solutions that frequently bear unclear physical meaning and make interpretation of results difficult if not impossible.Other techniques such as ridge regression produce unbiased estimates of parameters, but their expected values are not at all equal to the true values.Generally, they tend to be underestimated (sometimes grossly) even if the variance of these estimates can be so much lower than that of the least-squares estimator and, of course, the total expected mean squared error is also less, which makes it (in a certain sense) a "better" estimator for some but surely not all intended purposes. 8onte Carlo (MC) methods, on the other hand, are rather less known as an efficient method for dealing with the type of problem here addressed, although they are widely used in computer simulations.][11][12] The Monte Carlo method assumes that the data set is a sample of all possible sets, randomly drawn by the experimental method within experimental uncertainty.Using this procedure one can simulate as many synthetic data sets as one needs (n), drawn from a particular model, just by introducing the appropriate random noise (δ i ): Random numbers are generated using any algorithm that can ensure long non-repeating sequences with the appropriate standard deviation and zero mean value.Usually, since it is assumed that coefficient uncertainties follow a Gaussian distribution, the generated numbers are normally distributed in order to later allow the calculation of confidence intervals for each equation parameter.Each synthetic data collection is then subjected to a regression analysis leading to n "synthetic" parameter sets (a 1 , a 2 ,..., a n ), (b 1 , b 2 ,..., b n ), etc. 12 Alternatively, using the method proposed by Alper and Gelb, 9,10 confidence intervals (given from the adequate number of computer runs) can be easily and accurately calculated without the usual assumption of normality for uncertainty distribution by the simple exclusion of the appropriate number of high and low solution values from the generated sets by ordering the calculated values and excluding the (100 -x)% / 2 top and bottom ones, thus obtaining x% CIs.The minimum number of sets to be generated for each confidence level has also been established by these authors.
Finally, another usual MLR assumption is that uncertainties for the whole experimental data set are equal.In reality experimental results may have very different standard deviations (s i ).In this instance, it is a well-known (but scarcely used) practice to use individual correcting weights (w i ).This type of individual weighting accounts for the accuracy of each value and, in most situations, it is considered to be inversely proportional to s i 2 , the variance of each set of replicate measurements, or the estimate of a single measurement. 13,14n the case of kinetic measurements, the accuracy of k i values may vary significantly due to experimental errors and especially if values from different authors and/or different measurement techniques have been used.
][18][19] When experimental k i values are converted to Y i (e.g., ln k, 1 / k), assuming the usual hypothesis that ΔY i and Δk i are relatively small, we can write, following de Levie, 18 (4) The global weighting factor (w' i ) dictated by the mathematical transformation of the data, is given by, (5)     Fighting Collinearity in QSPR Equations for Solution Kinetics with the Monte Carlo Method and Total Weighting J. Braz.Chem.Soc.2072   In the present case, the transformation of the rate constant (k i ) into ln (k i ) yields, (6)   Both types of weighing factors should be used together: their combination (W i ), the total weight, is achieved by multiplying the two weighing factors.In the present case, this produces the expression, (7)   In this paper, we show the advantages of the combined use of these methodology changes in multiparametric QSPR applied to kinetics, through its application to a data set of 20 rate constants for the Menschutkin reaction of Et 3 N with EtI in mono-and di-alcohols at 25.00 ºC (Table 1) that were previously correlated with two sets of four solvent descriptors from two multiparametric QSPR models by Calado et al. 22 In this work, the authors applied two suitable models to their experimental data set, one predominantly solvatochromic, the Taft-Abboud-Kamlet-Abraham equation (TAKA), 23 and one based prevalently on macroscopic descriptors, the Gonçalves-Albuquerque-Simões equation (GAS). 24or the TAKA equation we have where the solvatochromic descriptors are p * , a measure of the solvent's dipolarity/polarizability; a and b which are, respectively, measures of the solvent's ability to donate (Lewis acidity) or accept (Lewis basicity) a hydrogen bond from a given substrate; and C which is the cohesive energy density, a macroscopic descriptor intended to measure the solvent's contribution to the formation of a cavity to harbor the substrate's molecule.For the GAS equation we have In this equation, f(e) is the Kirkwood function of

Methodology
The procedure used was as follows: (i) classical regression was performed over the experimental data using a correlation equation to determine a set of coefficients; (ii) each k i value was scattered by adding a random number, δ i (equation 3), to construct a new "experimental" mathematical solution; (iii) the resulting MC data set was analyzed through classical regression and a new set of coefficients was obtained; (iv) this sequence was repeated n times, where n is the necessary total number of synthetic data sets; and (v) the CI for each coefficient was obtained as described above by ordering the coefficient values and excluding the (100 -68)% / 2 top and bottom ones.
The computational conditions were chosen using Alper and Gelb's criteria; 9,10 therefore, since we chose to use 68% confidence intervals (ca.1s), the number of MC simulations, n, equaled 200.Normally distributed random uncertainty was introduced through a routine written before.The number set had µ = 0 and s = 1.5 × 10 -4 .This latter value for the standard deviation can be considered as a scattering factor.It is usually chosen from information on the magnitude of the uncertainties affecting the dependent variable.In the present study, a different approach, previously developed by us, was used: 17 using the TAKA equation in each monodescriptor reduced form, where no collinearity problems can be present, the standard deviation of the random uncertainty numbers set was adjusted by repeating the MC procedure with different scattering factors until the averaged coefficients' confidence intervals were identical to those calculated through classical regression.The same scattering factor was then used for equations 7 and 8.This procedure has avoided the need for previous knowledge or pre-assumption of a given uncertainty level.
Goodness-of-fit has been analyzed in terms of the standard deviation of the fit (s fit ), defined in the usual manner, (10)   where n is the number of data points, p is the number of equation parameters and W i is the total weighting which is equal to unity if none of the weighting types are used.

Results and Discussion
The degree of collinearity among equation descriptors has been determined in terms of the determination coefficient (r 2 ) (Tables 2 and 3), showing that several descriptors are highly correlated in both models.
The correlation of each descriptor with linear combinations of pairs of the remaining descriptors (Tables 4  and 5) is also above the limit (R 2 > 0.8) in several cases.
Results above are not surprising since the observed multicollinearity is mostly a result of studying a specific family of solvents.Nevertheless, the best equations' subsets were the same for the classical and Monte Carlo methods.
Regarding the best equation for the TAKA model, the stepwise regression method indicates that only p * is significant, i.e., ln k = a 0 + a 1 p* (11)   and therefore this model will be used only to show the strict equivalence of results for the two approaches, as can be seen in Table 6 in which the comparison of results between classical regression and the MC method for equation 11   shows only very small differences in the confidence intervals and no difference in s fit .
Table 7 depicts the results for the GAS model, for which the best equation is: (12)   Results clearly indicate that, although the coefficients obtained in both cases are similar, the use of the Monte Carlo method leads to significantly less conservative estimates of the coefficients' confidence intervals, with relative differences ranging from -26 to -64%.The MC fit itself shows a somewhat larger standard deviation but with no statistical significance (the F-test on variances confirmed that they are equal up to a significance level of 98%).
The effect of weighting in addition to the use of the MC method is shown in Table 8, along with the same effects on a subset that excludes water.
On one hand, it is clear from the 1 st and 3 rd rows that the use of weights improves regression results, as can be through s fit values.Furthermore, global weighting also increases the internal consistency of the tested data set: if one uses non-weighted regression, the comparison of results for the original data set (20 points) and a subset excluding water shows a change on the  sign associated with a 3 .However, a similar comparison for the weighted regression shows an alteration on the magnitude of a 3 but with the sign remaining unchanged.This apparently small difference is rather significant in terms of both interpretation and prediction.It is also easily seen that the uncertainties affecting each k value are generally very small when compared with the residuals of each fit, causing the chi-square of each fit to become larger than the number of degrees of freedom in all cases (a "rule of thumb" indicates that they should be close). 12This is a common situation in this type of context and is due probably to the empirical or semi-empirical nature of solvent descriptors and to the formal procedure of considering the mere additivity of all included effects.

Conclusions
From the above analysis, we can conclude that the proposed MC method, especially when combined with global weighting, allows a significant improvement on the calculation of uncertainties on MLR coefficients, and increases the overall consistency of the regression process otherwise affected by the presence of multicollinearity.The interpretation and/or prediction process becomes, therefore, more reliable.

r 2 :
determination coefficient; f(e): Kirkwood function of the relative permittivity; g(n D ): function of the refractive index; E T N : normalized Dimroth and Reichardt parameter; C: cohesive energy density.

Table 1 .
19te constants (k) for the reaction of Et 3 N with EtI and select solvent properties at 25.00 ºC19 k: rate constant; e: relative permittivity; n D : refractive index; p*: TAKA dipolarity/polarizability parameter; a: TAKA Lewis acidity parameter; b: TAKA Lewis basicity parameter; C: cohesive energy density.

Table 3 .
Correlation among variables in equation 9 r 2

Table 2 .
Correlation among variables in equation 8

Table 4 .
Correlation between each variable in equation 8 and linear combinations of two of the remaining variables

Table 5 .
Correlation between each variable in equation 9 and linear combinations of two of the remaining variables

Table 6 .
Comparison between classical and Monte Carlo methods for equation 11

Table 7 .
Comparison between classical and Monte Carlo methods for equation 12

Table 8 .
The effect of weighting, using the MC approach, on the original data set and on a subset excluding water for equation 12 n: number of data points; a i ± s(a i ): coefficient confidence interval; s fit : standard deviation of the fit.