As noted in the lecture, multi-season occupancy depends on data collected under the Robust Design, so that S sites are sampled on each of T primary occasions and J secondary occasions for J*T total samples per site (in general J can vary among primary occasions but for now we assume all J's are the same). The parameters needed to model these data are now:
psi(1) - initial occupancy probability at the beginning of the study
ext(t) - probability that an occupied site becomes unoccupied between primary periods t and t+1
col(t) - probability that an unoccupied site becomes occupied between primary periods t and t+1
p(j, j=1, J*T) - probability of detection of an occupied side for each of the J*T samples
GIven the first 3 parameters, occupancy after time t=1 is modeled as
psi(t+1)=psi(t)*[1-ext(t)]+ [1-psi(t)]*col(t)
The parameters, in turn, can be modeled in terms of site and occasion-specific covariates or grouping factors.
Multi-season occupancy modeling – Grand Skink example
I illustrate multi-season occupancy modeling with an example from Grand Skinks, an endangered reptile endemic to the Central Otago region of the South Island of New Zealand. The data are provided as one of the example data sets with the installation of the program PRESENCE from the PWRC site, but we will instead use the colext() function in the library unmarked. The data are located in a csv file containing the presence/absence records (columns B-P) over 5 years (“seasons”) and with 3 surveys per season, at each of 352 rocky outcrops (“tors”). A single site-specific covariate is provided in column R, which is an variable indicating whether the site is in pasture (1) or native grassland (0).
First we read the data in, specify site-level covariates and time effects at the seasonal (yearly) and survey levels:
> data<-read.csv("Grand_Skinks.csv", na.strings="-")
> y<-data[,2:16]
> S <- nrow(data) # number of sites
> J <- 3 # number of secondary sampling occasions
> T <- 5 # number of primary periods
>
> #site level covariates
> site.cov<-data.frame(pasture=data[,18])
> #dummy variables to createseasonal effects (need 2 for 3 occasions)
> year<-as.factor(rep(c(1,2,3,4,5),S))
>
> yearly.cov<-data.frame(year)
> #additive time effects within season. to model effects change between seasons create an interaction term
> secondary<-as.factor(rep(c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3),S))
>
> second.cov<-data.frame(secondary)
These are then combined into a specialized dataframe:
> skink<-unmarkedMultFrame(y=y,siteCovs=site.cov, numPrimary=5,yearlySiteCovs=yearly.cov,obsCovs=second.cov)
We can run an initial, simple model and obtain estimates
>fm1<-colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
pformula = ~1, data=skink)
Warning message:
14 sites have been discarded because of missing data.
> #back transformations
> backTransform(fm1,"det")
Backtransformed linear combination(s) of Detection estimate(s)
Estimate SE LinComb (Intercept)
0.689 0.0231 0.796 1
Transformation: logistic
> backTransform(fm1,"psi")
Backtransformed linear combination(s) of Initial estimate(s)
Estimate SE LinComb (Intercept)
0.395 0.0315 -0.428 1
Transformation: logistic
> backTransform(fm1,"col")
Backtransformed linear combination(s) of Colonization estimate(s)
Estimate SE LinComb (Intercept)
0.0683 0.0119 -2.61 1
Transformation: logistic
> backTransform(fm1,"ext")
Backtransformed linear combination(s) of Extinction estimate(s)
Estimate SE LinComb (Intercept)
0.1 0.0183 -2.2 1
Transformation: logistic
Let's proceed to construct several more models that permit combinations of time and site covariates effects, and do AIC model comparisons
> #additional models
> #constant initial psi , gamma , and epsilon
> #secondary occasion specific detection
> fm2<-colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
+ pformula = ~secondary, data=skink)
Warning message:
14 sites have been discarded because of missing data.
>
> #constant initial psi , gamma , and epsilon
> #primary occasion specific detection
> fm3<-colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
+ pformula = ~year, data=skink)
Warning message:
14 sites have been discarded because of missing data.
>
> #constant initial psi , gamma , and epsilon
> #primary*secondary occasion specific detection
> fm4<-colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
+ pformula = ~year*secondary, data=skink)
Warning message:
14 sites have been discarded because of missing data.
>
>
> #initial psi as function of pasture
> fm5<-colext(psiformula = ~pasture, gammaformula = ~1, epsilonformula = ~1,
+ pformula = ~1, data=skink)
Warning message:
14 sites have been discarded because of missing data.
>
> #constant initial psi , gamma and epsilon time dependent
> #
> fm6<-colext(psiformula = ~1, gammaformula = ~year, epsilonformula = ~year,
+ pformula =~1, data=skink)
Warning message:
14 sites have been discarded because of missing data.
>
> fm7<-colext(psiformula = ~pasture, gammaformula = ~year, epsilonformula = ~year,
+ pformula =~year+pasture, data=skink,control=list(maxit=1e4))
Warning message:
14 sites have been discarded because of missing data.
>
> fm8<-colext(psiformula = ~pasture, gammaformula = ~1, epsilonformula = ~1,
+ pformula =~year+pasture, data=skink,control=list(maxit=1e4))
Warning message:
14 sites have been discarded because of missing data.
>
> fm9<-colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~pasture,
+ pformula =~year, data=skink,control=list(maxit=1e4))
Warning message:
14 sites have been discarded because of missing data.
>
> fm10<-colext(psiformula = ~pasture, gammaformula = ~1, epsilonformula = ~pasture,
+ pformula =~year, data=skink,control=list(maxit=1e4))
Warning message:
14 sites have been discarded because of missing data.
>
>
>
> fms <- fitList('psi0(.)col(.)ext(.)p(.)' = fm1,'psi0(.)col(.)ext(.)p(s)'=fm2, 'psi0(.)col(.)ext(.)p(t)'=fm3,'psi0(.)col(.)ext(.)p(t*s)'=fm4,'psi0(past)col(.)ext(.)p(.)'=fm5,'psi0(.)col(t)ext(t)p(.)'=fm6,'psi0(past)col(t)ext(t)p(t+past)'=fm7,'psi0(past)col(.)ext(.)p(t+past)'=fm8,'psi0(.)col(.)ext(past)p(t)'=fm9,'psi0(past)col(.)ext(past)p(t)'=fm10)
There were 11 warnings (use warnings() to see them)
> mod.sel<-modSel(fms)
Warning message:
In sqrt(diag(vcov(x, altNames = TRUE))) : NaNs produced
> mod.sel
nPars AIC delta AICwt cumltvWt
psi0(past)col(.)ext(past)p(t) 10 1744.71 0.00 9.0e-01 0.90
psi0(past)col(.)ext(.)p(t+past) 10 1750.05 5.34 6.2e-02 0.96
psi0(past)col(.)ext(.)p(.) 5 1751.54 6.83 3.0e-02 0.99
psi0(past)col(t)ext(t)p(t+past) 16 1753.52 8.81 1.1e-02 1.00
psi0(.)col(.)ext(past)p(t) 9 1766.66 21.95 1.5e-05 1.00
psi0(.)col(.)ext(.)p(s) 6 1769.05 24.34 4.6e-06 1.00
psi0(.)col(.)ext(.)p(t*s) 18 1769.82 25.11 3.2e-06 1.00
psi0(.)col(.)ext(.)p(t) 8 1771.98 27.27 1.1e-06 1.00
psi0(.)col(.)ext(.)p(.) 4 1775.02 30.31 2.4e-07 1.00
psi0(.)col(t)ext(t)p(.) 10 1777.43 32.71 7.1e-08 1.00
The warning messages simply indicated that there were missing values in the data and are not cause for alarm. Note that for the last 2 models I included a 'maxit' option-- this was because these 2 models failed to converge in the default number of iterations allowed.
The model selection procedure indicates that the top 3 models all contain an effect on initial occupancy of pasture. If we wish we can explore this relationship further, by model averaging this parameter (or better yet, predictions for specified levels of pasture =1 or 0). This time I am using a built in predict() function in unmarked that automatically performs the model averaging over models in the specified list of models that has already been run.
Here are data and R script for this example.