[BioC] Limma normalization error and loess.R segmentation type fault (windows)

[Gordon K Smyth]

Hi Marcus,

I have committed a fix to limma in Bioc devel repository, to work-around 
problems with lm.wfit() and loess() functions with zero weights. 
normalizeWithinArrays() now runs smoothly on the data example you sent me. 
Please try it out and see if it solves your problems.

My solution has been to disallow zero weights in loessFit(), instead to 
reset zero weights to a very small positive value so as to avoid 
instabilities in defining the fit.

My thinking on weights is that normalizeWithinArrays() should not 
introduce NA values for zero weight observations.  Rather, the philosophy 
is that weights can be used to define each probes influence on the 
normalizing curve, with the curve being extended to all observations 
regardless of their influence.  We use this sort of strategy for 
upweighting control probes, as for example in:

   http://genomebiology.com/2007/8/1/R2

Of course, the weights still enter into downstream functions such as 
lmFit(), and zero weights at this stage have the effect of removing 
observations entirely.

Best wishes
Gordon

> Date: Wed, 5 Sep 2012 16:16:40 +1200
> From: Marcus Davy <mdavy86 at gmail.com>
> To: Gordon K Smyth <smyth at wehi.edu.au>
> Cc: Bioconductor mailing list <bioconductor at r-project.org>
> Subject: Re: [BioC] Limma normalization error and loess.R segmentation
> 	type	fault (windows)
>
> Hi Gordon,
> some good points you make here.
>
> I have checked back to the original slide and there was a background smear
> over one corner covering an entire print tip region, so a researcher would
> have subjectively flagged all those spots as bad, which is why they were
> all allocated zero weight. For global loess normalization such as with
> agilent arrays, a researcher would also potentially remove the entire
> biological replicate from analysis.
>
> I am tending towards the solution should also be NA for the situation
> above. Maybe a conditional switch could be used to cover the options,
> ordinary unweighted normalization, and doing no normalization etc for
> historical/backward compatability reasons.
>
> Marcus
>
>
> On Wed, Sep 5, 2012 at 11:21 AM, Gordon K Smyth <smyth at wehi.edu.au> wrote:
>
>> Hi Marcus,
>>
>> I guess the fundamental question then is, what should loess normalization
>> return (for a array or print-tip group or whatever) when all the weights
>> are zero?
>>
>> Strict mathematics would suggest that the solution should be NA.
>>
>> Practical considerations suggest to me that normalizeWithinArrays() might
>> be better performing ordinary unweighted loess normalization in this case,
>> because no probe is weighted as more reliable than any other.
>>
>> Or alternatively, one might say that the loess curve can't be estimated,
>> so the raw expression values should be returned without adjustment.  So
>> loess normalization with zero weights is equivalent to no normalization.
>> That is what loessFit() and normalizeWithinArrays() have been doing up to R
>> 2.14.1.
>>
>> If you make the weights all zero for a print tip, do you want
>> normalizeWithinArrays to return NAs for all probes in that print tip group
>> on that array?  Or do ordinary unweighted normalization?  Or do no
>> normalization?
>>
>> Regards
>> Gordon
>>
>>
>> On Wed, 5 Sep 2012, Marcus Davy wrote:
>>
>>  It looks like this error is related to a particular print tip having all
>>> weights=0 as input into stats:::simpleLoess. I am trying to construct a
>>> simple reproducible example.
>>>
>>> Marcus
>>>

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}