 |
 |
You do also have always to consider why you are doing this calculation - usually to satisfy a sceptical and possibly ill-informed referee. A major reason for doing this is to justify including an outer resolution shell of data (see this BB passim), and for this I have come to prefer the random half-dataset correlation coefficient in shells. A CC has a more straightforward distribution than an R-factor (though not entirely without problems). It is independent of the SD estimates, and easy to understand. Phil > > 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 <[log in to unmask] >> <mailto:[log in to unmask]>> 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 >> [log in to unmask] >> <mailto:[log in to unmask]> >> 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> >> >> >> >> >> >> >> > >