The R2ucare package contains R functions to perform goodness-of-fit tests for capture-recapture models. It also has various functions to manipulate capture-recapture data.
First things first, load the R2ucare package:
There are 3 main data formats when manipulating capture-recapture data, corresponding to the 3 main computer software available to fit corresponding models: RMark, E-SURGE and Mark. With R2ucare, it is easy to work with any of these formats. We will use the classical dipper dataset, which is provided with the package (thanks to Gilbert Marzolin for sharing his data).
RMark files
# # read in text file as described at pages 50-51 in http://www.phidot.org/software/mark/docs/book/pdf/app_3.pdf
library(RMark)
data_path <- "/Library/Frameworks/R.framework/Versions/4.1/Resources/library/RMark/extdata/"
dipper <- import.chdata(paste0(data_path, "dipper.txt"),
field.names = c("ch", "sex"), header = FALSE)
dipper <- as.data.frame(table(dipper))
str(dipper)## 'data.frame': 64 obs. of 3 variables:
## $ ch : Factor w/ 32 levels "0000001","0000010",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ sex : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
## $ Freq: int 22 12 11 7 6 6 10 1 1 5 ...Get encounter histories, counts and groups:
dip.hist <- matrix(as.numeric(unlist(strsplit(as.character(dipper$ch),""))),
nrow = length(dipper$ch),
byrow = T)
dip.freq <- dipper$Freq
dip.group <- dipper$sex
head(dip.hist)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 0 0 0 0 0 1
## [2,] 0 0 0 0 0 1 0
## [3,] 0 0 0 0 0 1 1
## [4,] 0 0 0 0 1 0 0
## [5,] 0 0 0 0 1 1 0
## [6,] 0 0 0 0 1 1 1
head(dip.freq)## [1] 22 12 11 7 6 6
head(dip.group)## [1] Female Female Female Female Female Female
## Levels: Female MaleE-SURGE filesLet’s use the read_headed function.
dipper <- system.file("extdata", "ed.txt", package = "R2ucare")
dipper <- read_headed(dipper)Get encounter histories, counts and groups:
dip.hist <- dipper$encounter_histories
dip.freq <- dipper$sample_size
dip.group <- dipper$groups
head(dip.hist)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1 1 1 1 1 1 0
## [2,] 1 1 1 1 0 0 0
## [3,] 1 1 0 0 0 0 0
## [4,] 1 1 0 0 0 0 0
## [5,] 1 1 0 0 0 0 0
## [6,] 1 1 0 0 0 0 0
head(dip.freq)## [1] 1 1 1 1 1 1
head(dip.group)## [1] "male" "male" "male" "male" "male" "male"Mark filesLet’s use the read_inp function.
dipper <- system.file("extdata", "ed.inp", package = "R2ucare")
dipper <- read_inp(dipper, group.df = data.frame(sex = c('Male', 'Female')))Get encounter histories, counts and groups:
dip.hist <- dipper$encounter_histories
dip.freq <- dipper$sample_size
dip.group <- dipper$groups
head(dip.hist)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1 1 1 1 1 1 0
## [2,] 1 1 1 1 0 0 0
## [3,] 1 1 0 0 0 0 0
## [4,] 1 1 0 0 0 0 0
## [5,] 1 1 0 0 0 0 0
## [6,] 1 1 0 0 0 0 0
head(dip.freq)## [1] 1 1 1 1 1 1
head(dip.group)## [1] "Male" "Male" "Male" "Male" "Male" "Male"Split the dataset in females/males:
mask <- (dip.group == 'Female')
dip.fem.hist <- dip.hist[mask,]
dip.fem.freq <- dip.freq[mask]
mask <- (dip.group == 'Male')
dip.mal.hist <- dip.hist[mask,]
dip.mal.freq <- dip.freq[mask]Tadaaaaan, now you’re ready to perform Test.3Sr, Test3.Sm, Test2.Ct and Test.Cl for females:
test3sr_females <- test3sr(dip.fem.hist, dip.fem.freq)
test3sm_females <- test3sm(dip.fem.hist, dip.fem.freq)
test2ct_females <- test2ct(dip.fem.hist, dip.fem.freq)
test2cl_females <- test2cl(dip.fem.hist, dip.fem.freq)
# display results:
test3sr_females## $test3sr
## stat df p_val sign_test
## 4.985 5.000 0.418 1.428
##
## $details
## component stat p_val signed_test test_perf
## 1 2 0.858 0.354 0.926 Chi-square
## 2 3 3.586 0.058 1.894 Chi-square
## 3 4 0.437 0.509 0.661 Chi-square
## 4 5 0.103 0.748 -0.321 Chi-square
## 5 6 0.001 0.982 0.032 Chi-square
test3sm_females## $test3sm
## stat df p_val
## 2.041 3.000 0.564
##
## $details
## component stat df p_val test_perf
## 1 2 1.542 1 0.214 Fisher
## 2 3 0 1 1 Fisher
## 3 4 0.499 1 0.48 Fisher
## 4 5 0 0 0 None
## 5 6 0 0 0 None
test2ct_females## $test2ct
## stat df p_val sign_test
## 3.250 4.000 0.517 -0.901
##
## $details
## component dof stat p_val signed_test test_perf
## 1 2 1 0 1 0 Fisher
## 2 3 1 0 1 0 Fisher
## 3 4 1 0 1 0 Fisher
## 4 5 1 3.25 0.071 -1.803 Fisher
test2cl_females## $test2cl
## stat df p_val
## 0 0 1
##
## $details
## component dof stat p_val test_perf
## 1 2 0 0 0 None
## 2 3 0 0 0 None
## 3 4 0 0 0 NoneOr perform all tests at once:
overall_CJS(dip.fem.hist, dip.fem.freq)## chi2 degree_of_freedom p_value
## Gof test for CJS model: 10.276 12 0.592What to do if these tests are significant? If you detect a transient effect, then 2 age classes should be considered on the survival probability to account for this issue, which is straightforward to do in RMark (Cooch and White 2017; appendix C) or E-SURGE (Choquet et al. 2009). If trap dependence is significant, Cooch and White (2017) illustrate how to use a time-varying individual covariate to account for this effect in RMark, while Gimenez et al. (2003) suggest the use of multistate models that can be fitted with RMark or E-SURGE, and Pradel and Sanz (2012) recommend multievent models that can be fitted in E-SURGE only.
Now how to assess the fit of a model including trap-dependence and/or transience? For example, let’s assume we detected a significant effect of trap-dependence, we accounted for it in a model, and now we’d like to know whether our efforts paid off. Because the overall statistic is the sum of the four single components (Test.3Sr, Test3.Sm, Test2.Ct and Test.Cl), we obtain a test for the model with trap-dependence as follows:
overall_test <- overall_CJS(dip.fem.hist, dip.fem.freq) # overall test
twoct_test <- test2ct(dip.fem.hist, dip.fem.freq) # test for trap-dependence
stat_tp <- overall_test$chi2 - twoct_test$test2ct["stat"] # overall stat - 2CT stat
df_tp <- overall_test$degree_of_freedom - twoct_test$test2ct["df"] # overall dof - 2CT dof
pvalue <- 1 - pchisq(stat_tp, df_tp) # compute p-value for null hypothesis:
# model with trap-dep fits the data well
pvalue## stat
## 0.5338301Read in the geese dataset. It is provided with the package (thanks to Jay Hestbeck for sharing his data).
geese <- system.file("extdata", "geese.inp", package = "R2ucare")
geese <- read_inp(geese)Get encounter histories and number of individuals with corresponding histories
geese.hist <- geese$encounter_histories
geese.freq <- geese$sample_sizeAnd now, perform Test3.GSr, Test.3.GSm, Test3G.wbwa, TestM.ITEC and TestM.LTEC:
test3Gsr(geese.hist, geese.freq)## $test3Gsr
## stat df p_val
## 117.753 12.000 0.000
##
## $details
## occasion site stat df p_val test_perf
## 1 2 1 3.894777e-03 1 9.502378e-01 Chi-square
## 2 2 2 2.715575e-04 1 9.868523e-01 Chi-square
## 3 2 3 8.129814e+00 1 4.354322e-03 Chi-square
## 4 3 1 1.139441e+01 1 7.366526e-04 Chi-square
## 5 3 2 2.707742e+00 1 9.986223e-02 Chi-square
## 6 3 3 3.345916e+01 1 7.277633e-09 Chi-square
## 7 4 1 1.060848e+01 1 1.125702e-03 Chi-square
## 8 4 2 3.533332e-01 1 5.522323e-01 Chi-square
## 9 4 3 1.016778e+01 1 1.429165e-03 Chi-square
## 10 5 1 1.101349e+01 1 9.045141e-04 Chi-square
## 11 5 2 1.292013e-01 1 7.192616e-01 Chi-square
## 12 5 3 2.978513e+01 1 4.826802e-08 Chi-square
test3Gsm(geese.hist, geese.freq)## $test3Gsm
## stat df p_val
## 302.769 119.000 0.000
##
## $details
## occasion site stat df p_val test_perf
## 1 2 1 23.913378 14 4.693795e-02 Fisher
## 2 2 2 24.810007 16 7.324561e-02 Fisher
## 3 2 3 11.231939 8 1.889004e-01 Fisher
## 4 3 1 36.521484 14 8.712879e-04 Fisher
## 5 3 2 21.365358 17 2.103727e-01 Chi-square
## 6 3 3 23.072982 10 1.048037e-02 Fisher
## 7 4 1 55.338866 8 3.793525e-09 Fisher
## 8 4 2 17.172011 11 1.028895e-01 Fisher
## 9 4 3 45.089296 10 2.095523e-06 Chi-square
## 10 5 1 9.061514 3 2.848411e-02 Chi-square
## 11 5 2 5.974357 4 2.010715e-01 Chi-square
## 12 5 3 29.217786 4 7.060092e-06 Chi-square
test3Gwbwa(geese.hist, geese.freq)## $test3Gwbwa
## stat df p_val
## 472.855 20.000 0.000
##
## $details
## occasion site stat df p_val test_perf
## 1 2 1 19.5914428 2 5.568936e-05 Chi-square
## 2 2 2 37.8676763 2 5.986026e-09 Chi-square
## 3 2 3 4.4873614 1 3.414634e-02 Fisher
## 4 3 1 80.5903050 1 2.777187e-19 Chi-square
## 5 3 2 98.7610833 4 1.805369e-20 Chi-square
## 6 3 3 0.8071348 1 3.689687e-01 Fisher
## 7 4 1 27.7054638 1 1.412632e-07 Chi-square
## 8 4 2 53.6936048 2 2.190695e-12 Chi-square
## 9 4 3 25.2931602 1 4.924519e-07 Chi-square
## 10 5 1 43.6547442 1 3.917264e-11 Chi-square
## 11 5 2 50.9264976 2 8.738795e-12 Chi-square
## 12 5 3 29.4763896 2 3.974507e-07 Chi-square
testMitec(geese.hist, geese.freq)## $testMitec
## stat df p_val
## 68.222 27.000 0.000
##
## $details
## occasion stat df p_val test_perf
## 1 2 14.26651 9 0.1131556176 Chi-square
## 2 3 30.83749 9 0.0003156442 Chi-square
## 3 4 23.11762 9 0.0059375004 Chi-square
testMltec(geese.hist, geese.freq)## $testMltec
## stat df p_val
## 20.990 19.000 0.337
##
## $details
## occasion stat df p_val test_perf
## 1 2 14.103805 10 0.1683092 Chi-square
## 2 3 6.886466 9 0.6489387 Chi-squareOr all tests at once:
overall_JMV(geese.hist, geese.freq)## chi2 degree_of_freedom p_value
## Gof test for JMV model: 982.592 197 0If trap-dependence or transience is significant, you may account for these lacks of fit as in the Cormack-Jolly-Seber model example. If there are signs of a memory effect, it gets a bit trickier but you may fit a model to account for this issue using hidden Markov models (also known as multievent models when applied to capture-recapture data).
Several U-CARE functions to manipulate and process capture-recapture data can be mimicked with R built-in functions. For example, recoding multisite/state encounter histories in single-site/state ones is easy:
# Assuming the geese dataset has been read in R (see above):
geese.hist[geese.hist > 1] <- 1So is recoding states:
# Assuming the geese dataset has been read in R (see above):
geese.hist[geese.hist == 3] <- 2 # all 3s become 2sAlso, reversing time is not that difficult:
What about cleaning data, i.e. deleting empty histories, and histories with no individuals?
# Assuming the female dipper dataset has been read in R (see above):
mask = (apply(dip.fem.hist, 1, sum) > 0 & dip.fem.freq > 0) # select non-empty histories, and histories with at least one individual
sum(!mask) # how many histories are to be dropped?
dip.fem.hist[mask,] # drop these histories from dataset
dip.fem.freq[mask] # from countsSelecting or gathering occasions is as simple as that:
# Assuming the female dipper dataset has been read in R (see above):
dip.fem.hist[, c(1,4,6)] # pick occasions 1, 4 and 6 (might be a good idea to clean the resulting dataset)
gather_146 <- apply(dip.fem.hist[,c(1,4,6)], 1, max) # gather occasions 1, 4 and 6 by taking the max
dip.fem.hist[,1] <- gather_146 # replace occasion 1 by new occasion
dip.fem.hist <- dip.fem.hist[, -c(4,6)] # drop occasions 4 and 6Finally, suppressing the first encounter is achieved as follows:
# Assuming the geese dataset has been read in R (see above):
for (i in 1:nrow(geese.hist)){
occasion_first_encounter <- min(which(geese.hist[i,] != 0)) # look for occasion of first encounter
geese.hist[i, occasion_first_encounter] <- 0 # replace the first non zero by a zero
}
# delete empty histories from the new dataset
mask <- (apply(geese.hist, 1, sum) > 0) # select non-empty histories
sum(!mask) # how many histories are to be dropped?
geese.hist[mask,] # drop these histories from dataset
geese.freq[mask] # from countsBesides these simple manipulations, several useful U-CARE functions needed a proper R equivalent. I coded a few of them, see the list below and ?name-of-the-function for more details.
marray build the m-array for single-site/state capture-recapture data;multimarray build the m-array for multi-site/state capture-recapture data;group_data pool together individuals with the same encounter capture-recapture history.ungroup_data split encounter capture-recapture histories in individual ones.