How do I combine geometries in a shapefile based on a grouping variable?

You could also do this using dplyr's group_by() and summarize() functions:

states %>%
  group_by(group_var) %>% 
  summarize(geometry = st_union(geometry))

You were close; just cast the second argument in aggregate() as a list, like this:

x <- aggregate(states, 
               by = list(states$group_var),
               FUN = mean)

x
plot(x)