inst/doc/secrlinear-vignette.R

In secrlinear: Spatially Explicit Capture-Recapture for Linear Habitats

## ----setup, message = FALSE, warning = FALSE----------------------------------
library(secrlinear)   # also loads secr
options(digits = 4)   # for more readable output
inputdir <- system.file("extdata", package = "secrlinear")

## ----readarvicola, eval = TRUE------------------------------------------------
captfile <- paste0(inputdir, "/Jun84capt.txt")
trapfile <- paste0(inputdir, "/glymetrap.txt")
arvicola <- read.capthist(captfile, trapfile, covname = "sex")

## ----readglyme, eval = TRUE---------------------------------------------------
habitatmap <-  paste0(inputdir, "/glymemap.txt")
glymemask <- read.linearmask(file = habitatmap, spacing = 4)

## ----plotglyme, eval = TRUE, fig.width = 7, fig.height = 3.5------------------
par(mar = c(1,1,4,1))
plot(glymemask)
plot(arvicola, add = TRUE, tracks = TRUE)
plot(traps(arvicola), add = TRUE)

## ----fit1, eval = TRUE, warning = FALSE---------------------------------------
# 2-D habitat, Euclidean distance
fit2DEuc <- secr.fit(arvicola, buffer = 200, trace = FALSE)   

# 1-D habitat, Euclidean distance
fit1DEuc <- secr.fit(arvicola, mask = glymemask, trace = FALSE)

# 1-D habitat, river distance
fit1DNet <- secr.fit(arvicola, mask = glymemask, trace = FALSE,
                details = list(userdist = networkdistance))

## ----predict, eval = TRUE-----------------------------------------------------
predict(fit2DEuc)
predict(fit1DEuc)
predict(fit1DNet)

## ----silvermask, eval = TRUE--------------------------------------------------
habitatmap <- paste0(inputdir, "/silverstream.shp")
silverstreammask <- read.linearmask(file = habitatmap, spacing = 50)
par(mar = c(1,1,1,1))
plot(silverstreammask)

## ----networklength, eval = TRUE-----------------------------------------------
sldf <- attr(silverstreammask, "SLDF")
networklength <- sum(sp::SpatialLinesLengths(sldf)) / 1000   # km
discrepancy <- networklength - masklength(silverstreammask)  # km

## ----silvermask2, eval = FALSE------------------------------------------------
#  habitatmap <- paste0(inputdir, "/silverstream.shp")
#  silverstreamsf <- st_read(habitatmap)
#  silverstreamSLDF <- as(silverstreamsf, 'Spatial')
#  silverstreammask <- read.linearmask(data = silverstreamSLDF, spacing = 50)

## ----dataframemask, eval=TRUE-------------------------------------------------
x <- seq(0, 4*pi, length = 200)
xy <- data.frame(x = x*100, y = sin(x)*300)
linmask <- read.linearmask(data = xy, spacing = 20)

## ----plotlinmask, eval = TRUE-------------------------------------------------
plot(linmask)

## ----showpath, eval = FALSE---------------------------------------------------
#  # start interactive session and click on two points
#  showpath(silverstreammask, lwd = 3)

## ----makeline, eval = TRUE----------------------------------------------------
trps <- make.line(linmask, detector = "proximity", n = 40, startbuffer = 0, by = 300,
                  endbuffer = 80, cluster = c(0,40,80), type = 'randomstart')

## ----plotline, eval = TRUE, fig.width = 7, fig.height = 3.5-------------------
plot(linmask)
plot(trps, add = TRUE, detpar = list(pch = 16, cex = 1.2, col='red'))

## ----snappoints, eval = FALSE-------------------------------------------------
#  plot(silverstreammask)
#  loc <- locator(30)
#  xy <- snapPointsToLinearMask(data.frame(loc), silverstreammask)
#  tr <- read.traps(data = xy, detector = 'multi')
#  plot(tr, add = TRUE)

