Sighting-only camera density estimation

[doranmye]

Hi secr group, 

I'm looking for any help properly executing a SECR sighting-only model (secr version 3.1.0) on a remote camera grid. I am estimating density of lynx with a partially marked population for my master's thesis. Having a great time! But making some mistakes as I go. 

I have 75 "count" detectors (cameras) with varying effort over one session and 12 occasions (traphist). There were 26 known marked animals within the study area - 15 detected by cameras, 11 not detected. I have not attempted yet to model the marking process. There are 25 detections of marked and identified individuals (capthist); 9 detections of marked and unidentified individuals (Tm); and 267 unmarked detections (Tu). I have zero unresolved detections (Tn). 

I have concurrent genetic data from hair snares and had success with a simple SECR model, giving me density of 14 animals/ 100 km sq. This estimate makes a lot of biological sense (there are a lot of lynx right now) so I'm happy with that estimate for now. 

The "sighting-only" model for cameras - same population at the same time -  is giving me an estimate of <1 animal/ 100 km sq. This doesn't make sense. In addition, I am getting lots of warnings as I go regarding usage and sightings. As I described above, there are at least 26 animals (marked) in my ~300 km sq study area. After getting this result, I went back to try to find where I went wrong, and haven't been able to find anything. Anyone with experience with SECR mark-resight: can you spot where I went wrong? Here I'm including my data files, code, and a description of what I have tried. I am just after the null model for now. 

__________________________ 

#build mark-resight caphist object

>ResightJanApr17CH <- read.capthist('cams_janapr17_marked.txt', 'cams_janapr17_traps_use.txt', detector = "count", markocc = c(0,0,0,0,0,0,0,0,0,0,0,0))

            No errors found :-)

>ResightJanApr17CH <- addSightings(ResightJanApr17CH, unmarked = 'cams_janapr17_tu2.txt', nonID = 'cams_janapr17_tm2.txt', verify = TRUE)

            No errors found :-)

>summary(ResightJanApr17CH)

Object class      capthist 

Detector type     count 

Detector number   75 

Average spacing   1311.27 m 

x-range           183485.5 212318 m 

y-range           720452 739149.9 m 

 Usage range by occasion

    1 2 3 4 5 6 7 8 9 10 11 12

min 0 0 0 0 0 0 0 0 0  0  0  0

max 1 1 1 1 1 1 1 1 1  1  1  1

Marking occasions

 1 2 3 4 5 6 7 8 9 10 11 12

 0 0 0 0 0 0 0 0 0  0  0  0

Counts by occasion 

                   1  2  3  4  5  6  7  8  9 10 11 12 Total

n                  2  1  3  2  0  0  1  1  3  1  3  6    23

u                  2  0  2  2  0  0  1  1  3  0  2  2    15

f                 10  3  1  1  0  0  0  0  0  0  0  0    15

M(t+1)             2  2  4  6  6  6  7  8 11 11 13 15    15

losses             0  0  0  0  0  0  0  0  0  0  0  0     0

detections         2  1  3  2  0  0  1  1  3  2  3  7    25

detectors visited  1  1  3  2  0  0  1  1  3  2  3  5    22

detectors used    63 65 70 69 44 42 40 57 62 65 66 66   709

Empty histories :  11 

Sightings by occasion 

          1  2  3  4  5  6  7  8  9 10 11 12 Total

ID        2  1  3  2  0  0  1  1  3  2  3  7    25

Not ID    0  0  1  0  0  0  2  0  1  1  1  3     9

Unmarked  3 14 26 30 15 15  9 20 28 26 34 47   267

Uncertain 0  0  0  0  0  0  0  0  0  0  0  0     0

Total     5 15 30 32 15 15 12 21 32 29 38 57   301

#summary checks out 

##but when I call for usage or try to plot Tu, it shows that there is nothing

> usage(ResightJanApr17CH)

NULL

