---
title: "Mosa - Processing data"
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)
```
## **Processing MOSA data**

### **1.1 Check new data on errors**

This R-script reads the MOSA-CMR file as a txt  and checks if any obvious errors are present. Data from new captures should be added to the main excel and processed. The R-code will show  if any changes occured for recaptured animals (the same TOEcode) for species ("species_change"), NA values for X coordinates ("X_change"), NA values for X coordinates ("Y_change"), weight drops > half in comparison to the previous capture ("weight_changes"), sex conversion ("sex_changes"), or sexual reproduction changes from reproductively active to inactive ("sexcon_change") and, if animals with the same TOE code are captured longer than 2 years ago (not possible as Mastomys lives < 1 year in the wild). All changes were checked in the final dataset, except for errors against sexual condition reversion (active to inactive). This allows the researcher to decide for himself if he considers an animal that was previously sexually active but later inactive, to be active or not.


Read in all data from main txt MOSA file that includes all the captures. 
```{r,,warning=F,echo=T,message=F}
rm(list=ls())
# install.packages("readr")
library("readr")

mosa<-data.frame(read.csv(url("https://raw.githubusercontent.com/Joachimmarien/MOSA/main/mosa_Feb2023_flags.csv"), sep=';'))
# mosa <- read.delim("C:/Users/joach/OneDrive - Universiteit Antwerpen/Bureaublad/PhD selection/MOSA/CMR cleaning/since Feb 2022/mosa_Feb2023_flags2.txt")
mosa=mosa[,1:10] #select relevant column
mosa2=mosa[which(mosa[,3]=="MOSA"),]#Select field
mosa2$DATE=as.Date(as.character(mosa2$DATE), format="%d/%m/%Y") #Put date in correct form
ID=unique(mosa2$TOE) #select unique individuals
mosa2$WEIGHT=as.numeric(mosa2$WEIGHT)
```

Build data frame 'df3' that can be filled with all errors. 
```{r,,warning=F,echo=T,message=F}
i=1
df3=data.frame()
for (i in 1:length(ID)){
  df=mosa2[which(mosa2$TOE==ID[i]),]
  if (any(is.na(df))){
    df2b=cbind(df,matrix(NA,nrow=nrow(df),ncol=7))
    colnames(df2b)=colnames(df3)
    df3=rbind(df3,df2b)
  } else{
    if(nrow(df)>0){
      species="n";x="n";y="n";weight="n";sx="n";sexcon1="n"; date="n"
      if(nrow(df)>1){
        if(any(is.na(df$SPECIES)==T)){species="SP_NA"} #Check NA for species
        if(length(unique(df$SPECIES))>1){species="SP_CH"} 
        if(any(is.na(df$X)==T)){x="X_NA"}#Check NA for X coordinate
        if(any(is.na(df$Y)==T)){y="Y_NA"}#Check NA for Y coordinate
        if(any(is.na(df$WEIGHT)==T)){weight="WE_NA"}#Check NA for body weight
        if(any(is.na(df$SX)==T)){sx="SX_NA"}#Check NA for sex
        if(length(unique(df$SX))>1){sx="SX_CH"}#Check NA for sex
        if(any(is.na(df$SEXCON)==T)){sexcon1="CON_NA"}#Check NA for sexual condition
        if(any(is.na(df$DATE)==T)){date="dat_NA"}#Check NA for species
        j=2
        vec_weight=c()
        vec_weight2=c()
        sexconF=rep(0,nrow(df))
        sexconM=rep(0,nrow(df))
        sub=substring(df$SEXCON, 1, 1)
        dates=rep(0,nrow(df))
        for(j in 2:nrow(df)){
          vec_weight[j-1]=df$WEIGHT[j-1]/(df$WEIGHT[j]+1)
          if((sub[j-1]=="P") & (sub[j]!="P")){sexconF[j]=1} #check for sexual activity changes
          if((sub[j-1]=="S") & (sub[j]!="S")){sexconM[j]=1} #Check sexual activity change
          sexcon=c(sexconF,sexconM)
          dates[j-1]=as.numeric(df$DATE[j]-df$DATE[j-1]) # check if animals live linger than 2 year
        }
        if(any((vec_weight>2)|vec_weight==0)){weight="WE_CH"} # check body weight decrease for more than half of body weight
        if(any(sexcon>0)){sexcon1="CON_CH"} #check for sexual activity changes
        if(any(dates>400)){date="dat_long"} # check if animals live linger than 2 year
      } 
      df2=cbind(df,paste(species,x,y,weight,sx,sexcon1,date))
      df2b=cbind(df,species,x,y,weight,sx,sexcon1,date)
      if(any(df2[,11]!="n n n n n n n")){df3=rbind(df3,df2b)}
    }
  }
}
```

Show number of individual animals with an error. Most of these errors are changes from re productively active to inactive, which is not changed in the final data set as explained above. 
```{r,,warning=F,echo=T,message=F}
length(unique(df3[,8]))#number of wrong species
```

Show  individuals that change from species (e.g. first capture is *Mastomys natalensis* and second is *Lemniscomys rosalia*). 
```{r,,warning=F,echo=T,message=F}
species_change=df3[which(df3$species=="SP_CH" |df3$species=="SP_NA" ),]; length(unique(species_change[,8]))
species_change
```
Show individuals that have no X or Y-coordinate 
```{r,,warning=F,echo=T,message=F}
x_change=df3[which(df3$x=="X_NA"),]; length(unique(x_change[,8]))
x_change
y_change=df3[which(df3$y=="Y_NA"),]; length(unique(y_change[,8]))
y_change
```

Show individuals for which the body weight drops with more than half in comparison to the previous capture
```{r,,warning=F,echo=T,message=F}
weight_change=df3[which(df3$weight=="WE_CH" |df3$weight=="WE_NA" ),]; length(unique(weight_change[,8]))
weight_change
```

Show individuals for which the sex changes over captures
```{r,,warning=F,echo=T,message=F}
sex_change=df3[which(df3$sx=="SX_CH" |df3$sx=="SX_NA" ),]; length(unique(sex_change[,8]))
sex_change
```

Show individuals for which the sex_condition  changes over captures (from active to inactive)
```{r,,warning=F,echo=T,message=F}
sexcon_changes=df3[which(df3$sexcon1=="CON_CH"|df3$sexcon1=="CON_NA"),]; length(unique(sexcon_changes[,8]))
sexcon_changes[1:5,1:10]
```

Show individuals for which a recapture was more than two years ago. Some individuals might just meet these criteria, for others these are obvious mistakes in toe-code readings. 
```{r,,warning=F,echo=T,message=F}
date_changes=df3[which(df3$date=="dat_long"|df3$date=="dat_NA"),]; length(unique(date_changes[,8]))
date_changes[,1:10]
```



### **1.2 Modify the mosa data**  

The first lines convert the Y-axis letters to numbers. The FOR-loop is  used to convert the different reproductive states for male and female rodents to "Active" or "Inactive".There might be more elegant ways to code this, but this also works. 

```{r,,warning=F,echo=T,message=F}
Y2=rep(NA,nrow(mosa2))#change Y-axis letters to numbers
REPR=rep("AC",nrow(mosa2))
i=1
for(i in 1:nrow(mosa2)){
if(is.na(mosa2[i,6])){
  Y2[i]=NA
}else {
  if(mosa2[i,6]=="A"){Y2[i]=1}
  if(mosa2[i,6]=="B"){Y2[i]=2}
  if(mosa2[i,6]=="C"){Y2[i]=3}
  if(mosa2[i,6]=="D"){Y2[i]=4}
  if(mosa2[i,6]=="E"){Y2[i]=5}
  if(mosa2[i,6]=="F"){Y2[i]=6}
  if(mosa2[i,6]=="G"){Y2[i]=7}
  if(mosa2[i,6]=="H"){Y2[i]=8}
  if(mosa2[i,6]=="I"){Y2[i]=9}
  if(mosa2[i,6]=="J"){Y2[i]=10}
}
  if(is.na(mosa2[i,10])){
    REPR[i]=NA
  } else {
  if((mosa2[i,10]=="AN") | (mosa2[i,10]=="AV") | (mosa2[i,10]=="CSN")) {REPR[i]="IN"}
  }
}

