There is grandeur in this view of life

martins bioblogg

”How to draw the line” with ggplot2

with 6 comments

In a recent tutorial in the eLife journal, Huang, Rattner, Liu & Nathans suggested that researchers who draw scatterplots should start providing not one but three regression lines. I quote,

Plotting both regression lines gives a fuller picture of the data, and comparing their slopes provides a simple graphical assessment of the correlation coefficient. Plotting the orthogonal regression line (red) provides additional information because it makes no assumptions about the dependence or independence of the variables; as such, it appears to more accurately describe the trend in the data compared to either of the ordinary least squares regression lines.

eLife_regressions

Not that new, but I do love a good scatterplot, so I decided to try drawing some lines. I use the temperature and wind variables in the air quality data set (NA values removed). We will  need Hadley Wickham’s ggplot2 and Bendix Carstensen’s, Lyle Gurrin’s and Claus Ekstrom’s MethComp.

library(ggplot2)
library(MethComp)

data(airquality)
data <- na.exclude(airquality)

Let’s first make the regular old scatterplot with a regression (temperature as response; wind as predictor):

plot.y <- qplot(y=Temp, x=Wind, data=data)
model.y <- lm(Temp ~ Wind, data)
coef.y <- coef(model.y)
plot.y + geom_abline(intercept=coef.y[1],
  slope=coef.y[2])

drawing_the_line_y

And then, a regular old scatterplot of the other regression:

plot.x <- qplot(y=Wind, x=Temp, data=data)
model.x <- lm(Wind ~ Temp, data)
coef(model.x)
plot.x + geom_abline(intercept=coef(model.x)[1],
  slope=coef(model.x)[2])

drawing_the_line_x

So far, everything is normal. To put both lines in the same plot, we’ll need to rearrange the coefficients a little. From the above regression we get the equation x = a + by, and we rearrange it to y = – a / b + (1 / b) x.

rearrange.coef <- function(coef) {
  alpha <- coef[1]
  beta <- coef[2]
  new.coef <- c(-alpha/beta, 1/beta)
  names(new.coef) <- c("intercept", "slope")
  return(new.coef)
}
coef.x <- rearrange.coef(coef(model.x))

The third regression line is different: orthogonal, total least squares or Deming regression. There is a function for that in the MethComp package.

deming <- Deming(y=airquality$Temp, x=airquality$Wind)
deming

We can even use the rearrange.coef function above to see that the coefficients of Deming regression does not depend on which variable is taken as the response or predictor:

rearrange.coef(deming)
 intercept      slope 
24.8083259 -0.1906826
Deming(y=airquality$Wind, x=airquality$Temp)[1:2]
 Intercept      Slope 
24.8083259 -0.1906826

So, here is the final plot with all three lines:

plot.y + geom_abline(intercept=coef.y[1],
  slope=coef.y[2], colour="red") +
  geom_abline(intercept=coef.x[1],
  slope=coef.x[2], colour="blue") +
  geom_abline(intercept=deming[1],
  slope=deming[2], colour="purple")

drawing_the_line_3

Now for the dénouement of this post. Of course, there’s already an R function for doing this. It’s even in the same package as I used for the Deming regression; I just didn’t immediately rtfm. It uses base R graphics, though, and I tend to prefer ggplot2:

plot(x=data$Wind, y=data$Temp)
bothlines(x=data$Wind, y=data$Temp, Dem=T, col=c("red", "blue", "purple"))

This is the resulting plot (and one can even check out the source code of bothlines and see that I did my algebra correctly).

drawing_the_line_bothlines

In this case, as well as in a few others that I tried, it seems like one of the ordinary regression lines (in this case the one with wind as response and temperature as predictor) is much closer to the Deming regression line than the other. I wonder under what circumstances that is the case, and if it tells anything useful about the variables. I welcome any thoughts on this matter from you, dear reader.

Literature

Huang L, Rattner A, Liu H, Nathans J. (2013) Tutorial: How to draw the line in biomedical research. eLife e00638 doi:10.7554/eLife.00638

About these ads

Written by mrtnj

31 maj, 2013 at 00:28

6 svar

Subscribe to comments with RSS.

  1. Unfortunately it looks like MethComp has been removed from CRAN.

    Bryan Hanson

    31 maj, 2013 at 15:10

    • Oh, I didn’t notice, probably since I haven’t updated to R 3. I wonder why. :/

      Well, you can probably still install it from archives if you don’t have it. Or do orthogonal regression with the mcr package: http://cran.r-project.org/web/packages/mcr/index.html

      Like so,

      library(mcr)
      model <- mcreg(x=airquality$Temp, y=airquality$Wind, method.reg="Deming")
      getCoefficients(model)

      Cheers,
      martin

      mrtnj

      31 maj, 2013 at 16:18

  2. Found this through r-bloggers.com. Very interesting and useful tutorial, thank you!

    Jon L.

    31 maj, 2013 at 18:33

  3. But why not qr() – quantile regression?

    Drago

    2 juni, 2013 at 00:21

    • Thanks for the comments. Fun to hear from some fellow (I presume) R-blogges readers!

      You tell me, why quantile regression? :)

      (I also assume you mean rq, since qr is the QR-decomposition.)

      Cheers,

      m.

      mrtnj

      2 juni, 2013 at 00:47

  4. […] Neste post do RBloggers explica como! […]


Kommentera

Fyll i dina uppgifter nedan eller klicka på en ikon för att logga in:

WordPress.com Logo

Du kommenterar med ditt WordPress.com-konto. Logga ut / Ändra )

Twitter-bild

Du kommenterar med ditt Twitter-konto. Logga ut / Ändra )

Facebook-foto

Du kommenterar med ditt Facebook-konto. Logga ut / Ändra )

Google+ photo

Du kommenterar med ditt Google+-konto. Logga ut / Ändra )

Ansluter till %s

Följ

Få meddelanden om nya inlägg via e-post.

Gör sällskap med 1 168 andra följare

%d bloggare gillar detta: