Aggregate polygons based on area (weighted average)

Question

I’m working with geospatial data on districts and neighborhoods. I have a dataset with multiple features for each neighborhood. Now I want to aggregate these neighborhoods into the bigger districts, for this I have only the geometry. I tried the aggregate function from the SF library in R, but this takes the average of all polygons that intersects with the overlay polygon. Since some polygons only intersects partly I want to calculate the mean based on the area.

I have created an example below:

m1 = cbind(c(0, 0, 10, 0), c(0, 10, 10, 0))
m2 = cbind(c(0, 10, 10, 0), c(0, 0, 10, 0))

m3 = cbind(c(0, 0, 1, 1, 0), c(0, 1, 1, 0, 0))
m4 = cbind(c(1, 2, 2, 1, 1), c(9, 9, 8, 8, 9))
m5 = cbind(c(9, 9, 8, 8, 9), c(1, 2, 2, 1, 1))
pol = st_sfc(st_polygon(list(m1)), st_polygon(list(m2)))
pol2 = st_sfc(st_polygon(list(m3)), st_polygon(list(m4)), st_polygon(list(m5)))

plot(pol)
plot(pol2, add = TRUE)

d = st_sf(data.frame(value1=c(5, 8 ,7), value2=c(1, 4 , 9), geom=pol2))


(p_ag1 = aggregate(d, pol, mean))
#   value1 value2                       geometry
# 1    6.5    2.5 POLYGON ((0 0, 0 10, 10 10,...
# 2    6.0    5.0 POLYGON ((0 0, 10 0, 10 10,...

# Expected weigthed average based on area
#   value1 value2                       geometry
# 1    7.0    3.0   POLYGON ((0 0, 0 10, 10 10,...
# 2    6.33   6.33  POLYGON ((0 0, 10 0, 10 10,...

Can someone help me out?