Initialize maxnet repository

Initial set up of maxnet repository

.Rbuildignore

@@ -0,0 +1 @@
+^\.travis\.yml$

.gitignore

@@ -0,0 +1,56 @@
+# class files, jar files, zip files, emacs backup files, semantic cache and javadoc files
+*.class
+*.jar
+*.zip
+*~
+semantic.cache
+html/
+help.html
+
+# Windows image file caches
+Thumbs.db
+ehthumbs.db
+
+# Folder config file
+Desktop.ini
+
+# Recycle Bin used on file shares
+$RECYCLE.BIN/
+
+# Windows Installer files
+*.cab
+*.msi
+*.msm
+*.msp
+
+# Windows shortcuts
+*.lnk
+
+# =========================
+# Operating System Files
+# =========================
+
+# OSX
+# =========================
+
+.DS_Store
+.AppleDouble
+.LSOverride
+
+# Thumbnails
+._*
+
+# Files that might appear in the root of a volume
+.DocumentRevisions-V100
+.fseventsd
+.Spotlight-V100
+.TemporaryItems
+.Trashes
+.VolumeIcon.icns
+
+# Directories potentially created on remote AFP share
+.AppleDB
+.AppleDesktop
+Network Trash Folder
+Temporary Items
+.apdisk

.travis.yml

@@ -0,0 +1,5 @@
+# R for travis: see documentation at https://docs.travis-ci.com/user/languages/r
+
+language: R
+sudo: false
+cache: packages

DESCRIPTION

@@ -0,0 +1,12 @@
+Package: maxnet
+Type: Package
+Title: An implementation of Maxent using glmnet
+Version: 0.1
+Date: 2013-06-06
+Author: Steven Phillips
+Maintainer: Steven Phillips <mrmaxent@gmail.com>
+Imports:
+    glmnet
+Description: An implementation in R of the Maxent method for modeling species distributions from occurrence records, using glmnet for model fitting. It implements the same feature types and regularization options as the Maxent Java package.  See the Maxent website for more details.
+License: MIT + file LICENSE
+RoxygenNote: 5.0.1

LICENSE

@@ -0,0 +1,2 @@
+YEAR: 2016
+COPYRIGHT HOLDER: Steven Phillips

NAMESPACE

@@ -0,0 +1,12 @@
+# Generated by roxygen2: do not edit by hand
+
+S3method(plot,maxnet)
+S3method(predict,maxnet)
+export(categorical)
+export(hinge)
+export(maxnet)
+export(maxnet.default.regularization)
+export(maxnet.formula)
+export(response.plot)
+export(thresholds)
+import(stats)

R/categorical.R

@@ -0,0 +1,8 @@
+#' @export
+categorical <-
+function(x)
+{
+   f <- outer(x, levels(x), function(w,f) ifelse(w==f,1,0))
+   colnames(f) <- paste("", levels(x), sep=":")
+   f
+}

R/categoricalval.R

@@ -0,0 +1,5 @@
+categoricalval <-
+function(x, category)
+{
+   ifelse(x==category, 1, 0)
+}

R/data.R

@@ -0,0 +1,12 @@
+#' Occurrence records and background data for the brown-throated three-toed sloth, Bradypus variegatus
+#'
+#' A dataset containing environmental data at 116 Bradypus variegatus occurrence points
+#' and 10000 background points in South and Central America.  Occurrence data are from
+#' Anderson and Handley (2001); see Phillips et al. (2006) for descriptions of the 
+#' predictor variables.
+#'
+#' @references 
+#' Anderson, R. P. and Handley, Jr., C. O. (2001). A new species of three-toed sloth (Mammalia: Xenarthra) from Panama, with a review of the genus Bradypus. Proceedings of the Biological Society of Washington 114, 1-33.
+#'
+#' Phillips, S. J. et al. (2006). Maximum entropy modeling of species geographic distributions. Ecological Modelling 190, 231-259
+"bradypus"

R/hinge.R

