Merge pull request #99 from piLaboratory/dev

Sads 0.2.4

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 <[email protected]>
 Depends:

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)

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.
+

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]

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"))

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, ...))

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
   }

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.
+}
+

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.
+}
+

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