IELNE on telemetry+census mark-resight design

[cpenaherrera]

Dear all,

Thanks for allowing to be part of the phidot group. I have some technical issues with the design and model application I can't overcome and I would like to kindly request your expert advise.

My project is bridging acoustic telemetry and visual census for assessing aggregation size of scalloped hammerhead sharks in an oceanic island of the eastern tropical Pacific ocean. For this, I choose the season with the highest abundance on the year (September) and repeated my experiment during three years. In each yearly session, I tagged 10 to 30 sharks in the first two days by free diving, and then carried out visual census three times per day during five to eight days.

I assembled my data by using the counts of hammerheads (from visual census) and then pairing the resights with the amount of sharks detected in the acoustic receivers in the same site where census were carried out. Each census lasted only 20 minutes (time estimated from continuous tracking of the species in a contiguous site and that relates the amount of time sharks needed for coming back again to the same site).

For assessing this we choose (With my advisor) the Inmigration-Emigration Logit-Normal Estimator. We are confident that my design satisfies the models assumptions, but I'm assessing this in a yearly case as my sharks did not return to the same site in the following years, or even they did not return at all (We have many acoustic receivers placed around the island, and sister Islands to check their presence). We built the matrices with the idea that census in every day are secondary occasions, and days are primary ones, as follows (provided a dummy example):

Code: Select all

000101000100000 1;
...011110001111 1;
000111101000000 1;
000000...000000 1;
000000000010000 1;
100000100000000 1;
......000...000 1;
000001......011 1;
010000100100010 1;

I have run this model using mark and Rmark separately (I built my example code for R using a Poisson example posted here by Jesse Lewis and I adapted it to my IELNE model). Code as follows:

Code: Select all

# IELNE Example
# Created by Cesar Penaherrera 2014, using as basis Poisson exmaple from Jesse Lewis, Jeff Laake
# and Brett McClintock.
# Power computer version


# clear all list from system
# rm(list=ls())

# activate RMark
library(RMark)


# Set working directory and add the library of RMark
setwd("P:/@ Rmark/IELNE Example")

# import inp file. This data comes from IELNE Example 2011, counted individuals. 
IELNExample=convert.inp("P:/@ Rmark/IELNE Example/IELNExample.inp")
IELNExample

# calculates number of animals in input file
n = length(IELNExample[,1])
n

# Specify the desing of the study (create *.proc and *.ddl objects). This example is a
# robust design in which there are three secondary occassions per primary one. There are
# five primary ocassions, totalling 15 secondary ones. Time interval reflects this but omits
# the first seconday occassion of the first primary one.
IELNExample.proc=process.data(IELNExample,model="IELogitNormalMR",
   counts=list("Marked Superpopulation"=c(9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9),
      "Unmarked Seen"=c(129, 66, 26, 119, 59, 79, 41, 46, 28, 27, 3, 78, 54, 312, 138),
      "Marked Unidentified"=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)),
   time.intervals=c(0,0,1,0,0,1,0,0,1,0,0,1,0,0))

IELNExample.ddl=make.design.data(IELNExample.proc)

## Run Multiple Models
# FIRST: estimate starting values from alpha(.) N(.) sig(.) model. SIMANNEAL means Simulated
# annealing. This resembles Alt opt method from Mark GUI --> sidebar, page 39, chapter 10 Mark book.
all.dot = mark(IELNExample.proc, IELNExample.ddl,
   model.parameters=list(p=list(formula=~1, link='logit'),
      sigma = list(formula = ~1, link='logit'),
      Nbar = list(formula = ~1, link='log'),
      alpha = list(formula = ~1, link='log'),
      Nstar = list(formula = ~1, link='log')),
    threads=16,  options = "SIMANNEAL",filename= 'alldot3s') #options = "SIMANNEAL", Using this option considerably increases the processing time "CONSIDERABLY"
y

