Package 'unmarked'

Description

Fits hierarchical models of animal occurrence and abundance to data collected on species that may be detected imperfectly. Models include single- and multi-season site occupancy models, binomial N-mixture models, and multinomial N-mixture models. The data can arise from survey methods such as occurrence sampling, temporally replicated counts, removal sampling, double observer sampling, and distance sampling. Parameters governing the state and observation processes can be modeled as functions of covariates. General treatment of these models can be found in MacKenzie et al. (2006) and Royle and Dorazio (2008). The primary reference for the package is Fiske and Chandler (2011).

Details

Overview of Model-fitting Functions:

occu fits occurrence models with no linkage between abundance and detection (MacKenzie et al. 2002).

occuRN fits abundance models to presence/absence data by exploiting the link between detection probability and abundance (Royle and Nichols 2003).

occuFP fits occupancy models to data characterized by false negatives and false positive detections (e.g., Royle and Link [2006] and Miller et al. [2011]).

occuMulti fits multi-species occupancy model of Rota et al. [2016].

colext fits the mutli-season occupancy model of MacKenzie et al. (2003).

pcount fits N-mixture models (aka binomial mixture models) to repeated count data (Royle 2004a, Kery et al 2005).

distsamp fits the distance sampling model of Royle et al. (2004) to distance data recorded in discrete intervals.

gdistsamp fits the generalized distance sampling model described by Chandler et al. (2011) to distance data recorded in discrete intervals.

gpcount fits the generalized N-mixture model described by Chandler et al. (2011) to repeated count data collected using the robust design.

multinomPois fits the multinomial-Poisson model of Royle (2004b) to data collected using methods such as removal sampling or double observer sampling.

gmultmix fits a generalized form of the multinomial-mixture model of Royle (2004b) that allows for estimating availability and detection probability.

pcountOpen fits the open population model of Dail and Madsen (2011) to repeated count data. This is a genearlized form of the Royle (2004a) N-mixture model that includes parameters for recruitment and apparent survival.

Data: All data are passed to unmarked's estimation functions as a formal S4 class called an unmarkedFrame, which has child classes for each model type. This allows metadata (eg as distance interval cut points, measurement units, etc...) to be stored with the response and covariate data. See unmarkedFrame for a detailed description of unmarkedFrames and how to create them.

Model Specification: unmarked's model-fitting functions allow specification of covariates for both the state process and the detection process. For two-level hierarchical models, (eg occu, occuRN, pcount, multinomPois, distsamp) covariates for the detection process (at the site or observation level) and the state process (at the site level) are specified with a double right-hand sided formula, in that order. Such a formula looks like

$~ x1 + x2 + \ldots + x_n ~ x_1 + x_2 + \ldots + x_n$

where $x_1$ through $x_n$ are additive covariates of the process of interest. Using two tildes in a single formula differs from standard R convention, but it is informative about the model being fit. The meaning of these covariates, or what they model, is full described in the help files for the individual functions and is not the same for all functions. For models with more than two processes (eg colext, gmultmix, pcountOpen), single right-hand sided formulas (only one tilde) are used to model each parameter.

Utility Functions: unmarked contains several utility functions for organizing data into the form required by its model-fitting functions. csvToUMF converts an appropriately formated comma-separated values (.csv) file to a list containing the components required by model-fitting functions.

Author(s)

Ian Fiske, Richard Chandler, Andy Royle, Marc Kery, David Miller, and Rebecca Hutchinson

References

Chandler, R. B., J. A. Royle, and D. I. King. 2011. Inference about density and temporary emigration in unmarked populations. Ecology 92:1429-1435.

Dail, D. and L. Madsen. 2011. Models for estimating abundance from repeated counts of an open metapopulation. Biometrics 67:577-587.

Fiske, I. and R. B. Chandler. 2011. unmarked: An R package for fitting hierarchical models of wildlife occurrence and abundance. Journal of Statistical Software 43:1–23.

Kery, M., Royle, J. A., and Schmid, H. 2005 Modeling avian abundance from replicated counts using binomial mixture models. Ecological Applications 15:1450–1461.

MacKenzie, D. I., J. D. Nichols, G. B. Lachman, S. Droege, J. A. Royle, and C. A. Langtimm. 2002. Estimating site occupancy rates when detection probabilities are less than one. Ecology 83: 2248–2255.

MacKenzie, D. I., J. D. Nichols, J. E. Hines, M. G. Knutson, and A. B. Franklin. 2003. Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84:2200–2207.

MacKenzie, D. I., J. D. Nichols, J. A. Royle, K. H. Pollock, L. L. Bailey, and J. E. Hines. 2006. Occupancy Estimation and Modeling. Amsterdam: Academic Press.

Miller, D.A., J.D. Nichols, B.T. McClintock, E.H.C. Grant, L.L. Bailey, and L.A. Weir. 2011. Improving occupancy estimation when two types of observational error occur: non-detection and species misidentification. Ecology 92:1422-1428.

Rota, C.T., et al. 2016. A multi-species occupancy model for two or more interacting species. Methods in Ecology and Evolution 7: 1164-1173.

Royle, J. A. 2004a. N-Mixture models for estimating population size from spatially replicated counts. Biometrics 60:108–105.

Royle, J. A. 2004b. Generalized estimators of avian abundance from count survey data. Animal Biodiversity and Conservation 27:375–386.

Royle, J. A., D. K. Dawson, and S. Bates. 2004. Modeling abundance effects in distance sampling. Ecology 85:1591–1597.

Royle, J. A., and R. M. Dorazio. 2006. Hierarchical models of animal abundance and occurrence. Journal Of Agricultural Biological And Environmental Statistics 11:249–263.

Royle, J.A., and W.A. Link. 2006. Generalized site occupancy models allowing for false positive and false negative errors. Ecology 87:835-841.

Royle, J. A. and R. M. Dorazio. 2008. Hierarchical Modeling and Inference in Ecology. Academic Press.

Royle, J. A. and J. D. Nichols. 2003. Estimating Abundance from Repeated Presence-Absence Data or Point Counts. Ecology, 84:777–790.

Sillett, S. and Chandler, R.B. and Royle, J.A. and Kery, M. and Morrison, S.A. In Press. Hierarchical distance sampling models to estimate population size and habitat-specific abundance of an island endemic. Ecological Applications

Examples

## An example site-occupancy analysis

# Simulate occupancy data
set.seed(344)
nSites <- 100
nReps <- 5
covariates <- data.frame(veght=rnorm(nSites),
    habitat=factor(c(rep('A', 50), rep('B', 50))))

psipars <- c(-1, 1, -1)
ppars <- c(1, -1, 0)
X <- model.matrix(~veght+habitat, covariates) # design matrix
psi <- plogis(X %*% psipars)
p <- plogis(X %*% ppars)

y <- matrix(NA, nSites, nReps)
z <- rbinom(nSites, 1, psi)       # true occupancy state
for(i in 1:nSites) {
    y[i,] <- rbinom(nReps, 1, z[i]*p[i])
    }

# Organize data and look at it
umf <- unmarkedFrameOccu(y = y, siteCovs = covariates)
head(umf)
summary(umf)


# Fit some models
fm1 <- occu(~1 ~1, umf)
fm2 <- occu(~veght+habitat ~veght+habitat, umf)
fm3 <- occu(~veght ~veght+habitat, umf)


# Model selection
fms <- fitList(m1=fm1, m2=fm2, m3=fm3)
modSel(fms)

# Empirical Bayes estimates of the number of sites occupied
sum(bup(ranef(fm3), stat="mode"))     # Sum of posterior modes
sum(z)                                # Actual


# Model-averaged prediction and plots

# psi in each habitat type
newdata1 <- data.frame(habitat=c('A', 'B'), veght=0)
Epsi1 <- predict(fms, type="state", newdata=newdata1)
with(Epsi1, {
    plot(1:2, Predicted, xaxt="n", xlim=c(0.5, 2.5), ylim=c(0, 0.5),
        xlab="Habitat",
        ylab=expression(paste("Probability of occurrence (", psi, ")")),
        cex.lab=1.2,
        pch=16, cex=1.5)
    axis(1, 1:2, c('A', 'B'))
    arrows(1:2, Predicted-SE, 1:2, Predicted+SE, angle=90, code=3, length=0.05)
    })


# psi and p as functions of vegetation height
newdata2 <- data.frame(habitat=factor('A', levels=c('A','B')),
    veght=seq(-2, 2, length=50))
Epsi2 <- predict(fms, type="state", newdata=newdata2, appendData=TRUE)
Ep <- predict(fms, type="det", newdata=newdata2, appendData=TRUE)

op <- par(mfrow=c(2, 1), mai=c(0.9, 0.8, 0.2, 0.2))
plot(Predicted~veght, Epsi2, type="l", lwd=2, ylim=c(0,1),
    xlab="Vegetation height (standardized)",
    ylab=expression(paste("Probability of occurrence (", psi, ")")))
    lines(lower ~ veght, Epsi2, col=gray(0.7))
    lines(upper ~ veght, Epsi2, col=gray(0.7))
plot(Predicted~veght, Ep, type="l", lwd=2, ylim=c(0,1),
    xlab="Vegetation height (standardized)",
    ylab=expression(paste("Detection probability (", italic(p), ")")))
lines(lower~veght, Ep, col=gray(0.7))
lines(upper~veght, Ep, col=gray(0.7))
par(op)

Description