## ----transect, eval = FALSE---------------------------------------------------
#  transects <- read.traps('transectxy.txt', detector = 'transect')
#  capt <- read.table('capt.txt')
#  tempCH <- make.capthist(capt, transects, fmt = 'XY')
#  tempCH <- snip(tempCH, by = 100)   # for 100-m segments
#  CH <- reduce(tempCH, outputdetector = "count")

## ----silvertrps, eval = TRUE, echo = FALSE------------------------------------
trapfile <- paste0(inputdir, "/silverstreamtraps.txt")
tr <- read.traps(trapfile, detector = "multi")

## ----simCH, eval = TRUE, cache = TRUE-----------------------------------------
# simulate population of 2 animals / km
pop <- sim.linearpopn(mask = silverstreammask, D = 2)
# simulate detections using network distances
CH <- sim.capthist(traps = tr, popn = pop, noccasions = 4,
                   detectpar = list(g0 = 0.25, sigma = 500), 
                   userdist = networkdistance)
summary(CH)    # detector spacing uses Euclidean distances

## ----plotsim, eval=TRUE-------------------------------------------------------
# and plot the simulated detections...
par(mar = c(1,1,1,1))
plot(silverstreammask)
plot(CH, add = TRUE, tracks = TRUE, varycol = TRUE, rad = 100, cappar = list(cex = 2))
plot(tr, add = TRUE)

## ----sfit, eval = FALSE-------------------------------------------------------
#  userd <- networkdistance(tr, silverstreammask)
#  userd[!is.finite(userd)] <- 1e8  # testing
#  sfit <- secr.fit(CH, mask = silverstreammask, details = list(userdist = userd))
#  predict(sfit)

## ----regionN, eval = TRUE-----------------------------------------------------
region.N(fit2DEuc)
region.N(fit1DNet)

## ----plotregion, eval = TRUE, fig.width = 6.5, fig.height=3-------------------
par(mfrow = c(1,2), mar = c(1,1,1,1))
plot(fit2DEuc$mask)
plot(traps(arvicola), add = TRUE)
mtext(side = 3,line = -1.8, "fit2DEuc$mask", cex = 0.9)
plot(fit1DNet$mask)
plot(traps(arvicola), add = TRUE)
mtext(side = 3,line = -1.8,"fit1DNet$mask", cex = 0.9)

## ----derived, eval = TRUE-----------------------------------------------------
derived(fit2DEuc)
derived(fit1DNet)

## ----covariates, eval = FALSE-------------------------------------------------
#  # interactively obtain LineID for central 'spine' by clicking on
#  # each component line in plot
#  tmp <- getLineID(silverstreammask)
#  # extract coordinates of 'spine'
#  spine <- subset(silverstreammask, LineID = tmp$LineID)
#  # obtain network distances to spine and save for later use
#  netd <- networkdistance(spine, silverstreammask)  # matrix dim = c(nrow(spine), nrow(mask))
#  dfs <- apply(netd, 2, min) / 1000  # km
#  covariates(silverstreammask)$dist.from.spine <- dfs

## ----plotcovariate, eval = FALSE----------------------------------------------
#  par(mar=c(1,1,1,4))
#  plot(silverstreammask, covariate = 'dist.from.spine', col = topo.colors(13),
#       cex = 1.5, legend = FALSE)
#  strip.legend('right', legend = seq(0, 6.5, 0.5), col = topo.colors(13),
#               title = 'dist.from.spine km', height = 0.35)
#  plot(spine, add = TRUE, linecol = NA, cex = 0.3)

## ----checkmoves, eval = FALSE, strip.white = TRUE-----------------------------
#  # initially OK (no movement > 1000 m)--
#  checkmoves(arvicola, mask = glymemask, accept = c(0,1000))
#  # deliberately break graph of linear mask
#  attr(glymemask, 'graph')[200:203,201:204] <- NULL
#  # no longer OK --
#  out <- checkmoves(arvicola, mask = glymemask, accept = c(0,1000))
#  # display captures of animals 32 and 35 whose records span break
#  out$df

