Geomorphic Morphometrics on canid skulls

Running principal component and k-means clustering analyses

dir.create("./results", showWarnings = F)
dir.create("./results/figures", showWarnings = F)
source("./src/GM_stats.R")

Performing GPA

  |                                                                                      
  |                                                                                |   0%
  |                                                                                      
  |====================                                                            |  25%
  |                                                                                      
  |========================================                                        |  50%
  |                                                                                      
  |============================================================                    |  75%
  |                                                                                      
  |================================================================================| 100%

Making projections... Finished!

[1] -0.9999859
[1] 0.9999956
[1] -0.9999535
[1] -0.9999967
[1] -0.999442
[1] -0.9999388
[1] -1
[1] 0.9997545
[1] -0.9996835
[1] 0.9999969

Plotting PC and k-means results

K-means clustering plot:

Kennel-club groupings and PC results:

Discriminant Factors between groups

Comparison of landmarks between group means

# Extracting first 5 PCs and each kennel-club grouping scheme
skulls <- cbind(SlicerMorph.repc[, 1:5], 
                "UKC" = factor(SlicerMorph.repc$UKC, ordered = T, levels = UKC_group_order), 
                "AKC" = factor(SlicerMorph.repc$AKC, ordered = T, levels = AKC_group_order),
                "bitework" = factor(SlicerMorph.repc$Bitework),
                "scentwork" = factor(SlicerMorph.repc$Nosework))
skulls <- na.omit(skulls)

# Create a paired and ordered GPA-grouping object
mixdrop<-gpa$coords[,,order(dimnames(gpa$coords)[[3]])]
mixdrop<-mixdrop[,,-96]

reforge<-geomorph.data.frame(shape = mixdrop,
                             UKC = na.omit(factor(SlicerMorph.repc$UKC, ordered = T)),
                             AKC = na.omit(factor(SlicerMorph.repc$AKC, ordered = T)),
                             bitework = na.omit(factor(SlicerMorph.repc$Bitework[-96])),
                             scentwork = na.omit(factor(SlicerMorph.repc$Nosework[-96])))

# Create a function that subsets the mean landmarks for each group and does a pairwise subtraction and squaring of those means
meanforge <- function(target, group) {
  #sub-setting the means and sending them to a new object
  new.coords<-coords.subset(target[["shape"]], group = target[[group]])
  output_means<-lapply(new.coords, mshape)
  #perform the subtraction and squaring of the pairwise combinations
  result<- combn(output_means, 2, function(x) (x[[1]]-x[[2]])^2, simplify = FALSE)
  #rename the pairs so they're readable
  names(result)<-combn(names(output_means), 2, function(n) paste(n[1], "-", n[2]), simplify = TRUE)
  return(result)
}

#now run the function for each grouping scheme. Each member of the list will be named based on which two groups are compared
meanshape.UKC<-meanforge(reforge,"UKC")
meanshape.AKC<-meanforge(reforge,"AKC")
meanshape.nose<-meanforge(reforge,"scentwork")
meanshape.bite<-meanforge(reforge,"bitework")

Reformatting data for visualization:

calc_gpa_dist <- function(dat){
  require(tidyr)
  n <- length(names(dat))
  results <- matrix(NA, nrow = nrow(dat[[1]]), ncol = n)
  results <- as.data.frame(results)
  results <- cbind(seq(1,nrow(results)), results)
  colnames(results) <- c("landmark",names(dat))
  
  for(i in names(dat)){
    results[[i]] <- sqrt(dat[[i]][,"X"]^2 + dat[[i]][,"Y"]^2 + dat[[i]][,"Z"]^2)
  }
  results_long <- pivot_longer(results, cols = -"landmark", names_to = "comparison")
  results_long$landmark <- factor(results_long$landmark)
  return(results_long)
}

