Dear Boaz, You are quite correct, 'latter' and 'former' need to be switched in my email. Apologies to CCP4bb for the confusion caused!
Best wishes, George Prof. George M. Sheldrick FRS Dept. Structural Chemistry, University of Goettingen, Tammannstr. 4, D37077 Goettingen, Germany Tel. +49-551-39-3021 or -3068 Fax. +49-551-39-22582 On Sun, 13 Mar 2011, Boaz Shaanan wrote: > Dear George, > > While I agree with you I wonder whether in this statement: > > ..."The practice of quoting R-values both > for all data and for F>4sigma(F) seems to me to be useful. For example > if the latter is much larger than the former, maybe you are including a > lot of weak data"... > > Shouldn't it be: ..."former (i.e. R for all data) is much larger than the > latter (i.e. R for F>4sigma(F)"... ? > Just wondering, although it could be my late night misunderstanding. > > Best regards, > > Boaz > > Boaz Shaanan, Ph.D. > Dept. of Life Sciences > Ben-Gurion University of the Negev > Beer-Sheva 84105 > Israel > > Phone: 972-8-647-2220 Skype: boaz.shaanan > Fax: 972-8-647-2992 or 972-8-646-1710 > > > > ________________________________________ > From: CCP4 bulletin board [CCP4BB@JISCMAIL.AC.UK] On Behalf Of George M. > Sheldrick [gshe...@shelx.uni-ac.gwdg.de] > Sent: Sunday, March 13, 2011 12:11 AM > To: CCP4BB@JISCMAIL.AC.UK > Subject: Re: [ccp4bb] I/sigmaI of >3.0 rule > > Dear James, > > I'm a bit puzzled by your negative R-values and unstable behavior. In > practice, whether we refine against intensity or against |F|, it is > traditional to quote an R-factor (called R1 in small molecule > crystallography) R = Sum||Fo|-|Fc|| / Sum|Fo|. Reflections that have > negative measured intensities are either given F=0 or (e.g. using > TRUNCATE) F is set to a small positive value, both of which avoid having > to take the square root of a negative number which most computers don't > like doing. Then the 'divide by zero' catastropy and negative R-values > cannot happen because Sum|Fo| is always significantly greater than zero, > and in my experience there is no problem in calculating an R-value even > if the data are complete noise. The practice of quoting R-values both > for all data and for F>4sigma(F) seems to me to be useful. For example > if the latter is much larger than the former, maybe you are including a > lot of weak data. Similarly in calculating merging R-values, most > programs replace negative intensities by zero, again avoiding the > problems you describe. > > Best wishes, George > > Prof. George M. Sheldrick FRS > Dept. Structural Chemistry, > University of Goettingen, > Tammannstr. 4, > D37077 Goettingen, Germany > Tel. +49-551-39-3021 or -3068 > Fax. +49-551-39-22582 > > > On Sat, 12 Mar 2011, James Holton wrote: > > > > > The fundamental mathematical problem of using an R statistic on data with > > I/sigma(I) < 3 is that the assumption that the fractional deviates > > (I-<I>)/<I> > > obey a Gaussian distribution breaks down. And when that happens, the R > > calculation itself becomes unstable, giving essentially random R values. > > Therefore, including weak data in R calculations is equivalent to > > calculating > > R with a 3-sigma cutoff, and then adding a random number to the R value. > > Now, > > "random" data is one thing, but if the statistic used to evaluate the data > > quality is itself random, then it is not what I would call "useful". > > > > Since I am not very good at math, I always find myself approaching > > statistics > > by generating long lists of random numbers, manipulating them in some way, > > and > > then graphing the results. For graphing Rmerge vs I/sigma(I), one does find > > that "Bernhard's rule" of Rmerge = 0.8/( I/sigma(I) ) generally applies, but > > only for I/sigma(I) that is >= 3. It gets better with high multiplicity, > > but > > even with m=100, the Rmerge values for the I/sigma(I) < 1 points are all > > over > > the place. This is true even if you average the value of Rmerge over a > > million random number "seeds". In fact, one must do so much averaging, > > that I > > start to worry about the low-order bits of common random number generators. > > I > > have attached images of these "Rmerge vs I/sigma" graphs. The error bars > > reflect the rms deviation from the average of a large number of "Rmerge" > > values (different random number seeds). The "missing" values are actually > > points where the average Rmerge in 600000 trials (m=3) was still negative. > > > > The reason for this "noisy R factor problem" becomes clear if you consider > > the > > limiting case where the "true" intensity is zero, and make a histogram of ( > > I > > - <I> )/<I>. It is not a Gaussian. Rather, it is the Gaussian's evil > > stepsister: the Lorentzian (or "Cauchy distribution"). This distribution > > may > > look a lot like a Gaussian, but it has longer tails, and these tails give it > > the weird statistical property of having an undefined mean value. This is > > counterintuitive! Because you can clearly just look at the histogram and > > see > > that it has a central peak (at zero), but if you generate a million > > Lorentzian-distributed random numbers and take the average value, you will > > not > > get anything close to zero. Try it! You can generate a Lorentzian deviate > > from a uniform deviate like this: tan(pi*(rand()-0.5)), where "rand()" > > makes a > > random number from 0 to 1. > > > > Now, it is not too hard to understand how R could "blow up" when the "true" > > spot intensities are all zero. After all, as <I> approaches zero, the > > ratio ( > > I - <I> ) / <I> approaches a divide-by-zero problem. But what about when > > I/sigma(I) = 1? Or 2? If you look at these histograms, you find that they > > are a "cross" between a Gaussian and a Lorentzian (the so-called "Voigt > > function"), and the histogram does not become truly "Gaussian-looking" until > > I/sigma(I) = 3. At this point, the R factor behaves "Bernhard's rule" quite > > well, even with multiplicities as low as 2 or 3. This was the moment when I > > realized that the "early crystallographers" who first decided to use this > > 3-sigma cutoff, were smarter than I am. > > > > Now, you can make a Voigt function (or even a Lorentzian) look more like a > > Gaussian by doing something called "outlier rejection", but it is hard to > > rationalize why the "outliers" are being rejected. Especially in a > > simulation! Then again, the silly part of all this is all we really want is > > the "middle" of the histogram of ( I - <I> )/<I>. In fact, if you just pick > > the "most common" Rmerge, you would get a much better estimate of the "true > > Rmerge" in a given resolution bin than you would by averaging a hundred > > times > > more data. Such procedures are called "robust estimators" in statistics, > > and > > the "robust estimator" equivalents to the average and the "rms deviation > > from > > the average" are the median and the "median absolute deviation from the > > median". If you make a list of Lorentzian-random numbers as above, and > > compute the median, you will get a value very close to zero, even with > > modest > > multiplicity! And the "median absolute deviation from the median" rapidly > > converges to 1, which matches the "full width at half maximum" of the > > histogram quite nicely. > > > > So, what are the practical implications of this? Perhaps instead of the > > average Rmerge in each bin we should be looking at the median Rmerge? This > > will be the same as the average for the cases where I/sigma(I) > 3, but > > still > > be "well behaved" for the bins that contain only weak data. The last two > > graphs attached compare the average (red) Rmerge vs the median (blue). > > These > > were computed for the same data, and clearly the median Rmerge is a lot more > > stable. The blue error bars are the median absolute deviation from the > > median. True, the median Rmerge curves down a bit as I/sigma approaches > > zero > > at m=3, but at least it doesn't go negative! > > > > Just a suggestion. > > > > -James Holton > > MAD Scientist > > > > On 3/9/2011 8:33 AM, Graeme Winter wrote: > > > Hi James, > > > > > > May I just offer a short counter-argument to your case for not including > > > weak reflections in the merging residuals? > > > > > > Unlike many people I rather like Rmerge, not because it tells you how good > > > the data are, but because it gives you a clue as to how well the unmerged > > > measurements agree with one another. It's already been mentioned on this > > > thread that Rmerge is ~ 0.8 / <I/sigma> which means that the inverse is > > > also > > > true - an Rmerge of 0.8 indicates that the average measurement in the > > > shell > > > has an I/sigma of ~ 1 (presuming there are sufficient multiple > > > measurements > > > - if the multiplicity is < 3 or so this can be nonsense) > > > > > > This does not depend on the error model or the multiplicity. It just talks > > > about the average. Now, if we exclude all measurements with an I/sigma of > > > less than three we have no idea of how strong the reflections in the shell > > > are on average. We're just top-slicing the good reflections and asking how > > > well they agree. Well, with an I/sigma > > > > 3 I would hope they agree rather well if your error model is > > > reasonable. It would suddenly become rare to see an Rmerge > 0.3 in the > > > outer shell. > > > > > > I like Rpim. It tells you how good the average measurement should be > > > provided you have not too much radiation damage. However, without Rmerge I > > > can't get a real handle on how well the measurements agree. > > > > > > Personally, what I would like to see is the full contents of the Scala log > > > file available as graphs along with Rd from xdsstat and some other choice > > > statistics so you can get a relatively complete picture, however I > > > appreciate that this is unrealistic :o) > > > > > > Just my 2c. > > > > > > Cheerio, > > > > > > Graeme > > > > > > On 8 March 2011 20:07, James Holton <jmhol...@lbl.gov > > > <mailto:jmhol...@lbl.gov>> wrote: > > > > > > Although George does not mention anything about data reduction > > > programs, I take from his description that common small-molecule > > > data processing packages (SAINT, others?), have also been > > > modernized to record all data (no I/sigmaI > 2 or 3 cutoff). I > > > agree with him that this is a good thing! And it is also a good > > > thing that small-molecule refinement programs use all data. I > > > just don't think it is a good idea to use all data in R factor > > > calculations. > > > > > > Like Ron, I will probably be dating myself when I say that when I > > > first got into the macromolecular crystallography business, it was > > > still commonplace to use a 2-3 sigma spot intensity cutoff. In > > > fact, this is the reason why the PDB wants to know your > > > "completeness" in the outermost resolution shell (in those days, > > > the outer resolution was defined by where completeness drops to > > > ~80% after the 3 sigma spot cutoff). My experience with this, > > > however, was brief, as the maximum-likelihood revolution was just > > > starting to take hold, and the denzo manual specifically stated > > > that only bad people use sigma cutoffs > -3.0. Nevertheless, like > > > many crystallographers from this era, I have fond memories of the > > > REALLY low R factors you can get by using this arcane and now > > > reviled practice. Rsym values of 1-2% were common. > > > > > > It was only recently that I learned enough about statistics to > > > understand the wisdom of my ancestors and that a 3-sigma cutoff is > > > actually the "right thing to do" if you want to measure a > > > fractional error (like an R factor). That is all I'm saying. > > > > > > -James Holton > > > MAD Scientist > > > > > > > > > On 3/6/2011 2:50 PM, Ronald E Stenkamp wrote: > > > > > > My small molecule experience is old enough (maybe 20 years) > > > that I doubt if it's even close to representing current > > > practices (best or otherwise). Given George's comments, I > > > suspect (and hope) that less-than cutoffs are historical > > > artifacts at this point, kept around in software for making > > > comparisons with older structure determinations. But a bit of > > > scanning of Acta papers and others might be necessary to > > > confirm that. Ron > > > > > > > > > On Sun, 6 Mar 2011, James Holton wrote: > > > > > > > > > Yes, I would classify anything with I/sigmaI < 3 as > > > "weak". And yes, of course it is possible to get "weak" > > > spots from small molecule crystals. After all, there is no > > > spot so "strong" that it cannot be defeated by a > > > sufficient amount of background! I just meant that, > > > relatively speaking, the intensities diffracted from a > > > small molecule crystal are orders of magnitude brighter > > > than those from a macromolecular crystal of the same size, > > > and even the same quality (the 1/Vcell^2 term in Darwin's > > > formula). > > > > > > I find it interesting that you point out the use of a 2 > > > sigma(I) intensity cutoff for small molecule data sets! > > > Is this still common practice? I am not a card-carrying > > > "small molecule crystallographer", so I'm not sure. > > > However, if that is the case, then by definition there are > > > no "weak" intensities in the data set. And this is > > > exactly the kind of data you want for least-squares > > > refinement targets and computing "% error" quality metrics > > > like R factors. For likelihood targets, however, the > > > "weak" data are actually a powerful restraint. > > > > > > -James Holton > > > MAD Scientist > > > > > > On 3/6/2011 11:22 AM, Ronald E Stenkamp wrote: > > > > > > Could you please expand on your statement that > > > "small-molecule data has essentially no weak spots."? > > > The small molecule data sets I've worked with have > > > had large numbers of "unobserved" reflections where I > > > used 2 sigma(I) cutoffs (maybe 15-30% of the > > > reflections). Would you consider those "weak" spots > > > or not? Ron > > > > > > On Sun, 6 Mar 2011, James Holton wrote: > > > > > > I should probably admit that I might be indirectly > > > responsible for the resurgence of this I/sigma > 3 > > > idea, but I never intended this in the way > > > described by the original poster's reviewer! > > > > > > What I have been trying to encourage people to do > > > is calculate R factors using only hkls for which > > > the signal-to-noise ratio is > 3. Not refinement! > > > Refinement should be done against all data. I > > > merely propose that weak data be excluded from > > > R-factor calculations after the > > > refinement/scaling/mergeing/etc. is done. > > > > > > This is because R factors are a metric of the > > > FRACTIONAL error in something (aka a "% > > > difference"), but a "% error" is only meaningful > > > when the thing being measured is not zero. > > > However, in macromolecular crystallography, we > > > tend to measure a lot of "zeroes". There is > > > nothing wrong with measuring zero! An excellent > > > example of this is confirming that a systematic > > > absence is in fact "absent". The "sigma" on the > > > intensity assigned to an absent spot is still a > > > useful quantity, because it reflects how confident > > > you are in the measurement. I.E. a sigma of "10" > > > vs "100" means you are more sure that the > > > intensity is zero. However, there is no "R factor" > > > for systematic absences. How could there be! This > > > is because the definition of "% error" starts to > > > break down as the "true" spot intensity gets > > > weaker, and it becomes completely meaningless when > > > the "true" intensity reaches zero. > > > > > > Historically, I believe the widespread use of R > > > factors came about because small-molecule data has > > > essentially no weak spots. With the exception of > > > absences (which are not used in refinement), spots > > > from "salt crystals" are strong all the way out to > > > edge of the detector, (even out to the "limiting > > > sphere", which is defined by the x-ray > > > wavelength). So, when all the data are strong, a > > > "% error" is an easy-to-calculate quantity that > > > actually describes the "sigma"s of the data very > > > well. That is, sigma(I) of strong spots tends to > > > be dominated by things like beam flicker, spindle > > > stability, shutter accuracy, etc. All these > > > usually add up to ~5% error, and indeed even the > > > Braggs could typically get +/-5% for the intensity > > > of the diffracted rays they were measuring. > > > Things like Rsym were therefore created to check > > > that nothing "funny" happened in the measurement. > > > > > > For similar reasons, the quality of a model > > > refined against all-strong data is described very > > > well by a "% error", and this is why the > > > refinement R factors rapidly became popular. Most > > > people intuitively know what you mean if you say > > > that your model fits the data to "within 5%". In > > > fact, a widely used criterion for the correctness > > > of a "small molecule" structure is that the > > > refinement R factor must be LOWER than Rsym. This > > > is equivalent to saying that your curve (model) > > > fit your data "to within experimental error". > > > Unfortunately, this has never been the case for > > > macromolecular structures! > > > > > > The problem with protein crystals, of course, is > > > that we have lots of "weak" data. And by "weak", > > > I don't mean "bad"! Yes, it is always nicer to > > > have more intense spots, but there is nothing > > > shameful about knowing that certain intensities > > > are actually very close to zero. In fact, from > > > the point of view of the refinement program, isn't > > > describing some high-angle spot as: "zero, plus or > > > minus 10", better than "I have no idea"? Indeed, > > > several works mentioned already as well as the > > > "free lunch algorithm" have demonstrated that > > > these "zero" data can actually be useful, even if > > > it is well beyond the "resolution limit". > > > > > > So, what do we do? I see no reason to abandon R > > > factors, since they have such a long history and > > > give us continuity of criteria going back almost a > > > century. However, I also see no reason to punish > > > ourselves by including lots of zeroes in the > > > denominator. In fact, using weak data in an R > > > factor calculation defeats their best feature. R > > > factors are a very good estimate of the fractional > > > component of the total error, provided they are > > > calculated with strong data only. > > > > > > Of course, with strong and weak data, the best > > > thing to do is compare the model-data disagreement > > > with the magnitude of the error. That is, compare > > > |Fobs-Fcalc| to sigma(Fobs), not Fobs itself. > > > Modern refinement programs do this! And I say > > > the more data the merrier. > > > > > > > > > -James Holton > > > MAD Scientist > > > > > > > > > On 3/4/2011 5:15 AM, Marjolein Thunnissen wrote: > > > > > > hi > > > > > > Recently on a paper I submitted, it was the > > > editor of the journal who wanted exactly the > > > same thing. I never argued with the editor > > > about this (should have maybe), but it could > > > be one cause of the epidemic that Bart Hazes > > > saw.... > > > > > > > > > best regards > > > > > > Marjolein > > > > > > On Mar 3, 2011, at 12:29 PM, Roberto > > > Battistutta wrote: > > > > > > Dear all, > > > I got a reviewer comment that indicate the > > > "need to refine the structures at an > > > appropriate resolution (I/sigmaI of>3.0), > > > and re-submit the revised coordinate files > > > to the PDB for validation.". In the > > > manuscript I present some crystal > > > structures determined by molecular > > > replacement using the same protein in a > > > different space group as search model. > > > Does anyone know the origin or the > > > theoretical basis of this "I/sigmaI>3.0" > > > rule for an appropriate resolution? > > > Thanks, > > > Bye, > > > Roberto. > > > > > > > > > Roberto Battistutta > > > Associate Professor > > > Department of Chemistry > > > University of Padua > > > via Marzolo 1, 35131 Padova - ITALY > > > tel. +39.049.8275265/67 > > > fax. +39.049.8275239 > > > roberto.battistu...@unipd.it > > > <mailto:roberto.battistu...@unipd.it> > > > www.chimica.unipd.it/roberto.battistutta/ > > > > > > <http://www.chimica.unipd.it/roberto.battistutta/> > > > VIMM (Venetian Institute of Molecular > > > Medicine) > > > via Orus 2, 35129 Padova - ITALY > > > tel. +39.049.7923236 > > > fax +39.049.7923250 > > > www.vimm.it <http://www.vimm.it> > > > > > > > > > > > > > > > > > > > > > > > > > > >
