ggplot (in R)/matplotlib (in Python) with linear model?
I don't know what you mean by "with linear model" in the title, but here's code that generates a toy dataset and replicates your plot.

library(tidyverse)
x&lt;-crossing(year=paste("Year", 1:3), avg=c("Above 4.0", "Below 4.0"))
x$dat&lt;-replicate(6, tibble(wrkday=runif(100, 1000, 2000)))
x %&gt;% 
  unnest(dat) %&gt;% 
  ggplot(aes(dat)) + 
  geom_histogram(aes(y=..ncount..), bins=20) + 
  facet_grid(avg~year)




EDIT:

Given your comment I think I understand what you mean by "with linear model."  You want to see how wrkday varies as a function of year, faceted by avg.  Here's some code:

library(tidyverse)

x&lt;-tibble(year=sample(factor(paste("Year", 1:3)), 600, replace=T), 
          avg=sample(c("Above 4.0", "Below 4.0"), 600, replace=T)) %&gt;%
  mutate(wrkday=rnorm(600, mean=1000*(as.integer(year)-1)/2, sd=300))

ggplot(x, aes(x=as.integer(year), y=wrkday)) + 
  geom_jitter(width=0.2) + 
  geom_smooth(method="lm") + 
  facet_wrap(~avg)