Dear Ed, I've been reading this thread for a while, wondering if I should say anything. As Pavel knows, I've been suggesting that phenix.refine should include the experimental errors in the variance term ever since I realised they were being left out. On the other hand, Pavel has been asking (quite reasonably) for some evidence that it makes any difference. I feel that there ought to be circumstances where it does make a difference, but my attempts over the last few days to find a convincing test case have failed. Like you, I tried running Refmac with and without experimental sigmas, in my case on a structure where I know that we pushed the resolution limits (1bos), and I can't see any significant difference in the model or the resulting phases. Gabor Bunkoczi has suggested that it might be more relevant for highly anisotropic structures, so we'll probably look for a good example along those lines. In principle it should make a difference, but I think there's one point missing from your discussion below. If you leave out the experimental sigmas, then the refinement of sigmaA (or, equivalently, the alpha and beta parameters in phenix.refine) will lead to variances that already incorporate the average effect of experimental measurement error in each resolution shell. So if all reflections were measured to the same precision, it shouldn't matter at all to the likelihood target whether you add them explicitly or absorb them implicitly into the variances. The potential problems come when different reflections are measured with different levels of precision, e.g. because the data collection strategy gave widely varying redundancies for different reflections. In the statistics you give below, the key statistic is probably the standard deviation of sigf/sqrt(beta), which is actually quite small. So after absorbing the average effect of measurement error into the beta values, the residual variation is even less important to the total variance than you would think from the total value of sigf. I would still argue that it's relatively easy to incorporate the experimental error into the likelihood variances so it's worth doing even if we haven't found the circumstances where it turns out to matter! Regards, Randy Read On 12 Dec 2012, at 06:46, Ed Pozharski wrote:

On Tue, 2012-12-11 at 11:27 -0500, Douglas Theobald wrote:

...

What is the evidence, if any, that the exptl sigmas are actually negligible compared to fit beta (is it alluded to in Lunin 2002)? Is there somewhere in phenix output I can verify this myself?

Essentially, equation 4 in Lunin (2002) is the same as equation 14 in Murshudov (1997) or equation 1 in Cowtan (2005) or 12-79 in Rupp (2010). The difference is that instead of combination of sigf^2 and sigma_wc you have a single parameter, beta. One can do that assuming that sigf<sigf. I just pulled up a random dataset refined with phenix and here is what I see

min(sigf/sqrt(beta)) = 0.012 max(sigf/sqrt(beta)) = 0.851 mean(sigf/sqrt(beta)) = 0.144 std(sigf/sqrt(beta)) = 0.118

But there are two problems. First, in the highest resolution shell

min(sigf/sqrt(beta)) = 0.116 max(sigf/sqrt(beta)) = 0.851 mean(sigf/sqrt(beta)) = 0.339 std(sigf/sqrt(beta)) = 0.110

This is a bit more troubling. Notice that for acentrics it's 2sigf**2 +sigma_wc, thus the actual ratio should be increased by sqrt(2), getting uncomfortably close to 1/2. Still, given that one adds variances, this is at most 25% correction, and this *is* the high resolution shell.

Second, if one tries to interpret sqrt(beta) as a measure of model error in reciprocal space, one runs into trouble. This dataset was refined to R~18%. Assuming that sqrt(beta) should roughly predict discrepancy between Fo and Fc, it corresponds to R~30%. This suggests that for reasons I don't yet quite understand beta overestimates model variance. If it is simply doubled, then it becomes comparable to experimental error, at least in higher resolution shells.

...

And, in comparison, how does refmac handle the exptl sigmas? Maybe this last question is more appropriate for ccp4bb, but contrasting with phenix would be helpful for me. I know there's a box, checked by default, "Use exptl sigmas to weight Xray terms".

Refmac fits sigmaA to a certain resolution dependence and then adds experimental sigmas (or not as you noticed). I was told that the actual formulation is different from what is described in the original manuscript. But what's important that if one pulls out the sigma_wc as calculated by refmac it has all the same characteristics as sqrt(beta) - it is generally >>sigf and suggests model error in reciprocal space that is incompatible with (too large) observed R-values. Kevin Cowtan's spline approximation implemented in clipper libraries behaves much better, meaning that R-value expectations projected from sigma_wc are much closer to observed R-value.

