R/MoranI.R

   1 ## MoranI.R (2008-01-14)
   2
   3 ##   Moran's I Autocorrelation Index
   4
   5 ## Copyright 2004 Julien Dutheil, 2007-2008 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 ## code cleaned-up by EP (Dec. 2007)
  11
  12 Moran.I <- function(x, weight, scaled = FALSE, na.rm = FALSE,
  13                     alternative = "two.sided")
  14 {
  15     if(dim(weight)[1] != dim(weight)[2])
  16         stop("'weight' must be a square matrix")
  17     n <- length(x)
  18     if(dim(weight)[1] != n)
  19         stop("'weight' must have as many rows as observations in 'x'")
  20     ## Expected mean:
  21     ei <- -1/(n - 1)
  22
  23     nas <- is.na(x)
  24     if (any(nas)) {
  25         if (na.rm) {
  26             x <- x[!nas]
  27             n <- length(x)
  28             weight <- weight[!nas, !nas]
  29         } else {
  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))
  32         }
  33     }
  34
  35     ## normaling the weights:
  36     ## Note that we normalize after possibly removing the
  37     ## missing data.
  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
  42
  43     s <- sum(weight)
  44     m <- mean(x)
  45     y <- x - m # centre the x's
  46     cv <- sum(weight * y %o% y)
  47     v <- sum(y^2)
  48     obs <- (n/s) * (cv/v)
  49     ## Scaling:
  50     if (scaled) {
  51         i.max <- (n/s) * (sd(rowSums(weight) * y)/sqrt(v/(n - 1)))
  52         obs <- obs/i.max
  53     }
  54     ## Expected sd:
  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:
  58     ##S2 <- 0
  59     ##for (i in 1:n)
  60     ##    S2 <- S2 + (sum(weight[i, ]) + sum(weight[, i]))^2
  61
  62     s.sq <- s^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))
  67
  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)
  74 }
  75
  76 weight.taxo <- function(x)
  77 {
  78     d <- outer(x, x, "==")
  79     diag(d) <- 0 # implicitly converts 'd' into numeric
  80     d
  81 }
  82
  83 weight.taxo2 <- function(x, y)
  84 {
  85     d <- outer(x, x, "==") & outer(y, y, "!=")
  86     diag(d) <- 0
  87     d
  88 }
  89
  90 correlogram.formula <- function(formula, data = NULL, use = "all.obs")
  91 {
  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)
  95
  96     lhs <- formula[[2]]
  97     y.nms <- if (length(lhs) > 1)
  98         unlist(strsplit(as.character(as.expression(lhs)), " \\+ "))
  99     else as.character(as.expression(lhs))
 100
 101     rhs <- formula[[3]]
 102     gr.nms <- if (length(rhs) > 1)
 103         rev(unlist(strsplit(as.character(as.expression(rhs)), "/")))
 104     else as.character(as.expression(rhs))
 105
 106     if (is.null(data)) {
 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))
 110     } else {
 111         y <- data[y.nms]
 112         gr <- data[gr.nms]
 113     }
 114     if (use == "all.obs") {
 115         na.fail(y)
 116         na.fail(gr)
 117     }
 118     if (use == "complete.obs") {
 119         sel <- complete.cases(y, gr)
 120         y <- y[sel]
 121         gr <- gr[sel]
 122     }
 123     na.rm <- use == "pairwise.complete.obs"
 124
 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
 129             xx <- x[sel]
 130             g <- gr[sel, i]
 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
 135         }
 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"))
 140     }
 141
 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")
 145 }
 146
 147 plot.correlogram <-
 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, ...)
 151 {
 152     BG <- col[(x$p.values < test.level) + 1]
 153     if (pch > 20 && pch < 26) {
 154         bg <- col
 155         col <- CO <- "black"
 156     } else {
 157         CO <- BG
 158         BG <- bg <- NULL
 159     }
 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)
 163     if (legend)
 164         legend("top", legend = paste(c("P >=", "P <"), test.level),
 165                pch = pch, col = col, pt.bg = bg, pt.cex = cex, horiz = TRUE)
 166 }
 167
 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, ...)
 173 {
 174     n <- length(x)
 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]
 181     if (lattice) {
 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)
 188                         })
 189     } else {
 190         if (pch > 20 && pch < 26) {
 191             bg <- col
 192             CO <- rep("black", length(obs))
 193             col <- "black"
 194         } else {
 195             CO <- BG
 196             BG <- bg <- NULL
 197         }
 198         plot(as.numeric(gr), obs, type = "n", xlab = xlab,
 199              ylab = ylab, xaxt = "n")
 200         for (i in 1: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, ...)
 204         }
 205         axis(1, at = 1:length(x[[i]]$obs), labels = x[[i]]$labels)
 206         if (legend) {
 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)
 210         }
 211     }
 212 }