I'll use a modification of this example to ask my question about an apparent alternative way of presenting a partial regression plot, using the effects package.
First, we set up two, positively correlated predictor variables, x1 (which is of our primary interest) and x2 (a covariate of lesser interest), and their effects on the response variabe, y:
set.seed(1)
x <- runif(100, 30, 100)
set.seed(2)
x1 <- x + rnorm(100, sd=10)
set.seed(3)
c <- 1 # positive correlation between predictors
x2 <- c * x + rnorm(100, sd=10)
p1 <- 1 # positive effect of x1
p2 <- -1 # neg. effect of x2
set.seed(1)
y <- 50 + p1*x1 + p2*x2 + rnorm(100, sd=10) # independent effects of x1 and x2, plus noise)
plot(x1, y, main='y ~ x1')
To visualize the unique effect of x1 while accounting for x2, a partial regression plot is generally presented by plotting the residuals of x1 ~ x2 on the horizontal axis, against the residuals of y ~ x2 on the vertical axis.
# y, given x2
mod1 <- lm(y ~ x2 )
y.resid <- resid(mod1)
# x1, given x2
mod2 <- lm(x1 ~ x2)
x.resid <- resid(mod2)
plot(x.resid, y.resid, main='Partial regression: y|x2 ~ x1|x2')
However, the interpretation of this figure is not straight forward.
I found that an alternative partial regression plot can be presented using the effects package, after running the multiple linear regression model mod0 <- lm(y ~ x1 +x2):
library(effects)
effx1 <- effect("x1", mod0, partial.residuals=T)
plot(effx1, smooth.residuals=F)
Which is really more appealing, because the horizontal axis gives the actual values of x1. Also, the vertical axis appear to show the unique (positive) contribution of each x1 to y, but I don't quite understand how x2 is corrected for. I notice that it is not only a matter of scaling the vertical axis differently, as the pattern of the data points is different.
My questions:
- Is this effect plot really an equivalent presentation of the above partial regression plot?
- How is it computed? In other words, how would one produce it manually?
Many thanks!
