---
title: "Master thesis plots and analyses"
author: "Gaute Eiterjord"
date: "2023-12-15"
output:
  html_document: default
  word_document: default
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

This document has the scripts for all the analyses (and some more) that I performed as part of my master thesis in Natural Resource Management in 2023/2024. The datasets will hopefully be uploaded unto NMBUs webserver, if not it can be retrieved by sending an email to gaute.eiterjord@gmail.com 

## Loading data and packages

I start of by loading necessary packages and the dataset

```{r echo=T, results='hide',message=FALSE,warning=FALSE}
library(ggplot2)
library(plyr)
library(dplyr)
library(patchwork)
library(readxl)
library(ggsignif)
library(vegan)
library(knitr)
library(lmerTest)
library(reshape2)
library(scales)
coverage<-read_excel("env_data.xlsx")
env_data<-read_excel("env_data.xlsx")

```

# Cover

### Vizualising data

Making boxplots of the cover of functional groups and litter

```{r echo=T, results='hide',fig.show='hide'}
#1.1 Shrub
dwarf_shrub_box <- coverage %>% 
  ggplot(aes(x=Type,y=Shrub,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Dwarf shrub cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.2 Dryas
dryas_box <- coverage %>% 
  ggplot(aes(x=Type,y=Dryas,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Dryas cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.3 Other shrub
other_shrub_box <- coverage %>% 
  ggplot(aes(x=Type,y=Other_shrub,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Other cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.4 Forb
forb_box <- coverage %>% 
  ggplot(aes(x=Type,y=Forb,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Forb cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.5 Graminiods
gram_box <- coverage %>% 
  ggplot(aes(x=Type,y=Gram,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Graminoid cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.6 Bryophytes
bryo_box <- coverage %>% 
  ggplot(aes(x=Type,y=Bryo,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Bryophyte cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.7 Lichens
lich_box <- coverage %>% 
  ggplot(aes(x=Type,y=Lich,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Lichen cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.8 Dead Dryas
dead_dryas_box <- coverage %>% 
  ggplot(aes(x=Type,y=Dead_Dryas,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Dead Dryas cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.9 Litter
litter_box <- coverage %>% 
  ggplot(aes(x=Type,y=Litter,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Litter cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.10 Bare soil
bare_soil_box <- coverage %>% 
  ggplot(aes(x=Type,y=Bare_soil,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Bare soil cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.11 Biological crust
biol_crust_box <- coverage %>% 
  ggplot(aes(x=Type,y=Biol_crust,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Biological crust cover (%)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

```

Plot for Dryas and dead Dryas Cover

```{r,fig.show='hide',warning=FALSE}
wrap_plots(dryas_box,dead_dryas_box,ncol=2,nrow=1)
```

For dwarf shrub and litter

```{r,fig.show='hide',warning=FALSE}
wrap_plots(dwarf_shrub_box,litter_box,ncol=2,nrow=1)
```

For forbs and graminoids

```{r,fig.show='hide',warning=FALSE}
wrap_plots(gram_box,forb_box,ncol=2,nrow=1)
```

For lichen and bryophyte

```{r,fig.show='hide',warning=FALSE}
wrap_plots(bryo_box,lich_box,ncol=2,nrow=1)
```

Based on the boxplots, it looks like these are not different: Dryas,
Dwarf shrub, graminoid and forb While these look different: Dead Dryas,
litter, lichen and bryophyte

### Checking assumptions

The next steps are checking two of the assumptions for using a T-test:
Normally distributed data and equal variance. I will do this by
inspecting histograms, Shapiro-wilk tests and F-tests.

Histograms for checking the distribution

```{r echo=T, results='hide',warning=FALSE}
#2 Histograms for checking if data is normally distributed
#2.1 Shrub
shrub_hist<-coverage %>% 
  ggplot(aes(x=Shrub,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("Dwarf shrub cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.2 Dryas
dryas_hist<-coverage %>% 
  ggplot(aes(x=Dryas,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Dryas cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.3 Other shrub
other_shrub_hist<-coverage %>% 
  ggplot(aes(x=Other_shrub,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Other shrub cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.4 Forb
forb_hist<-coverage %>% 
  ggplot(aes(x=Forb,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Forb cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.5 Gram
gram_hist<-coverage %>% 
  ggplot(aes(x=Gram,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Graminoid cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.6 Lich
lich_hist<-coverage %>% 
  ggplot(aes(x=Lich,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Lichen cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.7 Bryo
bryo_hist<-coverage %>% 
  ggplot(aes(x=Bryo,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Bryophyte cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.8 Dead Dryas
dead_dryas_hist<-coverage %>% 
  ggplot(aes(x=Dead_Dryas,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Dead Dryas cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.9 Litter
litter_hist<-coverage %>% 
  ggplot(aes(x=Litter,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Litter cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.10 Bare soil
bare_soil_hist<-coverage %>% 
  ggplot(aes(x=Bare_soil,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Bare soil cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.11 Biological crust
biol_crust_hist<-coverage %>% 
  ggplot(aes(x=Biol_crust,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Biological crust cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.12 organizing the plots in the same image

wrap_plots(shrub_hist,forb_hist,
           gram_hist,lich_hist,
           ncol=2,nrow=2)

wrap_plots(bryo_hist,litter_hist,
           dryas_hist,dead_dryas_hist,
           ncol=2,nrow=2)
```

The following looks normally distributed: Shrub, Dryas The following
looks left skewed: Forb, Gram, Lich, Bryo, Litter, Dead Dryas, Bare soil

Testing the assumptions of normally distributed data and equal variance
with shapiro-wilk and F-tests. In both tests the p-value should be over
0,05 for the test to pass

```{r echo=T, results='hide'}
#3.1shapiro-wilk-test to check normality distribution
#first filtering between C and OTC populations
coverage_c<-coverage[coverage$Type=="C",]
#View(coverage_c)
shapiro.test(coverage_c$Shrub)
shapiro.test(coverage_c$Dryas)
shapiro.test(coverage_c$Other_shrub)
shapiro.test(coverage_c$Forb)
shapiro.test(coverage_c$Gram)
shapiro.test(coverage_c$Lich)
shapiro.test(coverage_c$Bryo)
shapiro.test(coverage_c$Litter)
shapiro.test(coverage_c$Dead_Dryas)
shapiro.test(coverage_c$Bare_soil)
shapiro.test(coverage_c$Biol_crust)

coverage_o<-coverage[coverage$Type=="OTC",]
#View(coverage_o)
shapiro.test(coverage_o$Shrub)
shapiro.test(coverage_o$Dryas)
shapiro.test(coverage_o$Other_shrub)
shapiro.test(coverage_o$Forb)
shapiro.test(coverage_o$Gram)
shapiro.test(coverage_o$Lich)
shapiro.test(coverage_o$Bryo)
shapiro.test(coverage_o$Litter)
shapiro.test(coverage_o$Dead_Dryas)
shapiro.test(coverage_o$Bare_soil)
shapiro.test(coverage_o$Biol_crust)

#3.2 Checking equal variances with F test
var.test(Shrub~Type,data=coverage)
var.test(Dryas~Type,data=coverage)
var.test(Other_shrub~Type,data=coverage)
var.test(Forb~Type,data=coverage)
var.test(Gram~Type,data=coverage)
var.test(Lich~Type,data=coverage)
var.test(Bryo~Type,data=coverage)
var.test(Dead_Dryas~Type,data=coverage)
var.test(Litter~Type,data=coverage)
var.test(Bare_soil~Type,data=coverage)
var.test(Biol_crust~Type,data=coverage)
```

The results are in, only dwarf shrub and Dryas has normal distribution,
while the following has equal variance: Dwarf shrub, Dryas, forb. The
others have not.

### Transforming data

