############################################################################
# R code belonging to
# Schuster, R., Germain R. R. & Roemer H. (2013) Multi-scale distribution and movement 
# effects of a montane highway on a large bodied- mammal community. PeerJ.
# Creative Commons Attribution license 3.0
###########################################################################
# forw.pscl.modlst: a function for "stepwise" regression using zeroinflated models of package pscl
# start.model: initial model e.g. zeroinfl(y~1, dist="poisson", link="logit", data=data)
# blocks: vector of covariate names e.g. var4 <- as.vector(c("woodZ","shrubZ")) 
# max.iter: number of iterations the function will run. If length(blocks) is e.g. 20
#           setting this parameter to 10 will result in models with max. 10 covariates
# AICcut: up to what delta AICc value should models to be cosidered as potential new base models
#
# Parts based on Forward.lmer by Rense Nieuwenhuis (http://www.rensenieuwenhuis.nl/r-sessions-32/)
# with some additions by for Nick Isaac
#
# A modified version of this code has been used in:
# Schuster, R., & Arcese, P. (2013). Using bird species community occurrence to 
# prioritize forests for old growth restoration. Ecography, 36(4), 499-507.
#
# Author: Richard Schuster (mail@richard-schuster.com)
# 12 July 2013
###########################################################################

# code needed for pscl package
addvar <- function(form, variable) formula(update(Formula(form),
as.formula(paste(". ~ . +", variable))))


forw.pscl.comb.modlst <- function(start.model, blocks, max.iter=NULL, AICcut = 1, dat = data, print.log=TRUE){

  modlst <- list(start.model)
  x <- 2
  coeff <- 1
	model.basis <- start.model
  keep <- list(start.model)
  AICmin <- AIClst <- AICc(start.model)
  #cutoff for when to exclude values
  cutoff <- 20
  #continue flag for debugging
  cont <- 0

	# Maximum number of iterations cannot exceed number of blocks, but this is also the default
	if (is.null(max.iter) | max.iter > length(blocks)) max.iter <- length(blocks)

	# Setting up the outer loop
	for(iter in 1:max.iter){
		models <- list()
    coeff <- coeff + 1
    cnt = 1
    for (jj in 1:length(keep)) {

  		# Iteratively updating the model with addition of one block of variable(s)
  		for(kk in 1:length(blocks)) {

        form <- as.formula(paste("~.+", blocks[kk]))

        if (class(dummy <- try(update(keep[[jj]], form))) == "zeroinfl") {
          flag <- 0 
          if (dummy$dist == "negbin")
            if (dummy$SE.logtheta == "NaN" || log(dummy$theta) > cutoff) flag <- 1
  
          if (flag == 0) {
            for (cc in 1:length(sqrt(diag(vcov(dummy))))) {
              if (diag(vcov(dummy))[[cc]] < 0 || sqrt(diag(vcov(dummy))[[cc]]) > cutoff) {
                flag <- 1
                break
              }
            }
          }
          if (flag == 0) { 
            for (bb in 1:length(AIClst)) {
              if (round(AICc(dummy),digits=6) == round(AIClst[bb], digits=6)) {
                flag <- 1
                break
              }
            }
          }
          if (flag == 0) {
            models[[cnt]] <- dummy
            modlst[[x]] <- models[[cnt]]
            AIClst <- c(AIClst, AICc(models[[cnt]]))
            x <- x + 1
            cnt <- cnt + 1
          }
  
        }
       }
    }

    if(length(LL <- unlist(lapply(models, function(x) {AICc(x)}))) == 0) break
    keep <- list()
    k = 1
    cont <- 0

    #check for improvement in AIC, if none stop loop
		for (mm in order(LL, decreasing=FALSE))
			{

      if (LL[mm]< AICmin + AICcut )
        {
        if (length(models[[mm]]$coefficients$count) == coeff) {
          keep[[k]]<- models[[mm]]
          k <- k + 1
          if (LL[mm]<AICmin) {
            AICmin <- LL[mm]
            cont <- 1
          }
        }
        }
     else break
     }
     
     if (length(keep) == 0) break                                      

    #for debugging
    #keep
    #cont
    
  }

#########################################################################
############## ZERO INTERCEPT ONLY ######################################
#########################################################################

  coeff <- 1
  tmp<-addvar (formula(start.model),"1|1")	
  if (start.model$dist == "negbin") dis="negbin"    
  if (start.model$dist == "poisson") dis="poisson"
  
  model.basis <- zeroinfl(tmp, dist = dis, link = "logit", data = dat)
  keep <- list(model.basis)
  AICmin <- AIClst <- AICc(model.basis)
  #continue flag for debugging
  cont <- 0

# Maximum number of iterations cannot exceed number of blocks, but this is also the default
	if (is.null(max.iter) | max.iter > length(blocks)) max.iter <- length(blocks)

	# Setting up the outer loop
	for(iter in 1:max.iter){
		models <- list()
    coeff <- coeff + 1
    cnt = 1
    for (jj in 1:length(keep)) {

  		# Iteratively updating the model with addition of one block of variable(s)
  		for(kk in 1:length(blocks)) {

        #form <- as.formula(paste("~.+", blocks[kk]))
    		tmp <- addvar(formula(keep[[jj]]), blocks[kk])
    
        if (class(dummy <- try(zeroinfl(tmp, dist = dis, link = "logit", data = dat))) == "zeroinfl") {
          flag <- 0 
          if (dummy$dist == "negbin")
            if (dummy$SE.logtheta == "NaN" || log(dummy$theta) > cutoff) flag <- 1
  
          if (flag == 0) {
            for (cc in 1:length(sqrt(diag(vcov(dummy))))) {
              if (diag(vcov(dummy))[[cc]] < 0 || sqrt(diag(vcov(dummy))[[cc]]) > cutoff) {
                flag <- 1
                break
              }
            }
          }
          if (flag == 0) { 
            for (bb in 1:length(AIClst)) {
              if (round(AICc(dummy),digits=6) == round(AIClst[bb], digits=6)) {
                flag <- 1
                break
              }
            }
          }
          if (flag == 0) {
            models[[cnt]] <- dummy
            modlst[[x]] <- models[[cnt]]
            AIClst <- c(AIClst, AICc(models[[cnt]]))
            x <- x + 1
            cnt <- cnt + 1
          }
  
        }
       }
    }

    if(length(LL <- unlist(lapply(models, function(x) {AICc(x)}))) == 0) break
    keep <- list()
    k = 1
    cont <- 0

    #check for improvement in AIC, if none stop loop
		for (mm in order(LL, decreasing=FALSE))
			{

      if (LL[mm]< AICmin + AICcut )
        {
        if (length(models[[mm]]$coefficients$count) == coeff) {
          keep[[k]]<- models[[mm]]
          k <- k + 1
          if (LL[mm]<AICmin) {
            AICmin <- LL[mm]
            cont <- 1
          }
        }
        }
     else break
     }
     
     if (length(keep) == 0) break                                      

    #for debugging
    #keep
    #cont
    
  }

  ## Create Model List
  modsel <- znblst(modlst)
	## Return the gathered output
  return(list(model=model.basis, modlst=modlst, modsel=modsel))
}

