Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Constraining slope in stat_smooth with ggplot (plotting ANCOVA)

Tags:

plot

r

ggplot2

Using ggplot(), I am trying to plot the results of an ANCOVA in which slopes of the two linear components are equal: i.e., lm(y ~ x + A). The default behavior for geom_smooth(method = "lm") is to plot separate slopes and intercepts for each level of each factor. For example, with two levels of A

library(ggplot2)
set.seed(1234)

n <- 20

x1 <- rnorm(n); x2 <- rnorm(n)
y1 <- 2 * x1 + rnorm(n)
y2 <- 3 * x2 + (2 + rnorm(n))
A <- as.factor(rep(c(1, 2), each = n))
df <- data.frame(x = c(x1, x2), y = c(y1, y2), A = A)

p <- ggplot(df, aes(x = x, y = y, color = A))
p + geom_point() + geom_smooth(method = "lm")

Default ggplot()

I can fit the ANCOVA separately with lm() and then use geom_abline() to manually add the lines. This approach has a couple of drawbacks like having the lines extend beyond the range of the data and manually specify the colors.

fm <- lm(y ~ x + A, data = df)
summary(fm)

a1 <- coef(fm)[1]
b <- coef(fm)[2]
a2 <- a1 + coef(fm)[3]

p + geom_point() + 
  geom_abline(intercept = a1, slope = b) + 
  geom_abline(intercept = a2, slope = b)

With geom_abline()

I know ancova() in the HH package automates the plotting, but I don't really care for lattice graphics. So I am looking for a ggplot()-centric solution.

library(HH)
ancova(y ~ x + A, data = df)

Is there a method to accomplish this using ggplot()? For this example, A has two levels, but I have situations with 3, 4, or more levels. The formula argument to geom_smooth() doesn't seem to have the answer (as far as I can tell).

like image 570
kmm Avatar asked Nov 22 '10 23:11

kmm


1 Answers

For completeness, this works:

library(ggplot2)
set.seed(1234)

n <- 20

x1 <- rnorm(n); x2 <- rnorm(n)
y1 <- 2 * x1 + rnorm(n)
y2 <- 3 * x2 + (2 + rnorm(n))
A <- as.factor(rep(c(1, 2), each = n))
df <- data.frame(x = c(x1, x2), y = c(y1, y2), A = A)
fm <- lm(y ~ x + A, data = df)

p <- ggplot(data = cbind(df, pred = predict(fm)),
  aes(x = x, y = y, color = A))
p + geom_point() + geom_line(aes(y = pred))
like image 98
kmm Avatar answered Sep 30 '22 10:09

kmm