# There is grandeur in this view of life

martins bioblogg

## ”How to draw the line” with ggplot2

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.

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])
```

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])
```

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")
```

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

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

Written by mrtnj

31 maj, 2013 at 00:28

### 6 svar

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! […]