all.dot$results$AICc
all.dot$results$real # $AICc [1] 907.8979 $AICc.unadjusted [1] 905.7327
#                       estimate      se             lcl          ucl      fixed
# p g1 s1 t1            0.282618  2.324560e-02    0.2393522    0.3303086
# sigma g1 a0 s1 t0     1.000000  9.760711e-19    1.0000000    1.0000000
# Nbar g1 a0 s1 t0    276.838130  2.034415e+01  242.2663300  322.6572300
# alpha g1 a0 s1 t1    18.019517  4.580464e+00    9.0418077   26.9972260
# Nstar g1 a0 s1 t0   335.804830  1.302565e+01  315.4364800  367.5089400

# Set this values of alpha, sigma and N as initial start values for alpha, U, sig for multiple models.
# Initial parameters must be provided in the log scale (pag 36, chapter 18, Mark book)
p = 0.282618
alpha = 18.019517
sigma = log(1.00000)
Nbar = log(276.838130)
Nstar = log(335.804830)

##### Running different models of interest
##  all.variable
all.var = mark(IELNExample.proc, IELNExample.ddl,
   model.parameters=list(p=list(formula=~session:time, link='logit'),
      sigma = list(formula = ~session, link='log'),
      Nbar = list(formula = ~session, link='log'),
      alpha = list(formula = ~session:time, link='logit'),
      Nstar = list(formula = ~session, link='log')),
   initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
    threads=16, options = "SIMANNEAL", filename='allvariable')
all.var$results$real

all.varLogit = mark(IELNExample.proc, IELNExample.ddl,
                      model.parameters=list(p=list(formula=~session:time, link='logit'),
                                                       sigma = list(formula = ~session, link='logit'),
                                                       Nbar = list(formula = ~session, link='logit'),
                                                       alpha = list(formula = ~session:time, link='logit'),
                                                       Nstar = list(formula = ~session, link='logit')),
                      initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                      threads=18, options = "SIMANNEAL", filename='allvariableLogit')
all.varLogit$results$real


##  allvar.sigmafixLogit
allvar.sigmafix = mark(IELNExample.proc, IELNExample.ddl,
   model.parameters=list(p=list(formula=~session:time, link='logit'),
      sigma = list(formula = ~1, fixed = 0, link='logit'),
      Nbar = list(formula = ~session, link='logit'),
      alpha = list(formula = ~session:time, link='logit'),
      Nstar = list(formula = ~session, link='logit')),
   initial = c(Nbar = Nbar, Nstar=Nstar, sigma = 0),
   threads=18, options = "SIMANNEAL", filename='allvar.sigmafixLogit')
#allvar.sigmafixLogit$results$real
#allvar.sigmafixLogit$results$AICc

## alpha.dot
#alpha.dot = mark(IELNExample.proc, IELNExample.ddl,
#   model.parameters=list(p=list(formula=~session:time),
#      sigma = list(formula = ~session),
#      Nbar = list(formula = ~session),
#      alpha = list(formula = ~1),
#      Nstar = list(formula = ~session)),
#   initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
#   threads=16, filename='alphadot')
#alpha.dot$results$real

## alpha.dotLogit
alpha.dotLogit = mark(IELNExample.proc, IELNExample.ddl,
                         model.parameters=list(p=list(formula=~session:time, link='logit'),
                                                          sigma = list(formula = ~session, link='logit'),
                                                          Nbar = list(formula = ~session, link='logit'),
                                                          alpha = list(formula = ~1, link='logit'),
                                                          Nstar = list(formula = ~session, link='logit')),
                         initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                         threads=16, filename='alphadotLogit')
#alpha.dot$results$real
#alpha.dot$results$AICc

## alpha.sessionLogit
alpha.sessionLogit = mark(IELNExample.proc, IELNExample.ddl,
                                       model.parameters=list(p=list(formula=~session:time, link='logit'),
                                                                        sigma = list(formula = ~session, link='logit'),
                                                                        Nbar = list(formula = ~session, link='logit'),
                                                                        alpha = list(formula = ~session, link='logit'),
                                                                        Nstar = list(formula = ~session, link='logit')),
                                       initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                                       threads=16,  options = "SIMANNEAL", filename='alphasessionLogit')
