#### 1 - linearni model pomoci permutaci
# Data z řeky Doubs
# druhova data jsou v d.spe
# promenne prostredi v d.env

library(vegan)
sp.rich<-specnumber(d.spe)
lm.rich.0<-lm(sp.rich ~ 1, data=d.env)
lm.step<-step(lm.rich.0, scope=~ dfs + alt + slo + flo + pH + har + phos + nit + amm + 
    oxy + bdo)
anova(lm.step)
# Response: sp.rich
#           Df  Sum Sq Mean Sq F value    Pr(>F)    
# alt        1 1017.17 1017.17 93.6541 2.184e-09 ***
# bdo        1  338.37  338.37 31.1546 1.305e-05 ***
# amm        1  245.79  245.79 22.6310 9.500e-05 ***
# phos       1   47.04   47.04  4.3312   0.04927 *  
# pH         1   70.92   70.92  6.5301   0.01804 *  
# har        1   35.80   35.80  3.2963   0.08309 .  
# nit        1   35.47   35.47  3.2661   0.08442 .  
# Residuals 22  238.94   10.86   
lm.step<-update(lm.step, .~.-har-nit)# Jeste vyradime prediktory, ktere nejsou prukazne
# i kdyz prosly pres AIC vyber
summary(lm.step)
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 97.508701  33.579732   2.904 0.007790 ** 
# alt         -0.015433   0.003719  -4.150 0.000361 ***
# bdo         -0.212005   0.029667  -7.146 2.19e-07 ***
# amm          2.316272   0.669092   3.462 0.002025 ** 
# phos         2.839591   1.057043   2.686 0.012904 *  
# pH          -1.000205   0.426993  -2.342 0.027784 *  
# 
# Residual standard error: 3.595 on 24 degrees of freedom
# Multiple R-squared:  0.8471,	Adjusted R-squared:  0.8153 
# F-statistic:  26.6 on 5 and 24 DF,  p-value: 4.696e-09


rda.rich.0<-rda(sp.rich ~ 1, data=d.env)
rda.rich.f<-rda(sp.rich ~ dfs + alt + slo + flo + pH + har + phos + nit + 
    amm + oxy + bdo, data=d.env)
rda.step<-ordistep(rda.rich.0, rda.rich.f)

rda.step
# Call: rda(formula = sp.rich ~ alt + bdo + amm, data = d.env)
# 
#               Inertia Proportion Rank
# Total          69.983      1.000     
# Constrained    55.218      0.789    1
# Unconstrained  14.765      0.211    1
# Inertia is variance 
# 
# Eigenvalues for constrained axes:
#  RDA1 
# 55.22 

anova(rda.step, by="terms")
# Permutation test for rda under reduced model
# Terms added sequentially (first to last)
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = sp.rich ~ alt + bdo + amm, data = d.env)
#          Df Variance      F Pr(>F)    
# alt       1   35.075 61.765  0.001 ***
# bdo       1   11.668 20.547  0.001 ***
# amm       1    8.476 14.925  0.006 ** 
# Residual 26   14.765   

RsquareAdj(rda.step)
# $r.squared
# [1] 0.789024
# 
# $adj.r.squared
# [1] 0.7646806

### Vyber pomoci RDA vychazi podobne. Prediktory na hrane signifikance
# ale neprojdou. R2 vychazi velmi podobne (kdyby lm neobsahoval o 
# jeden prediktor navic, vychazelo by uplne stejne).


##### Rhinanthus

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

to.rem<-colSums(rhin.spe>0)<2 #odstraneni druhu, ktere se vyskytuji pouze jedenkrat
rhin.spe<-rhin.spe[,-to.rem]

rhin.env$Year<-rhin.env$Year-2012#Preskalovani roku na 0-3
names(rhin.env)

# Specifikace modelu
rda.1<-rda(log(rhin.spe+1)~
             Rhinanthus.sown*Year+Mown.twice*Year+
             Condition(Rhinanthus.sown+Year+Mown.twice+
                         Block), data=rhin.env)
#specifikace permutacniho testu pomoci how:
# zasadni je specifikace plotu (whole-plotu) a jejich permutace
# a vypnuti permutace split-plotu
anova(rda.1, permutations = how(blocks = rhin.env$Block,
                               plots=Plots(strata=rhin.env$Plot,
                                   type="free"),
                               within=Within("none"),nperm=999), by="terms")

