weta: Occupancy data for Mahoenui Giant Weta

In RMark: R Code for Mark Analysis

Occupancy data for Mahoenui Giant Weta

Description

An occupancy data set for modelling presence/absence data for salamanders.

Format

A data frame with 72 observations (sites) on the following 7 variables.

ch

a character vector containing the presence (1) and absence (0), or (.) not visited for each of 5 visits to the site

Browse

0/1 dummy variable to indicate browsing

Obs1

observer number for visit 1; . used when site not visited

Obs2

observer number for visit 2; . used when site not visited

Obs3

observer number for visit 3; . used when site not visited

Obs4

observer number for visit 4; . used when site not visited

Obs5

observer number for visit 5; . used when site not visited

Details

This is a data set that accompanies program PRESENCE and is explained on pages 116-122 of MacKenzie et al. (2006).

References

MacKenzie, D.I., Nichols, J. D., Royle, J.A., Pollock, K.H., Bailey, L.L., and Hines, J.E. 2006. Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurence. Elsevier, Inc. 324p.

Examples

#  The data can be imported with the following command using the
#  tab-delimited weta.txt file in the data subdirectory.
#   weta=import.chdata("weta.txt",field.types=c(rep("f",6)))
#  Below is the first few lines of the data file that was constructed
#  from the .xls file that accompanies PRESENCE.
#ch	Browse	Obs1	Obs2	Obs3	Obs4	Obs5
#0000.	1	1	3	2	3	.
#0000.	1	1	3	2	3	.
#0001.	1	1	3	2	3	.
#0000.	0	1	3	2	3	.
#0000.	1	1	3	2	3	.
#0000.	0	1	3	2	3	.
#

# This example is excluded from testing to reduce package check time
# retrieve weta data
data(weta)
# Create function to fit the 18 models in the book
fit.weta.models=function()
{
#  use make.time.factor to create time-varying dummy variables Obs1 and Obs2
#  observer 3 is used as the intercept
   weta=make.time.factor(weta,"Obs",1:5,intercept=3)
#  Process data and use Browse covariate to group sites; it could have also
#  been used an individual covariate because it is a 0/1 variable.
   weta.process=process.data(weta,model="Occupancy",groups="Browse")
   weta.ddl=make.design.data(weta.process)
#  time factor variable copied to Day to match names used in book
   weta.ddl$p$Day=weta.ddl$p$time
# Define p models
   p.dot=list(formula=~1)
   p.day=list(formula=~Day)
   p.obs=list(formula=~Obs1+Obs2)
   p.browse=list(formula=~Browse)
   p.day.obs=list(formula=~Day+Obs1+Obs2)
   p.day.browse=list(formula=~Day+Browse)
   p.obs.browse=list(formula=~Obs1+Obs2+Browse)
   p.day.obs.browse=list(formula=~Day+Obs1+Obs2+Browse)
# Define Psi models
   Psi.dot=list(formula=~1)
   Psi.browse=list(formula=~Browse)
# Create model list
   cml=create.model.list("Occupancy")
# Run and return marklist of models
   return(mark.wrapper(cml,data=weta.process,ddl=weta.ddl,delete=TRUE))
}
weta.models=fit.weta.models()
# Modify the model table to show -2lnl and use AIC rather than AICc
weta.models$model.table=model.table(weta.models,use.AIC=TRUE,use.lnl=TRUE)
# Show new model table which duplicates the results except they have
# some type of error with the model Psi(.)P(Obs+Browse) which should have
# 5 parameters rather than 4 and the -2lnl also doesn't agree with the results here
weta.models
#
# display beta vcv matrix of the Psi parameters (intercept + browse=1)
# matches what is shown on pg 122 of Occupancy book
weta.models[[7]]$result$beta.vcv[8:9,8:9]
# compute variance-covariance matrix of Psi0(6; unbrowsed) ,Psi1(7; browsed)
vcv.psi=get.real(weta.models[[7]],"Psi",vcv=TRUE)$vcv.real
vcv.psi
# Compute proportion unbrowsed and browsed
prop.browse=c(37,35)/72
prop.browse
# compute std error of overall estimate as shown on pg 121-122
sqrt(sum(prop.browse^2*diag(vcv.psi)))
# compute std error and correctly include covariance between Psi0 and Psi1
sqrt( t(prop.browse) %*% vcv.psi %*% prop.browse )
# show missing part of variance 2 times cross-product of prop.browse * covariance
2*prod(prop.browse)*vcv.psi[1,2]
sqrt(sum(prop.browse^2*diag(vcv.psi))+2*prod(prop.browse)*vcv.psi[1,2])

RMark documentation built on Aug. 14, 2022, 1:05 a.m.