---
title: "Mosa: ensity estimation and survival analysis"
author: "Joachim Marien"
date: "31 maart 2022"
output: html_document
---

<style>
body {
text-align: justify}
</style>

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## **2. Estimation of population size**

Population dynamics are essentially concerned with all processes that play a role in the determination of population size. Consequently, if we want our study to help in understanding fluctuations and outbreak phenomena in *Mastomys*, we must first be sure that also our own study population displays  variations in population densities.
Estimating population size is an issue that has kept ecologist busy since long and much of the knowledge in this field has been gathered and discussed in handbooks. Direct counting of small mammals is impossible and one has to resort to a variety of assessment techniques to obtain a reasonably accurate estimate of population numbers. Some of these methods are enumeration methods, others have a statistical basis, providing confidence interval for estimates. The choice of which estimation technique is the most appropriate is not always evident. Here, we will discuss three possibilities and their application using the MOSA data. 

### **2.1 Enumeration methods (Minimum Number Alive)**

The easiest way of getting an idea of the population size is simply to count the number of captured animals. However, it is never certain, and even rather unlikely, that all individuals that were present were effectively captured.  

#### **2.1.1 Read capture-recapture data**
```{r,,warning=F,echo=T,message=F}
rm(list=ls()) #Remove previous R sessions

###if reading from github
# install.packages("readr")
library("readr")

prim_sec1<-data.frame(read.csv(url("https://raw.githubusercontent.com/Joachimmarien/MOSA/main/C_primary_secondary_capture_3consecutive_o_1_10000b.csv"), sep=';'))
prim_sec2<-data.frame(read.csv(url("https://raw.githubusercontent.com/Joachimmarien/MOSA/main/C_primary_secondary_capture_3consecutive_o_10000-20000b.csv"), sep=';'))
prim_sec=rbind(prim_sec1,prim_sec2)

#if reading from computer
# prim_sec <- read.delim("C:/Users/joach/OneDrive - Universiteit Antwerpen/Bureaublad/PhD selection/MOSA/CMR cleaning/since Feb 2022/R-scripts/markdown/Final markdown files/mosa_04_03_2023/C_primary_secondary_capture_3consecutive_o.txt")
prim_sec_mat=prim_sec[,10:ncol(prim_sec)]#Select capture-recapture matrix
prim_sec_char=prim_sec[,1:9]#select host characteristics

timeprim_sec<-data.frame(read.csv(url("https://raw.githubusercontent.com/Joachimmarien/MOSA/main/C_primary_secondary_time_variables_3consecutive_m2.csv"), sep=';'))
# timeprim_sec <- read.delim("C:/Users/joach/OneDrive - Universiteit Antwerpen/Bureaublad/PhD selection/MOSA/CMR cleaning/since Feb 2022/R-scripts/markdown/Final markdown files/mosa_04_03_2023/C_primary_secondary_time_variables_3consecutive_m.txt")
```

#### **2.2.2 Create primary capture-recapture matrix**

```{r,,warning=F,echo=T,message=F}
prim_cap_ses=unique(timeprim_sec$prim_capt) # find all primary capture sessions
prim_cap_numb=rep(NA,length(prim_cap_ses)) # create empty vector that will be filled with the number of captures in the FOR-loop
prim_mat=matrix(0,nrow=nrow(prim_sec_mat),ncol=length(prim_cap_ses)) # create empty primary capture matrix

i=1
for (i in 1:length(prim_cap_ses)){
  trap_sess=prim_sec_mat[,which(timeprim_sec$prim_capt==prim_cap_ses[i])] #select the capture-recapture matrix of secondary trapping sessions from primary session [i]
  prim_cap_numb[i]=length(which(apply(trap_sess,1,sum)>0)) # Count number of unique capture per secondary trapping session
  prim_mat[which(apply(trap_sess,1,sum)>0),i]=1 # Fill primary capture matrix based on secondary matrix
}
prim_mat[1:5,1:5]
```

To partially overcome this problem, Krebbs (1966) introduced the **Minimum Number Alive method (MNA)** which counts all the animals whose presence is documented. This is done by adding the number of captured animals at time *t* to the number of animals that were not captured at time *t* but had been captured before and were captured again later. Several problems should be taken into consideration with this technique. First of all, it is impossible to know whether an individual has escaped capture by chance, or whether it was temporarily dispersing during the trapping session and returned later. Similarly, animals may enter the population, stay there for some time and then disappear without ever being captured. Also, animals that skip capture once accidentally and then disappear before the following trapping session will be considered 'not present' one month too early. Nevertheless, the MNA method has the advantage to be an enumeration technique resulting in an exact number, not in an estimate associated with a given probability. Therefore, the MNA method is very valuable, albeit that it always underestimates  real population size.  

