colext.RdEstimate 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, ...)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).
unmarkedFitColExt object describing model fit.
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.
Right-hand sided formula for the initial probability of occupancy at each site.
Right-hand sided formula for colonization probability.
Right-hand sided formula for extinction probability.
Right-hand sided formula for detection probability.
unmarkedMultFrame object that supplies the data (see unmarkedMultFrame).
optionally, initial values for parameters in the optimization.
Optimization method used by optim.
logical specifying whether or not to compute standard errors.
Additional arguments to optim, such as lower and upper bounds
# 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)))
} # }