Methods for bracket extraction [ in Package ‘unmarked’

Usage

## S4 method for signature 'unmarkedEstimateList,ANY,ANY,ANY'
x[i, j, drop]
## S4 method for signature 'unmarkedFit,ANY,ANY,ANY'
x[i, j, drop]
## S4 method for signature 'unmarkedFrame,numeric,numeric,missing'
x[i, j]
## S4 method for signature 'unmarkedFrame,list,missing,missing'
x[i, j]
## S4 method for signature 'unmarkedMultFrame,missing,numeric,missing'
x[i, j]
## S4 method for signature 'unmarkedMultFrame,numeric,missing,missing'
x[i, j]
## S4 method for signature 'unmarkedFrameGMM,numeric,missing,missing'
x[i, j]
## S4 method for signature 'unmarkedFrameGDS,numeric,missing,missing'
x[i, j]
## S4 method for signature 'unmarkedFramePCO,numeric,missing,missing'
x[i, j]

Arguments

Methods

x = "unmarkedEstimateList", i = "ANY", j = "ANY", drop = "ANY"

Extract a unmarkedEstimate object from an unmarkedEstimateList by name (either 'det' or 'state')

x = "unmarkedFit", i = "ANY", j = "ANY", drop = "ANY"

Extract a unmarkedEstimate object from an unmarkedFit by name (either 'det' or 'state')

x = "unmarkedFrame", i = "missing", j = "numeric", drop = "missing"

Extract observations from an unmarkedFrame.

x = "unmarkedFrame", i = "numeric", j = "missing", drop = "missing"

Extract rows from an unmarkedFrame

x = "unmarkedFrame", i = "numeric", j = "numeric", drop = "missing"

Extract rows and observations from an unmarkedFrame

x = "unmarkedMultFrame", i = "missing", j = "numeric", drop = "missing"

Extract primary sampling periods from an unmarkedMultFrame

x = "unmarkedFrame", i = "list", j = "missing", drop = "missing"

List is the index of observations to subset for each site.

x = "unmarkedMultFrame", i = "numeric", j = "missing", drop = "missing"

Extract rows (sites) from an unmarkedMultFrame

x = "unmarkedGMM", i = "numeric", j = "missing", drop = "missing"

Extract rows (sites) from an unmarkedFrameGMM object

x = "unmarkedGDS", i = "numeric", j = "missing", drop = "missing"

Extract rows (sites) from an unmarkedFrameGDS object

x = "unmarkedPCO", i = "numeric", j = "missing", drop = "missing"

Extract rows (sites) from an unmarkedFramePCO object

Examples

data(mallard)
mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site,
	obsCovs = mallard.obs)
summary(mallardUMF)

mallardUMF[1:5,]
mallardUMF[,1:2]
mallardUMF[1:5, 1:2]

Description

Methods for function backTransform in Package ‘unmarked’. This converts from link-scale to original-scale

Usage

## S4 method for signature 'unmarkedFit'
backTransform(obj, type)
## S4 method for signature 'unmarkedEstimate'
backTransform(obj)

Arguments

Methods

obj = "unmarkedEstimate"

Typically done internally

obj = "unmarkedFit"

Back-transform a parameter from a fitted model. Only possible if no covariates are present. Must specify argument type as one of the values returned by names(obj).

obj = "unmarkedLinComb"

Back-transform a predicted value created by linearComb. This is done internally by predict but can be done explicitly by user.

Examples

## Not run: 

data(mallard)
mallardUMF <- unmarkedFramePCount(mallard.y, siteCovs = mallard.site, 
    obsCovs = mallard.obs)

(fm <- pcount(~ 1 ~ forest, mallardUMF))    # Fit a model
backTransform(fm, type="det")               # This works because there are no detection covariates
#backTransform(fm, type="state")             # This doesn't work because covariates are present
lc <- linearComb(fm, c(1, 0), type="state") # Estimate abundance on the log scale when forest=0
backTransform(lc)                           # Abundance on the original scale

## End(Not run)

Description

Data frames for 2 species from the breeding bird survey (BBS). Each data frame has a row for each site and columns for each sampling event. There is a point count and occurrence–designated by .bin– version for each species.

Usage

data(birds)

Format

catbird

A data frame of point count observations for the catbird.

catbird.bin

A data frame of occurrence observations for the catbird.

woodthrush

A data frame of point count observations for the wood thrush.

woodthrush.bin

A data frame of point count observations for the wood thrush.

Source

Royle J. N-mixture models for estimating population size from spatially replicated counts. Biometrics. 2004. 60(1):108–115.

Examples

data(birds)

Description

Extract coefficients

Usage

## S4 method for signature 'unmarkedFit'
coef(object, type, altNames = TRUE, fixedOnly=TRUE)
## S4 method for signature 'unmarkedEstimate'
coef(object, altNames = TRUE, fixedOnly=TRUE, ...)
## S4 method for signature 'linCombOrBackTrans'
coef(object)

Arguments

Value

A named numeric vector of parameter estimates.

Methods

object = "linCombOrBackTrans"

Object from linearComb

object = "unmarkedEstimate"

unmarkedEstimate object

object = "unmarkedFit"

Fitted model

Description

Estimate parameters of the colonization-extinction model, including covariate-dependent rates and detection process.

Usage

colext(psiformula= ~1, gammaformula =  ~ 1, epsilonformula = ~ 1,
    pformula = ~ 1, data, starts, method="BFGS", se=TRUE, ...)

Arguments

Details

This function fits the colonization-extinction model of MacKenzie et al (2003). The colonization and extinction rates can be modeled with covariates that vary yearly at each site using a logit link. These covariates are supplied by special unmarkedMultFrame yearlySiteCovs slot. These parameters are specified using the gammaformula and epsilonformula arguments. The initial probability of occupancy is modeled by covariates specified in the psiformula.

The conditional detection rate can also be modeled as a function of covariates that vary at the secondary sampling period (ie., repeat visits). These covariates are specified by the first part of the formula argument and the data is supplied via the usual obsCovs slot.

The projected and smoothed trajectories (Weir et al 2009) can be obtained from the smoothed.mean and projected.mean slots (see examples).

Value

unmarkedFitColExt object describing model fit.

References

MacKenzie, D.I. et al. (2002) Estimating Site Occupancy Rates When Detection Probabilities Are Less Than One. Ecology, 83(8), 2248-2255.

MacKenzie, D. I., J. D. Nichols, J. E. Hines, M. G. Knutson, and A. B. Franklin. 2003. Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84:2200–2207.

MacKenzie, D. I. et al. (2006) Occupancy Estimation and Modeling.Amsterdam: Academic Press.

Weir L. A., Fiske I. J., Royle J. (2009) Trends in Anuran Occupancy from Northeastern States of the North American Amphibian Monitoring Program. Herpetological Conservation and Biology. 4(3):389-402.

See Also

nonparboot, unmarkedMultFrame, and formatMult

Examples

# Fake data
R <- 4 # number of sites
J <- 3 # number of secondary sampling occasions
T <- 2 # number of primary periods

y <- matrix(c(
   1,1,0,  0,0,0,
   0,0,0,  0,0,0,
   1,1,1,  1,1,0,
   1,0,1,  0,0,1), nrow=R, ncol=J*T, byrow=TRUE)
y

site.covs <- data.frame(x1=1:4, x2=factor(c('A','B','A','B')))
site.covs

yearly.site.covs <- list(
   year = matrix(c(
      'year1', 'year2',
      'year1', 'year2',
      'year1', 'year2',
      'year1', 'year2'), nrow=R, ncol=T, byrow=TRUE)
      )
yearly.site.covs

obs.covs <- list(
   x4 = matrix(c(
      -1,0,1,  -1,1,1,
      -2,0,0,  0,0,2,
      -3,1,0,  1,1,2,
      0,0,0,   0,1,-1), nrow=R, ncol=J*T, byrow=TRUE),
   x5 = matrix(c(
      'a','b','c',  'a','b','c',
      'd','b','a',  'd','b','a',
      'a','a','c',  'd','b','a',
      'a','b','a',  'd','b','a'), nrow=R, ncol=J*T, byrow=TRUE))
obs.covs

umf <- unmarkedMultFrame(y=y, siteCovs=site.covs,
    yearlySiteCovs=yearly.site.covs, obsCovs=obs.covs,
    numPrimary=2)                  # organize data
umf                                # look at data
summary(umf)                       # summarize
fm <- colext(~1, ~1, ~1, ~1, umf)  # fit a model
fm



## Not run: 
# Real data
data(frogs)
umf <- formatMult(masspcru)
obsCovs(umf) <- scale(obsCovs(umf))

## Use 1/4 of data just for run speed in example
umf <- umf[which((1:numSites(umf)) %% 4 == 0),]

## constant transition rates
(fm <- colext(psiformula = ~ 1,
gammaformula = ~ 1,
epsilonformula = ~ 1,
pformula = ~ JulianDate + I(JulianDate^2), umf, control = list(trace=1, maxit=1e4)))

## get the trajectory estimates
smoothed(fm)
projected(fm)

# Empirical Bayes estimates of number of sites occupied in each year
re <- ranef(fm)
modes <- colSums(bup(re, stat="mode"))
plot(1:7, modes, xlab="Year", ylab="Sites occupied", ylim=c(0, 70))

## Find bootstrap standard errors for smoothed trajectory
fm <- nonparboot(fm, B = 100)  # This takes a while!
fm@smoothed.mean.bsse

