Regression visualization

Hi all,
Trying to figure out the best way to visualize a particular set of findings and hoping someone more knowledgeable than I can weigh in as I am way out of my depth here. We have a multiple regression equation that includes covariates of no interest and some interactions. For example's sakes, let's say: Y = Intercept + Categorical + NuisanceVar + InterestingVar + NuisanceVar*InterestingVar + Categorical*InterestingVar.

The latter is the main outcome of interest, essentially showing that the category (Drug A versus Drug B) reverses the relationship between InterestingVar and Y. The analysis was relatively straightforward, but I'm struggling to figure out the best way to visualize these findings. Its somewhat apparent from looking at a simple scatterplot of Y and InterestingVar, but the regression line for the full model is gnarly (i.e. jagged) because NuisanceVar adds an additional dimension that isn't in the visualiztion. I "think" I want a partial residual plot, but I've never done one and have only seen these looking at single IVs. So essentially my question is:
1) Am I being a dunce and is there an obvious easy to do this I'm just not thinking of?
2) If no, Is it kosher to essentially create a multiple-partial-residual plot that distills the above to Y = Categorical + Categorical*InterestingVar so we can see the relationship?
3) If we go with would one go about doing this?

For the R aficionados (of which I am not, but I'm trying to make the switch from SPSS!), I'm coping some sample code below to illustrate the issue.

Thanks to any advice anyone can offer.


catvar <- c(rep("Drug A", 25), rep("Drug B", 25))
xvals <- rnorm(50, mean=10, sd=2)
yvals <- rnorm (50, mean=20, sd=3)
covar <- rnorm(50, mean = 0, sd = 5)

scatterdata <- data.frame("CatVar"=catvar, "XVals"=xvals, "YVals" = yvals, "Covar" = covar)
colorpal <- c("red", "black")

sorteddata <- scatterdata[order(xvals, yvals),]

##Original plot
linefit <- lm(yvals~Covar+XVals+CatVar+Covar*XVals+CatVar*XVals, data=sorteddata)
fitvals <- data.frame("YVals" = fitted(linefit), "CatVar" = sorteddata$CatVar, "XVals" = sorteddata$XVals)
DrugAfitvals <- fitvals[which(sorteddata$CatVar=='Drug A'),]
DrugBfitvals <- fitvals[which(sorteddata$CatVar=='Drug B'),]

newplot <- plot_ly() %>%
add_trace(data=sorteddata, x= ~XVals, y = ~YVals, color = ~CatVar,
colors=colorpal, type = "scatter", mode = "markers", marker = list(size=10)) %>%
add_trace(data = DrugAfitvals,
x = ~XVals,
y = ~YVals, mode = "lines") %>%
add_trace(data = DrugBfitvals,
x = ~XVals,
y = ~YVals, mode = "lines")