[Vegan-commits] r521 - pkg/R

[noreply at r-forge.r-project.org]

Author: jarioksa
Date: 2008-10-06 13:04:20 +0200 (Mon, 06 Oct 2008)
New Revision: 521

Modified:
   pkg/R/rad.lognormal.R
   pkg/R/rad.null.R
   pkg/R/rad.preempt.R
   pkg/R/rad.zipf.R
   pkg/R/rad.zipfbrot.R
Log:
radfit: improved and more consisntent handling of communities with 0, 1, or 2 species

Modified: pkg/R/rad.lognormal.R
===================================================================
--- pkg/R/rad.lognormal.R	2008-10-05 17:59:22 UTC (rev 520)
+++ pkg/R/rad.lognormal.R	2008-10-06 11:04:20 UTC (rev 521)
@@ -6,10 +6,10 @@
     rnk <- -qnorm(ppoints(n))
     fam <- family(link = "log")
     ## Must be > 2 species to fit a model
-    if (length(x) > 2)
+    if (length(x) > 1)
         ln <- try(glm(x ~ rnk, family = fam))
-    if (length(x) <= 2) {
-        aic <- 4
+    if (length(x) < 2) {
+        aic <- NA
         dev <- rdf <-  0
         ln <- nl <- NA
         p <- rep(NA, 2)

Modified: pkg/R/rad.null.R
===================================================================
--- pkg/R/rad.null.R	2008-10-05 17:59:22 UTC (rev 520)
+++ pkg/R/rad.null.R	2008-10-06 11:04:20 UTC (rev 521)
@@ -7,14 +7,17 @@
     x <- as.rad(x)
     nsp <- length(x)
     wt <- rep(1, nsp)
-    if (nsp > 1) 
+    if (nsp > 0) { 
         fit <- rev(cumsum(1/nsp:1)/nsp) * sum(x)
-    else
-        fit <- x
+        aic <- aicfun(x, nsp, fit, wt, deviance)
+    }
+    else {
+        fit <- NA
+        aic <- NA
+    }
     res <- dev.resids(x, fit, wt)
     deviance <- sum(res)
     residuals <- x - fit
-    aic <- aicfun(x, wt, fit, wt, deviance)
     rdf <- nsp
     p <- NA
     names(p) <- "S"

Modified: pkg/R/rad.preempt.R
===================================================================
--- pkg/R/rad.preempt.R	2008-10-05 17:59:22 UTC (rev 520)
+++ pkg/R/rad.preempt.R	2008-10-06 11:04:20 UTC (rev 521)
@@ -17,6 +17,7 @@
     linkinv <- fam$linkinv
     dev.resids <- fam$dev.resids
     x <- as.rad(x)
+    nsp <- length(x)
     rnk <- seq(along = x) - 1
     wt <- rep(1, length(x))
     logJ <- log(sum(x))
@@ -28,17 +29,20 @@
         p <- rep(NA, 1)
         fit <- residuals <- prior.weights <- rep(NA, length(x))
     } else {
-        if (length(x) > 1) {
+        if (nsp > 1) {
             p <- plogis(canon$estimate)
             fit <- exp(logJ + log(p) + log(1 - p) * rnk)
         } else {
-            p <- 1
+            p <- if (nsp > 0) 1 else NA
             fit <- x
         }
         res <- dev.resids(x, fit, wt)
         deviance <- sum(res)
         residuals <- x - fit
-        aic <- aicfun(x, rep(1, length(x)), fit, wt, deviance) + 2
+        if (nsp > 0)
+            aic <- aicfun(x, rep(1, length(x)), fit, wt, deviance) + 2
+        else
+            aic <- NA
         rdf <- length(x) - 1
     }
     names(p) <- c("alpha")

Modified: pkg/R/rad.zipf.R
===================================================================
--- pkg/R/rad.zipf.R	2008-10-05 17:59:22 UTC (rev 520)
+++ pkg/R/rad.zipf.R	2008-10-06 11:04:20 UTC (rev 521)
@@ -5,10 +5,10 @@
     rnk <- seq(along = x)
     off <- rep(log(sum(x)), length(x))
     fam <- family(link = "log")
-    if (length(x) > 2)
+    if (length(x) > 1)
         ln <- try(glm(x ~ log(rnk) + offset(off), family = fam))
-    if (length(x) <= 2) {
-        aic <- 4
+    if (length(x) < 2) {
+        aic <- NA
         dev <- rdf <-  0
         ln <- nl <- NA
         p <- rep(NA, 2)

Modified: pkg/R/rad.zipfbrot.R
===================================================================
--- pkg/R/rad.zipfbrot.R	2008-10-05 17:59:22 UTC (rev 520)
+++ pkg/R/rad.zipfbrot.R	2008-10-06 11:04:20 UTC (rev 521)
@@ -11,11 +11,11 @@
     off <- rep(log(sum(x)), length(x))
     p <- 0
     fam <- family(link = "log")
-    if (length(x) > 3) 
+    if (length(x) > 2) 
         nl <- try(nlm(mandelfun, p = p, x = x, rnk = rnk, off = off, 
                       family = fam, hessian = TRUE, ...))
-    if (length(x) <= 3) {
-        aic <- 6
+    if (length(x) < 3) {
+        aic <- NA
         dev <- rdf <-  0
         ln <- nl <- NA
         p <- rep(NA, 3)