```{r,,warning=F,echo=T,message=F}
MNA=matrix(0,nrow=nrow(prim_mat),ncol=length(prim_cap_ses)) # create empty primary capture matrix
i=1
for(i in 1:nrow(MNA)){
  cap_hist=prim_mat[i,] #select each trapped animal seperately
  start=which(cap_hist==1)[1]; stop=which(cap_hist==1)[length(which(cap_hist==1))] # find first and last primary capture session in which animal was trapped
  MNA[i,start:stop]=1
}
MNA_cap_numb=apply(MNA,2,sum)
MNA[1:5,1:5]
```

### **2.2 Petersen estimate and Schnabel method**
 The well-known **Petersen estimate** is a basis for most estimation methods that use CMR data. The basic idea is that the proportion of marked animals in a recapture sample is equal to the proportion of all marked animals in the whole population. After correcting for bias, this gives an estimator of population size N:
 
 $N =\frac{(M+1)(C+1)}{R+1} -1$  where
 
N= number of animals in the population at the time of marking
 
M= number of animals marked in the first sample and released

C= number of animals caught in the second sample

R=number of marked animals caught in the second sample

One of the assumptions of this method is that one has to deal with a closed population, i.e. the population size N does not change between the first and the second sample. This is clearly not the case in our CMR-study of we should consider each month as a separate sampling session. In the period between two samples, immigrant and newborn juveniles will enter the population, animals may die or emigrate. However, if we consider each month independent, then we have three samples in each month (the three secondary trapping sessions). If we assume that the proportion of new animal that enter or leave the population between two consecutive nights is so small that it comes near to zero, we have a virtually closed population and the Petersen method can be used. An extension of the Peterson estimate for circumstances with more than two samples is the **Schnabel method** where new animals are marked in   every sample, each time resulting in an increased number of marked animals in the population. Here, we used the poisson distribution  to set confidence intervals if the sum of all recaptures is <50, and the t-distribution is used if the sum of all recaptures is >=50.
```{r,,warning=F,echo=T,message=F}
#install.packages("fishmethods")#When using first time, make sure to remove the # and install the package!!!
library("fishmethods") #install and open new R-package
schnabel=rep()
i=1
prim_capt_ses=unique(timeprim_sec$prim_capt) #vector with unique primary capture sessions
schnabel=rep(NA,length(prim_capt_ses)) #empty vector for Schnabbel median estimates
schnabel_LC=rep(NA,length(prim_capt_ses))#empty vector for Schnabbel lower confidence intervals
schnabel_UC=rep(NA,length(prim_capt_ses))#empty vector for Schnabbel higher confidence intervals

for(i in 1:length(prim_capt_ses)){
  df= prim_sec_mat[,which(timeprim_sec$prim_capt==prim_capt_ses[i])] # select secondary trapping session 
  df2=df[,1]+df[,2];t1=length(which(df2>1)) # find number of animals recaptured at the second secondary session
  df3=df[,2]+df[,3]; df4=df[,1]+df[,3];df5=df[,1]+df[,3]+df[,2];
  t2=length(which(df3>1))+length(which(df4>1))-length(which(df5>2)) #the third secondary trapping session
  c=c(apply(df,2,sum)) #A vector containing the number of animals caught in each mark-recapture experiment.
  r=c(0,t1,t2) #A vector containing the number of animal recaptured in each mark-recapture experiment.
  m=c(c[1],c[2]-t1,c[3]-t2) #A vector containing the newly marked animals in each mark-recapture experiment.
  x=schnabel(catch=c,recaps=r, newmarks=m, alpha=0.10) #calculate density estimates with confidence intervals 
  schnabel[i]=x[1,1]; schnabel_LC[i]=x[1,3]; schnabel_UC[i]=x[1,4] #Fill vectors in the FOR-loop
}
```

#### **2.2.1 Time intervals between primary captures**
Here, we select the dates and time intervals between primary capture sessions. 
```{r,,warning=F,echo=T,message=F}
i=2
time=data.frame(timeprim_sec[1,])
for (i in 2:length(prim_capt_ses)){
  time=rbind(time,timeprim_sec[which(timeprim_sec$prim_capt==prim_capt_ses[i]),][1,])
}
```

