Help in fitting a multi-session SECR model in R Studio

[Lachlan McRae]

Hi everyone!

I'm hoping to get some assistance with trying to fit a multi-session SECR model in R Studio please. I am trying to work out the density of Short-eared possums (Trichosurus caninus) in different habitat types. I have four trapping sessions (session 1 and 3 in rainforest, session 2 in wet-sclerophyll forest and session 4 in dry-sclerophyll forest) and there are 7-9 occasions in each session. I have been trying to follow the "Multi-session models in secr 4.1" document and this is where I am up to:

> setwd("C:/Users/lmcra/OneDrive/Desktop/Career/Uni/Honours/Data/Combined")

> library(secr)

> PossumCH <- read.capthist('CHistory_4sessions.txt', c("session1.txt", "session2.txt", "session3.txt", "session4.txt"), detector = "multi")

Session 4 

No live releases

Duplicated or missing row names (animal ID)

> summary(PossumCH, terse = TRUE)

            1  2  3  4

Occasions   8  7  8  9

Detections 63 10 51  0

Animals    17  5 18  0

Detectors  13 13 13 14

> masks <- make.mask(traps(PossumCH), buffer = 275, nx = 250, type = 'trapbuffer')

> fit <- secr.fit(PossumCH, buffer = masks, trace = FALSE, verify = FALSE)

Error in -buffer : invalid argument to unary operator

> fit

secr.fit(capthist = PossumCH, buffer = 4 * 75, trace = FALSE)

secr 4.3.1, 14:01:49 28 Feb 2021

$`1`

Detector type      multi 

Detector number    53 

Average spacing    76.87578 m 

x-range            425346.8 449006.4 m 

y-range            6662975 6734869 m 

 Usage range by occasion

    1 2 3 4 5 6 7 8 9

min 0 0 0 0 0 0 0 0 0

max 1 1 1 1 1 1 1 1 1

$`2`

Detector type      multi 

Detector number    53 

Average spacing    76.87578 m 

x-range            425346.8 449006.4 m 

y-range            6662975 6734869 m 

 Usage range by occasion

    1 2 3 4 5 6 7 8 9

min 0 0 0 0 0 0 0 0 0

max 1 1 1 1 1 1 1 1 1

            1  2

Occasions   9  9

Detections 73 51

Animals    22 18

Detectors  53 53

Model           :  D~1 g0~1 sigma~1 

Fixed (real)    :  none 

Detection fn    :  halfnormal

Distribution    :  poisson 

N parameters    :  3 

Log likelihood  :  -461.8017 

AIC             :  929.6033 

AICc            :  930.27 

Beta parameters (coefficients) 

           beta    SE.beta       lcl       ucl

D     -2.529221 0.16012344 -2.843057 -2.215385

g0    -1.630940 0.16913232 -1.962433 -1.299446

sigma  5.179691 0.06869625  5.045049  5.314334

Variance-covariance matrix of beta parameters 

                  D            g0        sigma

D      0.0256395147 -0.0009782534 -0.001286739

g0    -0.0009782534  0.0286057404 -0.005656104

sigma -0.0012867390 -0.0056561040  0.004719174

Fitted (real) parameters evaluated at base levels of covariates 

 session = 1 

       link     estimate SE.estimate          lcl         ucl

D       log   0.07972111  0.01284748   0.05824733   0.1091115

g0    logit   0.16370168  0.02315480   0.12320400   0.2142582

sigma   log 177.62799372 12.21678703 155.25194826 203.2290384

 session = 2 

       link     estimate SE.estimate          lcl         ucl

D       log   0.07972111  0.01284748   0.05824733   0.1091115

g0    logit   0.16370168  0.02315480   0.12320400   0.2142582

sigma   log 177.62799372 12.21678703 155.25194826 203.2290384

For some reason I am only getting density estimates for two sessions, not all four? The density estimates for the two sessions are also identical which can't be right. Can someone please take a look at the for me and let me know where I have gone wrong? Can I also get some assistance with exactly what I need to type into R to assign a habitat type to each session please? I have attached the four detector layouts for each session as text files and I have also attached my capture history text file.

