ENV567_ProblemSet_2_MGalindo
Maureen_Galindo
3/26/2020
PRACTICE PROBLEMS – Time series and spatial analysis
- Write a function to calculate the 3-point moving average for an input vector. The formula is: y_i^’= (y_(i-1)+y_i+y_(i+1))/3
moveave3pt <- function (x) {
y<-numeric(length(x)-2)
for (i in 2:(length(x)-1)) {
y[i]<-(x[i-1]+x[i]+x[i+1])/3
}
y
}- Read the temp.txt file. The data correspond to monthly average temperatures.
library(TTR)
temp<-read.table("temp.txt", header=T)
temptimeseries <- ts(temp$temps, frequency=12)2a) Plot the data
#plot "temp" as timeseries
plot(temptimeseries, ylab="Temp Over Months")
2b) Calculate the 5-point moving average and plot it together with the time series
plot(ts(temp$temps, frequency=1)) #monthly
tempmov<-SMA(ts(temp$temps, frequency=1), n = 5)
lines(tempmov, col="red", lwd=2)
2c) Decompose the time series into seasonal, trend and residual error components assuming the seasonal trend is yearly.
#use TTR function to decompose time series
tempdec<-decompose(temptimeseries)
plot(tempdec)
2d) Generate a temporal correlogram to assess the autocorrelation of the time series
#correlogram
acf(temptimeseries) 
2e) Generate a new correlogram but removing the trend and seasonal variation
acf(tempdec$random[!is.na(tempdec$random)])
2f) Generate a partial correlogram plot
pacf(temptimeseries) 
2g) Find the best ARMA model using the forecast package
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zootempfit<-auto.arima(temptimeseries)
summary(tempfit)## Series: temptimeseries
## ARIMA(0,0,1)(2,1,0)[12]
##
## Coefficients:
## ma1 sar1 sar2
## 0.2456 -0.5212 -0.3233
## s.e. 0.0774 0.0823 0.0827
##
## sigma^2 estimated as 3.772: log likelihood=-300.76
## AIC=609.52 AICc=609.81 BIC=621.4
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.2384692 1.846384 1.498262 0.2682136 11.90101 0.8327783
## ACF1
## Training set 0.0054130472h) Estimate future values using the previous ARMA model and plot the results
tempfcast<-forecast(tempfit)
plot(tempfcast)
3. Read the bac_dust.txt dataset. The data corresponds to proteobacterial abundance and bacterial species in dust samples collected around the globe.
bacdust<-read.table("bac_dust.txt",header=T)3a) Make a map of the sampling points using the R libraries maps and ggplot2. Color the points based on the continent they were collected.
library(maps)
library(ggplot2)
attach(bacdust)
world <- map_data("world")
ggplot(world)+
geom_polygon(aes(x = long, y = lat, group = group), fill = "lightgray", colour = "black", size = 0.1)+
geom_point(data=bacdust, aes(x=Longitude, y=Latitude, color=Continent), size=2)+
scale_color_manual(values = c("SouthAmerica" = 'green', "NorthAmerica" = 'purple', "Africa" = 'goldenrod1', "Oceania" = 'cornflowerblue', "Europe" = 'red', "Asia" = 'orange'))+
theme_bw(base_size = 15)+
theme(legend.title = element_blank(), axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank())
theme_bw()## List of 93
## $ line :List of 6
## ..$ colour : chr "black"
## ..$ size : num 0.5
## ..$ linetype : num 1
## ..$ lineend : chr "butt"
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ rect :List of 5
## ..$ fill : chr "white"
## ..$ colour : chr "black"
## ..$ size : num 0.5
## ..$ linetype : num 1
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ text :List of 11
## ..$ family : chr ""
## ..$ face : chr "plain"
## ..$ colour : chr "black"
## ..$ size : num 11
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : num 0
## ..$ lineheight : num 0.9
## ..$ margin : 'margin' num [1:4] 0pt 0pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ title : NULL
## $ aspect.ratio : NULL
## $ axis.title : NULL
## $ axis.title.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.75pt 0pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 2.75pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.bottom : NULL
## $ axis.title.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 2.75pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.y.left : NULL
## $ axis.title.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 0pt 2.75pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey30"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.2pt 0pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 2.2pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.bottom : NULL
## $ axis.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 1
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 2.2pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.y.left : NULL
## $ axis.text.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 0pt 2.2pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.ticks :List of 6
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ axis.ticks.x : NULL
## $ axis.ticks.x.top : NULL
## $ axis.ticks.x.bottom : NULL
## $ axis.ticks.y : NULL
## $ axis.ticks.y.left : NULL
## $ axis.ticks.y.right : NULL
## $ axis.ticks.length : 'unit' num 2.75pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ axis.ticks.length.x : NULL
## $ axis.ticks.length.x.top : NULL
## $ axis.ticks.length.x.bottom: NULL
## $ axis.ticks.length.y : NULL
## $ axis.ticks.length.y.left : NULL
## $ axis.ticks.length.y.right : NULL
## $ axis.line : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.line.x : NULL
## $ axis.line.x.top : NULL
## $ axis.line.x.bottom : NULL
## $ axis.line.y : NULL
## $ axis.line.y.left : NULL
## $ axis.line.y.right : NULL
## $ legend.background :List of 5
## ..$ fill : NULL
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.margin : 'margin' num [1:4] 5.5pt 5.5pt 5.5pt 5.5pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ legend.spacing : 'unit' num 11pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ legend.spacing.x : NULL
## $ legend.spacing.y : NULL
## $ legend.key :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.key.size : 'unit' num 1.2lines
## ..- attr(*, "valid.unit")= int 3
## ..- attr(*, "unit")= chr "lines"
## $ legend.key.height : NULL
## $ legend.key.width : NULL
## $ legend.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.text.align : NULL
## $ legend.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.title.align : NULL
## $ legend.position : chr "right"
## $ legend.direction : NULL
## $ legend.justification : chr "center"
## $ legend.box : NULL
## $ legend.box.just : NULL
## $ legend.box.margin : 'margin' num [1:4] 0cm 0cm 0cm 0cm
## ..- attr(*, "valid.unit")= int 1
## ..- attr(*, "unit")= chr "cm"
## $ legend.box.background : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ legend.box.spacing : 'unit' num 11pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ panel.background :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.border :List of 5
## ..$ fill : logi NA
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.spacing : 'unit' num 5.5pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ panel.spacing.x : NULL
## $ panel.spacing.y : NULL
## $ panel.grid :List of 6
## ..$ colour : chr "grey92"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major : NULL
## $ panel.grid.minor :List of 6
## ..$ colour : NULL
## ..$ size : 'rel' num 0.5
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major.x : NULL
## $ panel.grid.major.y : NULL
## $ panel.grid.minor.x : NULL
## $ panel.grid.minor.y : NULL
## $ panel.ontop : logi FALSE
## $ plot.background :List of 5
## ..$ fill : NULL
## ..$ colour : chr "white"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ plot.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 5.5pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.title.position : chr "panel"
## $ plot.subtitle :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 5.5pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : num 1
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 5.5pt 0pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption.position : chr "panel"
## $ plot.tag :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.tag.position : chr "topleft"
## $ plot.margin : 'margin' num [1:4] 5.5pt 5.5pt 5.5pt 5.5pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ strip.background :List of 5
## ..$ fill : chr "grey85"
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ strip.background.x : NULL
## $ strip.background.y : NULL
## $ strip.placement : chr "inside"
## $ strip.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey10"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 4.4pt 4.4pt 4.4pt 4.4pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.x : NULL
## $ strip.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.switch.pad.grid : 'unit' num 2.75pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ strip.switch.pad.wrap : 'unit' num 2.75pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ strip.text.y.left :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## - attr(*, "class")= chr [1:2] "theme" "gg"
## - attr(*, "complete")= logi TRUE
## - attr(*, "validate")= logi TRUE3b) Make a map of the sampling points using the R libraries maps and ggplot2. Size the points based on bacterial richness.
world_map <- map_data("world")
ggplot(world_map)+
geom_polygon(aes(x = long, y = lat, group = group), fill = "lightgray", colour = "black", size = 0.1)+
geom_point(data=bacdust, aes(x=Longitude, y=Latitude, size=Richness, color=Continent), alpha=0.6)
theme_bw(base_size = 15)+
theme(legend.title = element_blank(), axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank())## List of 93
## $ line :List of 6
## ..$ colour : chr "black"
## ..$ size : num 0.682
## ..$ linetype : num 1
## ..$ lineend : chr "butt"
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ rect :List of 5
## ..$ fill : chr "white"
## ..$ colour : chr "black"
## ..$ size : num 0.682
## ..$ linetype : num 1
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ text :List of 11
## ..$ family : chr ""
## ..$ face : chr "plain"
## ..$ colour : chr "black"
## ..$ size : num 15
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : num 0
## ..$ lineheight : num 0.9
## ..$ margin : 'margin' num [1:4] 0pt 0pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ title : NULL
## $ aspect.ratio : NULL
## $ axis.title : NULL
## $ axis.title.x : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.title.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 3.75pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.bottom : NULL
## $ axis.title.y : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.title.y.left : NULL
## $ axis.title.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 0pt 3.75pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey30"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.text.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 3pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.bottom : NULL
## $ axis.text.y : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.text.y.left : NULL
## $ axis.text.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 0pt 3pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.ticks :List of 6
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ axis.ticks.x : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.ticks.x.top : NULL
## $ axis.ticks.x.bottom : NULL
## $ axis.ticks.y : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.ticks.y.left : NULL
## $ axis.ticks.y.right : NULL
## $ axis.ticks.length : 'unit' num 3.75pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ axis.ticks.length.x : NULL
## $ axis.ticks.length.x.top : NULL
## $ axis.ticks.length.x.bottom: NULL
## $ axis.ticks.length.y : NULL
## $ axis.ticks.length.y.left : NULL
## $ axis.ticks.length.y.right : NULL
## $ axis.line : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.line.x : NULL
## $ axis.line.x.top : NULL
## $ axis.line.x.bottom : NULL
## $ axis.line.y : NULL
## $ axis.line.y.left : NULL
## $ axis.line.y.right : NULL
## $ legend.background :List of 5
## ..$ fill : NULL
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.margin : 'margin' num [1:4] 7.5pt 7.5pt 7.5pt 7.5pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ legend.spacing : 'unit' num 15pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ legend.spacing.x : NULL
## $ legend.spacing.y : NULL
## $ legend.key :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.key.size : 'unit' num 1.2lines
## ..- attr(*, "valid.unit")= int 3
## ..- attr(*, "unit")= chr "lines"
## $ legend.key.height : NULL
## $ legend.key.width : NULL
## $ legend.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.text.align : NULL
## $ legend.title : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ legend.title.align : NULL
## $ legend.position : chr "right"
## $ legend.direction : NULL
## $ legend.justification : chr "center"
## $ legend.box : NULL
## $ legend.box.just : NULL
## $ legend.box.margin : 'margin' num [1:4] 0cm 0cm 0cm 0cm
## ..- attr(*, "valid.unit")= int 1
## ..- attr(*, "unit")= chr "cm"
## $ legend.box.background : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ legend.box.spacing : 'unit' num 15pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ panel.background :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.border :List of 5
## ..$ fill : logi NA
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.spacing : 'unit' num 7.5pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ panel.spacing.x : NULL
## $ panel.spacing.y : NULL
## $ panel.grid :List of 6
## ..$ colour : chr "grey92"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major : NULL
## $ panel.grid.minor :List of 6
## ..$ colour : NULL
## ..$ size : 'rel' num 0.5
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major.x : NULL
## $ panel.grid.major.y : NULL
## $ panel.grid.minor.x : NULL
## $ panel.grid.minor.y : NULL
## $ panel.ontop : logi FALSE
## $ plot.background :List of 5
## ..$ fill : NULL
## ..$ colour : chr "white"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ plot.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 7.5pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.title.position : chr "panel"
## $ plot.subtitle :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0pt 0pt 7.5pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : num 1
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 7.5pt 0pt 0pt 0pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption.position : chr "panel"
## $ plot.tag :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.tag.position : chr "topleft"
## $ plot.margin : 'margin' num [1:4] 7.5pt 7.5pt 7.5pt 7.5pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ strip.background :List of 5
## ..$ fill : chr "grey85"
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ strip.background.x : NULL
## $ strip.background.y : NULL
## $ strip.placement : chr "inside"
## $ strip.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey10"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 6pt 6pt 6pt 6pt
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.x : NULL
## $ strip.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.switch.pad.grid : 'unit' num 3.75pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ strip.switch.pad.wrap : 'unit' num 3.75pt
## ..- attr(*, "valid.unit")= int 8
## ..- attr(*, "unit")= chr "pt"
## $ strip.text.y.left :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## - attr(*, "class")= chr [1:2] "theme" "gg"
## - attr(*, "complete")= logi TRUE
## - attr(*, "validate")= logi TRUE3c) Calculate Moran’s I autocorrelation for the variable Richness and Proteobacteria
library(ape)
bacdustdist <- as.matrix(dist(cbind(Latitude, Longitude)))
bacdustdistinv <- 1/bacdustdist
diag(bacdustdistinv) <- 0
Moran.I(Proteobacteria, bacdustdistinv) ## $observed
## [1] 0.216501
##
## $expected
## [1] -0.003225806
##
## $sd
## [1] 0.03858986
##
## $p.value
## [1] 1.241692e-08There is a significant autocorrelation
3d) Make a simple linear regression of Richness as a function of Proteobacteria. Are the residuals spatially autocorrelated?
bacdustlm <- lm(Proteobacteria ~ Richness)
summary(bacdustlm)##
## Call:
## lm(formula = Proteobacteria ~ Richness)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36741 -0.08251 0.00027 0.07813 0.54164
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.693e-01 1.489e-02 24.797 <2e-16 ***
## Richness -2.165e-05 2.960e-05 -0.731 0.465
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1405 on 309 degrees of freedom
## Multiple R-squared: 0.001728, Adjusted R-squared: -0.001502
## F-statistic: 0.5349 on 1 and 309 DF, p-value: 0.46513e) Perform a generalized linear mixed model (GLMM) using the coordinates as random effect and using the Matérn covariance function (R package spaMM).
library(spaMM)## Registered S3 methods overwritten by 'registry':
## method from
## print.registry_field proxy
## print.registry_entry proxy## spaMM (Rousset & Ferdy, 2014, version 3.2.0) is loaded.
## Type 'help(spaMM)' for a short introduction,
## 'news(package='spaMM')' for news,
## and 'citation(spaMM)' for proper citation.fitme(Proteobacteria ~ Richness+Matern(1|Latitude+Longitude),data=bacdust)## formula: Proteobacteria ~ Richness + Matern(1 | Latitude + Longitude)
## ML: Estimation of corrPars, lambda and phi by ML.
## Estimation of fixed effects by ML.
## Estimation of lambda and phi by 'outer' ML, maximizing p_v.
## Family: gaussian ( link = identity )
## ------------ Fixed effects (beta) ------------
## Estimate Cond. SE t-value
## (Intercept) 3.609e-01 2.040e-02 17.6897
## Richness 1.482e-05 2.821e-05 0.5252
## --------------- Random effects ---------------
## Family: gaussian ( link = identity )
## --- Correlation parameters:
## 1.nu 1.rho
## 5.145130 1.072765
## --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian;
## Latitude . : 0.005976
## # of obs: 311; # of groups: Latitude ., 311
## ------------- Residual variance -------------
## phi estimate was 0.0143184
## ------------- Likelihood values -------------
## logLik
## p_v(h) (marginal L): 193.72593f) Report the nu and the rho values for this spatial regression model.
Nu = 5.145 rho = 1.0727
3g) What can you say about the estimated correlation of this particular spatial regression model across increasing distances between pairs of locations.
dd<-dist(bacdust[,c("Longitude", "Latitude")])
mm<-MaternCorr(dd, nu=5.1456, rho=1.0727)
plot(as.numeric(dd), as.numeric(mm), xlab="Distance between pairs of locations", ylab="Estimated correlation")
at a distance of about 10 the correlation is pretty small 4. Download the map of Spain. 4a) Make a basic map of Spain
spain<-map_data("world", region=c("Spain"))
ggplot(spain)+
geom_polygon(aes(x = long, y = lat, group=group))+
theme_bw()
4b) Get information on the population of Spanish cities and make a map of Spain including the locations and population of Spanish cities.
cities<-get('world.cities')
citiesSpain<-cities[cities$country.etc == 'Spain',]
ggplot(spain)+
geom_polygon(aes(x = long, y = lat, group=group, ), fill='lightgray', size=0.1)+
geom_point(data=citiesSpain, aes(x=long, y=lat, size=pop, color=pop))+
theme_bw()
4c) Visit Spain and enjoy >It is on the bucket list. I hear the beaches are amazing.
PRACTICE PROBLEMS – Diversity analysis
- Load the dune_bio.txt dataset. Species should be in columns and sites in rows.
library(vegan)## Loading required package: permute##
## Attaching package: 'permute'## The following object is masked from 'package:spaMM':
##
## how## Loading required package: lattice## This is vegan 2.5-6dunebio<-read.table('dune_bio.txt', header = T, row.names=1)1a) Calculate the total number of individuals of all species
sum(dunebio)## [1] 6851b) Calculate the total number of individuals for each species
colSums(dunebio)## Belper Empnig Junbuf Junart Airpra Elepal Rumace Viclat Brarut Ranfla Cirarv
## 13 2 13 18 5 25 18 4 49 14 2
## Hyprad Leoaut Potpal Poapra Calcus Tripra Trirep Antodo Salrep Achmil Poatri
## 9 54 4 48 10 9 47 21 11 16 63
## Chealb Elyrep Sagpro Plalan Agrsto Lolper Alogen Brohor
## 1 26 20 26 48 58 36 151c) Calculate the average number of individuals for each species
colMeans(dunebio)## Belper Empnig Junbuf Junart Airpra Elepal Rumace Viclat Brarut Ranfla Cirarv
## 0.65 0.10 0.65 0.90 0.25 1.25 0.90 0.20 2.45 0.70 0.10
## Hyprad Leoaut Potpal Poapra Calcus Tripra Trirep Antodo Salrep Achmil Poatri
## 0.45 2.70 0.20 2.40 0.50 0.45 2.35 1.05 0.55 0.80 3.15
## Chealb Elyrep Sagpro Plalan Agrsto Lolper Alogen Brohor
## 0.05 1.30 1.00 1.30 2.40 2.90 1.80 0.751d) Calculate the total number of individuals for each site
rowSums(dunebio)## 2 13 4 16 6 1 8 5 17 15 10 11 9 18 3 20 14 19 12 7
## 42 33 45 33 48 18 40 43 15 23 43 32 42 27 40 31 24 31 35 401e) Calculate the average number of individuals for each site
rowMeans(dunebio)## 2 13 4 16 6 1 8 5
## 1.4000000 1.1000000 1.5000000 1.1000000 1.6000000 0.6000000 1.3333333 1.4333333
## 17 15 10 11 9 18 3 20
## 0.5000000 0.7666667 1.4333333 1.0666667 1.4000000 0.9000000 1.3333333 1.0333333
## 14 19 12 7
## 0.8000000 1.0333333 1.1666667 1.33333331f) Write a function to report the median number of individuals for each species and the median number of individuals for each site
mediandunes<-function(d) {
sitemed<-apply(d, 1, median)
specmed<-apply(d, 2, median)
print(sitemed)
print(specmed)
}
mediandunes(dunebio)## 2 13 4 16 6 1 8 5 17 15 10 11 9 18 3 20 14 19 12 7
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## Belper Empnig Junbuf Junart Airpra Elepal Rumace Viclat Brarut Ranfla Cirarv
## 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0 0.0
## Hyprad Leoaut Potpal Poapra Calcus Tripra Trirep Antodo Salrep Achmil Poatri
## 0.0 2.0 0.0 3.0 0.0 0.0 2.0 0.0 0.0 0.0 4.0
## Chealb Elyrep Sagpro Plalan Agrsto Lolper Alogen Brohor
## 0.0 0.0 0.0 0.0 1.5 2.0 0.0 0.0specnumber(dunebio)## 2 13 4 16 6 1 8 5 17 15 10 11 9 18 3 20 14 19 12 7
## 10 10 13 8 11 5 12 14 7 8 12 9 13 9 10 8 7 9 9 131g) Use the vegan function decostand() to transform the dataset to relative abundances. How would you confirm the transformation worked?
dunebiodf<-decostand(dunebio, method = "total")
summary(dunebiodf)## Belper Empnig Junbuf Junart
## Min. :0.00000 Min. :0.000000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.000000 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.00000 Median :0.000000 Median :0.00000 Median :0.00000
## Mean :0.01665 Mean :0.003226 Mean :0.01752 Mean :0.02728
## 3rd Qu.:0.04496 3rd Qu.:0.000000 3rd Qu.:0.00000 3rd Qu.:0.02273
## Max. :0.07407 Max. :0.064516 Max. :0.11429 Max. :0.13043
## Airpra Elepal Rumace Viclat
## Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.00000 Median :0.00000 Median :0.00000 Median :0.00000
## Mean :0.01151 Mean :0.04278 Mean :0.02105 Mean :0.00614
## 3rd Qu.:0.00000 3rd Qu.:0.02500 3rd Qu.:0.01190 3rd Qu.:0.00000
## Max. :0.13333 Max. :0.24242 Max. :0.12500 Max. :0.06250
## Brarut Ranfla Cirarv Hyprad
## Min. :0.00000 Min. :0.00000 Min. :0.000000 Min. :0.00000
## 1st Qu.:0.03333 1st Qu.:0.00000 1st Qu.:0.000000 1st Qu.:0.00000
## Median :0.05000 Median :0.00000 Median :0.000000 Median :0.00000
## Mean :0.07213 Mean :0.02353 Mean :0.002222 Mean :0.01786
## 3rd Qu.:0.12216 3rd Qu.:0.05265 3rd Qu.:0.000000 3rd Qu.:0.00000
## Max. :0.22222 Max. :0.12903 Max. :0.044444 Max. :0.16129
## Leoaut Potpal Poapra Calcus
## Min. :0.00000 Min. :0.000000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.05536 1st Qu.:0.000000 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.06977 Median :0.000000 Median :0.07778 Median :0.00000
## Mean :0.08170 Mean :0.008514 Mean :0.06960 Mean :0.01772
## 3rd Qu.:0.09498 3rd Qu.:0.000000 3rd Qu.:0.10000 3rd Qu.:0.00000
## Max. :0.19355 Max. :0.086957 Max. :0.22222 Max. :0.16667
## Tripra Trirep Antodo Salrep
## Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.03816 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.00000 Median :0.05530 Median :0.00000 Median :0.00000
## Mean :0.01003 Mean :0.06625 Mean :0.03471 Mean :0.01846
## 3rd Qu.:0.00000 3rd Qu.:0.08772 3rd Qu.:0.05312 3rd Qu.:0.00000
## Max. :0.10417 Max. :0.25000 Max. :0.26667 Max. :0.16129
## Achmil Poatri Chealb Elyrep
## Min. :0.00000 Min. :0.00000 Min. :0.000000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.000000 1st Qu.:0.00000
## Median :0.00000 Median :0.09651 Median :0.000000 Median :0.00000
## Mean :0.02458 Mean :0.08232 Mean :0.001515 Mean :0.03711
## 3rd Qu.:0.04738 3rd Qu.:0.12054 3rd Qu.:0.000000 3rd Qu.:0.08992
## Max. :0.13333 Max. :0.27273 Max. :0.030303 Max. :0.22222
## Sagpro Plalan Agrsto Lolper
## Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.00000 Median :0.00000 Median :0.03571 Median :0.06085
## Mean :0.02714 Mean :0.03767 Mean :0.07145 Mean :0.08353
## 3rd Qu.:0.05265 3rd Qu.:0.09635 3rd Qu.:0.15396 3rd Qu.:0.12863
## Max. :0.11429 Max. :0.13333 Max. :0.21212 Max. :0.38889
## Alogen Brohor
## Min. :0.00000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.00000 Median :0.00000
## Mean :0.04824 Mean :0.01757
## 3rd Qu.:0.08387 3rd Qu.:0.01163
## Max. :0.22857 Max. :0.09524confirm it works by summing all rows. should equal 1
rowSums(dunebiodf)## 2 13 4 16 6 1 8 5 17 15 10 11 9 18 3 20 14 19 12 7
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 11g2) Use the vegan function decostand() to standardize the dataset into the range 0 to 1
decostand(dunebio,"range")## Belper Empnig Junbuf Junart Airpra Elepal Rumace Viclat Brarut
## 2 1.0000000 0 0.00 0.00 0.0000000 0.000 0.0000000 0.0 0.0000000
## 13 0.0000000 0 0.75 0.00 0.0000000 0.000 0.0000000 0.0 0.0000000
## 4 0.6666667 0 0.00 0.00 0.0000000 0.000 0.0000000 0.0 0.3333333
## 16 0.0000000 0 0.00 0.75 0.0000000 1.000 0.0000000 0.0 0.6666667
## 6 0.0000000 0 0.00 0.00 0.0000000 0.000 1.0000000 0.0 1.0000000
## 1 0.0000000 0 0.00 0.00 0.0000000 0.000 0.0000000 0.0 0.0000000
## 8 0.0000000 0 0.00 1.00 0.0000000 0.500 0.0000000 0.0 0.3333333
## 5 0.6666667 0 0.00 0.00 0.0000000 0.000 0.8333333 0.0 0.3333333
## 17 0.0000000 0 0.00 0.00 0.6666667 0.000 0.0000000 0.0 0.0000000
## 15 0.0000000 0 0.00 0.75 0.0000000 0.625 0.0000000 0.0 0.6666667
## 10 0.6666667 0 0.00 0.00 0.0000000 0.000 0.0000000 0.5 0.3333333
## 11 0.0000000 0 0.00 0.00 0.0000000 0.000 0.0000000 1.0 0.6666667
## 9 0.0000000 0 1.00 1.00 0.0000000 0.000 0.3333333 0.0 0.3333333
## 18 0.6666667 0 0.00 0.00 0.0000000 0.000 0.0000000 0.5 1.0000000
## 3 0.6666667 0 0.00 0.00 0.0000000 0.000 0.0000000 0.0 0.3333333
## 20 0.0000000 0 0.00 1.00 0.0000000 0.500 0.0000000 0.0 0.6666667
## 14 0.0000000 0 0.00 0.00 0.0000000 0.500 0.0000000 0.0 0.0000000
## 19 0.0000000 1 0.00 0.00 1.0000000 0.000 0.0000000 0.0 0.5000000
## 12 0.0000000 0 1.00 0.00 0.0000000 0.000 0.3333333 0.0 0.6666667
## 7 0.0000000 0 0.50 0.00 0.0000000 0.000 0.5000000 0.0 0.3333333
## Ranfla Cirarv Hyprad Leoaut Potpal Poapra Calcus Tripra Trirep Antodo
## 2 0.0 0 0.0 0.8333333 0 0.8 0.00 0.0 0.8333333 0.00
## 13 0.5 0 0.0 0.3333333 0 0.4 0.00 0.0 0.3333333 0.00
## 4 0.0 1 0.0 0.3333333 0 0.8 0.00 0.0 0.1666667 0.00
## 16 0.5 0 0.0 0.0000000 0 0.0 0.75 0.0 0.0000000 0.00
## 6 0.0 0 0.0 0.5000000 0 0.6 0.00 1.0 0.8333333 0.75
## 1 0.0 0 0.0 0.0000000 0 0.8 0.00 0.0 0.0000000 0.00
## 8 0.5 0 0.0 0.5000000 0 0.8 0.00 0.0 0.3333333 0.00
## 5 0.0 0 0.0 0.5000000 0 0.4 0.00 0.4 0.3333333 1.00
## 17 0.0 0 0.4 0.3333333 0 0.2 0.00 0.0 0.0000000 1.00
## 15 0.5 0 0.0 0.3333333 1 0.0 0.00 0.0 0.1666667 0.00
## 10 0.0 0 0.0 0.5000000 0 0.8 0.00 0.0 1.0000000 1.00
## 11 0.0 0 0.4 0.8333333 0 0.8 0.00 0.0 0.5000000 0.00
## 9 0.0 0 0.0 0.3333333 0 0.8 0.00 0.0 0.5000000 0.00
## 18 0.0 0 0.0 0.8333333 0 0.6 0.00 0.0 0.3333333 0.00
## 3 0.0 0 0.0 0.3333333 0 1.0 0.00 0.0 0.3333333 0.00
## 20 1.0 0 0.0 0.3333333 0 0.0 0.75 0.0 0.0000000 0.00
## 14 0.5 0 0.0 0.3333333 1 0.0 1.00 0.0 1.0000000 0.00
## 19 0.0 0 1.0 1.0000000 0 0.0 0.00 0.0 0.3333333 1.00
## 12 0.0 0 0.0 0.3333333 0 0.0 0.00 0.0 0.5000000 0.00
## 7 0.0 0 0.0 0.5000000 0 0.8 0.00 0.4 0.3333333 0.50
## Salrep Achmil Poatri Chealb Elyrep Sagpro Plalan Agrsto Lolper
## 2 0.0 0.75 0.7777778 0 0.6666667 0.0 0.0 0.000 0.7142857
## 13 0.0 0.00 1.0000000 1 0.0000000 0.4 0.0 0.625 0.0000000
## 4 0.0 0.00 0.5555556 0 0.6666667 1.0 0.0 1.000 0.7142857
## 16 0.0 0.00 0.2222222 0 0.0000000 0.0 0.0 0.875 0.0000000
## 6 0.0 0.50 0.4444444 0 0.0000000 0.0 1.0 0.000 0.8571429
## 1 0.0 0.25 0.2222222 0 0.6666667 0.0 0.0 0.000 1.0000000
## 8 0.0 0.00 0.4444444 0 0.0000000 0.4 0.0 0.500 0.5714286
## 5 0.0 0.50 0.6666667 0 0.6666667 0.0 1.0 0.000 0.2857143
## 17 0.0 0.50 0.0000000 0 0.0000000 0.0 0.4 0.000 0.0000000
## 15 0.0 0.00 0.0000000 0 0.0000000 0.0 0.0 0.500 0.0000000
## 10 0.0 1.00 0.4444444 0 0.0000000 0.0 0.6 0.000 0.8571429
## 11 0.0 0.00 0.0000000 0 0.0000000 0.4 0.6 0.000 1.0000000
## 9 0.0 0.00 0.5555556 0 1.0000000 0.4 0.0 0.375 0.2857143
## 18 0.6 0.00 0.0000000 0 0.0000000 0.0 0.6 0.000 0.2857143
## 3 0.0 0.00 0.6666667 0 0.6666667 0.0 0.0 0.500 0.8571429
## 20 1.0 0.00 0.0000000 0 0.0000000 0.0 0.0 0.625 0.0000000
## 14 0.0 0.00 0.0000000 0 0.0000000 0.0 0.0 0.500 0.0000000
## 19 0.6 0.00 0.0000000 0 0.0000000 0.6 0.0 0.000 0.0000000
## 12 0.0 0.00 0.4444444 0 0.0000000 0.8 0.0 0.500 0.0000000
## 7 0.0 0.50 0.5555556 0 0.0000000 0.0 1.0 0.000 0.8571429
## Alogen Brohor
## 2 0.250 1.00
## 13 0.625 0.00
## 4 0.250 0.75
## 16 0.500 0.00
## 6 0.000 0.00
## 1 0.000 0.00
## 8 0.625 0.00
## 5 0.000 0.50
## 17 0.000 0.00
## 15 0.000 0.00
## 10 0.000 1.00
## 11 0.000 0.00
## 9 0.375 0.00
## 18 0.000 0.00
## 3 0.875 0.00
## 20 0.000 0.00
## 14 0.000 0.00
## 19 0.000 0.00
## 12 1.000 0.00
## 7 0.000 0.501h) Use the vegan function decostand() to standardize the dataset to mean=1 and variance=1
decostand(dunebio,"standardize")## Belper Empnig Junbuf Junart Airpra Elepal Rumace
## 2 2.2596374 -0.2236068 -0.4686477 -0.5559102 -0.3179054 -0.5298423 -0.4990282
## 13 -0.6250061 -0.2236068 1.6943416 -0.5559102 -0.3179054 -0.5298423 -0.4990282
## 4 1.2980896 -0.2236068 -0.4686477 -0.5559102 -0.3179054 -0.5298423 -0.4990282
## 16 -0.6250061 -0.2236068 -0.4686477 1.2971238 -0.3179054 2.8611484 -0.4990282
## 6 -0.6250061 -0.2236068 -0.4686477 -0.5559102 -0.3179054 -0.5298423 2.8278264
## 1 -0.6250061 -0.2236068 -0.4686477 -0.5559102 -0.3179054 -0.5298423 -0.4990282
## 8 -0.6250061 -0.2236068 -0.4686477 1.9148018 -0.3179054 1.1656531 -0.4990282
## 5 1.2980896 -0.2236068 -0.4686477 -0.5559102 -0.3179054 -0.5298423 2.2733506
## 17 -0.6250061 -0.2236068 -0.4686477 -0.5559102 2.2253377 -0.5298423 -0.4990282
## 15 -0.6250061 -0.2236068 -0.4686477 1.2971238 -0.3179054 1.5895269 -0.4990282
## 10 1.2980896 -0.2236068 -0.4686477 -0.5559102 -0.3179054 -0.5298423 -0.4990282
## 11 -0.6250061 -0.2236068 -0.4686477 -0.5559102 -0.3179054 -0.5298423 -0.4990282
## 9 -0.6250061 -0.2236068 2.4153380 1.9148018 -0.3179054 -0.5298423 0.6099233
## 18 1.2980896 -0.2236068 -0.4686477 -0.5559102 -0.3179054 -0.5298423 -0.4990282
## 3 1.2980896 -0.2236068 -0.4686477 -0.5559102 -0.3179054 -0.5298423 -0.4990282
## 20 -0.6250061 -0.2236068 -0.4686477 1.9148018 -0.3179054 1.1656531 -0.4990282
## 14 -0.6250061 -0.2236068 -0.4686477 -0.5559102 -0.3179054 1.1656531 -0.4990282
## 19 -0.6250061 4.2485292 -0.4686477 -0.5559102 3.4969592 -0.5298423 -0.4990282
## 12 -0.6250061 -0.2236068 2.4153380 -0.5559102 -0.3179054 -0.5298423 0.6099233
## 7 -0.6250061 -0.2236068 0.9733452 -0.5559102 -0.3179054 -0.5298423 1.1643991
## Viclat Brarut Ranfla Cirarv Hyprad Leoaut Potpal
## 2 -0.3823007 -1.286103 -0.596107 -0.2236068 -0.3645567 1.4749716 -0.3248931
## 13 -0.3823007 -1.286103 1.107056 -0.2236068 -0.3645567 -0.4489044 -0.3248931
## 4 -0.3823007 -0.236223 -0.596107 4.2485292 -0.3645567 -0.4489044 -0.3248931
## 16 -0.3823007 0.813657 1.107056 -0.2236068 -0.3645567 -1.7314884 -0.3248931
## 6 -0.3823007 1.863537 -0.596107 -0.2236068 -0.3645567 0.1923876 -0.3248931
## 1 -0.3823007 -1.286103 -0.596107 -0.2236068 -0.3645567 -1.7314884 -0.3248931
## 8 -0.3823007 -0.236223 1.107056 -0.2236068 -0.3645567 0.1923876 -0.3248931
## 5 -0.3823007 -0.236223 -0.596107 -0.2236068 -0.3645567 0.1923876 -0.3248931
## 17 -0.3823007 -1.286103 -0.596107 -0.2236068 1.2556951 -0.4489044 -0.3248931
## 15 -0.3823007 0.813657 1.107056 -0.2236068 -0.3645567 -0.4489044 2.9240383
## 10 1.5292029 -0.236223 -0.596107 -0.2236068 -0.3645567 0.1923876 -0.3248931
## 11 3.4407065 0.813657 -0.596107 -0.2236068 1.2556951 1.4749716 -0.3248931
## 9 -0.3823007 -0.236223 -0.596107 -0.2236068 -0.3645567 -0.4489044 -0.3248931
## 18 1.5292029 1.863537 -0.596107 -0.2236068 -0.3645567 1.4749716 -0.3248931
## 3 -0.3823007 -0.236223 -0.596107 -0.2236068 -0.3645567 -0.4489044 -0.3248931
## 20 -0.3823007 0.813657 2.810219 -0.2236068 -0.3645567 -0.4489044 -0.3248931
## 14 -0.3823007 -1.286103 1.107056 -0.2236068 -0.3645567 -0.4489044 2.9240383
## 19 -0.3823007 0.288717 -0.596107 -0.2236068 3.6860728 2.1162636 -0.3248931
## 12 -0.3823007 0.813657 -0.596107 -0.2236068 -0.3645567 -0.4489044 -0.3248931
## 7 -0.3823007 -0.236223 -0.596107 -0.2236068 -0.3645567 0.1923876 -0.3248931
## Poapra Calcus Tripra Trirep Antodo Salrep Achmil
## 2 0.8663817 -0.4047136 -0.3645567 1.3951437 -0.6174222 -0.3943958 1.7746310
## 13 -0.2165954 -0.4047136 -0.3645567 -0.1842643 -0.6174222 -0.3943958 -0.6453204
## 4 0.8663817 -0.4047136 -0.3645567 -0.7107336 -0.6174222 -0.3943958 -0.6453204
## 16 -1.2995726 2.0235680 -0.3645567 -1.2372029 -0.6174222 -0.3943958 -0.6453204
## 6 0.3248931 -0.4047136 3.6860728 1.3951437 1.1466413 -0.3943958 0.9679805
## 1 0.8663817 -0.4047136 -0.3645567 -1.2372029 -0.6174222 -0.3943958 0.1613301
## 8 0.8663817 -0.4047136 -0.3645567 -0.1842643 -0.6174222 -0.3943958 -0.6453204
## 5 -0.2165954 -0.4047136 1.2556951 -0.1842643 1.7346624 -0.3943958 0.9679805
## 17 -0.7580840 -0.4047136 -0.3645567 -1.2372029 1.7346624 -0.3943958 0.9679805
## 15 -1.2995726 -0.4047136 -0.3645567 -0.7107336 -0.6174222 -0.3943958 -0.6453204
## 10 0.8663817 -0.4047136 -0.3645567 1.9216130 1.7346624 -0.3943958 2.5812814
## 11 0.8663817 -0.4047136 -0.3645567 0.3422051 -0.6174222 -0.3943958 -0.6453204
## 9 0.8663817 -0.4047136 -0.3645567 0.3422051 -0.6174222 -0.3943958 -0.6453204
## 18 0.3248931 -0.4047136 -0.3645567 -0.1842643 -0.6174222 1.7568540 -0.6453204
## 3 1.4078703 -0.4047136 -0.3645567 -0.1842643 -0.6174222 -0.3943958 -0.6453204
## 20 -1.2995726 2.0235680 -0.3645567 -1.2372029 -0.6174222 3.1910205 -0.6453204
## 14 -1.2995726 2.8329952 -0.3645567 1.9216130 -0.6174222 -0.3943958 -0.6453204
## 19 -1.2995726 -0.4047136 -0.3645567 -0.1842643 1.7346624 1.7568540 -0.6453204
## 12 -1.2995726 -0.4047136 -0.3645567 0.3422051 -0.6174222 -0.3943958 -0.6453204
## 7 0.8663817 -0.4047136 1.2556951 -0.1842643 0.5586201 -0.3943958 0.9679805
## Poatri Chealb Elyrep Sagpro Plalan Agrsto Lolper
## 2 1.3677199 -0.2236068 1.2980896 -0.6426846 -0.6668859 -0.8944272 0.7429511
## 13 2.0782237 4.2485292 -0.6250061 0.6426846 -0.6668859 0.9689628 -1.0259800
## 4 0.6572160 -0.2236068 1.2980896 2.5707383 -0.6668859 2.0869968 0.7429511
## 16 -0.4085397 -0.2236068 -0.6250061 -0.6426846 -0.6668859 1.7143188 -1.0259800
## 6 0.3019641 -0.2236068 -0.6250061 -0.6426846 1.8980600 -0.8944272 1.0967373
## 1 -0.4085397 -0.2236068 1.2980896 -0.6426846 -0.6668859 -0.8944272 1.4505235
## 8 0.3019641 -0.2236068 -0.6250061 0.6426846 -0.6668859 0.5962848 0.3891648
## 5 1.0124679 -0.2236068 1.2980896 -0.6426846 1.8980600 -0.8944272 -0.3184076
## 17 -1.1190435 -0.2236068 -0.6250061 -0.6426846 0.3590924 -0.8944272 -1.0259800
## 15 -1.1190435 -0.2236068 -0.6250061 -0.6426846 -0.6668859 0.5962848 -1.0259800
## 10 0.3019641 -0.2236068 -0.6250061 -0.6426846 0.8720816 -0.8944272 1.0967373
## 11 -1.1190435 -0.2236068 -0.6250061 0.6426846 0.8720816 -0.8944272 1.4505235
## 9 0.6572160 -0.2236068 2.2596374 0.6426846 -0.6668859 0.2236068 -0.3184076
## 18 -1.1190435 -0.2236068 -0.6250061 -0.6426846 0.8720816 -0.8944272 -0.3184076
## 3 1.0124679 -0.2236068 1.2980896 -0.6426846 -0.6668859 0.5962848 1.0967373
## 20 -1.1190435 -0.2236068 -0.6250061 -0.6426846 -0.6668859 0.9689628 -1.0259800
## 14 -1.1190435 -0.2236068 -0.6250061 -0.6426846 -0.6668859 0.5962848 -1.0259800
## 19 -1.1190435 -0.2236068 -0.6250061 1.2853692 -0.6668859 -0.8944272 -1.0259800
## 12 0.3019641 -0.2236068 -0.6250061 1.9280538 -0.6668859 0.5962848 -1.0259800
## 7 0.6572160 -0.2236068 -0.6250061 -0.6426846 1.8980600 -0.8944272 1.0967373
## Alogen Brohor
## 2 0.07610968 2.3056941
## 13 1.21775483 -0.5320832
## 4 0.07610968 1.5962497
## 16 0.83720645 -0.5320832
## 6 -0.68498709 -0.5320832
## 1 -0.68498709 -0.5320832
## 8 1.21775483 -0.5320832
## 5 -0.68498709 0.8868054
## 17 -0.68498709 -0.5320832
## 15 -0.68498709 -0.5320832
## 10 -0.68498709 2.3056941
## 11 -0.68498709 -0.5320832
## 9 0.45665806 -0.5320832
## 18 -0.68498709 -0.5320832
## 3 1.97885160 -0.5320832
## 20 -0.68498709 -0.5320832
## 14 -0.68498709 -0.5320832
## 19 -0.68498709 -0.5320832
## 12 2.35939999 -0.5320832
## 7 -0.68498709 0.8868054- Write a function to calculate the observed richness of a vector representing the abundances of different species in a community.
obsrich<- function(x){sum(x > 0)}
obsrich(dunebio)## [1] 197- Write a function to calculate the Shannon-Wiener diversity index of a vector representing the abundances of different species in a community.
H<- function(x){-sum(x[x>0]/sum(x[x>0])*log(x[x>0]/sum(x[x>0])))}
#pi = NumberofSpecies i /TotalNumberofSamples
H(dunebio)## [1] 5.181888- Write a function that calculates both the observed richness and the Shannon-Wiener diversity index for each community in a matrix, and writes an output table with the results. You can use the vegan built-in functions. Test your function with the dune_bio.txt dataset.
richdiv<- function(x) {
spec<-specnumber(x)
richdiv.df <- as.data.frame(spec)
richdiv.df$div<-diversity(x, index="shannon")
richdiv.df
}
richdiv(dunebio)## spec div
## 2 10 2.252516
## 13 10 2.099638
## 4 13 2.426779
## 16 8 1.959795
## 6 11 2.345946
## 1 5 1.440482
## 8 12 2.434898
## 5 14 2.544421
## 17 7 1.876274
## 15 8 1.979309
## 10 12 2.398613
## 11 9 2.106065
## 9 13 2.493568
## 18 9 2.079387
## 3 10 2.193749
## 20 8 2.048270
## 14 7 1.863680
## 19 9 2.134024
## 12 9 2.114495
## 7 13 2.471733- Make a rank-abundance curve using the second site (13) of the dune_bio.txt dataset. Fit a lognormal model to the data.
hist(as.numeric(dunebio[2,]), main = "Histogram")
plot(rad.lognormal(dunebio[2,]), lty=2, lwd=2)
plot(radfit(dunebio))## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
radfit(dunebio[2,])##
## RAD models, family poisson
## No. of species 10, total abundance 33
##
## par1 par2 par3 Deviance AIC BIC
## Null 3.91659 32.99179 32.99179
## Preemption 0.20971 2.02192 33.09712 33.39970
## Lognormal 0.99822 0.71061 1.13244 34.20764 34.81281
## Zipf 0.27866 -0.79375 0.79058 33.86578 34.47095
## Mandelbrot 0.40348 -0.95679 0.51608 0.75881 35.83401 36.74176rad.lognormal(dunebio[2,])##
## RAD model: Log-Normal
## Family: poisson
## No. of species: 10
## Total abundance: 33
##
## log.mu log.sigma Deviance AIC BIC
## 0.9982177 0.7106063 1.1324449 34.2076395 34.8128097- Create a Euclidean distance matrix after standardization using the first (A1), second (Moisture) and fifth (Manure) variables of dune_env.txt. Then create a Bray-Curtis distance matrix using dune_bio.txt. Perform a Mantel test on both distance matrices and plot the relationship.
duneenv<-read.table("dune_env.txt", sep="\t", header=T, row.names=1)
#standardized using the first (A1), second (Moisture) and fifth (Manure) variables of dune_env.txt.
duneenveuc<-vegdist(decostand(duneenv[,c(1,2,5)], method="standardize"), method="euclidean")
#create a Bray-Curtis distance matrix
dunebiobray<-vegdist(dunebio, method="bray")
#Perform a Mantel test
mantel(duneenveuc, dunebiobray)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = duneenveuc, ydis = dunebiobray)
##
## Mantel statistic r: 0.5496
## Significance: 0.001
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.142 0.184 0.221 0.265
## Permutation: free
## Number of permutations: 999#plot the relationship
plot(duneenveuc, dunebiobray)
PRACTICE PROBLEMS – Cluster analysis
- Make dendrograms of the results of hierarchical clustering using the single, average and complete methods. Use the dune_bio.txt dataset after creating a Bray-Curtis distance matrix.
plot(hclust(dunebiobray, method="single"))
plot(hclust(dunebiobray, method="average"))
plot(hclust(dunebiobray, method="complete"))
- Perform k-means partitions from 2 to 6 on the dune_bio.txt dataset and select the best partition using the Calinski criterion. Visualize the results.
dunebiok<-cascadeKM(dunebio, inf.gr=2, sup.gr=6)
plot(dunebiok, sortg=T)
dunebiok$results## 2 groups 3 groups 4 groups 5 groups 6 groups
## SSE 1218.500000 950.516667 777.333333 676.500000 593.750000
## calinski 5.611243 5.793253 5.633047 5.110033 4.737482PRACTICE PROBLEMS – Unconstrained ordination
- Perform a correspondence analysis (CA) on the dune_bio.txt dataset.
as.data.frame(dunebio)## Belper Empnig Junbuf Junart Airpra Elepal Rumace Viclat Brarut Ranfla Cirarv
## 2 3 0 0 0 0 0 0 0 0 0 0
## 13 0 0 3 0 0 0 0 0 0 2 0
## 4 2 0 0 0 0 0 0 0 2 0 2
## 16 0 0 0 3 0 8 0 0 4 2 0
## 6 0 0 0 0 0 0 6 0 6 0 0
## 1 0 0 0 0 0 0 0 0 0 0 0
## 8 0 0 0 4 0 4 0 0 2 2 0
## 5 2 0 0 0 0 0 5 0 2 0 0
## 17 0 0 0 0 2 0 0 0 0 0 0
## 15 0 0 0 3 0 5 0 0 4 2 0
## 10 2 0 0 0 0 0 0 1 2 0 0
## 11 0 0 0 0 0 0 0 2 4 0 0
## 9 0 0 4 4 0 0 2 0 2 0 0
## 18 2 0 0 0 0 0 0 1 6 0 0
## 3 2 0 0 0 0 0 0 0 2 0 0
## 20 0 0 0 4 0 4 0 0 4 4 0
## 14 0 0 0 0 0 4 0 0 0 2 0
## 19 0 2 0 0 3 0 0 0 3 0 0
## 12 0 0 4 0 0 0 2 0 4 0 0
## 7 0 0 2 0 0 0 3 0 2 0 0
## Hyprad Leoaut Potpal Poapra Calcus Tripra Trirep Antodo Salrep Achmil Poatri
## 2 0 5 0 4 0 0 5 0 0 3 7
## 13 0 2 0 2 0 0 2 0 0 0 9
## 4 0 2 0 4 0 0 1 0 0 0 5
## 16 0 0 0 0 3 0 0 0 0 0 2
## 6 0 3 0 3 0 5 5 3 0 2 4
## 1 0 0 0 4 0 0 0 0 0 1 2
## 8 0 3 0 4 0 0 2 0 0 0 4
## 5 0 3 0 2 0 2 2 4 0 2 6
## 17 2 2 0 1 0 0 0 4 0 2 0
## 15 0 2 2 0 0 0 1 0 0 0 0
## 10 0 3 0 4 0 0 6 4 0 4 4
## 11 2 5 0 4 0 0 3 0 0 0 0
## 9 0 2 0 4 0 0 3 0 0 0 5
## 18 0 5 0 3 0 0 2 0 3 0 0
## 3 0 2 0 5 0 0 2 0 0 0 6
## 20 0 2 0 0 3 0 0 0 5 0 0
## 14 0 2 2 0 4 0 6 0 0 0 0
## 19 5 6 0 0 0 0 2 4 3 0 0
## 12 0 2 0 0 0 0 3 0 0 0 4
## 7 0 3 0 4 0 2 2 2 0 2 5
## Chealb Elyrep Sagpro Plalan Agrsto Lolper Alogen Brohor
## 2 0 4 0 0 0 5 2 4
## 13 1 0 2 0 5 0 5 0
## 4 0 4 5 0 8 5 2 3
## 16 0 0 0 0 7 0 4 0
## 6 0 0 0 5 0 6 0 0
## 1 0 4 0 0 0 7 0 0
## 8 0 0 2 0 4 4 5 0
## 5 0 4 0 5 0 2 0 2
## 17 0 0 0 2 0 0 0 0
## 15 0 0 0 0 4 0 0 0
## 10 0 0 0 3 0 6 0 4
## 11 0 0 2 3 0 7 0 0
## 9 0 6 2 0 3 2 3 0
## 18 0 0 0 3 0 2 0 0
## 3 0 4 0 0 4 6 7 0
## 20 0 0 0 0 5 0 0 0
## 14 0 0 0 0 4 0 0 0
## 19 0 0 3 0 0 0 0 0
## 12 0 0 4 0 4 0 8 0
## 7 0 0 0 5 0 6 0 2dunebioca<-cca(dunebio)
summary(dunebioca)##
## Call:
## cca(X = dunebio)
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 2.115 1
## Unconstrained 2.115 1
##
## Eigenvalues, and their contribution to the scaled Chi-square
##
## Importance of components:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7
## Eigenvalue 0.5360 0.4001 0.2598 0.17598 0.14476 0.10791 0.09247
## Proportion Explained 0.2534 0.1892 0.1228 0.08319 0.06844 0.05102 0.04372
## Cumulative Proportion 0.2534 0.4426 0.5654 0.64858 0.71702 0.76804 0.81175
## CA8 CA9 CA10 CA11 CA12 CA13 CA14
## Eigenvalue 0.08091 0.07332 0.05630 0.04826 0.04125 0.03523 0.020529
## Proportion Explained 0.03825 0.03466 0.02661 0.02282 0.01950 0.01665 0.009705
## Cumulative Proportion 0.85000 0.88467 0.91128 0.93410 0.95360 0.97025 0.979955
## CA15 CA16 CA17 CA18 CA19
## Eigenvalue 0.014911 0.009074 0.007938 0.007002 0.003477
## Proportion Explained 0.007049 0.004290 0.003753 0.003310 0.001644
## Cumulative Proportion 0.987004 0.991293 0.995046 0.998356 1.000000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CA1 CA2 CA3 CA4 CA5 CA6
## Belper -0.50018 -0.35503 -0.15239 -0.704153 -0.058546 -0.07308
## Empnig -0.69027 3.26420 1.95716 -0.176936 -0.073518 0.16083
## Junbuf 0.08157 -0.68074 1.00545 1.078390 0.268360 -0.24168
## Junart 1.27580 0.09963 -0.09320 0.005536 0.289410 0.78247
## Airpra -1.00434 3.06749 1.33773 0.194305 -1.081813 0.53699
## Elepal 1.76383 0.34562 -0.57336 -0.002976 -0.332396 0.14688
## Rumace -0.65289 -0.25525 -0.59728 1.160164 0.255849 0.32730
## Viclat -0.61893 0.37140 -0.46057 -1.000375 1.162652 -1.44971
## Brarut 0.18222 0.26477 -0.16606 0.064009 0.576334 -0.07741
## Ranfla 1.55886 0.30700 -0.29765 0.046974 -0.008747 0.14744
## Cirarv -0.05647 -0.76398 0.91793 -1.175919 -0.384024 0.13985
## Hyprad -0.85408 2.52821 1.13951 -0.175115 -0.311874 -0.11177
## Leoaut -0.19566 0.38884 0.03975 -0.130392 0.141232 -0.23717
## Potpal 1.91690 0.52150 -1.18215 -0.021738 -1.359988 -1.31207
## Poapra -0.38919 -0.32999 -0.02015 -0.358371 0.079296 0.05165
## Calcus 1.95199 0.56743 -0.85948 -0.098969 -0.556737 -0.23282
## Tripra -0.88116 -0.09792 -1.18172 1.282429 0.325706 0.33388
## Trirep -0.07666 -0.02032 -0.20594 0.026462 -0.186748 -0.53957
## Antodo -0.96676 1.08361 -0.17188 0.459788 -0.607533 0.30425
## Salrep 0.61035 1.54868 0.04970 -0.607136 1.429729 0.55183
## Achmil -0.90859 0.08461 -0.58636 -0.008919 -0.660183 0.18877
## Poatri -0.18185 -0.53997 0.23388 0.178834 -0.155342 0.07584
## Chealb 0.42445 -0.84402 1.59029 1.248755 -0.207480 -0.87566
## Elyrep -0.37074 -0.74148 0.26238 -0.566308 -0.270122 0.72624
## Sagpro 0.00364 0.01719 1.11570 0.066981 0.186654 -0.32463
## Plalan -0.84058 0.24886 -0.78066 0.371149 0.271377 -0.11989
## Agrsto 0.93378 -0.20651 0.28165 0.024293 -0.139326 0.02256
## Lolper -0.50272 -0.35955 -0.21821 -0.474727 0.101494 0.01594
## Alogen 0.40088 -0.61839 0.85013 0.346740 0.016509 -0.10169
## Brohor -0.65762 -0.40634 -0.30685 -0.496751 -0.561358 -0.07004
##
##
## Site scores (weighted averages of species scores)
##
## CA1 CA2 CA3 CA4 CA5 CA6
## 2 -0.63268 -0.6958357 -0.09708 -1.18695 -0.97686 -0.06575
## 13 0.42445 -0.8440195 1.59029 1.24876 -0.20748 -0.87566
## 4 -0.05647 -0.7639784 0.91793 -1.17592 -0.38402 0.13985
## 16 2.00229 0.1090627 -0.33414 0.33760 -0.50097 0.76159
## 6 -0.85633 -0.0005408 -1.39735 1.59909 0.65494 0.19386
## 1 -0.81167 -1.0826714 -0.14479 -2.10665 -0.39287 1.83462
## 8 0.76268 -0.2968459 0.35648 -0.10772 0.17507 0.36444
## 5 -0.95293 -0.1846015 -0.95609 0.86853 -0.34552 0.98333
## 17 -1.47545 2.7724102 0.40859 0.75117 -2.59425 1.10122
## 15 1.91384 0.5079036 -0.96567 0.04227 -0.50681 -0.19370
## 10 -0.87885 -0.0353136 -0.82987 -0.68053 -0.75438 -0.81070
## 11 -0.64223 0.4440332 -0.17371 -1.09684 1.37462 -2.00626
## 9 0.09693 -0.7864314 0.86492 0.40090 0.28704 1.02783
## 18 -0.31241 0.6328355 -0.66501 -1.12728 2.65575 -0.97565
## 3 -0.10148 -0.9128732 0.68815 -0.68137 -0.08709 0.28678
## 20 1.94438 1.0688809 -0.66595 -0.55317 1.59606 1.70292
## 14 1.91996 0.5351062 -1.39863 -0.08575 -2.21317 -2.43044
## 19 -0.69027 3.2642026 1.95716 -0.17694 -0.07352 0.16083
## 12 0.28557 -0.6656161 1.64423 1.71496 0.65381 -1.17376
## 7 -0.87149 -0.2547040 -0.86830 0.90468 0.17385 0.034461a) How much variation is explained by the two first axes?
sum(0.2534,0.1892)## [1] 0.44261b) Make a screeplot of the CA results
screeplot(dunebioca)
1c) Plot the ordination results of the sites
plot(dunebioca, scaling=2, xlab="CA1 (0.2534)", ylab="CA2 (0.1892)") 
- Perform a nonmetric multidimensional scaling (NMDS) on the dune_bio.txt dataset after calculating Bray-Curtis distances among sites.
dunebionmds<-metaMDS(dunebio, distance="bray", k=2)## Run 0 stress 0.1192678
## Run 1 stress 0.1183186
## ... New best solution
## ... Procrustes: rmse 0.0202709 max resid 0.06495715
## Run 2 stress 0.1192678
## Run 3 stress 0.1192679
## Run 4 stress 0.1192686
## Run 5 stress 0.1183186
## ... Procrustes: rmse 7.720963e-05 max resid 0.0002342301
## ... Similar to previous best
## Run 6 stress 0.1192678
## Run 7 stress 0.1812932
## Run 8 stress 0.1183186
## ... New best solution
## ... Procrustes: rmse 1.616026e-05 max resid 5.037771e-05
## ... Similar to previous best
## Run 9 stress 0.1192683
## Run 10 stress 0.1183186
## ... Procrustes: rmse 2.586111e-05 max resid 8.049044e-05
## ... Similar to previous best
## Run 11 stress 0.1183186
## ... Procrustes: rmse 2.651651e-05 max resid 8.200842e-05
## ... Similar to previous best
## Run 12 stress 0.1192679
## Run 13 stress 0.1808915
## Run 14 stress 0.1809578
## Run 15 stress 0.1812936
## Run 16 stress 0.1192678
## Run 17 stress 0.1192687
## Run 18 stress 0.2075713
## Run 19 stress 0.1183186
## ... Procrustes: rmse 2.210693e-05 max resid 6.769228e-05
## ... Similar to previous best
## Run 20 stress 0.1809578
## *** Solution reached2a) What is the stress?
dunebionmds$stress## [1] 0.11831862b) Make a Shepard plot of the NMDS results
stressplot(dunebionmds)
2c) Plot the ordination results of the sites
plot(dunebionmds, type="t", display="sites")
rownames(duneenv)==rownames(dunebio)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUEduneenv$Axis01<-dunebionmds$points[,1]
duneenv$Axis02<-dunebionmds$points[,2]
duneenv$Richness<-specnumber(dunebio)2d) [Advanced – ggplot2] Plot the ordination results using ggplot2. Use the dune_env.txt dataset to make the size of the points proportional to the site richness (number of species) and the color to represent the variable Management.
ggplot(duneenv, aes(x=Axis01, y=Axis02))+
geom_point(aes(fill=Management, size=Richness), alpha=0.6, pch=21)
#PRACTICE PROBLEMS – Constrained ordination
- Perform a constrained correspondence analysis (CCA) on the dune_bio.txt dataset using the first (A1), second (Moisture) and fifth (Manure) variables of dune_env.txt as explanatory.
dunecca<-cca(dunebio ~ A1 + Moisture + Manure, data=duneenv)
summary(dunecca)##
## Call:
## cca(formula = dunebio ~ A1 + Moisture + Manure, data = duneenv)
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 2.1153 1.0000
## Constrained 0.7692 0.3637
## Unconstrained 1.3460 0.6363
##
## Eigenvalues, and their contribution to the scaled Chi-square
##
## Importance of components:
## CCA1 CCA2 CCA3 CA1 CA2 CA3 CA4
## Eigenvalue 0.4264 0.2347 0.10812 0.3124 0.2101 0.17518 0.13748
## Proportion Explained 0.2016 0.1110 0.05112 0.1477 0.0993 0.08282 0.06499
## Cumulative Proportion 0.2016 0.3125 0.36366 0.5114 0.6107 0.69348 0.75847
## CA5 CA6 CA7 CA8 CA9 CA10 CA11
## Eigenvalue 0.09442 0.09295 0.07380 0.06360 0.05541 0.04663 0.021117
## Proportion Explained 0.04464 0.04394 0.03489 0.03007 0.02620 0.02204 0.009983
## Cumulative Proportion 0.80311 0.84705 0.88194 0.91201 0.93820 0.96025 0.970230
## CA12 CA13 CA14 CA15 CA16
## Eigenvalue 0.019730 0.016183 0.012531 0.009960 0.004566
## Proportion Explained 0.009327 0.007651 0.005924 0.004709 0.002159
## Cumulative Proportion 0.979558 0.987208 0.993132 0.997841 1.000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## CCA1 CCA2 CCA3
## Eigenvalue 0.4264 0.2347 0.1081
## Proportion Explained 0.5543 0.3051 0.1406
## Cumulative Proportion 0.5543 0.8594 1.0000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## Belper -0.72508 0.09015 0.28678 -0.14997 -0.08941 0.64602
## Empnig 0.92386 -1.49359 -1.29239 2.87128 1.31651 0.74936
## Junbuf 0.50547 0.17496 -0.35707 0.33578 -1.58298 -0.40514
## Junart 1.08445 -0.22093 -0.33048 -0.84568 0.15861 -0.23372
## Airpra 0.32923 -1.52131 -0.75777 2.80560 1.49278 0.23227
## Elepal 1.39255 0.01967 0.39847 -0.90296 0.89564 -0.18234
## Rumace -0.56966 0.07756 0.44317 0.23147 -0.33959 -1.25112
## Viclat -0.95631 -1.05874 0.13224 -0.42277 -0.16218 0.84521
## Brarut 0.05861 -0.29093 0.14295 -0.17325 0.10170 -0.08514
## Ranfla 1.30695 -0.17793 0.06749 -0.84091 0.45935 -0.18497
## Cirarv -0.40316 1.49565 -0.09657 0.17279 0.57950 1.51909
## Hyprad 0.14257 -1.37899 -0.68907 2.13762 1.11990 0.54595
## Leoaut -0.10805 -0.39424 0.03343 0.21126 0.05153 0.18568
## Potpal 1.82080 -0.59463 2.50587 -0.53778 0.46936 0.51885
## Poapra -0.47298 0.19818 -0.15806 -0.08959 -0.09919 0.21588
## Calcus 1.32589 -0.43845 0.22643 -1.29121 0.99188 -0.18948
## Tripra -0.94282 0.14003 0.52488 0.16097 0.18159 -1.70499
## Trirep -0.03691 -0.24820 0.22610 -0.03138 -0.23464 0.04178
## Antodo -0.42848 -0.66783 0.01590 1.13407 0.58011 -0.56583
## Salrep 0.38937 -1.51269 -0.78060 -0.47447 0.69625 0.26139
## Achmil -0.84528 -0.28199 0.08034 0.23663 0.15351 -0.33457
## Poatri -0.09591 0.48672 -0.08195 0.10052 -0.39496 -0.08346
## Chealb 1.33134 1.08879 -0.17910 0.92504 -1.53277 -0.26877
## Elyrep -0.45995 0.49110 -0.07018 -0.09270 -0.32202 0.53794
## Sagpro 0.36673 0.22891 -0.41200 0.63169 -0.40176 0.48274
## Plalan -0.89463 -0.37036 0.37142 0.16434 0.15598 -0.65298
## Agrsto 0.80663 0.42851 0.03513 -0.31737 0.09257 0.16535
## Lolper -0.68266 0.25866 -0.13366 -0.17559 0.06653 0.20307
## Alogen 0.52884 0.74865 -0.28721 0.11872 -0.57814 0.07754
## Brohor -0.77674 0.11771 0.03195 -0.07993 -0.06235 0.29678
##
##
## Site scores (weighted averages of species scores)
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## 2 -0.8544 0.57195 -0.04460 -0.44645 -0.473383 1.10419
## 13 0.7656 1.43234 -1.04660 0.92504 -1.532773 -0.26877
## 4 -0.1565 1.36916 -0.67602 0.17279 0.579503 1.51909
## 16 2.0459 0.46834 0.70504 -0.99503 1.669766 -0.75493
## 6 -1.0571 -0.29557 1.50936 0.06506 0.216688 -2.09457
## 1 -1.2439 1.24492 -0.99278 -0.48730 0.708935 1.13256
## 8 0.8113 0.66762 -0.54772 -0.28201 0.298269 -0.31177
## 5 -1.1247 -0.04619 1.17587 0.61920 0.008872 -0.95018
## 17 -0.7722 -2.94478 -1.24411 2.70708 1.757175 -0.54337
## 15 2.0065 -0.48090 2.87633 -0.26585 0.401900 0.57133
## 10 -1.0958 -0.42412 0.49762 -0.26374 -0.332785 0.08742
## 11 -0.7816 -0.90608 -0.28150 -0.26599 -0.008900 1.12676
## 9 0.1817 0.86595 -0.91242 -0.34851 -2.059004 -0.10377
## 18 -0.6053 -1.51952 0.07324 -0.89534 -0.298138 1.03991
## 3 -0.2093 1.35237 -0.66038 0.06196 0.171454 0.84659
## 20 1.8996 -1.40266 -0.55715 -2.22940 0.920734 -0.49851
## 14 1.9462 -0.67177 3.54931 -0.80971 0.536813 0.46637
## 19 0.2690 -3.39538 -3.20303 2.87128 1.316513 0.74936
## 12 0.6252 1.06204 -0.88731 0.77476 -2.069361 -0.26842
## 7 -1.0498 -0.01449 0.64118 -0.05748 0.266550 -1.48584
##
##
## Site constraints (linear combinations of constraining variables)
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## 2 -1.0722 -0.15065 0.02248 -0.44645 -0.473383 1.10419
## 13 1.3313 1.08879 -0.17910 0.92504 -1.532773 -0.26877
## 4 -0.4032 1.49565 -0.09657 0.17279 0.579503 1.51909
## 16 1.2912 1.04854 -0.34917 -0.99503 1.669766 -0.75493
## 6 -0.9651 -0.04331 0.47601 0.06506 0.216688 -2.09457
## 1 -1.0995 1.27129 -0.50141 -0.48730 0.708935 1.13256
## 8 1.0904 0.84728 -1.19953 -0.28201 0.298269 -0.31177
## 5 -0.6973 0.22504 1.60982 0.61920 0.008872 -0.95018
## 17 -0.5627 -1.56290 0.04417 2.70708 1.757175 -0.54337
## 15 1.9681 -0.44704 3.12947 -0.26585 0.401900 0.57133
## 10 -0.6232 -0.89889 -0.41620 -0.26374 -0.332785 0.08742
## 11 -1.1054 -0.90857 0.08601 -0.26599 -0.008900 1.12676
## 9 0.4481 -0.77218 -0.96709 -0.34851 -2.059004 -0.10377
## 18 -0.9913 -1.51891 0.77314 -0.89534 -0.298138 1.03991
## 3 -0.3898 1.50907 -0.03988 0.06196 0.171454 0.84659
## 20 0.8971 -1.52042 -1.40577 -2.22940 0.920734 -0.49851
## 14 1.6735 -0.74222 1.88228 -0.80971 0.536813 0.46637
## 19 0.9239 -1.49359 -1.29239 2.87128 1.316513 0.74936
## 12 0.7625 0.26751 0.15988 0.77476 -2.069361 -0.26842
## 7 -1.1327 0.51336 -0.43788 -0.05748 0.266550 -1.48584
##
##
## Biplot scores for constraining variables
##
## CCA1 CCA2 CCA3 CA1 CA2 CA3
## A1 0.6048 0.04018 0.7954 0 0 0
## Moisture 0.9746 -0.06076 -0.2158 0 0 0
## Manure -0.2059 0.96205 -0.1791 0 0 01a) What is the variation explained by the two first constrained axes?
sum(0.2016, 0.1110)## [1] 0.312631.26%
1b) What is the adjusted R2 of the model?
RsquareAdj(dunecca)$adj.r.squared## [1] 0.24446311c) Use a permutation test to assess the overall statistical significance of the model
anova(dunecca)## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = dunebio ~ A1 + Moisture + Manure, data = duneenv)
## Df ChiSquare F Pr(>F)
## Model 3 0.76925 3.048 0.001 ***
## Residual 16 1.34602
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11d) Use a permutation test to assess the marginal statistical significance of each explanatory variable of the model. Which variables are significant?
anova(dunecca, by="margin")## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = dunebio ~ A1 + Moisture + Manure, data = duneenv)
## Df ChiSquare F Pr(>F)
## A1 1 0.13047 1.5508 0.107
## Moisture 1 0.30886 3.6714 0.001 ***
## Manure 1 0.23418 2.7837 0.001 ***
## Residual 16 1.34602
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11e) Make a triplot of the ordination results focusing on the sites
plot(dunecca, scaling=2, xlab="CCA1 (20%)", ylab="CCA2 (11%)")
- Perform a permutational multivariate analysis of variance (PERMANOVA) with 9,999 permutations using the dune_bio.txt dataset after calculating Bray-Curtis distances and the first (A1), second (Moisture) and fifth (Manure) variables of dune_env.txt as explanatory.
# class before knitting is dataframe, but after knitting one class variable changes to a function
class(dunebio)## [1] "data.frame"class(duneenv)## [1] "data.frame"#calculate Bray-Curtis distances using the first (A1), second (Moisture) and fifth (Manure) variables of duneenv as explanatory and run PERMANOVA
#adonis(vegdist(dunebio, method="bray") ~ duneenv$A1 + duneenv$Moisture + duneenv$Manure, by="margin",permutations=9999)
#error - could not knit adonisCall: adonis(formula = vegdist(dunebio, method = “bray”) ~ duneenv\(A1 + duneenv\)Moisture + duneenv\(Manure, permutations = 9999, by = "margin") Permutation: free Number of permutations: 9999 Terms added sequentially (first to last) Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) duneenv\)A1 1 0.5107 0.51069 5.8578 0.15417 1e-04 duneenv\(Moisture 1 0.7196 0.71958 8.2537 0.21722 1e-04 *** duneenv\)Manure 1 0.6874 0.68742 7.8849 0.20752 1e-04 Residuals 16 1.3949 0.08718 0.42109
Total 19 3.3126 1.00000
— Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’
Which explanatory variables are significant? Moisture and Manure are significant.
What is the explanatory power (R2) of the significant variables? Moisture R^2 is 0.20289 and Manure R^2 is 0.20752
Perform a new PERMANOVA analysis with 9,999 permutations with A1 as the explanatory variable but constrain the permutations within the Use variable.
#calculate Bray-Curtis distances using the first (A1) variable and constrain within the the Use variable of duneenv and run PERMANOVA
#adonis(vegdist(dunebio, method="bray") ~ duneenv$A1, strata = duneenv$Use, by="margin",permutations=9999)Call: adonis(formula = vegdist(dunebio, method = “bray”) ~ duneenv\(A1, permutations = 9999, strata = duneenv\)Use, by = “margin”) Blocks: strata Permutation: free Number of permutations: 9999 Terms added sequentially (first to last) Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
duneenv$A1 1 0.5107 0.51069 3.2808 0.15417 0.0118 * Residuals 18 2.8019 0.15566 0.84583
Total 19 3.3126 1.00000
— Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1