3_prot1.Rmd

```{r, echo = FALSE, message=FALSE, warning=FALSE}
library(knitr)
opts_chunk$set(collapse = TRUE, comment = "#>")
library(png)
library(grid)
library(gridExtra)
library(cerUB)
library(biplot2d3d)
```

```{r, echo=FALSE, message=FALSE, warning=FALSE}
# this performs the initalization procedures
data(amphorae)

varCode <- code_variables(amphorae)

cleanAmphorae <- clean_and_format(
  amphorae,
  completion_variable = c(
    # The variable with completion info
    "CHARAC", 
    # the value indicating completion
    "complete"
  ), 
  categorical_columns = 1:112, 
  numerical_columns = 113:ncol(amphorae),
  # values converted to NA
  as_na = c("NULL", "indeterminate", "unfired"),
  # method for replacing NAs
  method = NULL, 
  # don't use the following variables
  columns_to_exclude = c("VOID_VESIC_MEGA", "VOID_VUGH_MEGA",
                         "VOID_CHAN_MEGA", "VOID_PLAN_MEGA",
                         "COAR_R_DAC_AND", "COAR_R_EVAP",
                         "COAR_R_CONGBREC", "COAR_R_SERP",
                         "COAR_C_SPL", "COAR_C_OPX",
                         "COAR_C_OL", "COAR_C_SIL",
                         "COAR_C_ST", "COAR_C_ZRN",
                         "COAR_C_PY", "FINE_C_OPX",
                         "FINE_C_ZRN"),
  # don't use the following observations
  # (Italic amphorae from Port Vendres 4)
  rows_to_exclude = c("PV4033", # PV4-IND4
                      "PV4017", # PV4-CAMP
                      # PV4-ITT
                      "PV4021", "PV4023", "PV4024", 
                      "PV4025", "PV4035", "PV4037",
                      # PV4-NAP
                      "PV4022", "PV4026", "PV4027", 
                      "PV4028", "PV4029", "PV4030",
                      "PV4036")
)
isShipwreck <-
  cleanAmphorae$Site_Name=="Cap del Vol" |
  cleanAmphorae$Site_Name=="Ullastres I" |
  cleanAmphorae$Site_Name=="Port-Vendres 4"

ProvenanceGroup <- c()
isTrueIND <- c()

# coerce the original group variables (factors) into character vectors
# so we can use stringr package to operate on them.
cleanAmphorae$FabricGroup <- 
  as.character(cleanAmphorae$FabricGroup)
cleanAmphorae$ChemReferenceGroup <- 
  as.character(cleanAmphorae$ChemReferenceGroup)

for (i in 1:nrow(cleanAmphorae)){
  groupChem <-
    stringr::str_split(cleanAmphorae$ChemReferenceGroup[i], "-")[[1]]
  groupFabric <-
    stringr::str_split(cleanAmphorae$FabricGroup[i], "-")[[1]]
  group <- ""
  isATrueInd <- FALSE

  if (groupChem[2] == "IND" || groupFabric[2] == "IND") {
    group <- cleanAmphorae$ChemReferenceGroup[i]
    if (!isShipwreck[i]) isATrueInd <- TRUE
    index <- 1
    for (j in 1:length(ProvenanceGroup)){
      if (ProvenanceGroup[j] == paste(group, index, sep = ""))
        index <- index + 1
    }
    group <- paste(group, index, sep = "")
    cleanAmphorae$ChemReferenceGroup[i] <- group
    cleanAmphorae$FabricGroup[i] <- group
  }
  else {
    if (groupChem[1] == "ULL" || 
        groupChem[1] == "PV4" || 
        groupChem[1] == "CDV") {
      group <- cleanAmphorae$ChemReferenceGroup[i]
    }
    else if (groupChem[1] == groupFabric[1]){
      group <- groupChem[1]
    }
  }
  ProvenanceGroup <- c(ProvenanceGroup, group[1])
  isTrueIND <- c(isTrueIND, isATrueInd)
}

factor_list <-
  list(
    Site = factor(cleanAmphorae$Site_Name[!isShipwreck]),
    FabricGroup = factor(cleanAmphorae$FabricGroup[!isShipwreck]),
    ChemGroup = factor(cleanAmphorae$ChemReferenceGroup[!isShipwreck]),
    ProvGroup = factor(ProvenanceGroup[!isShipwreck])
  )

factor_list_Shipwreck <-
  list(
    Site = factor(cleanAmphorae$Site_Name[!isTrueIND]),
    FabricGroup = factor(cleanAmphorae$FabricGroup[!isTrueIND]),
    ChemGroup = factor(cleanAmphorae$ChemReferenceGroup[!isTrueIND]),
    ProvGroup = factor(ProvenanceGroup[!isTrueIND])
  )

labels_code <- as.character(row.names(cleanAmphorae)) # using row names
labels_cross <- rep("+", nrow(cleanAmphorae)) # using +
labels_x <- rep(4, nrow(cleanAmphorae)) # using pch code
labels_point <- rep(20, nrow(cleanAmphorae)) # using pch code

labels_list <- list(
  Code = labels_code[!isShipwreck],
  Cross = labels_cross[!isShipwreck],
  X = labels_x[!isShipwreck],
  Point = labels_point[!isShipwreck]
)

labels_list_Shipwreck <- list(
  Code = labels_code[!isTrueIND],
  Cross = labels_cross[!isTrueIND],
  X = labels_x[!isTrueIND],
  Point = labels_point[!isTrueIND]
)

color_list <- list()

for (i in 1:length(factor_list)){
  cv <- rainbow(nlevels(factor_list[[i]]), v=.8)
  color_list[[i]] <- cv
  names(color_list)[i] = names(factor_list)[i]
}

color_list_Shipwreck <- list()

for (i in 1:length(factor_list_Shipwreck)){
  cv <- rainbow(nlevels(factor_list_Shipwreck[[i]]), v=.8)
  color_list_Shipwreck[[i]] <- cv
  names(color_list_Shipwreck)[i] = names(factor_list_Shipwreck)[i]
}

excep_cols <- c("INCLUS_DISTRIB","INCLUS_ORIENT","COAR_ROUNDNESS",
                "COAR_FORM","COAR_SPACING","COAR_SORTING","FINE_FORM")

chemVars16 <- c("Fe2O3","Al2O3","TiO2","MgO","CaO","SiO2",
                "Th","Nb","Zr","Y","Ce","Ga","V","Zn","Ni","Cr")


```

