---
title: "Functional principal component analysis" 
subtitle: "Weddell seals dive shape" 
author: "Adélie ANTOINE & David NERINI"  
date: "(Last update 28-04-22)" 
output:
  html_document: 
    theme: readable 
    toc: yes 
    number_sections: yes 
    toc_float: yes 
    highlight: tango
  pdf_document:
    toc: yes 
    number_sections: yes 
    highlight: tango
header-includes: 
  - /renewcommand{/contentsname}{Table des matières} 
encoding: UTF-8 
---


```{r info, include=FALSE}
##############################################################################
# This code has been written by Adelie Antoine, with the help of David Nerini.
# 
# It is designed to perform functional principal component analysis and
# clustering analysis on dive shapes.
# 
# Last modification : April 2023
##############################################################################
```


```{r setup, include=FALSE}

knitr::opts_chunk$set(fig.align = "center", 
                      message = FALSE, 
                      warning = FALSE)
```

```{r packages, include=FALSE}
rm(list =ls()) 
library(fda)
library(rgl)
library(mclust)
library(factoextra)
library(cluster)
library(knitr)
library(kableExtra)
library(fpc)
library(tidyverse)
library (ggplot2)
library(cowplot)
```

/pagenumbering{arabic} 

# Data preparation

```{r data_upload, include=FALSE}
setwd("...")
data314a=read.csv(".../lw19_314a_summary_dives.csv",header=TRUE)
data314a$dive=paste(data314a$dive,"_1",sep="")
data314b=read.csv(".../lw19_314b_summary_dives.csv",header=TRUE)
data314b$dive=paste(data314b$dive,"_2",sep="")
data315a=read.csv(".../lw19_315a_summary_dives.csv",header=TRUE)
data315a$dive=paste(data315a$dive,"_3",sep="")


# dataframe preparation, called "plong"
numind=c(rep(1,length(unique(data314a$dive))),
         rep(2,length(unique(data314b$dive))),
         rep(4,length(unique(data315a$dive))))

data=rbind(data314a,data314b,data315a)
data[,1]=-data[,1]
plong=by(data[,2:1],data[,3],list)
```

Before applying fPCAs to the data, the 854 recorded dives must be visualised to make sure there is no sensor error or unexpected data.
For obvious reasons, the 854 dive profiles will not be provided here. A sample of 4 dives is presented below. 
```{r visu_test, echo= FALSE,out.width="70%", out.height="70%"}
nplong=length(plong)
rgeprof=matrix(NA,nplong,2)

for(k in (6:9)){ # iteration nb, chose any
  plong[[k]][,1]=plong[[k]][,1]-min(plong[[k]][,1])
  lonplong=length(plong[[k]][,1])
  plong[[k]][1,2]=0
  plong[[k]][lonplong,2]=0
  plot(plong[[k]],type="l",main=paste( 'Dive number', k, sep = ' ' ),
       xlab="Time (s)", ylab="Depth (m)") 
  locator(1)
}
```
These dive profiles illustrate the variety of dive types in the studied dataset, with duration, depth, skewness and overall shape differing among dives.

Following Godard _et al._, 2019 (in _Front. Mar. Sc._), fPCAs are lead on normalised dives. By reducing depth and duration to a 0 to 1 scale, dive structures are made comparable. The code is shown below.

```{r normalisation_dives, include= TRUE, eval = TRUE}
nbcoef=50
myb=create.bspline.basis(c(0,1),nbcoef)
penMat=eval.penalty(myb,2)
penMat=cbind(penMat,0,0)
penMat=rbind(penMat,0,0)

X.fd=fd(coef=NULL,basisobj=myb)
max_depth=NULL
duree=NULL
for(k in 1:nplong){
  cat("k=",k,fill=TRUE)
  tps=plong[[k]][,1]
  duree=c(duree,max(tps)-min(tps))
  tps=1/(max(tps)-min(tps))*(tps-min(tps))
  prof=plong[[k]][,2]
  max_depth=c(max_depth,min(prof))
  prof=-prof/min(prof)
  
  ### Compute the spline coefficients under constraints 
  k1=bsplineS(0,breaks=c(0,myb$params,1),nderiv=0)
  k1=t(k1)
  k1b=bsplineS(1,breaks=c(0,myb$params,1),nderiv=0)
  k1b=t(k1b)
  X=eval.basis(tps,myb)
  Mat=rbind(cbind(t(X)%*%X,k1,k1b),cbind(t(k1),0,0),cbind(t(k1b),0,0))
  ### Keep the spline coefficients, here KAPPA=1e-06
  achap=solve(Mat+0.000001*penMat)%*%rbind(t(X)%*%prof,0,0)
  X.fd$coefs=cbind(X.fd$coefs,as.matrix(achap[1:nbcoef]))
}
### Keep the coefficients
X.fd$coefs=X.fd$coefs[,-1]

# plot one normalised dive for visualisation
plot=plot(X.fd[68,])
```


