Re: [R] particle count probability
Somewhere, buried in the vast literature on the Wicksell problem, there is probably an answer, or at least a hint. > On Feb 20, 2019, at 11:16 AM, PIKAL Petr wrote: > > Dear all > > Sorry, this is probably the most off-topic mail I have ever sent to this help > list. However maybe somebody could point me to right direction or give some > advice. > > In microscopy particle counting you have finite viewing field and some > particles could be partly outside of this field. My problem/question is: > > Do bigger particles have also bigger probability that they will be partly > outside this viewing field than smaller ones? > > Saying it differently, although there is equal count of bigger (white) and > smaller (black) particles in enclosed picture (8), due to the fact that more > bigger particles are on the edge I count more small particles (6) than big > (4). > > Is it possible to evaluate this feature exactly i.e. calculate some bias > towards smaller particles based on particle size distribution, mean particle > size and/or image magnification? > > Best regards > Petr Pikal > Osobn? ?daje: Informace o zpracov?n? a ochran? osobn?ch ?daj? obchodn?ch > partner? PRECHEZA a.s. jsou zve?ejn?ny na: > https://www.precheza.cz/zasady-ochrany-osobnich-udaju/ | Information about > processing and protection of business partner's personal data are available > on website: https://www.precheza.cz/en/personal-data-protection-principles/ > D?v?rnost: Tento e-mail a jak?koliv k n?mu p?ipojen? dokumenty jsou d?v?rn? a > podl?haj? tomuto pr?vn? z?vazn?mu prohl??en? o vylou?en? odpov?dnosti: > https://www.precheza.cz/01-dovetek/ | This email and any documents attached > to it may be confidential and are subject to the legally binding disclaimer: > https://www.precheza.cz/en/01-disclaimer/ > > __ > R-help@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. __ R-help@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.
Re: [R] Siegel nonparametric regression / mblm package
A quick look at the code for Siegel in mblm reveals that it is extremely inefficient, but it seems to be correct. One “explanation” for this behavior, presuming that we haven’t overlooked something more basic, is that such high breakdown estimates sacrifice some efficiency, that is to say, they are more variable than other methods when the data is well behaved, and of course, the Galton data is famously “almost Gaussian”. > On Feb 11, 2019, at 12:47 PM, Marco Besozzi wrote: > > Thank you very much for your reply. > If I have well understood, unfortunately in this way I have lost the only > idea I had... > Do you believe that a problem in the R algorithm employed in the package mblm > for Siegel regression is possible? > And do you know if Siegel regression is available in a different package? I > was unable to find it. > Thanks again! > Best regards. > > P.S.: sorry for my bad english... > > Il giorno lun 11 feb 2019 alle ore 12:54 Roger Koenker <mailto:rkoen...@illinois.edu>> ha scritto: > My first thought was also that this was an artifact of the ties, but > dithering the data > n <- length(child) > child <- child + runif(n,-.5,.5) > parent <- parent + runif(n,-.5,.5) > > and rerunning yields the same discrepancy between the Siegel and other fits. > Curiously, both > lmsreg and ltsreg from MASS produce lines that are more steeply sloped than > those > of the other methods. Since I stupidly forgot to set.seed(), YMMV. > > > On Feb 11, 2019, at 10:24 AM, Marco Besozzi > <mailto:marco.bes...@gmail.com>> wrote: > > > > I employed the "galton" set of data included in the package "psych". With > > the package "mblm" I obtained the Theil-Sen nonparametric regression and > > the Siegel non parametric regression, and compared them with the ordinary > > least square regression line. > > The results of standard regression and Theil-Sen regression are practically > > identical. But the Siegel regression seems to have a bias that I cannot > > understand. May I ask for a possible explanation? The bias may be related > > to the number of ties in the set of data? Here's the code and the image. > > > > Best regards. > > > > Marco Besozzi > > # Theil-Sen and Siegel nonparametric regression with package mblm > > # comparison with ordinary least squares (parametric) regression > > # on galton set of data included in the package psych > > # > > library(psych) > > attach(galton) > > library(mblm) > > # > > reglin_yx <- lm(child ~ parent, data=galton) # ordinary least squares > > (parametric) regression > > a_yx <- reglin_yx$coefficients[1] # intercept a > > b_yx <- reglin_yx$coefficients[2] # slope b > > # > > regnonTS <- mblm(child ~ parent, data=galton, repeated=FALSE) # Theil-Sen > > nonparametric regression (wait a few minutes!) > > a_TS <- regnonTS$coefficients[1] # intercept a > > b_TS <- regnonTS$coefficients[2] # slope b > > # > > regnonS = mblm(child ~ parent, data=galton, repeated=TRUE) # Siegel > > nonparametric regression > > a_S <- regnonS$coefficients[1] # intercept a > > b_S <- regnonS$coefficients[2] # slope b > > # > > # xy plot of data and regression lines > > # > > windows() # open a new window > > plot(parent, child, xlim = c(60,80), ylim = c(60,80), pch=1, xlab="Parent > > heigt (inch)", ylab="Chile height (inch)", main="Regression lines > > comparison", cex.main = 0.9) # data plot > > abline(a_yx, b_yx, col="green", lty=1) # ordinary least squares > > (parametric) regression line > > abline(a_TS, b_TS, col="blue", lty=1) # Theil-Sen nonparametric regression > > line > > abline(a_S, b_S, col="red", lty=1) # Siegel nonparametric regression > > legend(60, 80, legend=c("Ordinary least squares regression", "Theil-Sen > > nonparametric regression","Siegel nonparametric regression"), > > col=c("green", "blue", "red"), lty=c(4,4,1), cex=0.8) # add a legend > > # > > __ > > R-help@r-project.org <mailto:R-help@r-project.org> mailing list -- To > > UNSUBSCRIBE and more, see > > https://stat.ethz.ch/mailman/listinfo/r-help > > <https://stat.ethz.ch/mailman/listinfo/r-help> > > PLEASE do read the posting guide > > http://www.R-project.org/posting-guide.html > > <http://www.r-project.org/posting-guide.html> > > and provide commented, minimal, self-contained, reproducible code. > [[alternative HTML version deleted]] __ R-help@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.
Re: [R] Siegel nonparametric regression / mblm package
My first thought was also that this was an artifact of the ties, but dithering the data n <- length(child) child <- child + runif(n,-.5,.5) parent <- parent + runif(n,-.5,.5) and rerunning yields the same discrepancy between the Siegel and other fits. Curiously, both lmsreg and ltsreg from MASS produce lines that are more steeply sloped than those of the other methods. Since I stupidly forgot to set.seed(), YMMV. > On Feb 11, 2019, at 10:24 AM, Marco Besozzi wrote: > > I employed the "galton" set of data included in the package "psych". With > the package "mblm" I obtained the Theil-Sen nonparametric regression and > the Siegel non parametric regression, and compared them with the ordinary > least square regression line. > The results of standard regression and Theil-Sen regression are practically > identical. But the Siegel regression seems to have a bias that I cannot > understand. May I ask for a possible explanation? The bias may be related > to the number of ties in the set of data? Here's the code and the image. > > Best regards. > > Marco Besozzi > # Theil-Sen and Siegel nonparametric regression with package mblm > # comparison with ordinary least squares (parametric) regression > # on galton set of data included in the package psych > # > library(psych) > attach(galton) > library(mblm) > # > reglin_yx <- lm(child ~ parent, data=galton) # ordinary least squares > (parametric) regression > a_yx <- reglin_yx$coefficients[1] # intercept a > b_yx <- reglin_yx$coefficients[2] # slope b > # > regnonTS <- mblm(child ~ parent, data=galton, repeated=FALSE) # Theil-Sen > nonparametric regression (wait a few minutes!) > a_TS <- regnonTS$coefficients[1] # intercept a > b_TS <- regnonTS$coefficients[2] # slope b > # > regnonS = mblm(child ~ parent, data=galton, repeated=TRUE) # Siegel > nonparametric regression > a_S <- regnonS$coefficients[1] # intercept a > b_S <- regnonS$coefficients[2] # slope b > # > # xy plot of data and regression lines > # > windows() # open a new window > plot(parent, child, xlim = c(60,80), ylim = c(60,80), pch=1, xlab="Parent > heigt (inch)", ylab="Chile height (inch)", main="Regression lines > comparison", cex.main = 0.9) # data plot > abline(a_yx, b_yx, col="green", lty=1) # ordinary least squares > (parametric) regression line > abline(a_TS, b_TS, col="blue", lty=1) # Theil-Sen nonparametric regression > line > abline(a_S, b_S, col="red", lty=1) # Siegel nonparametric regression > legend(60, 80, legend=c("Ordinary least squares regression", "Theil-Sen > nonparametric regression","Siegel nonparametric regression"), > col=c("green", "blue", "red"), lty=c(4,4,1), cex=0.8) # add a legend > # > __ > R-help@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. __ R-help@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.
Re: [R] Summarizing R script
I use R CMD BATCH foo which produces a file called foo.Rout and provided the script includes sessionInfo() constitutes a quite sufficient summary for my purposes, it isn’t exactly pretty, but it is informative. > On Sep 26, 2018, at 3:00 PM, Spencer Brackett > wrote: > > R users, > > Is anyone aware of the proper procedure for summarizing a script(your > complete list of functions, arguments , and error codes within your R > console for say a formal report or publication? > > Many thanks, > > Best wishes, > > Spencer Brackett > > -- Forwarded message - > From: CHATTON Anne via R-help > Date: Wed, Sep 26, 2018 at 6:03 AM > Subject: [R] Problems to obtain standardized betas in multiply-imputed data > To: r-help@r-project.org > > > Dear all, > > I am having problems in obtaining standardized betas on a multiply-imputed > data set. Here are the codes I used : > imp = mice(data, 5, maxit=10, seed=42, print=FALSE) > FitImp <- with(imp,lm(y ~ x1 + x2 + x3)) > Up to here everything is fine. But when I ask for the standardized > coefficients of the multiply-imputed regressions using this command : > sdBeta <- lm.beta(FitImp) > I get the following error message: > Error in b * sx : argument non numérique pour un opérateur binaire > > Can anyone help me with this please? > > Anne > > __ > R-help@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. > > [[alternative HTML version deleted]] > > __ > R-help@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. __ R-help@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.
Re: [R] Scaling - does it get any better results than not scaling?
In certain fields this sort of standardization has become customary based on some sort of (misguided) notion that it induces “normality.” For example, in anthropometric studies based on the international Demographic and Health Surveys (DHS) childrens’ heights are often transformed to Z-scores prior to subsequent analysis under the dubious presumption that variability around the Z-scores at various ages will be Gaussian. In my experience this is rarely justified, and analysts would be better off modeling the original data rather than doing the preliminary transformation. This is discussed in further detail here: https://projecteuclid.org/euclid.bjps/1313973394. > On Jul 17, 2018, at 5:53 AM, Michael Thompson > wrote: > > Hi, > I seem to remember from classes that one effect of scaling / standardising > data was to get better results in any analysis. But what I'm seeing when I > study various explanations on scaling is that we get exactly the same > results, just that when we look at standardised data it's easier to see > proportionate effects. > This is all very well for the data scientist to further investigate, but from > a practical point of view, (especially IF it doesn't improve the accuracy of > the result) surely it adds complication to 'telling the story' > of the model to non-DS people? > So, is scaling a technique for the DS to use to find effects, while > eventually delivering a non-scaled version to the users? > I'd like to be able to give the true story to my students, not some fairy > story based on my misunderstanding. Hope you can help with this. > Michael > > __ > R-help@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. __ R-help@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.
Re: [R] Lattice Histogram Scaling
My apologies, the data can now be found at: url <- "http://www.econ.uiuc.edu/~roger/research/ebayes/velo.d; x <- scan(url,skip = 1) If I could get each of the histograms to mimic what is produced by hist(x, 100, freq = FALSE) I’ve experimented with xlim, ylim, without success so far... url:www.econ.uiuc.edu/~roger Roger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 > On Aug 14, 2017, at 8:36 PM, Richard M. Heiberger <r...@temple.edu> wrote: > > Your example is not easily reproducible. > The REBayes requires Rmosek which requires a system command MOSEK. > Please try again with an example using data in base R. > > Meanwhile, my guess is that you will need to do something like > explicitly specifying xlim and ylim so all panels have the same > limits. > > On Mon, Aug 14, 2017 at 5:41 PM, Roger Koenker <rkoen...@illinois.edu> wrote: >> I am trying to do some comparisons of density estimators using lattice. >> The code below attempts to plot the same histogram in each panel and >> then overplots a kernel estimate with different bandwidths. Finding >> packet.number() was a bit arduous, but seems to do what I want. My >> concern now is that close examination of the 4 histograms reveals that >> they are different even though they use the same data, and use the >> same binning. Can someone explain this, or better yet suggest a fix? >> Admittedly, the kernel estimates are rather silly, they are just standing >> in for something else that I would like to think is less silly. >> >> Many thanks, >> Roger >> >> >> require(REBayes) # Source for the velo data >> require(lattice) >> x <- velo >> x <- x[!(x == 0)] >> bandwidths <- (1:4)*10 >> lp <- histogram(~x|bw, nint = 100, border = grey(.9), col = grey(.9), >>type = "density", panel = function(x, bw = bandwidths, ...){ >>panel.histogram(x, ...) >>f <- density(x, bw[packet.number()]) >>panel.lines(f$x, f$y, col = "blue", lwd = 1.5) >>}) >> print(lp) >> >> >> url:www.econ.uiuc.edu/~rogerRoger Koenker >> emailrkoen...@uiuc.eduDepartment of Economics >> vox: 217-333-4558University of Illinois >> fax: 217-244-6678Urbana, IL 61801 >> >> __ >> R-help@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. __ R-help@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.
[R] Lattice Histogram Scaling
I am trying to do some comparisons of density estimators using lattice. The code below attempts to plot the same histogram in each panel and then overplots a kernel estimate with different bandwidths. Finding packet.number() was a bit arduous, but seems to do what I want. My concern now is that close examination of the 4 histograms reveals that they are different even though they use the same data, and use the same binning. Can someone explain this, or better yet suggest a fix? Admittedly, the kernel estimates are rather silly, they are just standing in for something else that I would like to think is less silly. Many thanks, Roger require(REBayes) # Source for the velo data require(lattice) x <- velo x <- x[!(x == 0)] bandwidths <- (1:4)*10 lp <- histogram(~x|bw, nint = 100, border = grey(.9), col = grey(.9), type = "density", panel = function(x, bw = bandwidths, ...){ panel.histogram(x, ...) f <- density(x, bw[packet.number()]) panel.lines(f$x, f$y, col = "blue", lwd = 1.5) }) print(lp) url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 __ R-help@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.
Re: [R] Help with a specific quantile regression problem
What ensures that the tau-th quantile of the residuals is (nearly) zero, is that there IS an intercept in the model, this is one of the conditions required for the subgradient to contain 0 provided there is an intercept, when there is no intercept there is constraint to enforce this any more. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 > On Feb 24, 2017, at 5:10 AM, Helio Santana <helio.santana.1...@gmail.com> > wrote: > > Dear R community, > > I am a beginner in quantile regression and I have a question about a > specific problem. I have used the quantreg package to fit a QR with > inequality constrains: > > n <- 100 > p <- 5 > X <- matrix(rnorm(n*p),n,p) > y <- 0.95*apply(X,1,sum)+rnorm(n) > R <- cbind(0,rbind(diag(p),-diag(p))) > r <- c(rep(0,p),-rep(1,p)) > model <- rq(y~X,R=R,r=r,method="fnc") > > So, > >> quantile(model$residuals,0.5) >> -6.68836e-11 (It should be close to 0) > > > However, if I try to impose no intercept in the last problem: > > R <- cbind(0,rbind(diag(p),-diag(p))) > R <- R[,2:dim(R)[2]] > r <- c(rep(0,p),-rep(1,p)) > model <- rq(y~X-1,R=R,r=r,method="fnc") > > I obtain: > >> quantile(model$residuals,0.5) >> -0.03465427 > > As you can see, this quantile value is not close to 0 as I expected. Have I > an error in the formulation of the QR? > > Is it possible to fit a QR with inequality constrains and no intercept? > > Is there another alternative for solving this kind of problem? > > I would appreciate your comments. > > Best regards, > > Helio > > [[alternative HTML version deleted]] > > __ > R-help@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. __ R-help@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.
Re: [R] Functional programming?
Thanks again to all. This was highly educational for me. As it turned out Mikiko's suggestion turned out to be most easily adaptable to my real problem. I’m not at all able to evaluate relative efficiency at this point, but fortunately this isn’t an important factor (so far) in what I’m trying to do. As always the knowledge and generosity of the R-help community is inspiring. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 > On Mar 2, 2016, at 1:25 PM, Mikko Korpela <mikko.korp...@aalto.fi> wrote: > > On 02.03.2016 18:47, Roger Koenker wrote: >> I have a (remarkably ugly!!) code snippet (below) that, given >> two simple functions, f and g, generates >> a list of new functions h_{k+1} = h_k * g, k= 1, …, K. Surely, there are >> vastly >> better ways to do this. I don’t particularly care about the returned list, >> I’d be happy to have the final h_K version of the function, >> but I keep losing my way and running into the dreaded: >> >> Error in h[[1]] : object of type 'closure' is not subsettable >> or >> Error: evaluation nested too deeply: infinite recursion / >> options(expressions=)? >> >> Mainly I’d like to get rid of the horrible, horrible paste/parse/eval evils. >> Admittedly >> the f,g look a bit strange, so you may have to suspend disbelief to imagine >> that there is >> something more sensible lurking beneath this minimal (toy) example. >> >>f <- function(u) function(x) u * x^2 >>g <- function(u) function(x) u * log(x) >>set.seed(3) >>a <- runif(5) >>h <- list() >>hit <- list() >>h[[1]] <- f(a[1]) >>hit[[1]] <- f(a[1]) >>for(i in 2:5){ >> ht <- paste("function(x) h[[", i-1, "]](x) * g(", a[i], ")(x)") >> h[[i]] <- eval(parse(text = ht)) >> hit[[i]] <- function(x) {force(i); return(h[[i]] (x))} >> } >>x <- 1:99/10 >>plot(x, h[[1]](x), type = "l") >>for(i in 2:5) >> lines(x, h[[i]](x), col = i) > > Here is my (ugly?) suggestion: > > f <- function(u) function(x) u * x^2 > g <- function(u) function(x) u * log(x) > set.seed(3) > a <- runif(5) > h <- f(a[1]) > for (i in 2:5) { >body(h) <- call("*", body(h), >as.call(list(do.call("g", list(a[i])), quote(x > } > > -- > Mikko Korpela > Aalto University School of Science > Department of Computer Science __ R-help@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.
Re: [R] Functional programming?
Thanks, Duncan and Bert, Duncan’s version does replicate my result, Bert’s does something a bit different, now I just need some time to digest what you have done, and try to see how and why. Many thanks!!! Roger url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 > On Mar 2, 2016, at 12:23 PM, Duncan Murdoch <murdoch.dun...@gmail.com> wrote: > > On 02/03/2016 11:47 AM, Roger Koenker wrote: >> I have a (remarkably ugly!!) code snippet (below) that, given >> two simple functions, f and g, generates >> a list of new functions h_{k+1} = h_k * g, k= 1, …, K. Surely, there are >> vastly >> better ways to do this. I don’t particularly care about the returned list, >> I’d be happy to have the final h_K version of the function, >> but I keep losing my way and running into the dreaded: >> >> Error in h[[1]] : object of type 'closure' is not subsettable >> or >> Error: evaluation nested too deeply: infinite recursion / >> options(expressions=)? >> >> Mainly I’d like to get rid of the horrible, horrible paste/parse/eval evils. >> Admittedly >> the f,g look a bit strange, so you may have to suspend disbelief to imagine >> that there is >> something more sensible lurking beneath this minimal (toy) example. >> >> f <- function(u) function(x) u * x^2 >> g <- function(u) function(x) u * log(x) >> set.seed(3) >> a <- runif(5) >> h <- list() >> hit <- list() >> h[[1]] <- f(a[1]) >> hit[[1]] <- f(a[1]) >> for(i in 2:5){ >> ht <- paste("function(x) h[[", i-1, "]](x) * g(", a[i], ")(x)") >> h[[i]] <- eval(parse(text = ht)) >> hit[[i]] <- function(x) {force(i); return(h[[i]] (x))} >> } >> x <- 1:99/10 >> plot(x, h[[1]](x), type = "l") >> for(i in 2:5) >> lines(x, h[[i]](x), col = i) > > I don't understand what "hit" is for, but something like this should do it: > > > hlist <- function(maxk, f,g,a) { > h <- list() > h[[1]] <- f(a[1]) > for (j in 2:maxk) { > h[[j]] <- local({ > k <- j > function(x) { > result <- h[[1]](x) > for (i in 2:k) { > result <- result*g(a[i])(x) > } > result > } > }) > } > h > } > > f <- function(u) function(x) u * x^2 > g <- function(u) function(x) u * log(x) > set.seed(3) > a <- runif(5) > h <- hlist(5, f, g, a) __ R-help@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.
[R] Functional programming?
I have a (remarkably ugly!!) code snippet (below) that, given two simple functions, f and g, generates a list of new functions h_{k+1} = h_k * g, k= 1, …, K. Surely, there are vastly better ways to do this. I don’t particularly care about the returned list, I’d be happy to have the final h_K version of the function, but I keep losing my way and running into the dreaded: Error in h[[1]] : object of type 'closure' is not subsettable or Error: evaluation nested too deeply: infinite recursion / options(expressions=)? Mainly I’d like to get rid of the horrible, horrible paste/parse/eval evils. Admittedly the f,g look a bit strange, so you may have to suspend disbelief to imagine that there is something more sensible lurking beneath this minimal (toy) example. f <- function(u) function(x) u * x^2 g <- function(u) function(x) u * log(x) set.seed(3) a <- runif(5) h <- list() hit <- list() h[[1]] <- f(a[1]) hit[[1]] <- f(a[1]) for(i in 2:5){ ht <- paste("function(x) h[[", i-1, "]](x) * g(", a[i], ")(x)") h[[i]] <- eval(parse(text = ht)) hit[[i]] <- function(x) {force(i); return(h[[i]] (x))} } x <- 1:99/10 plot(x, h[[1]](x), type = "l") for(i in 2:5) lines(x, h[[i]](x), col = i) Thanks, Roger __ R-help@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.
Re: [R] algorithmic method quantile regression
Did you read item 1 in the quantreg FAQ()? url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 > On Oct 14, 2015, at 2:56 PM, T.Riedle <tr...@kent.ac.uk> wrote: > > Greetings R Community, > I am trying to run a quantile regression using the quantreg package. My > regression includes 7 independent variables with approx. 800 daily > observations each. Thus, I think that the Barrodale and Roberts algorithm > should do the trick. However, the Frisch-Newton after preprocessing returns > different results and more significant coefficients than the br method. Which > algorithmic method should I use now? Do the results mean that the > Frisch-Newton after preprocessing dominates the br method? > > [[alternative HTML version deleted]] > > __ > R-help@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. __ R-help@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.
Re: [R] quantile regression: warning message
see the output from the quantreg FAQ: FAQ() especially point 2. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 > On Oct 13, 2015, at 2:55 PM, T.Riedle <tr...@kent.ac.uk> wrote: > > Greetings R Community, > I am trying to run a quantile regression using the quantreg package. My code > looks as follows: > RegressionUtilitiesUK<-rq(ReturnUtilities~yield.spread.change+ReturnFTSE, > tau=0.01,data=State_variables_UK_calm) > > Unfortunately, the summary() function returns the results but also following > warning message: > Warning message: > In rq.fit.br(x, y, tau = tau, ...) : Solution may be nonunique > > My question now is if I should worry about the results. Are my results > correct and how can I avoid this message? I do not understand the message. > > Thanks a lot for your feedback. > > [[alternative HTML version deleted]] > > __ > R-help@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. __ R-help@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.
Re: [R] Quantile Regression without intercept
> On Oct 7, 2015, at 3:46 PM, Preetam Pal <lordpree...@gmail.com> wrote: > > So, what does weighted quantile regression even aim to achieve? Invariably, > this plane would not split the data set into the requisite fractions….. This officially ties this thread into a loop, since the OP wants to know why he wanted to do what he originally requested. If one “knows” that the intercept is zero, then it is always good to impose this, but like I said one does this at one’s peril. An example from economics might help: suppose you are estimating Engel curves in an unwelfare state, so 0 income implies 0 expenditure, then all (quantile) Engel curves pass through the origin and one might want to impose this. On the other hand maybe not... > From: Roger Koenker > Sent: 06-10-2015 07:09 PM > To: Lorenz, David > Cc: r-help@r-project.org > Subject: Re: [R] Quantile Regression without intercept > > > > On Oct 6, 2015, at 8:32 AM, Lorenz, David <lor...@usgs.gov> wrote: > > > > Thanks for the details, I suspected something like that. > > I think that begs the question: what is the meaning of quantile regression > > through the origin? If the tau=.5 line does not pass through 1/2 the data > > how do I interpret the line? > > As an estimate of the conditional median (quantile) function when constrained > to pass through > the origin… as with least squares fitting without an intercept, you do this > at your peril. > > > > > > On Tue, Oct 6, 2015 at 8:03 AM, Roger Koenker <rkoen...@illinois.edu> wrote: > > > > > On Oct 6, 2015, at 7:58 AM, Lorenz, David <lor...@usgs.gov> wrote: > > > > > > Did you verify that the correct percentages were above/below the > > > regression > > > lines? I did a quick check and for example did not consistently get 50% of > > > the observed response values greater than the tau=.5 line. I did when I > > > included the nonzero intercept term. > > > > Your "correct percentages" are only correct when you have an intercept in > > the model, > > without an intercept there is no gradient condition to ensure that. > > > > > > > > > > > >> Date: Mon, 5 Oct 2015 21:14:04 +0530 > > >> From: Preetam Pal <lordpree...@gmail.com> > > >> To: stephen sefick <ssef...@gmail.com> > > >> Cc: "r-help@r-project.org" <r-help@r-project.org> > > >> Subject: Re: [R] Quantile Regression without intercept > > >> Message-ID: <56129a41.025f440a.b1cf4.f...@mx.google.com> > > >> Content-Type: text/plain; charset="UTF-8" > > >> > > >> Yes..it works. Thanks ?? > > >> > > >> -Original Message- > > >> From: "stephen sefick" <ssef...@gmail.com> > > >> Sent: ?05-?10-?2015 09:01 PM > > >> To: "Preetam Pal" <lordpree...@gmail.com> > > >> Cc: "r-help@r-project.org" <r-help@r-project.org> > > >> Subject: Re: [R] Quantile Regression without intercept > > >> > > >> I have never used this, but does the formula interface work like lm? > > >> Y~X-1? > > >> > > >> > > >> On Mon, Oct 5, 2015 at 10:27 AM, Preetam Pal <lordpree...@gmail.com> > > >> wrote: > > >> > > >> Hi guys, > > >> > > >> Can you instruct me please how to run quantile regression without the > > >> intercept term? I only know about the rq function under quantreg package, > > >> but it automatically uses an intercept model. Icant change that, it > > >> seems. > > >> > > >> I have numeric data on Y variable (Gdp) and 2 X variables (Hpa and > > >> Unemployment). Their sizes are 125 each. > > >> > > >> Appreciate your help with this. > > >> > > >> Regards, > > >> Preetam > > >>[[alternative HTML version deleted]] > > >> > > >> __ > > >> R-help@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. > > >> > > >> > > >> > > >> > > >> > > >> > > >> -- > > >> > > >> St
Re: [R] Quantile Regression without intercept
> On Oct 6, 2015, at 8:32 AM, Lorenz, David <lor...@usgs.gov> wrote: > > Thanks for the details, I suspected something like that. > I think that begs the question: what is the meaning of quantile regression > through the origin? If the tau=.5 line does not pass through 1/2 the data how > do I interpret the line? As an estimate of the conditional median (quantile) function when constrained to pass through the origin… as with least squares fitting without an intercept, you do this at your peril. > > > On Tue, Oct 6, 2015 at 8:03 AM, Roger Koenker <rkoen...@illinois.edu> wrote: > > > On Oct 6, 2015, at 7:58 AM, Lorenz, David <lor...@usgs.gov> wrote: > > > > Did you verify that the correct percentages were above/below the regression > > lines? I did a quick check and for example did not consistently get 50% of > > the observed response values greater than the tau=.5 line. I did when I > > included the nonzero intercept term. > > Your "correct percentages" are only correct when you have an intercept in the > model, > without an intercept there is no gradient condition to ensure that. > > > > > > > >> Date: Mon, 5 Oct 2015 21:14:04 +0530 > >> From: Preetam Pal <lordpree...@gmail.com> > >> To: stephen sefick <ssef...@gmail.com> > >> Cc: "r-help@r-project.org" <r-help@r-project.org> > >> Subject: Re: [R] Quantile Regression without intercept > >> Message-ID: <56129a41.025f440a.b1cf4.f...@mx.google.com> > >> Content-Type: text/plain; charset="UTF-8" > >> > >> Yes..it works. Thanks ?? > >> > >> -Original Message- > >> From: "stephen sefick" <ssef...@gmail.com> > >> Sent: ?05-?10-?2015 09:01 PM > >> To: "Preetam Pal" <lordpree...@gmail.com> > >> Cc: "r-help@r-project.org" <r-help@r-project.org> > >> Subject: Re: [R] Quantile Regression without intercept > >> > >> I have never used this, but does the formula interface work like lm? Y~X-1? > >> > >> > >> On Mon, Oct 5, 2015 at 10:27 AM, Preetam Pal <lordpree...@gmail.com> > >> wrote: > >> > >> Hi guys, > >> > >> Can you instruct me please how to run quantile regression without the > >> intercept term? I only know about the rq function under quantreg package, > >> but it automatically uses an intercept model. Icant change that, it seems. > >> > >> I have numeric data on Y variable (Gdp) and 2 X variables (Hpa and > >> Unemployment). Their sizes are 125 each. > >> > >> Appreciate your help with this. > >> > >> Regards, > >> Preetam > >>[[alternative HTML version deleted]] > >> > >> __ > >> R-help@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. > >> > >> > >> > >> > >> > >> > >> -- > >> > >> Stephen Sefick > >> ** > >> Auburn University > >> Biological Sciences > >> 331 Funchess Hall > >> Auburn, Alabama > >> 36849 > >> ** > >> sas0...@auburn.edu > >> http://www.auburn.edu/~sas0025 > >> ** > >> > >> Let's not spend our time and resources thinking about things that are so > >> little or so large that all they really do for us is puff us up and make us > >> feel like gods. We are mammals, and have not exhausted the annoying little > >> problems of being mammals. > >> > >>-K. Mullis > >> > >> "A big computer, a complex algorithm and a long time does not equal > >> science." > >> > >> -Robert Gentleman > >>[[alternative HTML version deleted]] > >> > >> > >> > >> > > > > [[alternative HTML version deleted]] > > > > __ > > R-help@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. > > __ R-help@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.
Re: [R] Quantile Regression without intercept
> On Oct 6, 2015, at 7:58 AM, Lorenz, Davidwrote: > > Did you verify that the correct percentages were above/below the regression > lines? I did a quick check and for example did not consistently get 50% of > the observed response values greater than the tau=.5 line. I did when I > included the nonzero intercept term. Your "correct percentages" are only correct when you have an intercept in the model, without an intercept there is no gradient condition to ensure that. > > > >> Date: Mon, 5 Oct 2015 21:14:04 +0530 >> From: Preetam Pal >> To: stephen sefick >> Cc: "r-help@r-project.org" >> Subject: Re: [R] Quantile Regression without intercept >> Message-ID: <56129a41.025f440a.b1cf4.f...@mx.google.com> >> Content-Type: text/plain; charset="UTF-8" >> >> Yes..it works. Thanks ?? >> >> -Original Message- >> From: "stephen sefick" >> Sent: ?05-?10-?2015 09:01 PM >> To: "Preetam Pal" >> Cc: "r-help@r-project.org" >> Subject: Re: [R] Quantile Regression without intercept >> >> I have never used this, but does the formula interface work like lm? Y~X-1? >> >> >> On Mon, Oct 5, 2015 at 10:27 AM, Preetam Pal >> wrote: >> >> Hi guys, >> >> Can you instruct me please how to run quantile regression without the >> intercept term? I only know about the rq function under quantreg package, >> but it automatically uses an intercept model. Icant change that, it seems. >> >> I have numeric data on Y variable (Gdp) and 2 X variables (Hpa and >> Unemployment). Their sizes are 125 each. >> >> Appreciate your help with this. >> >> Regards, >> Preetam >>[[alternative HTML version deleted]] >> >> __ >> R-help@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. >> >> >> >> >> >> >> -- >> >> Stephen Sefick >> ** >> Auburn University >> Biological Sciences >> 331 Funchess Hall >> Auburn, Alabama >> 36849 >> ** >> sas0...@auburn.edu >> http://www.auburn.edu/~sas0025 >> ** >> >> Let's not spend our time and resources thinking about things that are so >> little or so large that all they really do for us is puff us up and make us >> feel like gods. We are mammals, and have not exhausted the annoying little >> problems of being mammals. >> >>-K. Mullis >> >> "A big computer, a complex algorithm and a long time does not equal >> science." >> >> -Robert Gentleman >>[[alternative HTML version deleted]] >> >> >> >> > > [[alternative HTML version deleted]] > > __ > R-help@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. __ R-help@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.
Re: [R] Quantile Regression without intercept
as for lm() or any other linear model fitting…. rq( y ~ x - 1, … ) url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 > On Oct 5, 2015, at 10:27 AM, Preetam Pal <lordpree...@gmail.com> wrote: > > Hi guys, > > Can you instruct me please how to run quantile regression without the > intercept term? I only know about the rq function under quantreg package, but > it automatically uses an intercept model. Icant change that, it seems. > > I have numeric data on Y variable (Gdp) and 2 X variables (Hpa and > Unemployment). Their sizes are 125 each. > > Appreciate your help with this. > > Regards, > Preetam > [[alternative HTML version deleted]] > > __ > R-help@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. __ R-help@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.
Re: [R] Quantile regression model with nonparametric effect and interaction
The main effect trend seems rather dangerous, why not just estimate the f’s in a loop? url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jun 11, 2015, at 1:57 AM, Waltl, Sofie (sofie.wa...@uni-graz.at) sofie.wa...@uni-graz.at wrote: Dear all, I would like to estimate a quantile regression model including a bivariate nonparametric term which should be interacted with a dummy variable, i.e., log p ~ year + f(a,b):year. I tried to use Roger Koenker's quantreg package and the functions rqss and qss but it turns out that interactions are not possible in this setting. Also weights are not implemented yet to build a work-around. I am looking for something like the by-statement in Simon Wood's package mgcv. Does anything comparable exist? I am grateful for any help on this issue. Kind regards, S. Waltl [[alternative HTML version deleted]] __ R-help@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. __ R-help@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.
Re: [R] Making my own graphics device
In ancient times, ie circa 1981, the S language certainly supported HP pen plotters so there should be code somewhere that could be resuscitated, he said naively. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Aug 17, 2014, at 2:58 PM, Thomas Levine _...@thomaslevine.com wrote: I want to make my own graphics device am thus looking for documentation about graphics devices. The only thing I've found so far is these directions for making graphics devices with the RGraphicsDevice package. http://www.omegahat.org/RGraphicsDevice/ Could someone point me to any other resources? Or just some documentation about how to edit base R? If I don't get anything, I'm just going to stare at the grDevices section of the R source code (src/library/grDevices/src) until I figure out how it works. In case you're curious, I want to make a graphics device that saves the graph in Hewlett-Packard Graphics Language. https://en.wikipedia.org/wiki/HPGL Thanks Tom __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Prediction intervals (i.e. not CI of the fit) for monotonic loess curve using bootstrapping
To follow up on David's suggestion on this thread, I might add that the demo(predemo) in my quantreg package illustrates a variety of approaches to prediction intervals for quantile regression estimates. Adapting this to monotone nonparametric estimation using rqss() or cobs would be quite straightforward, although the theory for such bands is rather difficult and still under construction. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Aug 13, 2014, at 9:59 AM, Jan Stanstrup jan.stanst...@fmach.it wrote: Thanks to all of you for your suggestions and comments. I really appreciate it. Some comments to Dennis' comments: 1) I am not concerned about predicting outside the original range. That would be nonsense anyway considering the physical phenomenon I am modeling. I am, however, concerned that the bootstrapping leads to extremely wide CIs at the extremes of the range when there are few data points. But I guess there is not much I can do about that as long as I rely on bootstrapping? 2) I have made a function that does the interpolation to the requested new x's from the original modeling data to get the residual variance and the model variance. Then it interpolates the combined SDs back the the new x values. See below. 3) I understand that. For this project it is not that important that the final prediction intervals are super accurate. But I need to hit the ballpark. I am only trying to do something that doesn't crossly underestimate the prediction error and doesn't make statisticians loose their lunch a first glance. I also cannot avoid that my data contains erroneous values and I will need to build many models unsupervised. But the fit should be good enough that I plan to eliminate values outside some multiple of the prediction interval and then re-calculate. And if the model is not good in any range I will throw it out completely. Based on the formula of my last message I have made a function that at least gives less optimistic intervals than what I could get with the other methods I have tried. The function and example data can be found here https://github.com/stanstrup/retpred_shiny/blob/master/retdb_admin/make_predictions_CI_tests.R in case anymore has any comments, suggestions or expletives to my implementation. -- Jan Stanstrup Postdoc Metabolomics Food Quality and Nutrition Fondazione Edmund Mach On 08/12/2014 05:40 PM, Bert Gunter wrote: PI's of what? -- future individual values or mean values? I assume quantreg provides quantiles for the latter, not the former. (See ?predict.lm for a terse explanation of the difference). Both are obtainable from bootstrapping but the details depend on what you are prepared to assume. Consult references or your local statistician for help if needed. -- Bert Bert Gunter Genentech Nonclinical Biostatistics (650) 467-7374 Data is not information. Information is not knowledge. And knowledge is certainly not wisdom. Clifford Stoll On Tue, Aug 12, 2014 at 8:20 AM, David Winsemius dwinsem...@comcast.net wrote: On Aug 12, 2014, at 12:23 AM, Jan Stanstrup wrote: Hi, I am trying to find a way to estimate prediction intervals (PI) for a monotonic loess curve using bootstrapping. At the moment my approach is to use the boot function from the boot package to bootstrap my loess model, which consist of loess + monoproc from the monoproc package (to force the fit to be monotonic which gives me much improved results with my particular data). The output from the monoproc package is simply the fitted y values at each x-value. I then use boot.ci (again from the boot package) to get confidence intervals. The problem is that this gives me confidence intervals (CI) for the fit (is there a proper way to specify this?) and not a prediction interval. The interval is thus way too optimistic to give me an idea of the confidence interval of a predicted value. For linear models predict.lm can give PI instead of CI by setting interval = prediction. Further discussion of that here: http://stats.stackexchange.com/questions/82603/understanding-the-confidence-band-from-a-polynomial-regression http://stats.stackexchange.com/questions/44860/how-to-prediction-intervals-for-linear-regression-via-bootstrapping. However I don't see a way to do that for boot.ci. Does there exist a way to get PIs after bootstrapping? If some sample code is required I am more than happy to supply it but I thought the question was general enough to be understandable without it. Why not use the quantreg package to estimate the quantiles of interest to you? That way you would not be depending on Normal theory assumptions which you apparently
Re: [R] goodness of fit for nonlinear quantile regression
Elsa, It is usual to write to package authors/maintainers about packages before trying R-help. Inference for nonlinear quantile regression is still a bit underdeveloped. To get logLik you can use: logLik.nlrq - function (object, ...) { n - length(object$m$resid()) p - length(object$m$getPars()) tau - object$m$tau() fid - object$m$rho val - n * (log(tau * (1 - tau)) - 1 - log(fid/n)) attr(val, n) - n attr(val, df) - p class(val) - logLik val } and from there AIC should work as usual for nls models. RK url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Oct 8, 2013, at 8:36 AM, Elsa Youngsteadt wrote: Hello, I am having trouble obtaining AIC or pseudo R^2 for nonlinear quantile regression fits. I would like to use one of these measures to compare non-nested models (each with a single, unique predictor variable). I am trying to fit the following gaussian peak model using the quantreg package: fit1.nlrq - nlrq(y ~ a*exp(-((x-b)/c)**2), data=data, start = list(a=.2,b=25.5,c=1), tau=0.5, trace=T); (and so on, for multiple tau; I would like a local goodness of fit measure at each tau, to help compare this model to a similar one that uses, say, x1 as a predictor instead of x.) Parameter estimates and predictions for the model looks as expected, but when I try to use AIC(fit1.nlrq) or AIC.nlrq(fit1.nlrq) I get the following output numeric(0) attr(,edf) [1] 0 Similarly, logLik(fit1.nlrq) yields 'log Lik.' (df=0) Can someone advise? (The output for deviance does exist and is a number...) As an alternative, I could perhaps calculate pseudo R2. I think this would be [1 - (deviance of fitted model/ deviance of null model)] but don't know how to obtain the deviance for the relevant null models. Can someone offer code or advice on this front? Perhaps someone has a more desirable alternative altogether for comparing among non-nested, nonlinear, quantile regression models. I am new to R and to R-help, please advise of posting mistakes or missing information! sessionInfo is as follows: R version 3.0.0 (2013-04-03) Platform: x86_64-w64-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C [5] LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] quantreg_4.98 SparseM_1.02 loaded via a namespace (and not attached): [1] tools_3.0.0 -- Elsa Youngsteadt Research Associate Department of Entomology North Carolina State University Campus Box 7613 Raleigh, NC 27695 919-515-1661 __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Matrix Multiplication using R.
In the event that these are moderately sparse matrices, you could try Matrix or SparseM. Roger Koenker rkoen...@illinois.edu On Aug 14, 2013, at 10:40 AM, Praveen Surendran wrote: Dear all, I am exploring ways to perform multiplication of a 9 x 4 matrix with it's transpose. As expected even a 4 x 100 %*% 100x4 didn't work on my desktop... giving the error Error: cannot allocate vector of length 16 However I am trying to run this on one node (64GB RAM; 2.60 GHz processor) of a high performance computing cluster. Appreciate if anyone has any comments on whether it's advisable to perform a matrix multiplication of this size using R and also on any better ways to handle this task. Kind Regards, Praveen. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Quantile regression for binary choice and heckit
This is a bit like asking how should I tweak my sailboat so I can explore the ocean floor. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On May 29, 2013, at 1:32 AM, Michal Kvasnička wrote: Hallo. Is there any package / code snippet to estimate quantile regression for a binary choice model (like probit) and selection model (like heckit)? I found that quantreg package can estimate tobit-like model, but I can't figure out how to tweak it for probit / heckit. Best wishes, Michal __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] updating observations in lm
The essential trick here is the Sherman-Morrison-Woodbury formula. My quantreg package has a lm.fit.recursive function that implements a fortran version for adding observations, but like biglm I don't remove observations at the other end either. Roger Koenker rkoen...@illinois.edu On May 27, 2013, at 2:07 PM, Greg Snow wrote: Look at the biglm package. It does 2 of the 3 things that you asked for: Construct an initial lm fit and add a new block of data to update that fit. It does not remove data, but you may be able to look at the code and figure out a way to modify it to do the final piece. On Mon, May 27, 2013 at 9:12 AM, ivo welch ivo.we...@anderson.ucla.eduwrote: dear R experts---I would like to update OLS regressions with new observations on the front of the data, and delete some old observations from the rear. my goal is to have a flexible moving-window regression, with a minimum number of observations and a maximum number of observations. I can keep (X' X) and (X' y), and add or subtract observations from these two quantities myself, and then use crossprod. strucchange does recursive residuals, which is closely related, but it is not designed for such flexible movable windows, nor primarily designed to produce standard errors of coefficients. before I get started on this, I just wanted to inquire whether someone has already written such a function. regards, /iaw Ivo Welch (ivo.we...@gmail.com) __ R-help@r-project.org mailing list 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. -- Gregory (Greg) L. Snow Ph.D. 538...@gmail.com [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] inverse for formula transformations on LHS
Paul, Inverting log(y) is just the beginning of the problem, after that you need to teach predict.lm() that E(y |x) = exp(x'betahat + .5*sigmahat^2) and then further lessons are required to get it to understand how to adapt its confidence and prediction bands… and then you need to generalize all of this to other transformations. Quite a project! Best, Roger Roger Koenker rkoen...@illinois.edu On May 17, 2013, at 12:21 PM, Paul Johnson wrote: This is an R formula handling question. It arose in class. We were working on the Animals data in the MASS package. In order to see a relationship, you need to log brain and body weight. It's a fun one for teaching regression, if you did not try it yet. There are outliers too! Students wanted to make a predicted value plot in the non-logged values of y, for comparison, and I wondered if I couldn't automate this somehow for them. It made me wonder how R manages formulae and if a transformation like log(y) can be be mechanically inverted. So we have something concrete to talk about, suppose x and y are variables in dat, a person fits m1 - lm(log(y) ~ log(x), data = dat) termplot shows log(y) on the vertical. What if I want y on the vertical? Similarly, predict gives values on the log(y) scale, there's no argument like type = untransformed. I want my solution to be a bit general, so that it would give back predicted y for formulae like sqrt(y) or exp(y) or log(y + d) or whatever other math people might throw in there. Here's what I can tell so far about R's insides. The formula handler makes a list out of the formula, I can get that from the terms object that the model generates. The formula list has ~ as element 1, and log(x) becomes element [[2]]. Where in the R source code can I see how R looks at the symbol log(y) and discerns that there is a variable y that needs to be logged? If I could understand that, and if R has a table of inverse functions, then maybe I could see what to do. If you have ideas, I'm very grateful if you share them. pj -- Paul E. Johnson Professor, Political Science Assoc. Director 1541 Lilac Lane, Room 504 Center for Research Methods University of Kansas University of Kansas http://pj.freefaculty.org http://quant.ku.edu [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] question
good usually means good relative to something else, in this case the comparison seems, as Michael has already said, f0 - rq(y ~ 1, tau = ?) and then one can compute the R1 version that I originally suggested. But since there is still no explicit way to evaluate this, it is all a bit pointless. Roger Koenker rkoen...@illinois.edu On Apr 24, 2013, at 6:37 PM, R. Michael Weylandt wrote: On Tue, Apr 23, 2013 at 2:54 PM, nafiseh hagiaghamohammadi n_hajiaghamohammadi2...@yahoo.com wrote: Hi I fit one linear quantile regression with package quantreg and I want to khow this model is good or not.Is there method for checking it? Thanks your advice I ask this question because there is 2 models,f0 and f1 in (R1 - 1 - f1$rho/f0$rho ), is it true? but I fit 1 model and I want to check goodness of fit for 1 model . Please keep your responses on list so you can get a quick reply even when I'm otherwise busy. I think you could -- for a rough and ready comparison -- compare against a constant (empirical quantile) model (not unlike how basic OLS models compare against the constant mean predictor) but someone else might know if there's any subtleties about quantile regression that should be noted here. MW __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] question
In quantreg if you do FAQ() item 4 is: 4. [R^2] I am currently trying to caculate the coefficient of determination for different quantile regression models. For example, how do you calculate the the sum of the weighted absolute deviations in the models ... R-squared is evil, that is why there isn't an automated way to compute something similar for quantile regression in the quantreg package. But if you insist use: R1 - 1 - f1$rho/f0$rho Provided that f0 is nested within f1 and the taus are the same, 0 = R1 = 1. If you want to test f1 vs f0 then use anova(f1,f0) For further details see: Koenker R. and Jose A.F. Machado. Goodness of Fit and Related Inference Processes for Quantile Regression J. of Am Stat. Assoc, (1999), 94, 1296-1310. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Apr 22, 2013, at 2:59 AM, nafiseh hagiaghamohammadi wrote: Hi Does anyone know if there is a method to calculate a goodness-of-fit statistic for quantile regressions with package quantreg? Tanks [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Singular design matrix in rq
Jonathan, This is not what we call a reproducible example... what is raw_data? Does it have something to do with mydata? what is i? Roger url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Apr 16, 2013, at 2:58 PM, Greenberg, Jonathan wrote: Quantreggers: I'm trying to run rq() on a dataset I posted at: https://docs.google.com/file/d/0B8Kij67bij_ASUpfcmJ4LTFEUUk/edit?usp=sharing (it's a 1500kb csv file named singular.csv) and am getting the following error: mydata - read.csv(singular.csv) fit_spl - rq(raw_data[,1] ~ bs(raw_data[,i],df=15),tau=1) Error in rq.fit.br(x, y, tau = tau, ...) : Singular design matrix Any ideas what might be causing this or, more importantly, suggestions for how to solve this? I'm just trying to fit a smoothed hull to the top of the data cloud (hence the large df). Thanks! --jonathan -- Jonathan A. Greenberg, PhD Assistant Professor Global Environmental Analysis and Remote Sensing (GEARS) Laboratory Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 607 South Mathews Avenue, MC 150 Urbana, IL 61801 Phone: 217-300-1924 http://www.geog.illinois.edu/~jgrn/ AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007 __ R-help@r-project.org mailing list 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.
Re: [R] bootstrapping quantile regression
There is no automatic clustering option for QR bootstrapping. You will have to roll your own. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Oct 31, 2012, at 1:38 AM, Kay Cichini wrote: sry, I forgot to replace rlm() - but actually I tried both and the question applies to both approaches.. Am 31.10.2012 00:19 schrieb Kay Cichini kay.cich...@gmail.com: HI everyone, I try to get some bootstrap CIs for coefficients obtained by quantile regression. I have influencial values and thus switched to quantreg.. The data is clustered and within clusters the variance of my DV = 0.. Is this sensible for the below data? And what about the warnings? Thanks in advance for any guidance, Kay dput(d) structure(list(Porenfläche = c(4990L, 7002L, 7558L, 7352L, 7943L, 7979L, 9333L, 8209L, 8393L, 6425L, 9364L, 8624L, 10651L, 8868L, 9417L, 8874L, 10962L, 10743L, 11878L, 9867L, 7838L, 11876L, 12212L, 8233L, 6360L, 4193L, 7416L, 5246L, 6509L, 4895L, 6775L, 7894L, 5980L, 5318L, 7392L, 7894L, 3469L, 1468L, 3524L, 5267L, 5048L, 1016L, 5605L, 8793L, 3475L, 1651L, 5514L, 9718L), P.Perimeter = c(2791.9, 3892.6, 3930.66, 3869.32, 3948.54, 4010.15, 4345.75, 4344.75, 3682.04, 3098.65, 4480.05, 3986.24, 4036.54, 3518.04, 3999.37, 3629.07, 4608.66, 4787.62, 4864.22, 4479.41, 3428.74, 4353.14, 4697.65, 3518.44, 1977.39, 1379.35, 1916.24, 1585.42, 1851.21, 1239.66, 1728.14, 1461.06, 1426.76, 990.388, 1350.76, 1461.06, 1376.7, 476.322, 1189.46, 1644.96, 941.543, 308.642, 1145.69, 2280.49, 1174.11, 597.808, 1455.88, 1485.58), P.Form = c(0.0903296, 0.148622, 0.183312, 0.117063, 0.122417, 0.167045, 0.189651, 0.164127, 0.203654, 0.162394, 0.150944, 0.148141, 0.228595, 0.231623, 0.172567, 0.153481, 0.204314, 0.262727, 0.200071, 0.14481, 0.113852, 0.291029, 0.240077, 0.161865, 0.280887, 0.179455, 0.191802, 0.133083, 0.225214, 0.341273, 0.311646, 0.276016, 0.197653, 0.326635, 0.154192, 0.276016, 0.176969, 0.438712, 0.163586, 0.253832, 0.328641, 0.230081, 0.464125, 0.420477, 0.200744, 0.262651, 0.182453, 0.200447), Durchlässigkeit = c(6.3, 6.3, 6.3, 6.3, 17.1, 17.1, 17.1, 17.1, 119, 119, 119, 119, 82.4, 82.4, 82.4, 82.4, 58.6, 58.6, 58.6, 58.6, 142, 142, 142, 142, 740, 740, 740, 740, 890, 890, 890, 890, 950, 950, 950, 950, 100, 100, 100, 100, 1300, 1300, 1300, 1300, 580, 580, 580, 580), Gebiete = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 6L, 6L, 6L, 6L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 3L, 7L, 7L, 7L, 7L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 5L, 5L, 5L, 5L, 12L, 12L, 12L, 12L, 8L, 8L, 8L, 8L), .Label = c(6.3, 17.1, 58.6, 82.4, 100, 119, 142, 580, 740, 890, 950, 1300), class = factor)), .Names = c(Porenfläche, P.Perimeter, P.Form, Durchlässigkeit, Gebiete), row.names = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48), class = data.frame) ## do quantile regression and bootstrap the coefficients, allowing for clustered data ## by putting Gebiet as strata argument (?), ## dv variation within clusters/Gebiet = 0! bs - function(formula, data, indices) { d - data[indices, ] # allows boot to select sample fit - rlm(formula, data = d) return(coef(fit)) } results - boot(data = d, statistic = bs, strata = d$Gebiete, R = 199, formula = Durchlässigkeit ~ P.Perimeter + P.Form) # get 99% confidence intervals boot.ci(results, type=bca, index=1, conf = .99) # intercept boot.ci(results, type=bca, index=2, conf = .99) # P.Perimeter boot.ci(results, type=bca, index=3, conf = .99) # P.Form -- Kay Cichini, MSc Biol Grubenweg 22, 6071 Aldrans E-Mail: kay.cich...@gmail.com -- [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] standard error for quantile
The rank test inversion option that you are trying to use won't work with only one coefficient, and therefore with univariate quantiles, if you use summary(rq(rnorm(50) ~ 1, tau = .9), cov = TRUE) you will have better luck. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Oct 30, 2012, at 9:42 AM, Koenker, Roger W wrote: Petr, You can do: require(quantreg) summary(rq(x ~ 1, tau = c(.10,.50,.99)) url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Oct 30, 2012, at 9:37 AM, Bert Gunter wrote: Petr: 1. Not an R question. 2. You want the distribution of order statistics. Search on that. It's basically binomial/beta. -- Bert On Tue, Oct 30, 2012 at 6:46 AM, PIKAL Petr petr.pi...@precheza.cz wrote: Dear all I have a question about quantiles standard error, partly practical partly theoretical. I know that x-rlnorm(10, log(200), log(2)) quantile(x, c(.10,.5,.99)) computes quantiles but I would like to know if there is any function to find standard error (or any dispersion measure) of these estimated values. And here is a theoretical one. I feel that when I compute median from given set of values it will have lower standard error then 0.1 quantile computed from the same set of values. Is it true? If yes can you point me to some reasoning? Thanks for all answers. Regards Petr PS. I found mcmcse package which shall compute the standard error but which I could not make to work probably because I do not have recent R-devel version installed Error in eval(expr, envir, enclos) : could not find function .getNamespace Error : unable to load R code in package 'mcmcse' Error: package/namespace load failed for 'mcmcse' Maybe I will also something find in quantreg package, but I did not went through it yet. __ R-help@r-project.org mailing list 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. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] standard error for quantile
Petr, You can do: require(quantreg) summary(rq(x ~ 1, tau = c(.10,.50,.99)) url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Oct 30, 2012, at 9:37 AM, Bert Gunter wrote: Petr: 1. Not an R question. 2. You want the distribution of order statistics. Search on that. It's basically binomial/beta. -- Bert On Tue, Oct 30, 2012 at 6:46 AM, PIKAL Petr petr.pi...@precheza.cz wrote: Dear all I have a question about quantiles standard error, partly practical partly theoretical. I know that x-rlnorm(10, log(200), log(2)) quantile(x, c(.10,.5,.99)) computes quantiles but I would like to know if there is any function to find standard error (or any dispersion measure) of these estimated values. And here is a theoretical one. I feel that when I compute median from given set of values it will have lower standard error then 0.1 quantile computed from the same set of values. Is it true? If yes can you point me to some reasoning? Thanks for all answers. Regards Petr PS. I found mcmcse package which shall compute the standard error but which I could not make to work probably because I do not have recent R-devel version installed Error in eval(expr, envir, enclos) : could not find function .getNamespace Error : unable to load R code in package 'mcmcse' Error: package/namespace load failed for 'mcmcse' Maybe I will also something find in quantreg package, but I did not went through it yet. __ R-help@r-project.org mailing list 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. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] SparseM buglet
Sam, Thanks for pointing this out, but I have to point out in turn that this isn't a SparseM function, it is part of the package e1071, maintained by David Meyer. Roger Roger Koenker rkoen...@illinois.edu On Aug 24, 2012, at 3:07 PM, Sam Steingold wrote: read.matrix.csr does not close the connection: library('SparseM') Package SparseM (0.96) loaded. read.matrix.csr(foo) ... Warning message: closing unused connection 3 (foo) -- Sam Steingold (http://sds.podval.org/) on Ubuntu 12.04 (precise) X 11.0.11103000 http://www.childpsy.net/ http://truepeace.org http://camera.org http://pmw.org.il http://think-israel.org http://dhimmi.com My other CAR is a CDR. __ R-help@r-project.org mailing list 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.
Re: [R] quantreg Wald-Test
On Jul 28, 2012, at 8:55 AM, stefan23 wrote: Dear all, I know that my question is somewhat special but I tried several times to solve the problems on my own but I am unfortunately not able to compute the following test statistic using the quantreg package. Well, here we go, I appreciate every little comment or help as I really do not know how to tell R what I want it to do^^ My situation is as follows: I have a data set containing a (dependent) vector Y and the regressor X. My aim is to check whether the two variables do not granger-cause each other in quantiles. I started to compute via quantreg for a single tau:= q: rq(Y_t~Y_(t-1)+Y_(t-2)+...+X_(t-1)+X_(t-2)+...,tau=q) This gives me the quantile regression coefficients. Now I want to check whether all the coefficients of X are equal to zero (for this specific tau). Can I do this by applying rq.anova ? I have already asked a similiar question but I am not sure if anova is really calculating this for me.. Currently I am calculating fitunrestricted=rq(Y_t~Y_(t-1)+Y_(t-2)+...+X_(t-1)+X_(t-2)+...,tau=q) fitrestrited=rq(Y_t~Y_(t-1)+Y_(t-2)+...,tau=q) anova(fitrestricted,fitunrestricted) If this is correct can you tell me how the test value is calculated in this case, Yes, this gives a Wald test if test = Wald is specified. You need to read the man page and the references mentioned there for further details. R-help isn't really intended to be a substitute for reading the literature. or in other words: My next step is going to replicate this procedure for a whole range of quantiles (say for tau in [a,b]). To apply a sup-Wald-test I am wondering if it is correct to choose the maximum of the different test values and to simulate the critical values by using the data tabulated in Andrees(1993) Do you mean Andrews? Such sup-Wald tests have been extensively discussed in various places, simulating asymptotic critical values is probably the most straightforward approach, but you need to be careful about tail behavior. (or simulate the vectors of independent Brownian Motions)...please feel free to comment, I am really looking forward to your help! Thank you very much Cheers Stefan -- View this message in context: http://r.789695.n4.nabble.com/quantreg-Wald-Test-tp4638198.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] The best solver for non-smooth functions?
There are obviously a large variety of non-smooth problems; for CVAR problems, if by this you mean conditional value at risk portfolio problems, you can use modern interior point linear programming methods. Further details are here: http://www.econ.uiuc.edu/~roger/research/risk/risk.html Roger Koenker rkoen...@illinois.edu On Jul 18, 2012, at 3:09 PM, Cren wrote: # Whoops! I have just seen there's a little mistake # in the 'sharpe' function, because I had to use # 'w' array instead of 'ead' in the cm.CVaR function! # This does not change the main features of my, # but you should be aware of it --- # The function to be minimized sharpe - function(w) { - (t(w) %*% y) / cm.CVaR(M, lgd, ead, N, n, r, rho, alpha, rating) } # This becomes... sharpe - function(w) { - (t(w) %*% y) / cm.CVaR(M, lgd, w, N, n, r, rho, alpha, rating) } # ...substituting 'ead' with 'w'. -- View this message in context: http://r.789695.n4.nabble.com/The-best-solver-for-non-smooth-functions-tp4636934p4636936.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Quantile Regression - Testing for Non-causalities in quantiles
Take a look at demo(Mel) in the quantreg package. Roger Koenker rkoen...@illinois.edu On Jul 14, 2012, at 6:55 AM, stefan23 wrote: Dear all, I am searching for a way to compute a test comparable to Chuang et al. (Causality in Quantiles and Dynamic Stock Return-Volume Relations). The aim of this test is to check wheter the coefficient of a quantile regression granger-causes Y in a quantile range. I have nearly computed everything but I am searching for an estimator of the density of the distribution at several points of the distribution. As the quantreg-package of Roger Koenker is also able to compute confidence intervalls for quantile regression (which also contain data concerning the estimated density) I wanted to ask wether someone could tell me if it is possible to extract the density of the underlying distribution by using the quantreg package. I hope my question is not to confusing, thank you very, very much in adavanve I appreciate every comment=) Cheers Stefan -- View this message in context: http://r.789695.n4.nabble.com/Quantile-Regression-Testing-for-Non-causalities-in-quantiles-tp4636511.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Quantile regression: Discrepencies Between optimizer and rq()
Optim() by default is using Nelder-Mead which is an extremely poor way to do linear programming, despite the fact that ?optim says that: It will work reasonably well for non-differentiable functions.I didn't check your coding of the objective function fully, but at the very least you should explicitly pass the arguments y, x, and quant. and you need to replace what you call sample.cond.quantile by 0 in the definition of w.less. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jun 7, 2012, at 1:49 PM, Kevin Chang wrote: Hello Everyone, I'm currently learning about quantile regressions. I've been using an optimizer to compare with the rq() command for quantile regression. When I run the code, the results show that my coefficients are consistent with rq(), but the intercept term can vary by a lot. I don't think my optimizer code is wrong and suspects it has something to do with the starting values. The results seems very sensitive to different starting values and I don't know how to make sense of it. Advice from the community would be greatly appreciated. Sincerely, Kevin Chang ## CODE Below ### library(quantreg) data(engel) y-cbind(engel[,2]) x-cbind(rep(1,length(engel[,1])),engel[,1]) x1-cbind(engel[,1]) nn-nrow(engel) nn bhat.ls-solve(t(x)%*%x)%*%t(x)%*%y #bhat.ls # QUANTILES quant=.25 fr.1=function(bhat.opt) { uu=y-x%*%bhat.opt sample.cond.quantile=quantile(uu,quant) w.less=rep(0,nn) for(ii in 1:nn){if(uu[ii]sample.cond.quantile) w.less[ii]=1} sum((quant-1)*sum((y-x%*%bhat.opt)*w.less) #negative residuals +quant*sum((y-x%*%bhat.opt)*(1-w.less))) #positive residuals } start-c(0,0) result=optim(start,fr.1) bhat.cond=result$par #Quantile Command Results fit.temp=rq(y~x1,tau=quant) fit.temp #OPTIMIZER Results bhat.cond #OLS Command Results mean=lm(y~x1) mean [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] R quantreg anova: How to change summary se-type
Stefan, You could try this: make a private version of anova.rqlist and change the call to lapply that computes summaries so that se = ker instead of se = nid. Please let me know if this does what you would like to do. This is about 20 lines into the code. Could you also explain what you mean by leads to mistakes below? Thanks, Roger url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On May 28, 2012, at 7:54 AM, stefan23 wrote: He folks=) I want to check whether a coefficient has an impact on a quantile regression (by applying the sup-wald test for a given quantile range [0.05,0.95]. Therefore I am doing the following calculations: a=0; for (i in 5:95/100){ fitrestricted=rq(Y~X1+X2,tau=i) tifunrestrited=rq(Y~X1+X2+X3,tau=i) a[i]=anova(fitrestricted,fitunrestricted)$table$Tn) #gives the Test-Value } supW=max(a) As anova is using the summary.rq function I want to change the Standard error method used (default: se=nid leads to mistakes, I prefer se=ker). Do you know how to handle this information in the anova syntax? Thank you very much Stefan -- View this message in context: http://r.789695.n4.nabble.com/R-quantreg-anova-How-to-change-summary-se-type-tp4631576.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] ANOVA in quantreg - faulty test for 'nesting'?
Yes, the models are nested, and yes I probably should have been more clever about parsing formulae, for these cases. I'll have a look at the code for glm, I presume that lm() is also capable of groking this? Meanwhile, I would like to point out that any of these linear restrictions can be reformulated as an exclusion restriction, see for example the remark on pages 5-6 of http://www.econ.uiuc.edu/~roger/research/ranks/ranktests.pdf or the references cited there. And using this sort of parameterization you can use anova.rq(). Roger Koenker rkoen...@illinois.edu On Apr 20, 2012, at 7:59 PM, Galen Sher wrote: Thanks to Brian's suggestion of using the logLik() function, I've dug a little deeper. I definitely think f1 and f2 are nested models. For example, by adding x2 to fmla1, we obtain a formula that quite clearly nests fmla1 and achieves the same log likelihood as that obtained for f2. Here is the extra code to show this: fmla3 = y~I(x1+x2)+x2 f3=glm(fmla3) logLik(f1); logLik(f2); logLik(f3) If f2=f3, as the log likelihoods would suggest, then this gives us a workaround: define the intermediate formula fmla3 and the fit f3 as above, and then conduct the analysis of variance between models f1 and f3 instead of f1 and f2. This doesn't offend anova.rq() any longer: f3.qr = rq(fmla3) anova(f1.qr,f3.qr) #this is actually anova(f1.qr, f2.qr) which resulted in an error earlier -Galen On Fri, Apr 20, 2012 at 6:47 PM, Brian S Cade ca...@usgs.gov wrote: Galen: Interesting, first time I tried it I couldn't get anova.glm to compute the p-values (got a warning) but I tried it again and it worked. Your larger (alternative hypothesis) model is y = B0 + B1*x1 + B2*x2 + e and your smaller (null hypotheisis) model is y = B0 + B3*(x1 + x2). So I guess I see that what you're trying to test in this case is that B1 = B2. I don't think Roger Koenker anticipated such a test with anova.rq. Other options besides using information criteria (AIC, BIC, etc) for comparing nonnested models include the Vuong test. But not sure how readily the theory of Vuong's test (like a paired t-test) extends to quantile regression. Brian Brian S. Cade, PhD U. S. Geological Survey Fort Collins Science Center 2150 Centre Ave., Bldg. C Fort Collins, CO 80526-8818 email: brian_c...@usgs.gov tel: 970 226-9326 From: Galen Sher galens...@gmail.com To: Brian S Cade ca...@usgs.gov Date: 04/20/2012 09:59 AM Subject: Re: [R] ANOVA in quantreg - faulty test for 'nesting'? -- Thanks Brian. I think anova.glm() requires the user to specify the appropriate distribution. In the example above, if I use either of the following commands anova(f1,f2,test=Chisq) #or anova(f1,f2,test=F) then anova.glm() will compute and display p-values associated with the deviance statistics. The reason I thought these models are nested is because the first model can be thought of as the second model estimated under an additional linear equality constraint. I suppose that's less of an R question and more of a regression question. Thanks for the logLik suggestion. In the absence of more information I'll have to do that - I'm just wary of conducting the test myself! Regards, Galen On Fri, Apr 20, 2012 at 4:31 PM, Brian S Cade *ca...@usgs.gov*ca...@usgs.gov wrote: Galen: Actually don't see how the models are nested (ask yourself what parameter in the model with 3 parameters is set to zero in the 2 parameter model?) and indeed if I try your code anova.glm will compute the difference in deviance between the two models but it does not compute a probability value for that difference in deviance as that would not make sense for models that aren't nested. Koenker's implementation of anova.rq immediately detects that the models aren't nested so doesn't even compute the deviance difference. You could use the logLik function on the rq objects to get their log likelihoods or use AIC (BIC) to compare the quantile regression models. Brian Brian S. Cade, PhD U. S. Geological Survey Fort Collins Science Center 2150 Centre Ave., Bldg. C Fort Collins, CO 80526-8818 email: *brian_c...@usgs.gov* brian_c...@usgs.gov tel: *970 226-9326* 970%20226-9326 From: galen *galens...@gmail.com* galens...@gmail.com To: * r-help@r-project.org* r-help@r-project.org Date: 04/19/2012 02:40 PM Subject: [R] ANOVA in quantreg - faulty test for 'nesting'? Sent by: * r-help-boun...@r-project.org* r-help-boun...@r-project.org -- I am trying to implement an ANOVA on a pair of quantile regression models in R. The anova.rq() function performs a basic check to see whether the models are nested, but I think this check is failing in my case. I think my models are nested despite the anova.rqlist() function saying otherwise. Here is an example where the GLM ANOVA regards the models as nested
Re: [R] non linear quantile regression - Median not plotting where it should
Dan, It is hard to say without being able to reproduce your example. If you send me the data I could try to advise something. Roger url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Feb 13, 2012, at 7:50 AM, Dan Morovitz wrote: Hi, I'm attempting to calculate the 0.25 and 0.97 quantiles for tree height (0-50 meters) against tree age (0-300 years) and I am running into some difficulty with the plotted grafic. I've run the examples in the quantreg help and can get those to work properly and by plugging in my data I can also get the lines plotted on my dataset. Unfortunately I'm running into a problem with the median and other regression lines with tree age younger than 50 years, basically the median is in this range overestimated and even comes out of the rage of oberservations. here is the code I'm using. # then fit the median using nlrq spruce.nlrq - nlrq(height ~ SSlogis(age, Asym, mid, scal),data=spruce, tau=0.5, trace=TRUE) lines(1:200, predict(spruce.nlrq, newdata=list(age=1:200)), col=2) I believe this has something to do with the SSlogis, as this gives the parameters for an S curve. My data set does not have the typical S curve shape, instead you could image it as starting at the inflection point of an S curve. This is what I expect the .025 quantile to be similar to: x - seq(1,100,1) plot(log(x)) the 0,975 should also have a logarithmic shape but no such a steep incline. Ive tried using different self starting models as found under: apropos(^ss) however I have not gotten them to work, or at least to fix the problem. Curves similar to mine look like these site index curves: http://www.extension.umn.edu/distribution/naturalresources/components/3473-13.html In the example from the nlrq help the lines of the median and the various quantiles all start from the same location, tha is basically (x=0, y=0) in the coodinate plane. With my problem, the lines start to be drawn from various different positions ( the lines always start at age(x)=0, but the height(y) can range between 0 and 15). Additionally, the data set is quite large. with about 50,000 oberservations on age and height. Does anyone have a suggestion on how to fix this problem? Thanks in advance Dan [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] about interpretation of anova results...
If you were to rtfm you'd see that your anova suggests that the slope coefficients are the same for your tau = .15 and tau = .30 models. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Dec 5, 2011, at 1:18 AM, narendarreddy kalam wrote: quantreg package is used. *fit1 results are* Call: rq(formula = op ~ inp1 + inp2 + inp3 + inp4 + inp5 + inp6 + inp7 + inp8 + inp9, tau = 0.15, data = wbc) Coefficients: (Intercept) inp1 inp2 inp3 inp4 inp5 -0.191528450 0.005276347 0.021414032 0.016034803 0.007510343 0.005276347 inp6 inp7 inp8 inp9 0.058708544 0.005224906 0.006804871 -0.003931540 Degrees of freedom: 673 total; 663 residual *fit2 results are* Call: rq(formula = op ~ inp1 + inp2 + inp3 + inp4 + inp5 + inp6 + inp7 + inp8 + inp9, tau = 0.3, data = wbc) Coefficients: (Intercept) inp1 inp2 inp3 inp4 -1.11e-01 5.776765e-19 4.635734e-18 1.874715e-18 2.099872e-18 inp5 inp6 inp7 inp8 inp9 -4.942052e-19 1.11e-01 2.205289e-18 4.138435e-18 9.300642e-19 Degrees of freedom: 673 total; 663 residual anova(fit1,fit2); Quantile Regression Analysis of Deviance Table Model: op ~ inp1 + inp2 + inp3 + inp4 + inp5 + inp6 + inp7 + inp8 + inp9 Joint Test of Equality of Slopes: tau in { 0.15 0.3 } Df Resid Df F value Pr(F) 1 9 1337 0.5256 0.8568 Warning messages: 1: In summary.rq(x, se = nid, covariance = TRUE) : 93 non-positive fis 2: In summary.rq(x, se = nid, covariance = TRUE) : 138 non-positive fis how to interpret the above results?? what is the use of anova function?? will it give the best among fit1 fit2.. -- View this message in context: http://r.789695.n4.nabble.com/about-interpretation-of-anova-results-tp4159510p4159510.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] linear against nonlinear alternatives - quantile regression
Roger Koenker rkoen...@illinois.edu On Nov 5, 2011, at 1:02 PM, Julia Lira wrote: Dear David, Indeed rq() accepts a vector fo tau. I used the example given by Frank to run fitspl4 - summary(rq(b1 ~ rcs(x,4), tau=c(a1,a2,a3,a4))) and it works. I even can use anova() to test equality of slopes jointly across quantiles. however, it would be interesting to test among different specifications, e.g. rcs(x,4) against rcs(x,3). but it does not work. Probably because the models aren't nested... Thanks for all suggestions! Julia From: dwinsem...@comcast.net Date: Sat, 5 Nov 2011 13:42:34 -0400 To: f.harr...@vanderbilt.edu CC: r-help@r-project.org Subject: Re: [R] linear against nonlinear alternatives - quantile regression I suppose this constitutes thread drift, but your simple example, Frank, made wonder if Rq() accepts a vector argument for tau. I seem to remember that Koencker's rq() does.. Normally I would consult the help page, but the power is still out here in Central Connecticut and I am corresponding with a less capable device. I am guessing that if Rq() does accept such a vector that the form of the nonlinearity would be imposed at all levels of tau. -- David On Nov 5, 2011, at 10:43 AM, Frank Harrell f.harr...@vanderbilt.edu wrote: Just to address a piece of this - in the case in which you are currently focusing on only one quantile, the rms package can help by fitting restricted cubic splines for covariate effects, and then run anova to test for nonlinearity (sometimes a dubious practice because if you then remove nonlinear terms you are mildly cheating). require(rms) f - Rq(y ~ x1 + rcs(x2,4), tau=.25) anova(f) # tests associations and nonlinearity of x2 Frank Julia Lira wrote: Dear all, I would like to know whether any specification test for linear against nonlinear model hypothesis has been implemented in R using the quantreg package. I could read papers concerning this issue, but they haven't been implemented at R. As far as I know, we only have two specification tests in this line: anova.rq and Khmaladze.test. The first one test equality and significance of the slopes across quantiles and the latter one test if the linear specification is model of location or location and scale shift. Do you have any suggestion? Thanks a lot! Best regards, Julia [[alternative HTML version deleted]] __ R-help@ mailing list 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. - Frank Harrell Department of Biostatistics, Vanderbilt University -- View this message in context: http://r.789695.n4.nabble.com/linear-against-nonlinear-alternatives-quantile-regression-tp3993327p3993416.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Question on estimating standard errors with noisy signals using the quantreg package
On Oct 31, 2011, at 7:30 AM, Thorsten Vogel wrote: Dear all, My question might be more of a statistics question than a question on R, although it's on how to apply the 'quantreg' package. Please accept my apologies if you believe I am strongly misusing this list. To be very brief, the problem is that I have data on only a random draw, not all of doctors' patients. I am interested in the, say, median number of patients of doctors. Does it suffice to use the nid option in summary.rq? More specifically, if the model generating the number of patients, say, r_i, of doctor i is r_i = const + u_i, then I think I would obtain the median of the number of doctors' patients using rq(r~1, ...) and plugging this into summary.rq() using the option se=iid. How big are the r_i? I presume that they are big enough so that you don't want to worry about the integer features of the data? Are there really no covariates? If so then you are fine with the iid option, but if not, probably better to use nid. If the r_i can be small, it is worth considering the dithering approach of Machado and Santos-Silva (JASA, 2005). Unfortunately, I don't observe r_i in the data but, instead, in the data I only have a fraction p of these r_i patients. In fact, with (known) probability p a patient is included in the data. Thus, for each doctor i the number of patients IN THE DATA follows a binomial distribution with parameters r_i and p. For each i I now have s_i patients in the data where s_i is a draw from this binomial distribution. That is, the problem with the data is that I don't observe r_i but s_i. Is it reasonable to assume that the p is the same across doctors? This seems to be some sort of compound Poisson problem to me, but I may misunderstand your description. Simple montecarlo experiments confirm my intuition that standard errors should be larger when using the noisy information s_i/p instead of (the unobserved) r_i. My guess is that I can consistently estimate any quantile of the number of doctors' patients AND THEIR STANDARD ERRORS using the quantreg's rq command: rq(I(s/p)~1, ...) and the summary.rq() command with option se=nid. Am I correct? I am greatful for any help on this issue. Best regards, Thorsten Vogel __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] nlrq {quantreg}
The model _is_ linear in parameters, after the log transformation of the response, so you don't need nlrq. If you really want something like: y = exp(a + b x) + u then you need to make a token effort to look at the documentation. Here is another example: x - exp(rnorm(50)) y - exp(1 + .5*x) + rnorm(50) nlrq(y ~ exp(a + b * x), start = list(a = 2, b = 1)) Nonlinear quantile regression model: y ~ exp(a + b * x) data: parent.frame tau: 0.5 deviance: 15.39633 a b 1.0348673 0.4962638 Roger Koenker rkoen...@illinois.edu On Oct 16, 2011, at 3:59 AM, Julia Lira wrote: Dear all, I sent an email on Friday asking about nlrq {quantreg}, but I haven't received any answer. I need to estimate the quantile regression estimators of a model as: y = exp(b0+x'b1+u). The model is nonlinear in parameters, although I can linearise it by using log.When I write: fitnl - nlrq(y ~ exp(x), tau=0.5) I have the following error: Error in match.call(func, call = cll) : invalid 'definition' argument Is there any way to estimate this model, or should I accept the following change: fitnl - rq(log(y) ~ x, tau=0.5) ? Thanks in advance! Best, Julia [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] How to keep a coefficient fixed when using rq {quantreg}?
Ah yes offsets, I've meant to look into this, but never quite understood why something like: rq((y - xk - bk) ~ x1 + x2) wasn't just as convenient Roger Koenker rkoen...@illinois.edu On Oct 14, 2011, at 6:55 PM, Tal Galili wrote: Hello all, I would like to compute a quantile regression using rq (from the quantreg package), while keeping one of the coefficients fixed. Is it possible to set an offset for rq in quantreg? (I wasn't able to make it to work) Thanks, Tal Contact Details:--- Contact me: tal.gal...@gmail.com | 972-52-7275845 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | www.r-statistics.com (English) __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] How to keep a coefficient fixed when using rq {quantreg}?
oops *bk not - bk. Roger Koenker rkoen...@illinois.edu On Oct 14, 2011, at 8:40 PM, Roger Koenker wrote: Ah yes offsets, I've meant to look into this, but never quite understood why something like: rq((y - xk - bk) ~ x1 + x2) wasn't just as convenient Roger Koenker rkoen...@illinois.edu On Oct 14, 2011, at 6:55 PM, Tal Galili wrote: Hello all, I would like to compute a quantile regression using rq (from the quantreg package), while keeping one of the coefficients fixed. Is it possible to set an offset for rq in quantreg? (I wasn't able to make it to work) Thanks, Tal Contact Details:--- Contact me: tal.gal...@gmail.com | 972-52-7275845 Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) | www.r-statistics.com (English) __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Change Variable Labels in Quantile Plot
It should be possible to say plot(..., xlab = foo, ylab = bar) ) if not please let me know. Roger Koenker rkoen...@illinois.edu On Aug 22, 2011, at 9:00 PM, Kitty Lee wrote: I have spent hours on this ---looked through the quantreg manual and r-help site--- still couldn't figure out the answer. Can someone please help me on this? I plot the result from quantile regression and want to change the variable labels: temp-rq(dep~inc+age50, data=newdata, tau=1:9/10) temp2-plot(summary(temp)) dimnames(temp2)[[1]]-c(Intercept, Per Capita Income, % Age 50 and Above) But after I manually change the dimnames, I can't replot the graph (temp2). Any suggestions? K. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] quantile regression: out of memory error
Paul, Yours is NOT a large problem, but it becomes a large problem when you ask for ALL the distinct QR solutions by specifying tau = -1. You probably don't want to see all these solutions, I suspect that only tau = 1:19/20 or so would suffice. Try this, and see how it goes. Roger url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jul 11, 2011, at 12:39 PM, Prew, Paul wrote: Hello, I’m wondering if anyone can offer advice on the out-of-memory error I’m getting. I’m using R2.12.2 on Windows XP, Platform: i386-pc-mingw32/i386 (32-bit). I am using the quantreg package, trying to perform a quantile regression on a dataframe that has 11,254 rows and 5 columns. object.size(subsetAudit.dat) 450832 bytes str(subsetAudit.dat) 'data.frame': 11253 obs. of 5 variables: $ Satisfaction : num 0.64 0.87 0.78 0.75 0.83 0.75 0.74 0.8 0.89 0.78 ... $ Return : num 0.84 0.92 0.91 0.89 0.95 0.81 0.9 0.87 0.95 0.88 ... $ Recommend: num 0.53 0.64 0.58 0.58 0.62 0.6 0.56 0.7 0.64 0.65 ... $ Cust.Clean : num 0.75 0.85 0.72 0.72 0.81 0.79 0.79 0.8 0.78 0.75 ... $ CleanScore.Cycle1: num 96.7 83.3 93.3 86.7 96.7 96.7 90 80 81.7 86.7 ... rq(subsetAudit.dat$Satisfaction ~ subsetAudit.dat$CleanScore.Cycle1, tau = -1) ERROR: cannot allocate vector of size 2.8 Gb I don’t know much about computers – software, hardware, algorithms – but does this mean that the quantreg package is creating some massive intermediate vector as it performs the rq function? Quantile regression is something I’m just starting to explore, but believe it involves ordering data prior to the regression, which could be a huge job when using 11,000 records. Does bigmemory have functionality to help me with this? Thank you, Paul Paul Prew ▪ Statistician 651-795-5942 ▪ fax 651-204-7504 Ecolab Research Center ▪ Mail Stop ESC-F4412-A 655 Lone Oak Drive ▪ Eagan, MN 55121-1560 CONFIDENTIALITY NOTICE: This e-mail communication and any attachments may contain proprietary and privileged information for the use of the designated recipients named above. Any unauthorized review, use, disclosure or distribution is prohibited. If you are not the intended recipient, please contact the sender by reply e-mail and destroy all copies of the original message. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] rqss help in Quantreg
This is _not_ a reproducible example so one can only guess, but the fact that with returns the analysis goes suggests that there is something amiss with your price variable. Roger Koenker rkoen...@illinois.edu On May 28, 2011, at 1:47 PM, Sergius Cerice wrote: Dear All, I,m trying to fulfill a constraint nonparametric quantile regression analysis for monthly stock index and gdp (159 cases of data) using rqss function of quantreg package. I need to specify that stock prices are nondecreasing with growing gdp. I tried the following simple code fit1-rqss(stock~gdp) fit2-rqss(stock~qss(gdp,constraint=I)+time) but R produces an error message for the firsts line of the code Error in rqss.fit(X, Y, tau = tau, rhs = rhs, nsubmax = nsubmax, nnzlmax = nnzlmax, : object 'rhs' is not found for the second line of the code Error in D %*% B : NA/NaN/Inf when calling external function (argument 7) If I use returns instead of prices, the analysis goes. But I need to regress prices. What is wrong in my specification? Are there any restrictions in the rqss approach? -- View this message in context: http://r.789695.n4.nabble.com/rqss-help-in-Quantreg-tp3392770p3557884.html Sent from the R help mailing list archive at Nabble.com. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] statistical question
The default rank test in the quantreg package would look like this summary(rq(y ~ d, tau = .5)) where d is a factor variable indicating which sample the elements of y belonged to. Summary returns a confidence interval for the coef of the factor variable -- if this interval excludes zero at the chosen alpha level then the difference in medians is significant. Roger Koenker rkoen...@illinois.edu On Mar 31, 2011, at 4:15 AM, Anna Lee wrote: Dear List! I want to compare medians of non normal distributed data. Is it possible and usefull to calculate 95% confidence intervals for medians? And if so - how can this be achieved in R? Thanks a lot! Anna -- Der Inhalt dieser E-Mail ist vertraulich. Sollte Ihnen die E-Mail irrtümlich zugesandt worden sein, bitte ich Sie, mich unverzüglich zu benachrichtigen und die E-Mail zu löschen. This e-mail is confidential. If you have received it in error, please notify me immediately and delete it from your system. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] rqss help in Quantreg
On Mar 20, 2011, at 10:24 PM, Man Zhang wrote: Dear All, I'm trying to construct confidence interval for an additive quantile regression model. In the quantreg package, vignettes section: Additive Models for Conditional Quantiles http://cran.r-project.org/web/packages/quantreg/index.html It describes how to construct the intervals, it gives the covariance matrix for the full set of parameters, \theta is given by the sandwich formula V = \tau (1-\tau) (t(X) \phi X)^(-1) (t(X) X)^(-1) (t(X) \phi X)^(-1) but gives no derivation of this result. Does anyone know how to obtain the above results? I need to construct these intervals with several adjustments, so I need to understand how V is constructed. Any proofs? As the author of the aforementioned vignette I suppose that I am meant to feel chastised by this request. I do not. There is a reasonable expectation in mathematics and most other fields that readers have some responsibility to familiarize themselves with the related literature, especially that cited in the work that they are currently reading. There seems to be a growing, deplorable expectation that learning a new subject is like being spoon fed some sort of pablum -- it isn't, sometimes you need to chew. I would remind would be R-posters of the admonition of Edward Davenant (quoted by Patrick Billingsley in his magisterial PM) I would have a man knockt on the head that should write anything in Mathematiques that had been written of before. Thanks, any help is greatly appreciated Man Zhang [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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.
Re: [R] Confidence interval on quantile regression predictions
You need to add explicitly newdata = list(x=x) and if you want percentile method you also need se = boot. Roger Koenker rkoen...@illinois.edu On Jan 11, 2011, at 3:54 AM, Davey, Andrew wrote: I am using the quantreg package to build a quantile regression model and wish to generate confidence intervals for the fitted values. After fitting the model, I have tried running predict() and predict.rq(), but in each case I obtain a vector of the fitted values only. For example: library(quantreg) y-rnorm(50,10,2) x-seq(1,50,1) model-rq(y~x,tau=0.5,method=br,model=TRUE) model$fitted predict(model,type=c(percentile),interval=c(confidence),level=0.95) #returns same output as model$fitted Am I missing something obvious? I have tried altering the settings for the type and interval arguments but without success. (I am running R v 2.10.1). Thank you for your help. Regards, Andrew Davey Senior Consultant Tel: + 44 (0) 1793 865023 Mob:+44 (0) 7747 863190 Fax: + 44 (0) 1793 865001 E-mail: andrew.da...@wrcplc.co.uk Website: www.wrcplc.co.uk P Before printing, please think about the environment. --- To read our latest newsletter visit http://www.wrcplc.co.uk/default.aspx?item=835 - Keeping our clients up-to-date with WRc's business developments --- Visit our websites www.wrcplc.co.uk and www.wrcnsf.com, as well as www.waterportfolio.com for collaborative research projects. --- The Information in this e-mail is confidential and is intended solely for the addressee. Access to this e-mail by any other party is unauthorised. If you are not the intended recipient, any disclosure, copying, distribution or any action taken in reliance on the information contained in this e-mail is prohibited and maybe unlawful. When addressed to WRc Clients, any opinions or advice contained in this e-mail are subject to the terms and conditions expressed in the governing WRc Client Agreement. --- WRc plc is a company registered in England and Wales. Registered office address: Frankland Road, Blagrove, Swindon, Wiltshire SN5 8YF. Company registration number 2262098. VAT number 527 1804 53. --- [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
[R] lty dots pdf issue
I'm trying to redo an old plot with: sessionInfo() R version 2.11.0 Under development (unstable) (2010-02-09 r51113) x86_64-apple-darwin9.8.0 When I do: pdf(lty.pdf,height = 6, width = 8) u - 1:100/100 y - matrix(rep(1:10,each = 100),100) matplot(u,y,lwd = 2,type =l) dev.off() the line types that have dots are difficult to distinguish because the dots that should appear are rendered as very thin vertical lines and it appears that the dashes are rendered without the lend round feature even though par() reports that that is the setting. I'm viewing all this with Acrobat 9 pro, but my printer produces something quite consistent with what is viewable on the screen. Apologies in advance if this has an obvious resolution, I didn't see anything in the r-help archive. Roger url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 __ R-help@r-project.org mailing list 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.
Re: [R] sparseM and kronecker product_R latest version
I can't reproduce this on my mac: sessionInfo() R version 2.10.1 (2009-12-14) x86_64-apple-darwin9.8.0 locale: [1] C attached base packages: [1] graphics grDevices datasets utils stats methods base other attached packages: [1] fortunes_1.3-7 quantreg_4.44 SparseM_0.83 loaded via a namespace (and not attached): [1] tools_2.10.1 url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jan 11, 2010, at 8:00 AM, David Winsemius wrote: On Jan 11, 2010, at 7:55 AM, Peter Ehlers wrote: Do you have the same problem with the example on the help page? ?'%x%-methods' Works for me on Windows Vista (32-bit OS) and R version 2.10.1 Patched (2010-01-05 r50896). -Peter Ehlers alessia matano wrote: Dear all, I just installed the new version of R, 2.10.1, and I am currently using the package sparseM. (I also use a 64 bit windows version) SparseM I got a problem that I never had: when I try to multiply with a kronecker product (%x%) two sparse matrixes I get the following message: Error in dim(x) - length(x) : invalid first argument I never had this problem with previous versions of R. I get the same error as you do when trying the example on the cited help page using SparseM 0.83 in a 64 bit Mac version of 2.10.1 A.csr %x% matrix(1:4,2,2) Error in dim(x) - length(x) : invalid first argument sessionInfo() R version 2.10.1 RC (2009-12-09 r50695) x86_64-apple-darwin9.8.0 locale: [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] grid splines stats graphics grDevices utils datasets methods base other attached packages: [1] SparseM_0.83 Matrix_0.999375-32 Epi_1.1.10 plotrix_2.7-2 ROCR_1.0-4 [6] gplots_2.7.4 caTools_1.10 bitops_1.0-4.1 gdata_2.6.1 gtools_2.6.1 [11] lattice_0.17-26Design_2.3-0 Hmisc_3.7-0survival_2.35-7 loaded via a namespace (and not attached): [1] cluster_1.12.1 tools_2.10.1 May you help me? thanks alessia __ -- Peter Ehlers University of Calgary 403.202.3921 __ David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list 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.
[R] Predicting the Present with R at Google
I thought some R-help readers might enjoy seeing the paper by Hal Varian on predicting the present using R and Google Trends that is linked via the following blog comment: http://googleresearch.blogspot.com/2009/04/predicting-present-with-google-trends.html Apologies in advance if this has already been mentioned, I didn't see it in the achives. -- View this message in context: http://www.nabble.com/Predicting-the-Present-with-R-at-Google-tp25513381p25513381.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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.
Re: [R] Browser and Debug?
At points of total desperation you can always consider the time-honored, avuncular advice -- RTFM, in this case Section 4.4 of Writing R Extensions. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Aug 13, 2009, at 10:20 AM, rkevinbur...@charter.net wrote: This may be asking more than can be reasonably expected. But, any tips on debugging the 'C' and Fortran code without trying to put together all the tools to compile from source? Kevin Erik Iverson eiver...@nmdp.org wrote: This article might help: http://www.biostat.jhsph.edu/~rpeng/docs/R-debug-tools.pdf -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org ] On Behalf Of Inchallah Yarab Sent: Thursday, August 13, 2009 9:40 AM To: r-help@r-project.org Subject: [R] Browser and Debug? Hi, Someone can explain to me how use Browser and Debug ? thank you [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Logit Model... GLM or GEE or ??
You could take a look at: M West, PJ Harrison, HS Migon - Journal of the American Statistical Association, 1985 - jstor.org Page 1. Dynamic Generalized Linear Models and Bayesian Forecasting and the subsequent literature it has generated... or along the same lines the literature on chess ratings. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Aug 6, 2009, at 5:00 PM, Noah Silverman wrote: Posted about this earlier. Didn't receive any response But, some further research leads me to believe that MAYBE a GLMM or a GEE function will do what I need. Hello, I have a bit of a tricky puzzle with trying to implement a logit model as described in a paper. The particular paper is on horseracing and they explain a model that is a logit trained per race, yet somehow the coefficients are combined across all the training races to come up with a final set of coefficients. My understanding is that they maximize log likelihood across the entire set of training races. Yet this isn't just as standard logit model as they are looking at data per race. This is a bit hard to explain, so I've attached a tiny pdf of the paragraph from the paper explaining this. Like everything else in the data/stat/econ world, there is probably a library in R that does this kind of thing, but after 3 days of heavy google research, I've been unable to find it. Does anyone have any suggestions?? Thanks. -N Attached is a jpg of the book page describing what I'm trying to do... __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
[R] Fwd: Making rq and bootcov play nice
John, You can make a local version of bootcov which either: deletes these arguments from the call to fitter, or modify the switch statement to include rq.fit, the latter would need to also modify rq() to return a fitFunction component, so the first option is simpler. One of these days I'll incorporate clustered se's into summary.rq, but meanwhile this seems to be a good alternative. Roger url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jul 23, 2009, at 8:10 PM, John Gardner wrote: I have a quick question, and I apologize in advance if, in asking, I expose my woeful ignorance of R and its packages. I am trying to use the bootcov function to estimate the standard errors for some regression quantiles using a cluster bootstrap. However, it seems that bootcov passes arguments that rq.fit doesn't like, preventing the command from executing. Here is an example: e-bootcov(rq(y~x),clust,B=10,fitter=rq.fit) (where clust is my clustering variable) results in Error in rq.fit.br(x, y, tau = tau, ...) : unused argument(s) (maxit = 15, penalty.matrix = NULL) In contrast, the lm.fit function seems to just ignore these arguments, resulting in the following warning: 10: In fitter(X[obs, , drop = FALSE], Y[obs, , drop = FALSE], maxit = maxit, : extra arguments maxitpenalty.matrix are just disregarded. Is there a way that I can either (a) modify bootcov so that it doesn't pass these arguments or (b) modify rq so that it ignores them? Thanks in advance, John Gardner __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] How to pass a character argument which contains expressions to arg.names in barplot?
Surely a fortune, because fortuna can be cruel. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jul 24, 2009, at 3:24 PM, jcano wrote: Dear Uwe Thank you very much for your unvaluable time and effort Cheers Javier Cano __ R-help@r-project.org mailing list 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.
[R] Constructing lists (yet, again)
This is an attempt to rescue an old R-help question that apparently received no response from the oblivion of collective silence, and besides I'm also curious about the answer From: Griffith Feeney (gfee...@hawaii.edu) Date: Fri 28 Jan 2000 - 07:48:45 EST wrote (to R-help) Constructing lists with list(name1=name1, name2=name2, ...) is tedious when there are many objects and names are long. Is there an R function that takes a character vector of object names as an argument and returns a list with each objected tagged by its name? The idiom lapply(ls(pat = ^name), function(x) eval(as.name(x))) makes the list, but (ironically) doesn't assign the names to the components. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 __ R-help@r-project.org mailing list 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.
Re: [R] Constructing lists (yet, again)
Marc, Thanks, sapply(ls(pat = ^name),get) was exactly what I was after. The default behavior for vectors of equal length is nice too, but I was most interested in the ragged case, which produces a list. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jul 23, 2009, at 2:41 PM, Marc Schwartz wrote: On Jul 23, 2009, at 9:19 AM, roger koenker wrote: This is an attempt to rescue an old R-help question that apparently received no response from the oblivion of collective silence, and besides I'm also curious about the answer From: Griffith Feeney (gfee...@hawaii.edu) Date: Fri 28 Jan 2000 - 07:48:45 EST wrote (to R-help) Constructing lists with list(name1=name1, name2=name2, ...) is tedious when there are many objects and names are long. Is there an R function that takes a character vector of object names as an argument and returns a list with each objected tagged by its name? The idiom lapply(ls(pat = ^name), function(x) eval(as.name(x))) makes the list, but (ironically) doesn't assign the names to the components. Roger, How about something like this: name1 - 1:3 name2 - 1:5 name3 - 1:9 name4 - 1:7 ls(pat = ^name) [1] name1 name2 name3 name4 sapply(ls(pat = ^name), get, simplify = FALSE) $name1 [1] 1 2 3 $name2 [1] 1 2 3 4 5 $name3 [1] 1 2 3 4 5 6 7 8 9 $name4 [1] 1 2 3 4 5 6 7 Is that what you are after? With sapply(), you can take advantage of the USE.NAMES argument, which defaults to TRUE and then set simplify to FALSE to force the result to be a list rather than a matrix. Of course, in the case I have above, when there are uneven length elements, the result will be a list anyway. HTH, Marc Schwartz __ R-help@r-project.org mailing list 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.
Re: [R] package quantreg behaviour in weights in function rq,
Václav, I think that the technical term for this is snafu. You are right of course that the magnitude of the weights should have no impact on the fitted values. The good news is that the coefficient estimates satisfy this obvious principle, but unfortunately the fitted.values component is being computed with the wrong form of the X matrix, so your plots are misaligned. The immediate work around is to compute fitted values from the original data and the coef component of the fit rather than relying on the fitted.values component. Thanks for pointing this out, I'll try to submit a repaired version of the package later in the day. Roger url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jul 21, 2009, at 7:03 AM, Václav Varvařovský wrote: Dear all, I am having v.4.36 of Quantreg package and I noticed strange behaviour when weights were added. Could anyone please explain me what if the results are really strange or the behavioiur is normal. As an example I am using dataset Engel from the package and my own weights. x-engel[1:50,1] y-engel[1:50,2] w-c(0.00123, 0.00050, 0.00126, 0.00183, 0.00036, 0.00100, 0.00122, 0.00133, 0.01208, 0.00126, 0.00102, 0.00183, 0.00063, 0.00134, 0.00084, 0.00087, 0.00118, 0.00894, 0.00105, 0.00154, 0.02829, 0.00095, 0.05943, 0.07003, 0.00692, 0.03610, 0.00316, 0.06862, 0.00439, 0.08974, 0.01960, 0.00185, 0.00348, 0.03597, 0.00210, 0.03929, 0.03535, 0.01463, 0.02254, 0.00089, 0.01495, 0.00178, 0.00351, 0.10338, 0.13662, 0.00157, 0.07689, 0.07304, 0.00194, 0.00142) windows(width = 8, height = 7) u1-rq(y~x,tau=c(1:10/10-0.05),weights=w) #note that by taking weights = 500*w the points #are all moving down (i thought it should have been invariant to the magnitude of weights) when all the weights remaint the same plot(x,y,ylim=c(-1000,max(y))) points(x,u1$fitted.values[,3],col=blue,cex=0.5) #weighted - fitted values nearly match original values windows(width = 8, height = 7) u1-rq(y~x,tau=c(1:10/10-0.05)) plot(x,y,ylim=c(-1000,max(y))) points(x,u1$fitted.values[,3],col=blue,cex=0.5) #unweighted - fitted values form a line Why the weighted quantile regression does not produce a line? Thank you for answers. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Matrix multiplication precision
A good discussion of this is provided by Gill, Murray and Wright Num Lin Alg and Opt, section 4.7.2. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jul 15, 2009, at 9:27 AM, Steve Lianoglou wrote: Hi Douglas, And the big lesson, of course, is the first rule of numerical linear algebra., Just because a formula is written in terms of the inverse of a matrix doesn't mean that is a good way to calculate the result; in fact, it is almost inevitably the worst way of calculating the result. You only calculate the inverse of a matrix after you have investigated all other possible avenues for calculating the result and found them to be fruitless. As a relative noob to the nitty gritty details of numerical linear algebra, could you elaborate on this a bit (even if it's just a link to a book/reference that explains these issues in more detail)? Where/what else should we be looking to do that would be better? Or are there really no general rules, and the answer just depends on the situation? Thanks, -steve -- Steve Lianoglou Graduate Student: Physiology, Biophysics and Systems Biology Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Fortran function for quantiles
You could look at the function kuantile() in the package quantreg and the associated fortran code dsel05.f which implements a version of the Floyd and Rivest (CACM, 1975) algorithm. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jul 13, 2009, at 6:24 PM, DHIMAN BHADRA wrote: Hi, I was wondering whether there is any Fortran function or associated library for evaluating the quantiles of a set of values (something which the R-function quantile() does). Any help will be much appreciated. Thanks and regards, Dhiman Bhadra [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Adaptive kernel density was ... Fw: (no subject)
On Jul 12, 2009, at 2:56 PM, David Winsemius wrote: On Jul 12, 2009, at 3:21 PM, maram salem wrote: Dear group, Thank u so much 4 ur help. I've tried the link, http://finzi.psych.upenn.edu/R/library/quantreg/html/akj.html for adaptive kernel density estimation. But since I'm an R beginer and the topic of adaptive estimation is new for me, i still can't figure out some of the arguments of akj(x, z =, p =, h = -1, alpha = 0.5, kappa = 0.9, iker1 = 0) I've a vector of 1000 values (my X), but I don't know how to get the Z That does seem rather trivial. According to the help page, those are just the points at which the density should be estimated. The example in the help page shows you how to create a suitable vector. and what's Kappa? Not so obvious. Experimentation shows that reducing kappa makes the estimates less smooth. I'm sorry if the question is trivial but I hope u could recommend some refrence if u know one. Koenker gives two references and apparently you have some other material you are reading. Your university should have access to the Project Euclid Annals of Statistics copies that are found with the obvious Google search strategy. Maybe you should be questioning the overall strategy of using a function you don't understand. Why, for instance, do you even have an interest in this function? The Silverman book on density estimation is as close to a canonical reference as one is likely to encounter in statistics, the akj() function implements Silverman's adaptive kernel proposal so it would be quite helpful to have this reference at hand. That is why it was cited in the man page for the function. __ R-help@r-project.org mailing list 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.
Re: [R] Dantzig Selector
There is an experimental version available here: http://www.econ.uiuc.edu/~roger/research/sparse/sfn.html that uses the interior point code in the package quantreg. There is an option to exploit possible sparsity of the X matrix. Comments would be welcome. url:www.econ.uiuc.edu/~rogerRoger Koenker email rkoen...@uiuc.edu Department of Economics vox:217-333-4558University of Illinois fax:217-244-6678Champaign, IL 61820 On Jul 8, 2009, at 6:35 PM, tzygmund mcfarlane wrote: Hi, I was wondering if there was an R package or routines for the Dantzig Selector (Candes Tao, 2007). I know Emmanuel Candes has Matlab routines to do this but I was wondering if someone had ported those to R. Thanks, T ---Reference--- @article{candes2007dantzig, title={{The Dantzig selector: statistical estimation when p is much larger than n}}, author={Candes, E. and Tao, T.}, journal={Annals of Statistics}, volume={35}, number={6}, pages={2313--2351}, year={2007}, publisher={Hayward, Calif.[etc] Institute of Mathematical Statistics [etc]} } __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Color of ecdf plots
having been bitten by this behavior many times, could I register a plea for allowing a col argument to plot.stepfun that would deal with both the horizontal and vertical segments -- I rather doubt that it is often desirable to have different colors for these. Note that verticals = TRUE is the default here so not strictly necessary in the example below. On Jul 3, 2009, at 12:34 PM, David Winsemius wrote: Does not work for me with rnorm(20) as an argument. (And supplying a dummy value for DIM.) The error after fixing the missing DIM issue was: Error in plot.stepfun(x, ..., ylab = ylab, verticals = verticals, pch = pch) : argument 4 matches multiple formal arguments The help page for plot.ecdf refers me to plot.stepfun (and the error message did too) which does not have a col= argument but rather col.hor= and col.vert= . Using those arguments instead of col=red, it now works. plot( ecdf(rnorm(20)), do.points=FALSE, verticals=TRUE, main=paste(Ecdf of distances ), col.hor=red, col.vert=red ) On Mac OS X 10.5.7, 64bit R 2.8.1 On Jul 3, 2009, at 12:24 PM, Uwe Ligges wrote: Lars Bergemann wrote: Hi. I have the following two ecdf plots in one graph: plot( ecdf(), do.points=FALSE, verticals=TRUE, main=paste(Ecdf of distances ,DIM,sep=), col=red ); lines( ecdf(), do.points=FALSE, verticals=TRUE ); How do I change the color of the resulting graph? Adding col=red to either plot or lines results in an error ... ... but works for me. Hence please specify your version of R, OS, exactly reproducible code (above the argument to ecdf are missing, obviously). Best, Uwe Ligges David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] odd behaviour in quantreg::rq
It's not clear to me whether you are looking for an exploratory tool or something more like formal inference. For the former, it seems that estimating a few weighted quantiles would be quite useful. at least it is rather Tukeyesque. While I'm appealing to authorities, I can't resist recalling that Galton's invention of correlation paper: Co- relations and their measurement, Proceedings of the Royal Society, 1888-89, used medians and interquartile ranges. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jul 1, 2009, at 2:48 PM, Dylan Beaudette wrote: Thanks Roger. Your comments were very helpful. Unfortunately, each of the 'groups' in this example are derived from the same set of data, two of which were subsets-- so it is not that unlikely that the weighted medians were the same in some cases. This all leads back to an operation attempting to answer the question: Of the 2 subsetting methods, which one produces a distribution most like the complete data set? Since the distributions are not normal, and there are area-weights involved others on the list suggested quantile- regression. For a more complete picture of how 'different' the distributions are, I have tried looking at the differences between weighted quantiles: (0.1, 0.25, 0.5, 0.75, 0.9) as a more complete 'description' of each distribution. I imagine that there may be a better way to perform this comparison... Cheers, Dylan On Tuesday 30 June 2009, roger koenker wrote: Admittedly this seemed quite peculiar but if you look at the entrails of the following code you will see that with the weights the first and second levels of your x$method variable have the same (weighted) median so the contrast that you are estimating SHOULD be zero. Perhaps there is something fishy about the data construction that would have allowed us to anticipate this? Regarding the fn option, and the non-uniqueness warning, this is covered in the (admittedly obscure) faq on quantile regression available at: http://www.econ.uiuc.edu/~roger/research/rq/FAQ # example: library(quantreg) # load data x - read.csv(url('http://169.237.35.250/~dylan/temp/test.csv')) # with weights summary(rq(sand ~ method, data=x, weights=area_fraction, tau=0.5), se='ker') #Reproduction with more convenient notation: X0 - model.matrix(~method, data = x) y - x$sand w - x$area_fraction f0 - summary(rq(y ~ X0 - 1, weights = w),se = ker) #Second reproduction with orthogonal design: X1 - model.matrix(~method - 1, data = x) f1 - summary(rq(y ~ X1 - 1, weights = w),se = ker) #Comparing f0 and f1 we see that they are consistent!! How can that be?? #Since the columns of X1 are orthogonal estimation of the 3 parameters are separable #so we can check to see whether the univariate weighted medians are reproducible. s1 - X1[,1] == 1 s2 - X1[,2] == 1 g1 - rq(y[s1] ~ X1[s1,1] - 1, weights = w[s1]) g2 - rq(y[s2] ~ X1[s2,2] - 1, weights = w[s2]) #Now looking at the g1 and g2 objects we see that they are equal and agree with f1. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jun 30, 2009, at 3:54 PM, Dylan Beaudette wrote: Hi, I am trying to use quantile regression to perform weighted- comparisons of the median across groups. This works most of the time, however I am seeing some odd output in summary(rq()): Call: rq(formula = sand ~ method, tau = 0.5, data = x, weights = area_fraction) Coefficients: ValueStd. Error t value Pr(|t|) (Intercept)45.44262 3.64706 12.46007 0.0 methodmukey-HRU 0.0 4.671150.0 1.0 ^ When I do not include the weights, I get something a little closer to a weighted comparison of means, along with an error message: Call: rq(formula = sand ~ method, tau = 0.5, data = x) Coefficients: ValueStd. Error t value Pr(|t|) (Intercept)44.91579 2.46341 18.23318 0.0 methodmukey-HRU 9.57601 9.293481.03040 0.30380 Warning message: In rq.fit.br(x, y, tau = tau, ...) : Solution may be nonunique I have noticed that the error message goes away when specifying method='fn' to rq(). An example is below. Could this have something to do with replication in the data? # example: library(quantreg) # load data x - read.csv(url('http://169.237.35.250/~dylan/temp/test.csv')) # with weights summary(rq(sand ~ method, data=x, weights=area_fraction, tau=0.5), se='ker') # without weights # note error message summary(rq(sand ~ method, data=x, tau=0.5), se='ker') # without weights, no error
Re: [R] odd behaviour in quantreg::rq
Admittedly this seemed quite peculiar but if you look at the entrails of the following code you will see that with the weights the first and second levels of your x$method variable have the same (weighted) median so the contrast that you are estimating SHOULD be zero. Perhaps there is something fishy about the data construction that would have allowed us to anticipate this? Regarding the fn option, and the non-uniqueness warning, this is covered in the (admittedly obscure) faq on quantile regression available at: http://www.econ.uiuc.edu/~roger/research/rq/FAQ # example: library(quantreg) # load data x - read.csv(url('http://169.237.35.250/~dylan/temp/test.csv')) # with weights summary(rq(sand ~ method, data=x, weights=area_fraction, tau=0.5), se='ker') #Reproduction with more convenient notation: X0 - model.matrix(~method, data = x) y - x$sand w - x$area_fraction f0 - summary(rq(y ~ X0 - 1, weights = w),se = ker) #Second reproduction with orthogonal design: X1 - model.matrix(~method - 1, data = x) f1 - summary(rq(y ~ X1 - 1, weights = w),se = ker) #Comparing f0 and f1 we see that they are consistent!! How can that be?? #Since the columns of X1 are orthogonal estimation of the 3 parameters are separable #so we can check to see whether the univariate weighted medians are reproducible. s1 - X1[,1] == 1 s2 - X1[,2] == 1 g1 - rq(y[s1] ~ X1[s1,1] - 1, weights = w[s1]) g2 - rq(y[s2] ~ X1[s2,2] - 1, weights = w[s2]) #Now looking at the g1 and g2 objects we see that they are equal and agree with f1. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jun 30, 2009, at 3:54 PM, Dylan Beaudette wrote: Hi, I am trying to use quantile regression to perform weighted- comparisons of the median across groups. This works most of the time, however I am seeing some odd output in summary(rq()): Call: rq(formula = sand ~ method, tau = 0.5, data = x, weights = area_fraction) Coefficients: ValueStd. Error t value Pr(|t|) (Intercept)45.44262 3.64706 12.46007 0.0 methodmukey-HRU 0.0 4.671150.0 1.0 ^ When I do not include the weights, I get something a little closer to a weighted comparison of means, along with an error message: Call: rq(formula = sand ~ method, tau = 0.5, data = x) Coefficients: ValueStd. Error t value Pr(|t|) (Intercept)44.91579 2.46341 18.23318 0.0 methodmukey-HRU 9.57601 9.293481.03040 0.30380 Warning message: In rq.fit.br(x, y, tau = tau, ...) : Solution may be nonunique I have noticed that the error message goes away when specifying method='fn' to rq(). An example is below. Could this have something to do with replication in the data? # example: library(quantreg) # load data x - read.csv(url('http://169.237.35.250/~dylan/temp/test.csv')) # with weights summary(rq(sand ~ method, data=x, weights=area_fraction, tau=0.5), se='ker') # without weights # note error message summary(rq(sand ~ method, data=x, tau=0.5), se='ker') # without weights, no error message summary(rq(sand ~ method, data=x, tau=0.5, method='fn'), se='ker') -- Dylan Beaudette Soil Resource Laboratory http://casoilresource.lawr.ucdavis.edu/ University of California at Davis 530.754.7341 __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Memory issues on a 64-bit debian system (quantreg)
Jonathan, Take a look at the output of sessionInfo(), it should say x86-64 if you have a 64bit installation, or at least I think this is the case. Regarding rqss(), my experience is that (usually) memory problems are due to the fact that early on the processing there is a call to model.matrix() which is supposed to create a design, aka X, matrix for the problem. This matrix is then coerced to matrix.csr sparse format, but the dense form is often too big for the machine to cope with. Ideally, someone would write an R version of model.matrix that would permit building the matrix in sparse form from the get-go, but this is a non-trivial task. (Or at least so it appeared to me when I looked into it a few years ago.) An option is to roll your own X matrix: take a smalller version of the data, apply the formula, look at the structure of X and then try to make a sparse version of the full X matrix. This is usually not that difficult, but usually is based on a rather small sample that may not be representative of your problems. Hope that this helps, Roger url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jun 24, 2009, at 4:07 PM, Jonathan Greenberg wrote: Rers: I installed R 2.9.0 from the Debian package manager on our amd64 system that currently has 6GB of RAM -- my first question is whether this installation is a true 64-bit installation (should R have access to 4GB of RAM?) I suspect so, because I was running an rqss() (package quantreg, installed via install.packages() -- I noticed it required a compilation of the source) and watched the memory usage spike to 4.9GB (my input data contains 500,000 samples). With this said, after 30 mins or so of processing, I got the following error: tahoe_rq - rqss(ltbmu_4_stemsha_30m_exp.img~qss(ltbmu_eto_annual_mm.img),tau=. 99,data=boundary_data) Error: cannot allocate vector of size 1.5 Gb The dataset is a bit big (300mb or so), so I'm not providing it unless necessary to solve this memory problem. Thoughts? Do I need to compile either the main R by hand or the quantreg package? --j __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Memory issues on a 64-bit debian system (quantreg)
my earlier comment is probably irrelevant since you are fitting only one qss component and have no other covariates. A word of warning though when you go back to this on your new machine -- you are almost surely going to want to specify a large lambda for the qss component in the rqss call. The default of 1 is likely to produce something very very rough with such a large dataset. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jun 24, 2009, at 5:04 PM, Jonathan Greenberg wrote: Yep, its looking like a memory issue -- we have 6GB RAM and 1GB swap -- I did notice that the analysis takes far less memory (and runs) if I: tahoe_rq - rqss(ltbmu_4_stemsha_30m_exp.img~ltbmu_eto_annual_mm.img,tau=. 99,data=boundary_data) (which I assume fits a line to the quantiles) vs. tahoe_rq - rqss(ltbmu_4_stemsha_30m_exp.img~qss(ltbmu_eto_annual_mm.img),tau=. 99,data=boundary_data) (which is fitting a spline) Unless anyone else has any hints as to whether or not I'm making a mistake in my call (beyond randomly subsetting the data -- I'd like to run the analysis on the full dataset to begin with) -- I'd like to fit a spline to the upper 1% of the data, I'll just wait until my new computer comes in next week which has more RAM. Thanks! --j roger koenker wrote: Jonathan, Take a look at the output of sessionInfo(), it should say x86-64 if you have a 64bit installation, or at least I think this is the case. Regarding rqss(), my experience is that (usually) memory problems are due to the fact that early on the processing there is a call to model.matrix() which is supposed to create a design, aka X, matrix for the problem. This matrix is then coerced to matrix.csr sparse format, but the dense form is often too big for the machine to cope with. Ideally, someone would write an R version of model.matrix that would permit building the matrix in sparse form from the get-go, but this is a non-trivial task. (Or at least so it appeared to me when I looked into it a few years ago.) An option is to roll your own X matrix: take a smalller version of the data, apply the formula, look at the structure of X and then try to make a sparse version of the full X matrix. This is usually not that difficult, but usually is based on a rather small sample that may not be representative of your problems. Hope that this helps, Roger url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Urbana, IL 61801 On Jun 24, 2009, at 4:07 PM, Jonathan Greenberg wrote: Rers: I installed R 2.9.0 from the Debian package manager on our amd64 system that currently has 6GB of RAM -- my first question is whether this installation is a true 64-bit installation (should R have access to 4GB of RAM?) I suspect so, because I was running an rqss() (package quantreg, installed via install.packages() -- I noticed it required a compilation of the source) and watched the memory usage spike to 4.9GB (my input data contains 500,000 samples). With this said, after 30 mins or so of processing, I got the following error: tahoe_rq - rqss(ltbmu_4_stemsha_30m_exp.img~qss(ltbmu_eto_annual_mm.img),tau=. 99,data=boundary_data) Error: cannot allocate vector of size 1.5 Gb The dataset is a bit big (300mb or so), so I'm not providing it unless necessary to solve this memory problem. Thoughts? Do I need to compile either the main R by hand or the quantreg package? --j __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Citing R/Packages Question
I've had an email exchange with the authors of a recent paper in Nature who also made a good faith effort to cite both R and the quantreg package, and were told that the Nature house style didn't allow such citations so they were dropped from the published paper and the supplementary material appearing on the Nature website. Since the CRAN website makes a special effort to make prior versions of packages available, it would seem to me to be much more useful to cite version numbers than access dates. There are serious questions about the ephemerality of url citations, not all of which are adequately resolved by the Wayback machine, and access dating, but it would be nice to have some better standards for such contingent citations rather than leave authors at the mercy of copy editors. I would also be interested in suggestions by other contributors. url:www.econ.uiuc.edu/~rogerRoger Koenker email rkoen...@uiuc.edu Department of Economics vox:217-333-4558University of Illinois fax:217-244-6678Champaign, IL 61820 On May 8, 2009, at 5:27 PM, Derek Ogle wrote: I used R and the quantreg package in a manuscript that is currently in the proofs stage. I cited both R and quantreg as suggested by citation() and noted the version of R and quantreg that I used in the main text as All tests were computed with the R v2.9.0 statistical programming language (R Development Core 2008). Quantile regressions were conducted with the quantreg v4.27 package (Koenker 2008) for R. The editor has asked me to also provide the date when the webpage was accessed for both R and quantreg. This does not seem like an appropriate request to me as both R and the quantreg package are versioned. This request seems to me to be the same as asking someone when they purchased commercial package X version Y (which I don't think would be asked). Am I thinking about this correctly or has the editor made a valid request? I would be interested in any comments or opinions. Dr. Derek H. Ogle Associate Professor of Mathematical Sciences and Natural Resources Northland College 1411 Ellis Avenue Box 112 Ashland, WI 715.682.1300 www.ncfaculty.net/dogle/ [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Question of Quantile Regression for Longitudinal Data
I was trying to resist responding to this question since the original questioner had already been admonished twice last october about asking questions on R-help about posted code that was not only not a part of R-base, but not even a part of an R package. But the quoted comment about Stata is too enticing a provocation to resist. First, it should be said that omitting intercepts in any regression setting should be undertaken at one's peril it is generally a very dangerous activity, somewhat akin to fitting interactions without main effects, but if there is a good rational for it, it is no different in principle for median regression than for mean regression. It may well be that Stata prohibits this sort of thing out of some sort of paternalistic motive, but in R the usual formula convention y ~ x1 + x2 -1 suffices. Of course it situations in which such a formula is used for several quantiles it should be understood that it is forcing each conditional quantile function through the origin effectively implies that the conditional distribution degenerates to a point mass at the origin. Second, I would like to remark that closed-form solutions are in the eye of the beholder, and many people who can recall the infamous formula: betahat = (X'X)^{-1} X'y would be hard pressed to dredge up enough linear algebra to use the formula for anything more than the bivariate case on the proverbial desert island without the aid of their trusty laptop Friday. Finally, cbind(1,x) does introduce an intercept in the code originally asked about, so if you don't want an intercept don't do that, but be sure that that is really want you want to do. url:www.econ.uiuc.edu/~rogerRoger Koenker email rkoen...@uiuc.edu Department of Economics vox:217-333-4558University of Illinois fax:217-244-6678Champaign, IL 61820 On Apr 26, 2009, at 6:35 AM, Tirthankar Chakravarty wrote: This is a nontrivial problem. This comes up often on the Statalist (-qreg- is for cross-section quantile regression): You want to fit a plane through the origin using the L-1 norm. This is not as easy as with L-2 norm (LS), as it is more than a matter of dropping a constant predictor yet otherwise using the same criterion of fit. You are placing another constraint on a problem that already does not have a closed-form solution, and it does not surprise me that -qreg- does not support this. (N.J. Cox) http://www.stata.com/statalist/archive/2007-10/msg00809.html You will probably have to program this by hand. Note also the degeneracy conditions in Koenker (2003, pg. 36--). I am not sure how this extends to panel data though. References: @book{koenker2005qre, title={{Quantile Regression; Econometric Society Monographs}}, author={Koenker, R.}, year={2005}, publisher={Cambridge University Press} } T On Sun, Apr 26, 2009 at 8:24 AM, Helen Chen 96258...@nccu.edu.tw wrote: Hi, I am trying to estimate a quantile regression using panel data. I am trying to use the model that is described in Dr. Koenker's article. So I use the code the that is posted in the following link: http://www.econ.uiuc.edu/~roger/research/panel/rq.fit.panel.R How to estimate the panel data quantile regression if the regression contains no constant term? I tried to change the code of rq.fit.panel by delect X=cbind(1,x) and would like to know is that correct ? Thanks I really would appreciate some suggestions. Best Helen Chen -- View this message in context: http://www.nabble.com/Question-of-%22Quantile-Regression-for-Longitudinal-Data%22-tp23239896p23239896.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. -- To every ω-consistent recursive class κ of formulae there correspond recursive class signs r, such that neither v Gen r nor Neg(v Gen r) belongs to Flg(κ) (where v is the free variable of r). __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] any other fast method for median calculation
There is a slightly faster algorithm in my quantreg package, see kuantile() but this is only significant when sample sizes are very large. In your case you really need a wrapper that keeps the loop over columns within some lower level language. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Apr 13, 2009, at 11:29 PM, Zheng, Xin (NIH) [C] wrote: Hi there, I got a data frame with more than 200k columns. How could I get median of each column fast? mapply is the fastest function I know for that, it's not yet satisfied though. It seems function median in R calculates median by sort and mean. I am wondering if there is another function with better algorithm. Any hint? Thanks, Xin Zheng __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
[R] extract argument names
I have a vector of character strings that look like R expressions: a - paste(qss(,paste(x,1:6,sep = ) ,, lambda =100), sep = ) a [1] qss(x1, lambda =100) qss(x2, lambda =100) qss(x3, lambda =100) [4] qss(x4, lambda =100) qss(x5, lambda =100) qss(x6, lambda =100) That I would like to operate on to obtain the names of the first argument, i.e. foo(a) [1] x1 x2 x3 x4 x5 x6 I thought there was some simple idiom involving deparse, but it is eluding my searches. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 __ R-help@r-project.org mailing list 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.
Re: [R] Help using smooth.spline with zoo object
For equally spaced observations this is quite simple and available in various packages but I like Gabor Grothendieck's version (which didn't come up immediately in my Rseek search: hpfilter - function(y,lambda=1600) eye - diag(length(y)) solve(eye+lambda*crossprod(diff(eye,d=2)),y)} url:www.econ.uiuc.edu/~rogerRoger Koenker email rkoen...@uiuc.edu Department of Economics vox:217-333-4558University of Illinois fax:217-244-6678Champaign, IL 61820 On Apr 4, 2009, at 2:19 PM, Rob Denniker wrote: Can someone please show me how to smooth time series data that I have in the form of a zoo object? I have a monthly economies series and all I really need is to see a less jagged line when I plot it. If I do something like s - smooth.spline(d.zoo$Y, spar = 0.2) plot(predict(s,index(d.zoo)), xlab = Year) # not defined for Date objects and if I do something like plot(predict(s,as.numeric(index(d.zoo))), xlab = Year) # one straight line, no matter the value of spar What am I doing wrong? (The unsmoothed series plots just fine - a noisy upward trend) Thanks in advance. P.S. I would really just like to keep varying the penalty parameter 'lambda' of a Hodrick-Prescott filter until I get a 'smooth enough' series, but an archive search suggests that if I ask for this, some smarty pants will reply saying that the HP filter is just the smooth.spline function applied with a particular choice of lambda and equally spaced observations. Receive Notifications of Incoming Messages Easily monitor multiple email accounts access them with a click. Visit http://www.inbox.com/notifier and check it out! __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] orthogonal/perpendicular distance from regression line
For the bivariate case: g - function(b,x,y) (abs(y - b[1] - b[2] * x))/sqrt(1 + crossprod(b)) url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Feb 20, 2009, at 11:12 AM, GAF wrote: Hi there, I am trying to measure orthogonal/perpendicular distances from regression lines (i.e. the shortest distance from a given point to my regression line). As it sounds rather easy (basically just an orthogonal/perpendicular residual) I hoped that there was some function in R that can do that. All efforts so far remained unsuccessful, however. Does anybody know? Thnx and cheers, Philipp -- View this message in context: http://www.nabble.com/orthogonal-perpendicular-distance-from-regression-line-tp22123179p22123179.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] alpha shape function
As it happens, I have also been looking into this. I began by considering Ken Clarkson's hull: http://www.netlib.org/voronoi/hull.html but eventually discovered that its alpha shapes don't seem to treat holes in regions, only simply connected regions. (I would be happy to hear to the contrary, if someone has other experience.) There is a nice matlab implementation http://www.mathworks.com/matlabcentral/fileexchange/6760 which I've also experimented with using R.matlab, but this is still very early days. The matlab version is only 2d, whereas Clarkson does moderate d which extends at least to 3d. Given tripack, it seems like alpha-shapes shouldn't be such a big enterprise, and might make a nice project for someone with an interest in computational geometry. Hint, hint. Nudge, Nudge. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Feb 18, 2009, at 7:20 AM, Pedroso MOACIR wrote: Hi all, I want to approximate te shape of an area defined by a set of points. The convex hull is not good enough, but I think that an alpha shape would be fine. I did an RSiteSearch(), google search, RSeek.org search, looked at the CRAN Views, but was unable do find a function in R that computes the alpha shape. Does anyone know if there is such a function in R? Theank you very much. Moacir Pedroso Embrapa - Brazil __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Percentiles/Quantiles with Weighting
url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Feb 17, 2009, at 1:58 PM, Brigid Mooney wrote: Thanks for pointing me to the quantreg package as a resource. I was hoping to ask be able to address one quick follow-up question... I get slightly different variants between using the rq funciton with formula = mydata ~ 1 as I would if I ran the same data using the quantile function. Example: mydata - (1:10)^2/2 pctile - seq(.59, .99, .1) quantile(mydata, pctile) 59%69%79%89%99% 20.015 26.075 32.935 40.595 49.145 rq(mydata~1, tau=pctile) Call: rq(formula = mydata ~ 1, tau = pctile) Coefficients: tau= 0.59 tau= 0.69 tau= 0.79 tau= 0.89 tau= 0.99 (Intercept)18 24.532 40.550 Degrees of freedom: 10 total; 9 residual Is it correct to assume this is due to the different accepted methods of calculating quantiles? If so, do you know where I would be able to see the algorithms used in these functions? I'm not finding it in the documentation for function rq, and am new enough to R that I don't know where those references would generally be. Yes, quantile() in base R documents 9 varieties of quantiles, 2 more than William Empson's famous 7 Types of Ambiguity. In quantreg the function rq() finds a solution to an underlying optimization problem and doesn't ask any further into the nature of the ambiguity -- it does often produce a warning indicating that there may be more than one solution. The default base R quantile is interpolated, while the default rq() with method = br using the simplex algorithm finds an order statistic, typically. If you prefer something more like interpolation, you can try rq() with method = fn which is using an interior point algorithm and when there are multiple solutions it tends to produce something more like the centroid of the solution set. I hope that this helps. On Tue, Feb 17, 2009 at 12:29 PM, roger koenker rkoen...@uiuc.edu wrote: http://www.nabble.com/weighted-quantiles-to19864562.html#a19865869 gives one possibility... url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Feb 17, 2009, at 10:57 AM, Brigid Mooney wrote: Hi All, I am looking at applications of percentiles to time sequenced data. I had just been using the quantile function to get percentiles over various periods, but am more interested in if there is an accepted (and/or R-implemented) method to apply weighting to the data so as to weigh recent data more heavily. I wrote the following function, but it seems quite inefficient, and not really very flexible in its applications - so if anyone has any suggestions on how to look at quantiles/percentiles within R while also using a weighting schema, I would be very interested. Note - this function supposes the data in X is time-sequenced, with the most recent (and thus heaviest weighted) data at the end of the vector WtPercentile - function(X=rnorm(100), pctile=seq(.1,1,.1)) { Xprime - NA for(i in 1:length(X)) { Xprime - c(Xprime, rep(X[i], times=i)) } print(Percentiles:) print(quantile(X, pctile)) print(Weighted:) print(Xprime) print(Weighted Percentiles:) print(quantile(Xprime, pctile, na.rm=TRUE)) } WtPercentile(1:10) WtPercentile(rnorm(10)) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Help with rgl
Why I love R [Number 6]: Chinese extend a helping hand to Russians who happen to be in Brazil about a package written in Germany. Trotsky would be proud -- and amazed! url:www.econ.uiuc.edu/~rogerRoger Koenker email rkoen...@uiuc.edu Department of Economics vox:217-333-4558University of Illinois fax:217-244-6678Champaign, IL 61820 On Feb 17, 2009, at 8:37 PM, Yihui Xie wrote: (1) you'll need ImageMagick installed to use the command convert to convert image sequences into GIF animations; see ?movie3d (2) viewport is read only!! see ?open3d carefully Regards, Yihui -- Yihui Xie xieyi...@gmail.com Phone: +86-(0)10-82509086 Fax: +86-(0)10-82509086 Mobile: +86-15810805877 Homepage: http://www.yihui.name School of Statistics, Room 1037, Mingde Main Building, Renmin University of China, Beijing, 100872, China On Tue, Feb 17, 2009 at 2:45 AM, Iuri Gavronski i...@ufrgs.br wrote: Hi, I don't know much about the RGL package, and I have read the documentation and tried some parameters, with no luck... I would like to generate a movie from a 3D object (code below), where the vortex A is closer to the observer, and then the object rotates and the B vortex gets closer. I would like to capture this movie to a file. By the way, I am not being able to insert unicode text with text3d. rgl 0.82, R 2.8.1, Windows Vista. Any help would be appreciated. Code follows: library(rgl) open3d() coord.1=c(0,100,0) coord.2=c(100,100,0) coord.3=c(100,0,0) coord.4=c(0,0,0) coord.5=c(50,50,70) pyrcolor=red triangles3d(rbind(coord.1,coord.4,coord.5),color=pyrcolor) triangles3d(rbind(coord.1,coord.2,coord.5),color=pyrcolor) triangles3d(rbind(coord.2,coord.3,coord.5),color=pyrcolor) triangles3d(rbind(coord.3,coord.4,coord.5),color=pyrcolor) quads3d(rbind(coord.1,coord.2,coord.3,coord.4),color=pyrcolor) vertices = LETTERS[1:5] text3d(coord.1,text=vertices[1],adj=1,color=blue) text3d(coord.2,text=vertices[2],adj=0,color=blue) text3d(coord.3,text=vertices[3],adj=0,color=blue) text3d(coord.4,text=vertices[4],adj=1,color=blue) text3d(coord.5,text=vertices[5],adj=0,color=blue) # couldn't make this work... #open3d(viewport=c(0,0,686,489)) #par3d(zoom = 1.157625) filename = piramide.png rgl.snapshot(filename) __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Concave Hull
Actually, I think that the survey on alpha shapes available from: http://www.cs.duke.edu/~edels/Surveys/ would be more closely aligned with what Michael was interested in... url:www.econ.uiuc.edu/~rogerRoger Koenker email rkoen...@uiuc.edu Department of Economics vox:217-333-4558University of Illinois fax:217-244-6678Champaign, IL 61820 On Jan 20, 2009, at 6:06 PM, David Winsemius wrote: The OP was asking whether concave hulls have been implemented. He wasn't very helpful with his link giving the example, since it was to the outside of a frame-based website. Perhaps this link (see the bottom of that page) will be more helpful: http://get.dsi.uminho.pt/local/results.html It has been discussed (briefly) in r-help: http://finzi.psych.upenn.edu/R/Rhelp02a/archive/75574.html Some of the material in Loader's Local Regression and Likelihood regarding classification looks potentially applicable. When first I saw this question I expected that one of the r-sig-Geo folks would have a ready answer. Perhaps a follow-up there would be a reasonable next step? -- David Winsemius On Jan 20, 2009, at 6:37 PM, Charles Geyer wrote: Message: 64 Date: Mon, 19 Jan 2009 15:14:34 -0700 From: Greg Snow greg.s...@imail.org Subject: Re: [R] Concave Hull To: Michael Kubovy kub...@virginia.edu, r-help r-help@r-project.org Message-ID: b37c0a15b8fb3c468b5bc7ebc7da14cc61c8360...@lp-exmbvs10.co.ihc.com Content-Type: text/plain; charset=us-ascii I don't know if it is the same algorithm or not, but there is the function chull that finds the convex hull. Also the R function redundant in the contributed package rcdd efficiently finds convex hulls in d-dimensional space for arbitrary d (chull only does d = 2). See Sections 4.2 and 5.2 of the rcdd package vignette. Hope this helps, -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- project.org] On Behalf Of Michael Kubovy Sent: Saturday, January 17, 2009 9:49 AM To: r-help Subject: [R] Concave Hull Dear Friends, Here is an algorithm for finding concave hulls: http://get.dsi.uminho.pt/local/ Has anyone implemented such an algorithm in R? RSiteSearch('concave hull') didn't reveal one (I think). _ Professor Michael Kubovy University of Virginia Department of Psychology Postal Address: P.O.Box 400400, Charlottesville, VA 22904-4400 Express Parcels Address: Gilmer Hall, Room 102, McCormick Road, Charlottesville, VA 22903 Office:B011;Phone: +1-434-982-4729 Lab:B019; Phone: +1-434-982-4751 WWW:http://www.people.virginia.edu/~mk9y/ Skype name: polyurinsane [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] lm: how are polynomial functions interpreted?
[E] Y = a + b*sin(d*x+phi) isn't a linear model and therefore can't be estimated with lm() -- you will need some heavier artillery. Linear as in lm() means linear in parameters. (As it happens, I'm adapting Gordon Smyth's pronyfreq S code for the above problem this afternoon, and have been wondering why someone else hasn't already done this? Any clues? url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Jan 12, 2009, at 4:57 PM, Carl Witthoft wrote: Well. *_* , I think it should have been clear that this was not a question for which any code exists. In fact, I gave two very specific examples of function calls. The entire point of my question was not what's up with my (putative) code and data but rather to try to understand the overarching philosophy of the way lm() treats the function it's given. I do understand the sneaky ways to make it do a linear fit with or without forcing the origin. And, sure, I could have run a data set thru a bunch of different quadratic-like functions to try to see what happens. Let me pick a more complicated example. The general case of a sin fit might be Y = a + b*sin(d*x+phi) .(where, to be pedantic, x is the only data input. All others are coefficients to be found) If I try y-lm(yin~I(sin(x))), what is the actual fit function? And so on. That's why I was hoping for a more general explanation of what lm() does. Charles C. Berry wrote: On Mon, 12 Jan 2009, c...@witthoft.com wrote: [nothing deleted] matplot(1:100, lm(rnorm(100)~poly(1:100,4),x=T)$x ) # for example __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html Ahem! and provide commented, minimal, self-contained, reproducible code. ..^ Charles C. Berry(858) 534-2098 Dept of Family/ Preventive Medicine E mailto:cbe...@tajo.ucsd.eduUC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] The R Inferno
I think that this continuation constitutes what Pat calls hijacking the thread at the end of his new and magnificent opus. The original thread should be reserved for kudos to Pat. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Jan 9, 2009, at 11:52 AM, Andrew Choens wrote: now if only i could get tips to sort a 5 column * 1 million rows dataset in less than ..eternity May I suggest mySQL, postgreSQL, etc.? If what you need to do is a basic sort, a database is going to be faster than R. -- Insert something humorous here. :-) __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Drawing from an empirical distribution
Sure, but it would be more 'fun' to modify ecdf() slightly to produce an ecqf() function -- essentially reversing the arguments to approxfun()-- and then use ecqf(runif(whatever)) no nit-picking about efficiency, please. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Jan 6, 2009, at 4:42 PM, Antonio, Fabio Di Narzo wrote: If the ecdf is 'ecdf(x)', do just: sample(x, size=whatever, replace=TRUE) HTH, Antonio. 2009/1/6 culpritNr1 ig2ar-s...@yahoo.co.uk: Hi All, Does anybody know if there is a simple way to draw numbers from an empirical distribution? I know that I can plot the empirical cumulative distribution function this easy: plot(ecdf(x)) Now I want to pick a number between 0 and 1 and go back to domain of x. Sounds simple to me. Any suggestion? Thank you, Your culprit (everybody needs a culprit) -- View this message in context: http://www.nabble.com/Drawing-from-an-empirical-distribution-tp21320810p21320810.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. -- Antonio, Fabio Di Narzo Ph.D. student at Department of Statistical Sciences University of Bologna, Italy __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Drawing from an empirical distribution
Nit-picking about syntax does seem needed, mea culpa, I intended something more like: Qn - ecqf(x) Qn(runif(whatever)) On Jan 6, 2009, at 5:06 PM, roger koenker wrote: Sure, but it would be more 'fun' to modify ecdf() slightly to produce an ecqf() function -- essentially reversing the arguments to approxfun()-- and then use ecqf(runif(whatever)) no nit-picking about efficiency, please. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Jan 6, 2009, at 4:42 PM, Antonio, Fabio Di Narzo wrote: If the ecdf is 'ecdf(x)', do just: sample(x, size=whatever, replace=TRUE) HTH, Antonio. 2009/1/6 culpritNr1 ig2ar-s...@yahoo.co.uk: Hi All, Does anybody know if there is a simple way to draw numbers from an empirical distribution? I know that I can plot the empirical cumulative distribution function this easy: plot(ecdf(x)) Now I want to pick a number between 0 and 1 and go back to domain of x. Sounds simple to me. Any suggestion? Thank you, Your culprit (everybody needs a culprit) -- View this message in context: http://www.nabble.com/Drawing-from-an-empirical-distribution-tp21320810p21320810.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. -- Antonio, Fabio Di Narzo Ph.D. student at Department of Statistical Sciences University of Bologna, Italy __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] plotting density for truncated distribution
Default kernel density estimation is poorly suited for this sort of situation. A better alternative is logspline -- see the eponymous package -- you can specify lower limits for the distribution as an option. url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Nov 25, 2008, at 10:43 AM, Jeroen Ooms wrote: I am using density() to plot a density curves. However, one of my variables is truncated at zero, but has most of its density around zero. I would like to know how to plot this with the density function. The problem is that if I do this the regular way density(), values near zero automatically get a very low value because there are no observed values below zero. Furthermore there is some density below zero, although there are no observed values below zero. This illustrated the problem: mydata - rnorm(10); mydata - mydata[mydata0]; plot(density(mydata)); the 'real' density is exactly the right half of a normal distribution, so truncated at zero. However using the default options, the line seems to decrease with a nice curve at the left, with some density below zero. This is pretty confusing for the reader. I have tried to decrease the bw, masks (but does not fix) some of the problem, but than also the rest of the curve loses smoothness. I would like to make a plot of this data that looks like the right half of a normal distribution, while keeping the curve relatively smooth. Is there any way to specify this truncation in the density function, so that it will only use the positive domain to calculate density? -- View this message in context: http://www.nabble.com/plotting-density-for-truncated-distribution-tp20684995p20684995.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Need help: Loading a large data set.
Both SparseM and Matrix have facilities for rbind and cbind that allow you to concatenate pieces of sparse matrices together. On Nov 23, 2008, at 2:19 PM, Atul Kulkarni wrote: Hi All, I am dealing with a large data set which translates in to a sparse matrix, I want to load this data that is spread in approximately 17000+ files each defining a row and each file has variable number of records that come with its column number and the value that they store. I wish to load this data in memory from these files one by one. Is there anyway I can do this in R, before I start processing? I am sure this is not the first time R or the community is confronted with this kind of a problem but I could not find the documentation for loading data in to sparse matrix I found quite a few packages for sparse matrix but they all were concentrating on how to do various operations with the matrix once the matrix is loaded. I need to first load the data in the system before I can think about analysing. Regards, Atul. Graduate Student, Department of Computer Science, University of Minnesota Duluth, Duluth, MN, 55812. www.d.umn.edu/~kulka053 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Dequantizing
I'm rather doubtful that you can improve on the uniform jittering strategy you originally considered. It would require intimate knowledge about the non-uniformity of the density in the spacings between your quantized version. But if you really _knew_ the parent distribution then something like the following might have been what you had in mind: # Toy dequantization example rate - 1 x - sort(rexp(100,rate)) xu - x + runif(100) y - floor(x) ty - table(y) p - c(0,cumsum(table(y)/length(y))) pup - p[-1] plo - p[-length(p)] fun - function(ty,plo,pup) qexp(runif(ty,plo,pup),rate) z - unlist(mapply(fun, ty = ty, plo = plo, pup = pup)) url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Nov 20, 2008, at 9:43 AM, Stavros Macrakis wrote: I have some data measured with a coarsely-quantized clock. Let's say the real data are q- sort(rexp(100,.5)) The quantized form is floor(q), so a simple quantile plot of one against the other can be calculated using: plot(q,type=l); points(floor(q),col=red) which of course shows the characteristic stair-step. I would like to smooth the quantized form back into an approximation of the underlying data. The simplest approach I can think of adds a uniform random variable of the size of the quantization: plot(q,type=l); points(floor(q),col=red); points(floor(q)+runif(100,0,1),col=blue) This gives pretty good results for uniform distributions, but less good for others (like exponential). Is there a better interpolation/smoothing function for cases like this, either Monte Carlo as above or deterministic? Thanks, -s __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Oja median
Cross posting is sociopathic. The Oja median _is_ robust in several reasonable senses, but admittedly not from a breakdown perspective. There are several other options: for a good review see, e.g. Multivariate Median, A. Niinimaa and H. Oja Encyclopedia of Statistical Sciences url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Nov 19, 2008, at 5:40 AM, [EMAIL PROTECTED] wrote: Hi Roger, As we know that The Oja median has (finite) breakdown point 2/n, i.e., is not robust in any reasonable sense, and is quite expensive to compute, so do we have some better methodology to compute multivariate median Rahul Agarwal Analyst Equities Quantitative Research UBS_ISC, Hyderabad On Net: 19 533 6363 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Quantile Regression fixed effects model
see: http://www.econ.uiuc.edu/~roger/research/panel/rq.fit.panel.R url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Nov 14, 2008, at 1:20 PM, dimitris kapetanakis wrote: Dear R users, I am trying to estimate a fixed effect quantile regression for different quantiles. As Dr. Koenker mention on his article (2004) the model should be estimated simultaneously so it is going to have the same fixed effects for all quantiles. The problem is that when I am using the following code R estimates the model for each quantile separately so for every quantile there is different fixed effects estimation. Can you help me to estimate the quantile fixed effects simultaneously for different quantiles? for(i in 1:7){ assign(qr.y.[i], rq(y~factor(year)+factor(state)+x+I(x^2)+I(x^3), tau=quant[i], data=file)) } Thank you in advance Dimitris -- View this message in context: http://www.nabble.com/Quantile-Regression-fixed-effects-model-tp20506811p20506811.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
[R] Election Maps
Those of you with an interest in the US election and/or statistical graphics may find the maps at: http://www-personal.umich.edu/~mejn/election/2008/ interesting. url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 __ R-help@r-project.org mailing list 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.
Re: [R] Quantile Regression for Longitudinal Data:error message
If you are going to insist on doing such things you will have to learn to read the documentation. In this case if you do a traceback() you will see that the error is occurring in rq.fit.slm and when you do ?rq.fit.slm you will see that there are several storage sizes that can be adjusted: nsubmax: upper bound of the dimension of lindx in Cholesky factorization; computed automatically inside the routine if missing. tmpmax: upper bound of the working array in Cholesky factorization; computed automatically inside the routine if missing. nnzlmax: upper bound of the non-zero entries in the Cholesky factor L; computed automatically inside the routine if missing. cachsz: size of the cache on the machine; default to 64. Since you don't give, as stipulated by the posting guide, a reproducible example, one can only speculate about what you are trying to do, but it appears that you are trying to simultaneously estimate 99 quantiles with a sample size of length(w) = 99 observations and this is Unsound. I would also like to repeat my earlier comment that this is NOT a proper R help question since it refers not to something that is in base R, nor even something in a package, but to code that I happen to have posted to complement a paper that was published several years ago. url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Oct 31, 2008, at 11:32 AM, Helen Chen wrote: Quantile Regression for Longitudinal Data. Hi, I am trying to estimate a quantile regression using panel data. I am trying to use the model that is described in Dr. Koenker's article. So I use the code the that is posted in the following link: http://www.econ.uiuc.edu/~roger/research/panel/rq.fit.panel.R I am trying to change the number quantiles being estimated. I change the codes about w and taus ,as following w=c(.01,.02,.03,.04,.05,.06,.07,.08,.09,.1,.11,.12,.13,.14,.15,.16,.17 ,.18,.19,.2,.21,.22,.23,.24,.25,.26,.27,.28,.29,.3,.31,.32,.33,.34,.35 ,.36,.37,.38,.39,.4,.41,.42,.43,.44,.45,.46,.47,.48,.49,.5,.49,.48,.47 ,.46,.45,.44,.43,.42,.41,.4,.39,.38,.37,.36,.35,.34,.33,.32,.31,.3,.29 ,. 28,.27,.26,.25,.24,.23,.22,.21,.2,.19,.18,.17,.16,.15,.14,.13,.12,.11 ,.1,.09,.08,.07,.06,.05,.04,.03,.02,.01) ,taus=(1:99)/100,lambda = 1 But I get error message: .local(x, pivot, ...) : Increase tmpmax So I am wondering if I am doing something wrong or I mistake w's meaning. Thanks I really would appreciate some suggestions. Best Helen Chen -- View this message in context: http://www.nabble.com/Quantile-Regression-for-Longitudinal-Data%3Aerror-message-tp20269386p20269386.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Question of Quantile Regression for Longitudinal Data.
This is really NOT an R-help question. You just need to read the code more carefully and everything will be revealed there. This isn't like the Pisan Cantos. url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Oct 26, 2008, at 8:59 PM, Helen Chen wrote: Hi, I am trying to estimate a quantile regression using panel data. I am trying to use the model that is described in Dr. Koenker's article. So I use the code the that is posted in the following link: http://www.econ.uiuc.edu/~roger/research/panel/rq.fit.panel.R This code run perfectly. Then I want to know what the result means.The result show $ierr,$it,and $time. What these estimators means in quantile panel data? And does $coef ,the coefficient estimators order by m*p + n elements? Here comes my questions What do these estimators' mean in quantile panel data? And does the coefficient, $coef, estimators order by m*p + n elements? Thanks I really would appreciate some suggestions. Best Helen Chen -- View this message in context: http://www.nabble.com/Question-of-%22Quantile-Regression-for-Longitudinal-Data%22.-tp20180665p20180665.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Error in Switch in KhmaladzeTest
Hey back, If you had RTFM you would know that the first argument of the function KhmaladzeTest was supposed to be a formula, so try: T - KhmaladzeTest(Logloss~AS+AM+CB+CF+RB+RBR, data=CPBP, nullH=location) url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Oct 15, 2008, at 10:39 AM, [EMAIL PROTECTED] wrote: Hey, My dataset has 1 dependent variable(Logloss) and 7 independent dummy variables(AS,AM,CB,CF,RB,RBR,TS) , it's attached in this email. The problem is I cant finish Khmaladze test because there's an error Error in switch(mode(x), NULL = structure(NULL, class = formula), : invalid formula which I really dont know how to fix. My R version is 2.7.2. The packages loaded are Quantreq and Sparse M, the process is as follows: CPBP-read.table(CPBP.txt) local({pkg - select.list(sort(.packages(all.available = TRUE))) + if(nchar(pkg)) library(pkg, character.only=TRUE)}) Loading required package: SparseM Package SparseM (0.78) loaded. To cite, see citation(SparseM) Package quantreg (4.22) loaded. To cite, see citation(quantreg) fit-rqProcess(Logloss~AS+AM+CB+CF+RB+RBR,data=CPBP,taus=-1) KhmaladzeTest(fit,nullH=location) Error in switch(mode(x), NULL = structure(NULL, class = formula), : invalid formula I would greatly appreciate if you can help me out, thank you very much (See attached file: CPBP.txt) Best Regards, Lulu Ren - ** This E-mail is confidential. It may also be legally privileged. If you are not the addressee you may not copy, forward, disclose or use any part of it. If you have received this message in error, please delete it and all copies from your system and notify the sender immediately by return E-mail. Internet communications cannot be guaranteed to be timely, secure, error or virus-free. The sender does not accept liability for any errors or omissions. ** SAVE PAPER - THINK BEFORE YOU PRINT! CPBP.txt__ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] plot - central limit theorem
Galton's 19th century mechanical version of this is the quincunx. I have a (very primitive) version of this for R at: http://www.econ.uiuc.edu/~roger/courses/476/routines/quincunx.R url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 Jörg Groß wrote: Hi, Is there a way to simulate a population with R and pull out m samples, each with n values for calculating m means? I need that kind of data to plot a graphic, demonstrating the central limit theorem and I don't know how to begin. So, perhaps someone can give me some tips and hints how to start and which functions to use. thanks for any help, joerg __ R-help@r-project.org mailing list 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.
Re: [R] Maximum number of pasted 'code' lines?
Nothing except: Nothing exceeds like excel. url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Oct 14, 2008, at 1:02 PM, Michael Just wrote: Hello, I write most of my R code in excel and then paste it into R. I am wondering if there is a limit to how much I can paste? I want to paste about 19,000 lines of code should this work? I am doing this because when I did it chunks it took about an hour and half. I thought if I could insert it all and leave it for that long that would be better time management. I am still waiting for 19,000 paste to process. Any thoughts or suggestions? Thanks, Michael [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Quantiles of weighted sample
see http://www.nabble.com/weighted-quantiles-to19864562.html#a19865869 url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Oct 8, 2008, at 12:37 PM, Giovanni Petris wrote: Hello! I am wondering if there is a function in R that returns quantiles of a weighted sample, i.e., a set of values, each coming with a different weight. The function quantile() does that only for the case when the weights are all equal. In other words, I am looking for a quantile function that applies to a discrete distribution specified by values and their probabilities. Thanks in advance. Best, Giovanni Petris -- Giovanni Petris [EMAIL PROTECTED] Associate Professor Department of Mathematical Sciences University of Arkansas - Fayetteville, AR 72701 Ph: (479) 575-6324, 575-8630 (fax) http://definetti.uark.edu/~gpetris/ __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
[R] partial matching and dots?
I'm writing a new predict method and would like to be able to pass an argument called se via the ... mechanism. However, predict has a se.fit argument that wants to interpret my specification of se as a partially matched version of se.fit. Surely there a standard treatment for this ailment, but I can't seem to find it. url:www.econ.uiuc.edu/~rogerRoger Koenker email [EMAIL PROTECTED] Department of Economics vox:217-333-4558University of Illinois fax:217-244-6678Champaign, IL 61820 __ R-help@r-project.org mailing list 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.
Re: [R] Quantile Regression for Longitudinal Data. Warning message: In rq.fit.sfn
This is a little esoteric for R-help. As the posting guide says, you should write the package maintainer with this sort of question. Without the data it is difficult to judge what is happening, a couple of possibilities are: o all is well and warning just conveys an exaggerated sense of skepticism about one or more of the Cholesky steps, o or, if the coefs really look like nonsense, you have a near singularity in the design, if so you can try to look in more detail at the dense form of design and try to see what the problem is, provided of course that it is small enough to look at with the usual tools, svd, etc. url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Sep 30, 2008, at 12:53 PM, dimitris kapetanakis wrote: Hi, I am trying to estimate a quantile regression using panel data. I am trying to use the model that is described in Dr. Koenker's article. So I use the code the that is posted in the following link: http://www.econ.uiuc.edu/~roger/research/panel/rq.fit.panel.R While this code run perfectly, it does not work for my data providing a warning message: In rq.fit.sfn(D, y, rhs = a) : tiny diagonals replaced with Inf when calling blkfct So I am wondering if I am doing something wrong or if it is data's problem (which runs perfectly at least for fixed and random effects as well as for quantile regression adding dummies of states and years just like fixed effects). Any help would be highly appreciate Thanks Dimitris ## The code that I use before using the link is: data$lex2-data$lex^2 data$lex3-data$lex^3 data$constant-rep(1,times=length(data$lex)) X-cbind(data$lex, data$lex2, data$lex3, data$constant) y-c(data$ban) s-c(data$state) -- View this message in context: http://www.nabble.com/Quantile-Regression-for-Longitudinal-Data.-Warning-message%3A-In-rq.fit.sfn-tp19747218p19747218.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Modality Test
the diptest package, perhaps? url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Sep 9, 2008, at 11:23 AM, Amin W. Mugera wrote: Dear Readers: I have two issues in nonparametric statistical analysis that i need help: First, does R have a package that can implement the multimodality test, e.g., the Silverman test, DIP test, MAP test or Runt test. I have seen an earlier thread (sometime in 2003) where someone was trying to write a code for the Silverman test of multimodality. Is there any other tests that can enable me to know how many modes are in a distribution? Second, i would like to test whether two distributions are equal. Does R have a package than can implement the Li (1996) test of the equality of two distributions? Is there any other test i can use rather than the Li test? Thank you in advance for your help. Amin Mugera Graduate Student AgEcon Dept. Kansas State University __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] binary order combinations
Does ?combos in the quantreg package do what you want? url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Sep 5, 2008, at 9:58 AM, Dimitri Liakhovitski wrote: Dear all! I have a vector of names names-(V1, V2, V3,., V15) I could create all possible combinations of these names (of all lengths) using R: combos-lapply(1:15,function(x) {combn(names,x) }) I get a list with all possible combinations of elements of 'names' that looks like this (just the very beginning of it): [[1]] - the first element contains all combinations of 1 name [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [, 13] [,14] [1,] V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 [,15] [1,] V15 [[2]] - the second element contains all possible combinations of 2 names [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [1,] V1 V1 V1 V1 V1 V1 V1 V1 V1 V1 V1 V1 V1 [2,] V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 . . . etc. My question is: Is there any way to re-arrange all sub-elements of the above list (i.e., all possible combinations of names such as V1, V1:V3, V1:V2:V4:V5) in a binary system order. More specifically, according to this system: V1=1 V2=2 V3=4 V4=8 V5=16, etc So, I'd like those combinations to be arranged in a vector in the following order: 1. V1 (because V1=1) 2. V2 (because V2=2) 3. V1:V2 (because V1=1 and V2=2 so that 1+2=3) 4. V3 (because V3=4) 5. V1:V3 (because V1=1 and V3=4 so that 1+4=5) 6. V2:V3 (because V2=2 and V3=4 so that 2+4=6) 7. V1:V2:V3 (because V1=1 and V2=2 and V3=4 so that 1+2+4=7) 8. V4 (because V4=8) etc. Is it at all possible? Or maybe there is a way to create the name combinations in such an order in the first place? Thank you very much! Dimitri Liakhovitski __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.