Any help is super appreciated!

Cheers,

Lachlan McRae

[Murray Efford]

The 'Duplicated or missing rownames' message is spurious (and won't appear in next version) - it's just saying there were no captures in that session.

'buffer' and 'mask' are alternatives in secr.fit. 'buffer' should be a number or omitted entirely, not given a list of masks. You probably meant

fit <- secr.fit(PossumCH, mask = masks, trace = FALSE, verify = FALSE)

[Lachlan McRae]

Thanks for getting back to me Murray.

I entered what you suggested but now this error message appears at the end:

> setwd("C:/Users/lmcra/OneDrive/Desktop/Career/Uni/Honours/Data/Combined")

> library(secr)

> PossumCH <- read.capthist('CHistory_4sessions.txt', c("session1.txt", "session2.txt", "session3.txt", "session4.txt"), detector = "multi")

Session 4 

No live releases

Duplicated or missing row names (animal ID)

> summary(PossumCH, terse = TRUE)

            1  2  3  4

Occasions   8  7  8  9

Detections 63 10 51  0

Animals    17  5 18  0

Detectors  13 13 13 14

> masks <- make.mask(traps(PossumCH), buffer = 275, nx = 250, type = 'trapbuffer')

> secr.fit(PossumCH, mask = masks, trace = FALSE, verify = FALSE)

Error in tmp[, 4] : subscript out of bounds

Any ideas on how to get past this error?

Cheers,

Lachlan

[Murray Efford]

It works for me. Could be a version problem - do you have latest (4.3.3)?

Also nx = 275 is way larger than it needs to be. I got almost same answer with nx = 32 in a fraction of the time.

[Lachlan McRae]

I did not have version 4.3.3. I just downloaded it and it now works! Thanks Murray!

However, the results produce the same density estimate for all four sessions which is not what I want. I went to the "Faster fitting of fully session-specific models" section of the "Multi-session models in secr 4.1" document and used the code provided there (fits <- lapply(PossumCH, secr.fit, buffer = 275, trace = FALSE)) to produce unique density estimates for each session but it has now produced another error message:

> setwd("C:/Users/lmcra/OneDrive/Desktop/Career/Uni/Honours/Data/Combined")

> library(secr)

This is secr 4.3.3. For overview type ?secr

Warning message:

package ‘secr’ was built under R version 4.0.4 

> PossumCH <- read.capthist('CHistory_4sessions.txt', c("session1.txt", "session2.txt", "session3.txt", "session4.txt"), detector = "multi")

Session 4 

No live releases

Duplicated or missing row names (animal ID)

> summary(PossumCH, terse = TRUE)

            1  2  3  4

Occasions   8  7  8  9

Detections 63 10 51  0

Animals    17  5 18  0

Detectors  13 13 13 14

> masks <- make.mask(traps(PossumCH), buffer = 275, nx = 32, type = 'trapbuffer')

> secr.fit(PossumCH, mask = masks, trace = FALSE, verify = FALSE)

secr.fit(capthist = PossumCH, mask = masks, verify = FALSE, trace = FALSE)

secr 4.3.3, 11:35:13 03 Mar 2021

$`1`

Detector type      multi 

Detector number    13 

Average spacing    77.62104 m 

x-range            448469.3 449006.4 m 

y-range            6662975 6663358 m 

 Usage range by occasion

    1 2 3 4 5 6 7 8

min 1 1 1 1 1 1 1 1

max 1 1 1 1 1 1 1 1

$`2`

Detector type      multi 

Detector number    13 

Average spacing    83.62952 m 

x-range            444601.2 445490 m 

y-range            6665718 6665971 m 

 Usage range by occasion

    1 2 3 4 5 6 7

min 1 1 1 1 1 1 1

max 1 1 1 1 1 1 1

$`3`

Detector type      multi 

Detector number    13 

Average spacing    74.26613 m 

x-range            425346.8 425663.6 m 

y-range            6730226 6730726 m 

 Usage range by occasion

    1 2 3 4 5 6 7 8

min 1 1 1 1 1 1 1 1

max 1 1 1 1 1 1 1 1

$`4`

Detector type      multi 

Detector number    14 

Average spacing    54.71638 m 

x-range            437823.5 438040 m 

y-range            6734524 6734869 m 

 Usage range by occasion

    1 2 3 4 5 6 7 8 9

min 0 0 0 0 0 0 0 0 0

max 1 1 1 1 1 1 1 1 1

            1  2  3  4

Occasions   8  7  8  9

Detections 63 10 51  0

Animals    17  5 18  0

Detectors  13 13 13 14

Model           :  D~1 g0~1 sigma~1 

Fixed (real)    :  none 

Detection fn    :  halfnormal

Distribution    :  poisson 

N parameters    :  3 

Log likelihood  :  -395.4022 

AIC             :  796.8044 

AICc            :  797.4711 

Beta parameters (coefficients) 

            beta    SE.beta       lcl        ucl

D     -1.2045132 0.19568197 -1.588043 -0.8209836

g0    -0.6374523 0.24942069 -1.126308 -0.1485967

sigma  4.3437422 0.09361131  4.160267  4.5272170

Variance-covariance matrix of beta parameters 

                 D           g0        sigma

D      0.038291432 -0.002331076 -0.009997925

g0    -0.002331076  0.062210681 -0.007014970

sigma -0.009997925 -0.007014970  0.008763077

Fitted (real) parameters evaluated at base levels of covariates 

 session = 1 

       link   estimate SE.estimate        lcl        ucl

D       log  0.2998379  0.05923905  0.2043251  0.4399987

g0    logit  0.3458227  0.05642628  0.2448431  0.4629190

sigma   log 76.9951337  7.22343437 64.0886593 92.5007744

 session = 2 

       link   estimate SE.estimate        lcl        ucl

D       log  0.2998379  0.05923905  0.2043251  0.4399987

g0    logit  0.3458227  0.05642628  0.2448431  0.4629190

sigma   log 76.9951337  7.22343437 64.0886593 92.5007744

 session = 3 

       link   estimate SE.estimate        lcl        ucl

D       log  0.2998379  0.05923905  0.2043251  0.4399987

g0    logit  0.3458227  0.05642628  0.2448431  0.4629190

sigma   log 76.9951337  7.22343437 64.0886593 92.5007744

 session = 4 

       link   estimate SE.estimate        lcl        ucl

D       log  0.2998379  0.05923905  0.2043251  0.4399987

g0    logit  0.3458227  0.05642628  0.2448431  0.4629190

sigma   log 76.9951337  7.22343437 64.0886593 92.5007744

> fits <- lapply(PossumCH, secr.fit, buffer = 275, trace = FALSE)

Session 4 

No live releases

Duplicated or missing row names (animal ID)

Error in FUN(X[[i]], ...) : 'verify' found errors in 'capthist' argument

In addition: Warning message:

In bufferbiascheck(output, buffer, biasLimit) :

  predicted relative bias exceeds 0.01 with buffer = 275

I can get passed the warning "In addition: Warning message: In bufferbiascheck(output, buffer, biasLimit) : predicted relative bias exceeds 0.01 with buffer = 275" by increasing the buffer width to 450 (even though I think this is far too excessive). I am hoping you can please help me out with the error that says "Error in FUN(X[[i]], ...) : 'verify' found errors in 'capthist' argument" ?

Cheers,

Lachlan

[Murray Efford]

You seem to be crashing around - may be time to slow down and do some homework. To model differing density among sessions call secr.fit with model = D ~ session. To avoid verify message, specify verify = FALSE as in your earlier example. Quite possible the estimate of sigma is large in one of your sessions, causing the bias warning (_warning_). And with no captures in session 4 it's silly to fit a model, so lapply(PossumCH[1:3]...) is better.

[Lachlan McRae]

Thanks for your help Murray. I apologise if the answers to my questions are obvious, I'm still learning the ins and outs of R and your guidance if very appreciated.

I have two more questions if that's okay?

Question 1) You said "To model differing density among sessions call secr.fit with model = D ~ session." however R produced a unique density estimate for each session just by adding " lapply(PossumCH[1:3]...)" like you mentioned above. See below: 