## try yearly transition rates
yearlySiteCovs(umf) <- data.frame(year = factor(rep(1:7, numSites(umf))))
(fm.yearly <- colext(psiformula = ~ 1,
gammaformula = ~ year,
epsilonformula = ~ year,
pformula = ~ JulianDate + I(JulianDate^2), umf,
	control = list(trace=1, maxit=1e4)))

## End(Not run)

Description

This function computes the weight for the MPLE penalty of Moreno & Lele (2010).

Usage

computeMPLElambda(formula, data, knownOcc=numeric(0), starts,
method="BFGS",engine=c("C", "R"))

Arguments

Details

See occuPEN for details and examples.

Value

The computed lambda.

Author(s)

Rebecca A. Hutchinson

References

Moreno, M. and S. R. Lele. 2010. Improved estimation of site occupancy using penalized likelihood. Ecology 91: 341-346.

See Also

unmarked, unmarkedFrameOccu, occu, occuPEN, occuPEN_CV, nonparboot

Description

Methods for function confint in Package ‘unmarked’

Usage

## S4 method for signature 'unmarkedBackTrans'
confint(object, parm, level)
## S4 method for signature 'unmarkedEstimate'
confint(object, parm, level)
## S4 method for signature 'unmarkedLinComb'
confint(object, parm, level)
## S4 method for signature 'unmarkedFit'
confint(object, parm, level, type, method)

Arguments

Value

A vector of lower and upper confidence intervals. These are asymtotic unless method='profile' is used on unmarkedFit objects in which case they are profile likelihood intervals.

See Also

unmarkedFit-class

Description

267 1-kmsq quadrats were surveyed 3 times per year during 1999-2007.

Usage

data(crossbill)

Format

A data frame with 267 observations on the following 58 variables.

id

Plot ID

ele

Elevation

forest

Percent forest cover

surveys

a numeric vector

det991

Detection data for 1999, survey 1

det992

Detection data for 1999, survey 2

det993

Detection data for 1999, survey 3

det001

Detection data for 2000, survey 1

det002

a numeric vector

det003

a numeric vector

det011

a numeric vector

det012

a numeric vector

det013

a numeric vector

det021

a numeric vector

det022

a numeric vector

det023

a numeric vector

det031

a numeric vector

det032

a numeric vector

det033

a numeric vector

det041

a numeric vector

det042

a numeric vector

det043

a numeric vector

det051

a numeric vector

det052

a numeric vector

det053

a numeric vector

det061

a numeric vector

det062

a numeric vector

det063

Detection data for 2006, survey 3

det071

Detection data for 2007, survey 1

det072

Detection data for 2007, survey 2

det073

Detection data for 2007, survey 3

date991

Day of the season for 1999, survey 1

date992

Day of the season for 1999, survey 2

date993

Day of the season for 1999, survey 3

date001

Day of the season for 2000, survey 1

date002

a numeric vector

date003

a numeric vector

date011

a numeric vector

date012

a numeric vector

date013

a numeric vector

date021

a numeric vector

date022

a numeric vector

date023

a numeric vector

date031

a numeric vector

date032

a numeric vector

date033

a numeric vector

date041

a numeric vector

date042

a numeric vector

date043

a numeric vector

date051

a numeric vector

date052

a numeric vector

date053

a numeric vector

date061

a numeric vector

date062

a numeric vector

date063

a numeric vector

date071

a numeric vector

date072

a numeric vector

date073

Day of the season for 2007, survey 3

Source

Schmid, H. N. Zbinden, and V. Keller. 2004. Uberwachung der Bestandsentwicklung haufiger Brutvogel in der Schweiz, Swiss Ornithological Institute Sempach Switzerland

See Also

Switzerland for corresponding covariate data defined for all 1-kmsq pixels in Switzerland. Useful for making species distribution maps.

Examples

data(crossbill)
str(crossbill)

Description

Test predictive accuracy of fitted models using several cross-validation approaches. The dataset is divided by site only into folds or testing and training datasets (i.e., encounter histories within sites are never split up).

Usage

## S4 method for signature 'unmarkedFit'
crossVal(
  object, method=c("Kfold","holdout","leaveOneOut"),
  folds=10, holdoutPct=0.25, statistic=RMSE_MAE, parallel=FALSE, ncores, ...)
## S4 method for signature 'unmarkedFitList'
crossVal(
  object, method=c("Kfold","holdout","leaveOneOut"),
  folds=10, holdoutPct=0.25, statistic=RMSE_MAE, parallel=FALSE, ncores, 
  sort = c("none", "increasing", "decreasing"), ...)

Arguments

Value

unmarkedCrossVal or unmarkedCrossValList object containing calculated statistic values for each fold.

Author(s)

Ken Kellner [email protected]

See Also

fitList, unmarkedFit

Examples

## Not run: 
#Get data
data(frogs)
pferUMF <- unmarkedFrameOccu(pfer.bin)
siteCovs(pferUMF) <- data.frame(sitevar1 = rnorm(numSites(pferUMF)))    
obsCovs(pferUMF) <- data.frame(obsvar1 = rnorm(numSites(pferUMF) * obsNum(pferUMF)))

#Fit occupancy model
fm <- occu(~ obsvar1 ~ 1, pferUMF)

#k-fold cross validation with 10 folds
(kfold = crossVal(fm, method="Kfold", folds=10))

#holdout method with 25
(holdout = crossVal(fm,method='holdout', holdoutPct=0.25))

#Leave-one-out method
(leave = crossVal(fm, method='leaveOneOut'))

#Fit a second model and combine into a fitList
fm2 <- occu(~1 ~1, pferUMF)
fl <- fitList(fm2,fm)

#Cross-validation for all fits in fitList using holdout method
(cvlist <- crossVal(fl, method='holdout'))


## End(Not run)

Description

Spatially-referenced elevation, forest cover, and vegetation data for Santa Cruz Island.

Usage

data(cruz)

Format

A data frame with 2787 observations on the following 5 variables.

x

Easting (meters)

y

Northing (meters)

elevation

a numeric vector, FEET (multiply by 0.3048 to convert to meters)

forest

a numeric vector, proportion cover

chaparral

a numeric vector, proportion cover

Details

The resolution is 300x300 meters.

The Coordinate system is EPSG number 26911

NAD_1983_UTM_Zone_11N Projection: Transverse_Mercator False_Easting: 500000.000000 False_Northing: 0.000000 Central_Meridian: -117.000000 Scale_Factor: 0.999600 Latitude_Of_Origin: 0.000000 Linear Unit: Meter GCS_North_American_1983 Datum: D_North_American_1983

Source

Brian Cohen of the Nature Conservancy helped prepare the data

References

Sillett, S. and Chandler, R.B. and Royle, J.A. and Kery, M. and Morrison, S.A. In Press. Hierarchical distance sampling models to estimate population size and habitat-specific abundance of an island endemic. Ecological Applications

Examples

## Not run: 
library(lattice)
data(cruz)
str(cruz)

levelplot(elevation ~ x + y, cruz, aspect="iso",
    col.regions=terrain.colors(100))

