1 ## MoranI.R (2008-01-14)
3 ## Moran's I Autocorrelation Index
5 ## Copyright 2004 Julien Dutheil, 2007-2008 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 ## code cleaned-up by EP (Dec. 2007)
12 Moran.I <- function(x, weight, scaled = FALSE, na.rm = FALSE,
13 alternative = "two.sided")
15 if(dim(weight)[1] != dim(weight)[2])
16 stop("'weight' must be a square matrix")
18 if(dim(weight)[1] != n)
19 stop("'weight' must have as many rows as observations in 'x'")
28 weight <- weight[!nas, !nas]
30 warning("'x' has missing values: maybe you wanted to set na.rm=TRUE?")
31 return(list(observed = NA, expected = ei, sd = NA, p.value = NA))
35 ## normaling the weights:
36 ## Note that we normalize after possibly removing the
38 ROWSUM <- rowSums(weight)
39 ## the following is useful if an observation has no "neighbour":
40 ROWSUM[ROWSUM == 0] <- 1
41 weight <- weight/ROWSUM # ROWSUM is properly recycled
45 y <- x - m # centre the x's
46 cv <- sum(weight * y %o% y)
51 i.max <- (n/s) * (sd(rowSums(weight) * y)/sqrt(v/(n - 1)))
55 S1 <- 0.5 * sum((weight + t(weight))^2)
56 S2 <- sum((apply(weight, 1, sum) + apply(weight, 2, sum))^2)
57 ## the above is the same than:
60 ## S2 <- S2 + (sum(weight[i, ]) + sum(weight[, i]))^2
63 k <- (sum(y^4)/n) / (v/n)^2
64 sdi <- sqrt((n*((n^2 - 3*n + 3)*S1 - n*S2 + 3*s.sq) -
65 k*(n*(n - 1)*S1 - 2*n*S2 + 6*s.sq))/
66 ((n - 1)*(n - 2)*(n - 3)*s.sq) - 1/((n - 1)^2))
68 alternative <- match.arg(alternative, c("two.sided", "less", "greater"))
69 pv <- pnorm(obs, mean = ei, sd = sdi)
70 if (alternative == "two.sided")
71 pv <- if (obs <= ei) 2*pv else 2*(1 - pv)
72 if (alternative == "greater") pv <- 1 - pv
73 list(observed = obs, expected = ei, sd = sdi, p.value = pv)
76 weight.taxo <- function(x)
78 d <- outer(x, x, "==")
79 diag(d) <- 0 # implicitly converts 'd' into numeric
83 weight.taxo2 <- function(x, y)
85 d <- outer(x, x, "==") & outer(y, y, "!=")
90 correlogram.formula <- function(formula, data = NULL, use = "all.obs")
92 err <- 'formula must be of the form "y1+...+yn ~ x1/x2/../xn"'
93 use <- match.arg(use, c("all.obs", "complete.obs", "pairwise.complete.obs"))
94 if (formula[[1]] != "~") stop(err)
97 y.nms <- if (length(lhs) > 1)
98 unlist(strsplit(as.character(as.expression(lhs)), " \\+ "))
99 else as.character(as.expression(lhs))
102 gr.nms <- if (length(rhs) > 1)
103 rev(unlist(strsplit(as.character(as.expression(rhs)), "/")))
104 else as.character(as.expression(rhs))
107 ## we 'get' the variables in the .GlobalEnv:
108 y <- as.data.frame(sapply(y.nms, get))
109 gr <- as.data.frame(sapply(gr.nms, get))
114 if (use == "all.obs") {
118 if (use == "complete.obs") {
119 sel <- complete.cases(y, gr)
123 na.rm <- use == "pairwise.complete.obs"
125 foo <- function(x, gr, na.rm) {
126 res <- data.frame(obs = NA, p.values = NA, labels = colnames(gr))
127 for (i in 1:length(gr)) {
128 sel <- if (na.rm) !is.na(x) & !is.na(gr[, i]) else TRUE
131 w <- if (i > 1) weight.taxo2(g, gr[sel, i - 1]) else weight.taxo(g)
132 o <- Moran.I(xx, w, scaled = TRUE)
133 res[i, 1] <- o$observed
134 res[i, 2] <- o$p.value
136 ## We need to specify the two classes; if we specify
137 ## only "correlogram", 'res' is coerced as a list
138 ## (data frames are of class "data.frame" and mode "list")
139 structure(res, class = c("correlogram", "data.frame"))
142 if (length(y) == 1) foo(y[[1]], gr, na.rm)
143 else structure(lapply(y, foo, gr = gr, na.rm = na.rm),
144 names = y.nms, class = "correlogramList")
148 function(x, legend = TRUE, test.level = 0.05,
149 col = c("grey", "red"), type = "b", xlab = "",
150 ylab = "Moran's I", pch = 21, cex = 2, ...)
152 BG <- col[(x$p.values < test.level) + 1]
153 if (pch > 20 && pch < 26) {
160 plot(1:length(x$obs), x$obs, type = type, xaxt = "n", xlab = xlab,
161 ylab = ylab, col = CO, bg = BG, pch = pch, cex = cex, ...)
162 axis(1, at = 1:length(x$obs), labels = x$labels)
164 legend("top", legend = paste(c("P >=", "P <"), test.level),
165 pch = pch, col = col, pt.bg = bg, pt.cex = cex, horiz = TRUE)
168 plot.correlogramList <-
169 function(x, lattice = TRUE, legend = TRUE,
170 test.level = 0.05, col = c("grey", "red"),
171 xlab = "", ylab = "Moran's I",
172 type = "b", pch = 21, cex = 2, ...)
175 obs <- unlist(lapply(x, "[[", "obs"))
176 pval <- unlist(lapply(x, "[[", "p.values"))
177 gr <- factor(unlist(lapply(x, "[[", "labels")),
178 ordered = TRUE, levels = x[[1]]$labels)
179 vars <- gl(n, nlevels(gr), labels = names(x))
180 BG <- col[(pval < test.level) + 1]
182 ## trellis.par.set(list(plot.symbol=list(pch=19)))
183 lattice::xyplot(obs ~ gr | vars, xlab = xlab, ylab = ylab,
184 panel = function(x, y) {
185 lattice::panel.lines(x, y, lty = 2)
186 lattice::panel.points(x, y, cex = cex, pch = 19, col = BG)
187 ##lattice::panel.abline(h = 0, lty = 3)
190 if (pch > 20 && pch < 26) {
192 CO <- rep("black", length(obs))
198 plot(as.numeric(gr), obs, type = "n", xlab = xlab,
199 ylab = ylab, xaxt = "n")
201 sel <- as.numeric(vars) == i
202 lines(as.numeric(gr[sel]), obs[sel], type = type, lty = i,
203 col = CO[sel], bg = BG[sel], pch = pch, cex = cex, ...)
205 axis(1, at = 1:length(x[[i]]$obs), labels = x[[i]]$labels)
207 legend("topright", legend = names(x), lty = 1:n, bty = "n")
208 legend("top", legend = paste(c("P >=", "P <"), test.level),
209 pch = pch, col = col, pt.bg = bg, pt.cex = cex, horiz = TRUE)