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

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

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.

Arguments

psiformula

Right-hand sided formula for the initial probability of occupancy at each site.

gammaformula

Right-hand sided formula for colonization probability.

epsilonformula

Right-hand sided formula for extinction probability.

pformula

Right-hand sided formula for detection probability.

data

unmarkedMultFrame object that supplies the data (see unmarkedMultFrame).

starts

optionally, initial values for parameters in the optimization.

method

Optimization method used by optim.

se

logical specifying whether or not to compute standard errors.

...

Additional arguments to optim, such as lower and upper bounds

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
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    1    1    0    0    0    0
#> [2,]    0    0    0    0    0    0
#> [3,]    1    1    1    1    1    0
#> [4,]    1    0    1    0    0    1

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

yearly.site.covs <- list(
   year = matrix(c(
      'year1', 'year2',
      'year1', 'year2',
      'year1', 'year2',
      'year1', 'year2'), nrow=R, ncol=T, byrow=TRUE)
      )
yearly.site.covs
#> $year
#>      [,1]    [,2]   
#> [1,] "year1" "year2"
#> [2,] "year1" "year2"
#> [3,] "year1" "year2"
#> [4,] "year1" "year2"
#> 

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
#> $x4
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]   -1    0    1   -1    1    1
#> [2,]   -2    0    0    0    0    2
#> [3,]   -3    1    0    1    1    2
#> [4,]    0    0    0    0    1   -1
#> 
#> $x5
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] "a"  "b"  "c"  "a"  "b"  "c" 
#> [2,] "d"  "b"  "a"  "d"  "b"  "a" 
#> [3,] "a"  "a"  "c"  "d"  "b"  "a" 
#> [4,] "a"  "b"  "a"  "d"  "b"  "a" 
#> 

umf <- unmarkedMultFrame(y=y, siteCovs=site.covs,
    yearlySiteCovs=yearly.site.covs, obsCovs=obs.covs,
    numPrimary=2)                  # organize data
#> Warning: obsCovs contains characters. Converting them to factors.
#> Warning: yearlySiteCovs contains characters. Converting them to factors.
umf                                # look at data
#> Data frame representation of unmarkedFrame object.
#>   y.1 y.2 y.3 y.4 y.5 y.6 x1 x2 year.1 year.2 x4.1 x4.2 x4.3 x4.4 x4.5 x4.6
#> 1   1   1   0   0   0   0  1  A  year1  year2   -1    0    1   -1    1    1
#> 2   0   0   0   0   0   0  2  B  year1  year2   -2    0    0    0    0    2
#> 3   1   1   1   1   1   0  3  A  year1  year2   -3    1    0    1    1    2
#> 4   1   0   1   0   0   1  4  B  year1  year2    0    0    0    0    1   -1
#>   x5.1 x5.2 x5.3 x5.4 x5.5 x5.6
#> 1    a    b    c    a    b    c
#> 2    d    b    a    d    b    a
#> 3    a    a    c    d    b    a
#> 4    a    b    a    d    b    a
summary(umf)                       # summarize
#> unmarkedFrame Object
#> 
#> 4 sites
#> Maximum number of observations per site: 6 
#> Mean number of observations per site: 6 
#> Number of primary survey periods: 2 
#> Number of secondary survey periods: 3 
#> Sites with at least one detection: 3 
#> 
#> Tabulation of y observations:
#>  0  1 
#> 14 10 
#> 
#> Site-level covariates:
#>        x1       x2   
#>  Min.   :1.00   A:2  
#>  1st Qu.:1.75   B:2  
#>  Median :2.50        
#>  Mean   :2.50        
#>  3rd Qu.:3.25        
#>  Max.   :4.00        
#> 
#> Observation-level covariates:
#>        x4         x5    
#>  Min.   :-3.000   a:10  
#>  1st Qu.: 0.000   b: 7  
#>  Median : 0.000   c: 3  
#>  Mean   : 0.125   d: 4  
#>  3rd Qu.: 1.000         
#>  Max.   : 2.000         
#> 
#> Yearly-site-level covariates:
#>     year  
#>  year1:4  
#>  year2:4  
fm <- colext(~1, ~1, ~1, ~1, umf)  # fit a model
fm
#> 
#> Call:
#> colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1, 
#>     pformula = ~1, data = umf)
#> 
#> Initial:
#>  Estimate   SE     z P(>|z|)
#>      1.16 1.21 0.956   0.339
#> 
#> Colonization:
#>  Estimate   SE      z P(>|z|)
#>     -6.61 28.5 -0.232   0.817
#> 
#> Extinction:
#>  Estimate   SE      z P(>|z|)
#>    -0.783 1.34 -0.583    0.56
#> 
#> Detection:
#>  Estimate   SE    z P(>|z|)
#>     0.615 0.58 1.06   0.289
#> 
#> AIC: 35.18469 



if (FALSE) { # \dontrun{
# 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)))
} # }