if(require(raster)) {
elev <- rasterFromXYZ(cruz[,1:3],
     crs="+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
elev
plot(elev)
}

## End(Not run)

Description

This function converts an appropriatedly formated comma-separated values file (.csv) to a format usable by unmarked's fitting functions (see Details).

Usage

csvToUMF(filename, long=FALSE, type, species, ...)

Arguments

Details

This function provides a quick way to take a .csv file with headers named as described below and provides the data required and returns of data in the format required by the model-fitting functions in unmarked. The .csv file can be in one of 2 formats: long or wide. See the first 2 lines of the examples for what these formats look like.

The .csv file is formatted as follows:

• col 1 is site labels.

• if data is in long format, col 2 is date of observation.

• next J columns are the observations (y) - counts or 0/1's.

• next is a series of columns for the site variables (one column per variable). The column header is the variable name.

• next is a series of columns for the observation-level variables. These are in sets of J columns for each variable, e.g., var1-1 var1-2 var1-3 var2-1 var2-2 var2-3, etc. The column header of the first variable in each group must indicate the variable name.

Value

an unmarkedFrame object

Author(s)

Ian Fiske [email protected]

Examples

# examine a correctly formatted long .csv
head(read.csv(system.file("csv","frog2001pcru.csv", package="unmarked")))

# examine a correctly formatted wide .csv
head(read.csv(system.file("csv","widewt.csv", package="unmarked")))

# convert them!
dat1 <- csvToUMF(system.file("csv","frog2001pcru.csv", package="unmarked"),
                 long = TRUE, type = "unmarkedFrameOccu")
dat2 <- csvToUMF(system.file("csv","frog2001pfer.csv", package="unmarked"),
                 long = TRUE, type = "unmarkedFrameOccu")
dat3 <- csvToUMF(system.file("csv","widewt.csv", package="unmarked"),
                 long = FALSE, type = "unmarkedFrameOccu")

Description

These functions represent the currently available detection functions used for modeling line and point transect data with distsamp. Detection functions begin with "g", and density functions begin with a "d".

Usage

gxhn(x, sigma)
gxexp(x, rate)
gxhaz(x, shape, scale)

dxhn(x, sigma)
dxexp(x, rate)
dxhaz(x, shape, scale)
drhn(r, sigma)
drexp(r, rate)
drhaz(r, shape, scale)

Arguments

See Also

distsamp for example of using these for plotting detection function

Examples

# Detection probabilities at 25m for range of half-normal sigma values.
round(gxhn(25, 10:15), 2)

# Plot negative exponential distributions
plot(function(x) gxexp(x, rate=10), 0, 50, xlab="distance",
    ylab="Detection probability")
plot(function(x) gxexp(x, rate=20), 0, 50, add=TRUE, lty=2)
plot(function(x) gxexp(x, rate=30), 0, 50, add=TRUE, lty=3)

# Plot half-normal probability density functions for line- and point-transects
par(mfrow=c(2, 1))
plot(function(x) dxhn(x, 20), 0, 50, xlab="distance",
    ylab="Probability density", main="Line-transect")
plot(function(x) drhn(x, 20), 0, 50, xlab="distance",
    ylab="Probability density", main="Point-transect")

Description

Fit the hierarchical distance sampling model of Royle et al. (2004) to line or point transect data recorded in discrete distance intervals.

Usage

distsamp(formula, data, keyfun=c("halfnorm", "exp",
  "hazard", "uniform"), output=c("density", "abund"),
  unitsOut=c("ha", "kmsq"), starts, method="BFGS", se=TRUE,
  engine=c("C", "R", "TMB"), rel.tol=0.001, ...)

Arguments

Details

Unlike conventional distance sampling, which uses the 'conditional on detection' likelihood formulation, this model is based upon the unconditional likelihood and allows for modeling both abundance and detection function parameters.

The latent transect-level abundance distribution $f(N | \mathbf{\theta})$ assumed to be Poisson with mean $\lambda$ (but see gdistsamp for alternatives).

The detection process is modeled as multinomial: $y_{ij} \sim Multinomial(N_i, \pi_{ij})$, where $\pi_{ij}$ is the multinomial cell probability for transect i in distance class j. These are computed based upon a detection function $g(x | \mathbf{\sigma})$, such as the half-normal, negative exponential, or hazard rate.

Parameters $\lambda$ and $\sigma$ can be vectors affected by transect-specific covariates using the log link.

Value

unmarkedFitDS object (child class of unmarkedFit-class) describing the model fit.

Note

You cannot use obsCovs.

Author(s)

Richard Chandler [email protected]

References

Royle, J. A., D. K. Dawson, and S. Bates (2004) Modeling abundance effects in distance sampling. Ecology 85, pp. 1591-1597.

Sillett, S. and Chandler, R.B. and Royle, J.A. and Kery, M. and Morrison, S.A. In Press. Hierarchical distance sampling models to estimate population size and habitat-specific abundance of an island endemic. Ecological Applications

See Also

unmarkedFrameDS, unmarkedFit-class fitList, formatDistData, parboot, sight2perpdist, detFuns, gdistsamp, ranef. Also look at vignette("distsamp").

Examples

## Line transect examples

data(linetran)

ltUMF <- with(linetran, {
   unmarkedFrameDS(y = cbind(dc1, dc2, dc3, dc4),
   siteCovs = data.frame(Length, area, habitat),
   dist.breaks = c(0, 5, 10, 15, 20),
   tlength = linetran$Length * 1000, survey = "line", unitsIn = "m")
   })

ltUMF
summary(ltUMF)
hist(ltUMF)

# Half-normal detection function. Density output (log scale). No covariates.
(fm1 <- distsamp(~ 1 ~ 1, ltUMF))

# Some methods to use on fitted model
summary(fm1)
backTransform(fm1, type="state")                # animals / ha
exp(coef(fm1, type="state", altNames=TRUE))     # same
backTransform(fm1, type="det")                  # half-normal SD
hist(fm1, xlab="Distance (m)")	# Only works when there are no det covars
# Empirical Bayes estimates of posterior distribution for N_i
plot(ranef(fm1, K=50))

# Effective strip half-width
(eshw <- integrate(gxhn, 0, 20, sigma=10.9)$value)

# Detection probability
eshw / 20 # 20 is strip-width


# Halfnormal. Covariates affecting both density and and detection.
(fm2 <- distsamp(~area + habitat ~ habitat, ltUMF))

# Hazard-rate detection function.
(fm3 <- distsamp(~ 1 ~ 1, ltUMF, keyfun="hazard"))

# Plot detection function.
fmhz.shape <- exp(coef(fm3, type="det"))
fmhz.scale <- exp(coef(fm3, type="scale"))
plot(function(x) gxhaz(x, shape=fmhz.shape, scale=fmhz.scale), 0, 25,
	xlab="Distance (m)", ylab="Detection probability")



## Point transect examples

# Analysis of the Island Scrub-jay data.
# See Sillett et al. (In press)

data(issj)
str(issj)

jayumf <- unmarkedFrameDS(y=as.matrix(issj[,1:3]),
 siteCovs=data.frame(scale(issj[,c("elevation","forest","chaparral")])),
 dist.breaks=c(0,100,200,300), unitsIn="m", survey="point")

(fm1jay <- distsamp(~chaparral ~chaparral, jayumf))




## Not run: 

data(pointtran)

ptUMF <- with(pointtran, {
	unmarkedFrameDS(y = cbind(dc1, dc2, dc3, dc4, dc5),
	siteCovs = data.frame(area, habitat),
	dist.breaks = seq(0, 25, by=5), survey = "point", unitsIn = "m")
	})

# Half-normal.
(fmp1 <- distsamp(~ 1 ~ 1, ptUMF))
hist(fmp1, ylim=c(0, 0.07), xlab="Distance (m)")

# effective radius
sig <- exp(coef(fmp1, type="det"))
ea <- 2*pi * integrate(grhn, 0, 25, sigma=sig)$value # effective area
sqrt(ea / pi) # effective radius

# detection probability
ea / (pi*25^2)


## End(Not run)

Description

Fit the model of Dail and Madsen (2011) and Hostetler and Chandler (2015) with a distance sampling observation model (Sollmann et al. 2015).

Usage

distsampOpen(lambdaformula, gammaformula, omegaformula, pformula,
    data, keyfun=c("halfnorm", "exp", "hazard", "uniform"),
    output=c("abund", "density"), unitsOut=c("ha", "kmsq"),
    mixture=c("P", "NB", "ZIP"), K,
    dynamics=c("constant", "autoreg", "notrend", "trend", "ricker", "gompertz"),
    fix=c("none", "gamma", "omega"), immigration=FALSE, iotaformula = ~1,
    starts, method="BFGS", se=TRUE, ...)

Arguments

Details

These models generalize distance sampling models (Buckland et al. 2001) by relaxing the closure assumption (Dail and Madsen 2011, Hostetler and Chandler 2015, Sollmann et al. 2015).

The models include two or three additional parameters: gamma, either the recruitment rate (births and immigrations), the finite rate of increase, or the maximum instantaneous rate of increase; omega, either the apparent survival rate (deaths and emigrations) or the equilibrium abundance (carrying capacity); and iota, the number of immigrants per site and year. Estimates of population size at each time period can be derived from these parameters, and thus so can trend estimates. Or, trend can be estimated directly using dynamics="trend".

When immigration is set to FALSE (the default), iota is not modeled. When immigration is set to TRUE and dynamics is set to "autoreg", the model will separately estimate birth rate (gamma) and number of immigrants (iota). When immigration is set to TRUE and dynamics is set to "trend", "ricker", or "gompertz", the model will separately estimate local contributions to population growth (gamma and omega) and number of immigrants (iota).

The latent abundance distribution, $f(N | \mathbf{\theta})$ can be set as a Poisson, negative binomial, or zero-inflated Poisson random variable, depending on the setting of the mixture argument, mixture = "P", mixture = "NB", mixture = "ZIP" respectively. For the first two distributions, the mean of $N_i$ is $\lambda_i$. If $N_i \sim NB$, then an additional parameter, $\alpha$, describes dispersion (lower $\alpha$ implies higher variance). For the ZIP distribution, the mean is $\lambda_i(1-\psi)$, where psi is the zero-inflation parameter.

For "constant", "autoreg", or "notrend" dynamics, the latent abundance state following the initial sampling period arises from a Markovian process in which survivors are modeled as $S_{it} \sim Binomial(N_{it-1}, \omega_{it})$, and recruits follow $G_{it} \sim Poisson(\gamma_{it})$. Alternative population dynamics can be specified using the dynamics and immigration arguments.

$\lambda_i$, $\gamma_{it}$, and $\iota_{it}$ are modeled using the the log link. $p_{ijt}$ is modeled using the logit link. $\omega_{it}$ is either modeled using the logit link (for "constant", "autoreg", or "notrend" dynamics) or the log link (for "ricker" or "gompertz" dynamics). For "trend" dynamics, $\omega_{it}$ is not modeled.

For the distance sampling detection process, half-normal ("halfnorm"), exponential ("exp"), hazard ("hazard"), and uniform ("uniform") key functions are available.

Value

An object of class unmarkedFitDSO

Warning

This function can be extremely slow, especially if there are covariates of gamma or omega. Consider testing the timing on a small subset of the data, perhaps with se=FALSE. Finding the lowest value of K that does not affect estimates will also help with speed.

Note

When gamma or omega are modeled using year-specific covariates, the covariate data for the final year will be ignored; however, they must be supplied.

If the time gap between primary periods is not constant, an M by T matrix of integers should be supplied to unmarkedFrameDSO using the primaryPeriod argument.

Secondary sampling periods are optional, but can greatly improve the precision of the estimates.

Optimization may fail if the initial value of the intercept for the detection parameter (sigma) is too small or large relative to transect width. By default, this parameter is initialized at log(average band width). You may have to adjust this starting value.

Author(s)

Richard Chandler, Jeff Hostetler, Andy Royle, Ken Kellner

References

Buckland, S.T., Anderson, D.R., Burnham, K.P., Laake, J.L., Borchers, D.L. and Thomas, L. (2001) Introduction to Distance Sampling: Estimating Abundance of Biological Populations. Oxford University Press, Oxford, UK.

Dail, D. and L. Madsen (2011) Models for Estimating Abundance from Repeated Counts of an Open Metapopulation. Biometrics. 67: 577-587.

Hostetler, J. A. and R. B. Chandler (2015) Improved State-space Models for Inference about Spatial and Temporal Variation in Abundance from Count Data. Ecology 96: 1713-1723.

Sollmann, R., Gardner, B., Chandler, R.B., Royle, J.A. and Sillett, T.S. (2015) An open-population hierarchical distance sampling model. Ecology 96: 325-331.

See Also

distsamp, gdistsamp, unmarkedFrameDSO

Examples

## Not run: 
  
  #Generate some data 
  set.seed(123)
  lambda=4; gamma=0.5; omega=0.8; sigma=25; 
  M=100; T=10; J=4
  y <- array(NA, c(M, J, T))
  N <- matrix(NA, M, T)
  S <- G <- matrix(NA, M, T-1)
  db <- c(0, 25, 50, 75, 100)

  #Half-normal, line transect
  g <- function(x, sig) exp(-x^2/(2*sig^2))

  cp <- u <- a <- numeric(J)
  L <-  1
  a[1] <- L*db[2]
  cp[1] <- integrate(g, db[1], db[2], sig=sigma)$value
  for(j in 2:J) {
    a[j] <-  db[j+1]  - sum(a[1:j])
    cp[j] <- integrate(g, db[j], db[j+1], sig=sigma)$value
  }
  u <- a / sum(a)
  cp <- cp / a * u
  cp[j+1] <- 1-sum(cp)

  for(i in 1:M) {
    N[i,1] <- rpois(1, lambda)
    y[i,1:J,1] <- rmultinom(1, N[i,1], cp)[1:J]

    for(t in 1:(T-1)) {
        S[i,t] <- rbinom(1, N[i,t], omega)
        G[i,t] <- rpois(1, gamma)
        N[i,t+1] <- S[i,t] + G[i,t]
        y[i,1:J,t+1] <- rmultinom(1, N[i,t+1], cp)[1:J]
        }
  }
  y <- matrix(y, M)
 
  #Make a covariate
  sc <- data.frame(x1 = rnorm(M))

  umf <- unmarkedFrameDSO(y = y, siteCovs=sc, numPrimary=T, dist.breaks=db, 
                          survey="line", unitsIn="m", tlength=rep(1, M))

  (fit <- distsampOpen(~x1, ~1, ~1, ~1, data = umf, K=50, keyfun="halfnorm"))
  
  #Compare to truth
  cf <- coef(fit)
  data.frame(model=c(exp(cf[1]), cf[2], exp(cf[3]), plogis(cf[4]), exp(cf[5])), 
             truth=c(lambda, 0, gamma, omega, sigma))

  #Predict
  head(predict(fit, type='lambda'))

  #Check fit with parametric bootstrap
  pb <- parboot(fit, nsims=15)
  plot(pb)
  
  # Empirical Bayes estimates of abundance for each site / year
  re <- ranef(fit)
  plot(re, layout=c(10,5), xlim=c(-1, 10))
 
  
## End(Not run)

Description

Organize models for model selection or model-averaged prediction.

Usage

fitList(..., fits, autoNames=c("object", "formula"))

Arguments

Note

Two requirements exist to conduct AIC-based model-selection and model-averaging in unmarked. First, the data objects (ie, unmarkedFrames) must be identical among fitted models. Second, the response matrix must be identical among fitted models after missing values have been removed. This means that if a response value was removed in one model due to missingness, it needs to be removed from all models.

Author(s)

Richard Chandler [email protected]

Examples

data(linetran)
(dbreaksLine <- c(0, 5, 10, 15, 20)) 
lengths <- linetran$Length * 1000

ltUMF <- with(linetran, {
	unmarkedFrameDS(y = cbind(dc1, dc2, dc3, dc4), 
	siteCovs = data.frame(Length, area, habitat), dist.breaks = dbreaksLine,
	tlength = lengths, survey = "line", unitsIn = "m")
	})

fm1 <- distsamp(~ 1 ~1, ltUMF)
fm2 <- distsamp(~ area ~1, ltUMF)
fm3 <- distsamp( ~ 1 ~area, ltUMF)

## Two methods of creating an unmarkedFitList using fitList()

# Method 1
fmList <- fitList(Null=fm1, .area=fm2, area.=fm3)

# Method 2. Note that the arugment name "fits" must be included in call.
models <- list(Null=fm1, .area=fm2, area.=fm3)
fmList <- fitList(fits = models)

# Extract coefficients and standard errors
coef(fmList)
SE(fmList)

# Model-averaged prediction
predict(fmList, type="state")

# Model selection
modSel(fmList, nullmod="Null")

Description

Extracted fitted values from a fitted model.

Usage

## S4 method for signature 'unmarkedFit'
fitted(object, na.rm = FALSE)
## S4 method for signature 'unmarkedFitColExt'
fitted(object, na.rm = FALSE)
## S4 method for signature 'unmarkedFitOccu'
fitted(object, na.rm = FALSE)
## S4 method for signature 'unmarkedFitOccuRN'
fitted(object, K, na.rm = FALSE)
## S4 method for signature 'unmarkedFitPCount'
fitted(object, K, na.rm = FALSE)
## S4 method for signature 'unmarkedFitDS'
fitted(object, na.rm = FALSE)

Arguments

Value

Returns a matrix of expected values

Methods

object = "unmarkedFit"

A fitted model

object = "unmarkedFitColExt"

A model fit by colext

object = "unmarkedFitOccu"

A model fit by occu

object = "unmarkedFitOccuRN"

A model fit by occuRN

object = "unmarkedFitPCount"

A model fit by pcount

object = "unmarkedFitDS"

A model fit by distsamp

Description

Convert individual-level distance data to the transect-level format required by distsamp or gdistsamp

Usage

formatDistData(distData, distCol, transectNameCol, dist.breaks,
                      occasionCol, effortMatrix)

Arguments

Details

This function creates a site (M) by distance interval (J) response matrix from a data.frame containing the detection distances for each individual and the transect names. Alternatively, if each transect was surveyed T times, the resulting matrix is M x JT, which is the format required by gdistsamp, seeunmarkedFrameGDS.

Value

An M x J or M x JT matrix containing the binned distance data. Transect names will become rownames and colnames will describe the distance intervals.

Note

It is important that the factor containing transect names includes levels for all the transects surveyed, not just those with >=1 detection. Likewise, if transects were visited more than once, the factor containing the occasion numbers should include levels for all occasions. See the example for how to add levels to a factor.

See Also

distsamp, unmarkedFrame

Examples

# Create a data.frame containing distances of animals detected
# along 4 transects.
dat <- data.frame(transect=gl(4,5, labels=letters[1:4]),
                  distance=rpois(20, 10))
dat

# Look at your transect names.
levels(dat$transect)

# Suppose that you also surveyed a transect named "e" where no animals were
# detected. You must add it to the levels of dat$transect
levels(dat$transect) <- c(levels(dat$transect), "e")
levels(dat$transect)

# Distance cut points defining distance intervals
cp <- c(0, 8, 10, 12, 14, 18)

# Create formated response matrix
yDat <- formatDistData(dat, "distance", "transect", cp)
yDat

# Now you could merge yDat with transect-level covariates and
# then use unmarkedFrameDS to prepare data for distsamp


## Example for data from multiple occasions

dat2 <- data.frame(distance=1:100, site=gl(5, 20),
                   visit=factor(rep(1:4, each=5)))
cutpt <- seq(0, 100, by=25)
y2 <- formatDistData(dat2, "distance", "site", cutpt, "visit")
umf <- unmarkedFrameGDS(y=y2, numPrimary=4, survey="point",
                        dist.breaks=cutpt, unitsIn="m")
 ## Example for datda from multiple occasions with effortMatrix
 
 dat3 <-  data.frame(distance=1:100, site=gl(5, 20), visit=factor(rep(1:4, each=5)))
 cutpt <- seq(0, 100, by=25)
 
 effortMatrix <- matrix(ncol=4, nrow=5, rbinom(20,1,0.8))
 
 y3 <- formatDistData(dat2, "distance", "site", cutpt, "visit", effortMatrix)

Description

This convenience function converts multi-year data in long format to unmarkedMultFrame Object. See Details for more information.

Usage

formatMult(df.in)

Arguments

Details

df.in is a data frame with columns formatted as follows:

Column 1 = year number
Column 2 = site name or number
Column 3 = julian date or chronological sample number during year
Column 4 = observations (y)
Column 5 – Final Column = covariates

Note that if the data is already in wide format, it may be easier to create an unmarkedMultFrame object directly with a call to unmarkedMultFrame.

Value

unmarkedMultFrame object

Description

Convert a data.frame between wide and long formats.

Usage

formatWide(dfin, sep = ".", obsToY, type, ...)
formatLong(dfin, species = NULL, type, ...)

Arguments

Details

Note that not all possible unmarkedFrame* classes have been tested with these functions. Multinomial data sets (e.g., removal, double-observer, capture-recapture) are almost certainly easier to enter directly to the constructor function and are not supported by formatLong or formatWide.

In order for these functions to work, the columns of dfin need to be in the correct order. formatLong requires that the columns are in the following scheme:

1. site name or number.

2. date or observation number.

3. response variable (detections, counts, etc).

4. The remaining columns are observation-level covariates.

formatWide requires particular names for the columns. The column order for formatWide is

1. (optional) site name, named “site”.

2. response, named “y.1”, “y.2”, ..., “y.J”.

3. columns of site-level covariates, each with a relevant name per column.

4. groups of columns of observation-level covariates, each group having the name form “someObsCov.1”, “someObsCov.2”, ..., “someObsCov.J”.

Value

A data.frame

See Also

csvToUMF

Description

frogs contains NAAMP data for Pseudacris feriarum (pfer) and Pseudacris crucifer (pcru) in 2001.

Usage

data(frogs)

Format

pcru.y

matrix of observed calling indices for pcru

pcru.bin

matrix of detections for pcru

pcru.data

array of covariates measured at the observation-level for pcru

pfer.y

matrix of observed calling indices for pfer

pfer.bin

matrix of detections for pfer

pfer.data

array of covariates measured at the observation-level for pfer

Details

The rows of pcru.y, pcru.bin, pfer.y, and pfer.bin correspond to sites and columns correspond to visits to each site. The first 2 dimensions of pfer.data and pcru.data are matrices of covariates that correspond to the observation matrices (sites $\times$ observation), with the 3rd dimension corresponding to separate covariates.

Source

https://www.pwrc.usgs.gov/naamp/

References

Mossman MJ, Weir LA. North American Amphibian Monitoring Program (NAAMP). Amphibian Declines: the conservation status of United States species. University of California Press, Berkeley, California, USA. 2005:307-313.

Examples

data(frogs)
str(pcru.data)

Description

Fit the model of Amundson et al. (2014) to point count datasets containing both distance and time of observation data. The Amundson et al. (2014) model is extended to account for temporary emigration by estimating an additional availability probability if multiple counts at a site are available. Abundance can be modeled as a Poisson, negative binomial, or Zero-inflated Poisson. Multiple distance sampling key functions are also available.

Usage

gdistremoval(lambdaformula=~1, phiformula=~1, removalformula=~1,
  distanceformula=~1, data, keyfun=c("halfnorm", "exp", "hazard", "uniform"),
  output=c("abund", "density"), unitsOut=c("ha", "kmsq"), mixture=c('P', 'NB', 'ZIP'), 
  K, starts, method = "BFGS", se = TRUE, engine=c("C","TMB"), threads=1, ...)

Arguments

Value

An object of class unmarkedFitGDR

Author(s)

Ken Kellner [email protected]

References

Amundson, C.L., Royle, J.A. and Handel, C.M., 2014. A hierarchical model combining distance sampling and time removal to estimate detection probability during avian point counts. The Auk 131: 476-494.

See Also

unmarkedFrameGDR, gdistsamp, gmultmix

Description

Extends the distance sampling model of Royle et al. (2004) to estimate the probability of being available for detection. Also allows abundance to be modeled using the negative binomial and zero-inflated Poisson distributions.

Usage

gdistsamp(lambdaformula, phiformula, pformula, data, keyfun =
c("halfnorm", "exp", "hazard", "uniform"), output = c("abund",
"density"), unitsOut = c("ha", "kmsq"), mixture = c("P", "NB", "ZIP"), K,
starts, method = "BFGS", se = TRUE, engine=c("C","R"), rel.tol=1e-4, threads=1, ...)

Arguments

Details

Extends the model of Royle et al. (2004) by estimating the probability of being available for detection $\phi$. To estimate this additional parameter, replicate distance sampling data must be collected at each transect. Thus the data are collected at i = 1, 2, ..., R transects on t = 1, 2, ..., T occassions. As with the model of Royle et al. (2004), the detections must be binned into distance classes. These data must be formatted in a matrix with R rows, and JT columns where J is the number of distance classses. See unmarkedFrameGDS for more information about data formatting.

The definition of availability depends on the context. The model is

$M_i \sim \text{Pois}(\lambda)$

$N_{i,t} \sim \text{Bin}(M_i, \phi)$

$y_{i,1,t}, \dots, y_{i,J,t} \sim \text{Multinomial}(N_{i,t}, \pi_{i,1,t}, \dots, \pi_{i,J,t})$

If there is no movement, then $M_i$ is local abundance, and $N_{i,t}$ is the number of individuals that are available to be detected. In this case, $\phi=g_0$. Animals might be missed on the transect line because they are difficult to see or detected. This relaxes the assumption of conventional distance sampling that $g_0=1$.

However, when there is movement in the form of temporary emigration, local abundance is $N_{i,t}$; it's the fraction of $M_i$ that are on the plot at time t. In this case, $\phi$ is the temporary emigration parameter, and we need to assume that $g_0=1$ in order to interpret $N_{i,t}$ as local abundance. See Chandler et al. (2011) for an analysis of the model under this form of temporary emigration.

If there is movement and $g_0<1$ then it isn't possible to estimate local abundance at time t. In this case, $M_i$ would be the total number of individuals that ever use plot i (the super-population), and $N_{i,t}$ would be the number available to be detected at time t. Since a fraction of the unavailable individuals could be off the plot, and another fraction could be on the plot, it isn't possible to infer local abundance and density during occasion t.

Value

An object of class unmarkedFitGDS.

Note

If you aren't interested in estimating $\phi$, but you want to use the negative binomial or ZIP distributions, set numPrimary=1 when formatting the data.

Note

You cannot use obsCovs, but you can use yearlySiteCovs (a confusing name since this model isn't for multi-year data. It's just a hold-over from the colext methods of formatting data upon which it is based.)

Author(s)

Richard Chandler [email protected]

References

Royle, J. A., D. K. Dawson, and S. Bates. 2004. Modeling abundance effects in distance sampling. Ecology 85:1591-1597.

Chandler, R. B, J. A. Royle, and D. I. King. 2011. Inference about density and temporary emigration in unmarked populations. Ecology 92:1429–1435.

See Also

distsamp

Examples

# Simulate some line-transect data

set.seed(36837)

R <- 50 # number of transects
T <- 5  # number of replicates
strip.width <- 50
transect.length <- 100
breaks <- seq(0, 50, by=10)

lambda <- 5 # Abundance
phi <- 0.6  # Availability
sigma <- 30 # Half-normal shape parameter

J <- length(breaks)-1
y <- array(0, c(R, J, T))
for(i in 1:R) {
    M <- rpois(1, lambda) # Individuals within the 1-ha strip
    for(t in 1:T) {
        # Distances from point
        d <- runif(M, 0, strip.width)
        # Detection process
        if(length(d)) {
            cp <- phi*exp(-d^2 / (2 * sigma^2)) # half-normal w/ g(0)<1
            d <- d[rbinom(length(d), 1, cp) == 1]
            y[i,,t] <- table(cut(d, breaks, include.lowest=TRUE))
            }
        }
    }
y <- matrix(y, nrow=R) # convert array to matrix

# Organize data
umf <- unmarkedFrameGDS(y = y, survey="line", unitsIn="m",
    dist.breaks=breaks, tlength=rep(transect.length, R), numPrimary=T)
summary(umf)


# Fit the model
m1 <- gdistsamp(~1, ~1, ~1, umf, output="density", K=50)

summary(m1)


backTransform(m1, type="lambda")
backTransform(m1, type="phi")
backTransform(m1, type="det")

## Not run: 
# Empirical Bayes estimates of abundance at each site
re <- ranef(m1)
plot(re, layout=c(10,5), xlim=c(-1, 20))

## End(Not run)

Description

Methods for function getB in Package ‘unmarked’. These methods return a matrix of probabilities detections were certain for occupancy models that account for false positives.

Description

Methods for function getFP in Package ‘unmarked’. These methods return a matrix of false positive detection probabilities.

Description

Methods for function getP in Package ‘unmarked’. These methods return a matrix of the back-transformed detection parameter ($p$ the detection probability or $\lambda$ the detection rate, depending on the model). The matrix is of dimension MxJ, with M the number of sites and J the number of sampling periods; or of dimension MxJT for models with multiple primary periods T.

Methods

signature(object = "unmarkedFit")

A fitted model object

signature(object = "unmarkedFitDS")

A fitted model object

signature(object = "unmarkedFitMPois")

A fitted model object

signature(object = "unmarkedFitGMM")

A fitted model object

signature(object = "unmarkedFitOccuCOP")

With unmarkedFitOccuCOP the object of a model fitted with occuCOP. Returns a matrix of $\lambda$ the detection rate.

Description

Multinomial calling index data.

Usage

data(gf)

Format

A list with 2 components

gf.data

220 x 3 matrix of count indices

gf.obs

list of covariates

References

Royle, J. Andrew, and William A. Link. 2005. A General Class of Multinomial Mixture Models for Anuran Calling Survey Data. Ecology 86, no. 9: 2505–2512.

Examples

data(gf)
str(gf.data)
str(gf.obs)

Description

A three level hierarchical model for designs involving repeated counts that yield multinomial outcomes. Possible data collection methods include repeated removal sampling and double observer sampling. The three model parameters are abundance, availability, and detection probability.

Usage

gmultmix(lambdaformula, phiformula, pformula, data, mixture = c("P", "NB", "ZIP"), K, 
         starts, method = "BFGS", se = TRUE, engine=c("C","R"), threads=1, ...)

Arguments

Details

The latent transect-level super-population abundance distribution $f(M | \mathbf{\theta})$ can be set as a Poisson, negative binomial, or zero-inflated Poisson random variable, depending on the setting of the mixture argument. mixture = "P", mixture = "NB", and mixture = "ZIP" select the Poisson, negative binomial, and zero-inflated Poisson distributions respectively. The mean of $M_i$ is $\lambda_i$. If $M_i \sim NB$, then an additional parameter, $\alpha$, describes dispersion (lower $\alpha$ implies higher variance). If $M_i \sim ZIP$, then an additional zero-inflation parameter $\psi$ is estimated.

The number of individuals available for detection at time j is a modeled as binomial: $N_{ij} \sim Binomial(M_i, \mathbf{\phi_{ij}})$.

The detection process is modeled as multinomial: $\mathbf{y_{it}} \sim Multinomial(N_{it}, \pi_{it})$, where $\pi_{ijt}$ is the multinomial cell probability for plot i at time t on occasion j.

Cell probabilities are computed via a user-defined function related to the sampling design. Alternatively, the default functions removalPiFun or doublePiFun can be used for equal-interval removal sampling or double observer sampling. Note that the function for computing cell probabilites is specified when setting up the data using unmarkedFrameGMM.

Parameters $\lambda$, $\phi$ and $p$ can be modeled as linear functions of covariates using the log, logit and logit links respectively.

Value

An object of class unmarkedFitGMM.

Note

In the case where availability for detection is due to random temporary emigration, population density at time j, D(i,j), can be estimated by N(i,j)/plotArea.

This model is also applicable to sampling designs in which the local population size is closed during the J repeated counts, and availability is related to factors such as the probability of vocalizing. In this case, density can be estimated by M(i)/plotArea.

If availability is a function of both temporary emigration and other processess such as song rate, then density cannot be directly estimated, but inference about the super-population size, M(i), is possible.

Three types of covariates can be supplied, site-level, site-by-year-level, and observation-level. These must be formatted correctly when organizing the data with unmarkedFrameGPC

Author(s)

Richard Chandler [email protected] and Andy Royle

References

Royle, J. A. (2004) Generalized estimators of avian abundance from count survey data. Animal Biodiversity and Conservation 27, pp. 375–386.

Chandler, R. B., J. A. Royle, and D. I. King. 2011. Inference about density and temporary emigration in unmarked populations. Ecology 92:1429-1435.

See Also

unmarkedFrameGMM for setting up the data and metadata. multinomPois for surveys where no secondary sampling periods were used. Example functions to calculate multinomial cell probabilities are described piFuns

Examples

# Simulate data using the multinomial-Poisson model with a
# repeated constant-interval removal design.

n <- 100  # number of sites
T <- 4    # number of primary periods
J <- 3    # number of secondary periods

lam <- 3
phi <- 0.5
p <- 0.3

#set.seed(26)
y <- array(NA, c(n, T, J))
M <- rpois(n, lam)          # Local population size
N <- matrix(NA, n, T)       # Individuals available for detection

for(i in 1:n) {
    N[i,] <- rbinom(T, M[i], phi)
    y[i,,1] <- rbinom(T, N[i,], p)    # Observe some
    Nleft1 <- N[i,] - y[i,,1]         # Remove them
    y[i,,2] <- rbinom(T, Nleft1, p)   # ...
    Nleft2 <- Nleft1 - y[i,,2]
    y[i,,3] <- rbinom(T, Nleft2, p)
    }

y.ijt <- cbind(y[,1,], y[,2,], y[,3,], y[,4,])


umf1 <- unmarkedFrameGMM(y=y.ijt, numPrimary=T, type="removal")

(m1 <- gmultmix(~1, ~1, ~1, data=umf1, K=30))

backTransform(m1, type="lambda")        # Individuals per plot
backTransform(m1, type="phi")           # Probability of being avilable
(p <- backTransform(m1, type="det"))    # Probability of detection
p <- coef(p)

# Multinomial cell probabilities under removal design
c(p, (1-p) * p, (1-p)^2 * p)

# Or more generally:
head(getP(m1))

# Empirical Bayes estimates of super-population size
re <- ranef(m1)
plot(re, layout=c(5,5), xlim=c(-1,20), subset=site%in%1:25)

Description

Fit multi-scale occupancy models as described in Nichols et al. (2008) to repeated presence-absence data collected using the robust design. This model allows for inference about occupancy, availability, and detection probability.

Usage

goccu(psiformula, phiformula, pformula, data, linkPsi = c("logit", "cloglog"),
      starts, method = "BFGS", se = TRUE, ...)

Arguments

Details

Primary periods could represent spatial or temporal sampling replicates. For example, you could have several spatial sub-units within each site, where each sub-unit was then sampled repeatedly. This is a frequent design for eDNA studies. Or, you could have multiple primary periods of sampling at each site (conducted at different times within a season), each of which contains several secondary sampling periods. In both cases the robust design structure can be used to estimate an availability probability in addition to detection probability. See Kery and Royle (2015) 10.10 for more details.

Value

An object of class unmarkedFitGOccu

Author(s)

Ken Kellner [email protected]

References

Kery, M., & Royle, J. A. (2015). Applied hierarchical modeling in ecology: Volume 1: Prelude and static models. Elsevier Science.

Nichols, J. D., Bailey, L. L., O'Connell Jr, A. F., Talancy, N. W., Campbell Grant, E. H., Gilbert, A. T., Annand E. M., Husband, T. P., & Hines, J. E. (2008). Multi-scale occupancy estimation and modelling using multiple detection methods. Journal of Applied Ecology, 45(5), 1321-1329.

See Also

occu, colext, unmarkedMultFrame, unmarkedFrameGOccu

Examples

set.seed(123)
M <- 100
T <- 5
J <- 4

psi <- 0.5
phi <- 0.3
p <- 0.4

z <- rbinom(M, 1, psi)
zmat <- matrix(z, nrow=M, ncol=T)

zz <- rbinom(M*T, 1, zmat*phi)
zz <- matrix(zz, nrow=M, ncol=T)

zzmat <- zz[,rep(1:T, each=J)]
y <- rbinom(M*T*J, 1, zzmat*p)
y <- matrix(y, M, J*T)
umf <- unmarkedMultFrame(y=y, numPrimary=T)

## Not run: 
  mod <- goccu(psiformula = ~1, phiformula = ~1, pformula = ~1, umf)
  plogis(coef(mod))

## End(Not run)

Description

Fit the model of Chandler et al. (2011) to repeated count data collected using the robust design. This model allows for inference about population size, availability, and detection probability.

Usage

gpcount(lambdaformula, phiformula, pformula, data,
mixture = c("P", "NB", "ZIP"), K, starts, method = "BFGS", se = TRUE,
engine = c("C", "R"), threads=1, ...)

Arguments

Details

The latent transect-level super-population abundance distribution $f(M | \mathbf{\theta})$ can be set as either a Poisson, negative binomial, or zero-inflated Poisson random variable, depending on the setting of the mixture argument. The expected value of $M_i$ is $\lambda_i$. If $M_i \sim NB$, then an additional parameter, $\alpha$, describes dispersion (lower $\alpha$ implies higher variance). If $M_i \sim ZIP$, then an additional zero-inflation parameter $\psi$ is estimated.

The number of individuals available for detection at time j is a modeled as binomial: $N_{ij} \sim Binomial(M_i, \mathbf{\phi_{ij}})$.

The detection process is also modeled as binomial: $y_{ikj} \sim Binomial(N_{ij}, p_{ikj})$.

Parameters $\lambda$, $\phi$ and $p$ can be modeled as linear functions of covariates using the log, logit and logit links respectively.

Value

An object of class unmarkedFitGPC

Note

In the case where availability for detection is due to random temporary emigration, population density at time j, D(i,j), can be estimated by N(i,j)/plotArea.

This model is also applicable to sampling designs in which the local population size is closed during the J repeated counts, and availability is related to factors such as the probability of vocalizing. In this case, density can be estimated by M(i)/plotArea.

If availability is a function of both temporary emigration and other processess such as song rate, then density cannot be directly estimated, but inference about the super-population size, M(i), is possible.

Three types of covariates can be supplied, site-level, site-by-year-level, and observation-level. These must be formatted correctly when organizing the data with unmarkedFrameGPC

Author(s)

Richard Chandler [email protected]

References

Royle, J. A. 2004. N-Mixture models for estimating population size from spatially replicated counts. Biometrics 60:108–105.

Chandler, R. B., J. A. Royle, and D. I. King. 2011. Inference about density and temporary emigration in unmarked populations. Ecology 92:1429-1435.

See Also

gmultmix, gdistsamp, unmarkedFrameGPC

Examples

set.seed(54)

nSites <- 20
nVisits <- 4
nReps <- 3

lambda <- 5
phi <- 0.7
p <- 0.5

M <- rpois(nSites, lambda) # super-population size

N <- matrix(NA, nSites, nVisits)
y <- array(NA, c(nSites, nReps, nVisits))
for(i in 1:nVisits) {
    N[,i] <- rbinom(nSites, M, phi) # population available during vist j
}
colMeans(N)

for(i in 1:nSites) {
    for(j in 1:nVisits) {
        y[i,,j] <- rbinom(nReps, N[i,j], p)
    }
}

ym <- matrix(y, nSites)
ym[1,] <- NA
ym[2, 1:nReps] <- NA
ym[3, (nReps+1):(nReps+nReps)] <- NA
umf <- unmarkedFrameGPC(y=ym, numPrimary=nVisits)

## Not run: 
fmu <- gpcount(~1, ~1, ~1, umf, K=40, control=list(trace=TRUE, REPORT=1))

backTransform(fmu, type="lambda")
backTransform(fmu, type="phi")
backTransform(fmu, type="det")

## End(Not run)

Description

Model abundance using a combination of distance sampling data (DS) and other similar data types, including simple point counts (PC) and occupancy/detection-nondetection (OC/DND) data.

Usage

IDS(lambdaformula = ~1, detformulaDS = ~1, detformulaPC = NULL, detformulaOC = NULL,
    dataDS, dataPC = NULL, dataOC = NULL, availformula = NULL,
    durationDS = NULL, durationPC = NULL, durationOC = NULL, keyfun = "halfnorm",
    maxDistPC, maxDistOC, K = 100, unitsOut = "ha", 
    starts = NULL, method = "BFGS", ...)

Arguments

Details

This function facilitates a combined analysis of distance sampling data (DS) with other similar data types, including simple point counts (PC) and occupancy/detection-nondetection (OC/DND) data. The combined approach capitalizes on the strengths and minimizes the weaknesses of each type. The PC and OC/DND data are viewed as latent distance sampling surveys with an underlying abundance model shared by all data types. All analyses must include some distance sampling data, but can include either PC or DND data or both. If surveys are of variable duration, it is also possible to estimate availability.

Input data must be provided as a series of separate unmarkedFrames: unmarkedFrameDS for the distance sampling data, unmarkedFramePCount for the point count data, and unmarkedFrameOccu for OC/DND data. See the help files for these objects for guidance on how to organize the data.

Value

An object of class unmarkedFitIDS

Note

Simulations indicated estimates of availability were very unreliable when including detection/non-detection data, so the function will not allow you to use DND data and estimate availability at the same time. In general estimation of availability can be difficult; use simulations to see how well it works for your specific situation.

Author(s)

Ken Kellner [email protected]

References

Kery M, Royle JA, Hallman T, Robinson WD, Strebel N, Kellner KF. 2024. Integrated distance sampling models for simple point counts. Ecology.

See Also

distsamp

Examples

## Not run: 

# Simulate data based on a real dataset

# Formulas for each model
formulas <- list(lam=~elev, ds=~1, phi=~1)

# Sample sizes
design <- list(Mds=2912, J=6, Mpc=506)

# Model parameters
coefs <- list(lam = c(intercept=3, elev=-0.5),
              ds = c(intercept=-2.5),
              phi = c(intercept=-1.3))

# Survey durations
durs <- list(ds = rep(5, design$Mds), pc=runif(design$Mpc, 3, 30))

set.seed(456)
sim_umf <- simulate("IDS", # name of model we are simulating for
                    nsim=1, # number of replicates
                    formulas=formulas, 
                    coefs=coefs,
                    design=design,
                    # arguments used by unmarkedFrameDS
                    dist.breaks = seq(0, 0.30, length.out=7),
                    unitsIn="km", 
                    # arguments used by IDS
                    # could also have e.g. keyfun here
                    durationDS=durs$ds, durationPC=durs$pc, durationOC=durs$oc,
                    maxDistPC=0.5, maxDistOC=0.5,
                    unitsOut="kmsq")

# Look at the results
lapply(sim_umf, head)

# Fit a model
(mod_sim <- IDS(lambdaformula = ~elev, detformulaDS = ~1,
                dataDS=sim_umf$ds, dataPC=sim_umf$pc,
                availformula = ~1, durationDS=durs$ds, durationPC=durs$pc,
                maxDistPC=0.5,
                unitsOut="kmsq"))

# Compare with known parameter values
# Note:  this is an unusually good estimate of availability
# It is hard to estimate in most cases
cbind(truth=unlist(coefs), est=coef(mod_sim))

# Predict density at each distance sampling site
head(predict(mod_sim, 'lam'))


## End(Not run)

Description

This function uses an ad-hoc averaging approach to impute missing entries in obsCovs. The missing entry is replaced by an average of the average for the site and the average for the visit number.

Usage

imputeMissing(umf, whichCovs = seq(length=ncol(obsCovs(umf))))

Arguments

Value

A version of umf that has the requested obsCovs imputed.

Author(s)

Ian Fiske

Examples

data(frogs)
pcru.obscovs <- data.frame(MinAfterSunset=as.vector(t(pcru.data[,,1])),
     Wind=as.vector(t(pcru.data[,,2])),
     Sky=as.vector(t(pcru.data[,,3])),
     Temperature=as.vector(t(pcru.data[,,4])))
pcruUMF <- unmarkedFrameOccu(y = pcru.bin, obsCovs = pcru.obscovs)
pcruUMF.i1 <- imputeMissing(pcruUMF)
pcruUMF.i2 <- imputeMissing(pcruUMF, whichCovs = 2)

Description

Data were collected at 307 survey locations ("point transects") on Santa Cruz Island, California during the Fall of 2008. The distance data are binned into 3 distance intervals [0-100], (100-200], and (200-300]. The coordinates of the survey locations as well as 3 habitat covariates are also included.

Usage

data(issj)

Format

A data frame with 307 observations on the following 8 variables.

issj[0-100]

Number of individuals detected within 100m

issj(100-200]

Detections in the interval (100-200m]

issj(200-300]

Detections in the interval (200-300m]

x

Easting (meters)

y

Northing (meters)

elevation

Elevation in meters

forest

Forest cover

chaparral

Chaparral cover

References

Sillett, S. and Chandler, R.B. and Royle, J.A. and Kery, M. and Morrison, S.A. In Press. Hierarchical distance sampling models to estimate population size and habitat-specific abundance of an island endemic. Ecological Applications

See Also

Island-wide covariates are also available cruz

Examples

data(issj)
str(issj)
head(issj)

umf <- unmarkedFrameDS(y=as.matrix(issj[,1:3]), siteCovs=issj[,6:8],
    dist.breaks=c(0,100,200,300), unitsIn="m", survey="point")
summary(umf)

Description

The Swiss breeding bird survey ("Monitoring Haufige Brutvogel" MHB) has monitored the populations of 150 common species since 1999. The MHB sample consists of 267 1-km squares that are laid out as a grid across Switzerland. Fieldwork is conducted by about 200 skilled birdwatchers, most of them volunteers. Avian populations are monitored using a simplified territory mapping protocol, where each square is surveyed up to three times during the breeding season (only twice above the tree line). Surveys are conducted along a transect that does not change over the years.

The list jay has the data for European Jay territories for 238 sites surveyed in 2002.

Usage

data("jay")

Format

jay is a list with 3 elements:

caphist

a data frame with rows for 238 sites and columns for each of the observable detection histories. For the sites visited 3 times, these are "100", "010", "001", "110", "101", "011", "111". Sites visited twice have "10x", "01x", "11x".

Each row gives the number of territories with the corresponding detection history, with NA for the detection histories not applicable: sites visited 3 times have NAs in the last 3 columns while those visited twice have NAs in the first 7 columns.

sitescovs

a data frame with rows for 238 sites, and the following columns:

1. elev : the mean elevation of the quadrat, m.

2. length : the length of the route walked in the quadrat, km.

3. forest : percentage forest cover.

covinfo

a data frame with rows for 238 sites, and the following columns:

1. x, y : the coordinates of the site.

2. date1, date2, date3 : the Julian date of the visit, with 1 April = 1. Sites visited twice have NA in the 3rd column.

3. dur1, dur2, dur3 : the duration of the survey, mins. For 10 visits the duration is not available, so there are additional NAs in these columns.

Note

In previous versions, jay had additional information not required for the analysis, and a data frame with essentially the same information as the Switzerland data set.

Source

Swiss Ornithological Institute

References

Royle, J.A., Kery, M., Gauthier, R., Schmid, H. (2007) Hierarchical spatial models of abundance and occurrence from imperfect survey data. Ecological Monographs, 77, 465-481.

Kery & Royle (2016) Applied Hierarachical Modeling in Ecology Section 7.9

Examples

data(jay)
str(jay)

# Carry out a simple analysis, without covariates:
# Create a customised piFun (see ?piFun for details)
crPiFun <- function(p) {
   p1 <- p[,1] # Extract the columns of the p matrix, one for 
   p2 <- p[,2] #   each of J = 3 sample occasions
   p3 <- p[,3]
   cbind(      # define multinomial cell probabilities:
      "100" = p1 * (1-p2) * (1-p3),
      "010" = (1-p1) * p2 * (1-p3),
      "001" = (1-p1) * (1-p2) * p3,
      "110" = p1 * p2 * (1-p3),
      "101" = p1 * (1-p2) * p3,
      "011" = (1-p1) * p2 * p3,
      "111" = p1 * p2 * p3,
      "10x" = p1*(1-p2),
      "01x" = (1-p1)*p2,
      "11x" = p1*p2)
}
# Build the unmarkedFrame object
mhb.umf <- unmarkedFrameMPois(y=as.matrix(jay$caphist),
  obsToY=matrix(1, 3, 10), piFun="crPiFun")
# Fit a model
( fm1 <- multinomPois(~1 ~1, mhb.umf) )