Functional PCAs are then performed. 
```{r func_pca}
xpca.fd=pca.fd(X.fd,3)
```

Three dive shapes are found, with each archetypal profile shown below. 
``` {r component_plot, echo=FALSE, out.width="70%", out.height="70%"}
plot(xpca.fd)
```

The following plot is used to determine the number of components to keep for classification. It appears that 5 dimensions are necessary here. 
```{r nb_dim_plot, echo=FALSE,out.width="70%", out.height="70%"}
plot(xpca.fd$values/sum(xpca.fd$values)*100,type="h")
```

# Clustering methods application

To build clusters according to the shape of each dive, several methods are tested and compared. The robustness of clusters is tested by estimating consistency of groups created following each method.

Three methods were used:

- **Model based clustering**, where classification and density are estimated based on finite Gaussian mixture modeling,

- **K means clustering**, which assigns values to the cluster whose centroid is closest using Euclidean distance,

- **Hierarchical clustering**, that calculates the pairwise dissimilarity between each observation in the dataset and fuses observations into clusters through a classification tree (dendrogram).

## Model based clustering 

### Cluster construction

The data matrix is provided to the function `Mclust` to perform the clustering. The data matrix must be standardised.
A summary displaying information on the best model is available, giving sample size and overall probability of belonging to each cluster. Estimation of the suitability of models is performed using the Bayesian information criterion (BIC hereafter).

The best suited model, with higher BIC, comprises 5 components. Group volume, shape and orientation can vary (VVV model).

The corresponding code and its evaluation is displayed below.
```{r mclust}
mod1 <- Mclust(xpca.fd[["scores"]])
summary(mod1, parameters = TRUE) 
```

Pairwise dimension clustering can be plotted, to illustrate classification.

```{r mclust_pairwise_classif}
plot(mod1, what ="classification")
```

Each model's BIC can be plotted, depending on the number of components chosen. The following plot illustrates the necessity of a 5-dimension mixture.

The BIC summary presents BIC value for the top three models, for complementary information.

The ICL criterion can also be used, to balance bias induced by BIC that approximates the density rather than the number of clusters as such. Here, it might not be relevant to resort to the ICL criterion since it is efficient only in cases where there is very little overlap.

```{r BIC}
BIC <- mclustBIC(xpca.fd[["scores"]])
plot(BIC) # 5 dimensions mixture, VVV model with top BIC value
summary(BIC) # summary of top 3 models
```

### Exploration of the resulting clusters

Scrucca _et al._, 2010 (in _Statistics and computing_) suggests reducing dimensions to project clusters into a suitable subspace, in order to visualise clustering structures and geometric characteristics. The reduction is performed using `MclustDR` function. The estimated directions are defined as a set of linear combinations of the original features, ordered by importance. The ordering is carried out through the quantification of associated eigenvalues.

The `summary`function displays information on resulting dimensions. Eigenvalue of each remaining dimension is plotted, showing the classification performed with `Mclust`function. 

```{r red_model,out.width="70%", out.height="70%"}
drmod <- MclustDR(mod1, lambda = 1) # reduction of dimensions
summary(drmod)
plot(drmod, what = "evalues") 
```

Several exploratory plots can be obtained. They are different representations of the same information. Each color represents a cluster.

Pairwise representation of the classification, for dimensions 1 to 4:
```{r explo_plots_pairs,out.width="70%", out.height="70%"}
plot(drmod, what = "pairs") 
```

Cluster scatterplot with contour depending on distance to centroid, for dimensions 1 and 2:
```{r explo_plots_contour,out.width="70%", out.height="70%"}
plot(drmod, what = "contour") 
```

Cluster limits according to the model, for dimensions 1 and 2:
```{r explo_plots_classif,out.width="70%", out.height="70%"}
plot(drmod, what = "classification") 
```

