diff --git a/DESCRIPTION b/DESCRIPTION index 20addf8..82a11e7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: sads Type: Package Title: Maximum Likelihood Models for Species Abundance Distributions -Version: 0.2.3 -Date: 2015-09-24 +Version: 0.2.4 +Date: 2015-11-02 Author: Paulo I. Prado, Murilo Dantas Miranda and Andre Chalom Maintainer: Paulo I. Prado Depends: diff --git a/NAMESPACE b/NAMESPACE index be1b584..1365b01 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -18,10 +18,12 @@ export(fitrad, fitsad, octav, octavpred, radpred, pprad, ppsad, qqrad, qqsad, ra ## Accessory functions export(plotprofmle, pred.logser, dtrunc, ptrunc, qtrunc, rsad, updatesad, updaterad, distr) ## Explicit classes and methods export -exportClasses(rad, octav, fitsad, fitrad) +exportClasses(rad, octav, fitsad, fitrad, summary.sads) exportMethods(plot, points, lines, octavpred, radpred, qqsad, qqrad, plotprofmle) # Methods that override bbmle counterparts: -exportMethods(show, nobs, AIC, AICc) +exportMethods(show, summary, nobs, AIC, AICc) +# Standard stats methods +exportMethods(coefficients, fitted, fitted.values, residuals) ## IMPORTS ## ## Import specific functions from other packages used by new methods importFrom("graphics", plot, points, lines, abline, axis, hist, par) diff --git a/NEWS b/NEWS index 4a2f918..3b94ce8 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,12 @@ +0.2.4 (2015-11-02) +* octavpred now accepts a 'preston' argument (as does 'octav') +* fixes a bug in fitvolkov under Windows +* added a friendlier error message when the octavs provided in octavpred do not cover the data +* deprecated the slot "distr" on fitsad/fitrad objects. Now users should use the `distr` function (see ?distr) +* fixed signature errors on several methods dealing with the `trunc` parameter +* now providing the coefficients, fitted, fitted.values, residuals methods on fitsad/fitrad methods +* providing a summary method for fitsad/fitrad classes that allows for fixed parameters (such as produced by fitls, fitvolkov, etc) + 0.2.3 (2015-09-24) * fitgeom, fitpoilog and fitnbinom now can get argument trunc=NULL to avoid zero-truncation (the default) * fitgeom uses default optimization method from mle2 instead of Brent's method, to avoid some overflow errors for large datasets (e.g. moths) @@ -34,3 +43,4 @@ error handling and large performance gains. 0.1.10 (2014-07-02) * Initial release on CRAN. + diff --git a/R/fitvolkov.R b/R/fitvolkov.R index 7a98330..5bb692a 100644 --- a/R/fitvolkov.R +++ b/R/fitvolkov.R @@ -1,9 +1,11 @@ fitvolkov <- function(x, trunc, start.value, ...){ dots <- list(...) if(missing(start.value)){ - sink("/dev/null") # as the following function outputs lots of garbage... - start.value = maxLikelihood.ESF(c(5, 0.5), x)$par + tmp <- tempfile() + sink(tmp) # as the following function outputs lots of garbage... + start.value <- maxLikelihood.ESF(c(5, 0.5), x)$par sink() + file.remove(tmp) ## as sink("dev/null") does not work in al OS' } thetahat <- start.value[1] mhat <-start.value[2] diff --git a/R/sads-classes.R b/R/sads-classes.R index 1f0ecf5..b266c96 100644 --- a/R/sads-classes.R +++ b/R/sads-classes.R @@ -21,6 +21,15 @@ setClass("rad", representation("data.frame"), validity = function(object) { setClass("fitsad", representation("mle2", sad="character", distr="character", trunc="numeric")) setClass("fitrad", representation("mle2", rad="character", distr="character", trunc="numeric", rad.tab="rad")) -#setClass("fitsadlist", representation("list")) distr.depr <- "The 'distr' slot of fitrad and fitsad objects have been deprecated. Please see ?distr" + +#' Summary for fitsad/fitrad calls +#' +#' This function works almost exactly as bbmle's summary.mle2, but it includes a "fixed parameters" +#' line for models with fixed parameters, such as \code{\link{fitls}} or \code{fitvolkov}. +#' NOTICE that the summary.mle2 is redefined in this package, as this class is not exported by bbmle. +#' @rdname summary.sads-class +setClass("summary.mle2", representation(call = "language", coef = "matrix",m2logL = "numeric", fixed="numeric")) +#' @rdname summary.sads-class +setClass("summary.sads", representation("summary.mle2", fixed="numeric")) diff --git a/R/sads-methods.R b/R/sads-methods.R index 6c226f7..d3c471a 100644 --- a/R/sads-methods.R +++ b/R/sads-methods.R @@ -233,6 +233,37 @@ showmle2 <- function(object) { setMethod("show", "fitsad", function(object){showmle2(object)}) setMethod("show", "fitrad", function(object){showmle2(object)}) +#### summary class dealing with fixed parameters (such as fitls, fitvolkov, etc) +#' @rdname summary.sads-class +#' @param object An object of class fitsad/fitrad is required to generate a summary.sads object. +setMethod("show", "summary.sads", function(object){ + cat("Maximum likelihood estimation\n\nCall:\n") + print(object@call) + cat("\nCoefficients:\n") + printCoefmat(object@coef) + if (length(object@fixed) > 0) { + cat("\nFixed parameters:\n") + print(object@fixed) + } + cat("\n-2 log L:", object@m2logL, "\n") + }) +sumle2 <- function(object, ...){ + cmat <- cbind(Estimate = object@coef, + `Std. Error` = sqrt(diag(object@vcov))) + zval <- cmat[,"Estimate"]/cmat[,"Std. Error"] + pval <- 2*pnorm(-abs(zval)) + coefmat <- cbind(cmat,"z value"=zval,"Pr(z)"=pval) + m2logL <- 2*object@min + fixed <- numeric() + if (! all(object@fullcoef %in% object@coef)) + fixed <- object@fullcoef [! object@fullcoef %in% object@coef] + new("summary.sads", call=object@call.orig, coef=coefmat, fixed=fixed, m2logL= m2logL) +} +#' @rdname summary.sads-class +setMethod("summary", "fitsad", function(object){sumle2(object)}) +#' @rdname summary.sads-class +setMethod("summary", "fitrad", function(object){sumle2(object)}) + ## radpred generic functions and methods ### setGeneric("radpred", def = function(object, sad, rad, coef, trunc , distr=NA, S, N) standardGeneric("radpred") @@ -441,7 +472,9 @@ setMethod("octavpred", signature(object="missing",sad="missing", rad="character" drad <- get(paste("d",rad,sep=""),mode="function") ab <- do.call(drad, c(list(x=1:S),coef,dots))*N } - Y = hist(ab, breaks=c(2^(min(oct)-2),n), plot=FALSE) + tryCatch({Y = hist(ab, breaks=c(2^(min(oct)-2),n), plot=FALSE)}, + error = function(cond) stop("Octaves do not span the entire range, try using a larger oct argument (maybe negative octaves?)") + ) res <- data.frame(octave = oct, upper = n, Freq = Y$count) if(preston) res <- prestonfy(res, ceiling(ab)) new("octav", res) @@ -680,3 +713,59 @@ setMethod("pprad", pprad(x=y, rad=rad, coef=coef, trunc=trunc, plot=plot, line=line, ...) } ) + +### Providing standard stats methods +#' Standard stats methods +#' +#' Provide the standard interface for fitted objects +#' +#' These methods are provided to allow for standard manipulation of \code{\link{fitsad}} +#' and \code{\link{fitrad}} objects using the generic methods defined in the "stats" package. +#' Please see the original man pages for each method. +#' +#' \code{coefficients} is an alias to \code{\link[stats]{coef}} (implemented in package "bbmle"). +#' +#' \code{fitted} and \code{fitted.values} provide an alternative interface to \code{\link{radpred}}; +#' these are also used to calcutate \code{residuals}. +#' +#' Notice that radpred is a preferred interface for most calculations, specially if there are several +#' ties. +#' +#' @param object An object from class fitsad or fitrad +#' @param \dots Other arguments to be forwarded for the lower level function +#' @rdname stats +setMethod("coefficients", signature(object="fitsad"), + function(object, ...) bbmle::coef(object, ...)) +#' @rdname stats +setMethod("coefficients", signature(object="fitrad"), + function(object, ...) bbmle::coef(object, ...)) +#' @rdname stats +setMethod("fitted.values", signature(object="fitsad"), + function(object, ...) fitted(object, ...)) +#' @rdname stats +setMethod("fitted.values", signature(object="fitrad"), + function(object, ...) fitted(object, ...)) +#' @rdname stats +setMethod("fitted", signature(object="fitsad"), + function(object, ...) { + rad <- radpred(object)$abund + rad <- rad[rev(order(object@data$x))] + names(rad) <- as.character(1:length(rad)) + return(rad) + } + ) +#' @rdname stats +setMethod("fitted", signature(object="fitrad"), + function(object, ...) { + rad <- radpred(object)$abund + rad <- rad[order(object@rad.tab$abund)] + names(rad) <- as.character(1:length(rad)) + return(rad) + } + ) +#' @rdname stats +setMethod("residuals", signature(object="fitsad"), + function(object, ...) object@data$x - fitted(object, ...)) +#' @rdname stats +setMethod("residuals", signature(object="fitrad"), + function(object, ...) object@rad.tab$abund - fitted(object, ...)) diff --git a/man/sads-package.Rd b/man/sads-package.Rd index 8316fd7..d325b7d 100644 --- a/man/sads-package.Rd +++ b/man/sads-package.Rd @@ -16,8 +16,8 @@ \tabular{ll}{ Package: \tab sads\cr Type: \tab Package\cr - Version: \tab 0.2.3\cr - Date: \tab 2015-09-24\cr + Version: \tab 0.2.4\cr + Date: \tab 2015-11-02\cr License: \tab gpl \cr LazyLoad: \tab yes\cr } diff --git a/man/stats.Rd b/man/stats.Rd new file mode 100644 index 0000000..45c1978 --- /dev/null +++ b/man/stats.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2 (4.1.1): do not edit by hand +% Please edit documentation in R/sads-methods.R +\docType{methods} +\name{coefficients,fitsad-method} +\alias{coefficients,fitrad-method} +\alias{coefficients,fitsad-method} +\alias{fitted,fitrad-method} +\alias{fitted,fitsad-method} +\alias{fitted.values,fitrad-method} +\alias{fitted.values,fitsad-method} +\alias{residuals,fitrad-method} +\alias{residuals,fitsad-method} +\title{Standard stats methods} +\usage{ +\S4method{coefficients}{fitsad}(object, ...) + +\S4method{coefficients}{fitrad}(object, ...) + +\S4method{fitted.values}{fitsad}(object, ...) + +\S4method{fitted.values}{fitrad}(object, ...) + +\S4method{fitted}{fitsad}(object, ...) + +\S4method{fitted}{fitrad}(object, ...) + +\S4method{residuals}{fitsad}(object, ...) + +\S4method{residuals}{fitrad}(object, ...) +} +\arguments{ +\item{object}{An object from class fitsad or fitrad} + +\item{\dots}{Other arguments to be forwarded for the lower level function} +} +\description{ +Provide the standard interface for fitted objects +} +\details{ +These methods are provided to allow for standard manipulation of \code{\link{fitsad}} +and \code{\link{fitrad}} objects using the generic methods defined in the "stats" package. +Please see the original man pages for each method. + +\code{coefficients} is an alias to \code{\link[stats]{coef}} (implemented in package "bbmle"). + +\code{fitted} and \code{fitted.values} provide an alternative interface to \code{\link{radpred}}; +these are also used to calcutate \code{residuals}. + +Notice that radpred is a preferred interface for most calculations, specially if there are several +ties. +} + diff --git a/man/summary.sads-class.Rd b/man/summary.sads-class.Rd new file mode 100644 index 0000000..feaee95 --- /dev/null +++ b/man/summary.sads-class.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2 (4.1.1): do not edit by hand +% Please edit documentation in R/sads-classes.R, R/sads-methods.R +\docType{class} +\name{summary.mle2-class} +\alias{show,summary.sads-method} +\alias{summary,fitrad-method} +\alias{summary,fitsad-method} +\alias{summary.mle2-class} +\alias{summary.sads-class} +\title{Summary for fitsad/fitrad calls} +\usage{ +\S4method{show}{summary.sads}(object) + +\S4method{summary}{fitsad}(object) + +\S4method{summary}{fitrad}(object) +} +\arguments{ +\item{object}{An object of class fitsad/fitrad is required to generate a summary.sads object.} +} +\description{ +This function works almost exactly as bbmle's summary.mle2, but it includes a "fixed parameters" +line for models with fixed parameters, such as \code{\link{fitls}} or \code{fitvolkov}. +NOTICE that the summary.mle2 is redefined in this package, as this class is not exported by bbmle. +} + diff --git a/vignettes/sads_intro.Rnw b/vignettes/sads_intro.Rnw index 92e592f..5a48d3d 100644 --- a/vignettes/sads_intro.Rnw +++ b/vignettes/sads_intro.Rnw @@ -367,7 +367,7 @@ sample with zeroes omitted: meanlog=5, sdlog=2)) @ -Once this is a Poisson sample of a lognormal community, the abundances +Since this is a Poisson sample of a lognormal community, the abundances in the sample should follow a Poisson-lognormal distribution with parameters $\mu + \log a $ and $\sigma$ \citep{grotan2008}. We can check this by fitting a Poisson-lognormal