[Vegan-commits] r516 - in pkg: R inst
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 2 16:06:00 CEST 2008
Author: jarioksa
Date: 2008-10-02 16:05:59 +0200 (Thu, 02 Oct 2008)
New Revision: 516
Modified:
pkg/R/plot.radfit.R
pkg/R/plot.radfit.frame.R
pkg/R/print.radfit.R
pkg/R/print.radfit.frame.R
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
pkg/R/radfit.data.frame.R
pkg/R/radfit.default.R
pkg/inst/ChangeLog
Log:
radfit: handling of rows with no species less or equal to the number or parameters
Modified: pkg/R/plot.radfit.R
===================================================================
--- pkg/R/plot.radfit.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/plot.radfit.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -1,7 +1,11 @@
"plot.radfit" <-
function (x, BIC = FALSE, legend = TRUE, ...)
{
+ if (length(x$y) == 0)
+ stop("No species, nothing to plot")
out <- plot(x$y, ...)
+ if (length(x$y) == 1)
+ return(invisible(out))
fv <- fitted(x)
if (BIC)
k = log(length(x$y))
Modified: pkg/R/plot.radfit.frame.R
===================================================================
--- pkg/R/plot.radfit.frame.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/plot.radfit.frame.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -19,7 +19,7 @@
}
Nhm <- length(x)
Abundance <- unlist(lapply(x, function(x) x$y))
- Rank <- unlist(lapply(x, function(x) 1:length(x$y)))
+ Rank <- unlist(lapply(x, function(x) if (length(x$y) > 0) 1:length(x$y) else NULL))
Site <- unlist(lapply(x, function(x) length(x$y)))
N <- Site
sitenames <- names(Site)
@@ -28,8 +28,9 @@
order.by <- 1:Nhm
else order.by <- order(order.by)
Site <- factor(Site, levels = sitenames[order.by])
- fit <- unlist(lapply(x, function(x) fitted(x)[, pickmod(x,
- pick, BIC)]))
+ fit <- unlist(lapply(x, function(x)
+ as.matrix(fitted(x))[, pickmod(x,
+ pick, BIC)]))
take <- sapply(x, function(x) pickmod(x, pick, BIC))
take <- rep(take, N)
cols <- trellis.par.get("superpose.line")$col
@@ -59,4 +60,3 @@
}, ...)
out
}
-
Modified: pkg/R/print.radfit.R
===================================================================
--- pkg/R/print.radfit.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/print.radfit.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -5,7 +5,8 @@
cat("No. of species ", length(x$y), ", total abundance ",
sum(x$y), "\n\n", sep = "")
p <- coef(x)
- p <- formatC(p, format="g", flag = " ", digits = digits)
+ if (any(!is.na(p)))
+ p <- formatC(p, format="g", flag = " ", digits = digits)
p <- apply(p, 2, function(x) gsub("NA", " ", x))
aic <- sapply(x$models, AIC)
bic <- sapply(x$models, AIC, k = log(length(x$y)))
@@ -15,4 +16,3 @@
print(out, quote=FALSE)
invisible(x)
}
-
Modified: pkg/R/print.radfit.frame.R
===================================================================
--- pkg/R/print.radfit.frame.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/print.radfit.frame.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -3,7 +3,6 @@
{
cat("\nDeviance for RAD models:\n\n")
out <- sapply(x, function(x) unlist(lapply(x$models, deviance)))
- colnames(out) <- colnames(out, do.NULL = FALSE, prefix = "row")
printCoefmat(out, na.print = "", ...)
invisible(x)
}
Modified: pkg/R/rad.lognormal.R
===================================================================
--- pkg/R/rad.lognormal.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/rad.lognormal.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -5,8 +5,19 @@
n <- length(x)
rnk <- -qnorm(ppoints(n))
fam <- family(link = "log")
- ln <- try(glm(x ~ rnk, family = fam))
- if (inherits(ln, "try-error")) {
+ ## Must be > 2 species to fit a model
+ if (length(x) > 2)
+ ln <- try(glm(x ~ rnk, family = fam))
+ if (length(x) <= 2) {
+ aic <- 4
+ dev <- rdf <- 0
+ ln <- nl <- NA
+ p <- rep(NA, 2)
+ fit <- x
+ res <- rep(0, length(x))
+ wts <- rep(1, length(x))
+ }
+ else if (inherits(ln, "try-error")) {
aic <- rdf <- ln <- nl <- dev <- NA
p <- rep(NA, 2)
fit <- res <- wts <- rep(NA, length(x))
Modified: pkg/R/rad.null.R
===================================================================
--- pkg/R/rad.null.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/rad.null.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -7,11 +7,14 @@
x <- as.rad(x)
nsp <- length(x)
wt <- rep(1, nsp)
- fit <- rev(cumsum(1/nsp:1)/nsp) * sum(x)
+ if (nsp > 1)
+ fit <- rev(cumsum(1/nsp:1)/nsp) * sum(x)
+ else
+ fit <- x
res <- dev.resids(x, fit, wt)
deviance <- sum(res)
residuals <- x - fit
- aic <- aicfun(x, wt, fit, wt, deviance)
+ aic <- aicfun(x, wt, fit, wt, deviance)
rdf <- nsp
p <- NA
names(p) <- "S"
@@ -21,4 +24,3 @@
class(out) <- c("radline", "glm")
out
}
-
Modified: pkg/R/rad.preempt.R
===================================================================
--- pkg/R/rad.preempt.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/rad.preempt.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -2,6 +2,8 @@
function (x, family = poisson, ...)
{
canfun <- function(p, x, ...) {
+ if (length(x) <= 1)
+ return(0)
p <- plogis(p)
if (p == 1)
p <- 1 - .Machine$double.eps
@@ -25,10 +27,14 @@
aic <- rdf <- deviance <- NA
p <- rep(NA, 1)
fit <- residuals <- prior.weights <- rep(NA, length(x))
- }
- else {
- p <- plogis(canon$estimate)
- fit <- exp(logJ + log(p) + log(1 - p) * rnk)
+ } else {
+ if (length(x) > 1) {
+ p <- plogis(canon$estimate)
+ fit <- exp(logJ + log(p) + log(1 - p) * rnk)
+ } else {
+ p <- 1
+ fit <- x
+ }
res <- dev.resids(x, fit, wt)
deviance <- sum(res)
residuals <- x - fit
@@ -42,4 +48,3 @@
class(out) <- c("radline", "glm")
out
}
-
Modified: pkg/R/rad.zipf.R
===================================================================
--- pkg/R/rad.zipf.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/rad.zipf.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -5,8 +5,18 @@
rnk <- seq(along = x)
off <- rep(log(sum(x)), length(x))
fam <- family(link = "log")
- ln <- try(glm(x ~ log(rnk) + offset(off), family = fam))
- if (inherits(ln, "try-error")) {
+ if (length(x) > 2)
+ ln <- try(glm(x ~ log(rnk) + offset(off), family = fam))
+ if (length(x) <= 2) {
+ aic <- 4
+ dev <- rdf <- 0
+ ln <- nl <- NA
+ p <- rep(NA, 2)
+ fit <- x
+ res <- rep(0, length(x))
+ wts <- rep(1, length(x))
+ }
+ else if (inherits(ln, "try-error")) {
aic <- rdf <- ln <- nl <- dev <- NA
p <- rep(NA, 2)
fit <- res <- wts <- rep(NA, length(x))
Modified: pkg/R/rad.zipfbrot.R
===================================================================
--- pkg/R/rad.zipfbrot.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/rad.zipfbrot.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -11,9 +11,19 @@
off <- rep(log(sum(x)), length(x))
p <- 0
fam <- family(link = "log")
- nl <- try(nlm(mandelfun, p = p, x = x, rnk = rnk, off = off,
- family = fam, hessian = TRUE, ...))
- if (inherits(nl, "try-error")) {
+ if (length(x) > 3)
+ nl <- try(nlm(mandelfun, p = p, x = x, rnk = rnk, off = off,
+ family = fam, hessian = TRUE, ...))
+ if (length(x) <= 3) {
+ aic <- 6
+ dev <- rdf <- 0
+ ln <- nl <- NA
+ p <- rep(NA, 3)
+ fit <- x
+ res <- rep(0, length(x))
+ wts <- rep(1, length(x))
+ }
+ else if (inherits(nl, "try-error")) {
aic <- rdf <- ln <- nl <- dev <- NA
p <- rep(NA, 3)
fit <- res <- wts <- rep(NA, length(x))
Modified: pkg/R/radfit.data.frame.R
===================================================================
--- pkg/R/radfit.data.frame.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/radfit.data.frame.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -1,6 +1,8 @@
"radfit.data.frame" <-
function(df, ...)
{
+ ## df *must* have rownames
+ rownames(df) <- rownames(df, do.NULL = TRUE)
out <- apply(df, 1, radfit, ...)
if (length(out) == 1)
out <- out[[1]]
Modified: pkg/R/radfit.default.R
===================================================================
--- pkg/R/radfit.default.R 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/R/radfit.default.R 2008-10-02 14:05:59 UTC (rev 516)
@@ -4,7 +4,7 @@
x <- as.rad(x)
NU <- rad.null(x, ...)
PE <- rad.preempt(x, ...)
- #BS <- rad.brokenstick(x, ...)
+ ##BS <- rad.brokenstick(x, ...)
LN <- rad.lognormal(x, ...)
ZP <- rad.zipf(x, ...)
ZM <- rad.zipfbrot(x, ...)
@@ -15,4 +15,3 @@
class(out) <- "radfit"
out
}
-
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2008-10-02 07:07:03 UTC (rev 515)
+++ pkg/inst/ChangeLog 2008-10-02 14:05:59 UTC (rev 516)
@@ -4,8 +4,11 @@
Version 1.16-1 (opened September 30, 2008)
- * print.radfit.frame: failed if the input data frame had no row
- names.
+ * radfit: Should work with empty sites (no species) or when the
+ number of species is less or equal the number of parameters
+ estimated (like may happen in sweeping analysis of data frames or
+ simulations). Takes care that input data frames have row names
+ which are necessary for displaying results.
* head.summary.cca & tail.summary.cca: shortcuts to
print(summary(x, ...), head, tail) for nicer Sweave tutorials.
More information about the Vegan-commits
mailing list