## ----showedges, eval = FALSE--------------------------------------------------
#  # problem shows up where voles recaptured either side of break:
#  showedges(glymemask, col = 'red', lwd = 6)
#  plot(out$CH, add = TRUE, tracks = TRUE, rad=8,cappar=list(cex=1.5))
#  pos <- traps(arvicola)['560.B',]
#  text(pos$x+5, pos$y+80, 'break', srt=90, cex=1.1)

## ----plotglymeedges, eval = FALSE---------------------------------------------
#  plot(glymemask)
#  replot(glymemask)  # click on corners to zoom in
#  showedges(glymemask, col = 'red', lwd = 2, add=T)
#  glymemask <- addedges(glymemask)

## ----linearHR, eval = FALSE---------------------------------------------------
#  par(mfrow = c(1,1), mar = c(1,1,1,5))
#  plot(silverstreammask)
#  centres <- data.frame(locator(4))
#  OK <- networkdistance(centres, silverstreammask) < 1000
#  for (i in 1:nrow(OK)) {
#      m1 <- subset(silverstreammask, OK[i,])
#      plot(m1, add = TRUE, col = 'red', cex = 1.7)
#      ml <- masklength(m1)
#      points(centres, pch = 16, col = 'yellow', cex = 1.4)
#      text (1406000, mean(m1$y), paste(ml, 'km'), cex = 1.2)
#  }
#  

## ----secrdesign, eval = TRUE, warning = FALSE---------------------------------
library(secrdesign)

# create a habitat geometry
x <- seq(0, 4*pi, length = 200)
xy <- data.frame(x = x*100, y = sin(x)*300)
linmask <- read.linearmask(data = xy, spacing = 5)

# define two possible detector layouts
trp1 <- make.line(linmask, detector = "proximity", n = 80, start = 200, by = 30)
trp2 <- make.line(linmask, detector = "proximity", n = 40, start = 200, by = 60)
trplist <- list(spacing30 = trp1, spacing60 = trp2)

# create a scenarios dataframe
scen1 <- make.scenarios(D = c(50,200), trapsindex = 1:2, sigma = 25, g0 = 0.2)

# we specify the mask, rather than construct it 'on the fly',
# we will use a non-Euclidean distance function for both
# simulating detections and fitting the model...
det.arg <- list(userdist = networkdistance)
fit.arg <- list(details = list(userdist = networkdistance))

# run the scenarios and summarise results
sims1 <- run.scenarios(nrepl = 50, trapset = trplist, maskset = linmask, 
    det.args = list(det.arg), fit.args = list(fit.arg),
    scenarios = scen1, seed = 345, fit = FALSE)
summary(sims1)

## ----sims2, eval = FALSE------------------------------------------------------
#  sims2 <- run.scenarios(nrepl = 5, trapset = trplist, maskset = linmask,
#       det.args = list(det.arg), scenarios = scen1, seed = 345, fit = TRUE)
#  summary(sims2)

