doubs.env<-read.delim2("doubs_env.csv", sep=";")
doubs.spe<-read.delim2("doubs_spe.csv", sep=";")

summary(doubs.env)
summary(doubs.spe)

library(vegan)

sha<-diversity(doubs.spe)
sp.rich<-specnumber(doubs.spe)
eve<-sha/log(sp.rich)

df.div<-data.frame(sha, sp.rich, eve)

library(psych)
pdf("doubs.pairs.pdf", 10,10)
pairs.panels(data.frame(doubs.env[,-1], df.div))
dev.off()

##
lm.0<-lm(sp.rich~+1, data=doubs.env)
names(doubs.env)
aa<-add1(lm.0, .~+alt+slo+pH+har+phos+nit+oxy+bdo, test="F")
# Model:
# sp.rich ~ +1
#        Df Sum of Sq    RSS    AIC F value    Pr(>F)    
# <none>              2029.5 128.43                      
# alt     1   1017.17 1012.3 109.56 28.1336 1.208e-05 ***
# slo     1    785.43 1244.1 115.75 17.6773 0.0002423 ***
# pH      1      1.79 2027.7 130.40  0.0247 0.8763483    
# har     1    335.88 1693.6 125.00  5.5529 0.0256792 *  
# phos    1      0.82 2028.7 130.42  0.0114 0.9158196    
# nit     1    409.17 1620.3 123.68  7.0706 0.0128114 *  
# oxy     1     23.68 2005.8 130.08  0.3306 0.5698964    
# bdo     1     42.75 1986.8 129.79  0.6025 0.4441517

cor(sp.rich, y=doubs.env)


#Uprava p hodnot Holmovou korekci
p.adjust(aa$`Pr(>F)`, method="holm")
# [1]           NA 9.665588e-05 1.696321e-03 1.000000e+00
# [5] 1.283960e-01 1.000000e+00 7.686811e-02 1.000000e+00
# [9] 1.000000e+00

lm.1<-update(lm.0, .~.+alt)
add1(lm.1, .~+alt+slo+pH+har+phos+nit+oxy+bdo, test="F")
lm.2<-update(lm.1, .~.+bdo)
add1(lm.2, .~+alt+slo+pH+har+phos+nit+oxy+bdo, test="F")
lm.3<-update(lm.2, .~.+oxy)
add1(lm.3, .~+alt+slo+pH+har+phos+nit+oxy+bdo, test="F")

anova(lm.3)
# Analysis of Variance Table
# 
# Response: sp.rich
#           Df  Sum Sq Mean Sq F value    Pr(>F)    
# alt        1 1017.17 1017.17  55.162 6.921e-08 ***
# bdo        1  338.37  338.37  18.350 0.0002224 ***
# oxy        1  194.54  194.54  10.550 0.0031970 ** 
# Residuals 26  479.43   18.44                      

plot(lm.3)

summary(lm.3)
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -6.5705 -3.3531  0.1593  2.2683  8.6743 
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 56.147340   8.078829   6.950 2.23e-07 ***
# alt         -0.024779   0.003159  -7.845 2.55e-08 ***
# bdo         -0.196552   0.038459  -5.111 2.51e-05 ***
# oxy         -0.220016   0.067737  -3.248   0.0032 ** 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 4.294 on 26 degrees of freedom
# Multiple R-squared:  0.7638,	Adjusted R-squared:  0.7365 
# F-statistic: 28.02 on 3 and 26 DF,  p-value: 2.638e-08