Ex5_solved

Solving Problems with R 205-2-7032 Lab assignment no.5

package vegan

The data set that we will be using is a real data set, contributed by Zehava Siegal, who works as an ecologist for the Nature and Parks Authority. The data contain numbers of annual plants from 40 species, collected near Kibutz Sede Boqer. In each plot there were three treatments: HA1-plowing the soil superficially; HA2-plowing the soil deeply; control-the soil was left untouched. This design was replicated in five plots. There are 5 areas, and in each one there were three subplots, with different treatments (HA1, HA2 and Control). So in total there were 15 sites, and in each one of them Zehava recorded the quantity of annual plants she found. The goal of the study was to determine which of the treatments results in a richer annual plant community. There could be many definitions of “rich plant community”, and we will explore several potential measures.

0. fill in your name

My_Name <- 'Gal Revivo'

1. Use ‘sp.csv’ data file and transform the data into the structure for the package, and call it VegSp. (20 points)

Have a look at the data file. It is in long format. Last week we learnt how to convert wide to long format. This week, we need to do the opposit, because the package we will use (vegan) for the analysis, requires data in wide format.

Structure of the data for the analysis should be: Each row is a location. Do not include location names in the table for analysis. sp1 sp2 sp3 sp4 sp5 sp6 sp7 sp8 6 5 1 3 12 50 0 0 12 21 5 25 2 17 9 0 The data frame sould not contain NA, please replace them to 0 (zero)

tip - try to look for function “dcast”" from package “reshape2”: https://seananderson.ca/2013/10/19/reshape/ .

# import data set and switch format
library(readr)
library(reshape2) #call in the different libraries that i will use for this exercise
library(vegan)
library(readr)
library(ggplot2)  #################################################################

Sp <- read_csv(file = 'C:/University/solve problems in R/Sp.csv')
Sp <- as.data.frame(Sp)#call in the Data set
View(Sp)
colnames(Sp)[2] <- "Area.num" #change the name of the column to be more easy to use
colnames(Sp)[5] <- "Obs" # dito

Vegsp <- dcast(Sp, Area.num + Treatment ~ Species, value.var = "Obs") # use dcast function from 'vegan' package to switch to wide format  

Vegsp[is.na(Vegsp)] <- 0 # switch from N/A to 0 cells

2. The first measure we will look at is species richness (a simple count of number of species). Plot the number of species that were observed for each place in bar graph, use ggplot. rotate the names of locations that will appear in the x axis to be vertical, using “theme(axis.text.x=element_text(angle = -90, hjust = 0))” 20 points

#calculate richness per location and per treatment

Vegsp$richness<-rowSums(Vegsp[,-1:-2] != 0) # add new column labeled 'richness' to sum every row (richness for location and treatment) without taking into account the first two columns

p<- ggplot(Vegsp,aes(Area.num, richness, fill= Treatment)) + # plot graph for the richness of different treatment and areas
  geom_bar(stat='identity',position = "dodge") +
  theme(axis.text.x=element_text(angle = -90, hjust = 0))

p+scale_fill_manual(values=c("#DA70D6", "#770737", "#673147")) #custom color pellet

3)species diversity with the package vegan- calculate three measures of diversity: Simpson, Shannon and Fisher’s alpha. Combine them with the species richness you calculated before in one data frame, and insert the place names. The different measures names should appear as the column names. 20 points

richness<-rowSums(Vegsp[,-1:-2] != 0) # two ways to calculate the richness - with rows summary and with 'specnumber' function

richness2<-specnumber(Vegsp2)

Vegsp2 <-Vegsp[,-1:-2] # new data frame without the first two columns

# different alpha diversities:
#simpson
Simpson <- diversity(Vegsp2,index = "simpson")
#Shannon
Shannon <- diversity(Vegsp2)
#Fisher
Fisher <- fisher.alpha(Vegsp2)

DiversityTable <- cbind(Vegsp[,1:2], Simpson,Shanon,Fisher, richness2) # take the four vectors from above and add them together 

colnames(DiversityTable) <- c("Area.num", "Treatment","Simpson", "Shannon", "Fisher", "Richness") # adding names for the columns

4. Based on last weeks’s lesson, make a box plot that summarizes the Shanon median for each treatment. Use diffrent colors for each treatment. Make sure that the axis have a meaning and add a title to the graph (10 points)

#boxplot the different treatment by the shannon median values (the line inside each of the boxplots)
ggplot (as.data.frame(DiversityTable), aes(x = Treatment, y = Shanon, fill = Treatment)) +  geom_boxplot(alpha=0.3) +
theme(legend.position="none")+ggtitle("Shannon diversity Median")

5. Calculate and compare in a graph the accumulation of species richness for treatment HA1 and HA2. Use the default method for the order of plots accumulation of species. Use function ‘specaccum’ from Vegan. USe “’”Package ‘vegan’" file to know how to use this function. You can find the file in the moodle.

10 points

H1 <- specaccum(Vegsp[Vegsp$Treatment=="HA1",-c(1:2)]) #calculate the accumulation of richness in HA1 treatment
H2 <- specaccum(Vegsp[Vegsp$Treatment=="HA2",-c(1:2)])#calculate the accumulation of richness in HA2 treatment

#plot graph on both treatments for comparison 
plot(H1, col="orange") 
plot(H2, col="purple", add = T) 

6. Cluster analysis - build single linkage cluster, do not forget to plot it. Use ‘vegan tutor for Ex5’, section 6. The tutor file is in moodle. 10 points

#Cluster analysis - single linkage cluster

distance<- vegdist(Vegsp[,-c(1:2,42)]) # use vegdist, to preform distance matrix from our original wide data on the 15 samples
clus1 <- hclust(distance, "single") #use single linkage cluster, this dendrogram is the one that the qusetion asked
clus2 <- hclust(distance, "complete") # the complete linkage cluster, just for me to see the differences
clus3 <- hclust(distance, "average") # avrage linkage cluster just for me to see the differences
#ploting the cluster
plot(clus1) # plot the cluster dendrograms
plot(clus2)
plot(clus3)

7. Display an interpretation of classes - use the cluster you built before and classify the places in 3 seprate groups. 10 points

plot(clus1) #plot the cluster from question 6

rect.hclust(clus1, 3) # classification into groups of threes
group3 <- cutree(clus1, 3)

Homework: 1. Complete all the tasks that you did not do in class. 2. Go over your chunks, after you made sure it runs, and make sure you wrote a comment for each line - a title that will explain what the line does. Try to use the terminology used in class whenever possible. 10% of every question will be given for comments. 3. Submit the RMD file and the HTML file in moodle. Both file names need to include your first and last name. If I had to submit the files I would name them: ex5_LiranSagi

                        Good Luck!