@@ -0,0 +1,13 @@
+#' @export
+hinge <-
+function(x, nknots=50)
+{
+   min <- min(x)
+   max <- max(x)
+   k <- seq(min, max, length=nknots)
+   lh <- outer(x, utils::head(k,-1), function(w,h) hingeval(w, h, max))
+   rh <- outer(x, k[-1], function(w,h) hingeval(w, min, h))
+   colnames(lh) <- paste("", utils::head(k,-1), max, sep=":")
+   colnames(rh) <- paste("", min, k[-1], sep=":")
+   cbind(lh, rh)
+}

R/hingeval.R

@@ -0,0 +1,5 @@
+hingeval <-
+function(x, min, max)
+{
+   pmin(1, pmax(0, (x-min)/(max-min)))
+}

R/maxnet.R

@@ -0,0 +1,33 @@
+#' @import stats
+#' @export
+maxnet <-
+function(p, data, f=maxnet.formula(p, data), regmult=1.0, 
+         regfun=maxnet.default.regularization, ...)
+{
+   mm <- model.matrix(f, data)
+   reg <- regfun(p,mm) * regmult
+   weights <- p+(1-p)*100
+   glmnet::glmnet.control(pmin=1.0e-8, fdev=0)  
+   model <- glmnet::glmnet(x=mm, y=as.factor(p), family="binomial", standardize=F, penalty.factor=reg, lambda=10^(seq(4,0,length.out=100))*sum(reg)/length(reg)*sum(p)/sum(weights), weights=weights, ...)
+   class(model) <- c("maxnet", class(model))
+   bb <- model$beta[,100]
+   model$betas <- bb[bb!=0]
+   model$alpha <- 0
+   rr <- predict.maxnet(model, data[p==0,], type="exponent", clamp=F)
+   raw <- rr / sum(rr)
+   model$entropy <- -sum(raw * log(raw))
+   model$alpha <- -log(sum(rr))
+   model$penalty.factor <- reg
+   model$featuremins <- apply(mm, 2, min)
+   model$featuremaxs <- apply(mm, 2, max)
+   vv <- (sapply(data, class)!="factor")
+   model$varmin <- apply(data[,vv], 2, min)
+   model$varmax <- apply(data[,vv], 2, max)
+   means <- apply(data[p==1,vv], 2, mean)
+   majorities <- sapply(names(data)[!vv], 
+      function(n) which.max(table(data[p==1,n])))
+   names(majorities) <- names(data)[!vv]
+   model$samplemeans <- c(means, majorities)
+   model$levels <- lapply(data, levels)
+   model
+}

R/maxnet.default.regularization.R

@@ -0,0 +1,40 @@
+#' @export
+maxnet.default.regularization <-
+function(p, m)
+{
+   isproduct <- function(x) grepl(":", x) & !grepl("\\(", x)
+   isquadratic <- function(x) grepl("^I\\(.*\\^2\\)", x)
+   ishinge <- function(x) grepl("^hinge\\(", x)
+   isthreshold <- function(x) grepl("^thresholds\\(", x)
+   iscategorical <- function(x) grepl("^categorical\\(", x)
+   regtable <- function(name, default) {
+      if (ishinge(name)) return(list(c(0,1), c(0.5,0.5)))
+      if (iscategorical(name)) return(list(c(0,10,17), c(0.65, 0.5, 0.25)))
+      if (isthreshold(name)) return(list(c(0,100), c(2.0, 1.0)))
+      default
+   }
+   lregtable <- list(c(0,10,30,100), c(1,1,0.2,0.05))
+   qregtable <- list(c(0,10,17,30,100), c(1.3,0.8,0.5,0.25,0.05))
+   pregtable <- list(c(0,10,17,30,100), c(2.6,1.6,0.9,0.55,0.05))
+   mm <- m[p==1,]
+   np <- nrow(mm)
+   lqpreg <- lregtable
+   if (sum(isquadratic(colnames(mm)))) lqpreg <- qregtable
+   if (sum(isproduct(colnames(mm)))) lqpreg <- pregtable
+   classregularization <- sapply(colnames(mm), function(n) {
+      t <- regtable(n, lqpreg)
+      approx(t[[1]], t[[2]], np, rule=2)$y
+   }) / sqrt(np)
+   # increase regularization for extreme hinges
+   ishinge <- grepl("^hinge\\(", colnames(mm))
+   hmindev <- sapply(1:ncol(mm), function(i) {
+      if (!ishinge[i]) return(0)
+      avg <- mean(mm[,i])
+      std <- max(sd(mm[,i]), 1/sqrt(np))
+      std*.5/sqrt(np)
+   })
+   # increase reg'n for threshold features that are all 1 or 0 on presences
+   tmindev <- sapply(1:ncol(mm), function(i) {
+      ifelse(isthreshold(colnames(mm)[i]) && (sum(mm[,i])==0 || sum(mm[,i])==nrow(mm)), 1,0)})
+   pmax(0.001 * (apply(m,2,max)-apply(m,2,min)), hmindev, tmindev, apply(as.matrix(mm), 2, sd) * classregularization)
+}