# Permutation test for rda under reduced model
# Terms added sequentially (first to last)
# Blocks:  rhin.env$Block 
# Plots: rhin.env$Plot, plot permutation: free
# Permutation: none
# Number of permutations: 999
# 
# Model: rda(formula = log(rhin.spe + 1) ~ Rhinanthus.sown * Year + Mown.twice * Year + Condition(Rhinanthus.sown + Year + Mown.twice + Block), data = rhin.env)
#                      Df Variance      F Pr(>F)    
# Rhinanthus.sown:Year  1   0.1256 1.0288  0.001 ***
# Year:Mown.twice       1   0.0776 0.6353  0.136    
# Residual             89  10.8675                  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


rda.2<-update(rda.1,.~.-Year:Mown.twice)## Jesne z modelu odstranime
# neprukaznou interakci casu a koseni dvakrat

#Finalni test modelu
anova(rda.2, permutations = how(blocks = rhin.env$Block,
                               plots=Plots(strata=rhin.env$Plot,
                                           type="free"),
                               within=Within("none"),nperm=999), by="terms")

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


### 3 Bruslarky
# Data jsou ulozena v df bru
bru.6<-bru[bru$instar==6,]#Selekce jedincu 6. instaru

pca.1<-rda(bru.6[,4:11], scale=T)# PCA se standardizaci dat
pca.1
# Call: rda(X = bru.6[, 4:11], scale = T)
# 
#               Inertia Rank
# Total               8     
# Unconstrained       8    8
# Inertia is correlations 
# 
# Eigenvalues for unconstrained axes:
#   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8 
# 6.831 0.790 0.203 0.061 0.058 0.027 0.020 0.011
screeplot(pca.1, bstick = T)# Screeplot s BS modelem

aa<-ordiplot(pca.1, type='n') # prazdny ordinacni diagram
text(aa, "sp", arrows = T, length=0.05) # pridani promennych ("druhu") do diagramu
# points(aa, "si", pch=".")# Zobrazeni vzorku teckami
# pridani obalek
ordispider(pca.1, groups=bru.6$species, # definice skupin
         col=2:7# specifikace barev 
       )
legend("topleft", fill=2:7, legend=levels(bru.6$species), 
       ncol=1, title="Species") # pridani legendy. Pro podrobnosti mrknete na ?legend

# otestovani pomoci RDA

rda.1<-rda(bru.6[,4:11]~bru.6$species, scale=T)
anova(rda.1)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = bru.6[, 4:11] ~ bru.6$species, scale = T)
#           Df Variance      F Pr(>F)    
# Model      5    6.352 87.881  0.001 ***
# Residual 114    1.648    

# Diskriminacni analyza
?lda
library(MASS)

lda.1<-lda(bru.6$species~antseg1 + antseg2 + antseg3 + antseg4 + midfem + midtib + hindfem + 
    hindtib, data=bru.6)

lda.1
# Call:
# lda(bru.6$species ~ antseg1 + antseg2 + antseg3 + antseg4 + midfem + 
    # midtib + hindfem + hindtib, data = bru.6)
# 
# Prior probabilities of groups:
#    canali    dissor    esakii     genit     notab      rufo 
# 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 
# 
# Group means:
#        antseg1 antseg2 antseg3 antseg4  midfem midtib hindfem hindtib
# canali  1.6750   1.130  1.0075  1.7540  5.4655 4.7855  5.9585  2.7160
# dissor  2.2655   1.668  1.2255  1.9735  8.8515 6.9790 10.4520  6.3635
# esakii  1.2390   0.915  0.8735  1.6270  3.9375 3.1845  4.2210  2.3860
# genit   2.2070   1.537  1.1055  1.7265  8.2565 6.5970 10.1785  6.1400
# notab   2.6310   1.918  1.3930  1.7770 10.1865 8.2175 11.5330  7.2900
# rufo    2.2790   1.694  1.2960  1.9385  8.3475 6.5755  9.9825  6.0275
# 
# Coefficients of linear discriminants:
#                LD1        LD2        LD3        LD4        LD5
# antseg1  0.6984262  4.2244975 -0.2882184 -4.2331457 -4.4952166
# antseg2 -0.1025288 -1.0157953  2.3909501  5.0926536  1.0403471
# antseg3 -5.3351297  0.4992989  1.8288249  6.1008294 -6.1111943
# antseg4  1.7987170 -1.4643478 -5.4587061  3.3955872  2.8454363
# midfem   0.7629668  1.0253636 -0.4539265  0.7160885  1.2935379
# midtib   0.4600042  2.6295965  0.5754982 -0.6621269  1.1138436
# hindfem -0.5072149 -2.0439907 -1.8292900 -0.7020821 -1.1301452
# hindtib  1.4956696 -1.3036616  2.2903072  0.1903232  0.4880557
# 
# Proportion of trace:
#    LD1    LD2    LD3    LD4    LD5 
# 0.7443 0.1482 0.0552 0.0424 0.0100 

