library(vegan)

library(corrplot)

library(berryFunctions)

library(packfor) ##alternatively adespatial can be used


##PAIRWISE CORRELATIONS BETWEEN ENV.TAB VARIABLES----------------


plot(env.tab$Altitude, env.tab$TRI)

plot(env.tab[,1:2], pch=16, col=densCols(env.tab[,1:2], colramp= colorRampPalette(c("gray80",
"gray20"))))


cor.test(env.tab$Altitude, env.tab$TRI)


cor(env.tab)


#plot correlation matrix

source("panelutils.R")

op <- par(mfrow=c(1,1), pty="s")

pairs(env.tab[1:4], lower.panel=panel.smooth, upper.panel=panel.cor, diag.panel=panel.hist,
main="Pearson correlation matrix")

par(op)


#another correlation matrix

cor.mat <- cor(env.tab)

head(cor.mat)


corrplot.mixed(cor.mat, lower = "number", upper = "ellipse", tl.pos = "lt",

               tl.cex=0.7, tl.col = "black", number.cex=0.7,

               title="Correlation matrix")



###EXPLORE RESPONSE VARIABLE----------------------------------

#map species richness

coord <- forest[, c("X", "Y")]


plot(coord$X, coord$Y)


colPoints(coord$X, coord$Y, forest$Diversity, add = F)


#specify my own color scale

my.col <- colorRampPalette(c("blue", "yellow", "red"))


colPoints(coord$X, coord$Y, forest$Diversity, add = F, col=my.col(100))


#statistics

summary(forest$Diversity)

boxplot(forest$Diversity)

hist(forest$Diversity)

plot(density(forest$Diversity))


hist(log(forest$Diversity))


#test normality

shapiro.test(forest$Diversity)

shapiro.test(log(forest$Diversity))

## N0 = výběr pochází normálně rozloženého souboru


###CORRELATION ANALYSIS OF DIVERSITY AND ENV.TAB VARIABLES-----------------

#correlate diversity and enviro variables

cor(forest$Diversity, env.tab)


plot(env.tab$Limestone, forest$Diversity, main="Limestone")

plot(env.tab$Arable.land, forest$Diversity, main="Arable land")


plot(env.tab$ForestAB, forest$Diversity, main="Natural forests")

cor(env.tab$ForestAB, forest$Diversity)



#loess smooth

l <- loess.smooth(env.tab$ForestAB, forest$Diversity)

lines(l, lwd=3, col="red")


##LINEAR REGRESSION-------------------------------------------

m <- lm(forest$Diversity ~ ForestAB, data=env.tab.s)

anova(m)

summary(m)


m <- lm(forest$Diversity ~ ForestAB + Altitude, data=env.tab.s)

anova(m)

summary(m)



m <- lm(forest$Diversity ~ ., data=env.tab.s)

m


anova(m)

summary(m)


m2 <- lm(forest$Diversity ~ ., data=env.tab.s[, sample(1:15)])

anova(m2)

summary(m2)



library("packfor")


####FORWARD SELECTION-------------------------------------------------


#forward selection based on significance level

fs <- forward.sel(forest$Diversity, env.tab.s, alpha=0.05)

fs


m <- lm(forest$Diversity ~ Limestone + Arable.land + Altitude, data=env.tab.s)


anova(m)

summary(m)


res <- resid(m)

pr <- predict(m)


#compare predicted and observed values

head(res)

head(pr)

head(forest$Diversity)


##calculate total variablity in response variable

TSS <- sum((mean(forest$Diversity) - forest$Diversity)^2)


#extract residuals

RSS <- sum(resid(m)^2)


#calculate R2

1-(RSS/TSS)



#evaluate model

# plot(m)

hist(res)


#map residuals

forest$res <- res


colPoints(coord$X, coord$Y, res, add=F, col=my.col(100), main="Residuals")


#map fitted values

forest$pr <- pr


colPoints(coord$X, coord$Y, pr, add=F, col=my.col(100), main="Fitted")


###VARIATION PARTITIONING-------------------------------------

library(vegan)

colnames(env.tab)


v <- varpart(forest$Diversity, env.tab.s[,5:6], env.tab.s[, 7:9], env.tab.s[,10:15])

plot(v)


showvarparts(3)