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.
Possible bug in rlm
2 messages · Francis Bursa, Richard Perry
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.