---
title: "FinalCode_BCBR_ANOVA_Biomass_2022_2_22"
author: "Jessica Rath"
date: "5/12/22"
output: html_document
---
Outline for ANOVA
## Step 1: Tidy your data and prepare it for analysis (load tidy, data, org by group using unite)
## Step 2: Assess what kind of variables we're dealing with. 
## Step 3: Visualize the data
## Step 4: Check assumptions
## Step 5: Conduct the statistical test
## Step 6: Assess the output. Does it make sense?
## Step 7: Post-Hoc Test
## Step 8: Plot and annotate to show significance

##work log 2022_5_14
TO DO Omit Degraded from all ANOVAS


##work log 3/9/22
DONE 
Export pdf of shrub size visualization, assumptions tables and ANOVA results for presentation

TO DO (at later date)
Create new df for by SPECIES for BiomassLogDataCS and BiomassLogDataREC
Assess Shrub Volume by top one to three species with the most biomass or abundance

##work log 3/2/2022
Done: 
corrected log data
Corrected the mean of the log for data
Created new dataframe for ANOVA BiomassLogData
Created new df for visuals BiomassLogSummary (includes se)

TO DO:
create ANOVA by Species 
ID Species to do top 2 species in study

##work log 2/28/22
Done:
calculated cylindrical volume with AH
relabeled biomass ggplots 
check that all is working
export biomass visuals for presentation

To Do:
lines 190-219 Ask AH to go over Kruskall Wallis test
Kruskall wallis: one-way ANOVA on ranks[1] is a non-parametric method for testing whether samples originate from the same distribution.

##Work Log 2/22/22
To Do: 
Calculate Biomass using cylinder formula
Check new totals working
write code for Kruskall Wallis test: one-way ANOVA on ranks[1] is a non-parametric method for testing whether samples originate from the same distribution


##Work Log 2/10/2022
JR Create MetaData, adjust annotations, replace data and file names with class requirements
JR Export files

##Work Log 01/29/2022
BH created Abundance as a function of SpeciesID with Treatment Fill (REC and CS separately)
We updated Biomass plot with log values of biomass and SE

To Do
JR attempts code for ANOVA assumptions (normal distribution and variance)


##Work log post 01/26/2022
Created summary table of biomass (AKA Shrub Size) with log10 values of biomass and standard error of logged biomass
Created ggplot of biomass ~ treeatment split up by species for each treatment

To-do:
Add standard error bars to ggplot^ 
widen columns to make error bars more visible
change biomass column to shrub size

3 Shrub size issues: Shrub size total, shrub size total for each species, and looking at one common species between intact and restored treatments to calculate effect size
Then move on to abundance and species diversity index (All the ANOVAs)

#Worklog 01/22/2022
What we did:
imported data into R as .csv
calculated biomass for each plant
filtered for only living CS and REC 
created new dataframe that includes treatment, category, SpeciesID, and Biomass
Created summary table for mean biomass
Did preliminary visualization of CS and REC Biomass by species and treatment
MALLAU biomass contains outliers that exceed 3e+05 cm^3, exlcuded from visualization 

Next Steps:
Going to split REC and CS plots for biomass into two separate plots to correct scaling
Calculate SEM for mean biomass
Compute summary statistics in summary table and create visualizations for species diversity and abundance
Create Summary tables for species diversity and abundance
Calculate Total Biomass/Shrub Size of all species ~ Treatment (X) (MALLAU yes/no)
Calculate Biomass/Shrub Size by Species  ~ Treatment (X)
Calculate Effect Size that is Calculate % difference from the global mean across all treatments of one species. ARTCAL Biomass by treatment /AVG ARTCAL in all treatments
Take MALLAU out
Find One Common Species Biomass/Shrub Size as a function of Treatment

### work log 1/19/2022
Done
created new rmd file for CS data
imported compiled CS data set
installed packages

Nest steps
Read CSV file into R Studio
Check Data format (num,char)
Create New Object with select columns
Summary statistics 

*****************************************************************************8
```{r Importing}
library(tidyverse)
library(fs)
library(readxl)
library(ggpubr)
library(rstatix)
```

#import data
```{r}
BCBR_PRIMARY_20220209 <- read_csv("BCBR_PRIMARY_20220209.csv")
View(BCBR_PRIMARY_20220209)
```

