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.package("vegan")

Load the libraries

library(ggplot2)
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4

Load 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   :27

We 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)

Principal components analysis - base R

Principal component analysis is typically done using base R functions using the prcomp() function.

Run the PCA.

# add prcomp()    prcomp() on a dataframe with NAs removed = a biplot
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()

# vegan package as rda() function -> better plotting for redundancy analysis
    # built specialized for ecological analysis
    # special because used WITH predictors -> if no predictors added, same as PCA
    # display = "sites" (corresponds to rows)     display = "species" (corresponds to cols)
    # type = "text" (plots row/col names as labels)     type = "points" (plots points)

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 with extracted PCAs + display="sites" creates full biplot WITHOUT vectors
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  -> uses polygons or "hulls" to group points together
          # ordiellipse() -> uses ellipses to group centroids (SD)
          # display onto respective biplot by running function in the same code chunk as biplot()
# 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

biplot(rda.out, display = "sites")

ordihull(rda.out,
         group = msleep_num$order,
         col = 1:7,
         lty = 1:7,
         lwd = c(3,6))