Plot regression adjusted for covariate I would like to plot values from a linear regression adjusted for a covariate.
Please consider the following dummy data in which y is predicted by x and the covariate a. There is an interaction between x and a.
# generate some dummy data
x1 <- rnorm(20,10,5)
x2 <- rnorm(20,20,5)

a1 <- rnorm(20,10,5)
a2 <- rnorm(20,20,5)

y1 <- x1*a1+rnorm(20,5,3)
y2 <- x2*a2+rnorm(20,5,3)

x <- c(x1,x2)
y <- c(y1,y2)
a <- c(a1,a2)

# create model with & without covariate
m1 <- lm(y~x*a)
m2 <- lm(y~a+a:x)
m3 <- lm(m2$residuals~x)

summary(m1)
summary(m2)
summary(m3)

# plot data
par(mfrow=c(1,2))
plot(x,y)
abline(lm(y~x, subset=c(1:20)))
abline(lm(y~x, subset=c(21:40)))
plot(x, m2$residuals)
abline(m3)



*

*Is this an appropriate way to visualize the data? 

*Any other ideas how to visualize a main effect when interactions are present?

*Why is in model m1 a main effect of x but when I analyze the residuals of model m2 in model m3 there is no effect of x? 

*Is that as inclusion of different predictors will always change the significance of other predictors?

 A: How about this:
# generate some dummy data
x1 <- rnorm(50,10,5)
x2 <- c(rep(0, 25), rep(1,25))
y <- x1*2 + x2*5 + x1*x2*2 + rnorm(50,0,10)

# create model with & without covariate

m1 <- lm(y ~ x1 + x2 + x1*x2)
summary(m1)

#plot
plot(y~x1, col = (x2+1))
abline(a = coef(m1)[1], b = coef(m1)[2])
abline(a = coef(m1)[1]+coef(m1)[3], b = coef(m1)[2] +
coef(m1)[4], col = 'red')

A: There is no interaction between x and a because the information is already in x about the dummy (0,1) variable - there is a step change in the mean of x. Perhaps you want all your x's to be the same? Combining m2 and m3 should definitely give the same results as m1.
A: It's a little hard to answer these questions because they are somewhat confused.  (I don't mean to be nasty about that, but sometimes you need to get clearer before you can ask the question that will get you what you need to know.)  One issue is with your example situation.  In your setup, x and a are confounded, and the true data generating process makes y a function of the x*a interaction only.  I will attempt to provide some general information that may be helpful for you:  


*

*Typically, to visualize an interaction between two continuous covariates, you might just plot y as a function of x at the mean of a, and at (say) the mean $\pm$1 SD of a.  If one of the variables in the interaction is categorical, you can plot them as @PeterFlom demonstrates.  

*Other possibilities for visualizing interactions include coplots, and interactive graphics.  

*You did not set the seed (see ?set.seed) for your example; when I ran it, I got different data than you did and did not find a main effect.  This is not too surprising, since there is no main effect in the data generating process.  That is, if you did find one earlier, it was a type I error.  

*The inclusion of more predictors will change the level of 'significance' of your variables by virtue of: consuming another degree of freedom, absorbing some of the (previous) residual variance, and, especially, if any of your covariates are correlated with each other (as yours are).  These topics are discussed extensively on CV; you can read around under multiple-regression and multicollinearity.  

A: To visualize the effects of interactions and changing variable values in a regression you may find the Predict.Plot and TkPredict functions in the TeachingDemos package for R useful.
These functions will plot the prediction line/curve relating one of the predictors to the response variable for given values of the other variables.  The TkPredict function then lets you interactively change the values of the other variables to see how that changes the relationship.  You can use Predict.Plot to overlay a few different relationships based on the values of the other variables to compare in a static plot.