meanshape.AKC[["plotting"]] <- calc_gpa_dist(meanshape.AKC)
meanshape.UKC[["plotting"]] <- calc_gpa_dist(meanshape.UKC)
meanshape.bite[["plotting"]] <- calc_gpa_dist(meanshape.bite)
meanshape.bite[["cutoff"]] <- data.frame("mean" = mean(meanshape.bite[["plotting"]]$value), 
                                        "sd" = sd(meanshape.bite[["plotting"]]$value))
meanshape.nose[["plotting"]] <- calc_gpa_dist(meanshape.nose)
meanshape.nose[["cutoff"]] <- data.frame("mean" = mean(meanshape.nose[["plotting"]]$value), 
                                        "sd" = sd(meanshape.nose[["plotting"]]$value))

meanshape.AKC[["plotting"]]$short <- ifelse(grepl("toy", meanshape.AKC[["plotting"]]$comparison, 
                                                  ignore.case = T),
                                           "toy", 
                                           ifelse(grepl("natural", meanshape.AKC[["plotting"]]$comparison,
                                                        ignore.case = T),
                                                  "natural", "other"))
meanshape.AKC[["plotting"]]$short <- factor(meanshape.AKC[["plotting"]]$short, 
                                            ordered = T, levels = c("toy", "natural", "other"))
meanshape.AKC[["cutoff"]] <- data.frame("mean" = mean(meanshape.AKC[["plotting"]]$value), 
                                        "sd" = sd(meanshape.AKC[["plotting"]]$value))
meanshape.UKC[["plotting"]]$short <- ifelse(grepl("companion", meanshape.UKC[["plotting"]]$comparison,
                                                      ignore.case = T),
                                           "companion", 
                                           ifelse(grepl("natural", meanshape.UKC[["plotting"]]$comparison,
                                                      ignore.case = T),
                                           "natural", "other"))
meanshape.UKC[["plotting"]]$short <- factor(meanshape.UKC[["plotting"]]$short, 
                                            ordered = T, levels = c("companion", "natural", "other"))
meanshape.UKC[["cutoff"]] <- data.frame("mean" = mean(meanshape.UKC[["plotting"]]$value), 
                                        "sd" = sd(meanshape.UKC[["plotting"]]$value))

AKC & UKC groups:

p_UKC_compare <- ggplot(meanshape.UKC[["plotting"]], aes(landmark, value, color = short)) + 
  geom_point(position = position_jitter(width = 0.1)) + 
  geom_hline(yintercept = meanshape.UKC[["cutoff"]]$mean+2*meanshape.UKC[["cutoff"]]$sd,
             lty = 2) +
  scale_color_viridis_d(name = "Comparison\nagainst:") +
  annotate("segment", x = 1, xend = 20, y = 0.0085, yend = 0.0085,
             lty = 1, color = "red") +
  annotate("text", x = 10, y = 0.0075, label = "Dorsal features", color = "red") +
  annotate("segment", x = 21, xend = 38, y = 0.0085, yend = 0.0085,
             lty = 1, color = "blue") +
  annotate("text", x = 30, y = 0.0075, label = "Ventral features", color = "blue") +
  ylab("Group-mean landmark\ndistance squared") +
  ggtitle("UKC group comparisons") +
  theme_bw() #+
  #theme(legend.position = "bottom")
  #theme(legend.position = "inside", legend.position.inside = c(0.75,0.65),
  #      legend.background = element_rect(color = "gray30", fill = "white", linetype="solid",
  #                                       linewidth = 0.2))

p_AKC_compare <- ggplot(meanshape.AKC[["plotting"]], aes(landmark, value, color = short)) + 
  geom_point(position = position_jitter(width = 0.1)) + 
  geom_hline(yintercept = meanshape.AKC[["cutoff"]]$mean+2*meanshape.AKC[["cutoff"]]$sd,
             lty = 2) +
  scale_color_viridis_d(name = "Comparison\nagainst:") +
  annotate("segment", x = 1, xend = 20, y = 0.0085, yend = 0.0085,
             lty = 1, color = "red") +
  annotate("text", x = 10, y = 0.0075, label = "Dorsal features", color = "red") +
  annotate("segment", x = 21, xend = 38, y = 0.0085, yend = 0.0085,
             lty = 1, color = "blue") +
  annotate("text", x = 30, y = 0.0075, label = "Ventral features", color = "blue") +
  ylab("Group-mean landmark\ndistance squared") +
  ggtitle("AKC group comparisons") +
  theme_bw() #+
  #theme(legend.position = "bottom")
    #theme(legend.position = "inside", legend.position.inside = c(0.75,0.65),
    #    legend.background = element_rect(color = "gray30", fill = "white", linetype="solid",
    #                                     linewidth = 0.2))