## ----appendix, eval = FALSE---------------------------------------------------
#  # It is efficient to pre-compute a matrix of distances between traps (rows)
#  # and mask points (columns)
#  distmat <- networkdistance (traps(arvicola), glymemask, glymemask)
#  
#  # Morning and evening trap checks as a time covariate
#  tcov <- data.frame(ampm = rep(c("am","pm"),3))
#  
#  glymefit1 <- secr.fit(arvicola, mask = glymemask, trace = FALSE,
#                        details = list(userdist = distmat),
#                        model = g0~1, hcov = "sex")
#  glymefit2 <- secr.fit(arvicola, mask = glymemask, trace = FALSE,
#                        details = list(userdist = distmat),
#                        model = g0~ampm, timecov = tcov, hcov = "sex")
#  glymefit3 <- secr.fit(arvicola, mask = glymemask, trace = FALSE,
#                        details = list(userdist = distmat),
#                        model = g0~ampm + h2, timecov = tcov, hcov = "sex")
#  glymefit4 <- secr.fit(arvicola, mask = glymemask, trace = FALSE,
#                        details = list(userdist = distmat),
#                        model = list(sigma~h2, g0~ampm + h2),
#                        timecov = tcov, hcov = "sex")
#  
#  fitlist <- secrlist(glymefit1, glymefit2, glymefit3, glymefit4)
#  # dropping the detectfn (halfnormal) column to save space...
#  AIC(fitlist)[,-2]
#  #                                       model npar logLik   AIC  AICc dAICc AICcwt
#  # glymefit4 D~1 g0~ampm + h2 sigma~h2 pmix~h2    7 -322.5 659.1 665.3  0.00      1
#  # glymefit3  D~1 g0~ampm + h2 sigma~1 pmix~h2    6 -347.3 706.7 711.1 45.80      0
#  # glymefit2       D~1 g0~ampm sigma~1 pmix~h2    5 -353.5 717.0 720.0 54.66      0
#  # glymefit1          D~1 g0~1 sigma~1 pmix~h2    4 -356.8 721.6 723.5 58.20      0
#  
#  # summaries of estimated density and sex ratio under different models
#  options(digits=3)
#  
#  # model does not affect density estimate
#  collate(fitlist, perm = c(2,3,1,4))[,,1,"D"]
#  #           estimate SE.estimate  lcl  ucl
#  # glymefit1     26.5        5.27 18.0 39.0
#  # glymefit2     26.4        5.26 18.0 38.9
#  # glymefit3     26.3        5.25 17.9 38.8
#  # glymefit4     27.2        5.45 18.5 40.2
#  
#  # model does affect the estimate of sex ratio (here proportion female)
#  collate(fitlist, perm=c(2,3,1,4))[,,1,"pmix"]
#  #           estimate SE.estimate   lcl   ucl
#  # glymefit1    0.615      0.0954 0.421 0.779
#  # glymefit2    0.615      0.0954 0.421 0.779
#  # glymefit3    0.634      0.0938 0.439 0.793
#  # glymefit4    0.669      0.0897 0.477 0.817
#  
#  # predictions from best model
#  newdata <- expand.grid(ampm = c("am", "pm"), h2 = c("F", "M"))
#  predict(glymefit4, newdata = newdata)
#  
#  # $`ampm = am, h2 = F`
#  #        link estimate SE.estimate    lcl    ucl
#  # D       log   27.239      5.4478 18.477 40.158
#  # g0    logit    0.218      0.0463  0.141  0.322
#  # sigma   log   13.624      1.8764 10.414 17.823
#  # pmix  logit    0.669      0.0897  0.477  0.817
#  #
#  # $`ampm = pm, h2 = F`
#  #        link estimate SE.estimate     lcl    ucl
#  # D       log   27.239      5.4478 18.4768 40.158
#  # g0    logit    0.116      0.0293  0.0694  0.186
#  # sigma   log   13.624      1.8764 10.4136 17.823
#  # pmix  logit    0.669      0.0897  0.4774  0.817
#  #
#  # $`ampm = am, h2 = M`
#  #        link estimate SE.estimate     lcl    ucl
#  # D       log   27.239      5.4478 18.4768 40.158
#  # g0    logit    0.153      0.0392  0.0908  0.246
#  # sigma   log   70.958     10.0551 53.8247 93.545
#  # pmix  logit    0.331      0.0897  0.1829  0.523
#  #
#  # $`ampm = pm, h2 = M`
#  #        link estimate SE.estimate     lcl    ucl
#  # D       log  27.2394      5.4478 18.4768 40.158
#  # g0    logit   0.0782      0.0201  0.0468  0.128
#  # sigma   log  70.9581     10.0551 53.8247 93.545
#  # pmix  logit   0.3311      0.0897  0.1829  0.523

## ----derivedapp, eval = FALSE-------------------------------------------------
#  derived(glymefit4, distribution = 'binomial')
#  #     estimate SE.estimate   lcl   ucl    CVn     CVa    CVD
#  # esa   0.9545          NA    NA    NA     NA      NA     NA
#  # D    27.2396       2.867 22.17 33.46 0.1038 0.01747 0.1053