```{r echo=T, results='hide',message=FALSE,warning=FALSE, fig.show='hide'}
#4.1 Transforming data
#4.1.1 Transformation of Gram
gram_coverage <- coverage[,c("Type","Gram")]
#View(gram_coverage)
#log10
gram_coverage$gram_log10 <- log10(gram_coverage$Gram)
#ln
gram_coverage$gram_ln <- log(gram_coverage$Gram)
#square root
gram_coverage$gram_square <- sqrt(gram_coverage$Gram)
#cube root
gram_coverage$gram_cube <- (sign(gram_coverage$Gram)*abs(gram_coverage$Gram)^(1/3))

#histograms
#untransformed
gram_coverage %>% 
  ggplot(aes(x=Gram,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Graminoid cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#log10
gram_coverage %>% 
  ggplot(aes(x=gram_log10,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Graminoid cover (log10)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#ln
gram_coverage %>% 
  ggplot(aes(x=gram_ln,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Graminoid cover (ln)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#square root
gram_coverage %>% 
  ggplot(aes(x=gram_square,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Graminoid cover (square root)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#cube root
gram_coverage %>% 
  ggplot(aes(x=gram_cube,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Graminoid cover (cube root)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#4.1.2 Shapiro-wilk-test of transformed data
gram_coverage_C <- gram_coverage[gram_coverage$Type == "C",]
#View(gram_coverage_C)

shapiro.test(gram_coverage_C$Gram)
shapiro.test(gram_coverage_C$gram_log10)
shapiro.test(gram_coverage_C$gram_ln)
shapiro.test(gram_coverage_C$gram_square)
shapiro.test(gram_coverage_C$gram_cube)

gram_coverage_OTC <- gram_coverage[gram_coverage$Type == "OTC",]
#View(gram_coverage_OTC)

shapiro.test(gram_coverage_OTC$Gram)
shapiro.test(gram_coverage_OTC$gram_log10)
shapiro.test(gram_coverage_OTC$gram_ln)
shapiro.test(gram_coverage_OTC$gram_square)
shapiro.test(gram_coverage_OTC$gram_cube)

#Cube transformation is good
var.test(gram_cube~Type,data=gram_coverage)

#4.2.1 Transformation of Bryo
bryo_coverage <- coverage[,c("Type","Bryo")]
#View(bryo_coverage)
#log10
bryo_coverage$bryo_log10 <- log10(bryo_coverage$Bryo+1)
#ln
bryo_coverage$bryo_ln <- log(bryo_coverage$Bryo+1)
#square root
bryo_coverage$bryo_square <- sqrt(bryo_coverage$Bryo)
#cube root
bryo_coverage$bryo_cube <- (sign(bryo_coverage$Bryo)*abs(bryo_coverage$Bryo)^(1/3))

#histograms
#untransformed
bryo_coverage %>% 
  ggplot(aes(x=Bryo,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Bryophyte cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#log10
bryo_coverage %>% 
  ggplot(aes(x=bryo_log10,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Bryophyte cover (log 10)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#ln
bryo_coverage %>% 
  ggplot(aes(x=bryo_ln,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Bryophyte cover (ln)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#square
bryo_coverage %>% 
  ggplot(aes(x=bryo_square,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Bryophyte cover (square root)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#cube
bryo_coverage %>% 
  ggplot(aes(x=bryo_cube,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Bryophyte cover (cube root)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#4.2.2 Shapiro-wilk-test of transformed data
bryo_coverage_C <- bryo_coverage[bryo_coverage$Type == "C",]
#View(bryo_coverage_C)

shapiro.test(bryo_coverage_C$Bryo)
shapiro.test(bryo_coverage_C$bryo_log10)
shapiro.test(bryo_coverage_C$bryo_ln)
shapiro.test(bryo_coverage_C$bryo_square)
shapiro.test(bryo_coverage_C$bryo_cube)

bryo_coverage_OTC <- bryo_coverage[bryo_coverage$Type == "OTC",]
#View(bryo_coverage_OTC)

shapiro.test(bryo_coverage_OTC$Bryo)
shapiro.test(bryo_coverage_OTC$bryo_log10)
shapiro.test(bryo_coverage_OTC$bryo_ln)
shapiro.test(bryo_coverage_OTC$bryo_square)
shapiro.test(bryo_coverage_OTC$bryo_cube)

#4.2.3 F-test to check equal variance
var.test(bryo_cube~Type,data=bryo_coverage)

#Cube-transformed data for bryophyte cover will be used as it is normally distributed and has equal variance

#4.3.1 Dead Dryas transformation
DD_coverage <- coverage[,c("Type","Dead_Dryas")]
#View(DD_coverage)
#log10
DD_coverage$DD_log10 <- log10(DD_coverage$Dead_Dryas+1)
#ln
DD_coverage$DD_ln <- log(DD_coverage$Dead_Dryas+1)
#square
DD_coverage$DD_square <- sqrt(DD_coverage$Dead_Dryas)
#cube
DD_coverage$DD_cube <- (sign(DD_coverage$Dead_Dryas)*abs(DD_coverage$Dead_Dryas)^(1/3))

#histograms
#untransformed
DD_coverage %>% 
  ggplot(aes(x=Dead_Dryas,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("DD cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#log10
DD_coverage %>% 
  ggplot(aes(x=DD_log10,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("DD cover (log 10)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#ln
DD_coverage %>% 
  ggplot(aes(x=DD_ln,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("DD cover (ln)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#square
DD_coverage %>% 
  ggplot(aes(x=DD_square,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("DD cover (square root)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#cube
DD_coverage %>% 
  ggplot(aes(x=DD_cube,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("DD cover (cube root)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#4.3.2 Shapiro-Wilk tests of transformed data
DD_C <- DD_coverage[DD_coverage$Type == "C",]
#View(DD_C)

shapiro.test(DD_C$Dead_Dryas)
shapiro.test(DD_C$DD_log10)
shapiro.test(DD_C$DD_ln)
shapiro.test(DD_C$DD_square)
shapiro.test(DD_C$DD_cube)

DD_O <- DD_coverage[DD_coverage$Type == "OTC",]
#View(DD_O)

shapiro.test(DD_O$Dead_Dryas)
shapiro.test(DD_O$DD_log10)
shapiro.test(DD_O$DD_ln)
shapiro.test(DD_O$DD_square)
shapiro.test(DD_O$DD_cube)

#4.3.3 F-test to check equal variance
var.test(DD_cube~Type,data=DD_coverage)

#Cube-transformed data will be used for Dead Dryas

#4.4.1 Lichen
lich_coverage <- coverage[,c("Type","Lich")]
#View(lich_coverage)
#log10
lich_coverage$lich_log10 <- log10(lich_coverage$Lich+1)
#ln
lich_coverage$lich_ln <- log(lich_coverage$Lich+1)
#square
lich_coverage$lich_square <- sqrt(lich_coverage$Lich)
#cube
lich_coverage$lich_cube <- (sign(lich_coverage$Lich)*abs(lich_coverage$Lich)^(1/3))

#histograms
#untransformed
lich_coverage %>% 
  ggplot(aes(x=Lich,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Lich cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#log10
lich_coverage %>% 
  ggplot(aes(x=lich_log10,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Lich cover (log 10)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#ln
lich_coverage %>% 
  ggplot(aes(x=lich_ln,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Lich cover (ln)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#square
lich_coverage %>% 
  ggplot(aes(x=lich_square,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Lich cover (square root)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#cube
lich_coverage %>% 
  ggplot(aes(x=lich_cube,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Lich cover (cube root)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#4.4.2 Shapiro-Wilk tests
lich_C <- lich_coverage[lich_coverage$Type == "C",]
#View(lich_C)

shapiro.test(lich_C$Lich)
shapiro.test(lich_C$lich_log10)
shapiro.test(lich_C$lich_ln)
shapiro.test(lich_C$lich_square)
shapiro.test(lich_C$lich_cube)

lich_O <- lich_coverage[lich_coverage$Type == "OTC",]
#View(lich_O)

shapiro.test(lich_O$Lich)
shapiro.test(lich_O$lich_log10)
shapiro.test(lich_O$lich_ln)
shapiro.test(lich_O$lich_square)
shapiro.test(lich_O$lich_cube)

#4.4.3 F test for equal variance
var.test(lich_cube~Type,data=lich_coverage)

#Cube-transformed data will be used for lichen coverage

#4.5.1 Litter
litter_coverage <- coverage[,c("Type","Litter")]
#View(litter_coverage)
#log10
litter_coverage$litter_log10 <- log10(litter_coverage$Litter)
#ln
litter_coverage$litter_ln <- log(litter_coverage$Litter)
#square
litter_coverage$litter_square <- sqrt(litter_coverage$Litter)
#cube
litter_coverage$litter_cube <- (sign(litter_coverage$Litter)*abs(litter_coverage$Litter)^(1/3))

#histograms
#untransformed
litter_coverage %>% 
  ggplot(aes(x=Litter,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Litter cover (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#log10
litter_coverage %>% 
  ggplot(aes(x=litter_log10,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Litter cover (log 10)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#ln
litter_coverage %>% 
  ggplot(aes(x=litter_ln,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Litter (ln)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#square
litter_coverage %>% 
  ggplot(aes(x=litter_square,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Litter cover (square root)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#cube
litter_coverage %>% 
  ggplot(aes(x=litter_cube,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Litter cover (cube root)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#4.5.2 Shapiro-Wilk tests
litter_C <- litter_coverage[litter_coverage$Type == "C",]
#View(litter_C)

shapiro.test(litter_C$Litter)
shapiro.test(litter_C$litter_log10)
shapiro.test(litter_C$litter_ln)
shapiro.test(litter_C$litter_square)
shapiro.test(litter_C$litter_cube)

litter_O <- litter_coverage[litter_coverage$Type == "OTC",]
#View(litter_O)

shapiro.test(litter_O$Litter)
shapiro.test(litter_O$litter_log10)
shapiro.test(litter_O$litter_ln)
shapiro.test(litter_O$litter_square)
shapiro.test(litter_O$litter_cube)

#Cube-transformed data will be used for litter

#4.5.3 F test for equal variance
var.test(litter_cube~Type,data=litter_coverage)

#4.6.1 Transforming forb
forb_coverage <- coverage[,c("Type","Forb")]
#View(forb_coverage)
#log10
forb_coverage$forb_log10 <- log10(forb_coverage$Forb)
#ln
forb_coverage$forb_ln <- log(forb_coverage$Forb)
#square
forb_coverage$forb_square <- sqrt(forb_coverage$Forb)
#cube
forb_coverage$forb_cube <- (sign(forb_coverage$Forb)*abs(forb_coverage$Forb)^(1/3))

#4.6.2 Shapiro-Wilk tests
forb_C <- forb_coverage[forb_coverage$Type == "C",]
#View(forb_C)

shapiro.test(forb_C$Forb)
shapiro.test(forb_C$forb_log10)
shapiro.test(forb_C$forb_ln)
shapiro.test(forb_C$forb_square)
shapiro.test(forb_C$forb_cube)

forb_O <- forb_coverage[forb_coverage$Type == "OTC",]
#View(forb_O)

shapiro.test(forb_O$Forb)
shapiro.test(forb_O$forb_log10)
shapiro.test(forb_O$forb_ln)
shapiro.test(forb_O$forb_square)
shapiro.test(forb_O$forb_cube)

#Cube will be used

#4.6.3 Equal variance test of transformed data
var.test(forb_cube~Type,data=forb_coverage)

#4.7.1 Transforming bare soil
bs_coverage <- coverage[,c("Type","Bare_soil")]
#View(bs_coverage)
#log10
bs_coverage$log10 <- log10(bs_coverage$Bare_soil+1)
#ln
bs_coverage$ln <- log(bs_coverage$Bare_soil+1)
#square
bs_coverage$square <- sqrt(bs_coverage$Bare_soil)
#cube
bs_coverage$cube <- (sign(bs_coverage$Bare_soil)*abs(bs_coverage$Bare_soil)^(1/3))

#4.7.2 Shapiro-Wilk tests
bs_C <- bs_coverage[bs_coverage$Type == "C",]
#View(forb_C)

shapiro.test(bs_C$Bare_soil)
shapiro.test(bs_C$log10)
shapiro.test(bs_C$ln)
shapiro.test(bs_C$square)
shapiro.test(bs_C$cube)

bs_O <- bs_coverage[bs_coverage$Type == "OTC",]
#View(forb_C)

shapiro.test(bs_O$Bare_soil)
shapiro.test(bs_O$log10)
shapiro.test(bs_O$ln)
shapiro.test(bs_O$square)
shapiro.test(bs_O$cube)

#transformation does not help, must use wilcox

```

The cube root transformation resulted in the following datasets
obtaining normal distribution and equal variance: Graminoids, forb,
lichen, bryophyte, litter and dead Dryas

Bare soil does not pass the test, and I must use Wilcox on it.

I can therefore perform a student's T-Test with the cube
root-transformed data, as well as the non-transformed data for the
variables that met the assumptions in the first place

```{r echo=T, results='hide',message=FALSE,warning=FALSE}
#5.5 Students T-test with transformed data
t.test(forb_cube ~ Type, data=forb_coverage, var.equal=TRUE, na.rm=TRUE)
t.test(gram_cube ~ Type, data=gram_coverage, var.equal=TRUE, na.rm=TRUE)
t.test(bryo_cube ~ Type, data=bryo_coverage, var.equal=TRUE, na.rm=TRUE)
t.test(DD_cube~Type,data=DD_coverage,var.equal=TRUE,na.rm=TRUE)
t.test(lich_cube~Type,data=lich_coverage,var.equal=TRUE,na.rm=TRUE)
t.test(litter_cube~Type,data=litter_coverage,var.equal=TRUE,na.rm=TRUE)

#5.6 Tests with the variables that did not need to be transformed
t.test(Shrub ~ Type, data=coverage, var.equal=TRUE, na.rm=TRUE)
t.test(Dryas ~ Type, data=coverage, var.equal=TRUE, na.rm=TRUE)

#5.7
wilcox.test(Bare_soil~Type,data=coverage)

```

### Results from T-tests

The results from the tests are Dwarf shrubs t = -0.076427, df = 40,
p-value = 0.9395

Dryas t = 0.54325, df = 40, p-value = 0.59

Forb t = -0.31737, df = 40, p-value = 0.7526

Gram t = -0.91931, df = 40, p-value = 0.3634