> fits <- lapply(PossumCH[1:3], secr.fit, buffer = 275, trace = FALSE)

Warning message:

In bufferbiascheck(output, buffer, biasLimit) :

  predicted relative bias exceeds 0.01 with buffer = 275

> class(fits) <- 'secrlist'

> predict(fits)

$`1`

       link   estimate SE.estimate        lcl        ucl

D       log  0.4785431   0.1358555  0.2772920  0.8258566

g0    logit  0.5126066   0.1227146  0.2865418  0.7336291

sigma   log 71.5382452   9.0716690 55.8507162 91.6321379

$`2`

       link     estimate SE.estimate         lcl

D       log   0.08781757  0.05349604  0.02920020

g0    logit   0.15430673  0.09784007  0.04028311

sigma   log 127.92925888 48.18079019 62.65560244

              ucl

D       0.2641053

g0      0.4423271

sigma 261.2040207

$`3`

       link   estimate SE.estimate        lcl        ucl

D       log  0.7315003  0.21801774  0.4129512  1.2957770

g0    logit  0.2542963  0.05903097  0.1563082  0.3856347

sigma   log 68.4717817  9.94469752 51.5857215 90.8853215

I can also get unique density estimates for each session by using model = D~session in secr.fit. See below:

> secr.fit(PossumCH, model = D~session, mask = masks, trace = FALSE, verify = FALSE)