R/maxnet.formula.R

@@ -0,0 +1,32 @@
+#' @export
+maxnet.formula <-
+function(p, data, classes="default")
+{
+   cont <- names(data)[sapply(data,class)!="factor"]
+   categorical <- names(data)[sapply(data,class)=="factor"]
+   np <- sum(p)
+   if (classes=="default") {
+      if (np < 10) classes <- "l"
+      else if (np < 15) classes <- "lq"
+      else if (np < 80) classes <- "lqh"
+      else classes <- "lqpht"
+   }
+   terms <- NULL
+   if (length(cont)) {
+      if (grepl("l", classes))
+         terms <- c(terms, paste(cont,collapse=" + "))
+      if (grepl("q", classes))
+         terms <- c(terms, paste("I(",cont,"^2)",sep="",collapse=" + "))
+      if (grepl("h", classes))
+         terms <- c(terms, paste("hinge(",cont,")",sep="",collapse=" + "))
+      if (grepl("t", classes))
+         terms <- c(terms, paste("thresholds(",cont,")",sep="",collapse=" + "))
+      if (grepl("p", classes)) {
+         m <- outer(cont, cont, function(x,y) paste(x,y,sep=":"))
+         terms <- c(terms, m[lower.tri(m)])
+      }
+   }
+   if (length(categorical))
+      terms <- c(terms, paste("categorical(",categorical,")",sep="",collapse=" + "))
+   formula(paste("~", paste(terms, collapse = " + "), "-1"))
+}

R/plot.maxnet.R

@@ -0,0 +1,17 @@
+#' @export
+plot.maxnet <-
+function(x, vars=names(x$samplemeans), common.scale=T, type=c("link","exponential","cloglog","logistic"), ylab=NULL, ...)
+{
+   type <- match.arg(type)
+   nc <- ceiling(sqrt(length(vars)))
+   nr <- ceiling(length(vars)/nc)
+   graphics::par(mfrow=c(nr,nc), mar=c(5,5,4,2)+.1)
+   ylim=NULL
+   if (common.scale && (type=="link" || type=="exponential")) {
+vals <- do.call(c, lapply(vars, function(v) 
+          response.plot(x, v, type, plot=F)))
+        ylim=c(min(vals), max(vals))
+   }
+   if (type=="cloglog" || type=="logistic") ylim=c(0,1)
+   for (v in vars) response.plot(x, v, type, ylim=ylim, ylab=ylab)
+}

R/predict.maxnet.R

@@ -0,0 +1,23 @@
+#' @export
+predict.maxnet <-
+function(object, newdata, clamp=T, type=c("link","exponential","cloglog","logistic"), ...)
+{
+   if (clamp) {
+      for (v in intersect(names(object$varmax), names(newdata))) {
+         newdata[,v] <- pmin(pmax(newdata[,v], object$varmin[v]), object$varmax[v])
+      }
+   }
+   terms <- sub("hinge\\((.*)\\):(.*):(.*)$", "hingeval(\\1,\\2,\\3)", names(object$betas))
+   terms <- sub("categorical\\((.*)\\):(.*)$", "categoricalval(\\1,\\2)", terms)
+   terms <- sub("thresholds\\((.*)\\):(.*)$", "thresholdval(\\1,\\2)", terms)
+   f <- formula(paste("~", paste(terms, collapse=" + "), "-1"))
+   mm <- model.matrix(f, data.frame(newdata))
+   if (clamp) mm <- t(pmin(pmax(t(mm), object$featuremins[names(object$betas)]), 
+                 object$featuremaxs[names(object$betas)]))
+   link <- (mm %*% object$betas) + object$alpha
+   type <- match.arg(type)
+   if (type=="link") return(link)
+   if (type=="exponential") return(exp(link))
+   if (type=="cloglog") return(1-exp(0-exp(object$entropy+link)))
+   if (type=="logistic") return(1/(1+exp(-object$entropy-link)))
+}