################################################################################
################################################################################


#Model list using the pscl package dataframe
# based on: setMethod("modSel", "unmarkedFitList", function(object, nullmod=NULL)
# from the R package unmarked
# and Model list using lm4 package dataframe
# Richard Schuster (mail@richard-schuster.com)
# 2012-08-09

znblst <- function(modlst)
{

    estList <- lapply(modlst, coef, altNames=TRUE)

    seList <- lapply(modlst, function(x) {
		se <- tryCatch(sqrt(diag(vcov(x, altNames=TRUE))),
			error=function(e) simpleError("Hessian is singular."))
        if(identical(class(se)[1], "simpleError")) {
            cat(se$message, fill=TRUE)
            se <- rep(NA, length(coef(x)))
            }
        names(se) <- names(sqrt(diag(vcov(x))))
        return(se)
        })
        

    eNames <- sort(unique(unlist(sapply(estList, names))))
    seNames <- paste("SE_", eNames, sep="")
    eseNames <- character(l <- length(c(eNames, seNames)))
    eseNames[seq(1, l, by=2)] <- eNames
    eseNames[seq(2, l, by=2)] <- seNames
    cNames <- c("formula", eseNames)
    out <- data.frame(matrix(NA, ncol=length(cNames), nrow=length(modlst)))
    colnames(out) <- cNames
    #out$model <- names(fits)
    out$formula <- sapply(modlst, function(x) {
          f <- as.character(formula(x))
          f <- paste(f[2], "~", f[3])
          f
        })
    for(i in 1:length(eNames)) {
        out[,eNames[i]] <- sapply(estList, function(x) x[eNames[i]])
        out[,seNames[i]] <- sapply(seList, function(x) x[eNames[i]])
        }

    if (modlst[[1]]$dist == "negbin") {
      out$theta <- sapply(modlst, function(x) x$theta)
      out$SE_theta <- sapply(modlst, function(x) x$SE.logtheta)
    }
    out$negLogLike <- sapply(modlst, function(x) logLik(x)[1])
    out$n <- sapply(modlst, function(x) x$n)
    out$converged <- sapply(modlst, function(x) x$converged)
    out$DF <- sapply(modlst, function(x) x$n - x$df.residual)    
    out$AICc <- sapply(modlst, function(x) AICc(x))
    out$deltaAIC <- out$AICc - min(out$AICc)
    out$AICwt <- exp(-out$delta / 2)
    out$AICwt <- out$AICwt / sum(out$AICwt)
    out <- out[order(out$AICc),]
    out$cumltvWt <- cumsum(out$AICwt)
    return(out)
}


################################################################################
################################################################################
# Function useage example
# Tried with R version 3.0.1 (2013-05-16)
# Test date 2013-07-12
library(pscl)
library(MuMIn)
library(Formula)

# data provided in Suppl. Material Data S1
hwy <- read.csv(file="Hwy_Modelling_Data.csv")
var1 <- as.vector(c("woodZ","shrubZ","herbZ","riparianZ","freshwZ"))

znbch<- zeroinfl(Carnivora ~1, dist = "poisson", link = "logit",  data = hwy)
carnh <- forw.pscl.comb.modlst(start.model=znbch,blocks=var1,max.iter=4, AICcut = 2,dat=hwy)
write.csv(carnh$modsel, "Carnivora_Hwy.csv", row.names = FALSE)