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.bu...@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@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]] ______________________________________________ R-help@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.