R/multispati.R

In ade4: Analysis of Ecological Data: Exploratory and Euclidean Methods in Environmental Sciences

"multispati" <- function(dudi, listw, scannf=TRUE, nfposi=2, nfnega=0) {
  
  .Deprecated(new="multispati", package="ade4", 
              msg="This function is now deprecated. Please use the 'multispati' function in the 'adespatial' package.")
  
  if(!inherits(dudi,"dudi")) stop ("object of class 'dudi' expected")
  if(!inherits(listw,"listw")) stop ("object of class 'listw' expected") 
  if(listw$style!="W") stop ("object of class 'listw' with style 'W' expected")
  NEARZERO <- 1e-14
  
  dudi$cw <- dudi$cw
  fun <- function (x) spdep::lag.listw(listw,x,TRUE)
  tablag <- apply(dudi$tab,2,fun)
  covar <- t(tablag)%*%as.matrix((dudi$tab*dudi$lw))
  covar <- (covar+t(covar))/2
  covar <- covar * sqrt(dudi$cw)
  covar <- t(t(covar) * sqrt(dudi$cw))
  covar <- eigen(covar, symmetric = TRUE)
  res <- list()
  res$eig <- covar$values[abs(covar$values)>NEARZERO]
  ndim <- length(res$eig)
  covar$vectors <- covar$vectors[, abs(covar$values)>NEARZERO]
  
  if (scannf) {
    barplot(res$eig)
    cat("Select the first number of axes (>=1): ")
    nfposi <- as.integer(readLines(n = 1))
    
    cat("Select the second number of axes (>=0): ")
    nfnega <- as.integer(readLines(n = 1))
  }
  if (nfposi <= 0)  nfposi <- 1
  if (nfnega<=0) nfnega <- 0       
  
  if(nfposi > sum(res$eig > 0)){
    nfposi <- sum(res$eig > 0)
    warning(paste("There are only",sum(res$eig>0),"positive factors."))
  }
  if(nfnega > sum(res$eig < 0)){
    nfnega <- sum(res$eig < 0)
    warning(paste("There are only",sum(res$eig< 0),"negative factors."))
  }
  res$nfposi <- nfposi
  res$nfnega <- nfnega
  agarder <- c(1:nfposi,if (nfnega>0) (ndim-nfnega+1):ndim else NULL)
  dudi$cw[which(dudi$cw == 0)] <- 1
  auxi <- data.frame(covar$vectors[, agarder] /sqrt(dudi$cw))
  names(auxi) <- paste("CS", agarder, sep = "")
  row.names(auxi) <- names(dudi$tab)
  res$c1 <- auxi                     
  auxi <- as.matrix(auxi)*dudi$cw
  auxi1 <- as.matrix(dudi$tab)%*%auxi
  auxi1 <- data.frame(auxi1)
  names(auxi1) <- names(res$c1)
  row.names(auxi1) <- row.names(dudi$tab)
  res$li <- auxi1
  auxi1 <- as.matrix(tablag)%*%auxi
  auxi1 <- data.frame(auxi1)
  names(auxi1) <- names(res$c1)
  row.names(auxi1) <-  row.names(dudi$tab)    
  res$ls <- auxi1
  auxi <- as.matrix(res$c1) * unlist(dudi$cw)
  auxi <- data.frame(t(as.matrix(dudi$c1)) %*% auxi)
  row.names(auxi) <- names(dudi$li)
  names(auxi) <- names(res$li)
  res$as <- auxi
  res$call <- match.call()
  class(res) <- "multispati"
  return(res)
}


"summary.multispati" <- function (object, ...) {
  
  .Deprecated(new="summary.multispati", package="ade4", 
              msg="This method is now deprecated. Please use the 'summary.multispati' method in the 'adespatial' package.")
  
  norm.w <- function(X, w) {
    f2 <- function(v) sum(v * v * w)/sum(w)
    norm <- apply(X, 2, f2)
    return(norm)
  }
  
  if (!inherits(object, "multispati")) stop("to be used with 'multispati' object")
  
  cat("\nMultivariate Spatial Analysis\n")
  cat("Call: ")
  print(object$call)
  
  appel <- as.list(object$call)
  dudi <- eval.parent(appel$dudi)
  listw <- eval.parent(appel$listw)
  
  ## les scores de l'analyse de base
  nf <- dudi$nf
  eig <- dudi$eig[1:nf]
  cum <- cumsum (dudi$eig) [1:nf]
  ratio <- cum/sum(dudi$eig)
  w <- apply(dudi$l1,2,spdep::lag.listw,x=listw)
  moran <- apply(w*as.matrix(dudi$l1)*dudi$lw,2,sum)
  res <- data.frame(var=eig,cum=cum,ratio=ratio, moran=moran)
  cat("\nScores from the initial duality diagramm:\n")
  print(res)
  
  ## les scores de l'analyse spatiale
  ## on recalcule l'objet en gardant tous les axes
  eig <- object$eig
  nfposi <- object$nfposi
  nfnega <- object$nfnega
  nfposimax <- sum(eig > 0)
  nfnegamax <- sum(eig < 0)
  
  ms <- multispati(dudi=dudi, listw=listw, scannf=FALSE,
                   nfposi=nfposimax, nfnega=nfnegamax)
  
  ndim <- dudi$rank
  nf <- nfposi + nfnega
  agarder <- c(1:nfposi,if (nfnega>0) (ndim-nfnega+1):ndim else NULL)
  varspa <- norm.w(ms$li,dudi$lw)
  moran <- apply(as.matrix(ms$li)*as.matrix(ms$ls)*dudi$lw,2,sum)
  res <- data.frame(eig=eig,var=varspa,moran=moran/varspa)
  
  cat("\nMultispati eigenvalues decomposition:\n")
  print(res[agarder,])
  return(invisible(res))
}