# Protocol 1 - Geochemical data {#prot1}

The following example applies protocol 1 to confirm the workshops' chemical reference groups.

Protocol 1 consist in: 

1. Select **_geochemical_** compositional data (CoDa);
2. Apply **_transformation_**;
3. Perform **_robust Principal Components Analysis_** (robPCA), implicitly using Euclidean distance;
4. Perform **_PERMANOVA & PERMDISP_** tests;

Last, search for outliers and re-do protocol excluding outliers.

NOTE: The [initial procedures](#init) must be ran at least once before any protocol can be applied.


## Ordination procedure

```{r}
prot1 <- apply_ordination(cleanAmphorae[!isShipwreck,], # no shipwrecks
                          "1", # select protocol 1
                          coda_override = chemVars16,
                          coda_transformation_method = "ILR")
```

The outcome is an *ordination object*. In this case, it is the output of **pcaCoDa** function in *robCompositions* package, in addition to several extra information, such as the transformed data ('transformed_data'), the distance matrix ('dist_matrix'), and the ready-to-plot texts indicating the fitness of the 2D/3D projections respect the distance matrix ('sub2d', 'sub3d'). The later are printed in the console once the object is created.

```{r}
class(prot1)

names(prot1)
```


## Simplify CoDa names

We may want to simplify the names of the transformed variables before plotting them in a biplot. The **transform_coda** function, which is called inside **apply_ordination** for protocol 1, generates composite names with format "transformationMethod-component" for all transformed variables (e.g., "CLR-Fe2O3"). The **simplify_coda_names** function replaces these names back to the shorter version (e.g., "Fe2O3"). However, you must always remember that the variables projected in biplots are not the originals but the transformed versions. This is particularly important when dealing with log-ratio variables since they contain information that goes beyond the original variable (i.e., divider).

```{r}
prot1 <- simplify_coda_names(prot1)
```


## Test the given chemical reference groups

Perform four tests (**anosim**, **betadisper**, **permdisp2**, and **permanova**) that assess the separation and uniformity of the given group factor. For more details on these tests, we refer to:

> Anderson, M.J., Walsh, D.C.I., 2013. PERMANOVA, ANOSIM, and the Mantel test in the face of heterogeneous dispersions: What null hypothesis are you testing? Ecol. Monogr. 83, 557-574. doi:10.1890/12-2010.1

The whole test batch may take several minutes depending on the size of the data matrix and the number of groups.

```{r}
prot1_tests <- test_groups(prot1$dist_matrix, factor_list$ChemGroup)
```

The tests outputs can be accessed by their name:

```{r}
names(prot1_tests)
```

The object also contains a "text" object, which is a function that generates a list of text lines for plotting the results of PERMANOVA and PERMDISP2 tests. It may feel confusing, but keep in mind that this "portable" function requires the same ordination object as an argument.

```{r, fig.width=3, fig.height=2, warning=FALSE}
displayTestText <- function(test_text) {
  par(mar = c(0, 0, 0, 0), fig = c(0.05, 0.9, 0.05, 0.9))
  plot.new()
  for (i in 1:length(test_text)) {
    
    # this is for calculating the vertical 
    # position of paragraphs and lines
    test_spacing_paragraph = 0.8
    test_spacing_line = 0.8
    
    first_line_pos_y <-
      1 - test_spacing_paragraph * ( (i - .9) / length(test_text) )
    
    pos_y <- first_line_pos_y
    
    if (length(test_text[[i]]) > 1) {
      
      next_paragraph_pos_y <-
        1 - test_spacing_paragraph * ( i / length(test_text) )
      
      for (j in 2:length((test_text[[i]])))
      {
        pos_y <-
          c(pos_y,
            first_line_pos_y - test_spacing_line * 
              ((j - 1) / (length((test_text[[i]])))) * 
              (first_line_pos_y - next_paragraph_pos_y)
          )
      }
    }
    # display a paragraph (a element of the list)
    text(0, pos_y, labels = test_text[[i]], cex = 0.8, pos = 4)
  }
}

displayTestText(prot1_tests$text(prot1_tests))
```

A rule-of-thumb for interpreting PERMANOVA and PERMDISP2 results is: if both p-values are low enough (e.g. < 0.05), the classification given is a good approximation of the data.

\pagebreak

## Biplots {#biplot}

Ordination objects are best represented in biplots, which simultaneously display the projections of observations (points) and variables (arrows) over the same space. There are several options for creating biplots in R, starting with the readily available **biplot** function:

```{r, fig.width=4.5, fig.height=4.5, fig.cap="Default biplot in R"}
biplot(prot1)
```


Although there are several options for customizing this default biplot function, we recommend the use of the *biplot2d3d* package. This package wraps a lot of possibilities in R.

```{r, message=FALSE}
library(biplot2d3d)
```

The *biplot2d3d* package use functions of other packages to allow the customization of virtually all aspects of a biplot. Another important feature of this package is the creation of three-dimensional interactive biplots using the *rgl* package.

\pagebreak


### Biplot 2D

You can consult all tuning options available in the **biplot_2d** function by calling up the help file:

```{r, eval=FALSE}
?biplot_2d
```

The default configuration will probably give you a much clearer picture than the **biplot** function, specially if your dataset contains more than 100 observations.

```{r, eval=FALSE}
biplot_2d(prot1)
```

```{r, echo=FALSE, fig.width=4.5, fig.height=4.5, fig.cap = "default 2D"}
biplot_2d(prot1,
          output_type = "preview")
```

\pagebreak

The default setting detaches the variables projections (arrows) and places them as a miniature in the bottom-right corner. In this format they may still be interpreted, much like the North when reading a map. Remember though: the more longer arrows you see, the less each one of them is reliable when comparing point values. Here, we were 'lucky' for getting two nearly orthogonal variables (CaO and MgO), which means, for instance, that those observations in the bottom-left corner are surely more calcareous than those in the top-right. See the [Appendix section](8_Appendix_biplot.html) for more details.

If detaching the arrows is not of your preference, you can disable this:

```{r, fig.width=4.5, fig.height=4.5, fig.cap = "detach_arrows = FALSE"}
biplot_2d(prot1, 
          detach_arrows = FALSE, 
          output_type = "preview")
```

\pagebreak

You can also visualize how the projection of points respond to a given typology (in this case, the chemical reference groups defined in previous studies).

```{r, fig.width=4.5, fig.height=4.5, fig.cap = "default with groups"}
biplot_2d(prot1, 
          groups = factor_list$ChemGroup,
          output_type = "preview")
```

\pagebreak

To get a prettier plot or match your research needs, you can play with the options given by the **biplot_2d**.  
Note that groups are by default marked using inertia ellipses. They can only be interpreted as confidence ellipses if each group can be assumed to be normally distributed in all variables considered (see more details on the scaling factor "group_ellipse_cex" in the help file). This is often not reasonable concerning groups that are either too small (< 10) or contain subgroups.  
In the argument "test_text" you can introduce the "text" function of the "prot1_tests" object.

```{r, fig.width=4.5, fig.height=4.5, fig.cap = "tunning appearance, groups with colors"}
biplot_2d(prot1, 
          groups = factor_list$ChemGroup, 
          group_color = color_list$ChemGroup,
          group_label_cex = 0.6,
          invert_coordinates = c(TRUE, TRUE),
          arrow_label_cex = 0.7,
          test_text = prot1_tests$text(prot1_tests),
          test_cex = 0.8,
          test_fig = c(0, 0.5, 0.65, .99),
          output_type = "preview")
```

\pagebreak

```{r, fig.width=4.5, fig.height=4.5, fig.cap = "labeled points"}
biplot_2d(prot1, 
          groups = factor_list$ChemGroup, 
          group_color = color_list$ChemGroup,
          group_label_cex = 0,
          group_star_cex = 0,
          group_ellipse_cex = 0,
          point_type = "label",
          point_label = row.names(cleanAmphorae)[!isShipwreck],
          point_label_cex = 0.5,
          invert_coordinates = c(TRUE, TRUE),
          arrow_label_cex = 0.7,
          test_text = prot1_tests$text(prot1_tests),
          test_cex = 0.8,
          test_fig = c(0, 0.5, 0.65, .99),
          output_type = "preview")
```

\pagebreak

It is also possible to save 2D biplots into various file formats (png, tiff, jpeg, eps):

```{r, eval=FALSE}
# better PNG version
biplot_2d(prot1,
          ordination_method = "PCA",
          invert_coordinates = c (TRUE,TRUE),
          grid_cex = 2.5,
          ylim = c(-.1,.1),
          point_type = "point",
          groups = factor_list$ChemGroup,
          group_color = color_list$ChemGroup,
          group_label_cex = 1.5,
          arrow_label_cex = 2,
          arrow_cex = 0.2,
          arrow_lwd = 2.5,
          arrow_fig = c(.6,.95,0,.35),
          subtitle_cex = 2.5,
          test_text = prot1_tests$text(prot1_tests),
          test_fig =c(0, 0.5, 0.62, .99),
          test_cex = 2,
          fitAnalysis_fig = c(0,.7,.05,.5),
          # saving settings
          file_name = "Prot1_Biplot2D",
          directory = directories$prot1,
          width = 1000, height = 1000,
          output_type = "png")
```

\pagebreak


### Biplot 3D

Most 2D options are also available when generating 3D biplots. Consult the help file for details.

```{r, eval=FALSE}
?biplot_3d
```

```{r, fig.width=4.5, fig.height=4.5}
biplot_3d(prot1,
          ordination_method = "PCA",
          groups = factor_list$ChemGroup,
          group_color = color_list$ChemGroup,
          point_type = "point",
          group_representation = "stars",
          star_centroid_radius = 0,
          star_label_cex = .8,
          test_text = prot1_tests$text(prot1_tests),
          test_cex = 1.25,
          test_fig = c(0, 0.5, 0.65, .99))

```
```{r, echo=FALSE}
rgl::rglwidget() # this allows you to see it in a .Rmd document.
```

\pagebreak

You can save an animated GIF and a PNG snapshot using the **animation** function. 

```{r, eval=FALSE}
biplot2d3d::animation(directory = directories$prot1,
                       file_name = "Prot1_Biplot3D")
```

You will need to install [ImageMagick](www.imagemagick.org) to be able to generate the GIF animation.

In these images, you get what you would see when running the biplot_3d function in a regular R session.

```{r Prot1Biplot3Danimation, eval=!params$ispdf, echo=FALSE, dependson='run.coevo.coeta-frames', fig.cap='Prot1_Biplot3D.gif'}
knitr::include_graphics("publication_appendices/Appendix_E_Animated_GIFs/E1_Prot1_Biplot3D.gif")
```

NOTE: Animated GIF will not be displayed in the pdf version of this document.

\pagebreak

## Comparing CoDa transformations

There is a lot of debate on which transformation is useful--or even *valid*--for analyzing geochemical compositions in Archaeometry. We show here how you can compare the results of applying different transformations to the same dataset. 

First, create different ordination objects for each type of CoDa transformation that you wish to compare:

```{r}
prot1_std <- apply_ordination(cleanAmphorae[!isShipwreck,],
                              "1", # select protocol 1
                              coda_override = chemVars16,
                              coda_transformation = "std")

prot1_log <- apply_ordination(cleanAmphorae[!isShipwreck,],
                              "1", # select protocol 1
                              coda_override = chemVars16,
                              coda_transformation = "log")

prot1_ALR <- apply_ordination(cleanAmphorae[!isShipwreck,],
                              "1", # select protocol 1
                              coda_override = chemVars16,
                              coda_transformation = "ALR",
                              # this is the divisor component
                              coda_alr_base = "Fe2O3")

prot1_CLR <- apply_ordination(cleanAmphorae[!isShipwreck,],
                              "1", # select protocol 1
                              coda_override = chemVars16,
                              coda_transformation = "CLR")
```

\pagebreak

Then, create the respective biplots:

```{r, warning=FALSE, message=FALSE, fig.width=4.5, fig.height=4.5}
biplot2d3d::biplot_2d(prot1_std,
                      groups = factor_list$ChemGroup,
                      group_color = color_list$ChemGroup,
                      group_label_cex = 0.6,
                      arrow_label_cex = 0.7,
                      ylim = c(-0.9, 0.6),
                      x_title = "Scaled and Centred",
                      x_title_cex = 1.5,
                      x_title_fig = c(0.05, 0.9, 0.85, 1),
                      output_type = "preview")
```

\pagebreak

```{r, warning=FALSE, message=FALSE, fig.width=4.5, fig.height=4.5}
biplot2d3d::biplot_2d(prot1_log,
                      groups = factor_list$ChemGroup,
                      group_color = color_list$ChemGroup,
                      group_label_cex = 0.6,
                      arrow_label_cex = 0.7,
                      ylim = c(-0.6, 0.6),
                      x_title = "Log-scaled",
                      x_title_cex = 1.5,
                      x_title_fig = c(0.05, 0.9, 0.85, 1),
                      output_type = "preview")
```

\pagebreak

```{r, warning=FALSE, message=FALSE, fig.width=4.5, fig.height=4.5}
biplot2d3d::biplot_2d(prot1_ALR,
                      groups = factor_list$ChemGroup,
                      group_color = color_list$ChemGroup,
                      group_label_cex = 0.6,
                      arrow_label_cex = 0.7,
                      ylim = c(-0.6, 0.6),
                      x_title = "Additive log-ratio",
                      x_title_cex = 1.5,
                      x_title_fig = c(0.05, 0.9, 0.85, 1),
                      output_type = "preview")
```

\pagebreak

```{r, warning=FALSE, message=FALSE, fig.width=4.5, fig.height=4.5}
biplot2d3d::biplot_2d(prot1_CLR,
                      groups = factor_list$ChemGroup,
                      group_color = color_list$ChemGroup,
                      group_label_cex = 0.6,
                      arrow_label_cex = 0.7,
                      ylim = c(-0.5, 0.4),
                      x_title = "Centred log-ratio",
                      x_title_cex = 1.5,
                      x_title_fig = c(0.05, 0.9, 0.85, 1),
                      output_type = "preview")
```