Re: Autoplot

[Jean-Olivier Irisson]

On 2012-Jun-18, at 13:11 , Lali wrote:
>
> I am currently working on some quantile regressions using the lmsqreg package and I was wondering if autoplot will work with this object class (lmsqreg.fit); which of course doesn't right off the bat, unless I define the method. The problem is that I am not an R expert and I have no idea how to define the method. I read this R bloggers article but a lot of it is just way over my head. Does anyone have any suggestions on how to go about this?

I have some still-in-development code to fortify output from rq models (quantile regression with the package quantreg, which seems to be the "classic" one). From then on, you would have to write the corresponding autoplot function (which is really:
1. calling fortify
2. calling ggplot() + geom_path() with the appropriate columns from the output of fortify

Hopefully, this, one day, will be in a package.

I paste the below:

fortify.rq <- function(model, data=model$model, ...) {
#
# Extract data from a requession quantile object into a data.frame
#

fitted <- predict(model, ...)
resid <- residuals(model)

# merge with data
res <- data.frame(data, .fitted=fitted, .resid=resid, row.names=NULL)

# identify quantile
res$.tau <- model$tau
res$.tau <- as.factor(res$.tau)

return(res)
}

fortify.rqs <- function(model, data=model$model, ...) {
#
# Extract data from a requession quantile object into a data.frame
#
suppressPackageStartupMessages(require("reshape2", quietly=TRUE))

# get fitted values for all quantiles
fit <- as.data.frame(predict(model))
# convert to long shape
fit <- melt(fit, measure.vars=1:ncol(fit), variable.name=".tau", value.name=".fitted")

# get residuals
resid <- as.data.frame(residuals(model))
resid <- melt(resid, measure.vars=1:ncol(resid), variable.name=".tau", value.name=".resid")
# merge with fit
fit$.resid <- resid$.resid

# identify quantiles
fit$.tau <- model$tau[as.numeric(fit$.tau)]
fit$.tau <- factor(fit$.tau)

# merge with data
res <- data.frame(data, fit, row.names=NULL)

return(res)
}


Jean-Olivier Irisson
---
Observatoire Océanologique
Station Zoologique, B.P. 28, Chemin du Lazaret
06230 Villefranche-sur-Mer
Tel: +33 04 93 76 38 04
Mob: +33 06 21 05 19 90
http://jo.irisson.com/
Send me large files at: http://jo.irisson.com/upload/

[Jean-Olivier Irisson]

On 2012-Jun-19, at 06:29 , Jean-Olivier Irisson wrote:
>
> From then on, you would have to write the corresponding autoplot function (which is really:
> 1. calling fortify
> 2. calling ggplot() + geom_path() with the appropriate columns from the output of fortify

Here are the types of plot you may want to do with the output of fortify

md <- fortify(m)

# check the distribution of residuals
ggplot(md) + geom_histogram(aes(x=.resid), binwidth=0.2) + facet_wrap(~.tau)
# a qqplot is also appropriate

# check homoscedasticity
ggplot(md) + geom_point(aes(x=.fitted, y=.resid))

# plot the fit
ggplot(m) + geom_point(aes(x=size, y=TG_Chol)) + geom_line(aes(x=size, y=.fitted, colour=.tau))
# here my x variable is size and my y variable is TC_Chol, this has to be adapted to your variables of course

[Lali]

Thanks so much for the code, I will try it out, when I figure how to extract the data from the lmsqreg.fit object (I will I just need to sit down and look at it carefully).

[Lali]

Hi Wrathematics,

I am afraid I couldn't follow your article to its whole extent, and that is entirely because I am still learning.

Here is what I used to generate my plot:

females <- read.csv2("females.csv")

female.weight <- lmsqreg.fit(log(females$heart),
                             log(females$weight),
                             pvec = (centiles <- c(0.03, 0.1, 0.5, 0.9, 0.97)),
                             lam.fixed=1)

###To plot, there are two options:

##1. For evaluation

plot.lmsqreg.fit(female.weight)

##2. for publication

plot(log(females$weight),log(females$heart),
     axes=FALSE,
     frame.plot=TRUE,
     xlab="Weight (kg)",
     ylab="Heart weight (g)",
     main="Women",
     col = "white",
     xlim = c(3.2,5.3))

axis(1,at=log(xat <- c(3,5,7,9,11,13,15,17,19)*10),
     labels=xat)

axis(2,at=log(yat <- (2:7)*100),labels=yat)

lpars <- list(lty=rep(c(2,1,2),
                      c(2,1,2)),
              col=c("#444444","#999999","#000000")
              [c(1,2,3,2,1)])

for(l in 1:nrow(female.weight$qsys$outmat) ){
  lines(female.weight$qsys$targetx,
        female.weight$qsys$outmat[l,],
        lty=lpars$lty[l],col=lpars$col[l])
 
  text(female.weight$qsys$targetx[lx <- length(female.weight$qsys$targetx)],
       female.weight$qsys$outmat[l,lx],
       paste(centiles[l]*100,"%",sep=""),adj=c(-0.5,0.5))
}

#After this, I overplot the pathological data, e.g.

#points(log(f.heart$weight[f.heart$cod.S == "CA"]),
#       log(f.heart$heart.w[f.heart$cod.S == "CA"]),
#       pch = 19, col = "red")


I hope it is ok with the group to attach a file, I couldn't find the posting rules.


On Tuesday, June 19, 2012 12:38:07 AM UTC+2, wrathematics wrote:

Hi Lali,

I wrote the article, so I might be able to help.  Can you give me some example code to generate the type of plot you want to reproduce (assuming plot() already works on lmsqreg.fit objects)?  I'm not familiar with this package.

On Tuesday, June 19, 2012 12:38:07 AM UTC+2, wrathematics wrote:

Hi Lali,

I wrote the article, so I might be able to help.  Can you give me some example code to generate the type of plot you want to reproduce (assuming plot() already works on lmsqreg.fit objects)?  I'm not familiar with this package.

[wrathematics]

Well, google decided to eat my original reply.  So let's try again.

My personal opinion is that autoplot methods should seek to replicate functionality of the original plot() functions.  So you really have two questions, with one (in my opinion) not having much in particular to do with autoplot.

I haven't had much time to work on this (and I burned some of that because I forgot about the aes/aes_string problem), but I've basically replicated the functionality of plot.lmsqreg.fit in ggplot.  I'm not proud of this solution, but I'm not going to have any more time to look over it until at least this weekend, so at least for now, it is what it is.

I had to make some insane choices to get things to work right, but it seems to work ok.  It definitely could stand to be significantly cleaned up.  See attached file.

As for your second plot, hopefully this and Jean-Oliver's posts will show you the way.  It looks like it will be a fairly similar translation going from autoplot() --> custom plot, as you did for  plot() --> custom plot...if that makes sense.

Best of luck.

[wrathematics]

And I realized the moment I posted that I forgot to change back the declarations of the dataframes to be local (set them global while testing).  So actually this doesn't work, unless you're comfortable making those dataframes global (I am not).

Not sure if there's a workaround to that or not.  If there's no resolution before this weekend, I'll try to take another look, but can't promise.

Sorry.  Best of luck.

[Lali]

Thank you so much, I will try it out (tomorrow at work) and tell you about it!

[Lali]

I am sorry it took me so long to get back to you, it works perfectly, thanks; this helped me move forward and now I am experimenting with my overplotting.

Thanks very much again for your time,

L

PS How do we mark this as solved?

On Tuesday, June 19, 2012 9:09:04 PM UTC+2, wrathematics wrote: