CODEWITHSUNDEEP
rplotstatistics
ChrKoenig
ChrKoenig

Reputation: 973

Partial residual plots for linear model including an interaction term

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

Upvotes: 4

Views: 11519

Answers (2)

jlhoward
jlhoward

Reputation: 59395

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

Upvotes: 8

MrFlick
MrFlick

Reputation: 206411

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.

Upvotes: 0

Related Questions