Have just checked with R 3.2.0 and MASS 7.3-40 and there still appears to
be a problem/strangeness at around k=2.5
On 23 April 2015 at 14:09, Francis Bursa <Francis.Bursa at quantics.co.uk>
wrote:
> Dear all,
>
> I believe I have found a bug in rlm in the MASS package. Specifically, the
> scale estimate can be wrong when there are no outliers. The following code
> snippet is an example:
>
> dose <- c(0,1,2,0,1,2)
> response <- c(0.659,1.633,3.621,1.803,3.093,4.424)
> line <- c(1,1,1,2,2,2)
> k2 <- seq(1.5,5,by=0.01)
> repNA <- rep(NA,length(k2))
> scale <- repNA
> niter <- repNA
> for (i in 1:length(k2)){
> rlm.fit <- rlm(response~dose+factor(line), psi=psi.huber,k=1.345,
> scale.est="proposal 2",k2=k2[i])
> scale[i] <- rlm.fit$s
> niter[i] <- length(rlm.fit$conv)
> }
> plot(k2,scale,type="b",col=niter)
>
> For this dataset there are no outliers, so I would expect the scale to be
> a smooth function of k2 once k2 is reasonably large, certainly for k2 >
2.
> However, there is a funny jump in the scale estimate around k2 = 2.4, just
> at the point where the number of iterations to convergence falls from 3 to
> 1.
>
> Looking at the source code, it appears that on each iteration, the scale
> is updated, then the parameters, and then a check for convergence is
> carried out just for the parameters, not the scale. So I would guess that
> in the range around k2=2.5 convergence is being reached when in fact the
> scale estimate hasn't converged.
>
> I am using MASS version 7.3-33 and R version 3.1.0 on Windows.
>
> I am not sure how common this issue is but there does not seem to be
> anything special about my dataset so it could be quite generic. Am I right
> that this is a bug?
>
> Many thanks,
> Francis Bursa
>
> ------------------
> Francis Bursa
> Statistician
>
> Quantics Consulting Ltd
> 28 Drumsheugh Gardens
> Edinburgh
> EH3 7RN
>
> Telephone: +44 (0) 131 440 2781 ext 207
>
> Quantics - complex data into clear results
> Quantics is an ISO 9001 registered company
>
> www.quantics.co.uk
>
> Please note that the contents of this e-mail (including any attachments)
> are confidential and may be legally privileged. If you are not the intended
> recipient you may not read, copy, distribute or make any other use of this
> email or its contents. If received in error, please tell us immediately by
> telephone on +44 (0) 131 440 2781 quoting the name of the sender and the
> intended recipient, then delete it from your system. Thank you.
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]