You can see the problem in the figure at the bottom of this question, but the output of gstat::idw for a stars input is a raster that doesn't line up with the input.

library(sf)
library(stars)
library(dplyr)
library(gstat)
library(ggplot2)


# Sample data
ex_points <- structure(list(
  lat = c(28.2911345, 28.2908078, 28.2910919, 28.2910929, 
          28.2906787, 28.2910489, 28.2908434, 28.2910878, 28.2910093, 28.2909689, 
          28.2906751, 28.2909233, 28.2908463, 28.2907666, 28.2912162, 28.2912135, 
          28.2912149, 28.2909714, 28.2910123, 28.2911411, 28.2909244, 28.29101, 
          28.2910492, 28.2908855, 28.290968, 28.2911294, 28.2911379, 28.2911784, 
          28.2911367, 28.2911737, 28.2910918, 28.2908019, 28.2910901, 28.290552, 
          28.2910951, 28.2911781, 28.2912249, 28.2906767, 28.290928, 28.2911758, 
          28.290637, 28.2906335, 28.2905937, 28.2906358, 28.2909653, 28.2911365, 
          28.2912134, 28.29122, 28.2910532, 28.291013, 28.2911337, 28.2911754, 
          28.2908402, 28.2908868, 28.2909677, 28.2910096, 28.2908444, 28.2908839, 
          28.2908062, 28.290758, 28.2907168, 28.2907154, 28.2907574, 28.2908823, 
          28.2911706, 28.291097, 28.2910515, 28.2910502, 28.2911762, 28.2911745, 
          28.2911319, 28.290596, 28.2905628, 28.2909705, 28.2912162), 
  lng = c(-112.4276232, 
          -112.4274911, -112.4274911, -112.4274213, -112.4276918, -112.4276217, 
          -112.4276875, -112.4276912, -112.4276255, -112.427619, -112.4276227, 
          -112.4276208, -112.427488, -112.4274902, -112.427423, -112.4275562, 
          -112.42769, -112.427356, -112.4273548, -112.4272236, -112.4276879, 
          -112.4276898, -112.4276892, -112.4274888, -112.4274902, -112.4275539, 
          -112.4273534, -112.4272874, -112.4272865, -112.427422, -112.4276212, 
          -112.4275576, -112.4275562, -112.4276889, -112.4272856, -112.4272232, 
          -112.4272207, -112.4275594, -112.4275579, -112.4273513, -112.4275586, 
          -112.4276227, -112.4276918, -112.4276929, -112.4276902, -112.4274217, 
          -112.4276253, -112.4273509, -112.4273544, -112.4274211, -112.4274899, 
          -112.4274853, -112.4275533, -112.4275553, -112.4275558, -112.4275565, 
          -112.4276235, -112.427625, -112.4276216, -112.4276247, -112.4276225, 
          -112.4276897, -112.427689, -112.4276905, -112.4276869, -112.4273561, 
          -112.4274207, -112.4275583, -112.4275556, -112.4276232, -112.4276855, 
          -112.4275558, -112.4276268, -112.42742, -112.4274891), 
  metric = c(0.744272, 
             0.724683, 0.727251, 0.720583, 0.751171, 0.747609, 0.723744, 0.720821, 
             0.72594, 0.68184, 0.644309, 0.7076, 0.787405, 0.673647, 0.711521, 
             0.738893, 0.735662, 0.728458, 0.717527, 0.704002, 0.750679, 0.72177, 
             0.760431, 0.734251, 0.729909, 0.762904, 0.682811, 0.714333, 0.725169, 
             0.715958, 0.711962, 0.775083, 0.739643, 0.686923, 0.721321, 0.545518, 
             0.724519, 0.657674, 0.750514, 0.679069, 0.68413, 0.70484, 0.739893, 
             0.660236, 0.73532, 0.66926, 0.720944, 0.687076, 0.66998, 0.741386, 
             0.741826, 0.686484, 0.74844, 0.673051, 0.605934, 0.760162, 0.661219, 
             0.738082, 0.672051, 0.731661, 0.552748, 0.740665, 0.658713, 0.75776, 
             0.753953, 0.752685, 0.709552, 0.698544, 0.701351, 0.759725, 0.731301, 
             0.628379, 0.611095, 0.754333, 0.733381)), 
  row.names = c(NA, -75L), 
  class = c("tbl_df", "tbl", "data.frame")
)

# Convert to sf
ex_sf <- 
  st_as_sf(ex_points, coords = c("lng", "lat"), remove = FALSE, crs = 4326)

# Construct raster grid
ex_grd <- st_bbox(ex_sf) %>% 
  st_as_stars(nx = 30, ny = 30) %>% 
  st_crop(
    st_convex_hull(st_combine(ex_sf))
  )

# Everything lines up as expected:
ggplot() +
  geom_stars(data = ex_grd) +
  geom_sf(data = ex_sf) +
  theme(axis.text.x = element_text(angle = 90))

# Fit a simple model
ex_model <- gstat::idw(metric~1, ex_sf, ex_grd)
#> [inverse distance weighted interpolation]

# Now points and raster are mis-aligned in the y-direction
ggplot() +
  geom_stars(data = ex_model) +
  geom_sf(data = ex_sf, aes(fill = metric), color = "black", shape = 21) +
  theme(axis.text.x = element_text(angle = 90)) +
  scale_fill_viridis_c()

^{Created on 2025-04-03 with reprex v2.0.2}

The points and raster predictions appear to be aligned vertically but offset by about 0.002º. This is dummy data at a similar latitude to my real points. When using my real data in a subset I get this offset. When I use my real data with my much larger dataset, I get a larger offset, but offset the other direction vertically (the raster appears north of the points by 0.01º instead of south).