I am new to GIS in R. I have two input rasters (an elevation DEM and a land use classified raster), and I am tasked with reprojecting, aggregating into courser resolution, and finally calculating the average land slope for each land use classification. I am having trouble finding the average slope for the land uses; this is what I need help with.

For ease assisting with this question, I have provided both of the rasters on a Google Drive for users to download. The elevation raster DEM can be downloaded here, whereas the land use classification raster can be downloaded here.

Here is what I have done so far:

#Load libraries
library(raster)
library(sf)

#Load data
cdl19 <- raster("CDL_2019_clip_20220329135051_1368678123.tif")
elev <- raster("elevation.tif")

#View information about the rasters
plot(cdl19)
hist(cdl19)

plot(elev)
hist(elev)

#Check coordinate reference systems
crs(cdl19)
crs(elev)

#Reproject each raster so that they will be in the same grid system
cor_ref = '+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'

cdl19_proj <- projectRaster(cdl19, crs = cor_ref, method = 'ngb') #Use 'ngb' for categorical data

elev_proj <- projectRaster(elev, crs = cor_ref, method = 'bilinear') #use 'bilinear' for continuous data

#The references should now be the same
#####################################################

#Resample to eliminate resolution problems
elev_resamp <- resample(elev_proj, cdl19_proj, method = 'bilinear')

#Aggregate to larger, 1km x 1km cell size
#Note: This is not exactly 1km x 1km, but projecting to a finer coordinate system would overwhelm my computer's memory capabilities.

cdl19_1k <- aggregate(cdl19_proj, fact = 33.333333333, fun = modal, na.rm=TRUE)

elev_1k <- aggregate(elev_resamp, fact = 33.333333333, fun = modal, na.rm=TRUE)


#Calculate slope of DEM
slope <- terrain(elev_1k, opt = 'slope', units = 'degrees')

#Crop and mask categorical raster to the DEM
cdl19_1k_crop <- crop(cdl19_1k, extent(slope))
mask <- mask(cdl19_1k_crop, mask = slope)

Here is where I need assistance. Calling the freq(mask) shows that there are 17 land use values (not including NA values). The output can be seen below:

      value count
 [1,]     1   137
 [2,]     4     9
 [3,]     5   230
 [4,]    24    16
 [5,]    28     1
 [6,]    36     7
 [7,]    37     1
 [8,]    61    13
 [9,]   111   136
[10,]   121     7
[11,]   122    52
[12,]   123     7
[13,]   124     2
[14,]   141   351
[15,]   143     1
[16,]   176  2021
[17,]   195    22
[18,]    NA  1342

How can I find the average slope of each land use value? I would presume that the solution involves combining the rasters and running a cellStats operation, but I'm not entirely sure.

Thanks in advance!