Raw Material Color

# compare colors
library(here)
library(colordistance)
library(recolorize)
library(dplyr)
library(vegan)
library(ape)

# get image paths
images <- colordistance::getImagePaths('images')

# generate color distance matrix
cdm <- imageClusterPipeline(images,
                            lower = FALSE,
                            upper = FALSE,
                            hist.bins = 3,
                            color.space = 'rgb')

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |======================================================================| 100%

# replace NAs with 0
cdm[is.na(cdm)] = 0
# export color distance matrix
write.csv(cdm, file = "color_distance_matrix.csv")

NMDS Ordination

# import color distance matrix
c.dat <- read.csv("color_distance_matrix.csv",
                   header = TRUE,
                   as.is = TRUE)

# omit 81
c.data <- c.dat %>% 
  select(c(X:X80)) %>% 
  filter(X %in% c(56:80))

# import qual data
e.data <- read.csv("qdata.csv",
                   header = TRUE,
                   as.is = TRUE)

# omit 81
e.data <- e.data %>% 
  filter(ColorGroup %in% c('A','B'))

# define add_image function
add_image <- function(obj,
                      x = NULL,
                      y = NULL,
                      width = NULL,
                      interpolate = TRUE,
                      angle = 0) {
  
  # get current plotting window parameters:
  usr <- graphics::par()$usr
  pin <- graphics::par()$pin
  
  # image dimensions and scaling factor:
  imdim <- dim(obj)
  sf <- imdim[1] / imdim[2]
  
  # set the width of the image (relative to x-axis)
  w <- width / (usr[2] - usr[1]) * pin[1]
  h <- w * sf
  hu <- h / pin[2] * (usr[4] - usr[3])
  
  # plot the image
  graphics::rasterImage(image = obj,
                        xleft = x - (width / 2), xright = x + (width / 2),
                        ybottom = y - (hu / 2), ytop = y + (hu/2),
                        interpolate = interpolate,
                        angle = angle)
}

# get image paths (sans img ID81)
images <- colordistance::getImagePaths('images2')

# cdm sans ID81
c.data2 <- as.matrix(c.dat %>% 
    filter(X %in% c(56:80)) %>% 
    select(c(X56:X80))
)

nmds_scores <- scores(metaMDS(comm = as.dist(c.data2)))

Run 0 stress 0.002359031 
Run 1 stress 0.002359087 
... Procrustes: rmse 9.379601e-05  max resid 0.0002291807 
... Similar to previous best
Run 2 stress 0.002359319 
... Procrustes: rmse 9.610158e-05  max resid 0.0002340066 
... Similar to previous best
Run 3 stress 0.002804939 
... Procrustes: rmse 0.00718958  max resid 0.01594965 
Run 4 stress 0.03008185 
Run 5 stress 0.006223433 
Run 6 stress 0.002760112 
... Procrustes: rmse 0.006642845  max resid 0.01468221 
Run 7 stress 0.002359052 
... Procrustes: rmse 7.494038e-05  max resid 0.0001833439 
... Similar to previous best
Run 8 stress 0.006182989 
Run 9 stress 0.006247482 
Run 10 stress 0.003441129 
Run 11 stress 0.002359048 
... Procrustes: rmse 7.519733e-05  max resid 0.0001856113 
... Similar to previous best
Run 12 stress 0.002372625 
... Procrustes: rmse 0.0008539385  max resid 0.001926306 
... Similar to previous best
Run 13 stress 0.03008237 
Run 14 stress 0.3473073 
Run 15 stress 0.006149061 
Run 16 stress 0.004521398 
Run 17 stress 0.006374702 
Run 18 stress 0.002359045 
... Procrustes: rmse 7.334519e-05  max resid 0.0001809859 
... Similar to previous best
Run 19 stress 0.03165429 
Run 20 stress 0.003066075 
*** Best solution repeated 6 times

# set plot parameters
plot(nmds_scores,
     xlim = c(-0.07, 0.12),
     ylim = c(-0.06, 0.06),
     cex = 1
)

#Create convex hulls that highlight point clusters based on grouping dataframe
ordihull(
  nmds_scores,
  e.data$ColorGroup,
  draw = c("polygon"),
  col = NULL,
  border = c("#798E87","#C27D38")
)

# add images
for (i in 1:length(images)) {
       
    # read image:
    img <- png::readPNG(images[i]) 
    # add image:
    add_image(img,
              x = nmds_scores[i, 1],
              y = nmds_scores[i, 2],
              width = 0.008)
}

Permutational MANOVA

# perMANOVA
# preliminaries
colnames(e.data) <- c('X','site','ColorGroup') # rename columns prior to join
elliptical <- left_join(e.data, c.data, by = "X") # left join by specimen number
elliptical.dist <- elliptical[,4:15] # distance matrix
set.seed(10) # make results reproducible

# model: biface color as a function of color group
el.biface.colour <- adonis2(elliptical.dist ~ ColorGroup,
                            data = elliptical,
                            permutations = 10000,
                            method = 'bray')

## does color differ by ColorGroup?
el.biface.colour

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 10000

adonis2(formula = elliptical.dist ~ ColorGroup, data = elliptical, permutations = 10000, method = "bray")
           Df SumOfSqs      R2      F Pr(>F)   
ColorGroup  1 0.017355 0.54754 12.102 0.0011 **
Residual   10 0.014341 0.45246                 
Total      11 0.031696 1.00000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1