2_init.Rmd

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

# Initial procedures {#init}

Load *cerUB* package.

```{r}
library(cerUB)
```


## Set directories for saving data

Set a list containing directories for easiness of reference:

```{r, results='hide', eval=FALSE}
directories <- list(
  # where you raw data is
  data = "data", 
  # where to save transformed compositional data (CoDa)
  transCoDa = "transformed_CoDa", 
  # where to save files concerning protocol workflow
  prot1 = "Protocol_1_geochemical_data",
  prot2 = "Protocol_2_petrographic_data",
  prot3 = "Protocol_3_geochemical_and_petrographic_data",
  prot4 = "Protocol_4_provenance_data",
  prot4_Shipwreck = "Protocols_4_provenance_data_with_shipwrecks"
)
```

Create the respective folders in the current R session working directory, if they do not exist:

```{r, results='hide', eval=FALSE}
lapply(directories, dir.create, showWarnings = FALSE)
```


## Read data

Load the amphorae dataset:

```{r}
data(amphorae)
```

Or, alternatively, your own dataset (e.g., CSV):

```{r, eval=FALSE}
dt <- cbind(read.csv(paste(directories$data,
                           "petrographic_data.csv",sep="/"),
                     # assuming the first column contains row names
                     row.names=1), 
            read.csv(paste(directories$data,
                           "geochemical_data.csv",sep="/"),
                     # assuming the first column contains row names
                     row.names=1)) 
```

Note that if you use your own dataset you must replace all references to "amphorae" with your data frame name (e.g. "dt").


## Codify petrographical variables

Create a two-column data frame containing the original names of petrographic variables and their respective codes:

```{r, warning=FALSE, message=FALSE}
varCode <- code_variables(amphorae)
```

 Petrographic variables must be named following *cerUB* naming system. Consult the documentation on the "amphorae" dataset by running:
 
```{r, eval=FALSE}
?amphorae
```

The result of the **code_variables** function is used in the **apply_ordination** function for protocols "2a", "2b", "3", and "4" (i.e., whenever there are ordinal input variables). The data frame containing the codification ('varCode') is attached as "ordination_object$variable_tags" to the resulting ordination object. 

```{r, echo=FALSE, results='asis'}
knitr::kable(rbind(head(varCode, 12),c("...","...")))
```


## Clean and format data

Cleaning and format procedures, including coercing variables as numeric or factor, excluding columns (constants, perturbed, unreliable) and rows (incomplete data, outliers).

```{r}
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")
)
```
```{r, echo=FALSE, results='asis'}
knitr::kable(cbind(" " = c("amphorae", "cleanAmphorae"),
                   Variables = c(ncol(amphorae), ncol(cleanAmphorae)),
                   Observations = c(nrow(amphorae), nrow(cleanAmphorae))))
```


## Subsetting criteria

Build vector indicating whether each observation is from a shipwreck:

```{r}
isShipwreck <-
  cleanAmphorae$Site_Name=="Cap del Vol" |
  cleanAmphorae$Site_Name=="Ullastres I" |
  cleanAmphorae$Site_Name=="Port-Vendres 4"
```
```{r, echo=FALSE, results='asis'}
knitr::kable(data.frame(cbind(Workshops = c(table(isShipwreck)["FALSE"]),
                   Shipwrecks = c(table(isShipwreck)["TRUE"])),
                   row.names = c(" ")))
```

Build vectors indicating provenance group and whether observations are true outliers (IND, observations with no group assigned). Also, reformat "FabricGroup" and "ChemReferenceGroup", so true outliers are singled out separately and not as a extra group.

```{r}
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)
}
```

```{r, echo=FALSE, results='asis'}
knitr::kable(data.frame(cbind(Assigned = c(table(isTrueIND)["FALSE"]),
                   Outliers = c(table(isTrueIND)["TRUE"])),
                   row.names = c(" ")))
```


## Organizing groups

Build lists of named group factors for easiness of reference.

First, create a list aiming to define workshops productions, so no shipwrecks:

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