mosa3=cbind(mosa2,REPR,Y2)# final data set
mosa3[1:5,]
```

### **1.3 Convert Mosa data to an encounter matrix**

This R-script reads in the corrected MOSA data and produces different output  files in which the encounter histories (100101) for each animal are given in a typical  RxK matrix, together with individual characteristics of the animal at first capture  (sex, reproductive state, body weight, Toe-ID, X-axis and Y-axis. The different files produced are:

  - primary_secondary_capture_3consecutive: presents the robust design matrix for the primary and secondary trappings sessions performed at the MOSA field that consisted for 3 consecutive days. It shows the encounter history for each animal at each trapping occasion.
  
  - primary_secondary_time_variables_3consecutive: shows the real date of the trapping occasion, primary and secondary capture occasion number, time interval between secondary capture occasions and time interval from starting date until the trapping occasion. Only for trappings sessions performed at the MOSA field that consisted for 3 consecutive days.


We will first preprocess the data. Each row in the mosa matrix (corresponding to an individual capture) is read and given a corresponding primary trapping session number. The output of the function is list with: a vector of dates for the different primary trapping sessions, the number of the primary trapping session, the day of the secondary trapping session in the primary trapping session. 

```{r,,warning=F,echo=T,message=F}
MN=mosa3
date=unique(MN$DATE)
ID=unique(MN$TOE)

