Построение графиков для множественной линейной модели, включающей непрерывные и дискретные предикторы, в R с помощью ggplot2

Делаю проект для онлайн-курса по статистике в R. Исходный датасет - “NELS88” из пакета “copulaData”. В соответствии с заданием и после некоторой валидации работаю с датафреймом “urpub2” и моделью “Model10”, включающей как непрерывные, так и дискретные предикторы (факторы), а также одно взаимодействие. Теперь пытаюсь построить три разных графика, отражающих связь зависимой переменной с предикторами, а также на одном из них визуализировать взаимодействие. В курсе учили изображать связь с наиболее важными предикторами при условии, что все остальные предикторы неизменны, например, их значения равны средним. Только как сделать неизменными факторы? Получается только "разбить" модель на три, что, как мне кажется, не совсем правильно. Можно ли добиться подобных нижеприведенным графиков, работая с одной моделью (Model10)? Код также прикрепляю.

library(copulaData)
library(ggplot2)
library(cowplot)
library(car)
library(dplyr)
data('NELS88')
urpub <- NELS88[NELS88$Urban == 1 & NELS88$Public == 1 & NELS88$Size > 100,]
urpub$Size_factor <- factor(urpub$Size)
urpub2 <- urpub[-19,]
Model10 <- lm(Science ~ Math + Reading + SES + Minority + Female + SES:Female, data = urpub2)

График для модели "Science ~ Math + Reading"

theme_set(theme_bw())
# Science ~ Math + Reading
Incorrect_M1 <- lm(Science ~ Math + Reading, data = urpub2)
new_data1 <- expand.grid(Math = seq(min(urpub2$Math), max(urpub2$Math), length.out = 100), Reading = seq(min(urpub2$Reading), max(urpub2$Reading), length.out = 100))
new_data1$predicted <- predict(Incorrect_M1, newdata = new_data1)
inc_Pl1 <- ggplot(new_data1, aes(x = Math, y = predicted, group = Reading)) + geom_line(aes(color = Reading)) + geom_point(data = urpub2, aes(x = Math, y = Science)) + scale_color_continuous(high = 'red', low = 'yellow') + ylab('Science')

График для модели "Science ~ Minority"

# Science ~ Minority
Incorrect_M2 <- lm(Science ~ Minority, data = urpub2)
new_data2 <- data.frame(Minority = factor(levels(urpub2$Minority), levels = levels(urpub2$Minority)))
new_data2$fit <- predict(Incorrect_M2, newdata = new_data2, se.fit = TRUE)$fit
new_data2$se <- predict(Incorrect_M2, newdata = new_data2, se.fit = TRUE)$se.fit
t_crit2 <- qt(0.975, df = nrow(urpub2) - length(coef(Incorrect_M2)))
new_data2$lwr <- new_data2$fit - t_crit2 * new_data2$se
new_data2$upr <- new_data2$fit + t_crit2 * new_data2$se
inc_Pl2 <- ggplot(new_data2, aes(x = Minority, y = fit)) + geom_bar(stat = 'identity', aes(fill = Minority)) + geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.2) + scale_x_discrete(labels = c('No', 'Yes')) + guides(fill = 'none') + labs(x = 'Minority', y = 'Science')

График для модели "Science ~ SES + Female + SES:Female"

# Science ~ SES + Female + SES:Female
Incorrect_M3 <- lm(Science ~ SES + Female + SES:Female, data = urpub2)
new_data3 <- urpub2 %>% group_by(Female) %>% do(data.frame(SES = seq(min(.$SES), max(.$SES), length.out = 100)))
Predictions3 <- predict(Incorrect_M3, newdata = new_data3, se.fit = TRUE)
new_data3$fit <- Predictions3$fit
new_data3$se <- Predictions3$se.fit
t_crit3 <- qt(0.975, df = nrow(urpub2) - length(coef(Incorrect_M3)))
new_data3$lwr <- new_data3$fit - t_crit3 * new_data3$se
new_data3$upr <- new_data3$fit + t_crit3 * new_data3$se
inc_Pl3 <- ggplot(new_data3, aes(x = SES, y = fit)) + geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr, group = Female))  + geom_line(aes(colour = Female), size = 1) + geom_point(data = urpub2, aes(x = SES, y = Science, colour = Female)) + scale_colour_discrete(name = 'Sex', labels = c('Male', 'Female')) + labs(y = 'Science')