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.