Addressing your exact output and expanding on the code you provided, you could query by date (the advantage of Posix) and index the rasters in a stack object. The use of "which" returns the raster stack/brick index meeting the query criteria. The use of "apply" is to sum rows, which represent pixels, to derive a seasonal sum to average. For large rasters that are not in memory you may need to use getValues.
library(raster)
r <- raster(ncol=10, nrow=10)
r[] <- runif(ncell(r))
s <- stack(r,r,r,r,r,r,r,r,r,r)
s <- stack(lapply(1:10, function(x) setValues(r, runif(ncell(r)))))
s <- setZ(s, as.Date('2000-1-1') + 0:9)
mean( apply(s[which(getZ(s) < "2000-01-04")], MARGIN=1, FUN=sum) )
mean( apply(s[which(getZ(s) >= "2000-01-04")], MARGIN=1, FUN=sum) )
# data.frame output
( s.mean <- data.frame(
s1.mean = mean( apply(s[which( getZ(s) < "2000-01-04")], MARGIN=2, FUN=sum) ),
s2.mean = mean( apply(s[which( getZ(s) <= "2000-01-04")], MARGIN=2, FUN=sum) ) ) )
I will say that this is not a valid way to summarize the data and you are collapsing a considerable amount of variation in to a single number but, I will leave it to you to figure out. However, please talk with somebody that is statistically savvy about the effect of skewed distributions on means and better ways to summarize spatially variable data.