version 1.0-7

DESCRIPTION

@@ -2,8 +2,8 @@ Package: FD
 Type: Package
 Title: Measuring functional diversity (FD) from multiple traits, and
         other tools for functional ecology
-Version: 1.0-5
-Date: 2009-11-30
+Version: 1.0-7
+Date: 2010-02-12
 Author: Etienne Lalibert�, Bill Shipley
 Maintainer: Etienne Lalibert� <[email protected]>
 Description: FD is a package to compute different multidimensional FD
@@ -16,6 +16,6 @@ LazyLoad: yes
 LazyData: yes
 Depends: ade4, ape, geometry, vegan
 Encoding: latin1
-Packaged: 2009-11-29 22:25:06 UTC; etienne
+Packaged: 2010-03-09 07:44:22 UTC; etienne
 Repository: CRAN
-Date/Publication: 2009-11-30 14:27:23
+Date/Publication: 2010-03-09 07:56:20

R/maxent.R

@@ -1,59 +1,88 @@
-maxent <- function (c.means, c.mat, prior, tol = 1e-08, lambda = F){
+maxent <- function(constr, states, prior, tol = 1e-07, lambda = FALSE){
 
      # check input
-    if (!is.numeric(c.means)) stop("c.means must be a numeric vector\n")
-    if (!is.matrix(c.mat)){
-       if (is.numeric(c.mat)){
-          s.names <- names(c.mat)
-          c.mat <- matrix(c.mat, nrow = 1, ncol = length(c.mat) )
-          dimnames(c.mat) <- list("constraint", s.names)
+    if (is.vector(constr)){
+        means.names <- names(constr)
+        constr <- matrix(constr, 1, length(constr) ) ; dimnames(constr) <- list("set1", means.names)
+    }
+    if (is.data.frame(constr)) constr <- as.matrix(constr)
+    if (!is.numeric(constr)) stop("constr must only contain numeric values\n")
+    if (!is.numeric(tol)) stop("tol must be numeric\n")
+    if (!is.logical(lambda)) stop("lambda must be logical\n")
+    if (is.vector(states)){
+          s.names <- names(states)
+          states <- matrix(states, nrow = 1, ncol = length(states) )
+          dimnames(states) <- list("constraint", s.names)
          }
-       else stop("if c.mat is not a matrix then it can only be a numeric vector\n")
-     }
-    s.names <- dimnames(c.mat)[[2]]
-    c.names <- dimnames(c.mat)[[1]]
-    dim.matrix <- dim(c.mat)
-    if (dim.matrix[2] == 1 && dim.matrix[1] > 1){
-        c.mat <- t(c.mat)
-        dim.matrix <- dim(c.mat)
-     }
-    n.species <- dim.matrix[2]
-    if (missing(prior) ) prob <- rep(1 / n.species, n.species) else prob <- prior
-    if (length(prob) != n.species) stop("number of states in prior not equal to number in c.mat\n")
-    if (any(is.na(c.means)) || any(is.na(c.mat)) || any(is.na(prob)) ) stop("no NA's allowed\n")
-    n.constraints <- length(c.means)
-    if (n.constraints != dim.matrix[1]) stop("number of constraint means not equal to number of constraints in c.mat\n")
-
-    # run algorithm
-    C.values <- rowSums(c.mat)
-    test <- 1e+10
-    iter <- 0
-    while (test > tol){
-        iter <- iter + 1
-        denom <- c.mat %*% prob
-        gamma.values <- log(c.means / denom) / C.values
-        if (any(is.na(gamma.values) ) ) stop("NA's in gamma.values\n")
-        unstandardized <- exp(t(gamma.values) %*% c.mat) * prob
-        new.prob <- as.vector(unstandardized / sum(unstandardized) )
-        if (any(is.na(new.prob) ) ) stop("NA's in new.prob\n")
-        test <- max(abs(prob - new.prob))
-        prob <- new.prob
-     }
+    if (is.data.frame(states)) states <- as.matrix(states)
+    if (!is.numeric(states)) stop("states must only contain numeric values\n")
+    if (dim(states)[2] == 1 && dim(states)[1] > 1) states <- t(states)
+    s.names <- dimnames(states)[[2]]
+    c.names <- dimnames(states)[[1]]
+    set.names <- dimnames(constr)[[1]]
+    n.states <- dim(states)[2]
+    n.traits <- dim(states)[1]
+    n.sets <- dim(constr)[1]
+    if (n.traits != dim(constr)[2]) stop("number of constraints in constr should be equal to number of constraints in states\n")
+    if (missing(prior) ) {prior <- matrix(1 / n.states, n.sets, n.states)   ; dimnames(prior) <- list(set.names, s.names)}
+    if (is.vector(prior)){
+        if (length(prior) != n.states) {stop("number of states in prior should be equal to number in states\n")}
+        if (n.sets == 1) {prior <- matrix(prior, 1, length(prior) ) ; dimnames(prior) <- list("set1", s.names)}
+        else {prior <- matrix(rep(prior, n.sets), n.sets, length(prior), byrow = T ) ; dimnames(prior) <- list(set.names, s.names)}
+    }     
+    if (is.data.frame(prior)) prior <- as.matrix(prior)
+    if (!is.numeric(prior)) stop("prior must only contain numeric values\n")
+    if (dim(prior)[2] == 1 && dim(prior)[1] > 1) prior <- t(prior)
+    if (dim(prior)[2] != n.states) stop("number of states in prior should be equal to number in states\n")
+    if (dim(prior)[1] > 1 && dim(prior)[1] != n.sets) stop("number of rows in prior should be 1 or equal to number of rows in constr\n")
+    if (dim(prior)[1] == 1) prior <- matrix(rep(prior[1,], n.sets), n.sets, ncol(prior), byrow = T )
+    if (any(prior < 0 | prior > 1)) stop("prior must contain probabilities between 0 and 1\n")
+    prior <- t(apply(prior, 1, function(x) x / sum(x))); dimnames(prior) <- list(set.names, s.names)
+    if (any(is.na(constr)) || any(is.na(states)) || any(is.na(prior)) ) stop("no NA's allowed\n")
+
+    # create storage objects
+    allprobs <- matrix(NA, n.sets, n.states); dimnames(allprobs) <- list(set.names, s.names)
+    moments <- matrix(NA, n.sets, n.traits) ; dimnames(moments) <- list(set.names, c.names)
+    entropy <- rep(NA, n.sets) ; names(entropy) <- set.names
+    iter <- rep(NA, n.sets) ; names(iter) <- set.names
+    if (lambda){
+       lambdas <- matrix(NA, n.sets, n.traits + 1)
+       dimnames(lambdas) <- list(set.names, c("intercept", c.names) )
+    }
+    # FORTRAN loop
+    for (i in 1:n.sets){
+       itscale <- .Fortran("itscale5", as.double(t(states)), as.integer(n.states), as.integer(n.traits), as.double(constr[i,]), as.double(prior[i,]), prob = double(n.states), entropy = double(1), niter = integer(1), as.double(tol), moments = double(n.traits), PACKAGE = "FD")
+       allprobs[i, ] <- itscale$prob
+       moments[i, ] <- itscale$moments
+       entropy[i] <- itscale$entropy
+       iter[i] <- itscale$niter
+       if (lambda) lambdas[i, ] <- coef(lm(log(itscale$prob) ~ t(states) ) )
+   }
 
   # output
   res <- list()
-  names(prob) <- s.names
-  res$prob <- prob
-  moments <- as.numeric(denom)
-  names(moments) <- rownames(denom)
-  res$moments <- moments
-  res$entropy <- -1 * sum(prob * log(prob) )
-  res$iter <- iter
-  if (lambda){
-    lambda <- coef(lm(log(prob) ~ t(c.mat) ) )
-    names(lambda) <- c("intercept", c.names)
-    res$lambda <- lambda
-   }
+  if (n.sets == 1){
+     res$prob <- allprobs[1, ]
+     res$moments <- moments[1, ]
+     names(entropy) <- NULL
+     res$entropy <- entropy
+     names(iter) <- NULL
+     res$iter <- iter
+     if (lambda) res$lambda <- lambdas[1, ]
+     res$constr <- constr[1, ]
+     res$states <- states
+     res$prior <- prior[1, ]
+  }
+  else{
+    res$prob <- allprobs
+    res$moments <- moments
+    res$entropy <- entropy
+    res$iter <- iter
+    if (lambda) res$lambda <- lambdas
+    res$constr <- constr
+    res$states <- states
+    res$prior <- prior
+  }
   return(res)
 }
 

R/maxent.test.R

@@ -0,0 +1,200 @@
+maxent.test <- function(model, obs, sub.c, nperm = 99, quick = FALSE, alpha = 0.05, plot = TRUE){
+
+    # check input
+    if (is.vector(obs)){
+        s.names <- names(obs)
+        obs <- matrix(obs, 1, length(obs) ) ; dimnames(obs) <- list("set1", s.names)
+    }
+    if (is.data.frame(obs)) obs <- as.matrix(obs)
+    obs.names <- dimnames(obs)
+    if (!is.numeric(obs)) stop("obs must only contain numeric values\n")
+    if (dim(obs)[2] == 1 && dim(obs)[1] > 1) obs <- t(obs)
+    if (any(is.na(obs) )) stop("no NA's allowed\n")
+    if (any(obs < 0 | obs > 1)) stop("obs must contain probabilities between 0 and 1\n")
+    obs <- t(apply(obs, 1, function(x) x / sum(x))) ; dimnames(obs) <- obs.names
+    states <- model$states
+    s.names <- dimnames(states)[[2]]
+    constr <- model$constr
+    prob <- model$prob
+    prior <- model$prior
+    n.states <- dim(states)[2]
+    if (is.vector(constr)){
+       means.names <- names(constr)
+       constr <- matrix(constr, 1, length(constr) ) ; dimnames(constr) <- list("set1", means.names)
+       prior <- matrix(prior, 1, length(prior) ) ; dimnames(prior) <- list("set1", s.names)
+       prob <- matrix(prob, 1, length(prob) ) ; dimnames(prob) <- list("set1", s.names)
+    }
+    n.sets <- dim(constr)[1]
+    if (n.sets != dim(obs)[1]) stop("number of rows in obs and constr should be equal\n")
+
+    # function to compute statistic
+    stat <- function(o, p, q) sum(o * log(p / q))
+
+    # vector to store null values
+    values <- rep(NA, nperm)
+
+    # test1 - compute p values and stat
+    if (missing(sub.c)){
+       obs.stat <- stat(obs, prob, prior)
+       count <- 0
+       within.ci <- FALSE
+       while(!within.ci && count < nperm){
+          count <- count + 1
+          prob.temp <- matrix(NA, n.sets, n.states)
+             for (j in 1:n.sets){
+               shuffled <- sample(1:n.states, n.states)
+               states.perm <- states[, shuffled, drop = F]; colnames(states.perm) <- s.names
+               constr.perm <- functcomp(t(states.perm), obs[j, , drop = F])
+               prob.temp[j, ] <- maxent(constr.perm, states.perm, prior[j, ])$prob
+             }
+       values[count] <- stat(obs, prob.temp, prior)
+       if (quick){
+          val.temp <- values[1:count]
+          p.temp <- (length(val.temp[val.temp >= obs.stat]) + 1) / (length(val.temp) + 1)
+          ci.hi <- p.temp + 1.96 * sqrt( (p.temp * (1 - p.temp)) / count )
+          within.ci <- ci.hi <= alpha
+          if (within.ci) nperm = count
+        }
+      }
+    }
+
+    # test2 - compute p values and stat
+    else{
+       if (length(sub.c) >= n.states) stop("sub.c contains as many or more elements than there are states\n") 
+       if (is.character(sub.c) && !all(sub.c %in% dimnames(states)[[1]]) ) stop("sub.c does not match the names of the state attributes\n")
+       if (is.character(sub.c) && !all(sub.c %in% dimnames(constr)[[2]]) ) stop("sub.c does not match the constraint names\n")
+       if (is.character(sub.c)) sub.c <- which(dimnames(states)[[1]] %in% sub.c)
+       if (!is.vector(sub.c) ) stop("sub.c must be a vector\n")
+       prob.a <- maxent(constr[, -sub.c, drop = F], states[-sub.c, , drop = F], prior)$prob
+       obs.stat <- stat(obs, prob, prob.a)
+       count <- 0
+       within.ci <- FALSE
+       while(!within.ci && count < nperm){
+          count <- count + 1
+          prob.temp <- matrix(NA, n.sets, n.states)
+             for (j in 1:n.sets){
+               shuffled <- sample(1:n.states, n.states)
+               states.perm <- states
+               states.perm[sub.c, ] <- states[sub.c, shuffled, drop = F]; colnames(states.perm) <- s.names
+               constr.perm <- functcomp(t(states.perm), obs[j, , drop = F])
+               prob.temp[j, ] <- maxent(constr.perm, states.perm, prior[j, ])$prob
+             }
+       values[count] <- stat(obs, prob.temp, prob.a)
+       if (quick){
+          val.temp <- values[1:count]
+          p.temp <- (length(val.temp[val.temp >= obs.stat]) + 1) / (length(val.temp) + 1)
+          ci.hi <- p.temp + 1.96 * sqrt( (p.temp * (1 - p.temp)) / count )
+          within.ci <- ci.hi <= alpha
+          if (within.ci) nperm = count
+        }
+      }
+    }
+
+     # permutational p value and ci
+     values <- values[!is.na(values)]
+     p.perm <- (length(values[values >= obs.stat]) + 1) / (length(values) + 1)
+     p.perm.hi <- p.perm + 1.96 * sqrt( (p.perm * (1 - p.perm)) / nperm )
+     p.perm.lo <- p.perm - 1.96 * sqrt( (p.perm * (1 - p.perm)) / nperm )
+
+    # measure of fit beween obs and pred relative to prior
+    opqfit <- function(o, p, q) 1 - (sum((o - p)^2) / sum((o - q)^2) )
+    r2.op <- cor(as.double(prob), as.double(obs) )^2
+    if (length(unique(as.double(prior))) != 1 )   r2.oq <- cor(as.double(prior), as.double(obs) )^2
+    else r2.oq = 0
+    fit <- opqfit(obs, prob, prior)
+    if (!missing(sub.c) ){
+         r2.opa <- cor(as.double(prob.a), as.double(obs) )^2
+         fit.a <- opqfit(obs, prob.a, prior)
+    }
+    # plot results
+    if (plot){
+        if (missing(sub.c) ){
+             par(las = 1, mfcol = c(2, 2), oma = c(0, 0, 3, 0))
+             fit.text <- bquote(fit[bold(o*","*p)*"|"*bold(q)] ==.(round(fit, 3)))
+             r2.op.text <- bquote(italic(r)^2 ==.(round(r2.op, 3)))
+             r2.oq.text <- bquote(italic(r)^2 ==.(round(r2.oq, 3)))
+             plot(as.double(obs), as.double(prob), xlim = c(0, 1), ylim = c(0, 1), xlab = "", ylab = "predicted probabilities", main = expression(bold(p)(bold(C)*","*~bold(q)) ) )
+             abline(0, 1, col = "grey25")
+             text(0.1, 0.9, fit.text, cex = 1.2, pos = 4) 
+             text(0.1, 0.75, r2.op.text, cex = 1.2, pos = 4) 
+             lines(x = c(0.05, 0.05), c(0, 1), lty = "dashed", col = "grey50")
+             lines(x = c(0, 1), c(0.05, 0.05), lty = "dashed", col = "grey50")
+             mtext("arithmetic scale", line = 0)
+             plot(as.double(obs) + 1e-4, as.double(prob) + 1e-4, xlim = c(1e-4, 1), ylim = c(1e-4, 1), xlab = "observed probabilities", ylab = "predicted probabilities", log = "xy", main = expression(bold(p)(bold(C)*","*~bold(q)) ))
+             abline(0, 1, col = "grey25")
+             lines(x = c(0.05 + 1e-4, 0.05 + 1e-4), c(0 + 1e-4, 1 + 1e-4), lty = "dashed", col = "grey50")
+             lines(x = c(1e-4, 1 + 1e-4), c(0.05+1e-4, 0.05 + 1e-4) , lty = "dashed", col = "grey50")
+             mtext("log10 scale, + 1e-4", line = 0)
+
+             plot(as.double(obs), as.double(prior), xlim = c(0, 1), ylim = c(0, 1), xlab = "", ylab = "prior", main = expression(bold(q) ) )
+             abline(0, 1, col = "grey25")
+             text(0.1, 0.9, r2.oq.text, cex = 1.2, pos = 4)
+             lines(x = c(0.05, 0.05), c(0, 1), lty = "dashed", col = "grey50")
+             lines(x = c(0, 1), c(0.05, 0.05), lty = "dashed", col = "grey50")
+             mtext("arithmetic scale", line = 0)
+             plot(as.double(obs) + 1e-4, as.double(prior) + 1e-4, xlim = c(1e-4, 1), ylim = c(1e-4, 1), xlab = "observed probabilities", ylab = "prior", log = "xy", main = expression(bold(q) ))
+             abline(0, 1, col = "grey25")
+             lines(x = c(0.05 + 1e-4, 0.05 + 1e-4), c(0 + 1e-4, 1 + 1e-4), lty = "dashed", col = "grey50")
+             lines(x = c(1e-4, 1 + 1e-4), c(0.05+1e-4, 0.05 + 1e-4) , lty = "dashed", col = "grey50")
+             mtext("log10 scale, + 1e-4", line = 0)
+             mtext(expression(H[0] %->% bold(p)(bold(C)*","*~bold(q))==bold(q) ),  cex = 1.2, outer = T, las = 1, line = 1)
+             p.text <- bquote(italic(P) ==.(round(p.perm, 3)))
+             mtext(p.text,  cex = 1.2, outer = T, las = 1, line = -0.5)
+         }
+
+
+
+        else{
+             par(las = 1, mfcol = c(2, 2), oma = c(0, 0, 3, 0))
+             fit.text <- bquote(fit[bold(o*","*p)*"|"*bold(q)] ==.(round(fit, 3)))
+             fit.a.text <- bquote(fit[bold(o*","*p)*"|"*bold(q)] ==.(round(fit.a, 3)))
+             r2.op.text <- bquote(italic(r)^2 ==.(round(r2.op, 3)))
+             r2.opa.text <- bquote(italic(r)^2 ==.(round(r2.opa, 3)))
+             plot(as.double(obs), as.double(prob), xlim = c(0, 1), ylim = c(0, 1), xlab = "", ylab = "predicted probabilities", main = expression(bold(p)(bold(A)~union(bold(B)*","*~bold(q))) ) )
+             abline(0, 1, col = "grey25")
+             text(0.1, 0.9, fit.text, cex = 1.2, pos = 4) 
+             text(0.1, 0.75, r2.op.text, cex = 1.2, pos = 4) 
+             lines(x = c(0.05, 0.05), c(0, 1), lty = "dashed", col = "grey50")
+             lines(x = c(0, 1), c(0.05, 0.05), lty = "dashed", col = "grey50")
+             mtext("arithmetic scale", line = 0)
+             plot(as.double(obs) + 1e-4, as.double(prob) + 1e-4, xlim = c(1e-4, 1), ylim = c(1e-4, 1), xlab = "observed probabilities", ylab = "predicted probabilities", log = "xy",  main = expression(bold(p)(bold(A)~union(bold(B)*","*~bold(q))) )  )
+             abline(0, 1, col = "grey25")
+             lines(x = c(0.05 + 1e-4, 0.05 + 1e-4), c(0 + 1e-4, 1 + 1e-4), lty = "dashed", col = "grey50")
+             lines(x = c(1e-4, 1 + 1e-4), c(0.05+1e-4, 0.05 + 1e-4) , lty = "dashed", col = "grey50")
+             mtext("log10 scale, + 1e-4", line = 0)
+
+             plot(as.double(obs), as.double(prob.a), xlim = c(0, 1), ylim = c(0, 1), xlab = "", ylab = "", main = expression(bold(p)(bold(A)*","*~bold(q)) ) )
+             abline(0, 1, col = "grey25")
+             text(0.1, 0.9, fit.a.text, cex = 1.2, pos = 4) 
+             text(0.1, 0.75, r2.opa.text, cex = 1.2, pos = 4) 
+             lines(x = c(0.05, 0.05), c(0, 1), lty = "dashed", col = "grey50")
+             lines(x = c(0, 1), c(0.05, 0.05), lty = "dashed", col = "grey50")
+             mtext("arithmetic scale", line = 0)
+             plot(as.double(obs) + 1e-4, as.double(prob.a) + 1e-4, xlim = c(1e-4, 1), ylim = c(1e-4, 1), xlab = "observed probabilities", ylab = "", log = "xy", main = expression(bold(p)(bold(A)*","*~bold(q)) ) )
+             abline(0, 1, col = "grey25")
+             lines(x = c(0.05 + 1e-4, 0.05 + 1e-4), c(0 + 1e-4, 1 + 1e-4), lty = "dashed", col = "grey50")
+             lines(x = c(1e-4, 1 + 1e-4), c(0.05+1e-4, 0.05 + 1e-4) , lty = "dashed", col = "grey50")
+             mtext("log10 scale, + 1e-4", line = 0)
+             mtext(expression(H[0] %->% bold(p)(bold(A)~union(bold(B))*","*~bold(q)) == bold(p)(bold(A)*","*~bold(q)) ),  cex = 1.2, outer = T, las = 1, line = 1)
+             p.text <- bquote(italic(P) ==.(round(p.perm, 3)))
+             mtext(p.text,  cex = 1.2, outer = T, las = 1, line = -0.5)
+        }
+
+    }
+
+    # return results
+    res <- list()
+    res$fit <- fit
+    if (!missing(sub.c) ){
+       res$fit.a <- fit.a
+       res$r2.a <- r2.opa 
+    }
+    res$r2 <- r2.op
+    if(missing(sub.c)) res$r2.q <- r2.oq
+    res$obs.stat <- obs.stat
+    res$nperm <- nperm
+    res$pval <- p.perm
+    res$ci.pval <- c(p.perm.lo, p.perm.hi)
+    return(res)
+}
+

inst/CITATION

@@ -4,9 +4,11 @@ citEntry(entry="article",
          title = "A distance-based framework for measuring functional diversity from multiple traits",
          author = personList(person("Etienne", "Lalibert�", email="[email protected]"), person("Pierre", "Legendre", email="[email protected]")),
  journal="Ecology",
- year="in press",
+ year="2010",
+ volume="91",
+ pages="299--305",
 
- textVersion="Lalibert�, E., and P. Legendre (in press) A distance-based framework for measuring functional diversity from multiple traits. Ecology.")
+ textVersion="Lalibert�, E., and P. Legendre (2010) A distance-based framework for measuring functional diversity from multiple traits. Ecology 91:299-305.")
 
 if(!exists("meta") || is.null(meta)) meta <- packageDescription("foo")
 year <- sub(".*(2[[:digit:]]{3})-.*", "\\1", meta$Date)
@@ -24,4 +26,4 @@ year,
 vers, ".", sep=""))
 
 
-citFooter("See http://www.elaliberte.info/publications for the final reference when the article is published.")
+citFooter("See http://www.elaliberte.info/publications for a PDF version of the article.")

inst/NEWS

@@ -3,6 +3,17 @@
 FD RELEASE VERSIONS
 ===================
 
+CHANGES IN FD 1.0-7
+
+	- added maxent.test, following Shipley (2010)
+        - changed tol to 1e-7 in maxent
+        - modified arguments of maxent
+        - maxent can now be run on multiple sites
+
+CHANGES IN FD 1.0-6
+
+	- added compiled FORTRAN code for loop in maxent(). About 10x boost in speed.
+
 CHANGES IN FD 1.0-5
 
 	- corrector minor error in functcomp()