Abstract
A Monte Carlo method is used in addition to functional and individual weighting to overcome multicollinearity problems in multiple linear regression equations applied as quantitativestructurepropertyrelationships, allowing the estimation of correct coefficient confidence intervals. The method was applied to rate constants for the Menschutkin reaction between Et_{3}N and EtI in mono and dialcohols, at 25.00 ºC. 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.
Keywords:
Monte Carlo; multicollinearity; correlation; MLR
Introduction
The search for quantitative structureproperty relationships (QSPR) to interpret and predict solvent effects upon reaction rates often resorts to linear model equations, of the general form,
where k stands for rate constant and x_{n} is a descriptor of a given type of substratesolvent or solventsolvent 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}1 Koppel, I. A.; Palm, V. A. In Advances in Linear Free Energy Relationships; Chapman, N. B.; Shorter, J., eds.; Plenum Press: New York, 1972, ch. 5.^{,}^{2}2 Kamlet, M. J.; Abboud, J.L. M.; Abraham, M. H.; Taft, R. W.; Prog. Phys. Org. Chem. 1981, 35, 485.
Over the years a large number of solvent descriptors have been proposed to describe substratesolvent interactions, some based on macroscopic properties such as dielectric constant, dipole moment and refractive index and others on microscopic properties, usually obtained from probes showing spectroscopic shifts due to different solvent effects upon the signalwave number (UVVis and infrared (IR)) or chemical shift (nuclear magnetic resonance (NMR)).^{3}3 Politzer, P.; Murray, J. S. In Quantitative Treatments of Solute/ Solvent Interactions; Politzer, P.; Murray, J. S., eds.; Elsevier: New York, 1994, ch. 1.
4 Abboud, J.L. M.; Notario, R.; Pure Appl. Chem.
1999, 71, 645.^{}^{5}5 Reichardt, C.; Welton, T.; Solvents and Solvent Effects in Organic Chemistry , 4th ed.; WileyVCH: Weinheim, 2011.
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.^{5}5 Reichardt, C.; Welton, T.; Solvents and Solvent Effects in Organic Chemistry , 4th ed.; WileyVCH: Weinheim, 2011.
The 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.^{6}6 Martins, F.; Santos, S.; Ventura, C.; ElvasLeitão, R.; Santos, L.; Vitorino, S.; Reis, M.; Miranda, V.; Correia, H. F.; AiresdeSousa, J.; Kovalishyn, V.; Latino, D. A. R. S.; Ramos, J.; Viveiros, M.; Eur. J. Med. Chem. 2014, 81, 119. However, 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, σ_{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}7 Milton, J. S.; Arnold, J. C.; Probability and Statistics in the Engineering and Computing Sciences; McGrawHill: New York, 1986.
A 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 leastsquares 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.^{8}8 Livingstone, D.; Data Analysis for Chemists; Oxford University Press: Oxford, 1995.
Monte 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. Unlike parametric statistics, they require no assumptions about the distribution of uncertainties or the collinearity among variables, thus allowing calculation of reliable coefficients and their confidence intervals.^{9}9 Alper, J. S.; Gelb, R. I.; J. Phys. Chem.
1990, 94, 4747.
10 Alper, J. S.; Gelb, R. I.; Talanta
1993, 40, 355.
11 Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T.; Numerical Recipes in Pascal; Cambridge University Press: Cambridge, 1989.^{}^{12}12 Buckland, S. T.; Biometrics
1984, 40, 811.
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 nonrepeating 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}12 Buckland, S. T.; Biometrics 1984, 40, 811.
Alternatively, using the method proposed by Alper and Gelb,^{9}9 Alper, J. S.; Gelb, R. I.; J. Phys. Chem. 1990, 94, 4747.^{,}^{10}10 Alper, J. S.; Gelb, R. I.; Talanta 1993, 40, 355. 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 (σ_{i}). In this instance, it is a wellknown (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 σ_{i}^{2}, the variance of each set of replicate measurements, or the estimate of a single measurement.^{13}13 de Levie, R.; J. Chem. Educ. 1986, 63, 10.^{,}^{14}14 Bevington, P. R.; Data Reduction and Error Analysis for the Physical Sciences; McGrawHill: New York, 1969.
In 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.
Much more uncommon, in the same context, is the use of global weighted regression to correct for effects due to transformations on the dependent variable, although this regression tool has been clearly shown to contribute towards correcting these same effects, mainly when significant compression/expansion of uncertainties is present as is the case with logarithms.^{14}14 Bevington, P. R.; Data Reduction and Error Analysis for the Physical Sciences; McGrawHill: New York, 1969.
15 Kalantar, A. H.; J. Phys. Chem.
1986, 90, 6301.
16 Shatkay, A.; Azor, M.; Anal. Chim. Acta
1981, 133, 183.
17 Gonçalves, R. M. C.; Martins, F. E. L.; Leitão, R. A. E.; J. Phys. Chem.
1994, 98, 9540.
18 de Levie, R.; Advanced Excel for Scientific Data Analysis; Oxford University Press: New York, 2004
19 Tellinghuisen, J.; J. Phys. Chem. A
2000, 104, 11829.
20 Sung, D. D.; Kim, J. Y.; Lee, I.; Chung, S. S.; Park, K. H.; Chem. Phys. Lett.
2004, 392, 378.^{}^{21}21 Bolster, C. H.; Tellinghuisen, J.; Soil Sci. Soc. Am. J.
2010, 74, 760.
For a logarithmic transformation such as Y_{i} = ln (k_{i}), the original assumption is, again, that k_{i} values are normally distributed, but we are actually minimizing the sum of squares of the deviations in ln (k_{i}), and these do not follow a Gaussian distribution, since there is a compression of the data for low k_{i} values and an expansion for the higher ones.^{16}16 Shatkay, A.; Azor, M.; Anal. Chim. Acta
1981, 133, 183.
17 Gonçalves, R. M. C.; Martins, F. E. L.; Leitão, R. A. E.; J. Phys. Chem.
1994, 98, 9540.
18 de Levie, R.; Advanced Excel for Scientific Data Analysis; Oxford University Press: New York, 2004^{}^{19}19 Tellinghuisen, J.; J. Phys. Chem. A
2000, 104, 11829.
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}18 de Levie, R.; Advanced Excel for Scientific Data Analysis; Oxford University Press: New York, 2004
The global weighting factor (w'_{i}) dictated by the mathematical transformation of the data, is given by,
In the present case, the transformation of the rate constant (k_{i}) into ln (k_{i}) yields,
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,
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 dialcohols 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}22 Calado, A. R. T.; Pinheiro, L. M. V.; Albuquerque, L. M. P. C.; Gonçalves, R. M. C.; Rosés, M.; Ràfols, C.; Bosch, E.; Collect. Czech. Chem. Commun. 1994, 59, 898. In this work, the authors applied two suitable models to their experimental data set, one predominantly solvatochromic, the TaftAbboudKamletAbraham equation (TAKA),^{23}23 Abraham, M. H.; Doherty, R. M.; Kamlet, M. J.; Harris, J. M.; Taft, R. W.; J. Chem. Soc., Perkin Trans. 2 1987, 913. and one based prevalently on macroscopic descriptors, the GonçalvesAlbuquerqueSimões equation (GAS).^{24}24 Gonçalves, R. M. C.; Simões, A. M. N.; Albuquerque, L. M. P. C.; J. Chem. Soc., Perkin Trans. 2 1990, 1379.
Rate constants (k) for the reaction of Et_{3}N with EtI and select solvent properties at 25.00 ºC^{19}
For the TAKA equation we have
where the solvatochromic descriptors are π^{*}, a measure of the solvent's dipolarity/polarizability; α and β 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(ε) is the Kirkwood function of the relative permittivity; e relates to dipolarity; g(n_{D}) is a function of the refractive index; n_{D} measures the polarizability effect; E_{T}^{N} is the normalized Dimroth and Reichardt parameter and measures a blend of dipolarity and Lewis acidity effects; and C has the same meaning as above.
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}9 Alper, J. S.; Gelb, R. I.; J. Phys. Chem. 1990, 94, 4747.^{,}^{10}10 Alper, J. S.; Gelb, R. I.; Talanta 1993, 40, 355. therefore, since we chose to use 68% confidence intervals (ca. 1σ), 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 σ = 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}17 Gonçalves, R. M. C.; Martins, F. E. L.; Leitão, R. A. E.; J. Phys. Chem. 1994, 98, 9540. 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 preassumption of a given uncertainty level.
Goodnessoffit has been analyzed in terms of the standard deviation of the fit (σ_{fit}), defined in the usual manner,
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.
Correlation between each variable in equation 8 and linear combinations of two of the remaining variables
Correlation between each variable in equation 9 and linear combinations of two of the remaining variables
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 π^{*} is significant, i.e.,
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 σ_{fit}.
Table 7 depicts the results for the GAS model, for which the best equation is:
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 Ftest 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.
The effect of weighting, using the MC approach, on the original data set and on a subset excluding water for equation 12
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 perceived through σ_{fit} values. Furthermore, global weighting also increases the internal consistency of the tested data set: if one uses nonweighted 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 chisquare 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).^{12}12 Buckland, S. T.; Biometrics 1984, 40, 811. This is a common situation in this type of context and is due probably to the empirical or semiempirical 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.
Supplementary Information
Supplementary information (Pascal code for MC confidence interval estimation) is available free of charge at http://jbcs.sbq.org.br as PDF file.
https://minio.scielo.br/documentstore/16784790/dSW8kLcvCHSyKz4DdP3NKXw/1a97c618653ae3d0dc361ea9f60a1ec1ba08f774.pdfAcknowledgments
Financial support by Fundação para a Ciência e Tecnologia, Portugal, under Projects PEST OE/QUI/UI10612/2013 and PEST UID/MULTI/00612/2013 is greatly appreciated.
References