Bryophyte t = 5.166, df = 40, p-value = 6.965e-06

Lichen t = 3.0217, df = 40, p-value = 0.004369

Dead Dryas t = -2.491, df = 40, p-value = 0.01699

Litter t = -3.2437, df = 40, p-value = 0.002385

Wilcox test for bare soil W=274,5, p=0.1708

It is also possible to run the non-parametric Wilcox test on all of
them. It gives the same results for which mean values are significantly
different

```{r echo=T, results='hide',message=FALSE,warning=FALSE}
#5.8 Mann-Whitney-Wilcoxon test, unparametric test for non-normally distributed data
wilcox.test(Shrub~Type,data=coverage)
wilcox.test(Dryas~Type,data=coverage)
wilcox.test(Other_shrub~Type,data=coverage)
wilcox.test(Forb~Type,data=coverage)
wilcox.test(Gram~Type,data=coverage)
wilcox.test(Lich~Type,data=coverage)
wilcox.test(Bryo~Type,data=coverage)
wilcox.test(Litter~Type,data=coverage)
wilcox.test(Dead_Dryas~Type,data=coverage)
wilcox.test(Bare_soil~Type,data=coverage)
wilcox.test(Biol_crust~Type,data=coverage)
```

The last two steps are producing 1) a table with mean values, standard
errors, and all the necessary information from the tests 2) boxplots
with p values

### Adding p values to boxplots

adding p values from the tests

```{r echo=T, results='hide', fig.show='hide',warning=FALSE}
#Dwarf shrub
text_height_dwarf_shrub <- max(coverage$Shrub) + 5
dwarf_shrub_box_with_p <- dwarf_shrub_box+
  geom_signif(y_position = text_height_dwarf_shrub, xmin = 1, xmax = 2, annotation = "p = 0.9395",textsize=5)+
  scale_y_continuous(limits=c(0,110))
dwarf_shrub_box_with_p

#Dryas
dryas_box_with_p <- dryas_box+
  geom_signif(y_position = 102, xmin = 1, xmax = 2, annotation = "p = 0.59",textsize=5)+
  scale_y_continuous(limits=c(0,110))
dryas_box_with_p

#Forb
text_height_forb <- max(coverage$Forb) + 1.5
forb_box_with_p <- forb_box+
  geom_signif(y_position = text_height_forb, xmin = 1, xmax = 2, annotation = "p = 0.7526",textsize=5)+
  scale_y_continuous(limits=c(0,30))
forb_box_with_p

#Gram
text_height_gram <- max(coverage$Gram) + 1.5
gram_box_with_p <- gram_box+
  geom_signif(y_position = text_height_gram, xmin = 1, xmax = 2, annotation = "p = 0.3634",textsize=5)+
  scale_y_continuous(limits=c(0,30))
gram_box_with_p

#Lichens
text_height_lich <- max(coverage$Lich) + 3
lich_box_with_p <- lich_box+
  geom_signif(y_position = text_height_lich, xmin = 1, xmax = 2, annotation = "p = 0.004369",textsize=5,fontface="bold")+
  scale_y_continuous(limits=c(0,60))
lich_box_with_p

#Bryophytes
text_height_bryo <- max(coverage$Bryo) + 3
bryo_box_with_p <- bryo_box+
  geom_signif(y_position = text_height_bryo, xmin = 1, xmax = 2, annotation = "p = 6.965e-06",textsize=5)+
  scale_y_continuous(limits=c(0,50))
bryo_box_with_p

#Dead Dryas
dead_dryas_box_with_p <- dead_dryas_box+
  geom_signif(y_position = 47, xmin = 1, xmax = 2, annotation = "p = 0.01699",textsize=5)+
  scale_y_continuous(limits=c(0,55))
dead_dryas_box_with_p

#Litter
text_height_litter <- max(coverage$Litter) + 3
litter_box_with_p <- litter_box+
  geom_signif(y_position = text_height_litter, xmin = 1, xmax = 2, annotation = "p = 0.002385",textsize=5)+
  scale_y_continuous(limits=c(0,60))
litter_box_with_p

```

### Final boxplots for cover with p values

Plot for Dryas and dead Dryas Cover

```{r,warning=FALSE}
wrap_plots(dryas_box_with_p,dead_dryas_box_with_p,ncol=2,nrow=1)
```

For dwarf shrub and litter

```{r}
wrap_plots(dwarf_shrub_box_with_p,litter_box_with_p,ncol=2,nrow=1)
```

For forbs and graminoids

```{r}
wrap_plots(gram_box_with_p,forb_box_with_p,ncol=2,nrow=1)
```

For lichen and bryophyte

```{r}
wrap_plots(bryo_box_with_p,lich_box_with_p,ncol=2,nrow=1)
```

Dryas and dead Dryas

```{r echo=FALSE}
wrap_plots(dryas_box_with_p,dead_dryas_box_with_p,ncol=2,nrow=1)
```

Shrub and litter

```{r echo=FALSE}
wrap_plots(dwarf_shrub_box_with_p,litter_box_with_p,ncol=2,nrow=1)
```

Forb and graminoids

```{r echo=FALSE}
wrap_plots(forb_box_with_p,gram_box_with_p,ncol=2,nrow=1)
```

Bryophytes and lichens

```{r echo=FALSE}
wrap_plots(lich_box_with_p,gram_box_with_p,ncol=2,nrow=1)
```

### Cover: More tidy plots

Alternatively, it looks more tidy with everyone with boxplots in the
same plot

```{r}
long_data<-melt(env_data,id.vars=c("Plot","Type"))
cover_data <- filter(long_data,variable%in%c("Shrub","Forb","Gram","Lich","Bryo","Litter","Bare_soil"))
#View(cover_data)

cover_box<-cover_data %>% 
  ggplot(aes(x = variable, y = value, fill = Type)) +
  geom_boxplot()+
  geom_signif(
    y_position = c(45, 60),
    xmin = c(3.8, 4.8),  
    xmax = c(4.2, 5.2),
    annotation = c("***", "***"),
    textsize=10)+
  geom_signif(
    y_position = c(60),
    xmin = c(5.8),
    xmax = c(6.2),
    annotation = c("***"),
    textsize=10)+
  ylab("Cover (%)")+
  xlab("")+
  theme_bw(base_size=30)+
  scale_fill_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"),name="Treatment",labels=c("Control","Warmed"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.text.x=element_text(angle=45,vjust=1,hjust=1),
        legend.position="right")+
  scale_x_discrete(labels=c("Shrub" = "Dwarf shrub", 
                            "Forb" = "Forb",
                            "Gram" = "Graminoid",
                            "Bryo" = "Bryophyte",
                            "Lich" = "Lichen",
                            "Bare_soil"="Bare soil"))
#* means p>0.05, 
#** means p>0.01
#*** means p>0.001 
cover_box
```

The same for Dryas and dead Dryas

### Dryas cover: More tidy plots

```{r}
Dryas_data <- filter(long_data,variable%in%c("Dryas","Dead_Dryas"))

dryas_box <- Dryas_data %>% 
  ggplot(aes(x = variable, y = value, fill = Type)) +
  geom_boxplot() +
  geom_signif(
    y_position = c(60),
    xmin = c(1.8),  
    xmax = c(2.2),
    annotation = c("*"),
    textsize = 10) +
  ylab("Cover (%)") +
  xlab("") +
  theme_bw(base_size = 30) +
  scale_fill_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"), name = "Treatment", labels = c("Control", "Warmed")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_x_discrete(labels = c("Dryas" = expression(italic("Dryas")), 
                              "Dead_Dryas" = expression("Dead" ~ italic("Dryas"))))+
  theme(legend.position = "none")
dryas_box
```

# Vegetation height

The process is repeated as above

```{r echo=F, results='hide',warning=FALSE,fig.show='hide'}
#1.1 Veg height
vegh_box <- coverage %>% 
  ggplot(aes(x=Type,y=Veg_height,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Vegetation height (cm)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.2 Dryas
dryash_box <- coverage %>% 
  ggplot(aes(x=Type,y=Dryas_height,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Dryas height (cm)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

wrap_plots(vegh_box,dryash_box,
           ncol=2,nrow=1)

```

Looks like both are different

```{r,warning=FALSE}
#2.1 Histogram
veg_hist<-coverage %>% 
  ggplot(aes(x=Veg_height,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("Vegetation height (cm)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.2 Dryas
dryash_hist<-coverage %>% 
  ggplot(aes(x=Dryas_height,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Dryas height(cm)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

wrap_plots(veg_hist,dryash_hist,
           ncol=2,nrow=1)

```

C looks normally distributed, OTC looks left-skewed

### Transforming general vegetation height

```{r echo=F, results='hide', fig.show='hide',warning=FALSE}
coverage_c<-coverage[coverage$Type=="C",]
#View(coverage_c)
shapiro.test(coverage_c$Veg_height)
shapiro.test(coverage_c$Dryas_height)

coverage_o<-coverage[coverage$Type=="OTC",]
#View(coverage_o)
shapiro.test(coverage_o$Veg_height)
shapiro.test(coverage_o$Dryas_height)

var.test(Veg_height~Type,data=coverage)
var.test(Dryas_height~Type,data=coverage)
```

The tests are not happy with the distribution in the OTC plots.
Vegetation height does not have equal variance

```{r echo=F, results='hide', fig.show='hide'}
#Transforming vegetation height
veg_height <- coverage[,c("Type","Veg_height")]
#View(veg_height)
#log10
veg_height$veg_log10 <- log10(veg_height$Veg_height)
#ln
veg_height$veg_ln <- log(veg_height$Veg_height)
#square
veg_height$veg_square <- sqrt(veg_height$Veg_height)
#cube
veg_height$veg_cube <- (sign(veg_height$Veg_height)*abs(veg_height$Veg_height)^(1/3))

veg_height %>% 
  ggplot(aes(x=Veg_height,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Veg height (cm)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

veg_height %>% 
  ggplot(aes(x=veg_log10,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Veg height (log10)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

veg_height %>% 
  ggplot(aes(x=veg_ln,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Veg height (ln)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

veg_height %>% 
  ggplot(aes(x=veg_cube,fill=Type))+
  geom_histogram(bins=10,color="black")+
  xlab("Veg height (cube)")+
  ylab("Number of plots")+
  facet_wrap(~Type)

#shapiro-wilk-tests

veg_C <- veg_height[veg_height$Type == "C",]
#View(veg_C)

shapiro.test(veg_C$Veg_height)
shapiro.test(veg_C$veg_log10)
shapiro.test(veg_C$veg_ln)
shapiro.test(veg_C$veg_square)
shapiro.test(veg_C$veg_cube)

veg_O <- veg_height[veg_height$Type == "OTC",]
#View(veg_O)

shapiro.test(veg_O$Veg_height)
shapiro.test(veg_O$veg_log10)
shapiro.test(veg_O$veg_ln)
shapiro.test(veg_O$veg_square)
shapiro.test(veg_O$veg_cube)

#cube passes

#Variance test
var.test(veg_cube~Type,data=veg_height)

#T-test
t.test(veg_cube ~ Type, data=veg_height, var.equal=TRUE, na.rm=TRUE)
```