bb<-ordiplot(lda.1, type='n') # prazdny ordinacni diagram

ordispider(lda.1, groups=bru.6$species, # definice skupin
         col=2:7)# specifikace barev 

#Z nejakeho duvodu nefunguje pridani sipek a nazvu druhu v jednom radku (asi nekompatibilita mezi veganem a MASS), tak je treba to 
# udelat samostatne.
text(bb, what="species", arrows = T, length=0.05) # pridani promennych ("druhu") do diagramu
text(bb, what="species", arrows = F) # pridani promennych ("druhu") do diagramu
# points(aa, "si", pch=".")# Zobrazeni vzorku teckami

legend("bottom", fill=2:7, legend=levels(bru.6$species), 
       ncol=6, title="Species") # pridani legendy. Pro podrobnosti mrknete na ?legend


## Sestaveni diskriminacni fuknce
# Potrebujeme postupny vyber, ten ale neni pro fuknci LDA k dispozici, takze vyuzijeme
# jeji pribuznost s CCA, pro kterou je postupny vyber implementovan

# CCA ovsem nepojme faktor jako odpoved, je proto treba prevest faktor na serii
# binarnich promennych (dummy variables)

install.packages("fastDummies")
library(fastDummies)
dummy_cols(bru.6$species)


cca.0<-cca(dummy_cols(bru.6$species)[,-1]~1, data=bru.6)
cca.step<-ordistep(cca.0, scope = ~antseg1 + antseg2 + antseg3 + antseg4 + midfem + midtib + hindfem + 
    hindtib)
##Postupny vyber vybral svechny parametry krome antseg 2
cca.step

# Call: cca(formula = dummy_cols(bru.6$species)[, -1] ~ hindtib + midtib + antseg4 +
# hindfem + antseg3 + midfem + antseg1, data = bru.6)
# 
#               Inertia Proportion Rank
# Total          5.0000     1.0000     
# Constrained    2.5904     0.5181    5
# Unconstrained  2.4096     0.4819    5
# Inertia is scaled Chi-square 
# 
# Eigenvalues for constrained axes:
#   CCA1   CCA2   CCA3   CCA4   CCA5 
# 0.9227 0.7026 0.4642 0.3642 0.1368 
# 
# Eigenvalues for unconstrained axes:
#    CA1    CA2    CA3    CA4    CA5 
# 0.8632 0.6358 0.5358 0.2974 0.0773 

# Rozdily mezi druhy odpovidaji 50% variability

anova(cca.step, by="margin")#Otestovani parcialnich efektu parametru
Model: cca(formula = dummy_cols(bru.6$species)[, -1] ~ hindtib + midtib + antseg4 + hindfem + antseg3 + midfem + antseg1, data = bru.6)
#           Df ChiSquare       F Pr(>F)    
# hindtib    1   0.29722 13.8147  0.001 ***
# midtib     1   0.21850 10.1559  0.001 ***
# antseg4    1   0.27792 12.9180  0.001 ***
# hindfem    1   0.29363 13.6483  0.001 ***
# antseg3    1   0.19664  9.1398  0.001 ***
# midfem     1   0.08133  3.7802  0.001 ***
# antseg1    1   0.08134  3.7805  0.004 ** 
# Residual 112   2.40961                   
  
### Kros-validace diskriminacni funkce
lda.2<-lda(bru.6$species~antseg1 + antseg3 + antseg4 + midfem + midtib + hindfem + 
    hindtib, data=bru.6, CV=T)# Fitovani LDA s klasifikacni funkci
lda.2$class #obsahuje klasifikaci vzorku pomoci funkce

table(bru.6$species==lda.2$class)# Porovnani mezi skutecnou druhovou prislusnosti a 
#klasifikaci na zaklade diskriminacni funkce
# FALSE  TRUE 
#    17   103 

# Proporce spravne urcenych jedincu:
1-17/(17+103)
# [1] 0.8583333
