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.