"print.multispati" <- function(x, ...) {
  
  .Deprecated(new="print.multispati", package="ade4", 
              msg="This method is now deprecated. Please use the 'print.multispati' method in the 'adespatial' package.")
  
  cat("Multispati object \n")
  cat("class: ")
  cat(class(x))
  cat("\n$call: ")
  print(x$call)
  cat("\n$nfposi:", x$nfposi, "axis-components saved")
  cat("\n$nfnega:", x$nfnega, "axis-components saved")
  #cat("\n$rank: ")
  #cat(x$rank)
  cat("\nPositive eigenvalues: ")
  l0 <- sum(x$eig >= 0)
  cat(signif(x$eig, 4)[1:(min(5, l0))])
  if (l0 > 5) 
    cat(" ...\n")
  else cat("\n")  
  cat("Negative eigenvalues: ")
  l0 <- sum(x$eig <= 0)
  cat(sort(signif(x$eig, 4))[1:(min(5, l0))])
  if (l0 > 5) 
    cat(" ...\n")
  else cat("\n")
  cat('\n')
  sumry <- array("", c(1, 4), list(1, c("vector", "length", 
                                        "mode", "content")))
  sumry[1, ] <- c('$eig', length(x$eig), mode(x$eig), 'eigen values')
  
  print(sumry, quote = FALSE)
  cat("\n")
  sumry <- array("", c(4, 4), list(1:4, c("data.frame", "nrow", "ncol", "content")))
  sumry[1, ] <- c("$c1", nrow(x$c1), ncol(x$c1), "column normed scores")
  sumry[2, ] <- c("$li", nrow(x$li), ncol(x$li), "row coordinates")
  sumry[3, ] <- c("$ls", nrow(x$ls), ncol(x$ls), 'lag vector coordinates')
  sumry[4, ] <- c("$as", nrow(x$as), ncol(x$as), 'inertia axes onto multispati axes')
  
  
  print(sumry, quote = FALSE)
  cat("other elements: ")
  if (length(names(x)) > 8) 
    cat(names(x)[9:(length(names(x)))], "\n")
  else cat("NULL\n")
}



"plot.multispati" <- function (x, xax = 1, yax = 2, ...) {
  
  .Deprecated(new="plot.multispati", package="ade4", 
              msg="This method is now deprecated. Please use the 'plot.multispati' method in the 'adespatial' package.")
  
  if (!inherits(x, "multispati")) 
    stop("Use only with 'multispati' objects")
  
  appel <- as.list(x$call)
  dudi <- eval.parent(appel$dudi)
  nf <- x$nfposi + x$nfnega
  if ((nf == 1) || (xax == yax)) {
    sco.quant(x$li[, 1], dudi$tab)
    return(invisible())
  }
  if (xax > nf) 
    stop("Non convenient xax")
  if (yax > nf) 
    stop("Non convenient yax")
  f1 <- function () 
  {
    opar <- par(mar = par("mar"))
    on.exit(par(opar))
    m <- length(x$eig)
    par(mar = c(0.8, 2.8, 0.8, 0.8))
    col.w <- rep(grey(1), m) # elles sont toutes blanches
    col.w[1:x$nfposi] <- grey(0.8)
    if (x$nfnega>0) col.w[m:(m-x$nfnega+1)] = grey(0.8)
    j1 <- xax
    if (j1>x$nfposi) j1 = j1-x$nfposi +m -x$nfnega
    j2 <- yax
    if (j2>x$nfposi) j2 = j2-x$nfposi +m -x$nfnega
    col.w[c(j1,j2)] = grey(0)
    barplot(x$eig, col = col.w)
    scatterutil.sub(cha ="Eigen values", csub = 2, possub = "topright")
  }
  
  def.par <- par(no.readonly = TRUE)
  on.exit(par(def.par))
  layout(matrix(c(3, 3, 1, 3, 3, 2), 3, 2))
  par(mar = c(0.2, 0.2, 0.2, 0.2))
  f1()
  s.arrow(x$c1, xax = xax, yax = yax, sub = "Canonical weights", 
          csub = 2, clabel = 1.25)
  s.match(x$li, x$ls, xax = xax, yax = yax, sub = "Scores and lag scores", csub = 2, clabel = 0.75) 
  
}