dblepr for 64B R
Hi David, Thanks, but this picture is exactly what I was getting in 64B..... and it is entirely bogus. Fortunately, this morning I discovered that ddot in fortran needs to be declared double in 64B, but apparently doesn't in 32B... and this resolves the problem. I've submitted an updated version of quantreg to CRAN that fixes the problem. The dblepr was just an intermediate debugging issue, and was entirely my stupidity. It works exactly as it should once variables are declared double. Why ddot behaves differently in 64B vs 32B is still mysterious to me, but so are many aspects of ordinary life, so I have resolved not to question further. Roger Koenker rkoenker at illinois.edu
On Jun 17, 2010, at 8:00 PM, David Winsemius wrote:
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 rkoenker at 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 at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-mac
David Winsemius, MD West Hartford, CT