Curiously, it does not make much difference in practice, i.e. refined model is not affected as much. For instance, with refmac there are no significant changes whether one uses experimental errors or not. I could think of several reasons for this, but haven't verified any.

Cheers,

Ed.

-- "I'd jump in myself, if I weren't so good at whistling." Julian, King of Lemurs _______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb

------ Randy J. Read Department of Haematology, University of Cambridge Cambridge Institute for Medical Research Tel: + 44 1223 336500 Wellcome Trust/MRC Building Fax: + 44 1223 336827 Hills Road E-mail: [email protected] Cambridge CB2 0XY, U.K. www-structmed.cimr.cam.ac.uk

On 12 Dec 2012, at 15:04, Douglas Theobald wrote:

On Dec 12, 2012, at 4:00 AM, Randy Read wrote:

...

Dear Ed,

I've been reading this thread for a while, wondering if I should say anything. As Pavel knows, I've been suggesting that phenix.refine should include the experimental errors in the variance term ever since I realised they were being left out. On the other hand, Pavel has been asking (quite reasonably) for some evidence that it makes any difference. I feel that there ought to be circumstances where it does make a difference, but my attempts over the last few days to find a convincing test case have failed. Like you, I tried running Refmac with and without experimental sigmas, in my case on a structure where I know that we pushed the resolution limits (1bos), and I can't see any significant difference in the model or the resulting phases. Gabor Bunkoczi has suggested that it might be more relevant for highly anisotropic structures, so we'll probably look for a good example along those lines.

In principle it should make a difference, but I think there's one point missing from your discussion below. If you leave out the experimental sigmas, then the refinement of sigmaA (or, equivalently, the alpha and beta parameters in phenix.refine) will lead to variances that already incorporate the average effect of experimental measurement error in each resolution shell. So if all reflections were measured to the same precision, it shouldn't matter at all to the likelihood target whether you add them explicitly or absorb them implicitly into the variances. The potential problems come when different reflections are measured with different levels of precision, e.g. because the data collection strategy gave widely varying redundancies for different reflections.

In the statistics you give below, the key statistic is probably the standard deviation of sigf/sqrt(beta), which is actually quite small. So after absorbing the average effect of measurement error into the beta values, the residual variation is even less important to the total variance than you would think from the total value of sigf.

So, IIUC, what you're suggesting is that the great majority of the variance of F/sigF can simply be explained by resolution, and that phenix betas (as currently implemented) already account for that.

I guess, to be more precise, I'm suggesting that this is probably true for my 1bos test case and the one Ed looked at, and maybe for lots of other examples. But I don't think it should cost much effort to avoid having to make this assumption, and there might be cases where it's not true, e.g. Gabor's example of highly anisotropic data. Randy ------ Randy J. Read Department of Haematology, University of Cambridge Cambridge Institute for Medical Research Tel: + 44 1223 336500 Wellcome Trust/MRC Building Fax: + 44 1223 336827 Hills Road E-mail: [email protected] Cambridge CB2 0XY, U.K. www-structmed.cimr.cam.ac.uk

To be honest, when I say "relatively easy" I'm thinking in terms of the single sigmaA parameter. Once you've got an estimate of sigmaA (e.g. by scaling the alpha parameter to apply to Es instead of Fs), then you can refine it using the likelihood target that has an explicit contribution from the sigf values. Or you could refine alpha and beta using the likelihood target expressed in terms of alpha and beta but adding the sigf contribution to beta. (But I believe my paper on SIGMAA showed that alpha and beta are not really independent anyway.) Randy On 12 Dec 2012, at 15:56, Ed Pozharski wrote:

Dear Randy,

On Wed, 2012-12-12 at 09:00 +0000, Randy Read wrote:

...

In the statistics you give below, the key statistic is probably the standard deviation of sigf/sqrt(beta), which is actually quite small. So after absorbing the average effect of measurement error into the beta values, the residual variation is even less important to the total variance than you would think from the total value of sigf.

You are absolutely right - this is one of the possible reasons (perhaps the main reason) why the effect on model upon incorporating sigf is not obvious. Indeed, in the example dataset that I used relative standard deviation of sigf in resolution shells ranges from 0.3 at high to 0.6 at low resolution. The incorporation of the shell-average sigf into beta is obscured by the fact that the two anti-correlate.

