setwd("/SET PATH TO DIRECTORY")
require(ade4)
require(FD)
require(cluster)

t <- read.csv("abiotic.csv", h = T)
head(t)
sea <- factor(t$Season)
tab <- t[,5:ncol(t)]

#PCA
pca <- dudi.pca(tab, scan = F, nf = 2)
s.corcircle(pca$co)
s.class(pca$li, sea)

#Mapping PCA axis scores
trans <- read.csv("Transect_INFO.csv", h = T, row.names = 1)
head(trans)
df <- cbind(t$Dis, t$Dep)
par(mfrow = c(2, 1))
plot(trans$dis, trans$dep, type = "l")
s.value(df, pca$li[,1], csi = 0.6, cleg = 0, add.p = T)
#Does seasonality depend on depth?
plot(trans$dis, trans$dep, type = "l")
s.value(df, pca$li[,2], csi = 0.6, cleg = 0, add.p = T)

#Clustering stations on the first axis
d <- gowdis(pca$li[1])
cl <- hclust(d, "ward.D2")
plot(cl, hang = -0.1)
g <- factor(cutree(cl, 4))
hab <- paste(rep("H", length(g)), g, sep = "")
hab <- factor(hab)
s.class(pca$li, hab)

par(mfrow = c(2, 1))
plot(trans$dis, trans$dep, type = "l")
s.value(df, pca$li[,1], csi = 0.6, cleg = 0, add.p = T)
plot(trans$dis, trans$dep, type = "l")
text(df, hab, font = 2)

#Masking seasonal effect
w <- wca(pca, sea, scan = F, nf = 2)
d <- gowdis(w$li[1])
cl <- hclust(d, "ward.D2")
plot(cl, hang = -0.1)
g <- factor(cutree(cl, 4))
hab <- paste(rep("H", length(g)), g, sep = "")
hab <- factor(hab)
levels(hab) <- c("L", "T", "C", "S")
s.class(w$li, hab)

par(mfrow = c(3, 1), mar = c(2, 2, 1, 1))
hab2 <- hab
levels(hab2) <- c(3, 4, 2, 1)
hab2 <- factor(as.character(hab2))
plot(trans$dis, trans$dep, type = "l")
s.value(df, w$li[,1], csi = 0.6, cleg = 0, add.p = T)
boxplot(t$Dep ~ hab2, col = 8, names = c("Stoss", "Crest", "Lee", "Through"))
plot(trans$dis, trans$dep, type = "l")
text(df, as.character(hab), font = 2, cex = 2)

