How do I plot the point estimates and the fixed effects for my LPM?

Question

My dataset looks like this:

birth_year = c(1945, 1950, 1960)
work = c(1,1,0)
city = c("NY", "London", "Paris")
education = c(5, 6, 7)

df <- data.frame(birth_year, work, city, education)

I have the following lm and would like to plot the fixed effects coefficients and standard errors for the birth_year variable using the sjPlot package.

model <- lm(work ~ education + as.factor(city) + as.factor(birth_year), data = df)

I have tried

plot_model(lm(work ~ as.factor(birth_year) + as.factor(city) + education. df, terms = "as.factor(birth_year)"))

But this is the error that I am getting:

ror in ans[ypos] <- rep(yes, length.out = len)[ypos] : 
  replacement has length zero
In addition: Warning messages:
1: In min(new_value, na.rm = TRUE) :
  no non-missing arguments to min; returning Inf
2: In min(dat$conf.low) : no non-missing arguments to min; returning Inf
3: In min(dat$estimate) : no non-missing arguments to min; returning Inf
4: In max(dat$conf.high) : no non-missing arguments to max; returning -Inf
5: In max(dat$estimate) : no non-missing arguments to max; returning -Inf
6: In min(ticks) : no non-missing arguments to min; returning Inf
7: In max(ticks) : no non-missing arguments to max; returning -Inf
8: In rep(yes, length.out = len) : 'x' is NULL so the result will be NULL