Another issue is that average value of beta does not matter that much and is only weakly controlled by data. This is obvious (to me but I may be wrong) from eq. (6-7) in Lunin (2002). It points out that near minimum (and we are talking here about effects on the *final model*) target may be approximated quadratically and the applied weight is essentially inversely proportional to beta (not exactly, of course, but it is the major effect beta has on the target). Thus, some inflation of beta over what it is expected to be (model variance in reciprocal space) will do very little to the target other than scaling it.

I think my main issue is with the idea that beta may be used as an estimate of model variance. Mathematically it probably does not matter, but we all tend to attach "physical" interpretations to model parameters, and here it does not work as it seems to suggest that crystallographic models are grossly overfitted.

...

I would still argue that it's relatively easy to incorporate the experimental error into the likelihood variances so it's worth doing even if we haven't found the circumstances where it turns out to matter!

I am not sure it is that easy. As I mentioned previously, the analytical expressions for alpha/beta break down when sigf is incorporated. Also, it is possible that this has already been done by Kevin Cowtan in his 2005 paper. My observation was that the spline coefficients return more reasonable estimates of model variance.

There are also very strange consequences in alpha/beta approach. Equations (6-7) in Lunin (2002) essentially set up the quadratic approximation. I already mentioned that target value Fs* oddly reduces to exact zero for "weak reflections", and if beta is overestimated those are not so weak anymore. In the example dataset that I used, in the low resolution shells average I/sigma of such zeroed reflections is as high as ~5-6. I can identify a reflection that according to eq. (6-7) will have a target value of zero and has I/sigma=22.5! I understand that this still has only minor effect on the final model because only ~12% of reflections are hit (and nobody minds that 5-10% of experimental data is routinely tossed for the sake of Rfree calculation). Still, it is puzzling and unexpected.

Also, the weights applied to individual reflections are approximated by ws* which has some peculiar properties, namely that it dips to zero around f/sqrt(beta)~1. Curiously, it recovers back to full weight for reflections that are weaker (Fig.3 in Lunin 2002). Again, I see no flaw in the math, but it is rather counterintuitive that reflections that roughly match model error are weighted down, while even weaker reflections are not. Given that these weaker reflections have their respective target Fs* reset to zero, there will be potential during minimization to simply run weak reflections down to zero across the board.

Oh, and by the way, phaser rocks :)

Ed.

-- Edwin Pozharski, PhD, Assistant Professor University of Maryland, Baltimore ---------------------------------------------- When the Way is forgotten duty and justice appear; Then knowledge and wisdom are born along with hypocrisy. When harmonious relationships dissolve then respect and devotion arise; When a nation falls to chaos then loyalty and patriotism are born. ------------------------------ / Lao Tse /

_______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb

------ Randy J. Read Department of Haematology, University of Cambridge Cambridge Institute for Medical Research Tel: + 44 1223 336500 Wellcome Trust/MRC Building Fax: + 44 1223 336827 Hills Road E-mail: [email protected] Cambridge CB2 0XY, U.K. www-structmed.cimr.cam.ac.uk

On Dec 12, 2012, at 1:46 AM, Ed Pozharski wrote:

On Tue, 2012-12-11 at 11:27 -0500, Douglas Theobald wrote:

...

What is the evidence, if any, that the exptl sigmas are actually negligible compared to fit beta (is it alluded to in Lunin 2002)? Is there somewhere in phenix output I can verify this myself?

Essentially, equation 4 in Lunin (2002) is the same as equation 14 in Murshudov (1997) or equation 1 in Cowtan (2005) or 12-79 in Rupp (2010). The difference is that instead of combination of sigf^2 and sigma_wc you have a single parameter, beta. One can do that assuming that sigf<

