Exploring the PCA published by Prentiss (1998) to understand the usefulness of the Sullivan and Rozen debitage typology

prentiss-1998-pca.R

#------------------------------------------------------
# Exploring the PCA published by Prentiss (1998) to understand the 
# usefulness of the Sullivan and Rozen debitage typology


# read in & tidy the data -----------------------------------------------
library(tidyverse)
# got these data from table 7 (p. 644) of https://www.jstor.org/stable/2694112
# OCR'd using https://tabula.technology/
prentiss <- readr::read_csv("tabula-Prentiss 1988.csv")

# the abbreviations are, from the paper:
# CF = Complete Flake; PF = Proximal Fragment; MDF = Medial/distal;Fragment; 
# NF = Nonorientable Fragment;SF = Split Flake; HH = Hard Hammer; SH = Soft Hammer; 
# PR = Pressure; PC = Prepared Core; UPC = Unprepared Core; BF = Biface; FL = Flake Tool.

# transpose to long format for analysis and plotting
prentiss_long <- 
  prentiss %>% 
  gather(variable, value, -var) %>% 
  spread(var, value) %>% 
  as.data.frame()

# split out the labels to use for plotting later,
# and convert them to factors because that's what the plotting functions need
prentiss_long_labels <- 
  prentiss_long %>% 
  separate(variable, 
           sep = "-", 
           into = c("technological_origin", "reduction_activity", "number"),
           remove = FALSE) %>% 
  mutate_at(vars(technological_origin, 
                 reduction_activity, 
                 number, 
                 variable), 
            as.factor)

# remove character cols for scatterplot panel
prentiss_long_rownames <- prentiss_long
row.names(prentiss_long_rownames) <- prentiss_long_rownames$variable
prentiss_long_rownames <- prentiss_long_rownames[,-1]

# scatterplot panel -----------------------------------------------

# show scatterplot panel to get an overview of how the variables relate
# to each other
library(GGally)
ggpairs(prentiss_long_rownames) +
  theme_bw()

# Compute & inspect the PCA -----------------------------------------------

# put the sample names as row names so 
# they show up later as labels on the plot
row.names(prentiss_long_labels) <- prentiss_long_labels$variable

library(FactoMineR)
# compute PCA
res.pca <- PCA(prentiss_long_labels, 
               quali.sup = 1:4,
               graph = FALSE)

# inspect eigenvalues, values >1 indicate that component captures more
# variability that any of the original measurement variables
eigenvalues <- res.pca$eig
head(eigenvalues[, 1:2])

# inspect distribution of PCs
library(factoextra)
fviz_screeplot(res.pca)

# Visualise output from the PCA ------------------------------------------

# plot variable loadings
fviz_pca_var(res.pca, 
             col.var="contrib") +
  scale_color_viridis_c() + 
  theme_bw()

# inspect PC1 vs PC2
fviz_pca_ind(res.pca, 
             geom = "text")  +
  theme_bw() +
  coord_equal()

# show ellipses for force application type
fviz_pca_biplot(res.pca,
                geom = "text",
                habillage = prentiss_long_labels$technological_origin,
                addEllipses=TRUE) +
  theme_bw() +
  coord_equal() +
  ggtitle("PCA showing artefacts grouped by technological origin")

# show ellipses for technology type
fviz_pca_biplot(res.pca,
                geom = "text",
                habillage = prentiss_long_labels$reduction_activity,
                addEllipses=TRUE) +
  theme_bw() +
  coord_equal() +
  ggtitle("PCA showing artefacts grouped by reduction activity")

tabula-Prentiss 1988.csv