Principal component analysis (PCA)
Anastasia Bernat
8/12/2020 - 12/16/2020
Winter 2020 Flight Trials
Flight Trials Winter 2020 were conducted from 2/17/2020 - 3/10/2020. Soapberry bugs were flight tested twice for multiple hours in the flight mill and observed from 8 AM to (5-8 PM) each day.
Size-adjusted CPCA (notes)
PCAs are good to…
Isolate changes in allometry: “Morphological relationships change with overall body size and body size often varies among populations. Therefore, quantitative analyses of individual traits from organisms in different populations or environments (e.g. in studies of phenotypic plasticity) often adjust for differences in body size to isolate changes in allometry.” (McCoy et al, 2006)
Help define scaling of traits with body size with PC1: “For most studies of morphological plasticity and size correction, it is only the first CPC (CPC1) that is of interest because it describes the general scaling of traits with body size.” (McCoy et al, 2006)
Define size patterns and size corrections: “If the first principal components are not shared (Fig. 1c), then the patterns of morphological variation are so fundamentally different that there can be no globally applied size correction because ‘size’ does not have a common meaning across groups. On the other hand, if all principal components are common to all groups (Fig. 1a), then the groups show identical patterns of within-group covariation, indicating the same allometries of traits with size in all groups.” (McCoy et al, 2006) AB: Looking at the figures….I would say differences in size but not shape.
Define common body size dimensions to evaluate group differences: “As long as the groups to be compared share CPC1 (Fig. 1a, b) and all the traits have strong loadings on CPC1, then it can be interpreted as a common body size dimension: i.e., a ‘size-axis’ that can be used to evaluate between-group differences in morphology that go beyond simple (isometric) changes in size” (McCoy et al, 2006)
Correct for the effects of a growth axis: “…where the main goal is to correct for the effects of a growth axis that explains a large fraction of the total variation (e.g., >90%)” (McCoy et al, 2006)
PCA Functions
Round <- function(number){
# for plotting
x <- round(number, 1)
if(x%%1 == 0){
return(paste(as.character(x), ".0", sep = ""))
}
else{
return(x)
}
}PCA_graphs <- function(dataset, PCA_title){
# cos2 and the alpha.var: alpha.var colours variables by cos2
# (importance of most important PC to variable),
# see https://personal.utdallas.edu/~herve/abdi-awPCA2010.pdf
GFpca <- PCA(dataset, scale.unit = TRUE, graph = TRUE, ncp = 10)
eig.val <- get_eigenvalue(GFpca)
var.val <- GFpca$var
print(eig.val) #will only show in in console
print(var.val)
scree <- fviz_eig(GFpca, addlabels = TRUE, ylim = c(0, 100))
print(scree)
labX <- paste("PC1 (", Round(eig.val[1, 2]), "%)", sep = "")
labY <- paste("PC1 (", Round(eig.val[2, 2]), "%)", sep = "")
leplot <- fviz_pca_biplot(GFpca, geom.id = c("point"),
geom.var = c("arrow", "text"),
alpha.var = "cos2",
label = "var", repel = T,
col.ind = "gray", col.var = "black")
print(leplot)
ggpubr::ggpar(leplot, title = PCA_title, xlab = labX, ylab = labY,
ggtheme = theme_classic(), font.main = c(20, "bold"),
font.x = 14, font.y = 14, font.tickslab = 12
)
D = cor(dataset)
test <- cor.mtest(dataset)$p
par(mfrow=c(1,2))
corrplot.mixed(D,lower.col = "black", number.cex = .7, p.mat=test, sig.level=0.05)
corrplot.mixed(D,lower.col = "black", number.cex = .7)
return(GFpca)
}Reading and Standardizaing Data
source_path = "~/Desktop/git_repositories/SBB-dispersal/avbernat_working_on/Rsrc/"
script_names = c("center_flight_data.R", # Re-centers data
"clean_flight_data.R" # Loads and cleans data
)
for (script in script_names) {
path = paste0(source_path, script)
source(path)
}##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, uniondata <- read_flight_data("data/all_flight_data-Winter2020.csv")
data_all <- data[[1]]
data_tested <- data[[2]]# Filter out duplicates (PCA applies to data tables where rows are considered as individuals)
data_tested <- data_tested[data_tested$trial_type == "T1",]
# Standardize variables
data_tested$thorax_s <- (data_tested$thorax-mean(data_tested$thorax, na.rm=TRUE)) / sd(data_tested$thorax, na.rm=TRUE)
data_tested$body_s <- (data_tested$body-mean(data_tested$body, na.rm=TRUE)) / sd(data_tested$body, na.rm=TRUE)
data_tested$wing_s <- (data_tested$wing-mean(data_tested$wing, na.rm=TRUE)) / sd(data_tested$wing, na.rm=TRUE)
data_tested$beak_s <- (data_tested$beak-mean(data_tested$beak, na.rm=TRUE)) / sd(data_tested$beak, na.rm=TRUE)
data_tested$mass_s <- (data_tested$mass-mean(data_tested$mass, na.rm=TRUE)) / sd(data_tested$mass, na.rm=TRUE)Morphology
Thorax, body, wing, and beak lengths only considered.
d <- data_tested %>%
select(thorax_s, body_s, wing_s, beak_s) %>%
filter(!is.na(body_s))
colnames(d) <- c("thorax", "body", "wing", "beak")
MorphPCA = PCA_graphs(d, "(a) ")
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 3.50143036 87.535759 87.53576
## Dim.2 0.37328713 9.332178 96.86794
## Dim.3 0.10977100 2.744275 99.61221
## Dim.4 0.01551152 0.387788 100.00000
## $coord
## Dim.1 Dim.2 Dim.3 Dim.4
## thorax 0.9593090 -0.03700249 -0.27928757 0.01885528
## body 0.9826965 -0.14281597 0.06407406 -0.09902406
## wing 0.9451532 -0.28258460 0.14715865 0.07194252
## beak 0.8497939 0.52121731 0.07751322 0.01321008
##
## $cor
## Dim.1 Dim.2 Dim.3 Dim.4
## thorax 0.9593090 -0.03700249 -0.27928757 0.01885528
## body 0.9826965 -0.14281597 0.06407406 -0.09902406
## wing 0.9451532 -0.28258460 0.14715865 0.07194252
## beak 0.8497939 0.52121731 0.07751322 0.01321008
##
## $cos2
## Dim.1 Dim.2 Dim.3 Dim.4
## thorax 0.9202737 0.001369184 0.078001546 0.0003555218
## body 0.9656923 0.020396401 0.004105485 0.0098057644
## wing 0.8933146 0.079854056 0.021655668 0.0051757262
## beak 0.7221497 0.271667486 0.006008299 0.0001745063
##
## $contrib
## Dim.1 Dim.2 Dim.3 Dim.4
## thorax 26.28279 0.3667911 71.058429 2.291986
## body 27.57994 5.4639979 3.740045 63.216018
## wing 25.51285 21.3921269 19.728041 33.366985
## beak 20.62442 72.7770840 5.473485 1.125011
vars <- MorphPCA$ind
par(mfrow=c(2,3))
#Correlations
m <- lm(d$thorax ~ vars$coord[,1])
plot(vars$coord[,1], d$thorax, xlab="PC1", ylab="thorax (mm)")
cor_val <- paste("R = ", round(cor(vars$coord[,1], d$thorax), 2))
text(-2,2, cor_val)
abline(m, lty=2)
m <- lm(d$body ~ vars$coord[,1])
plot(vars$coord[,1], d$body, xlab="PC1", ylab="body (mm)")
cor_val <- paste("R = ", round(cor(vars$coord[,1], d$body), 2))
text(-2,2, cor_val)
abline(m, lty=2)
m <- lm(d$wing ~ vars$coord[,1])
plot(vars$coord[,1], d$wing, xlab="PC1", ylab="wing (mm)")
cor_val <- paste("R = ", round(cor(vars$coord[,1], d$wing), 2))
text(-2,2, cor_val)
abline(m, lty=2)
m <- lm(d$beak ~ vars$coord[,1])
plot(vars$coord[,1], d$beak, xlab="PC1", ylab="beak (mm)")
cor_val <- paste("R = ", round(cor(vars$coord[,1], d$beak), 2))
text(-2,2, cor_val)
abline(m, lty=2)
m <- lm(d$beak ~ vars$coord[,2])
plot(vars$coord[,2], d$beak, xlab="PC2", ylab="beak (mm)")
cor_val <- paste("R = ", round(cor(vars$coord[,2], d$beak), 2))
text(-0.5,3, cor_val)
abline(m, lty=2)
#Plot by host
host <- data_tested %>%
select(host_plant)
d$host_plant <- host$host_plant
fviz_pca_ind(MorphPCA,
geom.ind = "point",
col.ind = d$host_plant,
legend.title = "Host Plant",
addEllipses=TRUE)
# Plot by sex
sex <- data_tested %>%
select(sex)
d$sex <- sex$sex
fviz_pca_ind(MorphPCA,
geom.ind = "point",
col.ind = d$sex,
legend.title = "Sex",
addEllipses=TRUE)
# Plot by distance from the sympatric zone
sym_dist <- data_tested %>%
select(sym_dist)
d$sym_dist <- sym_dist$sym_dist
fviz_pca_ind(MorphPCA,
geom.ind = "point",
col.ind = d$sym_dist,
legend.title = "Sympatric Zone Distance")
# Plot by population
pop <- data_tested %>%
select(population)
d$population <- pop$pop
fviz_pca_ind(MorphPCA,
geom.ind = "point",
col.ind = d$population,
legend.title = "Population",
addEllipses=TRUE)
Morphology w/ Mass
d <- data_tested %>%
select(thorax_s, body_s, wing_s, beak_s, mass_s) %>%
filter(!is.na(mass_s))
colnames(d) <- c("thorax", "body", "wing", "beak", "mass")
MorphPCA = PCA_graphs(d, "(b) ")
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 4.10248965 82.0497929 82.04979
## Dim.2 0.51045917 10.2091834 92.25898
## Dim.3 0.26872036 5.3744071 97.63338
## Dim.4 0.10447905 2.0895811 99.72296
## Dim.5 0.01385177 0.2770355 100.00000
## $coord
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## thorax 0.9544052 -0.1068106 0.01326175 -0.27806416 0.014372683
## body 0.9684503 -0.2190542 0.03179998 0.06598983 -0.093559921
## wing 0.9173345 -0.3613220 0.06701864 0.13672768 0.068976677
## beak 0.8643005 0.2856416 -0.41095930 0.04915868 0.009455980
## mass 0.8155896 0.4887957 0.30684556 0.04115444 0.006674031
##
## $cor
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## thorax 0.9544052 -0.1068106 0.01326175 -0.27806416 0.014372683
## body 0.9684503 -0.2190542 0.03179998 0.06598983 -0.093559921
## wing 0.9173345 -0.3613220 0.06701864 0.13672768 0.068976677
## beak 0.8643005 0.2856416 -0.41095930 0.04915868 0.009455980
## mass 0.8155896 0.4887957 0.30684556 0.04115444 0.006674031
##
## $cos2
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## thorax 0.9108894 0.01140850 0.0001758741 0.077319675 2.065740e-04
## body 0.9378959 0.04798472 0.0010112387 0.004354657 8.753459e-03
## wing 0.8415027 0.13055361 0.0044914975 0.018694460 4.757782e-03
## beak 0.7470154 0.08159111 0.1688875494 0.002416575 8.941556e-05
## mass 0.6651863 0.23892123 0.0941541956 0.001693688 4.454269e-05
##
## $contrib
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## thorax 22.20333 2.234949 0.06544875 74.004952 1.4913182
## body 22.86163 9.400306 0.37631637 4.167971 63.1937785
## wing 20.51200 25.575720 1.67143927 17.893022 34.3478195
## beak 18.20883 15.983866 62.84881144 2.312976 0.6455171
## mass 16.21421 46.805159 35.03798418 1.621079 0.3215667
Morphology w/ Wing2Body
# No mass
d <- data_tested %>%
select(thorax_s, wing2body_s, beak_s) %>%
filter(!is.na(beak_s))
colnames(d) <- c("thorax", "wing2body", "beak")
MorphPCA = PCA_graphs(d, "(c) ")
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 1.7833365 59.44455 59.44455
## Dim.2 1.0244727 34.14909 93.59364
## Dim.3 0.1921908 6.40636 100.00000
## $coord
## Dim.1 Dim.2 Dim.3
## thorax 0.9282759 0.21904146 -0.30054053
## wing2body -0.1423525 0.98594791 0.08742133
## beak 0.9494082 -0.06633477 0.30695877
##
## $cor
## Dim.1 Dim.2 Dim.3
## thorax 0.9282759 0.21904146 -0.30054053
## wing2body -0.1423525 0.98594791 0.08742133
## beak 0.9494082 -0.06633477 0.30695877
##
## $cos2
## Dim.1 Dim.2 Dim.3
## thorax 0.86169623 0.047979161 0.090324609
## wing2body 0.02026423 0.972093280 0.007642489
## beak 0.90137601 0.004400302 0.094223687
##
## $contrib
## Dim.1 Dim.2 Dim.3
## thorax 48.31933 4.6833028 46.997367
## wing2body 1.13631 94.8871785 3.976512
## beak 50.54436 0.4295187 49.026121
# Mass
d <- data_tested %>%
select(thorax_s, wing2body_s, beak_s, mass_s) %>%
filter(!is.na(mass_s))
colnames(d) <- c("thorax", "wing2body", "beak", "mass")
MorphPCA = PCA_graphs(d, "(d) ")
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.4968987 62.422467 62.42247
## Dim.2 1.0313420 25.783551 88.20602
## Dim.3 0.2854844 7.137109 95.34313
## Dim.4 0.1862749 4.656873 100.00000
## $coord
## Dim.1 Dim.2 Dim.3 Dim.4
## thorax 0.9003120 0.270103200 -0.10531713 -0.32463972
## wing2body -0.1826277 0.976798037 0.05512816 0.09733239
## beak 0.9193385 -0.007102094 -0.30132658 0.25292027
## mass 0.8987775 -0.064818606 0.42491865 0.08626527
##
## $cor
## Dim.1 Dim.2 Dim.3 Dim.4
## thorax 0.9003120 0.270103200 -0.10531713 -0.32463972
## wing2body -0.1826277 0.976798037 0.05512816 0.09733239
## beak 0.9193385 -0.007102094 -0.30132658 0.25292027
## mass 0.8987775 -0.064818606 0.42491865 0.08626527
##
## $cos2
## Dim.1 Dim.2 Dim.3 Dim.4
## thorax 0.81056162 7.295574e-02 0.011091699 0.105390947
## wing2body 0.03335289 9.541344e-01 0.003039114 0.009473594
## beak 0.84518319 5.043974e-05 0.090797710 0.063968664
## mass 0.80780099 4.201452e-03 0.180555857 0.007441697
##
## $contrib
## Dim.1 Dim.2 Dim.3 Dim.4
## thorax 32.462735 7.07386456 3.885221 56.578179
## wing2body 1.335773 92.51386761 1.064547 5.085813
## beak 33.849318 0.00489069 31.804791 34.341000
## mass 32.352174 0.40737714 63.245442 3.995008
