Averaging variables of smaller polygons to larger polygons in R using the SF package

Still getting used to the language of GIS with R so please bear with me.

I have this public dataset.

It has two levels of political boundaries: ADM1 and ADM2.

I am currently doing everything (i.e. plotting, calculations, etc.) for ADM2 which is each of the 2,440+ municipalities in Mexico because that is how the geometry is determined. However, I would like to calculate averages per states (ADM1) and start working with those by sort of "dissolving", if you will, the geometry of municipalities (ADM2) into geometries of (ADM1).

So far I have imported it:

library (sf)
library (dplyr)

mexi_sf<-read_sf( dsn = getwd()
, layer = "MexicoCaseStudy"
, stringsAsFactors = FALSE )

Fixed some labels and changed some classes:

c("M.xico" = "México", "Nuevo Le.n" = "Nuevo León","San Luis Potos."="San Luis Potosí","Quer.taro de Arteaga"="Querétaro de Arteaga"))

Changed some outliers:

for (row in 1:nrow(mexi_sf)){
if (mexi_sf\$FOODEXP[row]==-9999){
mexi_sf\$FOODEXP[row]=-1
}}

Here is my attempt to average the variable FOODEXP by ADM1NAME.

average_food_exp<-mexi_sf%>%
mutate(meanFOODEXP=mean(FOODEXP))

And obvs something is wrong because I can't plot it:

big_food_map_sf_trial<-tm_shape(average_food_exp,style="quantiles",n=8)+
tm_fill(col="meanFOODEXP",palette="Reds",title="mean FOODEXP (MXN 2004)")+
tm_layout(legend.outside =FALSE)+
tm_borders(col="black")+
tm_credits("Source: CIMMYT, June 2004, Mexican municipalities below food poverty line, year 2000.",position=c("left","bottom"),size=0.75)

save_tmap(tm=big_food_map_sf_trial,filename="big_food_map_sf_trial.png",width = 17,height=7)

Any ideas on how to proceed?

I think you need summarize(meanFOODEXP=mean(FOODEXP)) (not mutate) after the group_by.

E.g.

library(sf)