#Check the data
```{r}
view(BCBR_PRIMARY_20220209)
head(BCBR_PRIMARY_20220209)
```

#Removing the Notes or other non-necessary Columns by selecting columns of interest and renaming object as "CrownSprout_tidy"
```{r}
#view (BCBR_PRIMARY_20220209)

CrownSprout_tidy <- BCBR_PRIMARY_20220209 %>% select(1:11)
CrownSprout_tidy
```

#############################################################################
                     BIOMASS/SHRUB SIZE
############################################################################
#You are creating two data sets to use. 
#One for ANOVA and one for your PLOTS.
#Raw Data = BiomassLogData > for ANOVA
#Summarized Data - BiomassLogSummary > for PLOTS

#Create BiomassLogData (Total, REC, CS) for ANOVA 
```{r}
BiomassLogDataTotal <- CrownSprout_tidy %>% 
filter(Status == 'living', Category == c('CS','REC')) %>%
filter(Treatment != "Degraded") %>% 
drop_na() %>% mutate(Log_Biomass = log(biomass))

BiomassLogDataREC <- CrownSprout_tidy %>% 
filter(Status == 'living', Category == c('REC')) %>%
filter(Treatment != "Degraded") %>% 
drop_na() %>% mutate(Log_Biomass = log(biomass))

BiomassLogDataCS <- CrownSprout_tidy %>% 
filter(Status == 'living', Category == c('CS')) %>%
filter(Treatment != "Degraded") %>% 
drop_na() %>% mutate(Log_Biomass = log(biomass))


#view(BiomassLogDataTotal)
#view(BiomassLogDataREC)
#view(BiomassLogDataCS)

```

#Create BiomassLogSummary for PLOTS
```{r}
BiomassLogSummary <- CrownSprout_tidy %>% 
filter(Status == 'living', Category == c('CS','REC')) %>%
drop_na() %>% 
filter(Treatment != "Degraded") %>% 
select(Treatment, Category, SpeciesID, biomass) %>% 
group_by(Treatment, SpeciesID, Category) %>% 
#MeanLogBiomass_first calculation mean of biomass, untransformed mean
#se_biomass_second calculation is mean of the logbiomass
#se_logbiomass_third se of the logbiomass is calculated
summarise(Mean_biomass = mean(biomass),
          MeanLogBiomass = mean(log(biomass)), 
          se_biomass = sd(biomass)/sqrt(n()),
          se_logbiomass = sd(log(biomass))/sqrt(n()))
view(BiomassLogSummary)
```


#Create Mean Biomass Visualization for CS&REC across Species + Treatment w/ error bars
#facet_grid by category (REC and CS - two graphs one picture)
```{r}

ggplot(BiomassLogSummary) + (aes(x = Treatment, y = MeanLogBiomass, fill = SpeciesID)) +
  geom_col (position = position_dodge(width = 0.9, preserve = "single")) + 
  geom_errorbar(aes(ymin = (MeanLogBiomass - se_logbiomass),
                    ymax = (MeanLogBiomass + se_logbiomass) ),
                    position = position_dodge(width = 0.9, preserve = "single"),
                    width = 0.1) +
                    facet_grid(rows = vars(Category) ) +
                    theme(axis.text.x = element_text(angle = -90, size = 12, vjust = 0.5) ) + 
  labs(y = "Shrub Volume")

```
#Do not use for report
```{r}
jpeg("shrubsizeplot1.jpg", width = 700, height = 350) 
ggplot(BiomassLogSummary) + (aes(x = Treatment, y = MeanLogBiomass, fill = SpeciesID)) +
  geom_col (position = position_dodge(width = 0.9, preserve = "single")) + 
  geom_errorbar(aes(ymin = (MeanLogBiomass - se_logbiomass),
                    ymax = (MeanLogBiomass + se_logbiomass) ),
                    position = position_dodge(width = 0.9, preserve = "single"),
                    width = 0.1) +
                    facet_grid(rows = vars(Category) ) +
                    theme(axis.text.x = element_text(angle = -90, size = 12, vjust = 0.5) ) + 
  labs(y = "Shrub Volume")
dev.off() 
```

#normality
#assessing visually as first pass with histogram, plotting separately for each treatment. 
#They are evenly spaced from middle hump. Low number of observations n=12.
```{r}
ggplot(BiomassLogSummary, aes(x = MeanLogBiomass)) + geom_histogram() + facet_grid(Treatment~.)
```