From the log file: |-----------------------------------------------------------------------------| |R-free likelihood based estimates for figures of merit, absolute phase error,| |and distribution parameters alpha and beta (Acta Cryst. (1995). A51, 880-887)| | | | Bin Resolution No. Refl. FOM Phase Scale Alpha Beta | | # range work test error factor | | 1: 44.4859 - 3.0705 14086 154 0.93 12.12 1.00 0.98 118346.13| | 2: 3.0705 - 2.4372 13777 149 0.91 15.26 1.00 0.99 58331.77| | 3: 2.4372 - 2.1291 13644 148 0.94 11.42 1.00 0.99 23216.31| it appears that phenix estimates alpha and beta from the R-free set rather than from the working set (I might be misreading that). Is that correct?

On Dec 12, 2012, at 10:47 AM, Randy Read wrote:

On 12 Dec 2012, at 15:36, Douglas Theobald wrote:

...

On Dec 12, 2012, at 1:46 AM, Ed Pozharski wrote:

...

On Tue, 2012-12-11 at 11:27 -0500, Douglas Theobald wrote:

...

What is the evidence, if any, that the exptl sigmas are actually negligible compared to fit beta (is it alluded to in Lunin 2002)? Is there somewhere in phenix output I can verify this myself?

Essentially, equation 4 in Lunin (2002) is the same as equation 14 in Murshudov (1997) or equation 1 in Cowtan (2005) or 12-79 in Rupp (2010). The difference is that instead of combination of sigf^2 and sigma_wc you have a single parameter, beta. One can do that assuming that sigf<

From the log file:

|-----------------------------------------------------------------------------| |R-free likelihood based estimates for figures of merit, absolute phase error,| |and distribution parameters alpha and beta (Acta Cryst. (1995). A51, 880-887)| | | | Bin Resolution No. Refl. FOM Phase Scale Alpha Beta | | # range work test error factor | | 1: 44.4859 - 3.0705 14086 154 0.93 12.12 1.00 0.98 118346.13| | 2: 3.0705 - 2.4372 13777 149 0.91 15.26 1.00 0.99 58331.77| | 3: 2.4372 - 2.1291 13644 148 0.94 11.42 1.00 0.99 23216.31|

it appears that phenix estimates alpha and beta from the R-free set rather than from the working set (I might be misreading that). Is that correct?

Yes, using the cross-validation data was a key step in getting maximum likelihood refinement to work. A long time ago (a few years before our first paper on ML refinement) I implemented a first version of the MLF target we put into CNS, but the sigmaA values were estimated from the working data. What happened was that the data would be over-fit, then the sigmaA estimates would go up (with part of the increase being a result of the overfitting), then in the next cycle the pressure to fit the data compared to the restraints would be higher, and so on. The best I could claim for this at the time was that the resulting models were at least as good as the ones from least-squares refinement, but the R-factors were higher (indicating that there was still less over-fitting). It would have been hard to sell the advantage of higher R-factors to the protein crystallography community so it was good that, when we started using cross-validated sigmaA values, the convergence radius improved and we could get significantly better models with lower R-factors. I think you'll find that all the programs use just the cross-validation data to estimate the variance parameters for the likelihood target, not just phenix.refine.

Thanks Randy. I must say I'm quite surprised by this, and coming from a likelihoodist/Bayesian POV it seems very wrong :). I realize it works in practice, and works quite well evidently, but there's a very odd marriage of statistical philosophies going on here. From a likelihood POV, the joint ML estimates really should come from the same data set (i.e. the working set --- and there's the issue that there's already been a likelihood "compromise" of sorts by excluding some of the data from estimation anyway, a practice which violates the likelihood principle). And from a frequentist cross-validation POV, you certainly should not be refining against the test set --- that violates the very rationale of a test set.

Randy

...

_______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb

------ Randy J. Read Department of Haematology, University of Cambridge Cambridge Institute for Medical Research Tel: + 44 1223 336500 Wellcome Trust/MRC Building Fax: + 44 1223 336827 Hills Road E-mail: [email protected] Cambridge CB2 0XY, U.K. www-structmed.cimr.cam.ac.uk

_______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb

On Dec 12, 2012, at 12:17 PM, Randy Read wrote:

Hi,

I'm not too worried about refining sigmaA values against the test set, since there's only one parameter per resolution shell, and it doesn't (directly) affect the Fcalc values.

That's all good, and pragmatics should probably take priority. But it is somewhat concerning that you have to do it in the first place --- the fact that a parameter blows up in a likelihood analysis is usually a good indication that something is fundamentally wrong, usually non-identifiability of some sort. Or that there's an error in implementation. I wonder if anyone has recently re-visited the issue with the latest software --- sometimes these problems just "disappear" when other issues get worked out.