```{r, echo=FALSE, results='asis'}
knitr::kable(data.frame(factor_list,
                        row.names = row.names(cleanAmphorae)[!isShipwreck])[seq(1, nrow(cleanAmphorae[!isShipwreck,]), length.out = 12),])
```

Then, create a second list, aiming to assign shipwreck observations to workshop productions, so with shipwreck samples but no true outliers:

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


## Helper objects for plotting 

Build lists of named point types vectors for easiness of reference.

Create point type vectors for the whole dataset:

```{r}
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
```

Create a list aiming to define workshops productions:

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

Create a list aiming to assign shipwreck observations to workshop productions:

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

Build two other lists containing named group color vectors, picking different colors within the **rainbow** palette:

```{r}
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]
}
```

To visualize the colors:

```{r, eval=FALSE}
par(mar = c(0,0,0,0))
barplot(rep(100, nlevels(factor_list$ProvGroup)), 
        col = color_list$ProvGroup, 
        axes = FALSE)
```
```{r, echo=FALSE, fig.width=7, fig.height=1}
par(mar = c(0,0,0,0))
barplot(rep(100, nlevels(factor_list$ProvGroup)), 
        col = color_list$ProvGroup, 
        axes = FALSE)
```


## Enunciate exception columns

Create a vector that enunciate which ordinal variables have "none" as a exceptional value when calculating the distance between values.

```{r}
excep_cols <- c("INCLUS_DISTRIB","INCLUS_ORIENT","COAR_ROUNDNESS",
                "COAR_FORM","COAR_SPACING","COAR_SORTING","FINE_FORM")
```

In order to understand this "exceptional value" feature, compare the levels of regular and exceptional variables. However, to do that at this point you must re-assure the order of petrographic variables (i.e. format factors levels):

```{r}
cleanAmphorae <- order_petro(cleanAmphorae)
```

This step is not necessary for applying the protocols because the apply_ordination function already does it internally, before calculating distances.

```{r, echo=FALSE, results='asis'}
pander::pander(
  cbind(
    "Variable" = 
      c("INCLUS_DISTRIB", "TEMP", "COAR_FREQ", "COAR_ROUNDNESS",
        "COAR_R_CALS", "FINE_FORM", "FINE_C_QTZ"),
    "Values" = 
      c(paste(levels(cleanAmphorae$INCLUS_DISTRIB), collapse = ", "),
        paste(levels(cleanAmphorae$TEMP), collapse = ", "),
        paste(levels(cleanAmphorae$COAR_FREQ), collapse = ", "),
        paste(levels(cleanAmphorae$COAR_ROUNDNESS), collapse = ", "),
        paste(levels(cleanAmphorae$COAR_R_CALS), collapse = ", "),
        paste(levels(cleanAmphorae$FINE_FORM), collapse = ", "),
        paste(levels(cleanAmphorae$FINE_C_QTZ), collapse = ", ")
      )
  ), split.cells = 30, split.table = Inf
)
```


## Choose geochemical data

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


## Save transformed geochemical data to file (optional)

There is no need to save it in the environment, because **apply_ordination** will transform
the data internally and save the results in "ordination_object$transformed_data",
when applicable. 

```{r, eval=FALSE}
write(transform_coda(cleanAmphorae,
                     coda_variables = chemVars16,
                     method = c("CLR")),
      file = paste(directories$transCoDa, 
                   "transAmphorae_clr.csv", 
                   sep = "/"))
```

In the output table, columns will be ordered as:

1. variables not transformed,
2. Raw version of the selected variables,
3. Transformed version of the selected variables.

## Other CoDa packages

Be aware that compositional data (CoDa) analysis can be much more complex than what cerUB currently allows for. For more possibilities, you may explore other R packages: [*compositions*](https://cran.r-project.org/web/packages/compositions/index.html), [*zCompositions*](https://cran.r-project.org/web/packages/zCompositions/index.html), and [*robCompositions*](https://cran.r-project.org/web/packages/robCompositions/index.html).

However, before jumping into using more complex techniques, we do recommend a deeper introduction
to CoDa:

> Pawlowsky-Glahn, V., Buccianti, A., 2011. Compositional data analysis: theory and applications. Wiley.