alpha.session$results$real

## Nbar.dot
Nbar.dot = mark(IELNExample.proc, IELNExample.ddl,
                        model.parameters=list(p=list(formula=~session:time, link='logit'),
                                                         sigma = list(formula = ~session, link='logit'),
                                                         Nbar = list(formula = ~1, link='logit'),
                                                         alpha = list(formula = ~session:time, link='logit'),
                                                         Nstar = list(formula = ~session, link='logit')),
                        initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                        threads=16, options = "SIMANNEAL", filename='Nbardot')

## Nstar.dot
Nstar.dot = mark(IELNExample.proc, IELNExample.ddl,
                         model.parameters=list(p=list(formula=~session:time, link='logit'),
                                                          sigma = list(formula = ~session, link='logit'),
                                                          Nbar = list(formula = ~session, link='logit'),
                                                          alpha = list(formula = ~session:time, link='logit'),
                                                          Nstar = list(formula = ~1, link='logit')),
                         initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                         threads=16, options = "SIMANNEAL", filename='Nstardot')
Nstar.dot$results$AICc
Nstar.dot$results$real

## psigma.dot
psigma.dot = mark(IELNExample.proc, IELNExample.ddl,
                           model.parameters=list(p=list(formula=~1, link='logit'),
                                                            sigma = list(formula = ~1, link='logit'),
                                                            Nbar = list(formula = ~session, link='logit'),
                                                            alpha = list(formula = ~session:time, link='logit'),
                                                            Nstar = list(formula = ~session, link='logit')),
                           initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                           threads=16, options = "SIMANNEAL", filename='psigmadot')
#psigma.dot$results$AICc
## psigmaNstar.dot
psigmaNstar.dot = mark(IELNExample.proc, IELNExample.ddl,
                                  model.parameters=list(p=list(formula=~1, link='logit'),
                                                                   sigma = list(formula = ~1, link='logit'),
                                                                   Nbar = list(formula = ~session, link='log'),
                                                                   alpha = list(formula = ~session:time, link='logit'),
                                                                   Nstar = list(formula = ~1, link='log')),
                                  initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma, link='log'),
                                  threads=16, options = "SIMANNEAL", filename='psigmaNstardot')

## what if the model I want is fixed in N but varies p
## NbarNstar.dot

NbarStar.dot = mark(IELNExample.proc, IELNExample.ddl,
                              model.parameters=list(p=list(formula=~session:time),
                                                               sigma = list(formula = ~session),
                                                               Nbar = list(formula = ~1),
                                                               alpha = list(formula = ~session:time),
                                                               Nstar = list(formula = ~1)),
                              initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                              threads=16, options = "SIMANNEAL", filename='Nbarstardot')

## or if the Nbar varies per session but Nstar is fixed
## NbarNstar.dot

NbarStar.dot = mark(IELNExample.proc, IELNExample.ddl,
                              model.parameters=list(p=list(formula=~session),
                                                               sigma = list(formula = ~session),
                                                               Nbar = list(formula = ~session),
                                                               alpha = list(formula = ~session),
                                                               Nstar = list(formula = ~1)),
                              initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                              threads=16, options = "SIMANNEAL", filename='Nbarstardot')

## all.session but not time.

all.session = mark(IELNExample.proc, IELNExample.ddl,
                            model.parameters=list(p=list(formula=~session, link='logit'),
                                                             sigma = list(formula = ~session, link='log'),
                                                             Nbar = list(formula = ~session, link='log'),
                                                             alpha = list(formula = ~session, link='log'),
                                                             Nstar = list(formula = ~session, link='log')),
                            initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                            threads=16,
                            options = "SIMANNEAL",
                            filename='allsession')