prim_cap=function(MN,date,ID){
i=2
new_dates=as.Date(rep(NA,nrow(MN)))
trap_sess=rep(0,nrow(MN))
j=0;k=0

diff_days=MN$DATE-date[i-1]
new_dates[which((diff_days==0))]=date[i-1]
j=j+1; k=k+1
trap_sess[which((diff_days==0))]=j
cons_day=rep(NA,length(date))
cons_day[i-1]=k

for(i in 2:length(date)){
  time_interval=as.numeric(date[i]-date[i-1])
  if(time_interval==1){
    diff_days=MN$DATE-date[i]
    new_dates[which((diff_days==0))]=rep(new_dates[which((diff_days==0))[1]-1],length(which((diff_days==0))))
    trap_sess[which((diff_days==0))]=trap_sess[which((diff_days==0))[1]-1]
    k=k+1
    cons_day[i]=k
    
  } else {
    diff_days=MN$DATE-date[i]
    new_dates[which((diff_days==0))]=date[i]
    j=j+1
    trap_sess[which((diff_days==0))]=j
    k=1
    cons_day[i]=k
  }
}
data3=list(new_dates,trap_sess,cons_day)
}
data3=prim_cap(MN,date,ID)
```

Trapping occasions that were not included in the standard  28 day robust design were removed from the dataset to not interfere with the further survival and density estimates. 

```{r,,warning=F,echo=T,message=F}
xx=data.frame(data3[[3]],date)
MN=cbind(MN,data3[[1]],data3[[2]])
MN=MN[-(which(MN[,2]==xx[7,2])),]
MN=MN[-(which(MN[,2]==xx[8,2])),]
MN=MN[-(which(MN[,2]==xx[9,2])),]
MN=MN[-(which(MN[,2]==xx[10,2])),]
MN=MN[-(which(MN[,2]==xx[32,2])),]
MN=MN[-(which(MN[,2]==xx[33,2])),]
MN=MN[-(which(MN[,2]==xx[34,2])),]
MN=MN[-(which(MN[,2]==xx[35,2])),]
MN=MN[-(which(MN[,2]==xx[39,2])),]
MN=MN[-(which(MN[,2]==xx[40,2])),]
MN=MN[-(which(MN[,2]==xx[41,2])),]
MN=MN[-(which(MN[,2]==xx[57,2])),]
MN=MN[-(which(MN[,2]==xx[58,2])),]
MN=MN[-(which(MN[,2]==xx[59,2])),]
MN=MN[-(which(MN[,2]==xx[60,2])),]
MN=MN[-(which(MN[,2]==xx[61,2])),]
MN=MN[-(which(MN[,2]==xx[62,2])),]
MN=MN[-(which(MN[,2]==xx[63,2])),]
MN=MN[-(which(MN[,2]==xx[67,2])),]
MN=MN[-(which(MN[,2]==xx[89,2])),]
MN=MN[-(which(MN[,2]==xx[90,2])),]
MN=MN[-(which(MN[,2]==xx[91,2])),]
MN=MN[-(which(MN[,2]==xx[92,2])),]
MN=MN[-(which(MN[,2]==xx[96,2])),]
MN=MN[-(which(MN[,2]==xx[97,2])),]
MN=MN[-(which(MN[,2]==xx[113,2])),]
MN=MN[-(which(MN[,2]==xx[117,2])),]
MN=MN[-(which(MN[,2]==xx[118,2])),]
MN=MN[-(which(MN[,2]==xx[122,2])),]
MN=MN[-(which(MN[,2]==xx[123,2])),]
MN=MN[-(which(MN[,2]==xx[235,2])),]
MN=MN[-(which(MN[,2]==xx[236,2])),]
MN=MN[-(which(MN[,2]==xx[240,2])),]
MN=MN[-(which(MN[,2]==xx[241,2])),]
MN=MN[-(which(MN[,2]==xx[245,2])),]#
MN=MN[-(which(MN[,2]==xx[246,2])),]
MN=MN[-(which(MN[,2]==xx[247,2])),]
MN=MN[-(which(MN[,2]==xx[248,2])),]
MN=MN[-(which(MN[,2]==xx[252,2])),]
MN=MN[-(which(MN[,2]==xx[253,2])),]
MN=MN[-(which(MN[,2]==xx[254,2])),]
MN=MN[-(which(MN[,2]==xx[291,2])),]
MN=MN[-(which(MN[,2]==xx[301,2])),]
MN=MN[-(which(MN[,2]==xx[371,2])),]
MN=MN[-(which(MN[,2]==xx[372,2])),]
MN=MN[-(which(MN[,2]==xx[373,2])),]
MN=MN[-(which(MN[,2]==xx[374,2])),]
MN=MN[-(which(MN[,2]==xx[375,2])),]
MN=MN[-(which(MN[,2]==xx[376,2])),]
MN=MN[-(which(MN[,2]==xx[377,2])),]
MN=MN[-(which(MN[,2]==xx[381,2])),]
MN=MN[-(which(MN[,2]==xx[382,2])),]
MN=MN[-(which(MN[,2]==xx[383,2])),]
MN=MN[-(which(MN[,2]==xx[384,2])),]
MN[(which(MN[,2]==xx[469,2])),2]=as.Date("2004-04-13")
MN=MN[-(which(MN[,2]==xx[506,2])),]