Using a cube root transformation, the vegetation height can be compared
with a students T-test. The results are:

t = -3.374, df = 40, p-value = 0.001656

## Dryas height: Non-parametric

Transforming the Dryas height data did not help, and a non-parametric
wilcox.test must be used

```{r echo=F, results='hide', fig.show='hide',warning=FALSE}
wilcox.test(Dryas_height~Type,data=coverage,conf.int=TRUE)
```

The p-value for Dryas height is 0.00909

## Final plots for height

Final plots for height

```{r echo=F, results='hide',warning=FALSE}
text_height_veg.height <- max(coverage$Veg_height) + 0.5
vegh_box_with_p <- vegh_box+
  geom_signif(y_position = text_height_veg.height, xmin = 1, xmax = 2, annotation = "p = 0.001656",textsize=5)+
  scale_y_continuous(limits=c(0,9))
#vegh_box_with_p

dryash_box_with_p <- dryash_box+
  geom_signif(y_position = 6, xmin = 1, xmax = 2, annotation = "p = 0.00909",textsize=5)+
  scale_y_continuous(limits=c(0,7))
#dryash_box_with_p

wrap_plots(vegh_box_with_p,dryash_box_with_p,ncol=2,nrow=1)

```

### Height: More tidy plots

All within a boxplot

```{r}
Dryas_data <- filter(long_data,variable%in%c("Dryas","Dead_Dryas"))

dryas_box <- Dryas_data %>% 
  ggplot(aes(x = variable, y = value, fill = Type)) +
  geom_boxplot() +
  geom_signif(
    y_position = c(60),
    xmin = c(1.8),  
    xmax = c(2.2),
    annotation = c("*"),
    textsize = 10) +
  ylab("Cover (%)") +
  xlab("") +
  theme_bw(base_size = 30) +
  scale_fill_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"), name = "Treatment", labels = c("Control", "Warmed")) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_x_discrete(labels = c("Dryas" = expression(italic("Dryas")), 
                              "Dead_Dryas" = expression("Dead" ~ italic("Dryas"))))+
  theme(legend.position = "none")

height_data <- filter(long_data,variable%in%c("Veg_height","Dryas_height"))

height_box<-height_data %>% 
  ggplot(aes(x = variable, y = value, fill = Type)) +
  geom_boxplot() +
  ylim(0,9)+
  geom_signif(
    y_position = c(8,6),
    xmin = c(0.8,1.8),  
    xmax = c(1.2,2.2),
    annotation = c("***","****"),
    textsize=10)+
  ylab("Height (cm)")+
  xlab("")+
  theme_bw(base_size=30)+
  scale_fill_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"),name="Treatment",labels=c("Control","Warmed"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("Veg_height" = "Vegetation", 
                            "Dryas_height" = expression(italic("Dryas"))))
height_box

wrap_plots(dryas_box,height_box,ncol=2,nrow=1)

nested <- (dryas_box|height_box)+ plot_annotation(tag_levels = 'A')
nested #view multi-panel figure
```

# Dryas flowers

#Flowers

### Boxplot and histogram

```{r echo=F, results='hide', fig.show='hide'}
#All flowering species

Flower_box<- coverage %>% 
  ggplot(aes(x=Type,y=Flowering_species,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Flowering species")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))
print(Flower_box)

Flowering_species_hist<-coverage %>% 
  ggplot(aes(x=Flowering_species,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("Flowering species")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
print(Flowering_species_hist)

#Dryas flowers
coverage$Dryas_flowering_percentage <- coverage$Dryas_flowering*100

Dryasflower_box<- coverage %>% 
  ggplot(aes(x=Type,y=Dryas_flowering_percentage,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Flowering Dryas (% of subplots)")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))
print(Dryasflower_box)

Dryasflower_hist<-coverage %>% 
  ggplot(aes(x=Dryas_flowering,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("Dryas flower subplots")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
print(Dryasflower_hist)
```

Looks like flowering species are similar, but subplots with flowering
Dryas are different. Distribution of flowering species could be normally
distributed. Distribution of Dryas flowers is left-skewed. Testing the
assumptions:

```{r,results='hide',,warning=FALSE}
#checking normality by shapiro-wil test
coverage_c<-coverage[coverage$Type=="C",]
#View(coverage_c)
shapiro.test(coverage_c$Dryas_flowering)
shapiro.test(coverage_c$Flowering_species)

coverage_o<-coverage[coverage$Type=="OTC",]
#View(coverage_o)
shapiro.test(coverage_o$Dryas_flowering)
shapiro.test(coverage_o$Flowering_species)

var.test(Dryas_flowering~Type,data=coverage)
var.test(Flowering_species~Type,data=coverage)

```

Must transform Dryas flowers. Flowering species can be tested as is.

```{r,results='hide',,warning=FALSE}
#does cube root transformation help?
#Transforming data
Dryas_flowering <- coverage[,c("Type","Dryas_flowering")]
#View(Dryas_flowering)
#log10
Dryas_flowering$log10 <- log10(Dryas_flowering$Dryas_flowering+1)
#ln
Dryas_flowering$ln <- log(Dryas_flowering$Dryas_flowering+1)
#square
Dryas_flowering$square <- sqrt(Dryas_flowering$Dryas_flowering)
#cube
Dryas_flowering$cube <- (sign(Dryas_flowering$Dryas_flowering)*abs(Dryas_flowering$Dryas_flowering)^(1/3))

#shapiro-wilk-tests of transformed data

flowers_C <- Dryas_flowering[Dryas_flowering$Type == "C",]
#View(flowers_C)

shapiro.test(flowers_C$Dryas_flowering)
shapiro.test(flowers_C$log10)
shapiro.test(flowers_C$ln)
shapiro.test(flowers_C$square)
shapiro.test(flowers_C$cube)

flowers_O <- Dryas_flowering[Dryas_flowering$Type == "OTC",]
#View(flowers_O)

shapiro.test(flowers_O$Dryas_flowering)
shapiro.test(flowers_O$log10)
shapiro.test(flowers_O$ln)
shapiro.test(flowers_O$square)
shapiro.test(flowers_O$cube)

#Cube root does help. What about variance?

var.test(cube~Type,data=Dryas_flowering)
```

Cube root transformation helps, now testing this with Welch's t.test

```{r echo=F, results='hide', fig.show='hide'}
t.test(Flowering_species~Type,data=coverage,var.equal=TRUE)
t.test(cube~Type,data=Dryas_flowering,var.equal=FALSE)
```

Flowering species: t = -0.46625, df = 40, p-value = 0.6436

Dryas flowers: t = -1.8048, df = 32.446, p-value = 0.08039

Final boxplot with p-value

```{r echo=F, results='hide'}
text_height_Flowering_species <- max(coverage$Flowering_species)+0.5
Flower_box_with_p <- Flower_box+
  geom_signif(y_position=text_height_Flowering_species,xmin=1,xmax=2,annotation="0.6436",textsize=5)+
  scale_y_continuous(limits=c(1,11))
  

Dryasflower_box_with_p <- Dryasflower_box+
  geom_signif(y_position = 93, xmin = 1, xmax = 2, annotation = "p = 0.08039",textsize=5)+
  scale_y_continuous(limits=c(0,100))

wrap_plots(Flower_box_with_p,Dryasflower_box_with_p,ncol=2,nrow=1)

```

### Flowers: More tidy plots

```{r}
coverage$Dryas_flowering_percentage <- coverage$Dryas_flowering*100

Flower_box2<- coverage %>% 
  ggplot(aes(x=Type,y=Flowering_species,fill=Type))+
  geom_boxplot()+
  ylab("Flowering species")+
  xlab("")+
  theme_bw(base_size=30)+
  theme(legend.position="none")+
  scale_fill_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"),name="Treatment",labels=c("Control","Warmed"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

Dryasflower_box2<- coverage %>% 
  ggplot(aes(x=Type,y=Dryas_flowering_percentage,fill=Type))+
  geom_boxplot()+
  xlab("")+
  theme_bw(base_size=30)+
  theme(legend.position="right")+
  scale_fill_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"),name="Treatment",labels=c("Control","Warmed"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))+
  ylab(expression(paste(italic("Dryas"), " flowers (% of subplots)")))
  #theme(legend.position="none")

wrap_plots(Flower_box2,Dryasflower_box2,nrow=1,ncol=2)
nested <- (Flower_box2|Dryasflower_box2)+ plot_annotation(tag_levels = 'A')
nested #view multi-panel figure

```

# Biodiversity

## Species richness

Boxplots

```{r echo=F, results='hide', fig.show='hide'}
richness <- read_excel("env_data.xlsx")

#1. Boxplots
#1.1 All species
species_box <- richness %>% 
  ggplot(aes(x=Type,y=Species,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("All species")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.2 Vascular plants
vasc_box <- richness %>% 
  ggplot(aes(x=Type,y=Vascular_plants,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Vascular plant species")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.3 Lichens
lichens_box <- richness %>% 
  ggplot(aes(x=Type,y=Lichens,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Lichen species")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.4 Bryophytes
bryophytes_box <- richness %>% 
  ggplot(aes(x=Type,y=Bryophytes,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Bryophyte species")+
  xlab("")+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))

#1.5 combining
wrap_plots(species_box,vasc_box,
           ncol=2,nrow=1)
wrap_plots(lichens_box,bryophytes_box,
           ncol=2,nrow=1)

#2 Histogram
#2.1 All species
species_hist<-richness %>% 
  ggplot(aes(x=Species,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("All species")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.2 Vascular plants
vasc_hist<-richness %>% 
  ggplot(aes(x=Vascular_plants,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("Vascular plant species")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.3 Lichens
lichens_hist<-richness %>% 
  ggplot(aes(x=Lichens,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("Lichen species")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.4 Bryophytes
bryophytes_hist<-richness %>% 
  ggplot(aes(x=Bryophytes,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("Bryophyte species")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#2.5 Combining
wrap_plots(species_hist,vasc_hist,
           lichens_hist,bryophytes_hist,
           ncol=2,nrow=2)

```

Checking assumptions

```{r echo=F, results='hide', fig.show='hide'}
#3 Checking normality distribution and equal variances
#3.1 Shapiro-wilk tests
richness_c<-richness[richness$Type=="C",]
#View(richness_c)

shapiro.test(richness_c$Species)
shapiro.test(richness_c$Vascular_plants)
shapiro.test(richness_c$Lichens)
shapiro.test(richness_c$Bryophytes)

richness_o<-richness[richness$Type=="OTC",]
#View(richness_o)
shapiro.test(richness_o$Species)
shapiro.test(richness_o$Vascular_plants)
shapiro.test(richness_o$Lichens)
shapiro.test(richness_o$Bryophytes)

#F test, Equal variances test
var.test(Species~Type,data=richness)
var.test(Vascular_plants~Type,data=richness)
var.test(Lichens~Type,data=richness)
var.test(Bryophytes~Type,data=richness)

```

