Welcome back to the forum.
There is a procedure (already implemented) which enables calculations of the mean value of a target variable over a polygon.
I must first define what is a POLYGON in RGeoS. A POLYGON is a set of POLYSET(s). The number of polysets in a polygon is simple positive.
Then for each polyset, the mean value of the variable is calculated using a global Kriging (rather than a block kriging). In particular, it also gives access to the variance of the estimation error. As far as the data are concerned, it offers the possibility to select all the samples included inside the polyset, or to use all available data instead.
Finally this global estimation is carried out using a discretization of the polygon for the calculation of the covariance terms. Thsi discretization can be performed regularly (discertization points located on the nodes of a regular grid, with a mesh fixed by the user) or a stratified grid (where samples are generated randomly within cells of a regular grid defined by the user).
Hope this will serve.