```

Make the encounter matrix for the secondary trapping session by finding all captures for each animal. Add the characteristics (body weight, sex, reproductive state, x-y capture) of the animal when it was first captured. 
```{r,,warning=F,echo=T,message=F}
date=unique(MN$DATE) #create vector with unique dates
ID=unique(MN$TOE)#create vector with unique individuals
df=matrix(0,nrow=length(ID),ncol=length(date)) #empty matrix for encounter data
character=data.frame() #empty data.frame for host characteristic data

i=1
for(i in 1:length(ID)){ #for loop over the different unique individuals
  ID2=MN[which(MN$TOE==ID[i]),]
  df[i,pmatch(ID2$DATE,date)]=1 #when is the animal captured
  character=rbind(character,ID2[1,4:12]) # what are its characteristics
  print(i)
  }
```

Delete rows when no animals were captured and keep *Mastomys natalensis* only. Line 293 obviously has to change if someone wants to check the encouner histories of *Tatera*,  *Lemniscomys* or *Mus*. 
```{r,,warning=F,echo=T,message=F}
df=df[-(which(character$TOE=="NO_CAPTURE")),] #remove rows for which no animals were captured
character=character[-(which(character$TOE=="NO_CAPTURE")),] #remove rows for which no animals were captured
df=df[-(which(character$SPECIES!="MN")),] #keep mastomys natalensis only
character=character[-(which(character$SPECIES!="MN")),]#keep mastomys natalensis only
```

Make data frame and write to txt file in your activated directory.
```{r,,warning=F,echo=T,message=F}
mat=data.frame(character,df)
setwd("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")
write.table(mat, file = "C_primary_secondary_capture_3consecutive_o.txt", sep="\t")
```

**1.4 Create time variable matrix**

This FOR loop gives the different time variables needed for the later density and survival analyses. It produces for each secondary trapping sessions: the date, the time since previous secondary trapping session, time since start of the study, number of primary capture session. 

```{r,,warning=F,echo=T,message=F}
i=2;j=1;k=1
time_interval=rep(NA,length(date)) #time since previous secondary trapping session
time_interval2=rep(NA,length(date))#time since start of the study
prim_capt=rep(NA,length(date)) # number of primary capture session
sec_capt=rep(NA,length(date)) # number of secondary capture session within primary session
time_interval[1]=1
for (i in 2:length(date)){  #start FOR loop 
  time_interval[i]=as.numeric(date[i]-date[i-1]) #time since previous secondary trapping session
  time_interval2[i]=as.numeric(date[i]-date[1])#time since start of the study
  if(time_interval[i]==1){ #if days are consecutive
    prim_capt[i]=j # number of primary capture session
    k=k+1
    sec_capt[i]=k# number of secondary capture session within primary session
  } else { #if days are not consecutive
    j=j+1
    prim_capt[i]=j# number of primary capture session
    k=1
    sec_capt[i]=k# number of secondary capture session within primary session
  }
} 
prim_capt[1]=1; sec_capt[1]=1; time_interval2[1]=0 #Indicate first capture
all_capt=seq(1,length(date),1)
```


Make data frame for time variables and write to txt file in your activated directory.

```{r,,warning=F,echo=T,message=F}
mat2=data.frame(date,prim_capt,sec_capt,all_capt,time_interval,time_interval2)
setwd("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")
write.table(mat2, file = "C_primary_secondary_time_variables_3consecutive_m.txt", sep="\t")
```
