inst/doc/unmarked.R

In unmarked: Models for Data from Unmarked Animals

## ----echo=FALSE---------------------------------------------------------------
options(rmarkdown.html_vignette.check_title = FALSE)

## ----echo=FALSE---------------------------------------------------------------
tab1 <- data.frame(
  Model=c("Occupancy", "Royle-Nichols", "Point Count", "Distance-sampling",
          "Generalized distance-sampling", "Arbitrary multinomial-Poisson",
          "Colonization-extinction", "Generalized multinomial-mixture"),
  `Fitting Function`=c("occu","occuRN","pcount","distsamp","gdistsamp",
                       "multinomPois","colext","gmultmix"),
  Data=c("unmarkedFrameOccu","unmarkedFrameOccu","unmarkedFramePCount",
         "unmarkedFrameDS","unmarkedFrameGDS","unmarkedFrameMPois",
         "unmarkedMultFrame","unmarkedFrameGMM"),
  Citation=c("@mackenzie_estimating_2002","@royle_estimating_2003",
             "@royle_n-mixture_2004","@royle_modeling_2004",
             "@chandlerEA_2011","@royle_generalized_2004",
             "@mackenzie_estimating_2003","@royle_generalized_2004"),
        check.names=FALSE)

knitr::kable(tab1, format='markdown', align="lccc",
             caption="Table 1. Models handled by unmarked.")

## -----------------------------------------------------------------------------
library(unmarked)
wt <- read.csv(system.file("csv","widewt.csv", package="unmarked"))
y <- wt[,2:4]
siteCovs <-  wt[,c("elev", "forest", "length")]
obsCovs <- list(date=wt[,c("date.1", "date.2", "date.3")],
    ivel=wt[,c("ivel.1",  "ivel.2", "ivel.3")])
wt <- unmarkedFrameOccu(y = y, siteCovs = siteCovs, obsCovs = obsCovs)
summary(wt)

## -----------------------------------------------------------------------------
pcru_raw <- read.csv(system.file("csv","frog2001pcru.csv", package="unmarked"))

sites <- unique(pcru_raw$RouteNumStopNum)
M <- length(sites)
J <- max(table(pcru_raw$RouteNumStopNum))
y <- JulianDate <- MinAfterSunset <- Wind <- Sky <- Temperature <- matrix(NA, M, J)
for (i in 1:length(sites)){
  datsub <- pcru_raw[pcru_raw$RouteNumStopNum == sites[i],]
  n <- nrow(datsub)
  y[i,1:n] <- datsub$Pcru
  JulianDate[i,1:n] <- datsub$JulianDate
  MinAfterSunset[i,1:n] <- datsub$MinAfterSunset
  Wind[i,1:n] <- datsub$Wind
  Sky[i,1:n] <- datsub$Sky
  Temperature[i,1:n] <- datsub$Temperature
}
obscovs <- list(JulianDate = JulianDate, MinAfterSunset = MinAfterSunset,
                Wind = Wind, Sky = Sky, Temperature = Temperature)
pcru <- unmarkedFrameOccu(y = y, obsCovs = obscovs)

## -----------------------------------------------------------------------------
fm1 <- occu(~1 ~1, pcru)
fm2 <- occu(~ scale(MinAfterSunset) + scale(Temperature) ~ 1, pcru)
fm2

## -----------------------------------------------------------------------------
predict(fm2, 'state')[1,]

## -----------------------------------------------------------------------------
nd <- data.frame(MinAfterSunset = 0, Temperature = 0)
round(predict(fm2, type = 'det', newdata = nd, appendData = TRUE), 2)

## -----------------------------------------------------------------------------
nd <- data.frame(MinAfterSunset = 0, Temperature = -2:2)
round(predict(fm2, type = 'det', newdata = nd, appendData=TRUE), 2)

## ----eval=FALSE---------------------------------------------------------------
# confint(fm2, type='det')
# confint(fm2, type='det', method = "profile")

## ----echo=FALSE---------------------------------------------------------------
confint(fm2, type='det')
nul <- capture.output(ci <- confint(fm2, type='det', method = "profile"))
ci

## -----------------------------------------------------------------------------
fms <- fitList('psi(.)p(.)' = fm1, 'psi(.)p(Time+Temp)' = fm2)
modSel(fms)
predict(fms, type='det', newdata = nd)

## ----warning=FALSE------------------------------------------------------------
chisq <- function(fm) {
    umf <- fm@data
    y <- umf@y
    y[y>1] <- 1
    fv <- fitted(fm)
    sum((y-fv)^2/(fv*(1-fv)), na.rm=TRUE)
    }

set.seed(1)
(pb <- parboot(fm2, statistic=chisq, nsim=100, parallel=FALSE))

## -----------------------------------------------------------------------------
re <- ranef(fm2)
EBUP <- bup(re, stat="mode")
sum(EBUP) / numSites(pcru)