#### *2.2.2 Figure with all the density estimates**
```{r,,warning=F,echo=T,message=F}
plot(time$time_interval2,prim_cap_numb,xaxt='n',xlab="Years",ylab="Number of animals",pch=19,col="black",ylim=c(0,1400),cex.lab=1.5,main="MOSA - Mastomys natalensis")
lines(time$time_interval2,prim_cap_numb,lwd=3)
lines(time$time_interval2,MNA_cap_numb,col="red")
lines(time$time_interval2,schnabel,col="blue")
# lines(time$time_interval2,schnabel_LC,col="blue",lty=2) # Confidence interval on Schnabbel
# lines(time$time_interval2,schnabel_UC,col="blue",lty=2) # Confidence interval on Schnabbel
year=seq(292,10512,by=730); years2=seq(1995,2023,by=2)
axis(1,year,years2)
legend("topright",legend=c("Captures","Minimum number alive","Schnabel"),pch=c(19,19), col=c("black","red","blue"), pt.cex=c(1, 1), cex=1)
```

### **2.3 Closed Log-linear Capture-Recapture models**
A closed model can still incorporate sources of variation among capture probabilities: a temporal effect (subscript *t*), a heterogeneity between units (subscript *h*) and a behavioural effect (subscript *b*) (Baillargeon and Rivest, 2007). A temporal effect causes the capture probabilities to vary among capture occasions; heterogeneity causes the capture probabilities to vary among units. A behavioural effect means that the first capture changes the behaviour of a unit, so the capture probability differs before and after the first capture. These sources of variation lead to eight fundamental closed population models: M0 (no source of variation), Mt, Mh, Mth, Mb, Mtb, Mbh, Mtbh. 