## p.dot
p.dot = mark(IELNExample.proc, IELNExample.ddl,
                   model.parameters=list(p=list(formula=~1, link='logit'),
                                                    sigma = list(formula = ~session, link='logit'),
                                                    Nbar = list(formula = ~session, link='logit'),
                                                    alpha = list(formula = ~session:time, link='logit'),
                                                    Nstar = list(formula = ~session, link='logit')),
                   initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                   threads=16, options = "SIMANNEAL", filename='pdot')

## p.session
p.dot = mark(IELNExample.proc, IELNExample.ddl,
                   model.parameters=list(p=list(formula=~session, link='logit'),
                                                    sigma = list(formula = ~session, link='logit'),
                                                    Nbar = list(formula = ~session, link='logit'),
                                                    alpha = list(formula = ~session:time, link='logit'),
                                                    Nstar = list(formula = ~session, link='logit')),
                   initial = c(Nbar = Nbar, Nstar=Nstar, sigma = sigma),
                   threads=16, options = "SIMANNEAL", filename='pdot')

## Sigma.dot
sigma.dot = mark(IELNExample.proc, IELNExample.ddl,
                         model.parameters=list(p=list(formula=~session:time),
                                                          sigma = list(formula = ~1),
                                                          Nbar = list(formula = ~session),
                                                          alpha = list(formula = ~session:time),
                                                          Nstar = list(formula = ~session)),
                         initial = c(Nbar = Nbar,Nstar=Nstar,  sigma = sigma),
                         threads=16, options = "SIMANNEAL", filename='sigmadot')

### collect the results of all models from above
multiple.model.results = collect.models()
multiple.model.results$model.table=model.table(multiple.model.results,model.name=F)

### call results to see summary table, but first creates wider window to display results so
### that they do not wrap in R console
options(width = 200)
multiple.model.results

The first eight best models this example is at is follows:

#____model______________npar_____AICc_____DeltaAICc____weight_____Deviance
11___Nstar.dot____________38_____ 334.9172___0___________0.6631_____225.9839
5____allvar.sigmafix________37_____336.9673___2.050076____0.2379_____232.0662
8____alpha.sessionLogit_____36_____338.8858___3.968598____0.0912_____237.9293
4____all.varLogit___________42_____345.4997___10.582527___0.0033_____219.4997
9____Nbar.dot_____________38_____ 346.0984___11.18125____0.0025_____237.1651
3____all.var_______________42_____346.5746___11.657447___0.0020_____220.5746
15___sigma.dot ____________38_____352.0299___17.11271____0.0001_____243.0966

And the outcome of the best model is this:

