Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Calculation of Moran's I with 4000 records

I have 4000 records of volume in trees plantation. I need to calculate the Moran's I to the whole plantation. I use ape library because spdep is said to be slower. My code is this:

# Modified from http://www.ats.ucla.edu/stat/r/faq/morans_i.htm
require(ape)
df <- data.frame(
     x = 1:2000,
     y = 1:2000,
     v = rnorm(4000, mean=4) )
df.dists <- as.matrix(dist(cbind(df$x, df$y)))
df.dists.inv <- 1/df.dists
diag(df.dists.inv) <- 0
Moran.I(df$v, df.dists.inv)

When I run the code, I get overflow-like errors.

*Error in if (obs <= ei) 2 * pv else 2 * (1 - pv) : 
  missing value where TRUE/FALSE needed*

Using ff library

require(ape)
require(ff)
ffdf <- as.ffdf(df)
ffdf.dists <- as.matrix(dist(cbind(ffdf$x, ffdf$y)))
ffdf.dists.inv <- 1/df.dists
diag(ffdf.dists.inv) <- 0
Moran.I(ffdf$v, ffdf.dists.inv)

More error messages:

*Error in x - m : non-numeric argument to binary operator
In addition: Warning message:
In mean.default(x) : argument is not numeric or logical: returning NA*
  • How can I get the calculation to the whole plantation?

  • Should I use sdep instead of ape library?

  • How could parallel library solve this problem?

Thanks in advance Juan

like image 304
jlopez Avatar asked Sep 21 '13 11:09

jlopez


People also ask

How to calculate Moran's I?

Calculations for Moran's I are based on a weighted matrix, with units i and j. Similarities between units is calculated as the product of the differences between yi and yj with the overall mean. The Moran's statistic is calculated using the basic form, which is divided by the sample variance:s2 = (Σ(yi– ̄y)2)/n).


1 Answers

You have some infinite values in your matrix. This should work in the 2 cases ( with and without ff package)

df.dists.inv[is.infinite(df.dists.inv)] <- 0

Applying this with a small example:

require(ape)
set.seed(1)
df <- data.frame(
  x = 1:10,
  y = 1:10,
  v = rnorm(20, mean=4) )
.....

df.dists.inv[is.infinite(df.dists.inv)] <- 0
Moran.I(df$v, df.dists.inv)

$observed
[1] -0.02246154

$expected
[1] -0.05263158

$sd
[1] 0.05399303

$p.value
[1] 0.5763143
like image 114
agstudy Avatar answered Oct 03 '22 05:10

agstudy