> sightingPlot(ResightJanApr17CH, type = "Tu", mean = FALSE, dropunused = FALSE, fill = TRUE)

#so, this indicates to me that it might be something to do with my data files or with the first two lines of code creating the caphist object

##then again, the summary looks right, so the data apparently read in fine 

###here is the model fitting

#determine buffer

> 3*RPSV(ResightJanApr17CH)

[1] 11960.2

> resight17null<- secr.fit(ResightJanApr17CH, buffer=3*RPSV(ResightJanApr17CH))

Warning messages:

1: In autoini(ch, msk, binomN = details$binomN, ignoreusage = details$ignoreusage) :

  'autoini' failed to find g0; setting initial g0 = 0.1

2: In bias.D(buffer, temptrps, detectfn = output$detectfn, detectpar = dpar,  :

  bias.D() does not allow for variable effort (detector usage)

3: In bufferbiascheck(output, buffer, biasLimit) :

  predicted relative bias exceeds 0.01 with buffer = 11960.1986316536

##maybe there is the answer in these warnings?

#buffer again

>esa.plot(resight17null)

>suggest.buffer(resight17null)

[1] 16727

1: In bias.D(w, traps, detectfn, detectpar, noccasions, binomN,  ... :

  bias.D() does not allow for variable effort (detector usage)

##this error x31

#rerun with new buffer

>resight17null<- secr.fit(ResightJanApr17CH, buffer=16727)

#same warnings as above, incl. buffer warning 

>derived(resight17null)

 estimate  SE.estimate          lcl          ucl       CVn       CVa       CVD

  esa 1.285257e+05           NA           NA           NA        NA        NA        NA

  D   1.167082e-04 4.318877e-05 5.783072e-05 0.0002355288 0.2581989 0.2650964 0.3700578

Warning messages:

1: In esa(object, sessnum) : use of allsighting here untested

2: In esa(object, sessnum, beta) : use of allsighting here untested

>verify(ResightJanApr17CH)

No errors found :-)

#adjust overdispersion

>resight17null2<- secr.fit(ResightJanApr17CH, buffer= 16727, start = resight17null, details = list (nsim=2000))

Error in integrate(integrand2, lower = trapspacing/2^0.5, upper = Inf) : 

  roundoff error is detected in the extrapolation table

In addition: Warning message:

In bias.D(buffer, temptrps, detectfn = output$detectfn, detectpar = dpar,  :

  bias.D() does not allow for variable effort (detector usage)

>predict(resight17null2, newdata = data.frame(ts = factor(c('M', 'S'))))

$` = `

       link     estimate  SE.estimate          lcl          ucl

D       log 2.348774e-04 5.713125e-05 1.468129e-04 3.757666e-04

g0    logit 9.232586e-03 3.260598e-03 4.612438e-03 1.839508e-02

sigma   log 1.868323e+04 4.283001e+03 1.198983e+04 2.911326e+04

pID   logit 4.935568e-01 9.972226e-02 3.083782e-01 6.805203e-01

$` = `

       link     estimate  SE.estimate          lcl          ucl

D       log 2.348774e-04 5.713125e-05 1.468129e-04 3.757666e-04

g0    logit 9.232586e-03 3.260598e-03 4.612438e-03 1.839508e-02

sigma   log 1.868323e+04 4.283001e+03 1.198983e+04 2.911326e+04

pID   logit 4.935568e-01 9.972226e-02 3.083782e-01 6.805203e-01

> resight17null2$details$chat

        Tu        Tm       Tn

1 4.201612 0.9694753 4.221859

>derived(resight17null2)

  estimate SE.estimate          lcl          ucl       CVn       CVa       CVD

esa 4.986201e+05          NA           NA           NA        NA        NA        NA

D   3.008302e-05 9.42118e-06 1.651805e-05 5.478784e-05 0.2581989 0.1772299 0.3131727

Warning messages:

1: In esa(object, sessnum) : use of allsighting here untested

2: In esa(object, sessnum, beta) : use of allsighting here untested

________________

Hopefully it is something simple. I suspect there is something wrong in my data files, but I've checked them a hundred times and just can't see it. Any help/ suggestions are welcome. The next step that I can think of is simulating some similar mark-resight data and running that instead. 

Other things that I plan to do next: model marking process (live box trap marking); add in telemetry HRestimates; add in habitat mask; add in an all-zero Tn file. 

Attached are the data files that I am using. 

Thanks in advance for any help that you can offer and for your time reading. 

Darcy 

University of Alberta 

[Murray Efford]

Hi Darcy

It will take me some time to get to the bottom of this. In the meantime -

How were lynx marked? Are you relying on natural pelage differences?

Usage is an attribute of a detector layout, so try usage(traps( ResightJanApr17CH))

The warnings you're getting from autoini and bias.D can be ignored.

Murray

[Darcy Doran-Myers]

Lynx were marked with unique ear tags and collars in live box traps. I have not modeled the marking process yet, but intend to do so later. I am confident in the IDs included in the marked file. 

Ah, usage(traps( ResightJanApr17CH)) works. Thanks, Murray. That file is likely ok then. 

Darcy 

--
You received this message because you are subscribed to the Google Groups "secr" group.
To unsubscribe from this group and stop receiving emails from it, send an email to secrgroup+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Murray Efford]

Hi again

I don't see a problem when I fit these models. Are you using secr 3.1.0? The estimates come out at 14 /100 km^2 (95%CI 8-23/100 km^2), which is almost too good to be true (vs your previous estimate). (Multiply animals/ha by 10000 to get animals/100km^2).  It's a mistake to apply derived(). Admittedly the warning message "use of allsighting here untested" should be stronger - an outright error.

The arguments of sightingPlot need to be tuned to get a useful plot:

sightingPlot(ResightJanApr17CH, type = "Tu", mean = FALSE,

    gridsp = 1000, border = 4000, dropunused = FALSE, fill = FALSE, scale = 200)

As you probably realise, the newdata argument in predict specifies a predictor 'ts' that is not in the particular model (all occasions are sighting occasions). The sighting-only estimate is making a strong assumption: uniform distribution of marked animals across the habitat mask (buffered area) so I don't see why you would bother with this model if you have the marking data.

Murray

On Friday, October 27, 2017 at 8:57:28 AM UTC+13, doranmye wrote:

[Darcy Doran-Myers]

Thanks, Murray. I am using secr 3.1.0. The Tu plot turns out fine now - thanks for your suggestions. 

I thought that calling derived() gave me a better estimate for D. I checked out the help after you said it was a mistake and can see now that I was wrong. I am now getting the same answer that you got: 14 lynx/100 km sq (9-22/ 100 km sq). (Did you round differently for the upper and lower CI? If so, is there a reason you rounded down for the lower and up for the upper?) 

I'm still wondering about the buffer warnings that I outlined above. I thought that as long as the buffer warning showed up I should use the fitted model to suggest a new buffer and rerun: 

> suggest.buffer(resight17null)

[1] 25966

But the warnings don't let up no matter how big of a buffer I use. I see from your estimate of D that you stopped after fitting the model the first time and did not increase the buffer like I did. That's easier, so I'll run with it, but I would like to know why? 

Thanks for helping me with this and for all your work with the secr group. 

Darcy 

--
You received this message because you are subscribed to the Google Groups "secr" group.
To unsubscribe from this group and stop receiving emails from it, send an email to secrgroup+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

Darcy Doran-Myers

University of Alberta

Biological Sciences- Ecology

dora...@ualberta.ca

[Murray Efford]

You're right: 8.6 rounds to 9 not 8. I was in a hurry.

The estimates are definitely unstable, but then the assumptions are screwy (marked animals distributed evenly across mask).

Murray