How to plot interaction term in logistic regression in R

Question

I'm working on the logistic regression with interaction term between gender and education years as below. gender is a factor varibale with two value and eduyears2 is a numeric variable. and the depedndent variable is a dummy variable indicating if people have a job. How can I plot a line to visuallise the effects of education years on the dependent varibale in females and males

I'd be really appreciated to your help!

model <- glm(job~gender+eduyears2 + age2 + I(age2^2) + marriage +urban + time + gender*eduyears2, data = cfps1820, family = binomial)

Answer 1

I would generalize this to not assume age is quadratic and not assume years of education has a linear effect. The following fits both using restricted cubic spline functions, and plots predicted values by education with one curve for each gender, to show the amount of interactions. You can also use the contrast function to estimate gender differences as a function of education. The lrm function gives a lot of output that is special for logistic models. You may want to also allow time to be nonlinear.

require(rms)
# Compute good default plotting limits:
dd <- datadist(cfps1820); options(datadist='dd')
f <- lrm(job ~ gender * rcs(eduyears2, 4) + rcs(age, 4) + marriage + urban + time,
         data = cfps1820)
specs(f, long=TRUE)   # show spline knot locations etc.
ggplot(Predict(f, eduyears2, gender))

Answer 2

I would plot the education years --> job slopes separately for men and women. There are several ways to do this, in R perhaps the easiest is via SjPlot:

library(ggplot2)
library(SjPlot)
library(sjmisc)

p<-plot_model(model, type="pred", terms=c("eduyears2", "gender"))

#modify legend

p2<-p+scale_color_discrete(labels=c("Female", "Male"))

#Add axis titles and title

p_final<-p2+labs(title="Main title", x="Predictor variable", y="Dependent variable")