p_AKC_compare / p_UKC_compare

Nosework and Bitework:

p_bite_compare <-  ggplot(meanshape.bite[["plotting"]], aes(landmark, value, color = comparison)) + 
  geom_point(position = position_jitter(width = 0.1)) + 
  geom_hline(yintercept = meanshape.bite[["cutoff"]]$mean+2*meanshape.bite[["cutoff"]]$sd,
             lty = 2) +
  annotate("segment", x = 1, xend = 20, y = 0.007, yend = 0.007,
             lty = 1, color = "red") +
  annotate("text", x = 10, y = 0.0065, label = "Dorsal features", color = "red") +
    annotate("segment", x = 21, xend = 38, y = 0.007, yend = 0.007,
             lty = 1, color = "blue") +
  annotate("text", x = 30, y = 0.0065, label = "Ventral features", color = "blue") +
  scale_color_viridis_d(name = "Group-group\ncomparison") +
  ylab("Group-mean landmark\ndistance squared") +
  ggtitle("Bite-work group comparisons") +
  theme_bw()
p_nose_compare <- ggplot(meanshape.nose[["plotting"]], aes(landmark, value, color = comparison)) + 
  geom_point(position = position_jitter(width = 0.1)) + 
  geom_hline(yintercept = meanshape.nose[["cutoff"]]$mean+2*meanshape.nose[["cutoff"]]$sd,
             lty = 2) +
  annotate("segment", x = 1, xend = 20, y = 0.007, yend = 0.007,
             lty = 1, color = "red") +
  annotate("text", x = 10, y = 0.0065, label = "Dorsal features", color = "red") +
    annotate("segment", x = 21, xend = 38, y = 0.007, yend = 0.007,
             lty = 1, color = "blue") +
  annotate("text", x = 30, y = 0.0065, label = "Ventral features", color = "blue") +
  scale_color_viridis_d(name = "Group-group\ncomparison") +
  ylab("Group-mean landmark\ndistance squared") +
  ggtitle("Scent-work group comparisons") +
  theme_bw()

p_bite_compare / p_nose_compare + plot_layout(guides = "collect")

MANOVA

man_UKC <- manova(as.matrix(skulls[,1:5])~skulls$UKC)
summary.aov(man_UKC, test="Pillai")

 Response PC.1 :
             Df  Sum Sq  Mean Sq F value    Pr(>F)    
skulls$UKC    8 0.65020 0.081275  35.372 < 2.2e-16 ***
Residuals   107 0.24585 0.002298                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response PC.2 :
             Df   Sum Sq   Mean Sq F value    Pr(>F)    
skulls$UKC    8 0.073427 0.0091783  6.0866 1.902e-06 ***
Residuals   107 0.161352 0.0015080                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response PC.3 :
             Df   Sum Sq    Mean Sq F value  Pr(>F)  
skulls$UKC    8 0.016265 0.00203313  2.4906 0.01612 *
Residuals   107 0.087348 0.00081633                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response PC.4 :
             Df   Sum Sq    Mean Sq F value    Pr(>F)    
skulls$UKC    8 0.024974 0.00312178  5.7231 4.625e-06 ***
Residuals   107 0.058365 0.00054547                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response PC.5 :
             Df   Sum Sq    Mean Sq F value Pr(>F)
skulls$UKC    8 0.000574 0.00007179  0.1633 0.9951
Residuals   107 0.047029 0.00043952               