Vascular plants and bryophytes do not need to be transformed, but
bryophytes do not have equal variance and I must use Welch's instead of
Student's. And as transforming the others does not help, I must use a
non-parametric test wilcox test for them

```{r echo=F, results='hide', fig.show='hide',warning=FALSE}
wilcox.test(Species~Type,data=richness)
t.test(Vascular_plants~Type,data=richness,var.equal=TRUE)
wilcox.test(Lichens~Type,data=richness)
t.test(Bryophytes~Type,data=richness,var.equal=FALSE)
```

Results All species: p = 0.001687


Vascular plants: t = 0.49283, df = 40, p-value = 0.6248


Lichens: p = 0.0101

Bryophytes: t = 2.8749, df = 29.867, p-value = 0.007381

Final boxplots with p values

```{r echo=F, results='hide'}
text_height_all_species <- max(richness$Species) + 2
species_box_with_p <- species_box+
  geom_signif(y_position = text_height_all_species, xmin = 1, xmax = 2, annotation = "p = 0.001687",textsize=5)+
  scale_y_continuous(limits=c(15,55))

text_height_vasc <- max(richness$Vascular_plants) + 2
vasc_box_with_p <- vasc_box+
  geom_signif(y_position = text_height_vasc, xmin = 1, xmax = 2, annotation = "p = 0.6248",textsize=5)+
  scale_y_continuous(limits=c(10,30))

text_height_lich_species <- max(richness$Lichens) + 2
lichens_box_with_p <- lichens_box+
  geom_signif(y_position = text_height_lich_species, xmin = 1, xmax = 2, annotation = "p = 0.0101",textsize=5)+
  scale_y_continuous(limits=c(0,30))

text_height_bryophyte_species <- max(richness$Bryophytes) + 2
bryophytes_box_with_p <- bryophytes_box+
  geom_signif(y_position = text_height_bryophyte_species, xmin = 1, xmax = 2, annotation = "p = 0.007381",textsize=5)+
  scale_y_continuous(limits=c(0,30))

wrap_plots(species_box_with_p,vasc_box_with_p,ncol=2,nrow=1)
wrap_plots(lichens_box_with_p,bryophytes_box_with_p,ncol=2,nrow=1)

```

### Richness: More tidy plots

everyone in the same boxplot with stars

```{r}
long_data<-melt(env_data,id.vars=c("Plot","Type"))
species_data <- filter(long_data,variable%in%c("Species","Vascular_plants","Lichens","Bryophytes"))
#View(species_data)

species_box<- species_data %>% 
  ggplot(aes(x = variable, y = value, fill = Type)) +
  geom_boxplot() +
  ylim(0,55)+
  geom_signif(
    y_position = c(48, 20),
    xmin = c(0.8, 2.8),  
    xmax = c(1.2, 3.2),
    annotation = c("***", "**"),
    textsize=10)+
  geom_signif(
    y_position = c(17),
    xmin = c(3.8),
    xmax = c(4.2),
    annotation = c("***"),
    textsize=10)+
  ylab("Species richness")+
  xlab("")+
  theme_bw(base_size=30)+
  scale_fill_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"),name="Treatment",labels=c("Control","Warmed"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("Species" = "Total","Vascular_plants"="Vascular plants"))
species_box

```

## Diversity and evenness

The Shannon diversity index H and Pielou's evennes J is calculated with
the vegan package with

```{r, echo=FALSE, eval=FALSE}
species <- read_excel("subplotfreq.xlsx")
shannon<-read_excel("subplotfreq.xlsx")
#View(shannon)

shannonind<-diversity(shannon,index="shannon")
#View(shannonind)
shannon2<-as.data.frame(shannonind)
#View(shannon2)

J<-diversity(shannon,index="shannon")/log(specnumber(shannon))
J2<-as.data.frame(J)
#View(J2)

#env_data <- cbind(shannon2,J2,env_data)
```

I did these analyses from before and have added them to the main
env_data dataset

```{r echo=F, results='hide', fig.show='hide'}
env_data<-read_excel("env_data.xlsx")

shannon_box <- env_data %>% 
  ggplot(aes(x=Type,y=Shannon,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Diversity Shannon's H")+
  xlab("Plot type")+
  scale_x_discrete(labels=c('Control', 'Warmed')) +
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

J_box<-env_data %>% 
  ggplot(aes(x=Type,y=J,fill=Type))+
  geom_violin(width=1)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Evenness Pioelou's J")+
  xlab("Plot type")+
  scale_x_discrete(labels=c('Control','Warmed'))+
  theme_bw(base_size=20)+
  theme(legend.position="none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

wrap_plots(shannon_box,J_box,ncol=2,nrow=1)
```

They both look different, let's test some assumptions

```{r,results='hide',warning=FALSE}
shannon_hist<-env_data %>% 
  ggplot(aes(x=Shannon,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("Diversity Shannon's H")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

J_hist<-env_data %>% 
  ggplot(aes(x=J,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("Evenness Pioelou's J")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
wrap_plots(shannon_hist,J_hist,ncol=2,nrow=1)

#Shannon looks right skewed, J looks normally distributed

#normality distribution and equal variance tests
data_c<-env_data[env_data$Type=="C",]
#View(shannon_c)
shapiro.test(data_c$Shannon)
shapiro.test(data_c$J)

data_o<-env_data[env_data$Type=="OTC",]
shapiro.test(data_o$Shannon)
shapiro.test(data_o$J)

var.test(Shannon~Type,data=env_data)
var.test(J~Type,data=env_data)

#Shannon passes the test for normal distribution, Pielou does not. None have equal variance

```

Welch's T-test must be used for Shannon. Must try to transform data for Pielou:

```{r}
#Transforming data

J <- env_data[,c("Type","J")]

View(J)
#log10
J$log10 <- log10(J$J+1)
#ln
J$ln <- log(J$J+1)
#square
J$square <- sqrt(J$J)
#squared
J$squared<-(J$J)^2
#cube
J$cube <- (sign(J$J)*abs(J$J)^(1/3))

#shapiro-wilk-tests of transformed data

J_C <- J[J$Type == "C",]
View(J_C)

shapiro.test(J_C$J)
shapiro.test(J_C$log10)
shapiro.test(J_C$ln)
shapiro.test(J_C$square)
shapiro.test(J_C$squared)
shapiro.test(J_C$cube)

J_O <- J[J$Type == "OTC",]
View(J_O)

shapiro.test(J_O$J)
shapiro.test(J_O$log10)
shapiro.test(J_O$ln)
shapiro.test(J_O$square)
shapiro.test(J_O$squared)
shapiro.test(J_O$cube)

#transformation does not help. must use wilcox
```

Tests:

```{r, echo=F, results='hide'}
t.test(Shannon ~ Type, data=env_data, var.equal=FALSE)
wilcox.test(J~Type,data=env_data)
```

Shannon diversity index H: t = 4.051, df = 28.202, p-value = 0.0003626
Pielou's evennes index J W = 378, p-value = 2.882e-05

Boxplots with p values

```{r}
text_height_shannon <- max(env_data$Shannon) + 0.1
shannon_box_with_p <- shannon_box+
  geom_signif(y_position = text_height_shannon, xmin = 1, xmax = 2, annotation = "p = 0.0003626",textsize=5)+
  scale_y_continuous(limits=c(2,4))

text_height_J <- max(env_data$J) + 0.01
J_box_with_p <- J_box+
  geom_signif(y_position = text_height_J, xmin = 1, xmax = 2, annotation = "p = 2.882e-05",textsize=5)+
  scale_y_continuous(limits=c(0.7,1))

wrap_plots(shannon_box_with_p,J_box_with_p,ncol=2,nrow=1)

```

### Diversity and evenness: More tidy plots

boxplots with stars

```{r}
shannon_box2<-env_data %>% 
  ggplot(aes(x=Type,y=Shannon,fill=Type))+
  geom_boxplot()+
  ylim(2.3,3.8)+
  ylab("Diversity")+
  xlab("")+
  theme_bw(base_size=30)+
  theme(legend.position="none")+
  scale_fill_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"),name="Treatment",labels=c("Control","Warmed"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))+
          geom_signif(
            y_position = c(3.5),
            xmin = c(1),  
            xmax = c(2),
            annotation = c("****"),
            textsize=10)

J_box2<-env_data %>% 
  ggplot(aes(x=Type,y=J,fill=Type))+
  geom_boxplot()+
  ylim(0.8,1)+
  ylab("Evenness")+
  xlab("")+
  theme_bw(base_size=30)+
  #theme(legend.position="none")+
  scale_fill_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"),name="Treatment",labels=c("Control","Warmed"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_discrete(labels=c("C"="Control","OTC"="Warmed"))+
  geom_signif(
    y_position = c(0.97),
    xmin = c(1),  
    xmax = c(2),
    annotation = c("****"),
    textsize=10)

wrap_plots(shannon_box2,J_box2,ncol=2,nrow=1)

nested <- (shannon_box2|J_box2)+ plot_annotation(tag_levels = 'A')
nested #view multi-panel figure
```

#Soil moisture

```{r, warning=FALSE}
moisture <- read_excel("moisture3.xlsx")

moisture_box <- moisture %>% 
  ggplot(aes(x=Type,y=Moisture,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("Soil moisture (%)")+
  xlab("Plot type")+
  scale_x_discrete(labels=c('Control', 'Warmed'))+
  theme_bw(base_size=20)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
moisture_box

moisture_hist<-moisture %>% 
  ggplot(aes(x=Moisture,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("Soil moisture (%)")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
moisture_hist

```

Looks different, looks left-skewed

What about the differences between dates of measurements?

```{r,warning=FALSE}
ggplot(data=moisture,aes(x=Date,y=Moisture,color=Type))+
  geom_smooth(aes(color=Type),method="loess")+
  scale_color_manual(values=c("#00AFBB","#FC4E07"),name="Treatment",labels=c("Control","Warming"))+
  ylab("Soil moisture (%)")+
  xlab("Time")+
  theme_classic()

```

Seems to be variation between the dates. Best use a mixed linear model

```{r}
#with plot and date as crossed effect
mixed.moisture.lmer1<-lmer(Moisture~Type+(1|Date)+(1|Plot),data=moisture)
summary(mixed.moisture.lmer1)

#with date and type as fixed, plot as random
mixed.moisture.lmer2<-lmer(Moisture~Type+Date+(1|Plot),data=moisture)
summary(mixed.moisture.lmer2)
```

Which model to choose? I try using the AIC

```{r}
library(MASS)
AIC(mixed.moisture.lmer1,mixed.moisture.lmer2)
```

the first model has lowest AIC, but there is a sign effect between date

plot with p for model 1

```{r warning=FALSE}
moisture_box_with_p<-moisture_box+
  geom_signif(y_position = 95,xmin=1,xmax=2,annotation="p = 0.0613",textsize=5)+
  scale_y_continuous(limits=c(0,100))
moisture_box_with_p
  
```

plot for model 2