^{1}Koppel, I. A.; Palm, V. A. In Advances in Linear Free Energy Relationships; Chapman, N. B.; Shorter, J., eds.; Plenum Press: New York, 1972, ch. 5.

^{2}Kamlet, M. J.; Abboud, J.L. M.; Abraham, M. H.; Taft, R. W.; Prog. Phys. Org. Chem. 1981, 35, 485.

^{3}Politzer, P.; Murray, J. S. In Quantitative Treatments of Solute/ Solvent Interactions; Politzer, P.; Murray, J. S., eds.; Elsevier: New York, 1994, ch. 1.

^{4}Abboud, J.L. M.; Notario, R.; Pure Appl. Chem. 1999, 71, 645.

^{5}Reichardt, C.; Welton, T.; Solvents and Solvent Effects in Organic Chemistry , 4^{th} ed.; WileyVCH: Weinheim, 2011.

^{6}Martins, F.; Santos, S.; Ventura, C.; ElvasLeitão, R.; Santos, L.; Vitorino, S.; Reis, M.; Miranda, V.; Correia, H. F.; AiresdeSousa, J.; Kovalishyn, V.; Latino, D. A. R. S.; Ramos, J.; Viveiros, M.; Eur. J. Med. Chem. 2014, 81, 119.

^{7}Milton, J. S.; Arnold, J. C.; Probability and Statistics in the Engineering and Computing Sciences; McGrawHill: New York, 1986.

