F2 Markdown script
Bram Al-Sharifi, Sarina Kuipers, Juliette Debarge, Beau Strijker
Introduction to the data
We start with 5 different data sets. The first data set (animal_tracking) contains GPS information about the animals. Intrusion_tracking includes GPS information from the 5 different experiments. 1. Mimicking the movement of tourists by car (land-rover). 2. Mimicking the movement of tourists on foot . 3. Like exp. 1 but they get out of the car and split off into a small group that stands behind a rock or tree for 20 mins and then leaves. 4. Like exp 2 but they split off on foot into a small group, stand behind a rock or tree for 20 mins and then leaves. 5. Group of anti-poaching guards the area with the aim of tracking and killing the animal with the red spot on its back. There is no animal with the dot, but the poachers don’t know this, so this mimics real poaching behavior the best.
The intrusion_overview and animal_metadata hold information relating to experiments and the individual animals respectively. Finally, the weather_data data set holds information about the weather in that area, like wind, temperature and solar conditions.
Tidying Data
First we have to download the data and run the libraries
#downloading the packages for tidying and modeling
library(tidyverse)
library(momentuHMM)
library(tsibble)
library(corrplot)
library(ggfortify)
library(Rtsne)
library(raster)
#downloading the databases
animal_tracking <- read_rds("data/raw/animal_tracking.rds")
intrusion_tracking <- read_rds("data/raw/intrusion_tracking.rds")
intrusion_overview <- read_rds("data/raw/intrusion_overview.rds")
animal_metadata <- read_csv("data/raw/animal_metadata.csv")
weather <- read_csv("data/raw/weather_data.rds")
#Set the seed
set.seed(2024)First we will want to merge animal tracking, animal metadata, intrusion tracking and intrusion overview, as well as filter the 4 different animals so we only are working with zebra data.
At the same time, we must fix the coordinate system in animal
tracking so that it is the same as the rest of the data set. We do this
by multiplying the x and y value by 0.1, as it was originally in
decimeters. Then add the following values to the UTM35S
coordinates.
+ x –> 587800 E + y –> 7318800 N
The value in the temp column is also 100 times larger, which we mutate that so the values are accurate in degrees Celsius.
animal_tracking %>%
mutate(x = ((x*0.1)+587800),
y = ((y*0.1)+7318800)) %>% #adjusting the distances for animal tracking
left_join(animal_metadata, by="i")%>%
filter(species == "zebra") %>% #filtering for zebra data only
mutate(day = str_sub(t, start = 1, end = 10), # creating a column for day using t
temp = temp/100) %>% #making temp the real values
left_join((intrusion_tracking %>%
mutate(dayexp = str_sub(t,start = 1, end = 10))), by = "t") %>%
left_join(intrusion_overview, by = c("expid", "exptype")) %>%
rename(ID = i,
dttm = t)%>%
mutate(ID == as.factor(ID))-> zebra The next step is to create additional columns that computes the distances between the animal and the poachers/experiments, using the Pythagorean theorem. There can be 2 experiments running at the same time, so we are calculating the distances for columns mode 1 and mode 2.
zebradist <- zebra %>%
mutate(dx1 = x - x_1,
dy1 = y - y_1,
dist1 = sqrt(dx1^2 + dy1^2),
dx2 = x - x_2,
dy2 = y - y_2,
dist2 = sqrt(dx2^2 + dy2^2)) %>%
dplyr::select(-c(dx1, dy1, dx2, dy2))| ID | dist1 | dist2 |
|---|---|---|
| 17 | 3986.463 | 3988.380 |
| 17 | 4038.696 | 4043.050 |
| 17 | 4035.221 | 4022.997 |
Next we created a data set to specify the time of day, and filter the data set from 5am to 2pm.
This results in 2 different data sets. 1. Data for the morning on days that have experiments (even if those experiments aren’t in the morning) 2. Data for the morning on days without any experiments at all that day
#filter out times between 0500 and 1400
zebradist %>%
mutate(morning_start = 5,
morning_end = 14,
hour = hour(dttm),
timeoftheday = case_when(hour >= morning_start &
hour <= morning_end ~ "morning")) %>%
drop_na(timeoftheday) %>%
#this removes the data that is not related to the morning
mutate(burst = as.factor(day)) -> morningdata
#create lists of data for the morning
listdayexp <- unique(morningdata$dayexp)
listall <- unique(morningdata$day)
listdaywithoutexp <- listall[!listall %in% listdayexp]
#Data with experiments on that day
morningdata %>%
filter(day %in% listdayexp) %>%
mutate(burst = as.factor(day)) -> zebraexpmorn
#Data without experiments that day
morningdata %>%
filter(day %in% listdaywithoutexp) %>%
mutate(burst = as.factor(day)) -> zebrawithoutexpThe final step for tidying is to add a column of data that states whether the data is disturbed or not.
- Disturbed is data where the experiment is within 500m of the animals or 30 mins after they are no longer within 500 mins of the zebra
- Undisturbed is the rest of the data where the experiment is further than 500m of the animals
#creating a column that states when dist 1 or 2 are smaller than 500, then label500 is 1, otherwise 0
zebmornlab<- morningdata %>%
mutate(label500 = case_when(dist1 <= 500 ~ 1,
dist2 <= 500 ~ 1,
.default = 0))
#creating a column for label that states if label500 = 1, then the next 30 minutes is also 1
zebmornlab %>%
group_by(ID, burst) %>%
arrange(ID, burst, dttm) %>%
mutate(disturbedlabel = case_when(
lag(label500, 1) == 1~"Disturbed",
lag(label500, 2) == 1~"Disturbed",
lag(label500, 3) == 1~"Disturbed",
lag(label500, 4) == 1~"Disturbed",
lag(label500, 5) == 1~"Disturbed",
lag(label500, 6) == 1~"Disturbed",
lag(label500, 7) == 1~"Disturbed",
lag(label500, 8) == 1~"Disturbed",
lag(label500, 9) == 1~"Disturbed",
lag(label500, 10) == 1~"Disturbed",
lag(label500, 11) == 1~"Disturbed",
lag(label500, 12) == 1~"Disturbed",
lag(label500, 13) == 1~"Disturbed",
lag(label500, 14) == 1~"Disturbed",
lag(label500, 15) == 1~"Disturbed",
lag(label500, 16) == 1~"Disturbed",
lag(label500, 17) == 1~"Disturbed",
lag(label500, 18) == 1~"Disturbed",
lag(label500, 19) == 1~"Disturbed",
lag(label500, 20) == 1~"Disturbed",
lag(label500, 21) == 1~"Disturbed",
lag(label500, 22) == 1~"Disturbed",
lag(label500, 23) == 1~"Disturbed",
lag(label500, 24) == 1~"Disturbed",
lag(label500, 25) == 1~"Disturbed",
lag(label500, 26) == 1~"Disturbed",
lag(label500, 27) == 1~"Disturbed",
lag(label500, 28) == 1~"Disturbed",
lag(label500, 29) == 1~"Disturbed",
lag(label500, 30) == 1~"Disturbed",
label500 == 1 ~ "Disturbed",
label500 == 0 ~ "Undisturbed")) %>%
dplyr::select(-label500)->zebmornlabPCA, t-SNE, and further Exploration
The environmental variables must be tested, to see what effect they could have on the movement of disturbed and undisturbed zebras.
First we will have to prepare the data for the modelling by dplyr::selecting the useful columns for analysis
zebmornlab %>%
dplyr::select(ID,
burst,
dttm,
x,
y,
temp,
day,
dayexp,
start,
end,
disturbedlabel,
dist1,
dist2,
acs.mean.mean,
acs.mean.peak,
acs.mean.sd,
acs.sd.mean,
acs.sd.peak,
acs.fracmove) -> datThen to cut down on computation time, we selected 6 days, with 6 different individuals
daysSelect <- c("2017-11-07","2017-11-09","2017-11-10", "2017-11-13", "2017-11-22", "2017-12-04")
IDSelect <- c("17", "22", "31", "90", "91", "92")
datTS <- dat %>%
filter(day %in% daysSelect) %>% #filters to the days we want
filter(ID %in% IDSelect) %>% # filters to the individiuals we want to keep
dplyr::select(-day, -dayexp, -start, -end) %>% #we no longer need these columns
as_tsibble(key = c(ID, burst),
index = dttm)
datTS %>%
unite("ID", c(ID,burst)) %>%
as.data.frame(dat) -> datDF
datDF[5] <- as.data.frame(scale(datDF[5])) #standardize the integer columnsPCA
First we run a PCA with the following script, by first separating out the accelerometer data
datx <- prepData(data=datDF, covNames="temp")%>%
dplyr::select(step,
acs.mean.mean,
acs.mean.peak,
acs.mean.sd,
acs.sd.mean,
acs.sd.peak,
acs.fracmove,
disturbedlabel) %>% #select the columns we want to keep
drop_na()
#pull the labels
labelsPCA <- datx %>% dplyr::pull(disturbedlabel)
datx[1:7] <- as.data.frame(scale(datx[1:7]))
datx %>%
dplyr::select(-disturbedlabel) %>%
drop_na()-> datxPCAThe next step is to look at the corrplots to see if the ACS values and step are correlated
cor = cor(datxPCA)
cor_pc <- cor(datxPCA)
round(cor_pc, digits=2)
Lets have a closer look by running a PCA and looking at the output:
Here we can see that there is a strong correlation between acs.sd.peak
and sd.mean. And the rest have a strong affect on PC1 axis. The
distribution of disturbed is denser around the acs.mean.peak till step,
which sho that these variables have a strong correlation to increasing
of these values.
TSNE
Just for explorations sake, lets have a look at a non-linear model, the t-SNE
This model shows that there is a strong clustering when the values are
non-linear, that the Disturbed Zebras do have a value that could
determine poacher specific behaviors.
Animated model
Another exploration tool we implemented was running an animated diagram of a single day that had a poacher to visualize what it would look like.
Modeling
First we will have to prepare the data for an Hidden Markov Model (HMM).
datDF <- datDF %>%
dplyr::select(ID,
dttm,
x,
y,
temp,
disturbedlabel)
preppedData <- prepData(data=datDF, covNames="temp")
#Lets retrieve the labels
datxlabels <- preppedData %>%
pull(disturbedlabel)Model 0
The 0th model considers our original estimates that there would be 3 states of movement. However, we considered this to be too little as the variation in movement was too great.
# get the mean and sd values from the step column
hist(preppedData$step)
mean(preppedData$step, na.rm = TRUE) # 3.65
sd(preppedData$step, na.rm = TRUE) # 6.10
#plot(preppedData)
#Look at the distribution of angle
hist(preppedData$angle)
#Distance to use for the parameters
dist = list(step = "gamma",
angle = "wrpcauchy")
#starting parameters
Par0_m0 <- list(step=c(40, 4, 2, # means of 3states
8, 2, 2, # sd of first state and second state
0, 0, 0.3), # zero mass parameter specification for first and second state
angle=c(0.7, 0.3, 0.1)) # concentration high for first state, low for second.
m0 <- fitHMM(data = preppedData,
nbStates = 3,
dist = dist,
Par0 = Par0_m0,
estAngleMean = list(angle=FALSE)) # angle has a mean of zero## Value of the maximum log-likelihood: -22586.06
##
##
## step parameters:
## ----------------
## state 1 state 2 state 3
## mean 7.905359e+00 1.456498e+00 0.08171422
## sd 5.807768e+00 1.122567e+00 208.47152891
## zeromass 1.000132e-08 1.000015e-08 0.99999995
##
## angle parameters:
## -----------------
## state 1 state 2 state 3
## mean 0.0000000 0.000000 0.000000
## concentration 0.9637276 0.884906 0.216975
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -3.266252 -4.325395 -3.43482 -2.378091 -3.182405 -1.04517
##
## Transition probability matrix:
## ------------------------------
## state 1 state 2 state 3
## state 1 0.95113317 0.03628493 0.01258190
## state 2 0.02865101 0.88892161 0.08242738
## state 3 0.02977906 0.25240660 0.71781433
##
## Initial distribution:
## ---------------------
## state 1 state 2 state 3
## 0.53406404 0.40587611 0.06005985Model 1
The first model includes step, angle and no covariates. This is to see if the model improves including covariates later. It also considers there to be 4 states of movement.
# get the mean and sd values from the step column
#Look at the distribution of angle
hist(preppedData$angle)
#starting parameters
Par0_m1 <- list(step=c(150, 40, 4, 2, # means of 4 states
13, 8, 2, 2, # sd of first state and second state
0, 0, 0, 0.3), # zero mass parameter specification for first and second state
angle=c(0.99, 0.7, 0.3, 0.1)) # concentration high for first state, low for second.
m1 <- fitHMM(data = preppedData,
nbStates = 4,
dist = dist,
Par0 = Par0_m1,
estAngleMean = list(angle=FALSE)) # angle has a mean of zero## Value of the maximum log-likelihood: -19201.22
##
##
## step parameters:
## ----------------
## state 1 state 2 state 3 state 4
## mean 1.425397e+01 4.508344e+00 1.942493e+00 0.5616707
## sd 8.819541e+00 1.308504e+00 7.474716e-01 0.3512774
## zeromass 1.000241e-08 1.000187e-08 9.999848e-09 0.4414114
##
## angle parameters:
## -----------------
## state 1 state 2 state 3 state 4
## mean 0.0000000 0.0000000 0.000000 0.0000000
## concentration 0.9719025 0.9558443 0.921042 0.8012907
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 1 -> 3 1 -> 4 2 -> 1 2 -> 3 2 -> 4
## (Intercept) -2.933934 -4.456358 -3.520973 -3.476707 -2.439403 -3.556816
## 3 -> 1 3 -> 2 3 -> 4 4 -> 1 4 -> 2 4 -> 3
## (Intercept) -5.275932 -2.650842 -2.153787 -4.354006 -3.841725 -2.485444
##
## Transition probability matrix:
## ------------------------------
## state 1 state 2 state 3 state 4
## state 1 0.913773909 0.04860127 0.01060394 0.02702088
## state 2 0.026955905 0.87210461 0.07605877 0.02488072
## state 3 0.004290492 0.05923377 0.83910299 0.09737275
## state 4 0.011502508 0.01919878 0.07452450 0.89477421
##
## Initial distribution:
## ---------------------
## state 1 state 2 state 3 state 4
## 0.2077900 0.3355036 0.2682590 0.1884473Model 2
The second model includes step, angle and has temperature as its covariate. It still considers there to be 4 states of movement.
formula <- ~ temp
# initial parameter (obtained from m8)
Par0_m2 <- getPar0(model = m1, formula = formula)
# fit model
tic()
m2 <- fitHMM(data = preppedData,
nbStates = 4,
dist = dist,
Par0 = Par0_m2$Par,
formula = formula)## Value of the maximum log-likelihood: -19179.54
##
##
## step parameters:
## ----------------
## state 1 state 2 state 3 state 4
## mean 1.422352e+01 4.495562e+00 1.937739e+00 0.5616439
## sd 8.800235e+00 1.303349e+00 7.415505e-01 0.3510863
## zeromass 1.000459e-08 1.000319e-08 9.999689e-09 0.4411310
##
## angle parameters:
## -----------------
## state 1 state 2 state 3 state 4
## mean 0.0000000 0.0000000 0.0000000 0.000000
## concentration 0.9718758 0.9556607 0.9211487 0.801279
##
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
## 1 -> 2 1 -> 3 1 -> 4 2 -> 1 2 -> 3
## (Intercept) -2.8633693 -4.56169090 -3.508661611 -3.48139430 -2.441334673
## temp 0.2611145 -0.04810869 0.004056757 -0.07427068 -0.006451627
## 2 -> 4 3 -> 1 3 -> 2 3 -> 4 4 -> 1 4 -> 2
## (Intercept) -3.5423152 -5.295796 -2.6511962 -2.1736222 -4.393469 -3.8207720
## temp -0.1947821 -0.286750 -0.1475219 0.1611863 -0.485627 -0.4475546
## 4 -> 3
## (Intercept) -2.50306236
## temp 0.08456049
##
## Transition probability matrix (based on mean covariate values):
## ---------------------------------------------------------------
## state 1 state 2 state 3 state 4
## state 1 0.911197045 0.05200760 0.009516891 0.02727846
## state 2 0.026827420 0.87202561 0.075905118 0.02524186
## state 3 0.004214611 0.05933249 0.840799520 0.09565338
## state 4 0.011072269 0.01963160 0.073321219 0.89597491
##
## Initial distribution:
## ---------------------
## state 1 state 2 state 3 state 4
## 0.1997389 0.3397883 0.2662843 0.1941884Finally, lets have a look at the AIC values for both models to compare which is better
AIC(m1)## [1] 38464.43AIC(m2)## [1] 38445.07We can conclude that model 2 is better, even if the difference is little. If we had more time to run models, then we could include more environmental covariates, such as wind speed and solar, and other predictive variables such as different acs values.