```{r}
moisture_box2<- moisture %>% 
  ggplot(aes(x=factor(Date),y=Moisture,fill=Type))+
  geom_boxplot(alpha=0.8)+
  theme_classic(base_size=30)+
  ylab("Soil moisture (%)")+
  xlab("Date")+
  scale_fill_manual(values=c("Control"="#00AFBB","Warmed"="#FC4E07"))+
  theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))
moisture_box2
```

But I am not interested in the variation between dates, just the plot types. therefore model 1 is chosen with this plot

```{r}
moisture_box <- moisture %>% 
  ggplot(aes(x=Type,y=Moisture,fill=Type))+
  geom_boxplot()+
  ylab("Soil moisture (%)")+
  xlab("")+
  scale_x_discrete(labels=c('Control', 'Warmed'))+
  theme_bw(base_size=30)+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
moisture_box
```


# NDVI

```{r warning=FALSE}
#0.2 inspecting data
ndvi <- read_excel("ndvi2.xlsx")
#View(ndvi)

#1 boxplot
ndvi_box <- ndvi %>% 
  ggplot(aes(x=Type,y=NDVI,fill=Type))+
  geom_violin(width=1.0)+
  geom_boxplot(width=0.4,color="black",alpha=0.2)+
  ylab("NDVI")+
  xlab("Plot type")+
  theme_bw(base_size=20)+
  scale_x_discrete(labels=c('Control', 'Warmed'))+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ndvi_box

#histogram
#2 Histograms for checking if data is normally distributed
ndvi_hist<-ndvi %>% 
  ggplot(aes(x=NDVI,fill=Type),show.legend = FALSE)+
  geom_histogram(bins=10,color="black")+
  xlab("NDVI")+
  ylab("Number of plots")+
  facet_wrap(~Type)+
  theme_bw()+
  theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ndvi_hist
```

They do not look very different. The data looks somewhat normally
distributed, but its the distribution of the residuals I am interested
in here. What about between dates?

```{r warning=FALSE}
ggplot(data=ndvi,aes(x=Date,y=NDVI,color=Type))+
  geom_smooth(aes(color=Type),method="loess")+
  scale_color_manual(values=c("#00AFBB","#FC4E07"),name="Treatment",labels=c("Control","Warming"))+
  ylab("NDVI")+
  xlab("Time")+
  theme_classic()
```

Making mixed linear models

```{r}
mixed.ndvi.lmer<-lmer(NDVI~Type+(1|Date)+(1|Plot),data=ndvi)
summary(mixed.ndvi.lmer)

mixed.ndvi.lmer2<-lmer(NDVI~Type+Date+(1|Plot),data=ndvi)
summary(mixed.ndvi.lmer2)

AIC(mixed.ndvi.lmer,mixed.ndvi.lmer2)

```

There is sign. variation between dates, but not plot types

Adding the p value to the boxplot of model 1

```{r warning=FALSE}
ndvi_box_with_p<-ndvi_box+
  geom_signif(y_position = 0.9,xmin=1,xmax=2,annotation="p = 0.683",textsize=5)+
  scale_y_continuous(limits=c(0,1))
ndvi_box_with_p
```

plot for model 2

```{r}
ndvi_box2<-ndvi %>% 
  ggplot(aes(x=factor(Date),y=NDVI,fill=Type))+
  geom_boxplot(alpha=0.8)+
  theme_classic(base_size=30)+
  ylab("NDVI")+
  xlab("Date")+
  scale_fill_manual(values=c("Control"="#00AFBB","Warmed"="#FC4E07"))+
  theme(axis.text.x=element_text(angle=45,vjust=1,hjust=1))
ndvi_box2
```

But I am not interested in the variation between dates, just the plot types. therefore model 1 is chosen with this plot

```{r}
ndvi_box <- ndvi %>% 
  ggplot(aes(x=Type,y=NDVI,fill=Type))+
  geom_boxplot()+
  ylab("NDVI")+
  xlab("")+
  theme_bw(base_size=30)+
  scale_x_discrete(labels=c('Control', 'Warmed'))+
  #theme(legend.position = "none")+
  scale_fill_manual(values=c("#00AFBB","#FC4E07"),name="Treatment",labels=c("Control","Warmed"))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ndvi_box
```

### Moisture and NDVI: More tidy plots

```{r}
wrap_plots(moisture_box,ndvi_box,ncol=2,nrow=1)
nested <- (moisture_box|ndvi_box)+ plot_annotation(tag_levels = 'A')
nested #view multi-panel figure
```



# Tables for summarizing values

As I have not find a quick way to give me the mean values and standard
errors, Iam using dplyr and the summarise function to calculate these
values for the variables I have tested

```{r echo=F, fig.show='hide'}
summaries <- env_data %>% 
  group_by(Type) %>% 
  dplyr::summarise("Dryas cover (%)"=mean(Dryas,na.rm=T),
            "Dryas cover SE (%)"=sd(Dryas,na.rm=T)/sqrt(length(Dryas[!is.na(Dryas)])),
            "Dead Dryas cover (%)"=mean(Dead_Dryas,na.rm=T),
            "Dead Dryas cover SE (%)"=sd(Dead_Dryas,na.rm=T)/sqrt(length(Dead_Dryas[!is.na(Dryas)])),
            "Dwarf shrub cover (%)"=mean(Shrub),
            "Dwarf shrub cover SE (%)"=sd(Shrub)/sqrt(length(Shrub)),
            "Litter cover (%)"=mean(Litter),
            "Litter cover cover SE (%)"=sd(Litter)/sqrt(length(Litter)),
            "Forb cover (%)"=mean(Forb),
            "Forb cover SE (%)"=sd(Forb)/sqrt(length(Forb)),
            "Graminoid cover (%)"=mean(Gram),
            "Graminoid SE (%)"=sd(Gram)/sqrt(length(Gram)),
           "Lichen cover (%)"=mean(Lich),
            "Lichen cover SE (%)"=sd(Lich)/sqrt(length(Lich)),
            "Bryophyte cover (%)"=mean(Bryo),
            "Bryophyte cover SE (%)"=sd(Bryo)/sqrt(length(Bryo)),
            "Vegetation height (cm)"=mean(Veg_height),
            "Vegetation height SE (cm)"=sd(Veg_height)/sqrt(length(Veg_height)),
            "Dryas height (cm)"=mean(Dryas_height,na.rm=T),
            "Dryas heigth SE (cm)"=sd(Dryas_height,na.rm=T)/sqrt(length(Dryas_height[!is.na(Dryas_height)])),
            "Flowering species"=mean(Flowering_species),
            "Flowering species SE"=sd(Flowering_species)/sqrt(length(Flowering_species)),
            "Subplots with Dryas flowers (%)"=mean(Dryas_flowering*100,na.rm=T),
            "Subplots with Dryas flowers SE"=sd(Dryas_flowering*100,na.rm=T)/length(Dryas_flowering[!is.na(Dryas_flowering)]),
            "All species"=mean(Species),
            "All species SE"=sd(Species)/sqrt(length(Species)),
            "Vascular plant species"=mean(Vascular_plants),
            "Vascular plants species SE"=sd(Vascular_plants)/sqrt(length(Vascular_plants)),
            "Lichen species"=mean(Lichens),
            "Lichen species SE"=sd(Lichens)/sqrt(length(Lichens)),
            "Bryophyte species"=mean(Bryophytes),
            "Bryophyte species SE"=sd(Bryophytes)/sqrt(length(Bryophytes)),
            "Shannon diversity index H"=mean(Shannon),
            "Shannon_ diversity H SE"=sd(Shannon)/sqrt(length(Shannon)),
            "Pielou's evennes index J"=mean(J),
            "Pielou's evennes index J SE"=sd(J)/sqrt(length(J)),
           "Soil moisture mean"=mean(Moisture),
           "Soil moisture SE"=sd(Moisture)/sqrt(length(Moisture)),
           "NDVI mean"=mean(NDVI),
           "NDVI SE"=sd(NDVI)/sqrt(length(NDVI)),
            "Bare soil mean"=mean(Bare_soil),
            "Bare soil SE"=sd(Bare_soil)/sqrt(length(Bare_soil)))
summaries2<-round_df(summaries,digits=3)
summary_table<-t(summaries2)
#summary_table<-print(summaries,digits=3)
#kable(summary_table,format="html",digits=c(3,3))
kable(summary_table,format="html",digits=3)

```

###The final plots

```{r}
cover_box
dryas_box
height_box
wrap_plots(Flower_box2,Dryasflower_box2,nrow=1,ncol=2)
species_box
wrap_plots(shannon_box2,J_box2,ncol=2,nrow=1)

```

#Wilcox test for each species

```{r}
groups<-read_excel("groups.xlsx")
species2 <- read_excel("percentage2.xlsx")
env_data<-read_excel("env_data.xlsx")

merged_data3<-cbind(env_data,species2)  
names(merged_data3) <- gsub(" ", "_", names(merged_data3))

subset_data <- merged_data3[, c(1,2, 27:111)]

new_data <- subset_data

new_data[, 3:87] <- subset_data[, 3:87] * 100

df <- new_data

results_list <- list()

# Iterate over columns 3 to 87
for (i in 3:87) {
  # Calculate counts for values above 0 in each group
  counts_group_1 <- sum(df[df$Type == "C", i] > 0)  # Assuming "C" is group 1
  counts_group_2 <- sum(df[df$Type == "OTC", i] > 0)  # Assuming "OTC" is group 2
  
  # Calculate mean and standard error for each group
  group_means <- tapply(df[, i], df[, 2], mean)
  group_std_err <- tapply(df[, i], df[, 2], function(x) sd(x) / sqrt(length(x)))
  
  # Perform Wilcoxon rank sum test for each column based on groups in column 2
  result <- wilcox.test(df[, i] ~ df[, 2])
  
  # Determine significance level
  significance <- ifelse(result$p.value < 0.001, "***",
                         ifelse(result$p.value < 0.01, "**",
                                ifelse(result$p.value < 0.05, "*", "")))
  
  # Store the results in the list along with column names, counts for each group, means, standard errors, and significance
  results_list[[i]] <- c(
    Column_Name = colnames(df)[i],
    Count_Group_1 = counts_group_1,
    Count_Group_2 = counts_group_2,
    Mean_Group_1 = round(group_means[1], 3),
    Mean_Group_2 = round(group_means[2], 3),
    Std_Err_Group_1 = round(group_std_err[1], 3),
    Std_Err_Group_2 = round(group_std_err[2], 3),
    W = result$statistic,
    p_value = round(result$p.value, 3),
    Significance = significance
  )
}

results_df <- do.call(rbind, results_list)
results_df <- results_df[, c(1:4, 6, 5, 7:10)]
colnames(results_df)[1]<-"Species"
colnames(results_df)[2]<-"Counts (Control)"
colnames(results_df)[3]<-"Counts (Warmed)"
colnames(results_df)[4]<-"Mean (Control)"
colnames(results_df)[5]<-"SE (Control)"
colnames(results_df)[6]<-"Mean (Warmed)"
colnames(results_df)[7]<-"SE (Warmed)"
colnames(results_df)[8]<-"W"
colnames(results_df)[9]<-"p"
colnames(results_df)[10]<-" "

results_df<-cbind(results_df, groups[,2, drop=FALSE])
results_df <- results_df[, c(11,1:10)]
View(results_df)
```

