*** For details on how to be removed from this list visit the *** *** CCP4 home page http://www.ccp4.ac.uk ***
On 9/5/05 4:52 AM, "George M. Sheldrick" <[EMAIL PROTECTED]> a écrit : > I am now completely confused by the ML vs. LS discussion and hope that I > am not the only person with this problem! > > To clarify the issue, here is a simple clear-cut example of considerable > interest to small-molecule crystallographers: what is the best way to > determine the 'absolute configuration' using weak anomalous scattering > (e.g. for an organic compound containing only C, H, N and O). > > The conventional wisdom is that a so-called Flack parameter x is refined > as one of the LEAST-SQUARES parameters, where Iobs(hkl) is fitted by > Icalc(hkl) + x.Icalc(-h-k-l) with weights based on sigma(Iobs). Depending on the model assumed, this LS method can be interpreted as ML. For instance, let's assume that the 'real' intensity I(hkl) is subject to random Gaussian error. Then, one model is Iobs(hkl) = (1-x)I(hkl) + x.I(-h-k-l) + e.sigmaI(hkl) where e is a standard Gaussian random variable (with variance 1) and sigmaI(hkl) is the std deviation for each hkl. Here, x is a variable dependent on the current molecular model. Then, a LS fit of Iobs(hkl) to (1-x)I(hkl) + xI(-h-k-l) weighted by the inverse of sigmaI(hkl) is equivalent to the ML solution given by the above model. In this case there is a relatively simple analytical solution for x, if we are estimating I(hkl), I(-h-k-l), and sigmaI(hkl) by Iobs(hkl), Iobs(-h-k-l), and sigmaIobs(hkl) respectively. The LL(x) is given by -0.5[(Iobs(hkl) - (1-x)I(hkl) - x.I(-h-k-l))/sigmaI(hkl)]^2 up to an arbitrary constant. Note that this model has at least one theoretical problem, in that real, physically measurable intensities are always positive, and this model will give negative Iobs with high probability when sigmaI is large enough. IOW, intensities cannot be distributed by a Gaussian, except approximately. > Complex > scattering factors are employed to calculate the intensity Icalc. This > is similar to the way twins are refined. > > When x is 0 (with a 'small' esd from the full-matrix LS) the absolute > configuration is correct, when it is 1 (with a 'small' esd) it is wrong. > If there are no errors in the model then only these two values of x are > possible (racemic twinning, which would permit x to take any value in > the range 0 to 1, is rare but not unknown; however in some case it can > be ruled out by chemical evidence). Clearly the estimated esd of x is > just as important as the actual value of x. > > The following 'ML' approach was proposed in discussions at the Computing > School in Siena (contributions from Rob Hooft, Simon Parsons and David > Watkin are acknowledged, but I'm sure that there were others). We define > the quantity Qobs = [Iobs(hkl)-Iobs(-h-k-l)]/[Iobs(hkl)+Iobs(-h-k-l)] > and a corresponding Qcalc. This should cancel some systematic errors. > The esd of Qobs can be calculated by standard methods from the known > esds of Iobs(hkl) and Iobs(-h-k-l). Then we estimate the log likelihood > (LL) of a particular value of x by -0.5[(Qobs-Qcalc)/esd(Qobs)]^2 summed > over all reflections. This model is both LS and ML in the same way as the above model. The difference is that the implied model is a little different: here you assume that there is a 'real' Q(hkl), and that Qobs(hkl) is subject to random Gaussian errors as Qobs(hkl) = Q(hkl) + e.sigmaQ(hkl) > So my question is: should I estimate the value of x and its esd by > integrating over -infinity < x < +infinity (in practice this range can > be truncated to say -5 to +5) or should I integrate over 0 =< x =< 1 or > should I use only the values of the LL at x = 0 and x = 1? My opinion is to use 0 =< x =< 1 if you expect a reasonable possibility of racemic twinning; otherwise, constrain x to 0 and 1. > Preliminary tests suggest that these methods give appreciably lower esds for x > than the LS approach. As I said, the 'LS' approach above (the first model discussed) can also be thought of as ML. So the only real difference is in the assumed model, whether random Gaussian errors are applied to I(hkl) vs Q(hkl). This immediately begs the question as to which model is more appropriate. For one thing, it appears to me that Qobs(hkl) can take negative values in reality, so Gaussian errors are probably more justifiable for this model, and this may be the main reason it performs better. Note that, since both models (I and Q) have the same number of parameters, you can easily compare the LL from both models and see which one is largest (using, in effect, a trivial Akaike Info Criterion as a basis for model selection). If the the LL ratio is greater than about 3, then that is strong evidence in favor of the model with larger LL. Morris will probably kick me for suggesting it, but if you *must* have a P-value, there is an easy frequentist hypothesis test you can do. First, calculate the individual LL ratios for each hkl for the I vs Q models, and also calculate the overall LL ratio. The variance of the overall LL is equivalent (asymptotically) to the square root of the number of individual hkls times the variance of the individual LL ratios. If the models are equivalent, the overall LLR should be 0. So you can use the relations above to do a simple z-test for a significant difference from 0. A Bayesian in principle, a likelihoodist in practice, a frequentist in public. :) Cheers, Douglas > If I only consider the LL for x=0 and x=1 I have a > simple way of calculating the probability that x is 0 rather than 1 > (i.e. the probability that the absolute configuration of the model is > correct). Note that because of correlations between x and the atomic > positional parameters in polar space groups it might be necessary to > refine the structure to convergence for x=0 and x=1 to get the Icalc > values for this. > > I would appreciate clarification from the ML experts, with a view to > putting it into my programs (which are still used by some small molecule > crystallographers). > > George