#normality
#Amount of obsv is visualized below within bins.
```{r}
ggplot(BiomassLogSummary, aes(x = MeanLogBiomass)) + geom_histogram() + facet_grid(Treatment~SpeciesID)
```

##########################################################################
                        SHRUB VOLUME One Factor ANOVAs
##########################################################################
#Question 1
#Total CS Shrub Volume (Biomass) per Quad ~ Treatment (+Block)
#Total REC Shrub Volume (Biomass) per Quad ~ Treatment (+Block)

#First, group by Polygon, Quad, and Treatment and then calculate the sum of the biomass for each quad by Polygon. Use this biomass sum by quad for your ANOVAs. 
```{r}
BiomassLogDataRECQuad <- BiomassLogDataREC %>% group_by(Polygon, Quad, Treatment) %>% summarise(SumBiomass = sum(biomass))
BiomassLogDataCSQuad <- BiomassLogDataCS %>% group_by(Polygon, Quad, Treatment) %>% summarise(SumBiomass = sum(biomass))

view(BiomassLogDataRECQuad)
view(BiomassLogDataCSQuad )
```
#ANOVA, add + 1 so she does not get cranky about zero values in your data. 
#code: anova_results <- aov(values ~ ind, data=)
```{r}
RECVolANOVA <- aov(log(SumBiomass + 1) ~ Treatment + Polygon, data = BiomassLogDataRECQuad)
CSVolANOVA <- aov(log(SumBiomass + 1) ~ Treatment + Polygon, data = BiomassLogDataCSQuad)

summary(RECVolANOVA)
summary(CSVolANOVA)
```

#plot your ANOVA to check your assumptions. 
#The first Residual plot should be evenly distributed left to right and top to bottom. 
#Check the QQ plot too. Should follow along the line. A few outliers are OK bc ANOVA is robust analyses and can handle it. 
```{r}
plot(RECVolANOVA)
plot(CSVolANOVA)
```
#Write up: Results show no significant difference in total Crown Sprout biomass between Intact and Restored sites. There is also no significant difference in total Recruitment biomass between Intact and Restored sites. 


NEXT STEPS
Create new df for by SPECIES for BiomassLogDataCS and BiomassLogDataREC
Assess Shrub Volume by top one to three species with the most biomass or abundance
```{r}
#Per Species CS Shrub Volume (Biomass) per Quad ~ Treatment (+Block)
#Per Species REC Shrub Volume (Biomass) per Quad ~ Treatment (+Block)
```

Assess Shrub Volume by top one to three species with the most biomass or abundance
```{r}
#ARTCAL CS Shrub Volume (Biomass) per Quad ~ Treatment (+Block)
#ARTCAL REC Shrub Volume (Biomass) per Quad ~ Treatment (+Block)

```


*****************************************************
BONEYARD
*****************************************************
#attempt at showing total REC and CS Shrub size totals by Treatment. Need new DF with totals of MeanLogBiomass. 
```{r}
ggplot(BiomassLogSummary) +
    geom_col(aes(x = Treatment, y = MeanLogBiomass, fill = Treatment), position = position_dodge2(), width = 0.5) +
  scale_fill_manual(values = c("#5E4780","#4E6135")) +
  geom_errorbar(aes(x = Treatment, 
                    ymin = MeanLogBiomass - se_logbiomass,
                    ymax = MeanLogBiomass + se_logbiomass),
                    width = 0.1) +
                    labs(y = "Shrub Volume") +
                    theme(legend.position = "none") +
                    facet_grid(rows = vars(Category) ) +
                    
  #theme(axis.text.x = element_text(angle = -90, size = 12, vjust = 0.5) ) + 
  annotate("text", x = 1, y = 14, label = "A", size = 6) +
  annotate("text", x = 2, y = 14, label = "B", size = 6) 
``

#Kruskall Wallis test for Biomass
#Advised not to use Kruskall Wallis bc I plotted the ANOVAs and looked at Residuals and QQ plot there. Data was not normal so I logged the data then checked Residuals and QQ again--Homogeneity--enough.

#Advised not to use T-tests. T-tests in R is one of the most common tests in statistics. So, we use it to determine whether the means of two groups are equal to each other. The assumption for the test is that both groups are sampled from normal distributions with equal variances