The analysis of data from a closed population capture-recapture experiment amounts to finding the best fitting model and estimating the population size from the chosen model based on AIC or BIC criteria. Here, we follow steps proposed by [Baillargeon and Rivest](https://www.jstatsoft.org/article/view/v019i05) in the Rcapture package.  This package uses Poisson regressions to estimate parameters in a CMR experiment. The Poisson regressions are fitted with the glm function; then the loglinear parameters are transformed into demographic parameters. The analyis is performed on the secondary trapping occasion of each primary trapping session. 

```{r,,warning=F,echo=T,message=F}
library("Rcapture") #install and open new R-package

i=1
prim_capt_ses=unique(timeprim_sec$prim_capt) #vector with unique primary capture sessions
closedp_abundance=matrix(NA,ncol=length(prim_capt_ses),nrow=12) #abundance estimates for different models
closedp_stderr=matrix(NA,ncol=length(prim_capt_ses),nrow=12)#STDEV  for different models
closedp_AIC=matrix(NA,ncol=length(prim_capt_ses),nrow=12)#AIC for different models
closedp_BIC=matrix(NA,ncol=length(prim_capt_ses),nrow=12)#BIC for different models

for(i in 1:length(prim_capt_ses)){# For each primary trapping session
  df= prim_sec_mat[,which(timeprim_sec$prim_capt==prim_capt_ses[i])] # select secondary trapping session 
  if(any(apply(df,2,sum)<2) & sum(df<11)){
  } else {#Do not estimate for primary sessions with less than 11 captures or with zero at a secondary session
    closedp_abundance[,i]=closedp(df)$results[,1] # perform closed model abundance
    closedp_stderr[,i]=closedp(df)$results[,2]
    closedp_AIC[,i]=closedp(df)$results[,5]
    closedp_BIC[,i]=closedp(df)$results[,6]
    }
}

par(mfrow=c(4,3),mar=c(1,2,1,1))
i=1
model=c("M0","Mt","Mh Ch","Mh Pois","Mh Dar","Mh Gam","Mth Ch","Mth Pois","Mth Dar","Mth Gam","Mb","Mbh")
for(i in 1:12){
  closedp_abundance[i,which(closedp_abundance[i,]>2000)]=NA #remove estimates that don't make sense
  closedp_stderr[i,which(closedp_abundance[i,]>2000)]=NA
  closedp_AIC[i,which(closedp_abundance[i,]>2000)]=NA
  closedp_BIC[i,which(closedp_abundance[i,]>2000)]=NA
  
  plot(time$time_interval2,closedp_abundance[i,],xaxt='n',xlab="Years",ylab="Number of animals",pch=1,col="black",ylim=c(0,2000),cex.lab=1.5,main=model[i])
  lines(time$time_interval2,closedp_abundance[i,],lwd=0.8,col=1)
# lines(time$time_interval2,closedp_abundance[i,]+2*closedp_stderr[i,],lwd=1,lty=2) #if error bars are needed
# lines(time$time_interval2,closedp_abundance[i,]-2*closedp_stderr[i,],lwd=1,lty=2)#if error bars are needed
}
```

The figure shows the different fitted models for the MOSA data based on closed Loglinear models. The model with the best fit was different based on AIC or BIC criteria. 

```{r,,warning=F,echo=T,message=F}
par(mfrow=c(2,1),mar=c(4,5,1,1))
boxplot.matrix(closedp_AIC, use.cols = FALSE,xaxt='n',ylab="AIC",main="AIC")
axis(1,seq(1,12,1),model)
boxplot.matrix(closedp_BIC, use.cols = FALSE,xaxt='n',ylab="BIC",main="BIC")
axis(1,seq(1,12,1),model)
```

Model selection based on AIC

```{r,,warning=F,echo=T,message=F}
closedp_AIC2=closedp_AIC[,-(which(is.na(closedp_AIC[1,])))]
x=cbind(apply(closedp_AIC2,1,median),model)
x[which(x[,1]==min(x[,1])),]
```
Model selection based on BIC
```{r,,warning=F,echo=T,message=F}
closedp_BIC2=closedp_BIC[,-(which(is.na(closedp_BIC[1,])))]
x=cbind(apply(closedp_BIC2,1,median),model)
x[which(x[,1]==min(x[,1])),]
```

Models that included a time and heterogeneity effect were selected based on AIC criteria while models based on BIC criteria selected behavioural effect. Fortunately, the different methods show very similar patterns. Since we are mainly interested in changes in densities, we actually do not have to worry too much about which technique should be used and we can rely mainly on the MNA for other analyses, simply because they are no estimates but 'real' data. 


### **2.4 Open populations (the Cormack-Jolly-Seber model)**
In contrast to the earlier discussed models, the Cormack-Jolly Seber (CJS) method can be applied to open populations. This means that is not necessary to assume that population size is constant during the whole study period; in other words birth, dead and migration do not interfere with the estimation. However, survival should be age-dependent and the problem of equal catchability shows up here as well. Open population models apply when animals are released and recaptured or resighted at future capture occasions (Baillargeon and Rivest, 2007). Typically the capture occasions are distant in time and mortality occurs between them. 

Cormack-Jolly-Seber models are mark-recapture models used to estimate two parameters:

 - Detection probability (pt, the probability of encountering a live animal at time *t*)

 - Apparent survival ($\phi$t, the probability of an animal surviving and remaining in the study area between time *t* and *t* + 1).

We will use the RMark package to perform such a survival analysis on the MOSA data in R. An excellent introduction to the package on the famous dipper data is written by [James Peterson](https://jamesepaterson.github.io/jamespatersonblog/2020-04-26_introduction_to_CJS.html) . We also refer people to which come in contact with capture-recapture/resight types of analysis in the future to the excellent and lucidly written [MARK](http://www.phidot.org/software/mark/docs/book/) handbook. 
 


#### **2.4.1. Goodness-of-fit tests for the Cormack-Jolly-Seber model**

We will calculate apparent survival and separate it from recapture probability. Note that the CJS model is conditional on the first capture.

To start the analysis and perform model selection according to the principles of the Akaike Information Criterion (AIC) we first need to verify whether the initial (full) model does indeed provides a reasonable fit to the data at hand (ie capture histories of the birds). This Goodness of Fit test (GOF) is similar to performing a check of the residuals as in a linear regression model. It is however adapted to the specific type of capture-mark-recapture data. These test can be carried out in [U-CARE](https://www.cefe.cnrs.fr/fr/ressources/logiciels/34-french/recherche/bc/bbp/264-logiciels), a standalone program. Recently Olivier Giminez developed a package called R2ucare. This allows to access these GOF tests from within the R software. While the package is still under active development, we can now for the first time carry out an entire CMR analysis in R. By doing so the reproducibility of the entire analysis is ensured. 

A mandatory first step for this type of analysis is a GOF test to verify if your initial model (the one which contains most parameters) fits the data well enough. You can download, install and load the R2ucare package from the GitHub [repo](https://github.com/oliviergimenez/R2ucare) which hosts the development version of the package. 

We will first prepare the data for the GOF test by removing primary capture occasion with less than 5 captures. To run the test we need to get the encounter histories and their counts or frequencies.
```{r message=FALSE, warning=FALSE, include=T}
#install.packages("R2ucare")
library(R2ucare)
prim_mat2=prim_mat[,-(which(apply(prim_mat,2,sum)<5))] # delete prim capt occasion with less than 5 captures
prim_hist=prim_mat2[-(which(apply(prim_mat2,1,sum)<1)),] # delete rows with no capture occasion 
prim.freq=rep(1,nrow(prim_hist))
prim.freq=apply(prim_hist,1,sum)
```

Once this is done one can perform 4 different tests, named test3SR, test3SM, test2CT and test2CL. We will not cover those in detail, but you must be aware that the Goodness of Fit test is the starting point for every analysis of capture-recapture data. We will focus on the ones that test for an excess of transients (test3SR), and the one which tests for trap dependence (test2CT). To learn more about these test we refer to the user manual of U-CARE and to Choquet et al., (2009). 

Following the procedure described in Figure 1, we first assess the overall fit of the CJS model using the function  overall_CJS.

```{r message=FALSE, warning=FALSE, include=T}
overall_CJS(prim_hist, prim.freq)
```

Clearly, the CJS model does not fit the data well. We then test for an effect of trap-dependence using the test2ct function. The null hypothesis on 'trap-dependence' states that when individuals are caught, they become aware of the trap and will actively seek or avoid it at the next trapping occasion (e.g. M. natalensis becomes trap-happy in our dataset). 

```{r message=FALSE, warning=FALSE, include=T}
test2ct = test2ct(prim_hist, prim.freq,verbose=F)
```
The signed square root (sign_test) of the Pearson chi-square statistic is negative, suggesting  there is an excess of individuals encountered at a given occasion among the individuals encountered at the previous occasion. 

We will then test for an excess of transients using the test3sr function. The null hypothesis on 'transience' states that there is no difference in the re-encounter probability of newly trapped and recaptured individuals.

```{r message=FALSE, warning=FALSE, include=T}
test3sr = test3sr(prim_hist, prim.freq,verbose=F)
```
 The p signed square root (sign_test) of the Pearson chi-square statistic is positive, suggesting  there is an excess of never seen again among the newly marked animals.

Major deviations against assumptions on 'transience' and 'trap-dependence' were found (see Results). Because we are interested in survival  and recapture rates of resident animals only (not in migration), we can remove all transient animals from the CMR data set (Sluydts et al 2007).  Transient animals were defined as individuals that were captured only once during one secondary capture occasion. These individuals were (most likely) not re-encountered because they moved outside of the trapping grid, and not because they died shortly after release. 

```{r message=FALSE, warning=FALSE, include=T}
prim_hist=prim_mat[-(which(apply(prim_sec_mat,1,sum)<2)),] # delete animals captures once in primary capture session
char=prim_sec_char[-(which(apply(prim_sec_mat,1,sum)<2)),]
```



#### **2.4.2 Fit the FULL model**

Because running full time-dependent CJS models is computationally heavy, we will show such an analysis on the first two years of the MOSA data only. First we need to put all the data in a correct data.frame which can be read by the MARK function.

```{r message=FALSE, warning=FALSE, include=T}
prim_hist=prim_hist[1:1498,1:30] # selection for trying
char=char[1:1498,]

sex=as.factor(char$SX) #make factor from host characteristic data
reprod=as.factor(char$REPR)
weight=as.numeric(char$WEIGHT)

library(stringi) # Make a string from the encounter data
i=1
ch=rep(NA,nrow(prim_hist))
for(i in 1:nrow(prim_hist)){
 ch[i]=stri_paste(prim_hist[i,],collapse='')
}
masto=data.frame(ch,sex,reprod,weight) # make data.frame
```


Second, we need to add details about the model we want to fit. This step is equivalent to reading in the data in program MARK. MAKE SURE TO DOWLOAD THE MARK PROGRAM at http://www.phidot.org/software/mark/downloads/


```{r message=FALSE, warning=FALSE, include=T}
library(RMark)
masto.proc<-process.data(masto, model="CJS", begin.time=1)
masto.ddl<-make.design.data(masto.proc)
```

Now we can easily obtain a summary for the capture histories in the m-array format.
```{r message=FALSE, warning=FALSE, include=T}
summary_ch(masto.proc, bygroup=T, marray=T)
```

In R we will specify the different models using code equivalent to fitting a regression type of model. In order to fit a model with a survival and recapture probability varying over time, eg the FULL model proposed here, we can specify as:

```{r message=FALSE, warning=FALSE, include=T}
Phi.time=list(formula=~time)
p.time=list(formula=~time)
```

Hence running such a model does not require manual adaptation of the PIM matrices or the design matrix in MARK, which can become a cumbersome process, especially for long-term CMR data. 

Fitting the model is done by using the wrapper function for *mark*. You can ask for a summary which will provide beta and real estimates at once.

```{r message=FALSE, warning=FALSE, include=T}
M1<-mark(masto.proc, masto.ddl, model.parameters=list(Phi=Phi.time, p=p.time), output=F, delete=T)
```

```{r message=FALSE, warning=FALSE, include=T}
summary(M1)$beta[1:5,]
```

### **2.4.3  Fit *a priori* models and model selection using AIC**

When performing a model selection procedure the idea is to have a set of *a priori* candidate models. In the case of the MOSA data we will start with the different time and constant combinations for survival and recapture to illustrate how to make a proper AIC table.

First we specify the remaining formulas for constant survival and recapture:
```{r message=FALSE, warning=FALSE, include=T}
Phi.ct=list(formula=~1)
p.ct=list(formula=~1)
```

Now we will fit the remaining model combinations 

```{r message=FALSE, warning=FALSE, include=T}
M2<-mark(masto.proc, masto.ddl, model.parameters=list(Phi=Phi.time, p=p.ct), output=F, delete=T)
M3<-mark(masto.proc, masto.ddl, model.parameters=list(Phi=Phi.ct, p=p.time), output=F, delete=T)
M4<-mark(masto.proc, masto.ddl, model.parameters=list(Phi=Phi.ct, p=p.ct), output=F, delete=T)
```

Summarizing the data in a table can be done by using the collect.models function, which will look in your current workspace for all the models of a certain type you specify. The function will return a summary table as in program MARK.

```{r message=FALSE, warning=FALSE, include=T}
collect.models(lx=c("M1", "M2", "M3", "M4"),type="CJS", table=T)
```

The first model clearly fitted best based on AIC criteria. We can plot the survival estimates from the first model using the following code.

```{r message=FALSE, warning=FALSE, include=T}
phitable <- get.real(M1, "Phi", se=T)
library(plotrix)
par(mfrow=c(1,1))
plotCI(time$time_interval2[1:30], phitable$estimate[1:30], ui=phitable$ucl[1:30], li=phitable$lcl[1:30],ylab="Apparent survival probability",xaxt='n',xlab="Years", main="MOSA - Mastomys natalensis",lwd=1,ylim=c(0,1))
year=seq(292,10512,by=365); years2=seq(1995,2023,by=1)
axis(1,year,years2)
```

Given the estimated recapture probability at each occassion and the number of individuals captured/sighted/detected, we can now estimate abundances for all but the first occasion.  This code was inspired by an [Rpubs page](https://rpubs.com/ogimenez/201392) of Olivier Giminez. 

We first estimate the recapture probabilities
```{r message=FALSE, warning=FALSE, include=T}
ptable <- get.real(M3, "p", se=T)
recaptures=ptable[c("estimate","se","lcl","ucl")][1:29,]
```

Based on the recapture probablities, it is possible to estimate the number of resident animals in the population at each primary capture with confidence intervals around the estimates. We also tried to estimate the percentage of transient animals, although this was a rather back-of-the-enveloppe calculation.

```{r message=FALSE, warning=FALSE, include=T}
obs = gregexpr("1", masto$ch) # Count number per primary trapping session - transients
n_obs = summary(as.factor(unlist(obs))) # Count number per primary trapping session - transients
transients=prim_cap_numb[1:30]-n_obs #count transient animals
transients_perc=transients/n_obs
estim_abundance_residents = n_obs[-1]/recaptures$estimate
estim_abundance_residents_lcl = n_obs[-1]/recaptures$lcl#lower confidence intervals
estim_abundance_residents_ucl = n_obs[-1]/recaptures$ucl#upper confidence intervals
estim_abundance_residents_transients=estim_abundance_residents+estim_abundance_residents*transients_perc[-1]

plotCI(time$time_interval2[2:30], estim_abundance_residents, ui=estim_abundance_residents_lcl, li=estim_abundance_residents_ucl,xaxt='n',xlab="Years",ylab="Number of animals",pch=19,col="black",ylim=c(0,1000),cex.lab=1.5,main="MOSA - Mastomys natalensis densities")
points(time$time_interval2[2:30],prim_cap_numb[2:30],col="red",pch=19)
points(time$time_interval2[2:30],estim_abundance_residents_transients,col="blue",pch=19)
year=seq(292,10512,by=365); years2=seq(1995,2023,by=1)
axis(1,year,years2)
legend("topright",legend=c("CJS-residents","CJS-transients+residents","Number captured"),pch=c(19,19), col=c("black","blue","red"), pt.cex=c(1, 1), cex=1)
```

We can see that the density estimations of the resident animals approximate the number of captured animals per primary trapping session in the field. 

### **2.4.4 Survival analysis**
We continue with the survival analysis on the whole dataset. However, because fitting the models to the full time dependent model takes to long , we reduce the time variable to "year" and "season". Season was divided in a non-breeding season between December and May, and a breeding season between June and November. The objective of this analysis was the verify if survival and recapture rates differed over the years and over season. We also included a sex effect to check if survival differs between males and females. To not overcomplicate the models, we did not include any interactions between de the fixed effects. We only show the results of the full model including season + year + sex. 

We again needed to remove all the transient animals in the full dataset and prepare the data for processing by the RMark package. 

```{r message=FALSE, warning=FALSE, include=T}
prim_hist=prim_mat[-(which(apply(prim_sec_mat,1,sum)<2)),] # delete animals captures once in primary capture session
char=prim_sec_char[-(which(apply(prim_sec_mat,1,sum)<2)),]

sex=as.factor(char$SX)
reprod=as.factor(char$REPR)
weight=as.numeric(char$WEIGHT)

library(stringi)
i=1
ch=rep(NA,nrow(prim_hist))
for(i in 1:nrow(prim_hist)){
  ch[i]=stri_paste(prim_hist[i,],collapse='')
}
masto=data.frame(ch,sex,reprod,weight)

```

Second, we need to add details about the model we want to fit in the design matrix.
```{r message=FALSE, warning=FALSE, include=T}
# library(RMark)
time3=(time$time_interval+2)[-1]#find time difference between trapping session, add two days to interval, remove first day
diffsc=round(time3/30,1)# Scale all differences to 30 days.
masto.proc<-process.data(masto, model="CJS", groups=("sex"),time.intervals=diffsc)#correct for differences in time interval
masto.ddl<-make.design.data(masto.proc)
```

Here, we simplify the full time dependent model (estimating parameters for each primary trapping session) to a seasonal and year model. 

```{r message=FALSE, warning=FALSE, include=T}
time2=seq(1,nrow(time)) # number of primary trapping sessions
year=factor(c(substr(time$date, 7, 10)[-1],"2022"))#shift one month to let the breeding season end in November
year[52]="1998";year[65]="1999";year[209]="2009";year[267]="2015";year[348]="2022"; #Asure all is correct
month=substr(time$date, 4, 5) #select month
season=rep("breeding",length(month)) 
season[which((month=="12") | (month=="01") | (month=="02") | (month=="03") | (month=="04") | (month=="05"))]="non-breeding"
season=factor(season) #select which months are breeding and non-breeding seasons
env.data=data.frame(time=c(1,as.numeric(as.character(unique(masto.ddl$p$time)))),ses=season,yr=year)#put in one matrix
```

Merge covariate matrix to the design matrix voor P (recapture), Phi (Survival) parameters.
```{r message=FALSE, warning=FALSE, include=T}
masto.ddl$p=merge_design.covariates(masto.ddl$p,env.data)
masto.ddl$Phi=merge_design.covariates(masto.ddl$Phi,env.data)
```

Continue to build the different models you want to fit. Running each model takes a considerable amount of computation time. Therefore, we only activate the first full model.

```{r message=FALSE, warning=FALSE, include=T}
Phi.ses.yr.sex=list(formula=~ses+yr+sex)
p.ses.yr.sex=list(formula=~ses+yr+sex)
# Phi.ses.yr=list(formula=~ses+yr)
# p.ses.yr=list(formula=~ses+yr)
# Phi.ses=list(formula=~ses)
# p.ses=list(formula=~ses)
# Phi.sex=list(formula=~sex)
# p.sex=list(formula=~sex)
# Phi.yr=list(formula=~yr)
# p.yr=list(formula=~yr)
# Phi.ct=list(formula=~1)
# p.ct=list(formula=~1)
M1<-mark(masto.proc, masto.ddl, model.parameters=list(Phi=Phi.ses.yr.sex, p=p.ses.yr.sex), output=F, delete=T)
# M2<-mark(masto.proc, masto.ddl, model.parameters=list(Phi=Phi.ses.yr, p=p.ses.yr), output=F, delete=T)
# M3<-mark(masto.proc, masto.ddl, model.parameters=list(Phi=Phi.ses, p=p.ses), output=F, delete=T)
# M4<-mark(masto.proc, masto.ddl, model.parameters=list(Phi=Phi.sex, p=p.sex), output=F, delete=T)
# M5<-mark(masto.proc, masto.ddl, model.parameters=list(Phi=Phi.yr, p=p.yr), output=F, delete=T)
# M6<-mark(masto.proc, masto.ddl, model.parameters=list(Phi=Phi.ct, p=p.ct), output=F, delete=T)
# collect.models(lx=c("M1", "M2", "M3", "M4","M5,M6"),type="CJS", table=T)
```

Model 1 has the lowest AICc, but the difference with the second best model is only 1.5 (DeltaAICc). We therefore conclude that a model with 'sex' does not improve the model fit. The difference in AICc is much larger when removing 'year' (DeltaAICc=456) or 'season' (DeltaAICc=2077) from the models. While a survival analysis with model 2 (season + year) would make most sense to continue, we here continue with model 1 to visualize the different effects of season, year and sex on the survival and recapture probabilities. 

Plot survival estimates in function of year and season separately for males and females. 

```{r message=FALSE, warning=FALSE, include=T}
library(plotrix)
par(mfrow=c(1,2),mar=c(4,4,1,0.2))
phitable <- get.real(M1, "Phi", se=T) # get survival probabilities from model
surv=phitable[which(phitable$sex=="F"),][1:350,] # select first 350 rows and females only
surv_estimate=unique(surv$estimate) # estimates
surv_lcl=unique(surv$lcl)# lower confidence intervals
surv_ucl=unique(surv$ucl)# upper confidence intervals

plotCI(seq(1,length(surv_estimate)), surv_estimate, ui=surv_ucl, li=surv_lcl,xaxt='n',xlab="Years",ylab="Survival estimates",pch=19,col=rep(c(1,2),26),ylim=c(0,1),cex.lab=1.5,main="Mastomys (Females)")
year=seq(2,length(surv_estimate)+1,by=4); years2=seq(1995,2023,by=2)
axis(1,year,years2)
legend("bottomleft",legend=c("Non-breeding","Breeding"),pch=c(19,19), col=c("black","red"), pt.cex=c(1, 1), cex=0.8)

phitable <- get.real(M1, "Phi", se=T)
surv=phitable[which(phitable$sex=="M"),][1:350,]
surv_estimate=unique(surv$estimate)
surv_lcl=unique(surv$lcl)
surv_ucl=unique(surv$ucl)

plotCI(seq(1,length(surv_estimate)), surv_estimate, ui=surv_ucl, li=surv_lcl,xaxt='n',xlab="Years",ylab="Survival estimates",pch=19,col=rep(c(1,2),26),ylim=c(0,1),cex.lab=1.5,main="Mastomys (Males)")
year=seq(2,length(surv_estimate)+1,by=4); years2=seq(1995,2023,by=2)
axis(1,year,years2)
```

Do the same for the recapture estimates. 
```{r message=FALSE, warning=FALSE, include=T}
par(mfrow=c(1,2),mar=c(4,4,1,0.2))
ptable <- get.real(M1, "p", se=T)
surv=ptable[which(ptable$sex=="F"),][1:350,]
surv_estimate=unique(surv$estimate)
surv_lcl=unique(surv$lcl)
surv_ucl=unique(surv$ucl)

plotCI(seq(1,length(surv_estimate)), surv_estimate, ui=surv_ucl, li=surv_lcl,xaxt='n',xlab="Years",ylab="Recapture estimates",pch=19,col=rep(c(1,2),26),ylim=c(0,1),cex.lab=1.5,main="Mastomys natalensis (Females)")
year=seq(2,length(surv_estimate)+1,by=4); years2=seq(1995,2023,by=2)
axis(1,year,years2)
legend("bottomleft",legend=c("Non-breeding","Breeding"),pch=c(19,19), col=c("black","red"), pt.cex=c(1, 1), cex=0.8)

ptable <- get.real(M1, "p", se=T)
surv=ptable[which(ptable$sex=="M"),][1:350,]
surv_estimate=unique(surv$estimate)
surv_lcl=unique(surv$lcl)
surv_ucl=unique(surv$ucl)

plotCI(seq(1,length(surv_estimate)), surv_estimate, ui=surv_ucl, li=surv_lcl,xaxt='n',xlab="Years",ylab="Recapture estimates",pch=19,col=rep(c(1,2),26),ylim=c(0,1),cex.lab=1.5,main="Mastomys natalensis (Males)")
year=seq(2,length(surv_estimate)+1,by=4); years2=seq(1995,2023,by=2)
axis(1,year,years2)
```

**Conclusion:** Survival estimates are clearly higher in the breeding compared to the non-breeding season. This difference remains more or less constant over the years and does not differ between the sexes. Recapture estimates are also higher during the non-breeding season compared to the breeding season. Here, we do see a patter over the years. The recapture probability seems to decrease between 1994-2003, increase between 2004-2011 after which it starts to decrease again. No significant differences are detected between males and females. 


