library(vegan)
#1
ryby.rich<-specnumber(doubs.spe)
names(doubs.env)

lm.0<-lm(log(ryby.rich)~+1, data=doubs.env)
as.formula(doubs.env)
lm.step<-step(lm.0, .~.+ alt + slo + flo + pH + har + log1p(phos) + nit + log1p(amm) + 
    oxy + log1p(bdo))
anova(lm.step)
lm.step.1<-update(lm.step, .~.-slo-pH)
summary(lm.step.1)
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  5.4133730  0.8019458   6.750 3.01e-07 ***
# alt         -0.0024826  0.0004071  -6.098 1.64e-06 ***
# log1p(bdo)  -0.4974455  0.1860609  -2.674   0.0126 *  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.5504 on 27 degrees of freedom
# Multiple R-squared:  0.5803,	Adjusted R-squared:  0.5492 
# F-statistic: 18.67 on 2 and 27 DF,  p-value: 8.124e-06

rda.0<-rda(log(ryby.rich)~1, data=doubs.env)
as.formula(doubs.env)
rda.step<-ordistep(rda.0, .~.+ alt + slo + flo + pH + har + log1p(phos) + nit + log1p(amm) + 
    oxy + log1p(bdo))
anova(rda.step, by="terms")
# Model: rda(formula = log(ryby.rich) ~ alt + log1p(bdo), data = doubs.env)
#            Df Variance       F Pr(>F)    
# alt         1 0.315279 30.1824  0.001 ***
# log1p(bdo)  1 0.074666  7.1479  0.014 *  
# Residual   27 0.282036  

#2
rhin.env<-data.frame(t(read.delim2("clipboard", row.names = 1)))
rhin.spe<-data.frame(t(read.delim2("clipboard", header=F, row.names = 1)))

rhin.spe<-rhin.spe[,colSums(rhin.spe>0)>1] # Odstranění 
# singulárních výskytů druhů

rda.1<-rda(rhin.spe~Rhinanthus.sown*Year+Mown.twice*Year+
             Condition(as.factor(Block)+Year+
                         Rhinanthus.sown+Mown.twice),
           data=rhin.env)#Specifikace modelu

anova(rda.1, by="terms", permutations = how(
  blocks=as.factor(rhin.env$Block), 
  plots = Plots(strata=as.factor(rhin.env$Plot), type="free"),
  within=Within(type="none"), nperm=999
))# Specifikace testu - důležité je nastavení bloků a plotů v how()

# Permutation test for rda under reduced model
# Terms added sequentially (first to last)
# Blocks:  as.factor(rhin.env$Block) 
# Plots: as.factor(rhin.env$Plot), plot permutation: free
# Permutation: none
# Number of permutations: 999
# 
# Model: rda(formula = rhin.spe ~ Rhinanthus.sown * Year + Mown.twice * Year + Condition(as.factor(Block) + Year + Rhinanthus.sown + Mown.twice), data = rhin.env)
#                      Df Variance      F Pr(>F)   
# Rhinanthus.sown:Year  1    39.90 5.6298  0.002 **
# Year:Mown.twice       1    11.76 1.6591  0.128   
# Residual             85   602.35                 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


rda.2<-rda(rhin.spe~Rhinanthus.sown*Year+
             Condition(as.factor(Block)+Year+
                         Rhinanthus.sown+Mown.twice),
           data=rhin.env)# Stejný model jako rda.1 jen bez neprůkazné interakce kosení a času
anova(rda.2)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = rhin.spe ~ Rhinanthus.sown * Year + Condition(as.factor(Block) + Year + Rhinanthus.sown + Mown.twice), data = rhin.env)
#          Df Variance     F Pr(>F)   
# Model     1    39.90 5.587  0.002 **
# Residual 86   614.11                
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#3
bru.6<-bru[bru$instar==6,]
summary(bru.6)

#PCA
pca.bru<-rda(bru.6[,4:11], scale=T)

#PCA diagram
aa<-ordiplot(pca.bru, display=c("si", "sp"), scaling=3, type="none")
ordispider(pca.bru, groups=bru.6$species, col=1:6)
text(aa, what="sp", arrows=T)

rda.1<-rda(bru.6[,4:11]~bru.6$species, scale=T)
anova(rda.1)

library(MASS)# je potřeba pro LDA
as.formula(bru.6)# pomůcka pro generování prediktorů oddělených +
# species ~ instar + sex + antseg1 + antseg2 + antseg3 + antseg4 + 
#     midfem + midtib + hindfem + hindtib

lda.1<-lda(species~antseg1 + antseg2 + antseg3 + antseg4 + 
    midfem + midtib + hindfem + hindtib, data=bru.6)
#Nafitování LDA

aa<-ordiplot(lda.1, display=c("sp", "si"), type="n")
ordispider(lda.1, groups=bru.6$species, col=1:6)
text(aa, what="sp", arrows = T)
text(aa, what="sp", arrows = F)#V LDA druhy nelze zobrazit
#šipkou i popiskem zároveň (nekompatibilita MASS a veganu)
# Je to třeba rozdělit do dvou kroků


# Budování diskriminační funkce pomocí reverzní CCA
library(fastDummies)

cca.bru<-cca(dummy_cols(bru.6$species,
                        remove_selected_columns = T)~1, 
             data=bru.6)
# CCA potřebuje převést faktor odpověď na dummy proměnné

# Postupný výběr prediktorů
cca.step<-ordistep(cca.bru, scope=~antseg1 + antseg2 + antseg3 + antseg4 + 
    midfem + midtib + hindfem + hindtib)

anova(cca.step, by="terms")
#          Df ChiSquare       F Pr(>F)    
# hindtib    1   0.90412 42.4882  0.001 ***
# midtib     1   0.55496 26.0800  0.001 ***
# antseg4    1   0.37192 17.4781  0.001 ***
# hindfem    1   0.32815 15.4212  0.001 ***
# antseg3    1   0.25414 11.9432  0.001 ***
# antseg1    1   0.09575  4.4998  0.002 ** 
# midfem     1   0.08133  3.8220  0.001 ***
# antseg2    1   0.04760  2.2371  0.044 *  
# Residual 111   2.36201    
cca.step<-update(cca.step, .~.-antseg2)# Odstranění prediktoru na hranici signifikance

lda.cv<-lda(species~hindtib + midtib +
antseg4 + hindfem + antseg3 + antseg1 + midfem, data=bru.6,
CV=T)# Klasifikační LDA s krosvalidací
lda.cv

table(lda.cv$class==bru.6$species)# Tabulka správně a špatně klasifikovaných jedinců
# FALSE  TRUE 
#    17   103