Re: [R] particle count probability

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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?

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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?

[Roger Koenker]

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?

[Roger Koenker]

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?

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

> 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

[Roger Koenker]

> 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

[Roger Koenker]

> On Oct 6, 2015, at 7:58 AM, Lorenz, David  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 
>> 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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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.

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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?

[Roger Koenker]

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

[Roger Koenker]

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()

[Roger Koenker]

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

[Roger Koenker]

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'?

[Roger Koenker]

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

[Roger Koenker]

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...

[Roger Koenker]

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]

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

[Roger Koenker]

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}

[Roger Koenker]

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}?

[Roger Koenker]

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}?

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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

[Roger Koenker]

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?

[roger koenker]

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 ??

[roger koenker]

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

[roger koenker]

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?

[roger koenker]

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)

[roger koenker]

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)

[roger koenker]

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,

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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)

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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)

[roger koenker]

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)

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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?

[roger koenker]

[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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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.

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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.

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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?

[roger koenker]

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

[roger koenker]

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?

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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

[roger koenker]

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.