Robust Zonal stats in R and QGIS

QGIS and R's raster package use different methods to estimate zonal statistics. Briefly:

QGIS compares the centroid of each raster cell to the polygon boundary, initially considering cells to be wholly within or outside of the polygon based on the centroid. However, if fewer than two cell centroids fall within the polygon, an exact vector-based calculation is performed instead. That calculation is essentially the same as xunilk's answer.

R's raster package also uses a centroid test to determine if cells are inside or outside of the polygon. When weights=TRUE, it disaggregates the raster by a factor of 10 before performing the analysis, which reduces the error incurred by ignoring partially covered cells but reduces performance substantially.

If you're working in R, I would recommend the exactextractr package (I am the author). It provides a replacement for raster::extract, using an exact vector-based calculation that is hundreds of times faster than the original in most cases.

Here's an example using the test data provided:

library(sf)
library(raster)
library(exactextractr)

polygons <- st_read('testdata/Polygons.shp')
grid <- raster('testdata/Raster.tif')

polygons$result <- exact_extract(grid, polygons, fun=weighted.mean)

example <- polygons[c(1:4, 19:27), 'result']
st_geometry(example) <- NULL
print(example)

#       result
# 1  -63.95060
# 2  -67.00000
# 3  -67.00000
# 4  -62.37653
# 19 -99.94981
# 20 -92.79999
# 21 -92.79999
# 22 -98.04804
# 23 -91.76212
# 24 -91.76212
# 25 -87.48955
# 26 -86.00000
# 27 -86.00000

A Robust Zonal stats in QGIS can be implemented with PyQGIS. Following code was run with your layers by using a filter to select only ID_sub.pat between 19 and 27 (showed in your image).

import processing

registry = QgsMapLayerRegistry.instance()

Polygons = registry.mapLayersByName('Polygons')
request = QgsFeatureRequest().setFilterExpression (u'"ID_sub.pat" >= 19 AND "ID_sub.pat" <= 27')

iterator = Polygons[0].getFeatures(request)

feats_Polygons = [ feat for feat in iterator ]

Raster = registry.mapLayersByName('Raster')
rprovider = Raster[0].dataProvider()

xmin, ymin, xmax, ymax = Raster[0].extent().toRectF().getCoords()

extent = str(xmin)+ ',' + str(xmax)+ ',' +str(ymin)+ ',' +str(ymax) 

parameters = {"TYPE":1, 
              "EXTENT":extent, 
              "HSPACING":1000, 
              "VSPACING":1000, 
              "CRS":Raster[0].crs().authid(), 
              "OUTPUT":None}

path = processing.runalg('qgis:creategrid', 
                         parameters)

grid = QgsVectorLayer(path['OUTPUT'],
                      'grid',
                      'ogr')

feats_grid = [ feat for feat in grid.getFeatures() ]

for featg in feats_grid:
    for featp in feats_Polygons:
        if featg.geometry().intersects(featp.geometry()):
            geom = featg.geometry().intersection(featp.geometry())
            pt = geom.centroid().asPoint()
            value = rprovider.identify(pt,
                                       QgsRaster.IdentifyFormatValue).results()[1]

            print value, geom.area()/62500.0, featp.id()+1

After running code at Python Console of QGIS it was obtained this result:

-88.0 0.11367261864 22
-88.0 0.744772910144 25
-86.0 0.0389546011361 22
-86.0 0.152627219776 23
-86.0 0.152627219776 24
-86.0 0.255227089856 25
-86.0 1.0 26
-86.0 1.0 27
-102.399993896 0.744772910144 19
-102.399993896 0.631100291505 22
-92.799987793 0.255227089856 19
-92.799987793 1.0 20
-92.799987793 1.0 21
-92.799987793 0.21627248872 22
-92.799987793 0.847372780224 23
-92.799987793 0.847372780224 24  

where at first column are raster values, at second column weight factor and at third column id Polygons feature. You can observe that for polygons 20,21,26,27, its factors is always one (as expected). In other cases is proportional to intersection of Polygons feature with grid representing raster.

Results were corroborated with Value Tool and QuickWKT QGIS plugins and they were as expected.