secr.fit(capthist = PossumCH, model = D ~ session, mask = masks, 

    verify = FALSE, trace = FALSE)

secr 4.3.3, 16:38:18 04 Mar 2021

Model           :  D~session g0~1 sigma~1 

Fixed (real)    :  none 

Detection fn    :  halfnormal

Distribution    :  poisson 

N parameters    :  6 

Log likelihood  :  -379.9028 

AIC             :  771.8057 

AICc            :  774.3511 

Beta parameters (coefficients) 

                  beta      SE.beta         lcl         ucl

D           -0.7366118 2.675597e-01  -1.2610192  -0.2122045

D.session2  -1.3370696 5.088190e-01  -2.3343365  -0.3398027

D.session3   0.1746731 3.381917e-01  -0.4881705   0.8375167

D.session4 -15.7691603 4.652433e-09 -15.7691603 -15.7691603

g0          -0.6407638 2.496777e-01  -1.1301232  -0.1514045

sigma        4.3505726 9.520928e-02   4.1639658   4.5371793

Variance-covariance matrix of beta parameters 

                       D    D.session2    D.session3

D           7.158818e-02 -5.782777e-02 -5.850716e-02

D.session2 -5.782777e-02  2.588967e-01  5.883713e-02

D.session3 -5.850716e-02  5.883713e-02  1.143736e-01

D.session4 -9.667793e-10  5.289717e-10  6.744524e-10

g0         -2.066598e-03 -1.216788e-03  4.165053e-05

sigma      -9.988681e-03 -6.217177e-04 -2.538103e-04

              D.session4            g0         sigma

D          -9.667793e-10 -2.066598e-03 -9.988681e-03

D.session2  5.289717e-10 -1.216788e-03 -6.217177e-04

D.session3  6.744524e-10  4.165053e-05 -2.538103e-04

D.session4  2.164513e-17 -6.437811e-10  3.060474e-10

g0         -6.437811e-10  6.233897e-02 -7.200086e-03

sigma       3.060474e-10 -7.200086e-03  9.064807e-03

Fitted (real) parameters evaluated at base levels of covariates 

 session = 1 

       link   estimate SE.estimate        lcl        ucl

D       log  0.4787332  0.13041668  0.2833651  0.8087993

g0    logit  0.3450739  0.05642664  0.2441384  0.4622210

sigma   log 77.5228377  7.39765175 64.3261229 93.4269018

 session = 2 

       link   estimate SE.estimate         lcl        ucl

D       log  0.1257221  0.06154612  0.05068533  0.3118465

g0    logit  0.3450739  0.05642664  0.24413837  0.4622210

sigma   log 77.5228377  7.39765175 64.32612285 93.4269018

 session = 3 

       link   estimate SE.estimate        lcl        ucl