R/response.plot.R

@@ -0,0 +1,20 @@
+#' @export
+response.plot <-
+function(mod, v, type, mm=mod$samplemeans, min=mod$varmin[v], max=mod$varmax[v], levels=unlist(mod$levels[v]), plot=T, ylim=NULL, ylab=NULL) {
+   nr <- if (is.null(levels)) 100 else length(levels)
+   m <- data.frame(matrix(mm,nr,length(mm),byrow=T))
+   colnames(m) <- names(mm)
+   m[,v] <- if (!is.null(levels)) levels else 
+      seq(min - 0.1*(max-min), max+0.1*(max-min), length=100)
+   preds <- predict(mod, m, type=type)
+   if (is.null(ylab))
+      ylab <- paste(toupper(substring(type,1,1)), substring(type, 2), sep="")
+   if (plot) {
+      if (is.null(levels)) {
+         plot(m[,v], preds, xlab=v, ylab=ylab, type="l", ylim=ylim) 
+      } else {
+         graphics::barplot(as.vector(preds), names.arg=levels, xlab=v, ylab=ylab, ylim=ylim)
+      }
+   }
+   else return(preds)
+}

R/thresholds.R

@@ -0,0 +1,11 @@
+#' @export
+thresholds <-
+function(x, nknots=50)
+{
+   min <- min(x)
+   max <- max(x)
+   k <- seq(min, max, length=nknots+2)[2:nknots+1]
+   f <- outer(x, k, function(w,t) ifelse(w>=t,1,0))
+   colnames(f) <- paste("",k, sep=":")
+   f
+}

R/thresholdval.R

@@ -0,0 +1,5 @@
+thresholdval <-
+function(x, knot) 
+{
+   ifelse(x >= knot, 1, 0)
+}

README.md

@@ -0,0 +1,5 @@
+# maxnet
+Maxent is a stand-alone Java application for modelling species geographic distributions.  This open source repository contains an R package, called "maxnet", which implements much of the functionality of the Java application.  We welcome contributions to maxnet.
+The current release of maxnet is also available for download on the CRAN website.  
+
+For information on the Maxent application, please see the Maxent home page at the American Museum of Natural History.  

man/hinge.Rd

@@ -0,0 +1,32 @@
+\name{hinge}
+\alias{hinge}
+\alias{hingeval}
+\alias{thresholds}
+\alias{thresholdval}
+\alias{categorical}
+\alias{categoricalval}
+\title{Maxent feature classes}
+\description{Create and evaluate Maxent's feature classes}
+\usage{
+hinge(x, nknots = 50)
+thresholds(x, nknots=50)
+categorical(x)
+hingeval(x, min, max)
+thresholdval(x, knot)
+categoricalval(x, category)
+}
+\arguments{
+  \item{x}{ a predictor: a factor for categorical and categoricalval, otherwise numeric. }
+  \item{nknots}{ number of knots. }
+  \item{min}{ value of \code{x} at which hinge feature is 0. }
+  \item{max}{ value of \code{x} at which hinge feature is 1. }
+  \item{knot}{ value of \code{x} at which threshold feature changes from 0 to 1. }
+  \item{category}{ a level of \code{x} }
+}
+\value{
+  \code{hinge}, \code{threshold} and \code{categorical} return a matrix with a column for each feature of the specified type.  \code{hinge} creates \code{2*nknots} hinge features, half with \code{min=min(x)} and half with \code{max=max(x)}, and knots evenly spaced between \code{min(x)} and \code{max(x)}.  \code{hingeval} calculates a single hinge feature: 0 if the predictor is below min, 1 if the predictor is above max, and linearly interpolated inbetween.  \code{thresholdval} calculates a single threshold feature: 1 if the predictor is above the knot, 0 otherwise.  \code{categoricalval} calculates a categorical feature:  1 if the predictor matches the category.
+}
+\author{
+  Steven Phillips
+}
+\keyword{Maxent}