Fit the MacKenzie et al. (2002) Occupancy Model with the penalized likelihood methods of Hutchinson et al. (2015)

This function fits the occupancy model of MacKenzie et al (2002) with the penalized methods of Hutchinson et al (2015).

occuPEN(formula, data, knownOcc=numeric(0), starts, method="BFGS",
    engine=c("C", "R"), lambda=0, pen.type = c("Bayes","Ridge","MPLE"), ...)

Arguments

formula

Double right-hand side formula describing covariates of detection and occupancy in that order.

data

An unmarkedFrameOccu object

knownOcc

Vector of sites that are known to be occupied. These should be supplied as row numbers of the y matrix, eg, c(3,8) if sites 3 and 8 were known to be occupied a priori.

starts

Vector of parameter starting values.

method

Optimization method used by optim.

engine

Either "C" or "R" to use fast C++ code or native R code during the optimization.

lambda

Penalty weight parameter.

pen.type

Which form of penalty to use.

...

Additional arguments to optim, such as lower and upper bounds

Details

See unmarkedFrame and unmarkedFrameOccu for a description of how to supply data to the data argument.

occuPEN fits the standard occupancy model based on zero-inflated binomial models (MacKenzie et al. 2006, Royle and Dorazio 2008) using the penalized likelihood methods described in Hutchinson et al. (2015). See occu for model details. occuPEN returns parameter estimates that maximize a penalized likelihood in which the penalty is specified by the pen.type argument. The penalty function is weighted by lambda.

The MPLE method includes an equation for computing lambda (Moreno & Lele, 2010). If the value supplied does not equal match the one computed with this equation, the supplied value is used anyway (with a warning).

Value

unmarkedFitOccuPEN object describing the model fit.

References

Hutchinson, R. A., J. V. Valente, S. C. Emerson, M. G. Betts, and T. G. Dietterich. 2015. Penalized Likelihood Methods Improve Parameter Estimates in Occupancy Models. Methods in Ecology and Evolution. DOI: 10.1111/2041-210X.12368

MacKenzie, D. I., J. D. Nichols, G. B. Lachman, S. Droege, J. Andrew Royle, and C. A. Langtimm. 2002. Estimating Site Occupancy Rates When Detection Probabilities Are Less Than One. Ecology 83: 2248-2255.

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

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

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

Author

Rebecca A. Hutchinson

See also

unmarked, unmarkedFrameOccu, occu, computeMPLElambda, occuPEN_CV, nonparboot

Examples

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

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)
obsCovs(umf) <- covariates
head(umf)
#> Data frame representation of unmarkedFrame object.
#>   y.1 y.2      veght habitat    veght.1    veght.2 habitat.1 habitat.2
#> 1   0   0  1.0733096       A  1.0733096  3.3024986         A         A
#> 2   0   0  3.3024986       A -0.7308712  1.1855582         A         A
#> 3   0   0 -0.7308712       A  0.3037686  0.3758026         A         A
#> 4   0   0  1.1855582       A -1.0129665 -1.6261948         A         A
#> 5   0   1  0.3037686       A -0.7386708 -0.4478443         A         A
#> 6   1   0  0.3758026       A  0.2905940  0.8757997         A         A
#>  [ reached 'max' / getOption("max.print") -- omitted 4 rows ]
summary(umf)
#> unmarkedFrame Object
#> 
#> 100 sites
#> Maximum number of observations per site: 2 
#> Mean number of observations per site: 2 
#> Sites with at least one detection: 22 
#> 
#> Tabulation of y observations:
#>   0   1 
#> 167  33 
#> 
#> Site-level covariates:
#>      veght          habitat
#>  Min.   :-1.80746   A:50   
#>  1st Qu.:-0.75747   B:50   
#>  Median :-0.05057          
#>  Mean   : 0.07205          
#>  3rd Qu.: 0.62237          
#>  Max.   : 3.32295          
#> 
#> Observation-level covariates:
#>      veght          habitat
#>  Min.   :-1.80746   A:50   
#>  1st Qu.:-0.75747   B:50   
#>  Median :-0.05057          
#>  Mean   : 0.07205          
#>  3rd Qu.: 0.62237          
#>  Max.   : 3.32295          


# Fit some models
fmMLE <- occu(~veght+habitat ~veght+habitat, umf)
fm1pen <- occuPEN(~veght+habitat ~veght+habitat, umf,lambda=0.33,pen.type="Ridge")
fm2pen <- occuPEN(~veght+habitat ~veght+habitat, umf,lambda=1,pen.type="Bayes")

# MPLE:
fm3pen <- occuPEN(~veght+habitat ~veght+habitat, umf,lambda=0.5,pen.type="MPLE")
#> Warning: Supplied lambda does not match the computed value. Proceeding with the supplied lambda.
MPLElambda = computeMPLElambda(~veght+habitat ~veght+habitat, umf) 
fm4pen <- occuPEN(~veght+habitat ~veght+habitat, umf,lambda=MPLElambda,pen.type="MPLE")

# nonparametric bootstrap for uncertainty analysis:
fm1pen <- nonparboot(fm1pen,B=20) # should use more samples
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 50 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 50 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 50 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 50 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 49 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 60 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 49 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 53 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 47 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 46 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 38 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 46 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 44 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 55 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 47 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 51 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 46 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 48 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 62 sites have been discarded because of missing data.
#> Warning: Some observations have been discarded because corresponding covariates were missing.
#> Warning: 46 sites have been discarded because of missing data.
vcov(fm1pen,method="nonparboot")
#>                  psi(Int) psi(veght) psi(habitatB)     p(Int)    p(veght)
#> psi(Int)       0.54191131  0.8086550             0 -0.3296919  0.11124043
#> psi(veght)     0.80865501  1.4055070             0 -0.6409539  0.23168191
#> psi(habitatB)  0.00000000  0.0000000             0  0.0000000  0.00000000
#> p(Int)        -0.32969191 -0.6409539             0  0.4352372 -0.10749050
#> p(veght)       0.11124043  0.2316819             0 -0.1074905  0.11454844
#> p(habitatB)   -0.04670788  0.0337051             0 -0.1864737  0.01056793
#>               p(habitatB)
#> psi(Int)      -0.04670788
#> psi(veght)     0.03370510
#> psi(habitatB)  0.00000000
#> p(Int)        -0.18647374
#> p(veght)       0.01056793
#> p(habitatB)    0.37218475