Why covariance values are different in R and ArcGIS?

Question

I've calculated covariance between 7 rasters using R and ArcGIS. It gives me different results. Each times covariance values from R are approximately 1.874 times higher than ArcGIS values.

R - package raster, function layerStats(stack, 'cov', na.rm=TRUE)

ArcGIS - Band Collection Statistics tool (http://resources.arcgis.com/en/help/main/10.1/index.html#//009z000000q3000000)

Why covariance values are different in R and ArcGIS?

Answer 1

Create some simple test rasters in R:

> m=matrix(1:9,3,3)
> m2 = matrix(c(9,2,3,4,1,5,6,8,7),3,3)

Then we can trivially compute the covariance between these matrices:

> cov(c(m),c(m2))
[1] 2.125

and I would wager doing the computation by hand would get the same result. What does layerStats do?

> D = stack(raster(m),raster(m2))

> layerStats(D, "cov", na.rm=TRUE)
$covariance
        layer.1 layer.2
layer.1   7.500   2.125
layer.2   2.125   7.500

There's the same 2.125 in the cross-correlation.

Now try with an NA in there:

> m[2,2]=NA
> cov(c(m),c(m2),use="complete.obs")
[1] 2.428571
> D = stack(raster(m),raster(m2))
> layerStats(D, "cov", na.rm=TRUE)
$covariance
         layer.1  layer.2
layer.1 8.571429 2.428571
layer.2 2.428571 7.500000

Again, cov agrees with layerStats if we remove the NA value.

So what does ArcGIS do in each of these situations? Save the D, import to ArcGIS, and find out...

EDIT

The results from ArcGIS:

Raster m and m2:

Raster with an NA and m2:

[Now, the question is - Which way of calculating covariance between rasters is correct?]