Uncertainty zones (grey) according to the model, where adjacent points belong to different clusters, for dimensions 1 and 2:
```{r explo_plots_boundaries,out.width="70%", out.height="70%"}
plot(drmod, what = "boundaries") 
```

Distribution of observations on each dimension axis can be plotted, to compare cluster positions in the four dimensions subspace.
```{r density_plots,out.width="70%", out.height="70%"}
plot(drmod, dimens = 1, what = "density")
plot(drmod, dimens = 2, what = "density")
plot(drmod, dimens = 3, what = "density")
```

### Result extraction

Resulting clusters can be extracted from the model. Assessing probability of each data point to belong to its cluster is needed to estimate robustness of cluster predictions. 

In the following table, the first five columns represent the position of data points on each dimension axis. The last two columns (in bold) represents the cluster number to which the given dive belongs to, and probability that the said dive does belong to this cluster. Only the first 6 rows are presented.

```{r results_mclust}
cluster = mod1$classification
uncertainty = mod1$uncertainty
clust_results = data.frame(xpca.fd[["scores"]],cluster)
results = data.frame(clust_results,uncertainty)
kable(head(results)) %>%
  kable_styling("striped") %>%
  column_spec(4:5, bold = TRUE) %>%
  add_header_above(c("Dimension axis" = 3, "Clustering results" = 2))
```

Several threshold can be tested. Only dives with uncertainty below a given threshold are considered to belong to the corresponding cluster. For every tested threshold, the resulting number of classified dives is compared. 
```{r max_uncert_mclust, include = FALSE}
a1 = max(results$uncertainty)
```
The maximum uncertainty of a dive dive to belong to the cluster it was classified in is `r a1`.

The different thresholds tested are: 

- Uncertainty > `r a1`, called _a1_ in the following code chunk,
- Uncertainty > 0.5,
- Uncertainty > 0.4,
- Uncertainty > 0.3,
- Uncertainty > 0.2,
- Uncertainty > 0.1,
- Uncertainty > 0.05,
- Uncertainty > 0.04,
- Uncertainty > 0.03,
- Uncertainty > 0.02,
- Uncertainty > 0.01.

```{r thr_mclust}
results_thr1 = results[!(results$uncertainty>a1),]
results_thr2 = results[!(results$uncertainty>0.5),]
results_thr3 = results[!(results$uncertainty>0.4),]
results_thr4 = results[!(results$uncertainty>0.3),]
results_thr5 = results[!(results$uncertainty>0.2),]
results_thr6 = results[!(results$uncertainty>0.1),]
results_thr7 = results[!(results$uncertainty>0.05),]
results_thr8 = results[!(results$uncertainty>0.04),]
results_thr9 = results[!(results$uncertainty>0.03),]
results_thr10 = results[!(results$uncertainty>0.02),]
results_thr11 = results[!(results$uncertainty>0.01),]

thr_summary <- data.frame(var1=c(a1,0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.04,0.03, 0.02, 0.01),
                 var2=c(nrow(results_thr1),nrow(results_thr2),nrow(results_thr3),
                        nrow(results_thr4),nrow(results_thr5),nrow(results_thr6),
                        nrow(results_thr7),nrow(results_thr8),nrow(results_thr9),nrow(results_thr10),nrow(results_thr11)))

a2 <- c(a1,0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.04, 0.03, 0.02, 0.01)
clust1 <- numeric()
clust2 <- numeric()
clust3 <- numeric()
clust4 <- numeric()
clust5 <- numeric()
for(j in 1:11){ 
  k = a2[c(j)] 
  x1 = which(results$cluster == 1 & results$uncertainty<k)
  x2 = which(results$cluster == 2 & results$uncertainty<k)
  x3 = which(results$cluster == 3 & results$uncertainty<k)
  x4 = which(results$cluster == 4 & results$uncertainty<k)
  x5 = which(results$cluster == 5 & results$uncertainty<k)
clust1=rbind(clust1, length(x1))
clust2=rbind(clust2, length(x2))
clust3=rbind(clust3, length(x3))
clust4=rbind(clust4, length(x4))
clust5=rbind(clust5, length(x5))
}

thr_summary <- cbind(thr_summary,clust1,clust2,clust3,clust4,clust5)
names(thr_summary) <- c("Threshold value","Total number of dives", "Cluster 1",
                         "Cluster 2", "Cluster 3", "Cluster 4", "Cluster 5")

kable(thr_summary,align = 'c') %>%
  kable_styling("striped") %>%
  add_header_above(c(" " = 2, "Number of dives per cluster" = 5))
```