D       log  0.5701027  0.15231451  0.3407578  0.9538067

g0    logit  0.3450739  0.05642664  0.2441384  0.4622210

sigma   log 77.5228377  7.39765175 64.3261229 93.4269018

 session = 4 

       link     estimate  SE.estimate          lcl

D       log 6.786319e-08 1.848732e-08 4.016863e-08

g0    logit 3.450739e-01 5.642664e-02 2.441384e-01

sigma   log 7.752284e+01 7.397652e+00 6.432612e+01

              ucl

D     1.14652e-07

g0    4.62221e-01

sigma 9.34269e+01

There is quite a bit of difference in the density estimates in session 2 and 3 depending on which method I use to estimate the densities. I'm hoping you can please explain what method I should use and why?

Question 2)

The main goal of my research project is to determine the density of Short-eared Possums in different habitat types and work out whether habitat type is a predictor of density. I'm hypothesizing that rainforest habitat will have the highest density of possums. I'm planning to compare the AICc value of a model with habitat as a covariant to a model without the habitat covariate to determine whether habitat is a predictor of density. When trying to add the covariate "Habitat" to the mask I am getting an error which I'm hoping you can take a look at please? See below:

> setwd("C:/Users/lmcra/OneDrive/Desktop/Career/Uni/Honours/Data/Combined")

> library(secr)

> PossumCH <- read.capthist('CHistory_4sessions.txt', c("session1.txt", "session2.txt", "session3.txt", "session4.txt"), detector = "multi")

Session 4 

No live releases

Duplicated or missing row names (animal ID)

> summary(PossumCH, terse = TRUE)

            1  2  3  4

Occasions   8  7  8  9

Detections 63 10 51  0

Animals    17  5 18  0

Detectors  13 13 13 14

> masks <- make.mask(traps(PossumCH), buffer = 275, nx = 32, type = 'trapbuffer')

> covariatesource <- read.mask('HabitatCovariate.txt')

> masks <- addCovariates(masks, covariatesource)

> fitted <- secr.fit(PossumCH, mask = masks, detectfn = 'HN', 

+                    model = D~Habitat, trace = FALSE, verify = FALSE)

Error in D.designdata(mask, model$D, grouplevels, session(capthist), sessioncov) : 

  Habitat not found

It is saying that "Habitat" is not found so I went back into my HabitatCovariate text file and I removed the "#" symbol from the top row where the column headings are so that the word "Habitat" is read but then a different error message appears. See below:

> covariatesource <- read.mask('HabitatCovariate.txt')

Error in read.mask("HabitatCovariate.txt") : 

  non-numeric x or y coordinates

How do I get past this issue? I have attached my HabitatCovariate text file.

Cheers,

Lachlan

[Murray Efford]

1) These are different models - the lapply() approach ensures all parameters are estimated independently for each session. The other approach forces shared estimates of the detection parameters. And the density estimates are really not that different if you look at the CI.

2) Your 'covariatesource' has a column labelled 'V3' not 'Habitat' (this is a quirk of how read.mask works). You can fix this with

names(covariates(covariatesource)) <- 'Habitat'

Better still, recognise that habitat varies at the session level and use a session covariate instead of a mask covariate. Then

sessionhab <- data.frame(habitat = c('Rainforest', 'DrySclerophyll', 'Rainforest', 'DrySclerophyll')) # matching sessions 1-4
fitted2 <- secr.fit(PossumCH, mask = masks, detectfn = 'HN',
    model = D~habitat, sessioncov = sessionhab, trace = FALSE, verify = FALSE)
predict(fitted, newdata = data.frame(habitat = factor(c('DrySclerophyll', 'Rainforest'))))
$`habitat = DrySclerophyll`

       link    estimate SE.estimate         lcl        ucl

D       log  0.07404832   0.0362017  0.02988519  0.1834740
g0    logit  0.34630252   0.0565685  0.24506654  0.4636731
sigma   log 77.55578937   7.3919958 64.36768216 93.4459695

$`habitat = Rainforest`

       link   estimate SE.estimate        lcl        ucl