#Ordination analyses

#Loading packages, data and setting up ordination

```{r warning=FALSE}
library(readxl)
library(vegan)
library(ggplot2)
library(ggrepel)
library(patchwork)
#library(tidyverse)
#library(tidyr)

# Import the species data # 
species <- read_excel("subplotfreq.xlsx")
nmds <- metaMDS(species, distance = "bray", k = 2)
sppscores(nmds) <- species

# importing the environmental variables
env_data <- read_excel("env_data.xlsx")

#Group names
groups<- read_excel("groups.xlsx")

```

The most important axes as x and y

```{r warning=FALSE}
plot(nmds,"species",choices = c(1,2))
orditorp(nmds, "species") 
 
axes1and2 <- scores(nmds)
 dim(axes1and2)
 NMDS1<-axes1and2$sites[,1]
 NMDS2<-axes1and2$sites[,2]
 co=c("blue","red")
 
 plot(NMDS1,NMDS2,xlim = c(-1.5, 1.2), ylim = c(-1,1.2))
 points(NMDS1[env_data$Type == "C"],NMDS2[env_data$Type == "C"], col="blue" ,pch=15,cex=1.5)
 points(NMDS1[env_data$Type == "OTC"],NMDS2[env_data$Type == "OTC"], col="red" ,pch=17,cex=1.5)
 legend(0.6,-0.6, legend=c("Control", "Warming"),
        col=c("blue", "red"), pch=(c(15,17)), cex=1.2)
 
 ordiellipse(nmds,groups=env_data$Type,kind="e",conf=0.95,label=FALSE,cex=1,col=co)
 env_var<- envfit(nmds ~ Lich +Litter+Bryo+Shrub+Moisture+Veg_height+Dryas_height, choices = c(1,2), data = env_data,na.rm=TRUE)
 plot(env_var,col="black")
 
```

###Making the same plot but with ggplot

```{r}
scores.nmds <- as.data.frame(scores(nmds, "sites") ) 
 scores.nmds$Type <- env_data$Type 
species_scores <- as.data.frame(scores(nmds, "species"))  
 merged_data<-cbind(groups,species_scores)  
 
 env <- envfit(nmds ~ Lich +Litter+Bryo+Shrub+
                 Moisture+Veg_height+
                 NDVI+Forb+Gram+Bare_soil,
               choices = c(1,2), data = env_data,na.rm=TRUE)
 
 en_coord_cont = as.data.frame(scores(env, "vectors")) * ordiArrowMul(env)
 row.names(en_coord_cont) <- c("Lichens","Litter","Bryophytes","Dwarf shrubs",
                              "Moisture","Vegetation height","NDVI",
                              "Forbs", "Graminoids","Bare soil")

 axes1and2 <- scores(nmds)
 dim(axes1and2)
 NMDS1<-axes1and2$sites[,1]
 NMDS2<-axes1and2$sites[,2]
 
nmds.plot <- ggplot(data = scores.nmds, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(aes(color = Type, shape= Type), size = 4, alpha = 0.5) +
  scale_colour_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"),name="Treatment",labels=c("Control","Warmed")) +
  geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
               data = en_coord_cont, size = 1, alpha = 0.85, colour = "grey30",
               arrow = arrow(length = unit(0.5, "cm"))) +
  #geom_vline(xintercept = c(0), linetype = 2, linewidth = 1) +
  #geom_hline(yintercept = c(0), linetype = 2, linewidth = 1) + 
  geom_label_repel(data = en_coord_cont, aes(x = NMDS1, y = NMDS2), colour = "grey30", 
                   fontface = "bold", label = row.names(en_coord_cont),
                   alpha = 0.7,
                   label.padding = 0.15, label.size=NA,
                   size=7,
                   box.padding = 0.35,
                   point.padding = 0.5) + 
  stat_ellipse(aes(color = Type), geom = "polygon", level = 0.95, alpha = 0.1) +
  scale_shape_manual(values = c("C" = 16, "OTC" = 17),guide="none") +
  theme_classic(base_size=20)+
  theme(axis.title.y=element_text(size=30),
        axis.title.x=element_text(size=30),
        axis.text.x = element_text(size = 26),
        axis.text.y = element_text(size = 26),
        legend.text = element_text(size=20),
        legend.title = element_text(size=22),
        legend.position = "right")+
  #theme(plot.margin = grid::unit(c(0, 0, 0, 0), "mm"),
  #     axis.line = element_line(colour = "black", size = 1),
  #    axis.ticks.length = unit(.25, "cm"),
  #   panel.background = element_blank(),
  #  axis.title.y = element_text(size = 14, face = "bold"),
  # axis.text.y = element_text(size = 16),
  #axis.text.x = element_text(size = 16),
  #axis.title.x = element_text(size = 14, face = "bold"),
  #plot.title = element_text(size = 14),
  #legend.text = element_text(size = 12),
  #legend.title = element_text(size = 12, face = "bold"),
#legend.position = "none") +
scale_linetype_manual(values = c(1, 1, 1, 5))+
  xlim(c(-1.4, 1.1)) +  
  ylim(c(-1, 1.2))+
  xlab("")
print(nmds.plot)

#theme(legend.position = c(0.9,0.2))
#nudge_x = 0.2, nudge_y = 0.1, 
#legend.position = c (1.4,0.5)
 
```

ordination plot for species, with different shapes, colours and fonts

```{r}
species.plot<-ggplot(merged_data, aes(x = NMDS1, y = NMDS2, color = Group)) +
  geom_point(fill="darkblue") +  # Adjust size as needed
  geom_text_repel(aes(label = Species,fontface=Fontface), size = 5, box.padding = 0.1, point.padding = 0.1) +  
  xlim(c(-1.4, 1.1)) +  
  ylim(c(-1, 1.2)) +  
  labs(x = "NMDS1", y = "NMDS2", color = "Group") +  
  theme_classic(base_size=20)+
  scale_color_manual(name="Functional group",breaks=c("Bryophyte","Lichen","Forb","Graminoid","Shrub"),
                     values=c("Bryophyte"="darkblue","Lichen"="darkblue",
                              "Forb"="goldenrod","Graminoid"="goldenrod","Shrub"="goldenrod"),
                     labels=c("Bryophytes(italic)","Lichens",
                              "Forbs","Graminoids(italic)","Dwarf shrubs(bold)"))+
  scale_shape_manual(name="Functional group",breaks=c("Bryophyte","Lichen","Forb","Graminoid","Shrub"),
                     values=c("Bryophyte"=18,"Lichen"=18,"Forb"=18,"Graminoid"=18,"Shrub"=18),
                     labels=c("Bryophytes (italic)","Lichens",
                              "Forbs","Graminoids (italic)","Dwarf shrubs (bold)"))+
  theme(axis.title.y=element_text(size=30),
        axis.title.x=element_text(size=30),
        axis.text.y = element_text(size = 26),
        axis.text.x = element_text(size = 26),
        legend.title = element_text(size=22),
        legend.text = element_text(size=20),
        legend.position="right")


species.plot
```

```{r}
#wrap_plots(nmds.plot,species.plot,nrow=2,ncol=1)
(nmds.plot/species.plot)+ plot_annotation(tag_levels = 'A')
```

Checking non-linear relationships in the data

```{r}
#Moisture
moisture.ordi <- ordisurf(nmds ~ Moisture, data = env_data)
summary(moisture.ordi)

#NDVI
NDVI.ordi <- ordisurf(nmds ~ NDVI, data = env_data)
summary(NDVI.ordi)

#Bryo
bryo.ordi <- ordisurf(nmds ~ Bryo, data = env_data)
summary(bryo.ordi)

#Litter
litter.ordi <- ordisurf(nmds ~ Litter, data = env_data)
summary(litter.ordi)

#Lichen
lichen.ordi <- ordisurf(nmds ~ Lich, data = env_data)
summary(lichen.ordi)

#Veg. height
veg.height.ordi <- ordisurf(nmds ~ Veg_height, data = env_data)
summary(veg.height.ordi)

#Dryas height
dryas.height.ordi <- ordisurf(nmds ~ Dryas_height, data = env_data)
summary(dryas.height.ordi)

#Shrub
shrub.ordi <- ordisurf(nmds ~ Shrub, data = env_data)
summary(shrub.ordi)

```

Checking the number of dimensions

```{r}
k_vec <- 1:10
stress <- numeric(length(k_vec))
species_dij <- metaMDSdist(species, trace = FALSE)
set.seed(25)
for(i in seq_along(k_vec)) {
  sol <- metaMDSiter(species_dij, k = i,
                     trace = FALSE)
  stress[i] <- sol$stress
}


plot(k_vec, stress, type = "b", ylab = "Stress",
     xlab = "Dimensions")
```

Four dimensions looks to have reduced the stress enough, but it was also under 0.20 with two dimensions, so I choose two.

Creating DCA

```{r}
dta<-species #raw species data 

dca.av <-decorana(species)

```

Comparing DCA and NMDS

```{r}
plot(dca.av)
orditorp(dca.av, "species") 
ordiellipse(dca.av,groups=env_data$Type,kind="e",conf=0.95,label=FALSE,cex=1,col=co)

plot(nmds)
orditorp(nmds, "species")
ordiellipse(nmds,groups=env_data$Type,kind="e",conf=0.95,label=FALSE,cex=1,col=co)

```

```{r}
#Extracting DCA-plot scores:#Note that origin=FALSE implies that origo of the ordination diagram is moved 
#from the centroid to the lower end of each axis
dca1<-scores(dca.av,display="sites",origin=FALSE)[,1]
dca2<-scores(dca.av,display="sites",origin=FALSE)[,2]
dca3<-scores(dca.av,display="sites",origin=FALSE)[,3]
dca4<-scores(dca.av,display="sites",origin=FALSE)[,4]

#Calculation of correlation coefficients:####

kendall_cor <- function(x,y) {
  require(dplyr)
  
  z <- cor.test(as.numeric(x)
                ,as.numeric(y), method = "kendall")
  q <- cbind(z$estimate, z$p.value)
  q <- as.data.frame(q)
  names(q)[1] <- ( "tau")
  names(q)[2] <- ( "p.verdi")
  q$n <- nrow(file)
  return(q)
}

```

```{r}
# two axes
nmds.scores <- scores(nmds,tidy = T,display = "sites")

dca1nmds1 <- kendall_cor(dca1,nmds.scores[,1])
dca1nmds2 <- kendall_cor(dca1,nmds.scores[,2])


dc1.cor <- rbind(dca1nmds1,dca1nmds2)
print(dc1.cor)

dca2nmds1 <- kendall_cor(dca2,nmds.scores[,1])
dca2nmds2 <- kendall_cor(dca2,nmds.scores[,2])

dc2.cor <- rbind(dca2nmds1,dca2nmds2)
print(dc2.cor)

```

NMDS1 correlates with DCA1: tau=-0.8, p=9.5*10^-19
NMDS2 correlates with DCA2: tau=0.4, p=4.83*10^-5