## K means clustering 

### Clustering construction 

To perform K means clustering, the data matrix must be normalised. It is therefore necessary to reduce and center the data frame.

```{r scale_kmeans}
df <- scale(xpca.fd[["scores"]]) 
head(df)
```

Contrary to usual method, the number of cluster will not be computed through `fviz_nbclust` function, nor by comparing the number of clusters to the gap statistic (using `clusGap` function). As the aim of using K means clustering here is to compare robustness of clusters resulting from model-based clustering, the same number of dimensions must be used. 

In the following analysis, a fixed number of clusters is provided to the function (5 dimensions). The `nstart` parameter is fixed at 150, to ensure stability of the model.

As shown in cluster plots examples, clusters overlay each other a lot, therefore results must be taken carefully. 
```{r kmeans}
km <- kmeans(df, centers = 5,nstart = 150) 
km
fviz_cluster(km, data = df,ellipse.type = "norm",axes = c(1, 2),) # to plot results
fviz_cluster(km, data = df,ellipse.type = "norm",axes = c(1, 3),) # to plot results
```

### Results extraction
The number of dives per cluster is presented in the following data frame.
```{r kmeans_size}
kmeans_results = as.data.frame(t(km$size))
names(kmeans_results) <- c("Cluster 1", "Cluster 2","Cluster 3","Cluster 4","Cluster 5")
kmeans_results
```

## Hierarchical clustering

### Clustering construction

As in K means clustering, the data matrix must be normalised to perform the analysis.

The first step of the analysis is to define the linkage method. To do so, we look at the agglomerative coefficient of each method, chosing the one that is closer to 1.

```{r linkage_method}
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

#function to compute agglomerative coefficient
ac <- function(x) {
  agnes(df, method = x)$ac
}

#calculate agglomerative coefficient for each clustering linkage method
sapply(m, ac) # best: closest to 1 (here Ward)
```
Ward's method seems to be the best suited for the analysis.

Hierarchical clustering is performed using Ward's minimum variance and a dendrogram is plotted. The dendrogram is then cut in the number of cluster needed, here 5.
```{r hclust}
clust <- agnes(df, method = "ward")
pltree(clust, cex = 0.6, hang = -1, main = "Dendrogram") 
#compute distance matrix
d <- dist(df, method = "euclidean")
#perform hierarchical clustering using Ward's method
final_clust <- hclust(d, method = "ward.D2" )
#cut the dendrogram into 5 clusters
groups <- cutree(final_clust, k=5)
```

### Result extraction 

The number of dives per cluster is given below. 
```{r result_hclust}
hmeans_results <- table(groups)
names(hmeans_results) <- c("Cluster 1","Cluster 2","Cluster 3","Cluster 4","Cluster 5")
hmeans_results
```

# Clustering results comparison

To estimate robustness of clusters obtained through the different methods used here, a comparison of the intersecting data points among clusters is performed.

Regarding model-based clustering, only the first threshold is tested to ensure all dives are classified, allowing for comparison with hierarchical and k means clustering. After computing a distance matrix, the function `cluster.stats` is applied to the three methods.


Pairwise comparison is applied to all three clustering methods. Corrected for chance rand index (`corrected.rand`) and variation of information index (`VI`) allow for estimating the similarity of information of both classification. 

Comparison of model-based and hierarchical clustering:
```{r pairwise_comparison_mh}
data_dist = xpca.fd[["scores"]]
comp_mh = cluster.stats(data_dist, results_thr1$cluster, alt.clustering = groups) #comparison mclust & hclust
comp_mh$corrected.rand
comp_mh$vi
```

Comparison of model-based and K means clustering:

```{r pairwise_comparison_mk}
comp_mk = cluster.stats(data_dist, results_thr1$cluster, alt.clustering = km$cluster) #comparison mclust & kclust
comp_mk$corrected.rand
comp_mk$vi
```


Comparison of K means and hierarchical clusterings:
```{r pairwise_comparison_kh}
comp_kh = cluster.stats(data_dist, groups, alt.clustering = km$cluster) #comparison hclust & kclust
comp_kh$corrected.rand
comp_kh$vi
```

