Full_Name: Charles J. Geyer
Version: 1.3.1
OS: linux-gnu-i686
Submission from: (NULL) (134.84.86.22)
ltsreg gives incomprehensible (to me) warnings
A homework problem for nonparametrics
########## start example ##########
library(bootstrap)
data(cell)
names(cell)
attach(cell)
library(lqs)
plot(V1, V2)
fred <- ltsreg(V2 ~ V1 + I(V1^2))
curve(predict(fred, data.frame(V1 = x)), add = TRUE)
coefficients(fred)
beta1.hat <- coefficients(fred)[2]
beta2.hat <- coefficients(fred)[3]
n <- length(V1)
nboot <- 1000
beta1.star <- double(nboot)
beta2.star <- double(nboot)
for (i in 1:nboot) {
k <- sample(1:n, replace = TRUE)
x.star <- V1[k]
y.star <- V2[k]
sally <- ltsreg(y.star ~ x.star + I(x.star^2))
curve(predict(sally, data.frame(x.star = x)),
add = TRUE, col = 2)
beta1.star[i] <- coefficients(sally)[2]
beta2.star[i] <- coefficients(sally)[3]
}
points(V1, V2, pch = 16)
curve(predict(fred, data.frame(V1 = x)),
add = TRUE, lwd = 2)
beta1.hat
sd(beta1.star)
beta2.hat
sd(beta2.star)
beta2.hat / sd(beta2.star)
########## end example ##########
This produces at least 50 warnings. The problem seems to be line 101 in
lqs.default
s <- sqrt(z$crit/quantile)/sqrt(1 - 2 * n * dnorm(1/c1)/(quantile * c1))
where z$crit is sometimes negative (like -1e-20) and so sqrt complains.
Now I have no idea what ltsreg is doing, and what z$crit is. The letters
suggest
a z critical value, which can't be negative. Certainly, the programmer who
decided
to stuff it into sqrt() didn't think it would be negative.
So this suggests a problem in the lqs_fitlots function call, but I haven't
looked
at that.
I understand that the problem may be that bootstrapping such a small data set
is insane, but I'm just copying Efron and Tibshirani (who do this with
lmsreg).
So the problem may be only that more informative error messages are needed (or
failure codes, although failure codes are somewhat antithetical to the S/R
spirit).
Anyway, just stuffing negative numbers into sqrt() doesn't look good.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