> Nstar.dot$results$real
Label______________Estimate_______SE________LCI___________UCI
p_g1 s1 t1___________0.3342______0.1177_____0.1511_______0.5859
p_g1 s1 t2___________0.1740______0.0655_____0.0794_______0.3399
p_g1 s1 t3___________0.0666______0.0271_____0.0294_______0.1437
p_g1 s2 t1___________0.3066______0.2277_____0.0514_______0.7831
p_g1 s2 t2___________0.0909______0.1043_____0.0084_______0.5425
p_g1 s2 t3___________0.2050______0.1994_____0.0229_______0.7393
p_g1 s3 t1___________0.3498______0.2072_____0.0827_______0.7624
p_g1 s3 t2___________0.2669______0.1444_____0.0789_______0.6074
p_g1 s3 t3___________0.1540______0.1125_____0.0324_______0.4973
p_g1 s4 t1___________0.1532______0.0654_____0.0631_______0.3271
p_g1 s4 t2___________0.0183______0.0126_____0.0047_______0.0690
p_g1 s4 t3___________0.3825______0.1273_____0.1773_______0.6405
p_g1 s5 t1___________0.0000______0.0002_____0.0000_______0.9989
p_g1 s5 t2___________0.0251______0.1662_____0.0000_______0.9999
p_g1 s5 t3___________0.0011______0.0068_____0.0000_______0.9941
sigma_g1 a0 s1 t0_____0.0043______0.0225_____0.0001_______0.1547
sigma_g1 a0 s2 t0_____2.1711______1.4041_____0.6816_______6.9159
sigma_g1 a0 s3 t0_____1.0218______1.4164_____0.1343_______7.7756
sigma_g1 a0 s4 t0_____0.0000______0.0000_____0.0000_______0.0000
sigma_g1 a0 s5 t0_____6.3018______5.8937_____1.3325______29.8031
Nbar_g1 a0 s1 t0_____385.1947___133.2195___229.2734_____801.7629
Nbar_g1 a0 s2 t0_____292.0240___107.0800___178.6371_____654.3899
Nbar_g1 a0 s3 t0_____140.7404___49.5199_____84.3027_____298.9328
Nbar_g1 a0 s4 t0_____194.1692___65.7387____121.6960_____409.7739
Nbar_g1 a0 s5 t0____1021.9470__1113.4850___398.6975____6644.5960
alpha_g1 a0 s1 t1______3.3009____28.5526____-52.6621______59.2639
alpha_g1 a0 s1 t2_____-2.1959___30.1132______-61.2178______56.8261
alpha_g1 a0 s2 t1_____19.7256___29.7372______-38.5594______78.0106
alpha_g1 a0 s2 t2_____5.2427____35.6308______-64.5936______75.0790
alpha_g1 a0 s3 t1_____-20.2744__27.0576______-73.3073______32.7585
alpha_g1 a0 s3 t2_____13.3254___21.7585______-29.3214______55.9721
alpha_g1 a0 s4 t1_____-7.1110___28.4978______-62.9666______48.7446
alpha_g1 a0 s4 t2_____-4.9420___28.4393______-60.6830______50.7990
alpha_g1 a0 s5 t1_____5.7363____263.1271____-509.9928_____521.4654
alpha_g1 a0 s5 t2_____19.7256___242.8864____-456.3318_____495.7829
Nstar_g1 a0 s1 t0_____438.5744__17.2072______414.3128_____484.5421

As you can see, I used SIMANEAL option as numerical convergence was not been reached, but the computing time went to the sky, especially for the years when there are more than 20 tags available. As I'm super newbie on this, I've been reading and reviewing, and re-running my models in several ways using different link functions, yet my model outputs keep giving me results with alphas (and their corresponding s.e.) way below cero, or p/sigma equalling cero.

My main theory is that the natural variability in the abundance of sharks is so high that is messing up the models. I used GAMs to predict the values based on all the census conditions (divers, environmental conditions) and they are working (with the predicted data) a bit better, yet the problem with real parameters estimation is still there. I search the way to incorpore my environmental data into the models but unfortunately I was not able to do it. I know Mark uses a linear approach and I wonder if the variability of the data might be influencing in the estimation of the real parameters.

Said this, I would like to kindly request your help on any advise that could help to fine tune the models in order to improve real parameters estimation, or even if I need to change the model approach, any advise will be greatly appreciated. My main concerns are:

1. I understand log link function is used in Mark (GUI), but results were improved by setting all parameters with Logit link. Am I violating any assumption by changing the link for Nbar, Nstar and sigma?
2. Can I use the predicted data on abundance to correct the variability in shark abundances? or, is it better to incorporate census and environmental covariates as part of Mark analysis in order to reduce the variability?
3. Perhaps somebody may rise that IELNE could not be the most appropriate for this type of data set. Despite my efforts on satisfying model assumptions it could also be true. Does anybody foresee a problem of using IELNE with this type of data set or in the design?

Many, many thanks for the help given on this matter!

Best regards,
Cesar Penaherrera

[bmcclintock]

Hi Cesar,

Sounds like you might be jumping too quickly into a more advanced MARK analysis. Based on the information provided, I can't say much on the appropriateness of these data for the model. However, a quick glance at your code and questions brought up a few points that I hope may help you out:

1) The link functions should not be played with. I suggest leaving these to the defaults (p=logit,sigma=log,alpha=identity,Nbar=log,Nstar=log).
2) Use the beta estimates for initial values (not log-transformed real estimates; these are not equivalent)
3) Models with alpha(.) don't make any biological sense. For most systems, \alpha should be time- and session-dependent and both Nbar and Nstar should be session-dependent. To test whether the population within the study area is relatively constant, fix the \alpha's to zero. To test if the population is closed to temporary emigration (i.e., Nbar=Nstar), the \alpha's and Nstar should be fixed to zero.
4) Negative \alpha's are not a problem. In fact, some \alpha's MUST be negative (note that \sum_i \alpha_i = 0; see documentation for the model).
5) \sigma estimates close to zero can indeed be indicative of convergence issues, but they can also occur when there is very little evidence of individual heterogeneity. In the latter case, fix \sigma to zero.

I assume you've closely read chapters 1-8 of the MARK book. Besides Chapter 18 of the Book, it might be helpful to closely examine the following:

Neal, A. K., G. C. White, R. B. Gill, D. F. Reed, and J. H. Olterman. 1993. Evaluation of mark–resight model assumptions for estimating mountain sheep numbers. Journal of Wildlife Management 57:436–450.

McClintock, B.T. & White, G.C. 2012. From NOREMARK to MARK: software for estimating demographic parameters using mark-resight methodology. Journal of Ornithology, 152(Suppl. 2): S641–S650.

Once you've got a good understanding of the parameters, I'd suggest re-starting your analysis in MARK (not using Rmark for now) by running model p(.)sigma(0),alpha(0),Nbar(session)Nstar(0), where (0) indicates the parameters are fixed to 0. If this converges, build up models from there using simulated annealing as needed.

Hope this helps -- good luck!
Brett

[cpenaherrera]

Dear Dr. McClintock,

I really appreciate the time spent to answer my questions . I'm using Rmark cuz the interface of R is more familiar to me and it help me to re-run my models in a very timely and easy way, but I did used Mark-gui for a considerably amount of time. Also, I have read the cites you shared (thanks!) but I must say I have had my moments of confusion with Mark book (like changing link functions for improving model outputs).

I will re-run all the analysis using your suggestions, but I would like to please help me with some quick questions:

1 If I understood well, your suggestions is to leave covariates out of the analysis for the moment. But if not how could I add them into my matrices? I read the chapter on Mark book where it says about this and the only examples it uses are covariates that relate directly to the recaptured individuals: weight, length, etc. The data I have is environmental conditions at the time of sampling: sea temperature, thermocline presence and depth, lunar phase, categorical current strength, visibility, etc.

2 How will impact setting sigma to zero if I know that heterogeneity will tell the model that I identified all my tags (that is the interesting thing from acoustic telemetry )

3 in your last sentence:

Once you've got a good understanding of the parameters, I'd suggest re-starting your analysis in MARK (not using Rmark for now) by running model p(.)sigma(0),alpha(0),Nbar(session)Nstar(0), where (0) indicates the parameters are fixed to 0. If this converges, build up models from there using simulated annealing as needed.

Sorry for the silly questions, but did you mean by "build up from there" that I should use the beta estimates from these model to build up other models?

Many, many thanks for your help on this!!!

Cesar.

[cooch]

... If I understood well, your suggestions is to leave covariates out of the analysis for the moment. But if not how could I add them into my matrices? I read the chapter on Mark book where it says about this and the only examples it uses are covariates that relate directly to the recaptured individuals: weight, length, etc. The data I have is environmental conditions at the time of sampling: sea temperature, thermocline presence and depth, lunar phase, categorical current strength, visibility, etc. ..

environmental covariates -- Chapter 6 (in any number of places, but especially starting with section 6.8)

individual covariates -- Chapter 11

[cpenaherrera]

environmental covariates -- Chapter 6 (in any number of places, but especially starting with section 6.

individual covariates -- Chapter 11

Thanks Dr. Cooch! I did read that part, and got confused in the data input inside the matrix. I'm rechecking all once again. Many thanks for your time and help . Models are been re-run at the moment but without covariates. Will post results once I got them again.

Cheers,

Cesar P.