References to analyse comparison results:

- https://cran.r-project.org/web/packages/fpc/fpc.pdf 
- Meila, M. (2007) Comparing clusterings-an information based distance, Journal of Multivariate Analysis, 98, 873-895.
- Gordon, A. D. (1999) Classification, 2nd ed. Chapman and Hall.


# Characteristics of clusters 

```{r mean_quantile_per_cluster}

time = seq(0, 1, by=(1/49))

idx_1 <- results_thr5[results_thr5$cluster == 1,] # extract dives in cluster 1
idx_1$index <- row.names(idx_1) # get index of dives 
vec_index_1 = as.numeric(c(idx_1$index)) # create vector to index in X.fd$coefs
mean_cl1 <- data.frame(NA_col = rep(NA, 50)) # create empty df to store results
for (k in 1:50) { 
  data = mean(X.fd$coefs[k, vec_index_1 ]) # get the mean
  q_2_5 = quantile(X.fd$coefs[k, vec_index_1 ], probs = c(.025)) # get quantile = 2.5%
  q_25 = quantile(X.fd$coefs[k, vec_index_1 ], probs = c(.25)) # get quantile = 25%
  q_75 = quantile(X.fd$coefs[k, vec_index_1 ], probs = c(.75)) # get quantile = 75%
  q_97_5 = quantile(X.fd$coefs[k, vec_index_1 ], probs = c(.975)) # get quantile = 97.5%
  mean_cl1[k , 1] <- data.frame(data)
  mean_cl1[k , 2] <- data.frame(q_2_5)
  mean_cl1[k , 3] <- data.frame(q_25)
  mean_cl1[k , 4] <- data.frame(q_75)
  mean_cl1[k , 5] <- data.frame(q_97_5)
}
mean_cl1[, 6] <- time
plot(mean_cl1$NA_col)


idx_2 <- results_thr5[results_thr5$cluster == 2,]
idx_2$index <- row.names(idx_2)
vec_index_2 = as.numeric(c(idx_2$index)) # create vector to index in X.fd$coefs
mean_cl2 <- data.frame(NA_col = rep(NA, 50)) # create empty df to store results
for (k in 1:50) { # loop to get mean of data 
  data = mean(X.fd$coefs[k, vec_index_2 ])
  mean_cl2[k , ] <- data
  q_2_5 = quantile(X.fd$coefs[k, vec_index_2 ], probs = c(.025)) # get quantile = 2.5%
  q_25 = quantile(X.fd$coefs[k, vec_index_2 ], probs = c(.25)) # get quantile = 25%
  q_75 = quantile(X.fd$coefs[k, vec_index_2 ], probs = c(.75)) # get quantile = 75%
  q_97_5 = quantile(X.fd$coefs[k, vec_index_2 ], probs = c(.975)) # get quantile = 97.5%
  mean_cl2[k , 1] <- data.frame(data)
  mean_cl2[k , 2] <- data.frame(q_2_5)
  mean_cl2[k , 3] <- data.frame(q_25)
  mean_cl2[k , 4] <- data.frame(q_75)
  mean_cl2[k , 5] <- data.frame(q_97_5)}
plot(mean_cl2$q_25)
mean_cl2[, 6] <- time

idx_3 <- results_thr5[results_thr5$cluster == 3,]
idx_3$index <- row.names(idx_3)
vec_index_3 = as.numeric(c(idx_3$index)) # create vector to index in X.fd$coefs
mean_cl3 <- data.frame(NA_col = rep(NA, 50)) # create empty df to store results
for (k in 1:50) { # loop to get mean of data 
  data = mean(X.fd$coefs[k, vec_index_3 ])
  q_2_5 = quantile(X.fd$coefs[k, vec_index_3 ], probs = c(.025)) # get quantile = 2.5%
  q_25 = quantile(X.fd$coefs[k, vec_index_3 ], probs = c(.25)) # get quantile = 25%
  q_75 = quantile(X.fd$coefs[k, vec_index_3 ], probs = c(.75)) # get quantile = 75%
  q_97_5 = quantile(X.fd$coefs[k, vec_index_3 ], probs = c(.975)) # get quantile = 97.5%
  mean_cl3[k , 1] <- data.frame(data)
  mean_cl3[k , 2] <- data.frame(q_2_5)
  mean_cl3[k , 3] <- data.frame(q_25)
  mean_cl3[k , 4] <- data.frame(q_75)
  mean_cl3[k , 5] <- data.frame(q_97_5)}
plot(mean_cl3$NA_col)
mean_cl3[, 6] <- time

idx_4 <- results_thr5[results_thr5$cluster == 4,]
idx_4$index <- row.names(idx_4)
vec_index_4 = as.numeric(c(idx_4$index)) # create vector to index in X.fd$coefs
mean_cl4 <- data.frame(NA_col = rep(NA, 50)) # create empty df to store results
for (k in 1:50) { # loop to get mean of data 
  data = mean(X.fd$coefs[k, vec_index_4 ])
  q_2_5 = quantile(X.fd$coefs[k, vec_index_4 ], probs = c(.025)) # get quantile = 2.5%
  q_25 = quantile(X.fd$coefs[k, vec_index_4 ], probs = c(.25)) # get quantile = 25%
  q_75 = quantile(X.fd$coefs[k, vec_index_4 ], probs = c(.75)) # get quantile = 75%
  q_97_5 = quantile(X.fd$coefs[k, vec_index_4 ], probs = c(.975)) # get quantile = 97.5%
  mean_cl4[k , 1] <- data.frame(data)
  mean_cl4[k , 2] <- data.frame(q_2_5)
  mean_cl4[k , 3] <- data.frame(q_25)
  mean_cl4[k , 4] <- data.frame(q_75)
  mean_cl4[k , 5] <- data.frame(q_97_5)}
plot(mean_cl4$NA_col)
mean_cl4[, 6] <- time

idx_5 <- results_thr5[results_thr5$cluster == 5,]
idx_5$index <- row.names(idx_5)
vec_index_5 = as.numeric(c(idx_5$index)) # create vector to index in X.fd$coefs
mean_cl5 <- data.frame(NA_col = rep(NA, 50)) # create empty df to store results
for (k in 1:50) { # loop to get mean of data 
  data = mean(X.fd$coefs[k, vec_index_5 ])
  q_2_5 = quantile(X.fd$coefs[k, vec_index_5 ], probs = c(.025)) # get quantile = 2.5%
  q_25 = quantile(X.fd$coefs[k, vec_index_5 ], probs = c(.25)) # get quantile = 25%
  q_75 = quantile(X.fd$coefs[k, vec_index_5 ], probs = c(.75)) # get quantile = 75%
  q_97_5 = quantile(X.fd$coefs[k, vec_index_5 ], probs = c(.975)) # get quantile = 97.5%
  mean_cl5[k , 1] <- data.frame(data)
  mean_cl5[k , 2] <- data.frame(q_2_5)
  mean_cl5[k , 3] <- data.frame(q_25)
  mean_cl5[k , 4] <- data.frame(q_75)
  mean_cl5[k , 5] <- data.frame(q_97_5)}
plot(mean_cl5$NA_col)
mean_cl5[, 6] <- time

```

