Combine a polygon map on top of a raster map in R

Question

I would like to combine a shapefile map of some US states on top of a raster map obtained from using library ggmap.

This is the code I have

library(USAboundaries)
library(ggmap)
library(ggplot2)
library(sf)

us.state.longlat <- us_states()
us.state.longlat.cropped.large <- sf::st_crop(us.state.longlat, 
                                              c(xmin=-100, xmax=-80, ymin=30, ymax=50))

bbox <- setNames(st_bbox(us.state.longlat.cropped.large), 
                 c("left", "bottom", "right", "top"))

ggmap(get_map(bbox, source = "stamen", zoom = 5)) +
  geom_sf(data = us.state.longlat.cropped.large,
          inherit.aes = FALSE,
          color = "red", 
          fill = "transparent")

This is what I got:

Also got the following warning that I am not sure what it means:

Coordinate system already present. Adding new coordinate system, which will replace the existing one.

As you see there is a mismatch between the US state shapefile and the raster map. The coordinate reference system for the shapefile map us.state.longlat.cropped.large is this:

st_crs(us.state.longlat.cropped.large)
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

The issue seems to be less worrisome when cropping a small area map. For example:

us.state.longlat.cropped.small <- sf::st_crop(us.state.longlat, c(xmin=-92, xmax=-85, ymin=32, ymax=36))
bbox <- setNames(st_bbox(us.state.longlat.cropped.small), 
                 c("left", "bottom", "right", "top"))

ggmap(get_map(bbox, source = "stamen", zoom = 7)) +
  geom_sf(data = us.state.longlat.cropped.small,
          inherit.aes = FALSE,
          color = "red", 
          fill = "transparent")

This is what I got for a reduced area of the map:

In any case, I don't know what I am missing here but in both cases state borders should match between the shapefile and the raster map.