#using Geomorph and doing pairwise comparisons
manUKC<-procD.lm(shape ~ UKC, data = reforge, RRPP = TRUE)
manAKC<-procD.lm(shape ~ AKC, data = reforge, RRPP = TRUE)
manBite<-procD.lm(shape ~ bitework, data = reforge, RRPP = TRUE)
manScent<-procD.lm(shape ~ scentwork, data = reforge, RRPP = TRUE)

hold.ukc<-manova(manUKC)
hold.akc<-manova(manAKC)
hold.bite<-manova(manBite)
hold.nose<-manova(manScent)

ukcpairs<-pairwise(manUKC, groups = reforge$UKC)
akcpairs<-pairwise(manAKC, groups = reforge$AKC)
bitepairs<-pairwise(manBite, groups = reforge$bitework)
nosepairs<-pairwise(manScent, groups = reforge$scentwork)


pvals.UKC<-summary(ukcpairs, test.type= "dist", confidence = 0.95, stat.table = TRUE)
pvals.UKC<-pvals.UKC$pairwise.tables$P[lower.tri(pvals.UKC$pairwise.tables$P)]
pvals.AKC<-summary(akcpairs, test.type= "dist", confidence = 0.95, stat.table = TRUE)
pvals.AKC<-pvals.AKC$pairwise.tables$P[lower.tri(pvals.AKC$pairwise.tables$P)]
pvals.bite<-summary(bitepairs, test.type= "dist", confidence = 0.95, stat.table = TRUE)
pvals.bite<-pvals.bite$pairwise.tables$P[lower.tri(pvals.bite$pairwise.tables$P)]
pvals.nose<-summary(nosepairs, test.type= "dist", confidence = 0.95, stat.table = TRUE)
pvals.nose<-pvals.nose$pairwise.tables$P[lower.tri(pvals.nose$pairwise.tables$P)]

pvals.ukc.adjusted<-p.adjust(pvals.UKC, method = "bonf")
pvals.akc.adjusted<-p.adjust(pvals.AKC, method = "bonf")
pvals.bite.adjusted<-p.adjust(pvals.bite, method = "bonf")
pvals.nose.adjusted<-p.adjust(pvals.nose, method = "bonf")

lda_UKC <- MASS::lda(skulls$UKC~as.matrix(skulls[,1:5]))
lda_UKC

Call:
lda(skulls$UKC ~ as.matrix(skulls[, 1:5]))

Prior probabilities of groups:
       NATURAL        Herding        Terrier      Companion        Gun dog       Guardian 
    0.19827586     0.23275862     0.08620690     0.12068966     0.05172414     0.09482759 
Northern Breed     Sighthound     Scenthound 
    0.02586207     0.14655172     0.04310345 

Group means:
               as.matrix(skulls[, 1:5])PC.1 as.matrix(skulls[, 1:5])PC.2
NATURAL                         -0.11869015                 -0.019063433
Herding                         -0.12101159                  0.026436494
Terrier                         -0.05852371                  0.001880539
Companion                        0.09983210                  0.018101105
Gun dog                         -0.11168214                  0.041234328
Guardian                        -0.06180271                  0.065382213
Northern Breed                  -0.06497801                  0.049661060
Sighthound                      -0.15480497                  0.029842740
Scenthound                      -0.12915491                  0.049368307
               as.matrix(skulls[, 1:5])PC.3 as.matrix(skulls[, 1:5])PC.4
NATURAL                          0.04817564                 -0.016530397
Herding                          0.04577986                  0.013615782
Terrier                          0.04485686                  0.032376419
Companion                        0.02331451                  0.011348484
Gun dog                          0.04624008                  0.007520818
Guardian                         0.06595086                  0.003014570
Northern Breed                   0.06220728                 -0.017512289
Sighthound                       0.02964807                  0.018704155
Scenthound                       0.04321946                  0.010813419
               as.matrix(skulls[, 1:5])PC.5