Nonetheless, I agree that it would be better in principle to refine all parameters against the working data, and maybe we could do that if we understood how to correct the variance estimates in the likelihood targets for the number of degrees of freedom.

Randy

On 12 Dec 2012, at 16:04, Douglas Theobald wrote:

...

On Dec 12, 2012, at 10:47 AM, Randy Read wrote:

...

On 12 Dec 2012, at 15:36, Douglas Theobald wrote:

...

On Dec 12, 2012, at 1:46 AM, Ed Pozharski wrote:

...

On Tue, 2012-12-11 at 11:27 -0500, Douglas Theobald wrote:

...

What is the evidence, if any, that the exptl sigmas are actually negligible compared to fit beta (is it alluded to in Lunin 2002)? Is there somewhere in phenix output I can verify this myself?

Essentially, equation 4 in Lunin (2002) is the same as equation 14 in Murshudov (1997) or equation 1 in Cowtan (2005) or 12-79 in Rupp (2010). The difference is that instead of combination of sigf^2 and sigma_wc you have a single parameter, beta. One can do that assuming that sigf<

From the log file:

|-----------------------------------------------------------------------------| |R-free likelihood based estimates for figures of merit, absolute phase error,| |and distribution parameters alpha and beta (Acta Cryst. (1995). A51, 880-887)| | | | Bin Resolution No. Refl. FOM Phase Scale Alpha Beta | | # range work test error factor | | 1: 44.4859 - 3.0705 14086 154 0.93 12.12 1.00 0.98 118346.13| | 2: 3.0705 - 2.4372 13777 149 0.91 15.26 1.00 0.99 58331.77| | 3: 2.4372 - 2.1291 13644 148 0.94 11.42 1.00 0.99 23216.31|

it appears that phenix estimates alpha and beta from the R-free set rather than from the working set (I might be misreading that). Is that correct?

Yes, using the cross-validation data was a key step in getting maximum likelihood refinement to work. A long time ago (a few years before our first paper on ML refinement) I implemented a first version of the MLF target we put into CNS, but the sigmaA values were estimated from the working data. What happened was that the data would be over-fit, then the sigmaA estimates would go up (with part of the increase being a result of the overfitting), then in the next cycle the pressure to fit the data compared to the restraints would be higher, and so on. The best I could claim for this at the time was that the resulting models were at least as good as the ones from least-squares refinement, but the R-factors were higher (indicating that there was still less over-fitting). It would have been hard to sell the advantage of higher R-factors to the protein crystallography community so it was good that, when we started using cross-validated sigmaA values, the convergence radius ! improved and we could get significantly better models with lower R-factors. I think you'll find that all the programs use just the cross-validation data to estimate the variance parameters for the likelihood target, not just phenix.refine.

Thanks Randy. I must say I'm quite surprised by this, and coming from a likelihoodist/Bayesian POV it seems very wrong :). I realize it works in practice, and works quite well evidently, but there's a very odd marriage of statistical philosophies going on here. From a likelihood POV, the joint ML estimates really should come from the same data set (i.e. the working set --- and there's the issue that there's already been a likelihood "compromise" of sorts by excluding some of the data from estimation anyway, a practice which violates the likelihood principle). And from a frequentist cross-validation POV, you certainly should not be refining against the test set --- that violates the very rationale of a test set.

...

Randy

...

_______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb

------ Randy J. Read Department of Haematology, University of Cambridge Cambridge Institute for Medical Research Tel: + 44 1223 336500 Wellcome Trust/MRC Building Fax: + 44 1223 336827 Hills Road E-mail: [email protected] Cambridge CB2 0XY, U.K. www-structmed.cimr.cam.ac.uk

_______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb

_______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb

------ Randy J. Read Department of Haematology, University of Cambridge Cambridge Institute for Medical Research Tel: + 44 1223 336500 Wellcome Trust/MRC Building Fax: + 44 1223 336827 Hills Road E-mail: [email protected] Cambridge CB2 0XY, U.K. www-structmed.cimr.cam.ac.uk

_______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb