SEEC Toolbox on estimating survival from CMR data
Res Altwegg
28/02/2020
In this Toolbox Seminar, we are looking at how to estimate survival from capture-mark-recapture data (watch the seminar here). The survey design involved capturing a number of individuals, marking them so that they become uniquely identifyable and then release them again (unharmed!). At a later point in time, one captures another sample of individuals from the same population and checks them for marks. Unmarked individuals are marked and then all individuals (previously marked and newly marked) are released again. We do that a few times (at least 3 occasions are needed to get a survival estimate), each time leaving a time gap that corresponds to the interval over which we want to estimate survival. For example for annual survival, we would do one survey each year, ideally at about the same time (but the models can deal with unequal time intervals).
With this design, we end up with a capture history – the information on whether an individual was captured or not at each occasion – for each individual. So for example, one individual might have the capture history 1101000. That means, this individual was marked on the first occasion, recaptured on the second occasion, then not seen on the third occasion but recaptured again on the forth and then never seen again. If your study is long enough, all individuals will eventually disappear but you won’t know when each of them died. All you know is when it was last seen. That information gives us an indication of the individual’s life span but it is an underestimate. The challenge is therefore to estimate how long an individual lived after you last saw it. And that estimate depends on the recapture probability: always manage to capture all individuals that are alive and present in the population, then not seeing an individual means that it has died (or is no longer present in the population). On the other hand, if your recapture probability is low, individuals can go undetected for quite a long time. So estimating the recapture probability is key to obtaining an unbiased survival estimate.
The capture histories actually do contain information on the recapture probability: it is in these `internal’ zeros, i.e. when we miss an individual but recapture it again later. We then know that the individual must always have been alive but we failed to capture it. Getting capture histories with lot of internal zeros means that our recapture probability is relatively low, compare to the case where we rarely have internal zeros in our capture histories.
Capture-mark-recapture methods are designed to estimate survival while accounting for the recapture probability. Here, I will use a data set on hadedas (Bostrychia hagedash) to illustrate these methods and how they are implemented in program MARK. Click here for the slides and find the data for ringing here and resightings here. We will use MARK through the R package RMark. You can get program MARK from here. Make sure to read the detailed and excellent ‘Gentle Introduction to Program MARK’, also known as the ‘MARK book’. You also need to R package RMark, which you can install from CRAN. You need both to run the code below.
install.packages(RMark)
Data wrangling
The data are from a project that Doug Harebottle and I started in 2006. We went to every hadeda nest that we could find around Cape Town and ringed the nestlings (lots of tree climbing). We ringed them with a numbered SAFRING metal ring and a colour ring with a unique 2-character code. The colour ring can be read from a distance, e.g. with binoculars but also off photos, and sometimes the birds get close enough that their rings can be read by eye. So all our ‘recaptures’ are actually ‘re-sightings’. And all our birds were ringed as nestlings.
First load a few R packages that we are going to need:
library(R2ucare)
library(RMark)
library(tidyverse)
library(readxl)
library(lubridate)
Exploring the data
Our data are in two separate Excel tables: one called “ringingMay2017.xlsx” holds all the ringing data, i.e. the data we recorded each time we ringed a new individual; the other is called “resightingsMay2017.xlsx” and holds all the resightings that we collected or that were reported to us.
(As an aside, I don’t recommend storing your important data in Excel. Excel can do all sorts of damage to your data sets. But here we go…)
We first read the data in, check that evreything looks fine and then extract the year from the date variable.
ringing<-read_excel("ringingMay2017.xlsx") # table with resightings
ringing$type<-"r"
head(ringing); tail(ringing)
## # A tibble: 6 x 22
## Ring_No cring `Colour of Colo… Capture_Date nest brood Age Sex
## <chr> <chr> <lgl> <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 878509 - NA 2008-12-11 00:00:00 66 299 1 0
## 2 871551 AA NA 2006-08-31 00:00:00 1 1 1 0
## 3 871552 AB NA 2006-08-31 00:00:00 1 1 1 0
## 4 871553 AC NA 2006-09-13 00:00:00 2 4 1 0
## 5 871554 AD NA 2006-09-13 00:00:00 2 4 1 0
## 6 871555 AF NA 2006-09-26 00:00:00 4 19 1 0
## # … with 14 more variables: Mass <dbl>, Wing <dbl>, Tail <dbl>, Tarsus <dbl>,
## # Culmen <dbl>, Head <dbl>, `Parent Ring ID` <chr>, Notes <chr>,
## # Canada_ring_leg <chr>, blood <chr>, number <dbl>, F1 <chr>, ringer <chr>,
## # type <chr>
## # A tibble: 6 x 22
## Ring_No cring `Colour of Colo… Capture_Date nest brood Age Sex
## <chr> <chr> <lgl> <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 878581 VN NA 2009-11-02 00:00:00 187 384 1 0
## 2 878570 VP NA 2009-10-21 00:00:00 NA NA 1 0
## 3 878572 VS NA 2009-11-13 00:00:00 188 386 1 0
## 4 878573 VT NA 2009-11-25 00:00:00 189 388 1 0
## 5 878574 VU NA 2009-11-25 00:00:00 121 387 1 0
## 6 878575 VV NA 2009-11-25 00:00:00 121 387 1 0
## # … with 14 more variables: Mass <dbl>, Wing <dbl>, Tail <dbl>, Tarsus <dbl>,
## # Culmen <dbl>, Head <dbl>, `Parent Ring ID` <chr>, Notes <chr>,
## # Canada_ring_leg <chr>, blood <chr>, number <dbl>, F1 <chr>, ringer <chr>,
## # type <chr>
ringing$year <- as.numeric(year(ringing$Capture_Date))
resight<-read_excel("resightingsMay2017.xlsx") # table with resightings
head(resight); tail(resight)
## # A tibble: 6 x 20
## ID date time locality south east type
## <dbl> <dttm> <dttm> <chr> <chr> <chr> <chr>
## 1 5.62e7 2017-04-27 00:00:00 NA 61 Ring… <NA> <NA> l
## 2 5.62e7 2017-04-25 00:00:00 1899-12-31 08:00:00 Constan… <NA> <NA> l
## 3 5.62e7 2017-04-25 00:00:00 1899-12-31 08:00:00 Constan… <NA> <NA> l
## 4 5.62e7 2017-04-25 00:00:00 1899-12-31 08:00:00 Constan… <NA> <NA> l
## 5 5.62e7 2017-03-17 00:00:00 1899-12-31 16:00:00 lower C… <NA> <NA> l
## 6 5.62e7 2017-03-05 00:00:00 NA 61 Ring… <NA> <NA> l
## # … with 13 more variables: cring <chr>, comments <chr>, observer <chr>,
## # Field1 <chr>, Field2 <chr>, Tarsus <lgl>, Culmen <lgl>, Head <lgl>, `Parent
## # Ring ID` <lgl>, Notes <chr>, Canada_ring_leg <lgl>, blood <lgl>,
## # number <lgl>
## # A tibble: 6 x 20
## ID date time locality south east type cring
## <dbl> <dttm> <dttm> <chr> <chr> <chr> <chr> <chr>
## 1 7 2006-10-07 00:00:00 NA 2 Tanja… 3400… 1825… d AH
## 2 5 2006-09-30 00:00:00 NA 2 Tanja… 3400… 1825… d AF
## 3 4 2006-09-22 00:00:00 NA Die Oog… 3402… 1826… l AB
## 4 3 2006-09-22 00:00:00 NA Die Oog… 3402… 1826… l AA
## 5 1 2006-09-11 00:00:00 NA Die Oog… 3402… 1826… l AA
## 6 2 2006-09-11 00:00:00 NA Die Oog… 3402… 1826… l AB
## # … with 12 more variables: comments <chr>, observer <chr>, Field1 <chr>,
## # Field2 <chr>, Tarsus <lgl>, Culmen <lgl>, Head <lgl>, `Parent Ring
## # ID` <lgl>, Notes <chr>, Canada_ring_leg <lgl>, blood <lgl>, number <lgl>
resight$year <- as.numeric(year(resight$date))
Most birds were seen alive but we had a few cases where hadedas died or were killed (mostly by cars) and then someone found the body and reported the ring, called ‘dead recoveries’. This info is stored in the variable ‘type’ (‘l’ for live and ‘d’ for dead). How many observations do we have of each type?
table(resight$type)
##
## * d l
## 1 30 1418
So 1418 resightings of live birds (some individuals were seen many times) and 30 individuals were found dead. There was also one odd case where one of our birds was found living at a rehab place (’*’).
Clearly, the dead recoveries contain important information on survival and these data could be used in the analyses. But for now, we will concentrate on the live resighting data and ignore dead recoveries.
resight.l<-resight[resight$type=="l",] # selects live resightings (as opposed to dead recoveries)
When were these resightings made?
table(resight.l$year)
##
## 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017
## 43 78 154 451 282 134 102 77 45 31 11 10
We note that most resightings were made up to 2013 and relatively fewer after that. This gives us a rough idea of the effort but not of how many individuals were seen.
table(ringing$year)
##
## 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016
## 27 65 47 75 26 3 10 4 1 2 1
Tabulating the number of birds ringed each year shows that we did most of our effort between 2006 and 2010 with a few more birds ringed up to 2016.
So, even though the data include observations up to 2017, we have relatively little data to estimate survival in the later years. Also, we need to keep in mind that we have very few new birds added after 2010 and so little information on survival of young birds after that.
Creating capture histories
Next, we need to turn these data into capture histories, which is the format we will need for the analysis. We want annual survival and want to record for each individual whether it was seen in a particular year or not.
We first put the two data sets together. We only really need the ring number (or the code on the colour ring ‘cring’) and the year in which each sighting was made. So we only retain those two columns of each data frame, row-bind them and then check that all looks good.
data.live<-rbind(ringing[,c("cring","year")],resight.l[,c("cring","year")])
head(data.live); tail(data.live)
## # A tibble: 6 x 2
## cring year
## <chr> <dbl>
## 1 - 2008
## 2 AA 2006
## 3 AB 2006
## 4 AC 2006
## 5 AD 2006
## 6 AF 2006
## # A tibble: 6 x 2
## cring year
## <chr> <dbl>
## 1 AJ 2006
## 2 AI 2006
## 3 AB 2006
## 4 AA 2006
## 5 AA 2006
## 6 AB 2006
A quick way of getting something that is close to a capture history is tabulating the number of observations we have for each individual by year. Then, we just turn all numbers greater than 1 into 1, i.e. we ignore multiple resightings of the same individual in a specific year and just record that the individual was seen that year.
ch <- table(data.live$cring,data.live$year) # tabluate number of encounters per individual and year (this is close to the capture histories we need)
ch[ch>1] <- 1 # turn multiple resightings into '1'
Here we’ve got our capture history:
head(ch)
##
## 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017
## - 0 0 1 0 0 0 0 0 0 0 0 0
## AA 1 0 0 0 0 0 0 0 0 0 0 0
## AB 1 0 0 0 0 0 0 0 0 0 0 0
## AC 1 0 0 0 0 0 0 0 0 0 0 0
## AD 1 0 0 0 0 0 0 0 0 0 0 0
## AF 1 0 0 0 0 0 0 0 0 0 0 0
dim(ch)
## [1] 261 12
So we have data on 261 individuals for 12 years.
Input file for using MARK as standalone program
MARK has a user-friendly graphical user interface and I recommend that you use MARK via this interface until you are familiar with its main features. The downside of menu-driven software is of course reproducibility. Your final analysis should be code-driven, hence RMark.
To use MARK via the graphical user interface, you need to have your data in a certain format (see MARK book). The code below creates the .inp file you need:
write.table(cbind(ch," ",1,";"),"hadeda.inp",sep="",row.names=F,col.names=F,quote=F)
The MARK GUI likes to live in a Windows environment. You can get it to work on other platforms. It works quite ok for me on Ubuntu using wine.
However, some tweaks are sometimes needed. E.g. I had to open the file created with the above command in a text editor and re-save, making sure it uses UTF-8 encoding and the line ending says ‘Windows’ (not Linux/Unix).
Analysis using RMark
From here on, we’ll use RMark to communicate with MARK. The main advantage is that you will have reproducible code and therefore a much better chance to redo your exact analysis at a future point in time, e.g. when you get comments from examiners or reviewers.
We first need a little more data wrangling to get the data into the format RMark wants:
ch1 <- apply(ch, 1 , paste , collapse = "" ) # capture history needs to be single text column
hade <- data.frame(ch=ch1, group = 1)
hade$ch <- as.character(hade$ch) # the capture histories need to be a character and entitled 'ch'
Now, we are ready to run our first model. Let’s go with all the defaults and see what happens:
# run default model Phi()P()
model1 <- mark(hade)
##
## Output summary for CJS model
## Name : Phi(~1)p(~1)
##
## Npar : 2
## -2lnL: 964.9276
## AICc : 968.9553
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 0.5803299 0.1004607 0.3834269 0.7772329
## p:(Intercept) -0.3789708 0.1363474 -0.6462117 -0.1117300
##
##
## Real Parameter Phi
##
## 1 2 3 4 5 6 7
## 1 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433
## 2 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433
## 3 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433
## 4 0.6411433 0.6411433 0.6411433 0.6411433
## 5 0.6411433 0.6411433 0.6411433
## 6 0.6411433 0.6411433
## 7 0.6411433
## 8
## 9
## 10
## 11
## 8 9 10 11
## 1 0.6411433 0.6411433 0.6411433 0.6411433
## 2 0.6411433 0.6411433 0.6411433 0.6411433
## 3 0.6411433 0.6411433 0.6411433 0.6411433
## 4 0.6411433 0.6411433 0.6411433 0.6411433
## 5 0.6411433 0.6411433 0.6411433 0.6411433
## 6 0.6411433 0.6411433 0.6411433 0.6411433
## 7 0.6411433 0.6411433 0.6411433 0.6411433
## 8 0.6411433 0.6411433 0.6411433 0.6411433
## 9 0.6411433 0.6411433 0.6411433
## 10 0.6411433 0.6411433
## 11 0.6411433
##
##
## Real Parameter p
##
## 2 3 4 5 6 7 8
## 1 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751
## 2 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751
## 3 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751
## 4 0.4063751 0.4063751 0.4063751 0.4063751
## 5 0.4063751 0.4063751 0.4063751
## 6 0.4063751 0.4063751
## 7 0.4063751
## 8
## 9
## 10
## 11
## 9 10 11 12
## 1 0.4063751 0.4063751 0.4063751 0.4063751
## 2 0.4063751 0.4063751 0.4063751 0.4063751
## 3 0.4063751 0.4063751 0.4063751 0.4063751
## 4 0.4063751 0.4063751 0.4063751 0.4063751
## 5 0.4063751 0.4063751 0.4063751 0.4063751
## 6 0.4063751 0.4063751 0.4063751 0.4063751
## 7 0.4063751 0.4063751 0.4063751 0.4063751
## 8 0.4063751 0.4063751 0.4063751 0.4063751
## 9 0.4063751 0.4063751 0.4063751
## 10 0.4063751 0.4063751
## 11 0.4063751
We can see that RMark has made MARK fit a Cormack-Jolly-Seber (CJS) model with constand survival (Phi) and recapture (p) probabilities.
For now, this isn’t such a useful model and we don’t need it yet. MARK creates a number of files in your working directory and so it is a good idea to properly remove models you don’t need any more. With the rm command, we remove the model from our work space and with the cleanup command, we remove the associated files (now ‘orphaned’) from our working directory.
rm(model1) # remove a model
cleanup(ask=FALSE) # removes 'orphaned' files
To run the models we are interested in, we need to take a bit more control of the process, rather than going with all the defaults. It turns out that when we called the mark function above, it carried out a number of processing steps in the background. We now do these steps explicitly. The first step is to process the data and it is her that we can tell MARK what kind of analysis we want to do. In our case, we want to fit CJS-type of models. We can also tell MARK that our first occasion was the year 2006.
The next step is to create the design data, i.e. all the columns we might need for the design matrix of any of our models. The concept of design matrix is a bit special here and this is where it helps to first get to grips with how MARK works before trying to steer it via RMark. The concept is slightly different from what it is in linear models but the ideas are the same. In MARK, the design matrix maps an underlying flexible model structure (where every possible survival and recapture probability has its own parameter in the case of RMark) onto a more constrained model that follows linear model thinking. Essentially, we will constrain the underlying parameters to be a linear combination of a series of ‘beta’ parameters. This allows us to estimate survival and recapture probabilities as a function of covariates and to build additive models.
Let’s run these two steps and have a look at the design data for survival and recapture.
# take a bit more control
hade.process <- process.data(hade,model="CJS", begin.time=2006)
hade.ddl <- make.design.data(hade.process)
head(hade.ddl$Phi)
## par.index model.index group cohort age time occ.cohort Cohort Age Time
## 1 1 1 1 2006 0 2006 1 0 0 0
## 2 2 2 1 2006 1 2007 1 0 1 1
## 3 3 3 1 2006 2 2008 1 0 2 2
## 4 4 4 1 2006 3 2009 1 0 3 3
## 5 5 5 1 2006 4 2010 1 0 4 4
## 6 6 6 1 2006 5 2011 1 0 5 5
head(hade.ddl$p)
## par.index model.index group cohort age time occ.cohort Cohort Age Time
## 1 1 67 1 2006 1 2007 1 0 1 0
## 2 2 68 1 2006 2 2008 1 0 2 1
## 3 3 69 1 2006 3 2009 1 0 3 2
## 4 4 70 1 2006 4 2010 1 0 4 3
## 5 5 71 1 2006 5 2011 1 0 5 4
## 6 6 72 1 2006 6 2012 1 0 6 5
With the design data in place, we need to set up some model structures. Following standard notation, ‘dot’ refers to the constant model.
dot <- list(formula=~1)
time <- list(formula=~time)
age <- list(formula=~age)
timeage <- list(formula=~time + age)
Now, we are ready to run some models. Let’s first run the fully time-dependent model, i.e. survival and recapture probabilities are allowed to vary from year to year.
Phi.time.P.time <- mark(hade.process,hade.ddl,model.parameters=list(Phi=time, p=time))
##
## Note: only 21 parameters counted of 22 specified parameters
## AICc and parameter count have been adjusted upward
##
## Output summary for CJS model
## Name : Phi(~time)p(~time)
##
## Npar : 22 (unadjusted=21)
## -2lnL: 930.9508
## AICc : 977.4012 (unadjusted=975.18273)
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 1.0822087 1.3641367 -1.5914993 3.7559168
## Phi:time2007 -1.6754363 1.5170495 -4.6488534 1.2979809
## Phi:time2008 0.3786464 1.6348750 -2.8257086 3.5830014
## Phi:time2009 -0.3607545 1.4596715 -3.2217107 2.5002016
## Phi:time2010 -1.2869290 1.4071798 -4.0450014 1.4711434
## Phi:time2011 1.0321527 2.4062831 -3.6841623 5.7484676
## Phi:time2012 -0.8350509 1.4660599 -3.7085285 2.0384266
## Phi:time2013 1.2778288 2.8638457 -4.3353087 6.8909664
## Phi:time2014 -0.9928638 1.4787792 -3.8912710 1.9055435
## Phi:time2015 -0.5416363 1.7369783 -3.9461139 2.8628413
## Phi:time2016 -0.1193493 310.6679600 -609.0285500 608.7898500
## p:(Intercept) -0.4192546 0.7072455 -1.8054558 0.9669466
## p:time2008 -0.0033050 0.8439037 -1.6573562 1.6507463
## p:time2009 -0.3403823 0.7875442 -1.8839689 1.2032043
## p:time2010 0.0994832 0.7788217 -1.4270075 1.6259738
## p:time2011 0.3706798 0.7908849 -1.1794545 1.9208141
## p:time2012 -0.0777099 0.8165128 -1.6780750 1.5226553
## p:time2013 -0.1785822 0.8237765 -1.7931841 1.4360198
## p:time2014 0.6192760 0.8699408 -1.0858081 2.3243600
## p:time2015 1.5384861 1.0733044 -0.5651906 3.6421627
## p:time2016 0.2650998 1.0911840 -1.8736209 2.4038206
## p:time2017 0.7923174 210.4966800 -411.7811800 413.3658200
##
##
## Real Parameter Phi
##
## 2006 2007 2008 2009 2010 2011 2012
## 2006 0.7469117 0.3558947 0.8116634 0.6729272 0.4489979 0.8922912 0.5614768
## 2007 0.3558947 0.8116634 0.6729272 0.4489979 0.8922912 0.5614768
## 2008 0.8116634 0.6729272 0.4489979 0.8922912 0.5614768
## 2009 0.6729272 0.4489979 0.8922912 0.5614768
## 2010 0.4489979 0.8922912 0.5614768
## 2011 0.8922912 0.5614768
## 2012 0.5614768
## 2013
## 2014
## 2015
## 2016
## 2013 2014 2015 2016
## 2006 0.9137288 0.5223214 0.6319456 0.7236939
## 2007 0.9137288 0.5223214 0.6319456 0.7236939
## 2008 0.9137288 0.5223214 0.6319456 0.7236939
## 2009 0.9137288 0.5223214 0.6319456 0.7236939
## 2010 0.9137288 0.5223214 0.6319456 0.7236939
## 2011 0.9137288 0.5223214 0.6319456 0.7236939
## 2012 0.9137288 0.5223214 0.6319456 0.7236939
## 2013 0.9137288 0.5223214 0.6319456 0.7236939
## 2014 0.5223214 0.6319456 0.7236939
## 2015 0.6319456 0.7236939
## 2016 0.7236939
##
##
## Real Parameter p
##
## 2007 2008 2009 2010 2011 2012 2013
## 2006 0.3966951 0.3959044 0.3187251 0.4207315 0.4878587 0.3782543 0.3548388
## 2007 0.3959044 0.3187251 0.4207315 0.4878587 0.3782543 0.3548388
## 2008 0.3187251 0.4207315 0.4878587 0.3782543 0.3548388
## 2009 0.4207315 0.4878587 0.3782543 0.3548388
## 2010 0.4878587 0.3782543 0.3548388
## 2011 0.3782543 0.3548388
## 2012 0.3548388
## 2013
## 2014
## 2015
## 2016
## 2014 2015 2016 2017
## 2006 0.5498393 0.7538461 0.4615374 0.5921988
## 2007 0.5498393 0.7538461 0.4615374 0.5921988
## 2008 0.5498393 0.7538461 0.4615374 0.5921988
## 2009 0.5498393 0.7538461 0.4615374 0.5921988
## 2010 0.5498393 0.7538461 0.4615374 0.5921988
## 2011 0.5498393 0.7538461 0.4615374 0.5921988
## 2012 0.5498393 0.7538461 0.4615374 0.5921988
## 2013 0.5498393 0.7538461 0.4615374 0.5921988
## 2014 0.7538461 0.4615374 0.5921988
## 2015 0.4615374 0.5921988
## 2016 0.5921988
Now have a look at the parameter estimates:
Phi.time.P.time$results$real # parameter estimates
## estimate se lcl ucl fixed
## Phi g1 c2006 a0 t2006 0.7469117 0.2578690 1.691731e-01 0.9771551
## Phi g1 c2006 a1 t2007 0.3558947 0.0828906 2.138351e-01 0.5288445
## Phi g1 c2006 a2 t2008 0.8116634 0.1377444 4.242717e-01 0.9618368
## Phi g1 c2006 a3 t2009 0.6729272 0.1143159 4.264007e-01 0.8506185
## Phi g1 c2006 a4 t2010 0.4489979 0.0854461 2.928337e-01 0.6159092
## Phi g1 c2006 a5 t2011 0.8922912 0.1905111 1.454322e-01 0.9975264
## Phi g1 c2006 a6 t2012 0.5614768 0.1322385 3.088490e-01 0.7858049
## Phi g1 c2006 a7 t2013 0.9137288 0.1984970 7.073770e-02 0.9993219
## Phi g1 c2006 a8 t2014 0.5223214 0.1424394 2.631597e-01 0.7699980
## Phi g1 c2006 a9 t2015 0.6319456 0.2500971 1.726511e-01 0.9338937
## Phi g1 c2006 a10 t2016 0.7236939 62.1223730 9.299577e-265 1.0000000
## p g1 c2006 a1 t2007 0.3966951 0.1692637 1.411882e-01 0.7245105
## p g1 c2006 a2 t2008 0.3959044 0.1101137 2.099935e-01 0.6177114
## p g1 c2006 a3 t2009 0.3187251 0.0752285 1.917492e-01 0.4798634
## p g1 c2006 a4 t2010 0.4207315 0.0794853 2.770779e-01 0.5791927
## p g1 c2006 a5 t2011 0.4878587 0.0884429 3.224847e-01 0.6559340
## p g1 c2006 a6 t2012 0.3782543 0.0959629 2.147158e-01 0.5751272
## p g1 c2006 a7 t2013 0.3548388 0.0966958 1.937679e-01 0.5572570
## p g1 c2006 a8 t2014 0.5498393 0.1253811 3.115611e-01 0.7672546
## p g1 c2006 a9 t2015 0.7538461 0.1498104 3.862410e-01 0.9371212
## p g1 c2006 a10 t2016 0.4615374 0.2065083 1.439556e-01 0.8137416
## p g1 c2006 a11 t2017 0.5921988 50.8352120 9.599050e-180 1.0000000
## note
## Phi g1 c2006 a0 t2006
## Phi g1 c2006 a1 t2007
## Phi g1 c2006 a2 t2008
## Phi g1 c2006 a3 t2009
## Phi g1 c2006 a4 t2010
## Phi g1 c2006 a5 t2011
## Phi g1 c2006 a6 t2012
## Phi g1 c2006 a7 t2013
## Phi g1 c2006 a8 t2014
## Phi g1 c2006 a9 t2015
## Phi g1 c2006 a10 t2016
## p g1 c2006 a1 t2007
## p g1 c2006 a2 t2008
## p g1 c2006 a3 t2009
## p g1 c2006 a4 t2010
## p g1 c2006 a5 t2011
## p g1 c2006 a6 t2012
## p g1 c2006 a7 t2013
## p g1 c2006 a8 t2014
## p g1 c2006 a9 t2015
## p g1 c2006 a10 t2016
## p g1 c2006 a11 t2017
Phi.time.P.time$results$beta # parameter estimates on logit scale
## estimate se lcl ucl
## Phi:(Intercept) 1.0822087 1.3641367 -1.5914993 3.7559168
## Phi:time2007 -1.6754363 1.5170495 -4.6488534 1.2979809
## Phi:time2008 0.3786464 1.6348750 -2.8257086 3.5830014
## Phi:time2009 -0.3607545 1.4596715 -3.2217107 2.5002016
## Phi:time2010 -1.2869290 1.4071798 -4.0450014 1.4711434
## Phi:time2011 1.0321527 2.4062831 -3.6841623 5.7484676
## Phi:time2012 -0.8350509 1.4660599 -3.7085285 2.0384266
## Phi:time2013 1.2778288 2.8638457 -4.3353087 6.8909664
## Phi:time2014 -0.9928638 1.4787792 -3.8912710 1.9055435
## Phi:time2015 -0.5416363 1.7369783 -3.9461139 2.8628413
## Phi:time2016 -0.1193493 310.6679600 -609.0285500 608.7898500
## p:(Intercept) -0.4192546 0.7072455 -1.8054558 0.9669466
## p:time2008 -0.0033050 0.8439037 -1.6573562 1.6507463
## p:time2009 -0.3403823 0.7875442 -1.8839689 1.2032043
## p:time2010 0.0994832 0.7788217 -1.4270075 1.6259738
## p:time2011 0.3706798 0.7908849 -1.1794545 1.9208141
## p:time2012 -0.0777099 0.8165128 -1.6780750 1.5226553
## p:time2013 -0.1785822 0.8237765 -1.7931841 1.4360198
## p:time2014 0.6192760 0.8699408 -1.0858081 2.3243600
## p:time2015 1.5384861 1.0733044 -0.5651906 3.6421627
## p:time2016 0.2650998 1.0911840 -1.8736209 2.4038206
## p:time2017 0.7923174 210.4966800 -411.7811800 413.3658200
The summary function is also useful here:
summary(Phi.time.P.time)
## Output summary for CJS model
## Name : Phi(~time)p(~time)
##
## Npar : 22 (unadjusted=21)
## -2lnL: 930.9508
## AICc : 977.4012 (unadjusted=975.18273)
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 1.0822087 1.3641367 -1.5914993 3.7559168
## Phi:time2007 -1.6754363 1.5170495 -4.6488534 1.2979809
## Phi:time2008 0.3786464 1.6348750 -2.8257086 3.5830014
## Phi:time2009 -0.3607545 1.4596715 -3.2217107 2.5002016
## Phi:time2010 -1.2869290 1.4071798 -4.0450014 1.4711434
## Phi:time2011 1.0321527 2.4062831 -3.6841623 5.7484676
## Phi:time2012 -0.8350509 1.4660599 -3.7085285 2.0384266
## Phi:time2013 1.2778288 2.8638457 -4.3353087 6.8909664
## Phi:time2014 -0.9928638 1.4787792 -3.8912710 1.9055435
## Phi:time2015 -0.5416363 1.7369783 -3.9461139 2.8628413
## Phi:time2016 -0.1193493 310.6679600 -609.0285500 608.7898500
## p:(Intercept) -0.4192546 0.7072455 -1.8054558 0.9669466
## p:time2008 -0.0033050 0.8439037 -1.6573562 1.6507463
## p:time2009 -0.3403823 0.7875442 -1.8839689 1.2032043
## p:time2010 0.0994832 0.7788217 -1.4270075 1.6259738
## p:time2011 0.3706798 0.7908849 -1.1794545 1.9208141
## p:time2012 -0.0777099 0.8165128 -1.6780750 1.5226553
## p:time2013 -0.1785822 0.8237765 -1.7931841 1.4360198
## p:time2014 0.6192760 0.8699408 -1.0858081 2.3243600
## p:time2015 1.5384861 1.0733044 -0.5651906 3.6421627
## p:time2016 0.2650998 1.0911840 -1.8736209 2.4038206
## p:time2017 0.7923174 210.4966800 -411.7811800 413.3658200
##
##
## Real Parameter Phi
##
## 2006 2007 2008 2009 2010 2011 2012
## 2006 0.7469117 0.3558947 0.8116634 0.6729272 0.4489979 0.8922912 0.5614768
## 2007 0.3558947 0.8116634 0.6729272 0.4489979 0.8922912 0.5614768
## 2008 0.8116634 0.6729272 0.4489979 0.8922912 0.5614768
## 2009 0.6729272 0.4489979 0.8922912 0.5614768
## 2010 0.4489979 0.8922912 0.5614768
## 2011 0.8922912 0.5614768
## 2012 0.5614768
## 2013
## 2014
## 2015
## 2016
## 2013 2014 2015 2016
## 2006 0.9137288 0.5223214 0.6319456 0.7236939
## 2007 0.9137288 0.5223214 0.6319456 0.7236939
## 2008 0.9137288 0.5223214 0.6319456 0.7236939
## 2009 0.9137288 0.5223214 0.6319456 0.7236939
## 2010 0.9137288 0.5223214 0.6319456 0.7236939
## 2011 0.9137288 0.5223214 0.6319456 0.7236939
## 2012 0.9137288 0.5223214 0.6319456 0.7236939
## 2013 0.9137288 0.5223214 0.6319456 0.7236939
## 2014 0.5223214 0.6319456 0.7236939
## 2015 0.6319456 0.7236939
## 2016 0.7236939
##
##
## Real Parameter p
##
## 2007 2008 2009 2010 2011 2012 2013
## 2006 0.3966951 0.3959044 0.3187251 0.4207315 0.4878587 0.3782543 0.3548388
## 2007 0.3959044 0.3187251 0.4207315 0.4878587 0.3782543 0.3548388
## 2008 0.3187251 0.4207315 0.4878587 0.3782543 0.3548388
## 2009 0.4207315 0.4878587 0.3782543 0.3548388
## 2010 0.4878587 0.3782543 0.3548388
## 2011 0.3782543 0.3548388
## 2012 0.3548388
## 2013
## 2014
## 2015
## 2016
## 2014 2015 2016 2017
## 2006 0.5498393 0.7538461 0.4615374 0.5921988
## 2007 0.5498393 0.7538461 0.4615374 0.5921988
## 2008 0.5498393 0.7538461 0.4615374 0.5921988
## 2009 0.5498393 0.7538461 0.4615374 0.5921988
## 2010 0.5498393 0.7538461 0.4615374 0.5921988
## 2011 0.5498393 0.7538461 0.4615374 0.5921988
## 2012 0.5498393 0.7538461 0.4615374 0.5921988
## 2013 0.5498393 0.7538461 0.4615374 0.5921988
## 2014 0.7538461 0.4615374 0.5921988
## 2015 0.4615374 0.5921988
## 2016 0.5921988
We see that RMark thinks that our model should have 22 parameters. However, if you look at the output carefully, you’ll see that the last Phi and the last p parameters have a standard error of 0. Standard errors that are very large or 0 indicate a problem. Often, this is because the data are too sparse and that parameter has not been well estimated. MARK tends not to count such parameters. That’s why RMark tells us that the unadjusted parameter count is 21. It’s the number of parameters MARK thinks are estimable. Where parameters are poorly (or not) estimated because of data sparseness, it makes sens to still count those parameters because they are part of the model structure and contribute to model complexity. Therefore, RMark by default counts all parameters we have specified. However, in some cases, parameters are technically not estimable because they are confounded. In that case, they should not contribute to the parameter count. It turns out that in this case, MARK is right: with the fully time-dependent CJS model, the final Phi and p estimate are confounded. With no further data, there is no way to tell whether all individuals survived but we had a low recapture probability on the last occasion or whether they all died and we saw everyone who was still alive (or anything in between). In fact, only the product of those two parameters is estimable and we therefore should only count them as a single parameter.
So we need to tell RMark that we want to change the parameter count:
Phi.time.P.time <- adjust.parameter.count(Phi.time.P.time, 21) # adjust parameter count because last Phi and P are confounded
##
## Number of parameters adjusted from 22 to 21
## Adjusted AICc = 975.182724057971
## Unadjusted AICc = 975.18273
Now, we are ready to run a few more models:
Phi.dot.P.dot <- mark(hade.process,hade.ddl,model.parameters=list(Phi=dot, p=dot))
##
## Output summary for CJS model
## Name : Phi(~1)p(~1)
##
## Npar : 2
## -2lnL: 964.9276
## AICc : 968.9553
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 0.5803299 0.1004607 0.3834269 0.7772329
## p:(Intercept) -0.3789708 0.1363474 -0.6462117 -0.1117300
##
##
## Real Parameter Phi
##
## 2006 2007 2008 2009 2010 2011 2012
## 2006 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433
## 2007 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433
## 2008 0.6411433 0.6411433 0.6411433 0.6411433 0.6411433
## 2009 0.6411433 0.6411433 0.6411433 0.6411433
## 2010 0.6411433 0.6411433 0.6411433
## 2011 0.6411433 0.6411433
## 2012 0.6411433
## 2013
## 2014
## 2015
## 2016
## 2013 2014 2015 2016
## 2006 0.6411433 0.6411433 0.6411433 0.6411433
## 2007 0.6411433 0.6411433 0.6411433 0.6411433
## 2008 0.6411433 0.6411433 0.6411433 0.6411433
## 2009 0.6411433 0.6411433 0.6411433 0.6411433
## 2010 0.6411433 0.6411433 0.6411433 0.6411433
## 2011 0.6411433 0.6411433 0.6411433 0.6411433
## 2012 0.6411433 0.6411433 0.6411433 0.6411433
## 2013 0.6411433 0.6411433 0.6411433 0.6411433
## 2014 0.6411433 0.6411433 0.6411433
## 2015 0.6411433 0.6411433
## 2016 0.6411433
##
##
## Real Parameter p
##
## 2007 2008 2009 2010 2011 2012 2013
## 2006 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751
## 2007 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751
## 2008 0.4063751 0.4063751 0.4063751 0.4063751 0.4063751
## 2009 0.4063751 0.4063751 0.4063751 0.4063751
## 2010 0.4063751 0.4063751 0.4063751
## 2011 0.4063751 0.4063751
## 2012 0.4063751
## 2013
## 2014
## 2015
## 2016
## 2014 2015 2016 2017
## 2006 0.4063751 0.4063751 0.4063751 0.4063751
## 2007 0.4063751 0.4063751 0.4063751 0.4063751
## 2008 0.4063751 0.4063751 0.4063751 0.4063751
## 2009 0.4063751 0.4063751 0.4063751 0.4063751
## 2010 0.4063751 0.4063751 0.4063751 0.4063751
## 2011 0.4063751 0.4063751 0.4063751 0.4063751
## 2012 0.4063751 0.4063751 0.4063751 0.4063751
## 2013 0.4063751 0.4063751 0.4063751 0.4063751
## 2014 0.4063751 0.4063751 0.4063751
## 2015 0.4063751 0.4063751
## 2016 0.4063751
Phi.age.P.dot <- mark(hade.process,hade.ddl,model.parameters=list(Phi=age, p=dot))
##
## Note: only 9 parameters counted of 12 specified parameters
## AICc and parameter count have been adjusted upward
##
## Output summary for CJS model
## Name : Phi(~age)p(~1)
##
## Npar : 12 (unadjusted=9)
## -2lnL: 936.8426
## AICc : 961.5801 (unadjusted=955.26509)
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 0.0206908 2.123647e-01 -0.3955441 0.4369256
## Phi:age1 0.0625891 4.229103e-01 -0.7663151 0.8914932
## Phi:age2 1.1876704 6.385317e-01 -0.0638518 2.4391927
## Phi:age3 1.8487813 1.172957e+00 -0.4502145 4.1477772
## Phi:age4 0.6574010 6.286634e-01 -0.5747793 1.8895813
## Phi:age5 1.0416520 8.691471e-01 -0.6618763 2.7451804
## Phi:age6 2.2800434 2.223503e+00 -2.0780222 6.6381090
## Phi:age7 49.5993070 8.146174e-06 49.5992910 49.5993230
## Phi:age8 -0.6350147 8.334957e-01 -2.2686663 0.9986369
## Phi:age9 69.9390390 1.182904e-14 69.9390390 69.9390390
## Phi:age10 -12.2810440 2.845141e+03 -5588.7582000 5564.1961000
## p:(Intercept) -0.1565221 1.424354e-01 -0.4356955 0.1226514
##
##
## Real Parameter Phi
##
## 2006 2007 2008 2009 2010 2011 2012
## 2006 0.5051725 0.5208079 0.7700089 0.8663972 0.6633127 0.7431380 0.9089378
## 2007 0.5051725 0.5208079 0.7700089 0.8663972 0.6633127 0.7431380
## 2008 0.5051725 0.5208079 0.7700089 0.8663972 0.6633127
## 2009 0.5051725 0.5208079 0.7700089 0.8663972
## 2010 0.5051725 0.5208079 0.7700089
## 2011 0.5051725 0.5208079
## 2012 0.5051725
## 2013
## 2014
## 2015
## 2016
## 2013 2014 2015 2016
## 2006 1.0000000 0.3510735 1.0000000 4.735809e-06
## 2007 0.9089378 1.0000000 0.3510735 1.000000e+00
## 2008 0.7431380 0.9089378 1.0000000 3.510735e-01
## 2009 0.6633127 0.7431380 0.9089378 1.000000e+00
## 2010 0.8663972 0.6633127 0.7431380 9.089378e-01
## 2011 0.7700089 0.8663972 0.6633127 7.431380e-01
## 2012 0.5208079 0.7700089 0.8663972 6.633127e-01
## 2013 0.5051725 0.5208079 0.7700089 8.663972e-01
## 2014 0.5051725 0.5208079 7.700089e-01
## 2015 0.5051725 5.208079e-01
## 2016 5.051725e-01
##
##
## Real Parameter p
##
## 2007 2008 2009 2010 2011 2012 2013
## 2006 0.4609492 0.4609492 0.4609492 0.4609492 0.4609492 0.4609492 0.4609492
## 2007 0.4609492 0.4609492 0.4609492 0.4609492 0.4609492 0.4609492
## 2008 0.4609492 0.4609492 0.4609492 0.4609492 0.4609492
## 2009 0.4609492 0.4609492 0.4609492 0.4609492
## 2010 0.4609492 0.4609492 0.4609492
## 2011 0.4609492 0.4609492
## 2012 0.4609492
## 2013
## 2014
## 2015
## 2016
## 2014 2015 2016 2017
## 2006 0.4609492 0.4609492 0.4609492 0.4609492
## 2007 0.4609492 0.4609492 0.4609492 0.4609492
## 2008 0.4609492 0.4609492 0.4609492 0.4609492
## 2009 0.4609492 0.4609492 0.4609492 0.4609492
## 2010 0.4609492 0.4609492 0.4609492 0.4609492
## 2011 0.4609492 0.4609492 0.4609492 0.4609492
## 2012 0.4609492 0.4609492 0.4609492 0.4609492
## 2013 0.4609492 0.4609492 0.4609492 0.4609492
## 2014 0.4609492 0.4609492 0.4609492
## 2015 0.4609492 0.4609492
## 2016 0.4609492
Phi.age.P.time <- mark(hade.process,hade.ddl,model.parameters=list(Phi=age, p=time))
##
## Note: only 19 parameters counted of 22 specified parameters
##
## AICc and parameter count have been adjusted upward
##
## Output summary for CJS model
## Name : Phi(~age)p(~time)
##
## Npar : 22 (unadjusted=19)
## -2lnL: 923.1815
## AICc : 969.6318 (unadjusted=963.00839)
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 0.0775336 2.333981e-01 -0.3799268 0.5349939
## Phi:age1 -0.0083096 4.599517e-01 -0.9098150 0.8931957
## Phi:age2 1.0005729 6.294298e-01 -0.2331095 2.2342554
## Phi:age3 2.0621411 1.555769e+00 -0.9871658 5.1114479
## Phi:age4 0.4828085 6.170069e-01 -0.7265251 1.6921420
## Phi:age5 0.8679830 7.596263e-01 -0.6208846 2.3568506
## Phi:age6 2.0050019 1.830331e+00 -1.5824474 5.5924513
## Phi:age7 85.1746540 5.345200e+00 74.6980620 95.6512470
## Phi:age8 -0.5699461 9.663291e-01 -2.4639512 1.3240590
## Phi:age9 74.5221230 7.550666e-06 74.5221080 74.5221380
## Phi:age10 -10.3449230 1.543812e+03 -3036.2174000 3015.5276000
## p:(Intercept) 0.1664723 6.181851e-01 -1.0451705 1.3781151
## p:time2008 -0.9468654 6.983980e-01 -2.3157254 0.4219947
## p:time2009 -0.5616410 6.761540e-01 -1.8869029 0.7636210
## p:time2010 -0.0660631 6.529538e-01 -1.3458527 1.2137264
## p:time2011 -0.3746375 6.778284e-01 -1.7031811 0.9539062
## p:time2012 -0.3512292 7.058833e-01 -1.7347605 1.0323020
## p:time2013 -0.8913780 7.210137e-01 -2.3045650 0.5218089
## p:time2014 0.4255695 7.533746e-01 -1.0510448 1.9021837
## p:time2015 0.3966368 8.567353e-01 -1.2825644 2.0758380
## p:time2016 -0.7649789 8.712432e-01 -2.4726156 0.9426579
## p:time2017 -0.2348356 1.012909e+00 -2.2201372 1.7504660
##
##
## Real Parameter Phi
##
## 2006 2007 2008 2009 2010 2011 2012
## 2006 0.5193737 0.5172991 0.7461355 0.8947000 0.6365317 0.7202126 0.8891941
## 2007 0.5193737 0.5172991 0.7461355 0.8947000 0.6365317 0.7202126
## 2008 0.5193737 0.5172991 0.7461355 0.8947000 0.6365317
## 2009 0.5193737 0.5172991 0.7461355 0.8947000
## 2010 0.5193737 0.5172991 0.7461355
## 2011 0.5193737 0.5172991
## 2012 0.5193737
## 2013
## 2014
## 2015
## 2016
## 2013 2014 2015 2016
## 2006 1.0000000 0.3793254 1.0000000 3.474678e-05
## 2007 0.8891941 1.0000000 0.3793254 1.000000e+00
## 2008 0.7202126 0.8891941 1.0000000 3.793254e-01
## 2009 0.6365317 0.7202126 0.8891941 1.000000e+00
## 2010 0.8947000 0.6365317 0.7202126 8.891941e-01
## 2011 0.7461355 0.8947000 0.6365317 7.202126e-01
## 2012 0.5172991 0.7461355 0.8947000 6.365317e-01
## 2013 0.5193737 0.5172991 0.7461355 8.947000e-01
## 2014 0.5193737 0.5172991 7.461355e-01
## 2015 0.5193737 5.172991e-01
## 2016 5.193737e-01
##
##
## Real Parameter p
##
## 2007 2008 2009 2010 2011 2012 2013
## 2006 0.5415222 0.3142352 0.4024737 0.5250812 0.4481458 0.4539417 0.3263136
## 2007 0.3142352 0.4024737 0.5250812 0.4481458 0.4539417 0.3263136
## 2008 0.4024737 0.5250812 0.4481458 0.4539417 0.3263136
## 2009 0.5250812 0.4481458 0.4539417 0.3263136
## 2010 0.4481458 0.4539417 0.3263136
## 2011 0.4539417 0.3263136
## 2012 0.3263136
## 2013
## 2014
## 2015
## 2016
## 2014 2015 2016 2017
## 2006 0.6438335 0.6371716 0.3546854 0.4829158
## 2007 0.6438335 0.6371716 0.3546854 0.4829158
## 2008 0.6438335 0.6371716 0.3546854 0.4829158
## 2009 0.6438335 0.6371716 0.3546854 0.4829158
## 2010 0.6438335 0.6371716 0.3546854 0.4829158
## 2011 0.6438335 0.6371716 0.3546854 0.4829158
## 2012 0.6438335 0.6371716 0.3546854 0.4829158
## 2013 0.6438335 0.6371716 0.3546854 0.4829158
## 2014 0.6371716 0.3546854 0.4829158
## 2015 0.3546854 0.4829158
## 2016 0.4829158
Phi.timeage.P.dot <- mark(hade.process,hade.ddl,model.parameters=list(Phi=timeage, p=dot))
##
## Note: only 15 parameters counted of 22 specified parameters
##
## AICc and parameter count have been adjusted upward
##
## Output summary for CJS model
## Name : Phi(~time + age)p(~1)
##
## Npar : 22 (unadjusted=15)
## -2lnL: 903.5938
## AICc : 950.0441 (unadjusted=934.73664)
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 0.7564328 0.7807013 -0.7737418 2.2866075
## Phi:time2007 -1.4433637 0.8675556 -3.1437728 0.2570453
## Phi:time2008 -0.1747138 0.8230125 -1.7878184 1.4383908
## Phi:time2009 -0.1604063 0.8288442 -1.7849410 1.4641283
## Phi:time2010 -1.5109004 0.8081384 -3.0948518 0.0730510
## Phi:time2011 0.1469165 1.0519597 -1.9149245 2.2087575
## Phi:time2012 -1.1201940 0.9027083 -2.8895024 0.6491144
## Phi:time2013 14.7977310 0.0000000 14.7977310 14.7977310
## Phi:time2014 -25.3638640 333.9256600 -679.8581700 629.1304400
## Phi:time2015 -2.8014329 1.4846581 -5.7113628 0.1084971
## Phi:time2016 12.8482220 1093.5212000 -2130.4534000 2156.1499000
## Phi:age1 -0.1207040 0.4254795 -0.9546439 0.7132358
## Phi:age2 45.5827580 0.0000000 45.5827580 45.5827580
## Phi:age3 1.0131786 0.5695258 -0.1030919 2.1294492
## Phi:age4 0.5471955 0.6524114 -0.7315308 1.8259217
## Phi:age5 24.3098730 333.9263500 -630.1857800 678.8055300
## Phi:age6 24.9994900 333.9263500 -629.4961700 679.4951500
## Phi:age7 54.0310250 0.0000000 54.0310250 54.0310250
## Phi:age8 0.5917444 1.7290288 -2.7971521 3.9806410
## Phi:age9 12.6406500 0.0000000 12.6406500 12.6406500
## Phi:age10 -0.1700383 0.0000000 -0.1700383 -0.1700383
## p:(Intercept) -0.1936090 0.1352251 -0.4586502 0.0714322
##
##
## Real Parameter Phi
##
## 2006 2007 2008 2009 2010 2011 2012
## 2006 0.6805788 0.3083947 1.0000000 0.8333010 0.4483667 1.0000000 1.0000000
## 2007 0.3347162 0.6132549 1.0000000 0.5643194 0.8100823 1.0000000
## 2008 0.6414629 0.6166427 1.0000000 0.8717508 0.5457304
## 2009 0.6447467 0.2941793 1.0000000 0.6568792
## 2010 0.3198486 0.6862500 1.0000000
## 2011 0.7116373 0.3811983
## 2012 0.4100494
## 2013
## 2014
## 2015
## 2016
## 2013 2014 2015 2016
## 2006 1.0000000 3.716378e-11 0.9999750 0.9999985
## 2007 1.0000000 1.000000e+00 0.1895010 1.0000000
## 2008 1.0000000 5.967783e-01 1.0000000 0.9999993
## 2009 0.9999999 4.261545e-01 1.0000000 1.0000000
## 2010 0.9999999 3.554450e-11 1.0000000 1.0000000
## 2011 1.0000000 5.664302e-11 0.1827532 1.0000000
## 2012 0.9999998 1.000000e+00 0.2627311 0.9999993
## 2013 0.9999998 1.822667e-11 1.0000000 0.9999996
## 2014 2.056499e-11 0.1028728 1.0000000
## 2015 0.1145586 0.9999986
## 2016 0.9999988
##
##
## Real Parameter p
##
## 2007 2008 2009 2010 2011 2012 2013
## 2006 0.4517484 0.4517484 0.4517484 0.4517484 0.4517484 0.4517484 0.4517484
## 2007 0.4517484 0.4517484 0.4517484 0.4517484 0.4517484 0.4517484
## 2008 0.4517484 0.4517484 0.4517484 0.4517484 0.4517484
## 2009 0.4517484 0.4517484 0.4517484 0.4517484
## 2010 0.4517484 0.4517484 0.4517484
## 2011 0.4517484 0.4517484
## 2012 0.4517484
## 2013
## 2014
## 2015
## 2016
## 2014 2015 2016 2017
## 2006 0.4517484 0.4517484 0.4517484 0.4517484
## 2007 0.4517484 0.4517484 0.4517484 0.4517484
## 2008 0.4517484 0.4517484 0.4517484 0.4517484
## 2009 0.4517484 0.4517484 0.4517484 0.4517484
## 2010 0.4517484 0.4517484 0.4517484 0.4517484
## 2011 0.4517484 0.4517484 0.4517484 0.4517484
## 2012 0.4517484 0.4517484 0.4517484 0.4517484
## 2013 0.4517484 0.4517484 0.4517484 0.4517484
## 2014 0.4517484 0.4517484 0.4517484
## 2015 0.4517484 0.4517484
## 2016 0.4517484
Phi.timeage.P.time <- mark(hade.process,hade.ddl,model.parameters=list(Phi=timeage, p=time))
##
## Note: only 28 parameters counted of 32 specified parameters
##
## AICc and parameter count have been adjusted upward
##
## Output summary for CJS model
## Name : Phi(~time + age)p(~time)
##
## Npar : 32 (unadjusted=28)
## -2lnL: 894.2482
## AICc : 963.4889 (unadjusted=954.2384)
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 0.9063589 1.3015972 -1.6447716 3.4574894
## Phi:time2007 -1.6204430 1.4100245 -4.3840911 1.1432050
## Phi:time2008 -0.2834119 1.2306417 -2.6954697 2.1286459
## Phi:time2009 -0.4300336 1.3149494 -3.0073344 2.1472672
## Phi:time2010 -1.7081780 1.2458527 -4.1500493 0.7336933
## Phi:time2011 0.2941235 1.5617350 -2.7668772 3.3551242
## Phi:time2012 -0.8151621 1.5256484 -3.8054329 2.1751087
## Phi:time2013 0.1387575 1.9719770 -3.7263174 4.0038325
## Phi:time2014 -24.4092750 386.0679800 -781.1025400 732.2839900
## Phi:time2015 -2.9286185 1.8268539 -6.5092522 0.6520151
## Phi:time2016 13.6569320 1486.0650000 -2899.0305000 2926.3444000
## Phi:age1 -0.0541976 0.4895458 -1.0137073 0.9053122
## Phi:age2 25.1034840 386.0771800 -731.6078100 781.8147800
## Phi:age3 1.0997105 0.7207172 -0.3128951 2.5123162
## Phi:age4 0.2962080 0.7718978 -1.2167117 1.8091277
## Phi:age5 22.9610420 386.0683400 -733.7329100 779.6549900
## Phi:age6 23.5166430 386.0680800 -733.1768000 780.2100900
## Phi:age7 25.1239480 386.0746100 -731.5822900 781.8301900
## Phi:age8 0.6802225 1.7980687 -2.8439922 4.2044372
## Phi:age9 9.1406636 0.0000000 9.1406636 9.1406636
## Phi:age10 -0.1184524 1300.7567000 -2549.6016000 2549.3647000
## p:(Intercept) -0.3392224 0.7625291 -1.8337795 1.1553347
## p:time2008 0.0859601 0.8920092 -1.6623781 1.8342982
## p:time2009 -0.0643481 0.7903281 -1.6133911 1.4846950
## p:time2010 0.1998344 0.7939023 -1.3562141 1.7558829
## p:time2011 0.4396746 0.8222024 -1.1718422 2.0511913
## p:time2012 -0.0331316 0.8321631 -1.6641713 1.5979082
## p:time2013 -0.5388487 0.8787082 -2.2611167 1.1834194
## p:time2014 0.2614249 0.8593844 -1.4229685 1.9458182
## p:time2015 1.4919641 1.1041269 -0.6721247 3.6560529
## p:time2016 0.5121381 1.0079552 -1.4634541 2.4877303
## p:time2017 0.3313746 0.9808244 -1.5910412 2.2537905
##
##
## Real Parameter Phi
##
## 2006 2007 2008 2009 2010 2011 2012
## 2006 0.7122545 0.3168509 1.0000000 0.8286424 0.3762229 1.0000000 1.0000000
## 2007 0.3286970 0.6384746 1.0000000 0.5739270 0.8170803 1.0000000
## 2008 0.6508885 0.6039923 1.0000000 0.9088930 0.5956578
## 2009 0.6168798 0.2981722 1.0000000 0.7669033
## 2010 0.3096365 0.7588317 1.0000000
## 2011 0.7686106 0.5092487
## 2012 0.5227834
## 2013
## 2014
## 2015
## 2016
## 2013 2014 2015 2016
## 2006 1.0000000 1.225266e-10 0.9991906 0.9999995
## 2007 1.0000000 8.349374e-01 0.2071753 1.0000000
## 2008 1.0000000 5.034317e-01 1.0000000 0.9999998
## 2009 0.7927077 3.677517e-01 1.0000000 1.0000000
## 2010 0.8951844 8.345548e-11 1.0000000 1.0000000
## 2011 1.0000000 1.863853e-10 0.1510933 1.0000000
## 2012 0.7292694 8.320977e-01 0.2844388 0.9999996
## 2013 0.7398360 5.878624e-11 1.0000000 0.9999998
## 2014 6.206023e-11 0.1114062 1.0000000
## 2015 0.1168855 0.9999995
## 2016 0.9999995
##
##
## Real Parameter p
##
## 2007 2008 2009 2010 2011 2012 2013
## 2006 0.4159984 0.4370207 0.4004548 0.4652093 0.525092 0.4079724 0.2935777
## 2007 0.4370207 0.4004548 0.4652093 0.525092 0.4079724 0.2935777
## 2008 0.4004548 0.4652093 0.525092 0.4079724 0.2935777
## 2009 0.4652093 0.525092 0.4079724 0.2935777
## 2010 0.525092 0.4079724 0.2935777
## 2011 0.4079724 0.2935777
## 2012 0.2935777
## 2013
## 2014
## 2015
## 2016
## 2014 2015 2016 2017
## 2006 0.4805604 0.7600113 0.5431215 0.4980381
## 2007 0.4805604 0.7600113 0.5431215 0.4980381
## 2008 0.4805604 0.7600113 0.5431215 0.4980381
## 2009 0.4805604 0.7600113 0.5431215 0.4980381
## 2010 0.4805604 0.7600113 0.5431215 0.4980381
## 2011 0.4805604 0.7600113 0.5431215 0.4980381
## 2012 0.4805604 0.7600113 0.5431215 0.4980381
## 2013 0.4805604 0.7600113 0.5431215 0.4980381
## 2014 0.7600113 0.5431215 0.4980381
## 2015 0.5431215 0.4980381
## 2016 0.4980381
Now that we’ve run a few models, we need to have a way to tell which one of these describes the data the best. That is a model selection question and we can use Akaike’s Information Criterion to rank the models:
# model selection
hade.cjs.results <- collect.models()
hade.cjs.results
## model npar AICc DeltaAICc weight Deviance
## 5 Phi(~time + age)p(~1) 22 950.0441 0.00000 9.955527e-01 250.3671
## 1 Phi(~age)p(~1) 12 961.5801 11.53601 3.112092e-03 283.6159
## 6 Phi(~time + age)p(~time) 32 963.4889 13.44478 1.198304e-03 241.0215
## 3 Phi(~1)p(~1) 2 968.9553 18.91115 7.790406e-05 311.7009
## 2 Phi(~age)p(~time) 22 969.6318 19.58769 5.554582e-05 269.9548
## 4 Phi(~time)p(~time) 21 975.1827 25.13858 3.461711e-06 277.7241
The best mode, according to AIC, is the one with additive affects of time and age on survival. That means, this model estimates a separate survival rate for each year and each age (1st, 2nd, 3rd, etc year) but without an interaction, i.e. a year with relatively low survival for one age class is also a year with low survival for all others. The recapture probability is set constant in this model.
So let’s have a look at the parameter estimates, first on the ‘real’ scale – i.e. these are survival probabilities – and then on the ‘beta’ scale – i.e. the logit scale.
# look at the parameter estimates
Phi.timeage.P.dot$results$real # parameter estimates
## estimate se lcl ucl
## Phi g1 c2006 a0 t2006 6.805788e-01 1.697177e-01 3.156702e-01 9.077618e-01
## Phi g1 c2006 a1 t2007 3.083947e-01 9.050090e-02 1.625616e-01 5.060049e-01
## Phi g1 c2006 a2 t2008 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2006 a3 t2009 8.333010e-01 9.067600e-02 5.817067e-01 9.472813e-01
## Phi g1 c2006 a4 t2010 4.483667e-01 1.668896e-01 1.780244e-01 7.531058e-01
## Phi g1 c2006 a5 t2011 1.000000e+00 3.747005e-09 1.000000e+00 1.000000e+00
## Phi g1 c2006 a6 t2012 1.000000e+00 6.675569e-09 1.000000e+00 1.000000e+00
## Phi g1 c2006 a7 t2013 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2006 a8 t2014 3.716378e-11 1.241011e-08 -2.428665e-08 2.436098e-08
## Phi g1 c2006 a9 t2015 9.999750e-01 0.000000e+00 9.999750e-01 9.999750e-01
## Phi g1 c2006 a10 t2016 9.999985e-01 1.006900e-03 6.736624e-299 1.000000e+00
## Phi g1 c2007 a0 t2007 3.347162e-01 7.207090e-02 2.106045e-01 4.868590e-01
## Phi g1 c2007 a1 t2008 6.132549e-01 1.014143e-01 4.068330e-01 7.856840e-01
## Phi g1 c2007 a2 t2009 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2007 a3 t2010 5.643194e-01 1.299635e-01 3.148918e-01 7.849534e-01
## Phi g1 c2007 a4 t2011 8.100823e-01 1.269734e-01 4.583302e-01 9.555599e-01
## Phi g1 c2007 a5 t2012 1.000000e+00 1.330412e-08 1.000000e+00 1.000000e+00
## Phi g1 c2007 a6 t2013 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2007 a7 t2014 1.000000e+00 5.196732e-11 1.000000e+00 1.000000e+00
## Phi g1 c2007 a8 t2015 1.895010e-01 1.775388e-01 2.368680e-02 6.926111e-01
## Phi g1 c2007 a9 t2016 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2008 a0 t2008 6.414629e-01 1.010047e-01 4.306840e-01 8.088414e-01
## Phi g1 c2008 a1 t2009 6.166427e-01 1.012878e-01 4.098760e-01 7.883691e-01
## Phi g1 c2008 a2 t2010 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2008 a3 t2011 8.717508e-01 8.322750e-02 6.124113e-01 9.669331e-01
## Phi g1 c2008 a4 t2012 5.457304e-01 1.365697e-01 2.898106e-01 7.795718e-01
## Phi g1 c2008 a5 t2013 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2008 a6 t2014 5.967783e-01 2.003055e-01 2.245301e-01 8.832504e-01
## Phi g1 c2008 a7 t2015 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2008 a8 t2016 9.999993e-01 7.471547e-04 1.443047e-298 1.000000e+00
## Phi g1 c2009 a0 t2009 6.447467e-01 9.057620e-02 4.553574e-01 7.975579e-01
## Phi g1 c2009 a1 t2010 2.941793e-01 7.145170e-02 1.751382e-01 4.499921e-01
## Phi g1 c2009 a2 t2011 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2009 a3 t2012 6.568792e-01 1.457061e-01 3.503150e-01 8.717465e-01
## Phi g1 c2009 a4 t2013 9.999999e-01 0.000000e+00 9.999999e-01 9.999999e-01
## Phi g1 c2009 a5 t2014 4.261545e-01 2.189464e-01 1.138126e-01 8.111130e-01
## Phi g1 c2009 a6 t2015 1.000000e+00 3.586294e-08 9.999999e-01 1.000000e+00
## Phi g1 c2009 a7 t2016 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2010 a0 t2010 3.198486e-01 8.273370e-02 1.824455e-01 4.977327e-01
## Phi g1 c2010 a1 t2011 6.862500e-01 1.698180e-01 3.179459e-01 9.112115e-01
## Phi g1 c2010 a2 t2012 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2010 a3 t2013 9.999999e-01 0.000000e+00 9.999999e-01 9.999999e-01
## Phi g1 c2010 a4 t2014 3.554450e-11 1.186925e-08 -2.322818e-08 2.329927e-08
## Phi g1 c2010 a5 t2015 1.000000e+00 7.147309e-08 9.999999e-01 1.000000e+00
## Phi g1 c2010 a6 t2016 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2011 a0 t2011 7.116373e-01 1.583273e-01 3.523175e-01 9.180067e-01
## Phi g1 c2011 a1 t2012 3.811983e-01 1.387486e-01 1.628273e-01 6.611471e-01
## Phi g1 c2011 a2 t2013 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2011 a3 t2014 5.664302e-11 1.891459e-08 -3.701595e-08 3.712923e-08
## Phi g1 c2011 a4 t2015 1.827532e-01 2.127021e-01 1.353130e-02 7.847416e-01
## Phi g1 c2011 a5 t2016 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2012 a0 t2012 4.100494e-01 1.285591e-01 1.969644e-01 6.632595e-01
## Phi g1 c2012 a1 t2013 9.999998e-01 0.000000e+00 9.999998e-01 9.999998e-01
## Phi g1 c2012 a2 t2014 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2012 a3 t2015 2.627311e-01 2.355879e-01 3.181050e-02 7.944545e-01
## Phi g1 c2012 a4 t2016 9.999993e-01 7.811892e-04 1.380171e-298 1.000000e+00
## Phi g1 c2013 a0 t2013 9.999998e-01 0.000000e+00 9.999998e-01 9.999998e-01
## Phi g1 c2013 a1 t2014 1.822667e-11 6.086359e-09 -1.191104e-08 1.194749e-08
## Phi g1 c2013 a2 t2015 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2013 a3 t2016 9.999996e-01 4.902106e-04 2.199414e-298 1.000000e+00
## Phi g1 c2014 a0 t2014 2.056499e-11 6.867181e-09 -1.343911e-08 1.348024e-08
## Phi g1 c2014 a1 t2015 1.028728e-01 1.215394e-01 8.604000e-03 6.024003e-01
## Phi g1 c2014 a2 t2016 1.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## Phi g1 c2015 a0 t2015 1.145586e-01 1.308648e-01 1.021460e-02 6.186130e-01
## Phi g1 c2015 a1 t2016 9.999986e-01 1.523400e-03 7.077305e-299 1.000000e+00
## Phi g1 c2016 a0 t2016 9.999988e-01 1.350200e-03 7.985259e-299 1.000000e+00
## p g1 c2006 a1 t2007 4.517484e-01 3.349140e-02 3.873061e-01 5.178505e-01
## fixed note
## Phi g1 c2006 a0 t2006
## Phi g1 c2006 a1 t2007
## Phi g1 c2006 a2 t2008
## Phi g1 c2006 a3 t2009
## Phi g1 c2006 a4 t2010
## Phi g1 c2006 a5 t2011
## Phi g1 c2006 a6 t2012
## Phi g1 c2006 a7 t2013
## Phi g1 c2006 a8 t2014
## Phi g1 c2006 a9 t2015
## Phi g1 c2006 a10 t2016
## Phi g1 c2007 a0 t2007
## Phi g1 c2007 a1 t2008
## Phi g1 c2007 a2 t2009
## Phi g1 c2007 a3 t2010
## Phi g1 c2007 a4 t2011
## Phi g1 c2007 a5 t2012
## Phi g1 c2007 a6 t2013
## Phi g1 c2007 a7 t2014
## Phi g1 c2007 a8 t2015
## Phi g1 c2007 a9 t2016
## Phi g1 c2008 a0 t2008
## Phi g1 c2008 a1 t2009
## Phi g1 c2008 a2 t2010
## Phi g1 c2008 a3 t2011
## Phi g1 c2008 a4 t2012
## Phi g1 c2008 a5 t2013
## Phi g1 c2008 a6 t2014
## Phi g1 c2008 a7 t2015
## Phi g1 c2008 a8 t2016
## Phi g1 c2009 a0 t2009
## Phi g1 c2009 a1 t2010
## Phi g1 c2009 a2 t2011
## Phi g1 c2009 a3 t2012
## Phi g1 c2009 a4 t2013
## Phi g1 c2009 a5 t2014
## Phi g1 c2009 a6 t2015
## Phi g1 c2009 a7 t2016
## Phi g1 c2010 a0 t2010
## Phi g1 c2010 a1 t2011
## Phi g1 c2010 a2 t2012
## Phi g1 c2010 a3 t2013
## Phi g1 c2010 a4 t2014
## Phi g1 c2010 a5 t2015
## Phi g1 c2010 a6 t2016
## Phi g1 c2011 a0 t2011
## Phi g1 c2011 a1 t2012
## Phi g1 c2011 a2 t2013
## Phi g1 c2011 a3 t2014
## Phi g1 c2011 a4 t2015
## Phi g1 c2011 a5 t2016
## Phi g1 c2012 a0 t2012
## Phi g1 c2012 a1 t2013
## Phi g1 c2012 a2 t2014
## Phi g1 c2012 a3 t2015
## Phi g1 c2012 a4 t2016
## Phi g1 c2013 a0 t2013
## Phi g1 c2013 a1 t2014
## Phi g1 c2013 a2 t2015
## Phi g1 c2013 a3 t2016
## Phi g1 c2014 a0 t2014
## Phi g1 c2014 a1 t2015
## Phi g1 c2014 a2 t2016
## Phi g1 c2015 a0 t2015
## Phi g1 c2015 a1 t2016
## Phi g1 c2016 a0 t2016
## p g1 c2006 a1 t2007
Phi.timeage.P.dot$results$beta # parameter estimates on logit scale
## estimate se lcl ucl
## Phi:(Intercept) 0.7564328 0.7807013 -0.7737418 2.2866075
## Phi:time2007 -1.4433637 0.8675556 -3.1437728 0.2570453
## Phi:time2008 -0.1747138 0.8230125 -1.7878184 1.4383908
## Phi:time2009 -0.1604063 0.8288442 -1.7849410 1.4641283
## Phi:time2010 -1.5109004 0.8081384 -3.0948518 0.0730510
## Phi:time2011 0.1469165 1.0519597 -1.9149245 2.2087575
## Phi:time2012 -1.1201940 0.9027083 -2.8895024 0.6491144
## Phi:time2013 14.7977310 0.0000000 14.7977310 14.7977310
## Phi:time2014 -25.3638640 333.9256600 -679.8581700 629.1304400
## Phi:time2015 -2.8014329 1.4846581 -5.7113628 0.1084971
## Phi:time2016 12.8482220 1093.5212000 -2130.4534000 2156.1499000
## Phi:age1 -0.1207040 0.4254795 -0.9546439 0.7132358
## Phi:age2 45.5827580 0.0000000 45.5827580 45.5827580
## Phi:age3 1.0131786 0.5695258 -0.1030919 2.1294492
## Phi:age4 0.5471955 0.6524114 -0.7315308 1.8259217
## Phi:age5 24.3098730 333.9263500 -630.1857800 678.8055300
## Phi:age6 24.9994900 333.9263500 -629.4961700 679.4951500
## Phi:age7 54.0310250 0.0000000 54.0310250 54.0310250
## Phi:age8 0.5917444 1.7290288 -2.7971521 3.9806410
## Phi:age9 12.6406500 0.0000000 12.6406500 12.6406500
## Phi:age10 -0.1700383 0.0000000 -0.1700383 -0.1700383
## p:(Intercept) -0.1936090 0.1352251 -0.4586502 0.0714322
We can see that some survival rates are estimated at 1 with a standard error of 0. Others have a huge standard error and a confidence interval that goes from 0 to 1. On the logit scale, we can see some very large estimates in absolute terms (anything < -10 on the logit scale means essentially 0 on the probability scale and anything > 10 on the logit scale means essentially 1 on the probability scale) with enormous standard errors. We can also see some estimates with standard errors of 0. All these estimates indicate estimation problems and it looks like our model is perhaps too complex for the data (i.e. we don’t have enough data to estimate all these parameters properly).
Plotting the survival estimates for the first two cohorts shows the problems more clearly:
# plotting results
plot(1:11, Phi.timeage.P.dot$results$real$estimate[1:11], type='l', las=1, ylab = "Annual survival", xlab = "Year of study")
lines(2:11, Phi.timeage.P.dot$results$real$estimate[12:21], col='red')
Perhaps we have been a bit too ambitious about age: we tried to estimate a separate survival rate for each year of age but we have few individuals that reached more than a few years of age.
So we need to revisit the age variable. Often, young animals have a lower survival rate than older ones. Immature individuals (hadedas start to breed at around 2 at the earliest, as far as we know) often move around a lot and are more likely to emigrate from the study area than older individuals. Since we are only able to estimate ‘local’ or ‘apparent’ survival (i.e. the product of the true survival rate and the probability of staying in the study area), we also should allow this age class to have a different survival rate.
For our next model, we would therefore like to have three age classes: “juvenile” (0-1 year old), “immature” (1-2 years old) and “adult” (2 years or older). We need to add a new age variable to the design data:
head(hade.ddl$Phi)
## par.index model.index group cohort age time occ.cohort Cohort Age Time
## 1 1 1 1 2006 0 2006 1 0 0 0
## 2 2 2 1 2006 1 2007 1 0 1 1
## 3 3 3 1 2006 2 2008 1 0 2 2
## 4 4 4 1 2006 3 2009 1 0 3 3
## 5 5 5 1 2006 4 2010 1 0 4 4
## 6 6 6 1 2006 5 2011 1 0 5 5
ageclass <- NA
ageclass[hade.ddl$Phi$Age == 0] <- "juvenile"
ageclass[hade.ddl$Phi$Age ==1] <- "immature"
ageclass[hade.ddl$Phi$Age > 1] <- "adult"
hade.ddl$Phi$ageclass <- as.factor(ageclass)
head(hade.ddl$Phi, n=15)
## par.index model.index group cohort age time occ.cohort Cohort Age Time
## 1 1 1 1 2006 0 2006 1 0 0 0
## 2 2 2 1 2006 1 2007 1 0 1 1
## 3 3 3 1 2006 2 2008 1 0 2 2
## 4 4 4 1 2006 3 2009 1 0 3 3
## 5 5 5 1 2006 4 2010 1 0 4 4
## 6 6 6 1 2006 5 2011 1 0 5 5
## 7 7 7 1 2006 6 2012 1 0 6 6
## 8 8 8 1 2006 7 2013 1 0 7 7
## 9 9 9 1 2006 8 2014 1 0 8 8
## 10 10 10 1 2006 9 2015 1 0 9 9
## 11 11 11 1 2006 10 2016 1 0 10 10
## 12 12 12 1 2007 0 2007 2 1 0 1
## 13 13 13 1 2007 1 2008 2 1 1 2
## 14 14 14 1 2007 2 2009 2 1 2 3
## 15 15 15 1 2007 3 2010 2 1 3 4
## ageclass
## 1 juvenile
## 2 immature
## 3 adult
## 4 adult
## 5 adult
## 6 adult
## 7 adult
## 8 adult
## 9 adult
## 10 adult
## 11 adult
## 12 juvenile
## 13 immature
## 14 adult
## 15 adult
Then set up to new model structures that make use of this new age variable, and then run two extra models:
ageclass <- list(formula=~ageclass)
timeageclass <- list(formula=~time + ageclass)
Phi.ageclass.P.dot <- mark(hade.process,hade.ddl,model.parameters=list(Phi=ageclass, p=dot))
##
## Output summary for CJS model
## Name : Phi(~ageclass)p(~1)
##
## Npar : 4
## -2lnL: 941.7753
## AICc : 949.8681
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 1.2463452 0.1952534 0.8636486 1.6290419
## Phi:ageclassimmature -1.1532483 0.3757420 -1.8897027 -0.4167939
## Phi:ageclassjuvenile -1.2191000 0.2836041 -1.7749640 -0.6632360
## p:(Intercept) -0.1666705 0.1423216 -0.4456208 0.1122798
##
##
## Real Parameter Phi
##
## 2006 2007 2008 2009 2010 2011 2012
## 2006 0.5068109 0.5232574 0.7766666 0.7766666 0.7766666 0.7766666 0.7766666
## 2007 0.5068109 0.5232574 0.7766666 0.7766666 0.7766666 0.7766666
## 2008 0.5068109 0.5232574 0.7766666 0.7766666 0.7766666
## 2009 0.5068109 0.5232574 0.7766666 0.7766666
## 2010 0.5068109 0.5232574 0.7766666
## 2011 0.5068109 0.5232574
## 2012 0.5068109
## 2013
## 2014
## 2015
## 2016
## 2013 2014 2015 2016
## 2006 0.7766666 0.7766666 0.7766666 0.7766666
## 2007 0.7766666 0.7766666 0.7766666 0.7766666
## 2008 0.7766666 0.7766666 0.7766666 0.7766666
## 2009 0.7766666 0.7766666 0.7766666 0.7766666
## 2010 0.7766666 0.7766666 0.7766666 0.7766666
## 2011 0.7766666 0.7766666 0.7766666 0.7766666
## 2012 0.5232574 0.7766666 0.7766666 0.7766666
## 2013 0.5068109 0.5232574 0.7766666 0.7766666
## 2014 0.5068109 0.5232574 0.7766666
## 2015 0.5068109 0.5232574
## 2016 0.5068109
##
##
## Real Parameter p
##
## 2007 2008 2009 2010 2011 2012 2013
## 2006 0.4584286 0.4584286 0.4584286 0.4584286 0.4584286 0.4584286 0.4584286
## 2007 0.4584286 0.4584286 0.4584286 0.4584286 0.4584286 0.4584286
## 2008 0.4584286 0.4584286 0.4584286 0.4584286 0.4584286
## 2009 0.4584286 0.4584286 0.4584286 0.4584286
## 2010 0.4584286 0.4584286 0.4584286
## 2011 0.4584286 0.4584286
## 2012 0.4584286
## 2013
## 2014
## 2015
## 2016
## 2014 2015 2016 2017
## 2006 0.4584286 0.4584286 0.4584286 0.4584286
## 2007 0.4584286 0.4584286 0.4584286 0.4584286
## 2008 0.4584286 0.4584286 0.4584286 0.4584286
## 2009 0.4584286 0.4584286 0.4584286 0.4584286
## 2010 0.4584286 0.4584286 0.4584286 0.4584286
## 2011 0.4584286 0.4584286 0.4584286 0.4584286
## 2012 0.4584286 0.4584286 0.4584286 0.4584286
## 2013 0.4584286 0.4584286 0.4584286 0.4584286
## 2014 0.4584286 0.4584286 0.4584286
## 2015 0.4584286 0.4584286
## 2016 0.4584286
Phi.timeageclass.P.dot <- mark(hade.process,hade.ddl,model.parameters=list(Phi=timeageclass, p=dot))
##
## Note: only 12 parameters counted of 14 specified parameters
## AICc and parameter count have been adjusted upward
##
## Output summary for CJS model
## Name : Phi(~time + ageclass)p(~1)
##
## Npar : 14 (unadjusted=12)
## -2lnL: 918.8078
## AICc : 947.8054 (unadjusted=943.54535)
##
## Beta
## estimate se lcl ucl
## Phi:(Intercept) 2.3229658 0.8397174 0.6771197 3.9688120
## Phi:time2007 -1.4632642 0.8727986 -3.1739495 0.2474211
## Phi:time2008 -0.0406714 0.8469205 -1.7006356 1.6192928
## Phi:time2009 -0.1381503 0.8385447 -1.7816980 1.5053973
## Phi:time2010 -1.2773726 0.8265514 -2.8974134 0.3426682
## Phi:time2011 -0.6270499 1.0335261 -2.6527610 1.3986613
## Phi:time2012 -1.5135325 0.9254581 -3.3274304 0.3003654
## Phi:time2013 13.0047720 2297.1134000 -4489.3375000 4515.3470000
## Phi:time2014 -1.5419169 1.0304345 -3.5615685 0.4777346
## Phi:time2015 -2.2376788 1.0368213 -4.2698486 -0.2055089
## Phi:time2016 15.9741370 3000.8475000 -5865.6871000 5897.6354000
## Phi:ageclassimmature -1.7145926 0.4636716 -2.6233889 -0.8057963
## Phi:ageclassjuvenile -1.5535709 0.4192323 -2.3752663 -0.7318756
## p:(Intercept) -0.1815968 0.1417470 -0.4594210 0.0962273
##
##
## Real Parameter Phi
##
## 2006 2007 2008 2009 2010 2011 2012
## 2006 0.68339 0.2984079 0.9074000 0.8988776 0.7399278 0.8450006 0.6919887
## 2007 0.3331729 0.6382327 0.8988776 0.7399278 0.8450006 0.6919887
## 2008 0.6745251 0.6154365 0.7399278 0.8450006 0.6919887
## 2009 0.6527716 0.3387209 0.8450006 0.6919887
## 2010 0.3756677 0.4953310 0.6919887
## 2011 0.5355263 0.2879914
## 2012 0.3221000
## 2013
## 2014
## 2015
## 2016
## 2013 2014 2015 2016
## 2006 0.9999998 0.6859061 0.5213089 1.0000000
## 2007 0.9999998 0.6859061 0.5213089 1.0000000
## 2008 0.9999998 0.6859061 0.5213089 1.0000000
## 2009 0.9999998 0.6859061 0.5213089 1.0000000
## 2010 0.9999998 0.6859061 0.5213089 1.0000000
## 2011 0.9999998 0.6859061 0.5213089 1.0000000
## 2012 0.9999988 0.6859061 0.5213089 1.0000000
## 2013 0.9999990 0.2822063 0.5213089 1.0000000
## 2014 0.3159338 0.1639255 1.0000000
## 2015 0.1872036 0.9999999
## 2016 0.9999999
##
##
## Real Parameter p
##
## 2007 2008 2009 2010 2011 2012 2013
## 2006 0.4547251 0.4547251 0.4547251 0.4547251 0.4547251 0.4547251 0.4547251
## 2007 0.4547251 0.4547251 0.4547251 0.4547251 0.4547251 0.4547251
## 2008 0.4547251 0.4547251 0.4547251 0.4547251 0.4547251
## 2009 0.4547251 0.4547251 0.4547251 0.4547251
## 2010 0.4547251 0.4547251 0.4547251
## 2011 0.4547251 0.4547251
## 2012 0.4547251
## 2013
## 2014
## 2015
## 2016
## 2014 2015 2016 2017
## 2006 0.4547251 0.4547251 0.4547251 0.4547251
## 2007 0.4547251 0.4547251 0.4547251 0.4547251
## 2008 0.4547251 0.4547251 0.4547251 0.4547251
## 2009 0.4547251 0.4547251 0.4547251 0.4547251
## 2010 0.4547251 0.4547251 0.4547251 0.4547251
## 2011 0.4547251 0.4547251 0.4547251 0.4547251
## 2012 0.4547251 0.4547251 0.4547251 0.4547251
## 2013 0.4547251 0.4547251 0.4547251 0.4547251
## 2014 0.4547251 0.4547251 0.4547251
## 2015 0.4547251 0.4547251
## 2016 0.4547251
Look at the update model selection table:
# model selection
hade.cjs.results <- collect.models()
hade.cjs.results
## model npar AICc DeltaAICc weight Deviance
## 8 Phi(~time + ageclass)p(~1) 14 947.8054 0.000000 5.936613e-01 265.5811
## 3 Phi(~ageclass)p(~1) 4 949.8681 2.062723 2.116529e-01 288.5486
## 6 Phi(~time + age)p(~1) 22 950.0441 2.238758 1.938200e-01 250.3671
## 1 Phi(~age)p(~1) 12 961.5801 13.774764 6.058802e-04 283.6159
## 7 Phi(~time + age)p(~time) 32 963.4889 15.683540 2.332928e-04 241.0215
## 4 Phi(~1)p(~1) 2 968.9553 21.149909 1.516681e-05 311.7009
## 2 Phi(~age)p(~time) 22 969.6318 21.826448 1.081398e-05 269.9548
## 5 Phi(~time)p(~time) 21 975.1827 27.377339 6.739459e-07 277.7241
We can see that the best model is now the one with age class and time effects on survival.
Let’s have a look at the parameter estimates and plot the ones for the first cohort, this time with confidence intervals:
# look at the parameter estimates
Phi.timeageclass.P.dot$results$real # parameter estimates
## estimate se lcl ucl fixed
## Phi g1 c2006 a0 t2006 0.6833900 1.711638e-01 3.140795e-01 0.9105122
## Phi g1 c2006 a1 t2007 0.2984079 9.152490e-02 1.529391e-01 0.5004877
## Phi g1 c2006 a2 t2008 0.9074000 4.426220e-02 7.772687e-01 0.9649320
## Phi g1 c2006 a3 t2009 0.8988776 4.644160e-02 7.655581e-01 0.9603129
## Phi g1 c2006 a4 t2010 0.7399278 8.315860e-02 5.494874e-01 0.8690502
## Phi g1 c2006 a5 t2011 0.8450006 9.039040e-02 5.849791e-01 0.9547214
## Phi g1 c2006 a6 t2012 0.6919887 9.855110e-02 4.758127e-01 0.8475730
## Phi g1 c2006 a7 t2013 0.9999998 5.063260e-04 4.473163e-298 1.0000000
## Phi g1 c2006 a8 t2014 0.6859061 1.345752e-01 3.909604e-01 0.8813602
## Phi g1 c2006 a9 t2015 0.5213089 1.555544e-01 2.429675e-01 0.7870195
## Phi g1 c2006 a10 t2016 1.0000000 3.395574e-05 8.713513e-297 1.0000000
## Phi g1 c2007 a0 t2007 0.3331729 7.175460e-02 2.096730e-01 0.4847944
## Phi g1 c2007 a1 t2008 0.6382327 1.048223e-01 4.201601e-01 0.8111531
## Phi g1 c2008 a0 t2008 0.6745251 1.063283e-01 4.450860e-01 0.8426385
## Phi g1 c2008 a1 t2009 0.6154365 1.085291e-01 3.944722e-01 0.7972182
## Phi g1 c2009 a0 t2009 0.6527716 9.520550e-02 4.521408e-01 0.8106921
## Phi g1 c2009 a1 t2010 0.3387209 8.357390e-02 1.977675e-01 0.5155719
## Phi g1 c2010 a0 t2010 0.3756677 9.800770e-02 2.096557e-01 0.5771406
## Phi g1 c2010 a1 t2011 0.4953310 1.817767e-01 1.909375e-01 0.8032255
## Phi g1 c2011 a0 t2011 0.5355263 1.889593e-01 2.064246e-01 0.8363481
## Phi g1 c2011 a1 t2012 0.2879914 1.281584e-01 1.061993e-01 0.5792870
## Phi g1 c2012 a0 t2012 0.3221000 1.196184e-01 1.396892e-01 0.5816617
## Phi g1 c2012 a1 t2013 0.9999988 2.812300e-03 8.053351e-299 1.0000000
## Phi g1 c2013 a0 t2013 0.9999990 2.394100e-03 9.460355e-299 1.0000000
## Phi g1 c2013 a1 t2014 0.2822063 1.506892e-01 8.381610e-02 0.6282004
## Phi g1 c2014 a0 t2014 0.3159338 1.627296e-01 9.549400e-02 0.6689149
## Phi g1 c2014 a1 t2015 0.1639255 1.054251e-01 4.160720e-02 0.4696305
## Phi g1 c2015 a0 t2015 0.1872036 1.126080e-01 5.123090e-02 0.4955637
## Phi g1 c2015 a1 t2016 0.9999999 1.886042e-04 1.568755e-297 1.0000000
## Phi g1 c2016 a0 t2016 0.9999999 1.605538e-04 1.842833e-297 1.0000000
## p g1 c2006 a1 t2007 0.4547251 3.514620e-02 3.871232e-01 0.5240383
## note
## Phi g1 c2006 a0 t2006
## Phi g1 c2006 a1 t2007
## Phi g1 c2006 a2 t2008
## Phi g1 c2006 a3 t2009
## Phi g1 c2006 a4 t2010
## Phi g1 c2006 a5 t2011
## Phi g1 c2006 a6 t2012
## Phi g1 c2006 a7 t2013
## Phi g1 c2006 a8 t2014
## Phi g1 c2006 a9 t2015
## Phi g1 c2006 a10 t2016
## Phi g1 c2007 a0 t2007
## Phi g1 c2007 a1 t2008
## Phi g1 c2008 a0 t2008
## Phi g1 c2008 a1 t2009
## Phi g1 c2009 a0 t2009
## Phi g1 c2009 a1 t2010
## Phi g1 c2010 a0 t2010
## Phi g1 c2010 a1 t2011
## Phi g1 c2011 a0 t2011
## Phi g1 c2011 a1 t2012
## Phi g1 c2012 a0 t2012
## Phi g1 c2012 a1 t2013
## Phi g1 c2013 a0 t2013
## Phi g1 c2013 a1 t2014
## Phi g1 c2014 a0 t2014
## Phi g1 c2014 a1 t2015
## Phi g1 c2015 a0 t2015
## Phi g1 c2015 a1 t2016
## Phi g1 c2016 a0 t2016
## p g1 c2006 a1 t2007
Phi.timeageclass.P.dot$results$beta # parameter estimates on logit scale
## estimate se lcl ucl
## Phi:(Intercept) 2.3229658 0.8397174 0.6771197 3.9688120
## Phi:time2007 -1.4632642 0.8727986 -3.1739495 0.2474211
## Phi:time2008 -0.0406714 0.8469205 -1.7006356 1.6192928
## Phi:time2009 -0.1381503 0.8385447 -1.7816980 1.5053973
## Phi:time2010 -1.2773726 0.8265514 -2.8974134 0.3426682
## Phi:time2011 -0.6270499 1.0335261 -2.6527610 1.3986613
## Phi:time2012 -1.5135325 0.9254581 -3.3274304 0.3003654
## Phi:time2013 13.0047720 2297.1134000 -4489.3375000 4515.3470000
## Phi:time2014 -1.5419169 1.0304345 -3.5615685 0.4777346
## Phi:time2015 -2.2376788 1.0368213 -4.2698486 -0.2055089
## Phi:time2016 15.9741370 3000.8475000 -5865.6871000 5897.6354000
## Phi:ageclassimmature -1.7145926 0.4636716 -2.6233889 -0.8057963
## Phi:ageclassjuvenile -1.5535709 0.4192323 -2.3752663 -0.7318756
## p:(Intercept) -0.1815968 0.1417470 -0.4594210 0.0962273
# plotting results
plot(1:11, Phi.timeageclass.P.dot$results$real$estimate[1:11], type='l', las=1, ylab = "Annual survival", xlab = "Year of study", ylim=c(0,1))
for (i in 1:11) segments(i, Phi.timeageclass.P.dot$results$real$lcl[i],y1=Phi.timeageclass.P.dot$results$real$ucl[i])
This looks a lot better, even though two survival probabilities (for 2013 and 2016) are still poorly estimated. Since we had little data after 2012 or so, this is perhaps not too surprising. We need to interpret the estimates for these later years very carefully.
Overall, it looks like juvenile survival was roughly around 0.7, immature survival around 0.3 (remembering that this probably reflects a lot of permanent emigration rather than truely low survival) and adult survival around 0.8-0.9. That gels well with what we would expect for this type of bird.
Goodness of fit
But how good is our model, really? The modelling approach we use here makes several assumptions. Most importantly:
- That every marked individual in the population at time t has the same survival rate from t to t+1.
- That every marked individual in the population at time t has the same recapture probability at that time.
- That no marks are lost or misread.
- That the process of capturing, marking and releasing individuals is short relative to the survival interval.
The last assumption is clearly violated in our case: we receive resightings on a continuous basis and have lumped them into yearly bins for this analysis. It turns out that survival estimates are quite robust to this kind of lumping even though in our case, it is perhaps a bit extreme and we might want to explore what happens if we restrict our analysis to include only resightings that were made during part of the year.
Assumption 3 is probably realistic. We are not aware that the birds lose rings and we kept some spare colour rings exposed to full sun for many years without noticing any fading of the colour or any other problems. The rings are quite conspicuous and display the 2-letter code clearly. Often, people also send us photos. Mis-identification is therefore also not likely to be a problem in our case.
Assumptions 1 and 2 can be tested with classical goodness-of-fit tests for the fully time-dependent CJS model. I.e. in our case, it tests model ‘Phi.time.P.time’. We know that is not our best model but if that model fits reasonably well, we know our best model is also ok.
require(R2ucare) # load U-CARE
ch.gof <- as.matrix(ch) # capture histories need to be in matrix format for this procedure
freq <- rep(1, length=dim(ch.gof)[1]) # each capture history refers to 1 individual, therefore freq = 1 for all rows
overall_CJS(ch.gof,freq)
## chi2 degree_of_freedom p_value
## Gof test for CJS model: 83.231 31 0
This statistic tests the null hypothesis that the data conform to the assumptions. A small p value indicates evidence against H_0 and therefore lack of fit. So the result above indicates that all is not well with our data. It is an overall test, pooling several subtests.
To diagnose the problem, we need to look at the individual subtests separately:
test3sr(ch.gof,freq)
## $test3sr
## stat df p_val sign_test
## 17.894 10.000 0.057 2.513
##
## $details
## component stat p_val signed_test test_perf
## 1 2 0.182 0.67 -0.427 Fisher
## 2 3 4.027 0.045 2.007 Chi-square
## 3 4 1.728 0.189 1.315 Chi-square
## 4 5 1.357 0.244 1.165 Chi-square
## 5 6 0.351 0.553 0.592 Fisher
## 6 7 0.144 0.705 0.379 Chi-square
## 7 8 8.968 0.003 2.995 Fisher
## 8 9 0 1 0 Fisher
## 9 10 0.51 0.475 0.714 Fisher
## 10 11 0.627 0.429 -0.792 Fisher
test3sm(ch.gof,freq)
## $test3sm
## stat df p_val
## 5.614 6.000 0.468
##
## $details
## component stat df p_val test_perf
## 1 2 0.936 1 0.333 Fisher
## 2 3 3.155 1 0.076 Chi-square
## 3 4 1.129 1 0.288 Chi-square
## 4 5 0 1 1 Fisher
## 5 6 0 1 1 Fisher
## 6 7 0.394 1 0.53 Fisher
## 7 8 0 0 0 None
## 8 9 0 0 0 None
## 9 10 0 0 0 None
## 10 11 0 0 0 None
test2ct(ch.gof,freq)
## $test2ct
## stat df p_val sign_test
## 53.996 9.000 0.000 -6.705
##
## $details
## component dof stat p_val signed_test test_perf
## 1 2 1 3.452 0.063 -1.858 Fisher
## 2 3 1 8.155 0.004 -2.856 Chi-square
## 3 4 1 9.525 0.002 -3.086 Chi-square
## 4 5 1 3.109 0.078 -1.763 Chi-square
## 5 6 1 13.493 0 -3.673 Chi-square
## 6 7 1 6.238 0.013 -2.498 Chi-square
## 7 8 1 7.025 0.008 -2.65 Chi-square
## 8 9 1 0 1 0 Fisher
## 9 10 1 2.999 0.083 -1.732 Fisher
test2cl(ch.gof,freq)
## $test2cl
## stat df p_val
## 5.727 6.000 0.454
##
## $details
## component dof stat p_val test_perf
## 1 2 1 1.913 0.167 Fisher
## 2 3 1 0.748 0.387 Chi-square
## 3 4 1 1.269 0.26 Chi-square
## 4 5 1 0.64 0.424 Chi-square
## 5 6 1 0 1 Fisher
## 6 7 1 1.157 0.282 Fisher
## 7 8 0 0 0 None
## 8 9 0 0 0 None
The Tests 3sr and 3sm roughly deal with the first assumption, i.e. homogeneous survival probabilities. These tests look ok. Even though we know that we have rather strong age effect and these are not accounted for in the GOF test, these tests do not provide evidence for strong heterogeneity in survival.
Test 2ct and 2cl roughly test the second assumption, i.e. homogeneous recapture probabilities. Test 2ct in particular returns a very small P value, indicating that our assumption of homogeneous recapture is violated. This probably has to do with spatial heterogeneity: while we ringed birds over a large area, most resighting effort was concentrated in a relatively small area around Constantia. The next step with this analysis would need to take spatial heterogeneity in recapture into account.
It is also always useful to look at the raw data. A useful summary is the so-called M-array. The first column (a separate vector in the output below) shows the number of individuals ‘released’ each year. This includes all newly ringed birds but also those that were resighted (in traditional capture-mark-recapture experiments they would have been recaptured and re-released). Then, the matrix shows the how many of these individuals were re-sighted again for the first time in subsequent years. And the final vector shows how many individuals were never seen again.
E.g. from the output below, we can see that we released 27 birds during the first year of our study. Of those, 8 were resighted for the first time in the following year, one was resighted for the first time in year 3 and two more in year 4. Sixteen were never seen again.
marray(ch.gof, freq)
## $R
## [,1]
## [1,] 27
## [2,] 73
## [3,] 59
## [4,] 95
## [5,] 65
## [6,] 29
## [7,] 29
## [8,] 16
## [9,] 20
## [10,] 16
## [11,] 7
##
## $m
## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 8 0 1 2 0 0 0 0 0 0 0
## [2,] 0 12 0 2 1 0 1 0 1 0 1
## [3,] 0 0 19 3 3 2 0 1 1 0 0
## [4,] 0 0 0 32 5 0 0 2 1 0 0
## [5,] 0 0 0 0 17 2 1 3 0 0 0
## [6,] 0 0 0 0 0 15 1 1 0 0 0
## [7,] 0 0 0 0 0 0 9 1 3 0 0
## [8,] 0 0 0 0 0 0 0 11 0 0 0
## [9,] 0 0 0 0 0 0 0 0 8 0 1
## [10,] 0 0 0 0 0 0 0 0 0 6 1
## [11,] 0 0 0 0 0 0 0 0 0 0 3
##
##
## $never
## [,1]
## [1,] 16
## [2,] 55
## [3,] 30
## [4,] 55
## [5,] 42
## [6,] 12
## [7,] 16
## [8,] 5
## [9,] 11
## [10,] 9
## [11,] 4
If recapture probabilities are fairly homogeneous, we would expect to see the highest numbers on the diagonal (it is most likely to see individuals for the first time during the next occasion) and then drop. That seems to be the case in our data. We can clearly see the drop in sample sizes towards the end of the study. But there are no other obvious prolems visible.
We can calculate the ‘overdispersion factor’. The assumptions above imply that our data should follow a certain multivariated distribution. Think of flipping coins and counting heads and tails. In reality, no data set really follows that kind of distribution. Individual’s fates might not be independent or they vary in their probability to be resighted or to survive. These kinds of things cause extra variation that our model doesn’t account for. Our parameter estimates may therefore have confidence intervals that are too narrow. Estimating overdispersion is therefore another way of examining goodness of fit. Dividing the \(\chi^2\) value from our GOF test by the degrees of freedom gives us an estimate of overdispersion. It should be close to 1. If it is a lot larger (e.g. >3), then that can indicate that our model does not describe the structure in our data well.
# measure of overdispersion
chi2 <- overall_CJS(ch.gof,freq)$chi2 # test statistic
df <- overall_CJS(ch.gof,freq)$degree_of_freedom # degrees of freedom
chi2 / df
## [1] 2.684871
Overdispersion is not terrible in this data set and we could consider adjusting our analysis for it. This would imply wider confidence intervals and a tendency to favour simpler models. However, the downside is loss of power and in this case I would certainly try to find a model structure that fits the data better. One option would be to use a multi-state model where the states are different areas. One could then allow the resighting probability to vary among areas.