Partial residual plots for linear model including an interaction term

Question

My model includes one response variable, five predictors and one interaction term for predictor_1 and predictor_2. I would like to plot partial residual plots for every predictor variable which I would normally realize using the crPlots function from the package car. Unfortunately the function complains that it doesn't work with models that include interaction terms.

Is there another way of doing what I want?

EDIT: I created a small example illustrating the problem

require(car)
R <-  c(0.53,0.60,0.64,0.52,0.75,0.66,0.71,0.49,0.52,0.59)
P1 <- c(3.1,1.8,1.8,1.8,1.8,3.2,3.2,2.8,3.1,3.3)
P2 <- c(2.1,0.8,0.3,0.5,0.4,1.3,0.5,1.2,1.6,2.1)

lm.fit1 <- lm(R ~ P1 + P2)
summary(lm.fit1)
crPlots(lm.fit1) # works fine

lm.fit2 <- lm(R ~ P1*P2)
summary(lm.fit2)
crPlots(lm.fit2) # not available

If you really want help, you should create a [minimal, reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). Include sample data and the code to fit the model, and be as specific as you can about the desired result for the test data so we can test potential solutions. — MrFlick, Jul 25 '14 at 14:18
That is actually a perfect minimal reproducible example. Very helpful. So when you add the interaction term, are you expecting a third plot? It's hard to imagine what the x-axis would be given that them has both P1 and P2 varying in it. Or are you just trying to plot the results for the non-interaction terms from a model that has interaction terms? It's still not 100% clear to me what your desired output is. — MrFlick, Jul 25 '14 at 15:00
Yes, actually I only want the plots for the non-interaction terms in the second model. — ChrKoenig, Jul 25 '14 at 15:19

jlhoward · Accepted Answer · 2014-07-25T21:36:58.510

7

Another way to do this is to put the interaction term in as a separate variable (which avoids hacking the code for crPlot(...)).

df <- data.frame(R,P1,P2,P1.P2=P1*P2)
lm.fit1 <- lm(R ~ ., df)
summary(lm.fit1)
crPlots(lm.fit1)

Note that summary(lm.fit1) yeilds exactly the same result as summary(lm(R~P1*P2,df)).

edited Jul 25 '14 at 21:36

answered Jul 25 '14 at 21:21

jlhoward

58,004
7
97
140

Thanks! I finally came back to this topic and your method worked perfectly. – ChrKoenig Jul 31 '14 at 09:29

score 0 · Answer 2 · answered Jul 25 '14 at 20:07

I must admit i'm not that familiar with partial residual plots so i'm not entirely sure what the proper interpretation of them should be given an interaction term. But basically, the equivalent of

crPlot(lm.fit1, "P1")

is

x <- predict(lm.fit1, type="term", term="P1")
y <- residuals(lm.fit1, type="partial")[,"P1"]
plot(x, y)
abline(lm(y~x), col="red", lty=2)
loessLine(x,y,col="green3",log.x = FALSE, log.y = FALSE, smoother.args=list())

so really, there's no real reason the same idea couldn't work with an interaction term as well. We just leave the partial contribution from a variable due to the interaction as a separate entity and just focus on the non-interaction contribution. So what i'm going to do is just take out the check for the interaction term and then we can use the function. Assuming that

body(car:::crPlot.lm)[[11]]
# if (any(attr(terms(model), "order") > 1)) {
#     stop("C+R plots not available for models with interactions.")
# }

we can copy and modify to create a new function with out the check

crPlot2 <- car:::crPlot.lm
body(crPlot2) <- body(crPlot2)[-11]
environment(crPlot2) <- asNamespace("car")

And then we can run

layout(matrix(1:2, ncol=2))
crPlot2(lm.fit2, "P1")
crPlot2(lm.fit2, "P2")

to get

enter image description here

I'm sure the authors had a good reason for not incorporating models with interaction terms so use this hack at your own risk. It's just unclear to me what should happen to the residual from the interaction term when making the plot.

Partial residual plots for linear model including an interaction term

2 Answers2