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 #2...how 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.

library(plotly)

library(dplyr)

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

attach(scatterdata)

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

newplot