Now that we have extracted the mean and quantiles, let's make graphs.
```{r plot_per_cluster}

c1 <- ggplot(mean_cl1, aes(x=time, y=NA_col)) +
  geom_line(color="red", size=2)+
  geom_line(aes(y=q_2_5),linetype="dotted", size=1) +
  geom_line(aes(y=q_25),linetype="dotted", size=1) +
  geom_line(aes(y=q_75),linetype="dotted", size=1) +
  geom_line(aes(y=q_97_5),linetype="dotted", size=1) +
  ggtitle("Cluster 1") +
  scale_y_continuous(n.breaks=3) +
  scale_x_continuous(n.breaks=3) +
  xlab("Time") + ylab("Depth") +
  theme_classic()+
  theme(axis.text.x=element_text(size=17), axis.text.y=element_text(size=17),
        axis.title = element_text(size=19),
        plot.title = element_text(size=21, face="bold"))
c1

c2 <- ggplot(mean_cl2, aes(x=time, y=NA_col)) +
 geom_line(color="red", size=2)+
  geom_line(aes(y=q_2_5),linetype="dotted", size=1) +
  geom_line(aes(y=q_25),linetype="dotted", size=1) +
  geom_line(aes(y=q_75),linetype="dotted", size=1) +
  geom_line(aes(y=q_97_5),linetype="dotted", size=1) +
  ggtitle("Cluster 2") +
  scale_y_continuous(n.breaks=3) +
  scale_x_continuous(n.breaks=3) +
  xlab("Time") + ylab("Depth") +
  theme_classic()+
  theme(axis.text.x=element_text(size=17), axis.text.y=element_text(size=17),
        axis.title = element_text(size=19),
        plot.title = element_text(size=21, face="bold"))
c2

c3 <- ggplot(mean_cl3, aes(x=time, y=NA_col)) +
 geom_line(color="red", size=2)+
  geom_line(aes(y=q_2_5),linetype="dotted", size=1) +
  geom_line(aes(y=q_25),linetype="dotted", size=1) +
  geom_line(aes(y=q_75),linetype="dotted", size=1) +
  geom_line(aes(y=q_97_5),linetype="dotted", size=1) +
  ggtitle("Cluster 3") +
  scale_y_continuous(n.breaks=3) +
  scale_x_continuous(n.breaks=3) +
  xlab("Time") + ylab("Depth") +
  theme_classic()+
  theme(axis.text.x=element_text(size=17), axis.text.y=element_text(size=17),
        axis.title = element_text(size=19),
        plot.title = element_text(size=21, face="bold"))
c3

c4 <- ggplot(mean_cl4, aes(x=time, y=NA_col)) +
 geom_line(color="red", size=2)+
  geom_line(aes(y=q_2_5),linetype="dotted", size=1) +
  geom_line(aes(y=q_25),linetype="dotted", size=1) +
  geom_line(aes(y=q_75),linetype="dotted", size=1) +
  geom_line(aes(y=q_97_5),linetype="dotted", size=1) +
  ggtitle("Cluster 4") +
  scale_y_continuous(n.breaks=3) +
  scale_x_continuous(n.breaks=3) +
  xlab("Time") + ylab("Depth") +
  theme_classic()+
  theme(axis.text.x=element_text(size=17), axis.text.y=element_text(size=17),
        axis.title = element_text(size=19),
        plot.title = element_text(size=21, face="bold"))
c4

c5 <- ggplot(mean_cl5, aes(x=time, y=NA_col)) +
 geom_line(color="red", size=2)+
  geom_line(aes(y=q_2_5),linetype="dotted", size=1) +
  geom_line(aes(y=q_25),linetype="dotted", size=1) +
  geom_line(aes(y=q_75),linetype="dotted", size=1) +
  geom_line(aes(y=q_97_5),linetype="dotted", size=1) +
  ggtitle("Cluster 5") +
  scale_y_continuous(n.breaks=3) +
  scale_x_continuous(n.breaks=3) +
  xlab("Time") + ylab("Depth") +
  theme_classic()+
  theme(axis.text.x=element_text(size=17), axis.text.y=element_text(size=17),
        axis.title = element_text(size=19),
        plot.title = element_text(size=21, face="bold"))
c5

legend <- get_legend( c1 )
prow <- plot_grid(
  c1 + theme(legend.position="none"),
  c2 + theme(legend.position="none"),
  c3 + theme(legend.position="none"),
  c4 + theme(legend.position="none"),
  c5 + theme(legend.position="none"),
  labels = c("A", "B", "C"),
  nrow = 2
)
prow


```



# Save files
To assess seal's behavior among dive groups, analyses must be conducted with MATLAB. The result file, containing cluster number for each dive and dive index for each seal, can be saved in .csv format.


# References 

Godard M, Manté C, Guinet C, Picard B, Nerini D. 2020. 'Diving Behavior of Mirounga leonina: A Functional Data Analysis Approach.' Front Mar Sci. 7-595. 

Kassambara A., Mundt F.. 2020. ‘factoextra: Extract and Visualize the Results of Multivariate Data Analyses’. R package version 1.0.7.

Maechler M., Rousseeuw P., Struyf A., Hubert M., Hornik K. 2019. ‘cluster: Cluster Analysis Basics and Extensions’. R package version 2.1.0.

Scrucca, L. 2010. 'Dimension reduction for model-based clustering.' Stat Comput 20, 471–484. 

Scrucca L., Fop M., Murphy T.B., Raftery A.E. 2016. ‘mclust 5: clustering, classification and density estimation using Gaussian finite mixture models.’ The R Journal, 8(1): 289–317.


