###Conducting and plotting a basic PCA in Geomorph
#Load libraries
library(geomorph)
## Warning: package 'geomorph' was built under R version 4.1.2
## Loading required package: RRPP
## Warning: package 'RRPP' was built under R version 4.1.2
## Loading required package: rgl
## Warning: package 'rgl' was built under R version 4.1.2
## Loading required package: Matrix
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.2
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.1 v forcats 0.5.1
## Warning: package 'tidyr' was built under R version 4.1.2
## Warning: package 'readr' was built under R version 4.1.2
## Warning: package 'forcats' was built under R version 4.1.2
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x tidyr::pack() masks Matrix::pack()
## x tidyr::unpack() masks Matrix::unpack()
setwd("C:/1awinz/R_work/Geomorph")
#Import tps data file and put it in a matrix called "stickleback"
stickleback <- readland.tps("0_stickleback_evol_lab_data.tps", specID = "imageID")
##
## No curves detected; all points appear to be fixed landmarks.
#Perform the general procrustes analysis to aling landmarks
stickleback.gpa <- gpagen(stickleback, print.progress = FALSE)
#Provides a summary of the GPA and consensus configuration
summary(stickleback.gpa)
##
## Call:
## gpagen(A = stickleback, print.progress = FALSE)
##
##
##
## Generalized Procrustes Analysis
## with Partial Procrustes Superimposition
##
## 16 fixed landmarks
## 0 semilandmarks (sliders)
## 2-dimensional landmarks
## 2 GPA iterations to converge
##
##
## Consensus (mean) Configuration
##
## X Y
## 1 -0.41117662 0.000346648
## 2 -0.24207301 0.073780821
## 3 -0.13831304 0.088875817
## 4 -0.05843837 0.086127379
## 5 0.06281091 0.075125846
## 6 0.25195219 0.025702587
## 7 0.32144687 0.012590577
## 8 0.35235182 -0.004515709
## 9 0.31605382 -0.018456326
## 10 0.24999653 -0.021609268
## 11 0.12676763 -0.057565703
## 12 0.03455044 -0.073319393
## 13 -0.11310984 -0.076212332
## 14 -0.25414681 -0.067450865
## 15 -0.35861145 -0.049873849
## 16 -0.14006107 0.006453771
#To plot the aligned specimens
plotAllSpecimens(stickleback.gpa$coords)
#Perform a PCA on the gpa file
PCA_s <- gm.prcomp(stickleback.gpa$coords)
summary(PCA_s)
##
## Ordination type: Principal Component Analysis
## Centering by OLS mean
## Orthogonal projection of OLS residuals
## Number of observations: 90
## Number of vectors 28
##
## Importance of Components:
## Comp1 Comp2 Comp3 Comp4
## Eigenvalues 0.000643479 0.0006081554 0.0002788315 0.0001582414
## Proportion of Variance 0.314053935 0.2968140235 0.1360854457 0.0772307255
## Cumulative Proportion 0.314053935 0.6108679586 0.7469534043 0.8241841298
## Comp5 Comp6 Comp7 Comp8
## Eigenvalues 6.731505e-05 5.293737e-05 3.800709e-05 3.336165e-05
## Proportion of Variance 3.285353e-02 2.583642e-02 1.854960e-02 1.628236e-02
## Cumulative Proportion 8.570377e-01 8.828741e-01 9.014237e-01 9.177060e-01
## Comp9 Comp10 Comp11 Comp12
## Eigenvalues 0.0000247564 2.081827e-05 1.906948e-05 1.746122e-05
## Proportion of Variance 0.0120825168 1.016049e-02 9.306977e-03 8.522055e-03
## Cumulative Proportion 0.9297885449 9.399490e-01 9.492560e-01 9.577781e-01
## Comp13 Comp14 Comp15 Comp16
## Eigenvalues 1.572406e-05 1.164505e-05 9.851833e-06 8.243437e-06
## Proportion of Variance 7.674224e-03 5.683438e-03 4.808249e-03 4.023261e-03
## Cumulative Proportion 9.654523e-01 9.711357e-01 9.759440e-01 9.799672e-01
## Comp17 Comp18 Comp19 Comp20
## Eigenvalues 7.319103e-06 6.543606e-06 5.847578e-06 5.589374e-06
## Proportion of Variance 3.572134e-03 3.193648e-03 2.853947e-03 2.727929e-03
## Cumulative Proportion 9.835394e-01 9.867330e-01 9.895870e-01 9.923149e-01
## Comp21 Comp22 Comp23 Comp24
## Eigenvalues 3.796926e-06 3.345781e-06 2.471579e-06 2.284161e-06
## Proportion of Variance 1.853113e-03 1.632929e-03 1.206270e-03 1.114799e-03
## Cumulative Proportion 9.941680e-01 9.958009e-01 9.970072e-01 9.981220e-01
## Comp25 Comp26 Comp27 Comp28
## Eigenvalues 1.094531e-06 1.032518e-06 9.338544e-07 7.870060e-07
## Proportion of Variance 5.341927e-04 5.039268e-04 4.557735e-04 3.841032e-04
## Cumulative Proportion 9.986562e-01 9.991601e-01 9.996159e-01 1.000000e+00
#Plot the PCA results
plot(PCA_s, main = "Stickleback PCA")
#Save the matric of shape scores from the PCA_s object. The PCA_s objest consists of a series of lists. The PC scores are the 10th list in this file so can be referenced and saved via PCA_s[[10]]
pc_scores<-PCA_s[[10]]
#import grouping file into workspace
read.table("grp2.txt", header=T)
## grp type pop sex
## 1 AR_M anadromous AR male
## 2 AR_M anadromous AR male
## 3 AR_M anadromous AR male
## 4 AR_M anadromous AR male
## 5 AR_M anadromous AR male
## 6 AR_M anadromous AR male
## 7 AR_M anadromous AR male
## 8 AR_M anadromous AR male
## 9 AR_M anadromous AR male
## 10 AR_M anadromous AR male
## 11 AR_M anadromous AR male
## 12 AR_M anadromous AR male
## 13 AR_M anadromous AR male
## 14 AR_M anadromous AR male
## 15 AR_M anadromous AR male
## 16 AR_F anadromous AR female
## 17 AR_F anadromous AR female
## 18 AR_F anadromous AR female
## 19 AR_F anadromous AR female
## 20 AR_F anadromous AR female
## 21 AR_F anadromous AR female
## 22 AR_F anadromous AR female
## 23 AR_F anadromous AR female
## 24 AR_F anadromous AR female
## 25 AR_F anadromous AR female
## 26 AR_F anadromous AR female
## 27 AR_F anadromous AR female
## 28 AR_F anadromous AR female
## 29 AR_F anadromous AR female
## 30 AR_F anadromous AR female
## 31 Tern_M benthic Tern male
## 32 Tern_M benthic Tern male
## 33 Tern_M benthic Tern male
## 34 Tern_M benthic Tern male
## 35 Tern_M benthic Tern male
## 36 Tern_M benthic Tern male
## 37 Tern_M benthic Tern male
## 38 Tern_M benthic Tern male
## 39 Tern_M benthic Tern male
## 40 Tern_M benthic Tern male
## 41 Tern_M benthic Tern male
## 42 Tern_M benthic Tern male
## 43 Tern_M benthic Tern male
## 44 Tern_M benthic Tern male
## 45 Tern_M benthic Tern male
## 46 Tern_F benthic Tern female
## 47 Tern_F benthic Tern female
## 48 Tern_F benthic Tern female
## 49 Tern_F benthic Tern female
## 50 Tern_F benthic Tern female
## 51 Tern_F benthic Tern female
## 52 Tern_F benthic Tern female
## 53 Tern_F benthic Tern female
## 54 Tern_F benthic Tern female
## 55 Tern_F benthic Tern female
## 56 Tern_F benthic Tern female
## 57 Tern_F benthic Tern female
## 58 Tern_F benthic Tern female
## 59 Tern_F benthic Tern female
## 60 Tern_F benthic Tern female
## 61 Stormy_M limnetic Stormy male
## 62 Stormy_M limnetic Stormy male
## 63 Stormy_M limnetic Stormy male
## 64 Stormy_M limnetic Stormy male
## 65 Stormy_M limnetic Stormy male
## 66 Stormy_M limnetic Stormy male
## 67 Stormy_M limnetic Stormy male
## 68 Stormy_M limnetic Stormy male
## 69 Stormy_M limnetic Stormy male
## 70 Stormy_M limnetic Stormy male
## 71 Stormy_M limnetic Stormy male
## 72 Stormy_M limnetic Stormy male
## 73 Stormy_M limnetic Stormy male
## 74 Stormy_M limnetic Stormy male
## 75 Stormy_M limnetic Stormy male
## 76 Stormy_F limnetic Stormy female
## 77 Stormy_F limnetic Stormy female
## 78 Stormy_F limnetic Stormy female
## 79 Stormy_F limnetic Stormy female
## 80 Stormy_F limnetic Stormy female
## 81 Stormy_F limnetic Stormy female
## 82 Stormy_F limnetic Stormy female
## 83 Stormy_F limnetic Stormy female
## 84 Stormy_F limnetic Stormy female
## 85 Stormy_F limnetic Stormy female
## 86 Stormy_F limnetic Stormy female
## 87 Stormy_F limnetic Stormy female
## 88 Stormy_F limnetic Stormy female
## 89 Stormy_F limnetic Stormy female
## 90 Stormy_F limnetic Stormy female
classifier2=read.table("grp2.txt", header=T)
attach(classifier2)
names(classifier2)
## [1] "grp" "type" "pop" "sex"
#merge grouping matrix with matrix of pc scores
pc_scores2<-cbind(classifier2, pc_scores)
#Now you can plot with ggplot2
ggplot(pc_scores2, aes(x=Comp1, y=Comp2)) + geom_point(aes(color=type, shape=sex), size=3)+ xlab("Principal Component I") + ylab("Principal Component II")+theme_classic()
