On Jun 16, 2010, at 5:24 PM, Roger Koenker wrote:
I'm trying to understand why 64B R produces an insane answer/plot
to the first
example(crq)
in my quantreg package, whereas 32B R produces a sane answer
with the same data. This involves some fortran and in my usual
primitive fashion I've been trying to debug with fortran printing
using the tried and true call dblepr. However, in my test function
call dblepr prints values like
z
[1] 1.644287e-313
even though the function returns the sane, correct values. Is this
a known issue, and if so is there an alternative strategy for
printing from fortran in 64b R?
TIA for any enlightenment.
Roger
PS. In the hope that someone will tell me that all would be well
if I only followed directions and upgraded R, I'll sheepishly admit:
sessionInfo()
R version 2.11.0 Under development (unstable) (2010-02-18 r51149)
x86_64-apple-darwin9.8.0
With a freshly installed 2.11.1 just a couple of days ago, and I
only run the 64-bit version, and I do not see any such result. There
is a line on the first pot that "shoots up". I attached
it.<First.quantreg.example.plot.pdf>
Here is the console output I get:
> example(crq)
crq> # An artificial Powell example
crq> set.seed(2345)
crq> x <- sqrt(rnorm(100)^2)
crq> y <- -0.5 + x +(.25 + .25*x)*rnorm(100)
crq> plot(x,y, type="n")
Hit <Return> to see next plot:
crq> s <- (y > 0)
crq> points(x[s],y[s],cex=.9,pch=16)
crq> points(x[!s],y[!s],cex=.9,pch=1)
crq> yLatent <- y
crq> y <- pmax(0,y)
crq> yc <- rep(0,100)
crq> for(tau in (1:4)/5){
crq+ f <- crq(Curv(y,yc) ~ x, tau = tau, method = "Pow")
crq+ xs <- sort(x)
crq+ lines(xs,pmax(0,cbind(1,xs)%*%f$coef),col="red")
crq+ abline(rq(y ~ x, tau = tau), col="blue")
crq+ abline(rq(yLatent ~ x, tau = tau), col="green")
crq+ }
Loading required package: survival
Loading required package: splines
Attaching package: 'survival'
The following object(s) are masked from 'package:quantreg':
untangle.specials
crq> legend(.15,2.5,c("Naive QR","Censored QR","Omniscient QR"),
crq+ lty=rep(1,3),col=c("blue","red","green"))
crq> data(uis)
crq> #estimate the Peng and Huang model using log(TIME) AFT
specification
crq> fit <- crq(Surv(log(TIME), CENSOR) ~ ND1 + ND2 + IV3 +
crq+ TREAT + FRAC + RACE + AGE * SITE, method =
"Por", data = uis)
crq> Sfit <- summary(fit,1:19/20)
bootstrap roughly 10 percent complete
bootstrap roughly 20 percent complete
bootstrap roughly 30 percent complete
bootstrap roughly 40 percent complete
bootstrap roughly 50 percent complete
bootstrap roughly 60 percent complete
bootstrap roughly 70 percent complete
bootstrap roughly 80 percent complete
bootstrap roughly 90 percent complete
bootstrap roughly 100 percent complete
crq> PHit <- coxph(Surv(TIME, CENSOR) ~ ND1 + ND2 + IV3 +
crq+ TREAT + FRAC + RACE + AGE * SITE, data = uis)
crq> plot(Sfit, CoxPHit = PHit)
Hit <Return> to see next plot:
crq> formula <- ~ ND1 + ND2 + IV3 + TREAT + FRAC + RACE + AGE *
SITE -1
crq> X <- data.frame(model.matrix(formula,data=uis))
crq> newd <- as.list(apply(X,2,median))
crq> pred <- predict(fit, newdata=newd, type = "stepfun")
crq> plot(pred,do.points=FALSE,xlab = expression(tau), ylab =
expression(Q(tau)),main= "Quantiles at Median Covariate Values")
Hit <Return> to see next plot:
crq> plot(rearrange(pred),add=TRUE,do.points=FALSE,col.vert ="red",
col.hor="red")
crq> legend(.15,7,c("Raw","Rearranged"),lty =
rep(1,2),col=c("black","red"))
Warning messages:
1: In crq.fit.pow(X, y, cen, tau = taus, weights, left = left, ...) :
Solution may be nonunique
2: In crq.fit.pow(X, y, cen, tau = taus, weights, left = left, ...) :
Solution may be nonunique
3: In crq.fit.pow(X, y, cen, tau = taus, weights, left = left, ...) :
Solution may be nonunique
4: In crq.fit.pow(X, y, cen, tau = taus, weights, left = left, ...) :
Solution may be nonunique
>
None of the plots have weird limits. If you tell me what objects to
inspect I would be happy to deliver more details. I do not see that
dblepr() is in the NAMESPACE.
--
David.
PPS. The fact that I can run 32B R on the same machine by simply
typing R --arch i386 is one of the great advances of modern science
in my opinion, or at least I thought so, until I discovered this
discrepancy
in what was being computed by the two versions.
url: www.econ.uiuc.edu/~roger Roger Koenker
email rkoen...@uiuc.edu Department of Economics
vox: 217-333-4558 University of Illinois
fax: 217-244-6678 Urbana, IL 61801
_______________________________________________
R-SIG-Mac mailing list
R-SIG-Mac@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-mac
David Winsemius, MD
West Hartford, CT
