This notebook goes over a few examples of analyses which are typically done to discern patterns of benthic species. They include: multivariate analyses (community similarity and environmental drivers), and responses of individual species to environmental factors. All data are from the MAREANO Programme.
In [9]:
samples <- read.csv("tab_station.csv",header=T)
data <- read.csv("tab_videolog_files.csv",header=T)
In [1]:
# database connection details: marbunn@postgres9.imr.no:5432 pw: xxxxxxx
# run sql: select
# refstation_no,sample_no,latitude,longitude,location from tab_station where equipment=2 order by datetime;
# save as: "tab_station.csv"
In [11]:
samples$LAT <- as.numeric(as.character(substr(samples$latitude,1,2)) ) +
as.numeric(as.character(substr(samples$latitude,3,9999)) )/60
samples$LONG <- samples$longitude
for(i in 1:dim(samples)[1]){
if(samples$longitude[i] > 999){
samples$LONG[i] <- as.numeric(as.character(substr(samples$longitude[i],1,2)) ) +
as.numeric(as.character(substr(samples$longitude[i],3,9999)) )/60
}
else {
samples$LONG[i] <- as.numeric(as.character(substr(samples$longitude[i],1,1)) ) +
as.numeric(as.character(substr(samples$longitude[i],2,9999)) )/60
}
}
In [ ]:
# run sql: SELECT tvl.*, rkw.known_wrong, rkw.id_videospecies, rvs.id, rvs.valid_species_name
# FROM tab_video_logfiles tvl
# LEFT JOIN ref_known_wrong_videospecies rkw ON tvl.biology_observation ilike rkw.known_wrong
# LEFT JOIN ref_videospecies rvs ON tvl.biology_observation ilike rvs.valid_species_name OR rkw.id_videospecies = rvs.id
# WHERE tvl.biology_observation is not null AND tvl.originator = 'Biology'
# ORDER BY sample_no, date_time;
# save as: "tab_video_logfiles.csv"
In [6]:
output <- aggregate(!is.na(data$id),by=list(data$sample_no), sum)
colnames(output)<-c("sample_no", "count_valid_taxa")
output <- merge(output,samples, all.y=TRUE)
write.csv(output, "output1.csv", row.names=F)
In [8]:
coral_spp <- c("Alcyonacea", "Isidella lofotensis", "Paragorgia arborea", "Paramuricea placomus",
"Primnoa resedaeformis", "Radicipes sp.")
subdata <- data[data$valid_species_name %in% coral_spp,]
output <- aggregate(!is.na(subdata$id),by=list(subdata$sample_no), sum)
colnames(output)<-c("sample_no", "count_Hcorals")
output <- merge(output,samples, all.y=TRUE)
write.csv(output[!is.na(output$count_Hcorals),], "output2.csv", row.names=F)
In [7]:
one_spp <- c("Parastichopus tremulus")
subdata <- data[data$valid_species_name == one_spp,]
output <- aggregate(!is.na(subdata$id),by=list(subdata$sample_no), sum)
colnames(output)<-c("sample_no", "count_SppX")
output <- merge(output,samples, all.y=TRUE)
write.csv(output, "output4.csv", row.names=F)
In [4]:
library(vegan)
data.sel <- data[data$ship_lon_deg>12 & data$ship_lon_deg<16,]
sppdata <- data.sel[data.sel$valid_species_name!="",colnames(data.sel)==c("sample_no","valid_species_name")]
sppdata$Value <- 1
sppmatrix<-xtabs(sppdata$Value~sppdata$sample_no+sppdata$valid_species_name)
sppmatrix<-as.data.frame.matrix(sppmatrix)
#dca <- decorana(sppmatrix)
sedidata <- data.sel[data.sel$bottom_observation!="No observation",colnames(data.sel)%in%c("sample_no","bottom_observation")]
sedidata$Value <- 1
sedidata$bottom_observation <-gsub("Coral reef", "Coral_reef", sedidata$bottom_observation)
sedimatrix<-xtabs(sedidata$Value~sedidata$sample_no+sedidata$bottom_observation)
sedimatrix<-as.data.frame.matrix(sedimatrix)
nmds <- metaMDS(sppmatrix, trace = FALSE)
ef <- envfit(nmds, sedimatrix, permu = 99)
plot(nmds, display = "sites")
plot(ef, p.max = 0.1)
In [5]:
ef <- envfit(nmds ~ Coral_reef + Mud, sedimatrix)
plot(nmds, display = "sites")
plot(ef)
tmp <- with(sedimatrix, ordisurf(nmds, Mud, add = TRUE))
with(sedimatrix, ordisurf(nmds, Coral_reef + Mud, add = TRUE, col = "green4"))
In [12]:
library(raster)
minT <- raster("NorKyst800_Troms_stat_2010_minT.tif")
# convert samples to spatial data
coordinates(samples)<-~LONG+LAT
proj4string(samples) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
samples <- spTransform(samples,
CRS=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
v<-extract(minT, samples)
v<-cbind(v,samples@data$sample_no)
colnames(v)<- c("minT","sample_no")
v1<-merge(v,output)
In [13]:
head(v1)
In [16]:
In [16]:
# using least squares regression
library(nlme)
v1<-v1[,1:3]
v2<-v1[complete.cases(v1),]
model1<- gls(count_SppX~minT, v2)
model1
# test the significance of the term (first fit null model)
model0<-gls(count_SppX~1, v2)
model0.ml <- update(model0, . ~ ., method = "ML")
model1.ml <- update(model1, . ~ ., method = "ML")
anova(model0.ml, model1.ml)
# Make a plot -------------------------------------------------------------
plot(v2$count_SppX~v2$minT,
main="",
xlab="minimum Temp",
ylab="Count of Parastichopus tremulus")
lines(sort(fitted(model1))~sort(v2$minT))
In [19]:
# using gam
library(gam)
gam1 <- gam(count_SppX>0~s(minT),family=binomial)
In [ ]: