The vegan package for PCA
Nathan Brouwer
10/25/2020
Preliminaries
Download the vegan package
Only do this once, then comment out of the script. You may have already done this this for a previous assignment.
#install.packages("vegan")Load the libraries
library(ggplot2)
library(vegan)## Warning: package 'vegan' was built under R version 4.2.2## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-4Load the msleep package
The mammals dataset is a classic dataset in the MASS package. msleep is an updated version of the data in the ggplot2 package that includes more numeric data (e.g. hours of sleep) and categorical data (e.g. if a species is endangered).
#make sure you've loaded ggplot2 w/ library(ggplot2_
data(msleep)Subset numeric variables
msleep has a mixture of numeric and categorical variables. We want only the numeric data and a few key labels, which I’ll also have to set up a special way. Don’t worry about the details.
First, we will subset all the columns we may want:
msleep_num <- msleep[,c("sleep_total","sleep_rem",
"sleep_cycle","awake","brainwt" ,
"bodywt",
"vore","order")]All data should be continuous numeric that can take on decimal values, except the last 2
summary(msleep_num[,-c(7:8)])## sleep_total sleep_rem sleep_cycle awake
## Min. : 1.90 Min. :0.100 Min. :0.1167 Min. : 4.10
## 1st Qu.: 7.85 1st Qu.:0.900 1st Qu.:0.1833 1st Qu.:10.25
## Median :10.10 Median :1.500 Median :0.3333 Median :13.90
## Mean :10.43 Mean :1.875 Mean :0.4396 Mean :13.57
## 3rd Qu.:13.75 3rd Qu.:2.400 3rd Qu.:0.5792 3rd Qu.:16.15
## Max. :19.90 Max. :6.600 Max. :1.5000 Max. :22.10
## NA's :22 NA's :51
## brainwt bodywt
## Min. :0.00014 Min. : 0.005
## 1st Qu.:0.00290 1st Qu.: 0.174
## Median :0.01240 Median : 1.670
## Mean :0.28158 Mean : 166.136
## 3rd Qu.:0.12550 3rd Qu.: 41.750
## Max. :5.71200 Max. :6654.000
## NA's :27We can look at these data with a scatterplot matrix
# add pairs()
pairs(msleep_num[,-c(7:8)])
The data are in tibble format (tbl_df)
is(msleep_num)## [1] "tbl_df" "tbl" "data.frame" "list" "oldClass"
## [6] "vector"Convert to dataframe because I don’t like tibbles!
# add data.frame()
msleep_num <- data.frame(msleep_num)Remove missing values with na.omit()
#add na.omit()
msleep_num <-na.omit(msleep_num)
msleep_num## sleep_total sleep_rem sleep_cycle awake brainwt bodywt vore
## 4 14.9 2.3 0.1333333 9.1 0.00029 0.019 omni
## 5 4.0 0.7 0.6666667 20.0 0.42300 600.000 herbi
## 9 10.1 2.9 0.3333333 13.9 0.07000 14.000 carni
## 12 9.4 0.8 0.2166667 14.6 0.00550 0.728 herbi
## 14 12.5 1.5 0.1166667 11.5 0.00640 0.420 herbi
## 17 9.1 1.4 0.1500000 14.9 0.00014 0.005 omni
## 18 17.4 3.1 0.3833333 6.6 0.01080 3.500 carni
## 20 18.0 4.9 0.3333333 6.0 0.00630 1.700 omni
## 22 19.7 3.9 0.1166667 4.3 0.00030 0.023 insecti
## 23 2.9 0.6 1.0000000 21.1 0.65500 521.000 herbi
## 25 10.1 3.5 0.2833333 13.9 0.00350 0.770 omni
## 28 12.5 3.2 0.4166667 11.5 0.02560 3.300 carni
## 29 9.8 1.1 0.5500000 14.2 0.00500 0.200 omni
## 34 8.0 1.9 1.5000000 16.0 1.32000 62.000 omni
## 38 10.1 1.2 0.7500000 13.9 0.17900 6.800 omni
## 40 14.3 3.1 0.2000000 9.7 0.00100 0.120 herbi
## 42 12.5 1.4 0.1833333 11.5 0.00040 0.022 herbi
## 43 19.9 2.0 0.2000000 4.1 0.00025 0.010 insecti
## 48 8.4 0.9 0.4166667 15.6 0.01210 2.500 herbi
## 50 9.7 1.4 1.4166667 14.3 0.44000 52.200 omni
## 54 9.4 1.0 0.6666667 14.6 0.18000 25.235 omni
## 64 13.0 2.4 0.1833333 11.0 0.00190 0.320 herbi
## 67 8.4 2.1 0.1666667 15.6 0.00120 0.075 insecti
## 68 11.3 1.1 0.1500000 12.7 0.00118 0.148 herbi
## 71 13.8 3.4 0.2166667 10.2 0.00400 0.101 herbi
## 74 9.1 2.4 0.5000000 14.9 0.18000 86.250 omni
## 77 4.4 1.0 0.9000000 19.6 0.16900 207.501 herbi
## 79 8.9 2.6 0.2333333 15.1 0.00250 0.104 omni
## 83 9.8 2.4 0.3500000 14.2 0.05040 4.230 carni
## order
## 4 Soricomorpha
## 5 Artiodactyla
## 9 Carnivora
## 12 Rodentia
## 14 Rodentia
## 17 Soricomorpha
## 18 Cingulata
## 20 Didelphimorphia
## 22 Chiroptera
## 23 Perissodactyla
## 25 Erinaceomorpha
## 28 Carnivora
## 29 Primates
## 34 Primates
## 38 Primates
## 40 Rodentia
## 42 Rodentia
## 43 Chiroptera
## 48 Lagomorpha
## 50 Primates
## 54 Primates
## 64 Rodentia
## 67 Soricomorpha
## 68 Rodentia
## 71 Rodentia
## 74 Artiodactyla
## 77 Perissodactyla
## 79 Scandentia
## 83 CarnivoraPrincipal components analysis - base R
Principal component analysis is typically done using base R functions using the prcomp() function.
Run the PCA.
# add prcomp()
pca.out <- prcomp(msleep_num[,-c(7,8)], scale = TRUE)Plot the biplot with biplot()
# add biplot()
biplot(pca.out)
Principal components analysis - vegan
The base R PCA output isn’t very flexible. The R package vegan is has a function rda() which does PCA and has many more nice features.
For another example see https://rpubs.com/brouwern/veganpca
Run the PCA using vegan::rda()
rda.out <- vegan::rda(msleep_num[,-c(7,8)], scale = TRUE)Extract what are known as the PCA scores using the scores() function.
# add scores()
rda_scores <- scores(rda.out)vegan let’s you plot a scatterplot of the scores or a full biplot with the arrows (vectors). The displays the 2D scatterplot plot without the arrows. For more info on what this code does, see the RPubs document linked above.
# add biplot()
biplot(rda.out,
display = "sites")
vegan has some nice tools for grouping things.
In this dataset I don’t expect there to be any interesting groups, but I’ll check anyway.
PCA is an exploratory method. First I’ll see if there are any groupings based on diet (“vore”). Not really.
#make the plot with biplot()
# add biplot()
biplot(rda.out, display = "sites")
#add "hulls" with the ordihull()function
# add ordihull()
ordihull(rda.out,
group = msleep_num$vore,
col = 1:7,
lty = 1:7,
lwd = c(3,6))
Task
Copy the previous code chunk and change the group so that the taxonomic order is plotted (column “order”). Upload the image to this assignment. Consider if there are any meaningful groups.
## add your code here
#make the plot with biplot()
# add biplot()
biplot(rda.out, display = "sites")
#add "hulls" with the ordihull()function
# add ordihull()
ordihull(rda.out,
group = msleep_num$order,
col = 1:7,
lty = 1:7,
lwd = c(3,6))