The correlation should be above 0,4

CCA comparing species compositions with plot types

```{r}
#Type plot
cca.analysis <- cca(species ~ Type, data = env_data)

anova(cca.analysis)
plot(cca.analysis)
```

CCA-analysis gives a significant p value below 0.05

What about RDA?

```{r}
rda.analysis1 <- rda(species ~ Type, data = env_data)

anova(rda.analysis1)
plot(rda.analysis1)
```

the p-value is below 0.05

Correlation between the axes and environmental variables

```{r}
#Correlation between axes and env variables
#moisture
cor.test(NMDS1, env_data$Moisture, method="k")
cor.test(NMDS2, env_data$Moisture, method="k")

#Bryo
cor.test(NMDS1, env_data$Bryo, method="k")
cor.test(NMDS2, env_data$Bryo, method="k")

#Lich
cor.test(NMDS1, env_data$Lich, method="k")
cor.test(NMDS2, env_data$Lich, method="k")

#Litter
cor.test(NMDS1, env_data$Litter, method="k")
cor.test(NMDS2, env_data$Litter, method="k")

#Dryas height
cor.test(NMDS1, env_data$Dryas_height, method="k")
cor.test(NMDS2, env_data$Dryas_height, method="k")

#Veg height
cor.test(NMDS1, env_data$Veg_height, method="k")
cor.test(NMDS2, env_data$Veg_height, method="k")

#NDVI
cor.test(NMDS1, env_data$NDVI, method="k")
cor.test(NMDS2, env_data$NDVI, method="k")

#Shrub
cor.test(NMDS1, env_data$Shrub, method="k")
cor.test(NMDS2, env_data$Shrub, method="k")

```

There are some correlations:

-Moisture with NMDS1

-NDVI with both\
- Shrub with both\
- Gram with NMDS2\
- Lich with NMDS1\
-Bryo with NMDS1\
-Litter with NMDS2\
-Veg. height with NMDS1\

###Repeating with percentage data

```{r}
# Import the species data # 
species2 <- read_excel("percentage2.xlsx")
nmds2 <- metaMDS(species2, distance = "bray", k = 2)
sppscores(nmds2) <- species2

scores.nmds2 <- as.data.frame(scores(nmds2, "sites") ) 
 scores.nmds2$Type <- env_data$Type 
species_scores2 <- as.data.frame(scores(nmds2, "species"))



 merged_data2<-cbind(groups,species_scores2)  
 
 env2 <- envfit(nmds2 ~ Lich +Litter+Bryo+Shrub+
                 Moisture+Veg_height+
                 NDVI+Forb+Gram+Bare_soil,
               choices = c(1,2), data = env_data,na.rm=TRUE)
 
 en_coord_cont2 = as.data.frame(scores(env2, "vectors")) * ordiArrowMul(env2)
 row.names(en_coord_cont2) <- c("Lichens","Litter","Bryophytes","Dwarf shrubs",
                              "Moisture","Vegetation height","NDVI",
                              "Forbs", "Graminoids","Bare soil")
 
 axes1and2 <- scores(nmds2)
 dim(axes1and2)
 NMDS1<-axes1and2$sites[,1]
 NMDS2<-axes1and2$sites[,2]
 

```

```{r}
nmds.plot2 <- ggplot(data = scores.nmds2, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(aes(color = Type, shape= Type), size = 4, alpha = 0.5) +
  scale_colour_manual(values = c("C" = "#00AFBB", "OTC" = "#FC4E07"),name="Treatment",labels=c("Control","Warmed")) +
  geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
               data = en_coord_cont2, size = 1, alpha = 0.85, colour = "grey30",
               arrow = arrow(length = unit(0.5, "cm"))) +
  #geom_vline(xintercept = c(0), linetype = 2, linewidth = 1) +
  #geom_hline(yintercept = c(0), linetype = 2, linewidth = 1) + 
  geom_label_repel(data = en_coord_cont2, aes(x = NMDS1, y = NMDS2), colour = "grey30", 
                   fontface = "bold", label = row.names(en_coord_cont2),
                   alpha = 0.7,
                   label.padding = 0.15, label.size=NA,
                   size=7,
                   box.padding = 0.35,
                   point.padding = 0.5) + 
  stat_ellipse(aes(color = Type), geom = "polygon", level = 0.95, alpha = 0.1) +
  scale_shape_manual(values = c("C" = 16, "OTC" = 17),guide="none") +
  theme_classic(base_size=20)+
  theme(axis.title.y=element_text(size=30),
        axis.title.x=element_text(size=30),
        axis.text.x = element_text(size = 26),
        axis.text.y = element_text(size = 26),
        legend.text = element_text(size=20),
        legend.title = element_text(size=22),
        legend.position = "right")+
scale_linetype_manual(values = c(1, 1, 1, 5))+
  xlim(c(-1.4, 1.1)) +  
  ylim(c(-1, 1.2))+
  xlab("")
print(nmds.plot)
```


```{r}
species.plot2<-ggplot(merged_data2, aes(x = NMDS1, y = NMDS2, color = Group)) +
  geom_point(fill="darkblue") +  # Adjust size as needed
  geom_text_repel(aes(label = Species,fontface=Fontface), size = 5, box.padding = 0.1, point.padding = 0.1) +  
  xlim(c(-1.4, 1.1)) +  
  ylim(c(-1, 1.2)) +  
  labs(x = "NMDS1", y = "NMDS2", color = "Group") +  
  theme_classic(base_size=20)+
  scale_color_manual(name="Functional group",breaks=c("Bryophyte","Lichen","Forb","Graminoid","Shrub"),
                     values=c("Bryophyte"="darkblue","Lichen"="darkblue",
                              "Forb"="goldenrod","Graminoid"="goldenrod","Shrub"="goldenrod"),
                     labels=c("Bryophytes(italic)","Lichens",
                              "Forbs","Graminoids(italic)","Dwarf shrubs(bold)"))+
  scale_shape_manual(name="Functional group",breaks=c("Bryophyte","Lichen","Forb","Graminoid","Shrub"),
                     values=c("Bryophyte"=18,"Lichen"=18,"Forb"=18,"Graminoid"=18,"Shrub"=18),
                     labels=c("Bryophytes (italic)","Lichens",
                              "Forbs","Graminoids (italic)","Dwarf shrubs (bold)"))+
  theme(axis.title.y=element_text(size=30),
        axis.title.x=element_text(size=30),
        axis.text.y = element_text(size = 26),
        axis.text.x = element_text(size = 26),
        legend.title = element_text(size=22),
        legend.text = element_text(size=20),
        legend.position="right")

```

```{r}
(nmds.plot2/species.plot2)+ plot_annotation(tag_levels = 'A')
```


repeating the CCA, RDA and correlation analyses

```{r}

#CCA

#Type plot
cca.analysis1 <- cca(species2 ~ Type, data = env_data)

anova(cca.analysis1)
plot(cca.analysis1)

#sign difference in distance between the plot types when it comes to species compositions

#RDA
rda.analysis1 <- rda(species2 ~ Type, data = env_data)

anova(rda.analysis1)
plot(rda.analysis1)

#RDA does not show sign difference


```

```{r}
#Correlation between axes and env variables
#moisture
cor.test(NMDS1, env_data$Moisture, method="k")
cor.test(NMDS2, env_data$Moisture, method="k")
#moisture is positively correlated with nmds1

#Bryo
cor.test(NMDS1, env_data$Bryo, method="k")
cor.test(NMDS2, env_data$Bryo, method="k")
#Bryophyte cover is sign positively correlated with nmds1

#Lich
cor.test(NMDS1, env_data$Lich, method="k")
cor.test(NMDS2, env_data$Lich, method="k")
#Lichen cover is sign positively correlated with nmds1

#Litter
cor.test(NMDS1, env_data$Litter, method="k")
cor.test(NMDS2, env_data$Litter, method="k")
#Litter cover is sign positively correlated with nmds2

#Veg height
cor.test(NMDS1, env_data$Veg_height, method="k")
cor.test(NMDS2, env_data$Veg_height, method="k")
#Vegetation height is not sign correlated with either

#NDVI
cor.test(NMDS1, env_data$NDVI, method="k")
cor.test(NMDS2, env_data$NDVI, method="k")
#NDVI is negatively correlated with nmds1

#Dwarf hrub
cor.test(NMDS1, env_data$Shrub, method="k")
cor.test(NMDS2, env_data$Shrub, method="k")
#Shrub is negatively correlated with nmds1, and negatively with nmds2

#Forb
cor.test(NMDS1, env_data$Forb, method="k")
cor.test(NMDS2, env_data$Forb, method="k")
#Forb is positively correlated with nmds2

#Gram
cor.test(NMDS1, env_data$Gram, method="k")
cor.test(NMDS2, env_data$Gram, method="k")
#Gram is positively correlated with nmds1 and nmds2

#Bare soil
cor.test(NMDS1, env_data$Bare_soil, method="k")
cor.test(NMDS2, env_data$Bare_soil, method="k")
#Bare soil is not sign correlated to eiterh
```

Creating DCA

```{r}
dta2<-species2 #raw species data 

dca.av2 <-decorana(species2)

```

Comparing DCA and NMDS

```{r}
plot(dca.av2)
orditorp(dca.av2, "species") 
ordiellipse(dca.av2,groups=env_data$Type,kind="e",conf=0.95,label=FALSE,cex=1,col=co)

plot(nmds2)
orditorp(nmds2, "species")
ordiellipse(nmds2,groups=env_data$Type,kind="e",conf=0.95,label=FALSE,cex=1,col=co)


```

```{r}
#Extracting DCA-plot scores:#Note that origin=FALSE implies that origo of the ordination diagram is moved 
#from the centroid to the lower end of each axis
dca1<-scores(dca.av2,display="sites",origin=FALSE)[,1]
dca2<-scores(dca.av2,display="sites",origin=FALSE)[,2]
dca3<-scores(dca.av2,display="sites",origin=FALSE)[,3]
dca4<-scores(dca.av2,display="sites",origin=FALSE)[,4]

#Calculation of correlation coefficients:####

kendall_cor <- function(x,y) {
  require(dplyr)
  
  z <- cor.test(as.numeric(x)
                ,as.numeric(y), method = "kendall")
  q <- cbind(z$estimate, z$p.value)
  q <- as.data.frame(q)
  names(q)[1] <- ( "tau")
  names(q)[2] <- ( "p.verdi")
  q$n <- nrow(file)
  return(q)
}

```

```{r}
# two axes
nmds.scores2 <- scores(nmds2,tidy = T,display = "sites")

dca1nmds1 <- kendall_cor(dca1,nmds.scores2[,1])
dca1nmds2 <- kendall_cor(dca1,nmds.scores2[,2])


dc1.cor <- rbind(dca1nmds1,dca1nmds2)
print(dc1.cor)

dca2nmds1 <- kendall_cor(dca2,nmds.scores2[,1])
dca2nmds2 <- kendall_cor(dca2,nmds.scores2[,2])

dc2.cor <- rbind(dca2nmds1,dca2nmds2)
print(dc2.cor)

```