NATURAL                           0.1158473
Herding                           0.1190599
Terrier                           0.1152870
Companion                         0.1182048
Gun dog                           0.1204955
Guardian                          0.1207284
Northern Breed                    0.1170088
Sighthound                        0.1147247
Scenthound                        0.1131374

Coefficients of linear discriminants:
                                    LD1        LD2        LD3       LD4          LD5
as.matrix(skulls[, 1:5])PC.1  21.465233   2.072852   2.157478  1.232714   0.61490553
as.matrix(skulls[, 1:5])PC.2   2.130815 -20.151953  14.104686 -8.508683   1.48405135
as.matrix(skulls[, 1:5])PC.3 -13.131521   1.216403  21.217452 26.189753   4.17032936
as.matrix(skulls[, 1:5])PC.4  10.561161 -31.562092 -23.132983 17.258716  -0.02015756
as.matrix(skulls[, 1:5])PC.5   2.636558  -2.368104   6.558913  4.347214 -46.98480500

Proportion of trace:
   LD1    LD2    LD3    LD4    LD5 
0.7331 0.1664 0.0818 0.0172 0.0015 

Testing of Functional Groups

bite_yes_canine <- bite_data$canine[bite_data$Bitework == "yes"]
bite_no_canine <- bite_data$canine[bite_data$Bitework == "no"]
shapiro.test(bite_yes_canine)

    Shapiro-Wilk normality test

data:  bite_yes_canine
W = 0.95076, p-value = 0.5726

shapiro.test(bite_no_canine)

    Shapiro-Wilk normality test

data:  bite_no_canine
W = 0.97244, p-value = 0.08128

t.test(bite_yes_canine, bite_no_canine)

    Welch Two Sample t-test

data:  bite_yes_canine and bite_no_canine
t = -0.30306, df = 36.495, p-value = 0.7636
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.152725  6.771939
sample estimates:
mean of x mean of y 
 57.78286  58.97325 

bite_yes_carnassial <- bite_data$carnassial[bite_data$Bitework == "yes"]
bite_no_carnassial <- bite_data$carnassial[bite_data$Bitework == "no"]
shapiro.test(bite_yes_carnassial)

    Shapiro-Wilk normality test

data:  bite_yes_carnassial
W = 0.92608, p-value = 0.2687

shapiro.test(bite_no_carnassial)

    Shapiro-Wilk normality test

data:  bite_no_carnassial
W = 0.96684, p-value = 0.03577

t.test(bite_yes_carnassial, bite_no_carnassial)

    Welch Two Sample t-test

data:  bite_yes_carnassial and bite_no_carnassial
t = -0.52911, df = 42.805, p-value = 0.5995
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.417946  5.503589
sample estimates:
mean of x mean of y 
 57.12857  59.08575 

bite_natural_canine <- bite_data$canine[bite_data$Bitework == "NATURAL"]
shapiro.test(bite_natural_canine)

    Shapiro-Wilk normality test

data:  bite_natural_canine
W = 0.91337, p-value = 0.2672

bite_fox_canine <- bite_data$canine[bite_data$Bitework == "FOX"]
shapiro.test(bite_fox_canine)

    Shapiro-Wilk normality test

data:  bite_fox_canine
W = 0.82835, p-value = 0.03197

t.test(bite_yes_canine, bite_natural_canine)

    Welch Two Sample t-test

data:  bite_yes_canine and bite_natural_canine
t = -1.6371, df = 13.462, p-value = 0.1248
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -29.241240   3.979681
sample estimates:
mean of x mean of y 
 57.78286  70.41364 

t.test(bite_yes_canine, bite_fox_canine)

    Welch Two Sample t-test

data:  bite_yes_canine and bite_fox_canine
t = -0.85342, df = 10.859, p-value = 0.4119
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -29.77613  13.15585
sample estimates:
mean of x mean of y 
 57.78286  66.09300 

For the paper:

Figure 2: k-means cluster on PC plot

Figure 3: UKC comparison figure

Figure 4: AKC comparison figure

Figure 5: task-specific group results