^{8}Livingstone, D.; Data Analysis for Chemists; Oxford University Press: Oxford, 1995.

^{9}Alper, J. S.; Gelb, R. I.; J. Phys. Chem. 1990, 94, 4747.

^{10}Alper, J. S.; Gelb, R. I.; Talanta 1993, 40, 355.

^{11}Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T.; Numerical Recipes in Pascal; Cambridge University Press: Cambridge, 1989.

^{12}Buckland, S. T.; Biometrics 1984, 40, 811.

^{13}de Levie, R.; J. Chem. Educ. 1986, 63, 10.

^{14}Bevington, P. R.; Data Reduction and Error Analysis for the Physical Sciences; McGrawHill: New York, 1969.

^{15}Kalantar, A. H.; J. Phys. Chem. 1986, 90, 6301.

^{16}Shatkay, A.; Azor, M.; Anal. Chim. Acta 1981, 133, 183.

^{17}Gonçalves, R. M. C.; Martins, F. E. L.; Leitão, R. A. E.; J. Phys. Chem. 1994, 98, 9540.

^{18}de Levie, R.; Advanced Excel for Scientific Data Analysis; Oxford University Press: New York, 2004

^{19}Tellinghuisen, J.; J. Phys. Chem. A 2000, 104, 11829.

^{20}Sung, D. D.; Kim, J. Y.; Lee, I.; Chung, S. S.; Park, K. H.; Chem. Phys. Lett. 2004, 392, 378.

^{21}Bolster, C. H.; Tellinghuisen, J.; Soil Sci. Soc. Am. J. 2010, 74, 760.

^{22}Calado, A. R. T.; Pinheiro, L. M. V.; Albuquerque, L. M. P. C.; Gonçalves, R. M. C.; Rosés, M.; Ràfols, C.; Bosch, E.; Collect. Czech. Chem. Commun. 1994, 59, 898.

^{23}Abraham, M. H.; Doherty, R. M.; Kamlet, M. J.; Harris, J. M.; Taft, R. W.; J. Chem. Soc., Perkin Trans. 2 1987, 913.

^{24}Gonçalves, R. M. C.; Simões, A. M. N.; Albuquerque, L. M. P. C.; J. Chem. Soc., Perkin Trans. 2 1990, 1379.
Publication Dates

Publication in this collection
Nov 2016
History

Received
16 Dec 2015 
Accepted
31 Mar 2016