D       log  0.5209476   0.1073522  0.3493064  0.7769294
g0    logit  0.3463025   0.0565685  0.2450665  0.4636731
sigma   log 77.5557894   7.3919958 64.3676822 93.4459695

Which gives virtually identical estimates to the other way.

[Lachlan McRae]

Thanks so much Murray! 

I've now determined the AICc value for the model without Habitat as a covariate (AICc = 796.804) and the model with Habitat as a covariate (AICc =  773.35). Since the model that contains Habitat as a covariate produces a lower AICc value, can I now say that it is the most supported model and therefore habitat type is a predictor of density? Can I go a step further than that and say that not only is habitat a predictor of density, but the density of short-eared possums is higher in Rainforest habitat compared to Dry Sclerophyll? Below is my R script:

> setwd("C:/Users/lmcra/OneDrive/Desktop/Career/Uni/Honours/Data/Combined")

> library(secr)

> PossumCH <- read.capthist('CHistory_4sessions.txt', c("session1.txt", "session2.txt", "session3.txt", "session4.txt"), detector = "multi")

Session 4 

No live releases

Duplicated or missing row names (animal ID)

> summary(PossumCH, terse = TRUE)

            1  2  3  4

Occasions   8  7  8  9

Detections 63 10 51  0

Animals    17  5 18  0

Detectors  13 13 13 14

> masks <- make.mask(traps(PossumCH), buffer = 275, nx = 32, type = 'trapbuffer')

> fitted1 <- secr.fit(PossumCH, mask = masks, trace = FALSE, verify = FALSE)

> predict(fitted1)

$`session = 1`

       link   estimate SE.estimate        lcl        ucl

D       log  0.2998377  0.05923903  0.2043249  0.4399984

g0    logit  0.3458208  0.05642590  0.2448420  0.4629164

sigma   log 76.9952927  7.22346556 64.0887653 92.5010035

$`session = 2`

       link   estimate SE.estimate        lcl        ucl

D       log  0.2998377  0.05923903  0.2043249  0.4399984

g0    logit  0.3458208  0.05642590  0.2448420  0.4629164

sigma   log 76.9952927  7.22346556 64.0887653 92.5010035

$`session = 3`

       link   estimate SE.estimate        lcl        ucl

D       log  0.2998377  0.05923903  0.2043249  0.4399984

g0    logit  0.3458208  0.05642590  0.2448420  0.4629164

sigma   log 76.9952927  7.22346556 64.0887653 92.5010035

$`session = 4`

       link   estimate SE.estimate        lcl        ucl

D       log  0.2998377  0.05923903  0.2043249  0.4399984

g0    logit  0.3458208  0.05642590  0.2448420  0.4629164

sigma   log 76.9952927  7.22346556 64.0887653 92.5010035

> covariatesource <- read.mask('HabitatCovariate.txt')

> masks <- addCovariates(masks, covariatesource)

> sessionhab <- data.frame(habitat = c('Rainforest', 'DrySclerophyll', 'Rainforest', 'DrySclerophyll'))

> fitted2 <- secr.fit(PossumCH, mask = masks, detectfn = 'HN', model = D~habitat, sessioncov = sessionhab, trace = FALSE, verify = FALSE)

> predict(fitted2, newdata = data.frame(habitat = factor(c('DrySclerophyll', 'Rainforest')))) 

$`habitat = DrySclerophyll`

       link    estimate SE.estimate         lcl        ucl

D       log  0.07404788  0.03620163  0.02988491  0.1834735

g0    logit  0.34630334  0.05656865  0.24506708  0.4636742

sigma   log 77.55572360  7.39198307 64.36763799 93.4458752

$`habitat = Rainforest`

       link   estimate SE.estimate        lcl        ucl

D       log  0.5209477  0.10735225  0.3493065  0.7769295

g0    logit  0.3463033  0.05656865  0.2450671  0.4636742

sigma   log 77.5557236  7.39198307 64.3676380 93.4458752

> AIC(fitted1)[,3:5]

        npar    logLik     AIC

fitted1    3 -395.4022 796.804

> AIC(fitted2)[,3:5]

        npar    logLik    AIC

fitted2    4 -382.6751 773.35

Cheers,

Lachlan