Carryover_SoyResults.Rmd

---
title: "Carryover_SoyResults"
author: "Kolby Grint"
date: "3/4/2021"
output: pdf_document
---


```{r, include = FALSE}
library(car)
library(ggplot2)
library(emmeans)
library(lme4)
library(lmerTest)
library(patchwork)
library(glmmTMB)
library(dplyr)
library(tidyverse)
library(multcomp)
```



```{r, include = FALSE, echo = FALSE}
getwd()

Soybean = read.csv(file="Soybean_MasterData.csv")

SoybeanComp = read.csv(file="SBYieldComponents.csv")
Soybean$Rep= factor(Soybean$Rep)
Soybean$Year= factor(Soybean$Year)

CC= read.csv(file= "CoverCropData.csv")

order= c("Tillage", "No-Till", "Cover Crop")
order2= c("CTRL", "CL25", "CL50", "ME25", "ME50")
```


```{r, include = FALSE, echo = FALSE}

SoybeanComp= SoybeanComp %>%
   mutate(Location = fct_recode(Location,
                                "Arlington" = "ARL",
                                "Lancaster" = "LAN",
                                "Lincoln" = "HAV"),
          Site_Yr = fct_recode(Site_Yr,
                          "Arlington 2019" = "ARL_19",
                          "Arlington 2020" = "ARL_20",
                          "Lincoln 2019" = "HAV_19",
                          "Lincoln 2020" = "HAV_20",
                          "Lancaster 2019" = "LAN_19",
                          "Lancaster 2020" = "LAN_20"),
          Soil = fct_recode(Soil,
                            "Cover Crop" = "NT+CC",
                            "No-Till" = "NT",
                            "Tillage" = "Till"),
          Total = round(Total, 3),
          Hundred = round(Hundred, 4))

SoybeanComp1= SoybeanComp %>%
  filter(!is.na(Total)) %>%
  filter(!is.na(Hundred)) %>%
  filter(!is.na(Seeds))


Soybean1 <- Soybean %>% 
  mutate(Rep = as_factor(Rep),
         Year = as_factor(Year),
         Soil= fct_recode(Soil, 
                          "Cover Crop" = "NT+CC",
                          "No-Till" = "NT",
                          "Tillage" = "Till"),
         Site_crop_yr = fct_recode(Site_crop_yr,
                          "Arlington 2019" = "ARL_SB_19",
                          "Arlington 2020" = "ARL_SB_20",
                          "Lincoln 2019" = "HAV_SB_19",
                          "Lincoln 2020" = "HAV_SB_20",
                          "Lancaster 2019" = "LAN_SB_19",
                          "Lancaster 2020" = "LAN_SB_20"),
         Location = fct_recode(Location,
                           "Arlington" = "Arlington",
                           "Lancaster" = "Lancaster",
                           "Lincoln" = "Havelock"),
         Crop.Canopy = round(Crop.Canopy/100, 2),
         yield = round(yield, 4)) %>%
  filter(!is.na(yield)) %>%
  filter(!is.na(Crop.Canopy)) %>%
  filter(!is.na(Stand.Count)) %>%
  janitor::clean_names()

SBCC <- CC %>%
  filter(Crop == "Soybean") %>%
  mutate(Rep = as_factor(Rep),
         Year = as_factor(Year),
         Site_Crop_Yr = fct_recode(Site_Crop_Yr,
                          "Arlington 2019" = "ARL_SB_19",
                          "Arlington 2020"  = "ARL_SB_20",
                          "Lincoln 2019" = "HAV_SB_19",
                          "Lincoln 2020" = "HAV_SB_20",
                          "Lancaster 2019" = "LAN_SB_19",
                          "Lancaster 2020" = "LAN_SB_20"),
         Canopy = round(Canopy/100, 3)) %>%
  filter(!is.na(Biomass_kg)) %>%
  janitor::clean_names()
```

#Analysis Procedure

Prior to analysis I took the approach of plotting the response variables with box-plots to visualize treatment differences with soil management between locations. The intention of this was to visualize the differences between locations as well determine if I feel comfortable pooling things within a location (not testing for a site-year or year effect). This has been the desired direction to simplify results for publication. When I didn't feel the data allowed for this I tested for differences between site-years in a condensed model. This approach is up for more discussion, as it eliminates the potential to view anomalies in the data which might be insightful and scientifically interesting. We are assuming editors won't like complex findings.


I would also like to point out that for every linear-mixed-effects model I am testing to see that the model meets the visual assumptions for normal distribution of residual and evenly distributed variance. I will demonstrate for the first model, and only include it in the output for future analyses when it is suspected that transformations need to be made.

More analyses were made than are going to be included in this pdf. I chose not to include them all so as to condense the output to the approach and considerations that I think are most valuable for the publication.

# Early seaston stand Counts Analysis

```{r, fig.width= 12, fig.height= 6}
Soybean1 %>%
  ggplot(aes(x = soil, y = stand_count, color = location)) +
  geom_boxplot() +
  facet_grid( ~ year)
```

```{r, fig.width= 12, fig.height= 6}
Soybean1 %>%
  ggplot(aes(x = herb, y = stand_count, color = location)) +
  geom_boxplot() +
  facet_grid(soil~ year)
```


Based on these visual representations it doesn't appear to me that there are  any consistent patterns as a function of soil management, herbicide treatment, location, or year. There does appear to be differences within a location across seasons. Therefore, I think it is best if we test for site-year differences as a fixed effect in a condensed models for our initial approach.

## Condensed Stand Count Model
```{r}
sb_stand= lmer(stand_count~ site_crop_yr*soil*herb + (1|site_crop_yr:rep), data= Soybean1)

Soybean1 %>% 
  ggplot(aes(x = stand_count)) +
  geom_density()
#Looks a little skewed left but good to see the distribution is one hump

qqnorm(resid(sb_stand))
#Residuals are concentrated in a "straight" line. 
#I am satisfied that this response variable is normally distributed

plot(sb_stand)
#Residuals are spread out evenly verticall and horizontally.
#I am satisfied that this response variable has an evenly distributed variance.
```

```{r}
anova(sb_stand)
#site-year fixed effect significant.
```
Herbicide carryover and soil management had no effect on early season stand counts in soybean.


### Separation of means for each site-year based on soil management
### Wasn't significant, done for fun (and to look at usefulness of means for yield later?)!
```{r}
sb_cc_stand_soillsmeans= lsmeans(sb_stand ,~ soil|site_crop_yr, contr="pairwise", adjust="none", type="response")
sb_cc_stand_SoilCLD<- CLD(sb_cc_stand_soillsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)
sb_cc_stand_SoilCLD
```
I wouldn't pay attention to the groupings since this interaction wasn't significant (although close).
The reduced stand for soil management with a cover crop could help explain final grain yield.


### Separation of means by site-year based on ANOVA
I probably won't include mean separation for later analyses often, I just wanted to show you for the first one.
```{r}
sb_cc_stand_lsmeans= lsmeans(sb_stand ,~ site_crop_yr, contr="pairwise", adjust="none", type="response")
#Chris, you can ignore the warnings. It is just a suggestion!

sb_cc_standCLD<- CLD(sb_cc_stand_lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

sb_cc_standCLD
```

```{r, include= FALSE, echo = FALSE}
sb_cc_standCLD <- as_tibble(sb_cc_standCLD) %>% 
  rename(stand_count = lsmean) %>%
  mutate(site_crop_yr = fct_recode(site_crop_yr,
                          "Arlington 2019" = "ARL_SB_19",
                          "Arlington 2020" = "ARL_SB_20",
                          "Lincoln 2019" = "HAV_SB_19",
                          "Lincoln 2020" = "HAV_SB_20",
                          "Lancaster 2019" = "LAN_SB_19",
                          "Lancaster 2020" = "LAN_SB_20")) %>%
  arrange(site_crop_yr)
```

### Stand Count figure -- Not useful for paper but interesting
```{r, echo= FALSE}
ggplot(sb_cc_standCLD, aes(x= site_crop_yr, y= stand_count, color= site_crop_yr)) +
  geom_point(size= 3) +
  geom_errorbar(aes(ymin= lower.CL, ymax= upper.CL), width= .3) +
  geom_text(aes(label = .group), nudge_y = 15000, size= 5) +
  geom_jitter(data = Soybean1 ,mapping = aes(y = stand_count), alpha = 0.2) +
  coord_flip() +
  theme_bw() +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "none") +
  labs(title = "Soybean Stand Counts")
```
Soybean stand counts were planted at 370,000 seeds ha^-1 in NE and ~346,000 seeds ha^-1 in WI

# Crop Canopy

```{r}
Soybean1 %>%
  ggplot(aes(x = soil, y = crop_canopy, color = location)) +
  geom_boxplot() +
  facet_grid(~ year)
```

Based on the variability in canopy coverage between cropping seasons within a location, 
I created a condensed model, similar to the stand count model, with site-year as 
a fixed effect.

I also decided to perform the analysis using a separate model for each site-year. 
The first model and figures resulting from this analysis will also be displayed.


## Condensed Canopy model

a glmmTMB (generalized linear) model was used instead of a standard linear mixed-effects model 
because % canopy cover is a non-continuous variable. There are no assumption to meet for this 
type of model.

```{r, }
sb_canopy = glmmTMB(crop_canopy~ soil*herb*site_crop_yr + (1|rep:site_crop_yr), data= Soybean1, beta_family(link="logit"))

Anova(sb_canopy)
#all 3 main fixed effects significant and the soil:site-year interaction
```
Means were separated for the Soil management:site-year interaction and herbicide 
treatment fixed effect separately. I delayed presenting the results until comparison to the separated analysis later on.

```{r, include= FALSE, echo= FALSE}
sb_canopy_herbemmeans = emmeans(sb_canopy, ~ herb, contr= "pairwise", adjust= "none", type= "response")
sb_canopy_herbCLD = CLD(sb_canopy_herbemmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

sb_canopy_soilemmeans = emmeans(sb_canopy, ~ soil|site_crop_yr, contr= "pairwise", adjust= "none", type= "response")
sb_canopy_soilCLD = CLD(sb_canopy_soilemmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

sb_canopy_herbCLD <- as_tibble(sb_canopy_herbCLD) %>%
  rename(crop_canopy = response) %>%
  mutate(herb= factor(herb, levels= order2)) %>%
  arrange(herb)

sb_canopy_soilCLD <- as_tibble(sb_canopy_soilCLD) %>%
  rename(crop_canopy = response) %>%
  mutate(soil= factor(soil, levels= order)) %>%
  arrange(soil)
```

```{r}
sb_canopy_soilCLD = sb_canopy_soilCLD %>%
  separate(site_crop_yr, c("location", "year"), " ") %>%
  mutate(year = as.factor(year))

```


```{r, echo= FALSE, include= FALSE}
sb_canopy <- ggplot(sb_canopy_soilCLD, aes(x= soil, y= crop_canopy*100, color= soil)) +
  geom_point(size= 3) +
  geom_errorbar(aes(ymin= lower.CL*100, ymax= upper.CL*100), width= .3) +
  geom_text(aes(label = .group), nudge_y = 12, size= 6) +
  geom_jitter(data = Soybean1 ,mapping = aes(y = crop_canopy*100), alpha = 0.2) +
  #coord_flip() +
  facet_grid(year~location) +
  theme_bw() +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "none",
        plot.title = element_text(size = 25),
        #axis.title.x = element_text(size = 25),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 25),
        axis.text.y = element_text(size = 15),
        axis.text.x = element_text(angle= 45, hjust = 1, size = 15),
        strip.text = element_text(size= 15)) +
  labs(x= "Soil Management Practice", y = "% Canopy Coverage")

```

```{r, echo= FALSE, include= FALSE}
sb_canopyherb <- ggplot(sb_canopy_herbCLD, aes(x= herb, y= crop_canopy*100, color= herb)) +
  geom_point(size= 3) +
  geom_errorbar(aes(ymin= lower.CL*100, ymax= upper.CL*100), width= .3) +
  geom_text(aes(label = .group), nudge_y = 4, size= 6) +
  geom_jitter(data = Soybean1 ,mapping = aes(y = crop_canopy*100), alpha = 0.2) +
  #coord_flip() +
  #facet_grid(year~location) +
  theme_bw() +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "none",
        plot.title = element_text(size = 25),
        #axis.title.x = element_text(size = 25),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 25),
        axis.text.y = element_text(size = 15),
        axis.text.x = element_text(size = 15)) +
  labs(x= "Herbicide Treatment", y= "% Canopy Coverage")
#Models below in comparison
```

## Separated analysis for each site-year
Only first model displayed
```{r,}
#Beginning of analysis for separate site-year models

arl19_soy_can = glmmTMB(crop_canopy~ soil*herb + (1|rep:site_crop_yr), data= (filter(Soybean1, site_crop_yr == "Arlington 2019")), beta_family(link="logit"))

Anova(arl19_soy_can)
#Soil and herb fixed effects significant
```

```{r, include = FALSE, echo = FALSE}
arl20_soy_can = glmmTMB(crop_canopy~ soil*herb + (1|rep:site_crop_yr), data= (filter(Soybean1, site_crop_yr == "Arlington 2020")), beta_family(link="logit"))

Anova(arl20_soy_can)
#nothing significant
```

```{r, echo= FALSE, include= FALSE}
arl19_soy_can_lsmeans= emmeans(arl19_soy_can ,~ soil, contr="pairwise", adjust="none", type="response")
arl19_soy_canCLD<- CLD(arl19_soy_can_lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

arl19_soy_canCLD <- as_tibble(arl19_soy_canCLD) %>%
  rename(crop_canopy = response) %>%
  mutate(soil= factor(soil, levels= order)) %>%
  arrange(soil)

arl19_soy_can_herblsmeans= emmeans(arl19_soy_can ,~ herb, contr="pairwise", adjust="none", type="response")
arl19_soy_canherbCLD<- CLD(arl19_soy_can_herblsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

arl19_soy_canherbCLD <- as_tibble(arl19_soy_canherbCLD) %>%
  rename(crop_canopy = response)%>%
  mutate(herb= factor(herb, levels= order2)) %>%
  arrange(herb)

arl20_soy_can_lsmeans= emmeans(arl20_soy_can ,~ soil, contr="pairwise", adjust="none", type="response")
arl20_soy_canCLD<- CLD(arl20_soy_can_lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

arl20_soy_canCLD <- as_tibble(arl20_soy_canCLD) %>%
  rename(crop_canopy = response) %>%
  mutate(soil= factor(soil, levels= order)) %>%
  arrange(soil)

arl20_soy_can_herblsmeans= emmeans(arl20_soy_can ,~ herb, contr="pairwise", adjust="none", type="response")
arl20_soy_canherbCLD<- CLD(arl20_soy_can_herblsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

arl20_soy_canherbCLD <- as_tibble(arl20_soy_canherbCLD) %>%
  rename(crop_canopy = response) %>%
  mutate(herb= factor(herb, levels= order2)) %>%
  arrange(herb)
```

```{r, include= FALSE, echo= FALSE}
lan19_soy_can = glmmTMB(crop_canopy~ soil*herb + (1|rep:site_crop_yr), data= (filter(Soybean1, site_crop_yr == "Lancaster 2019")), beta_family(link="logit"))

Anova(lan19_soy_can)
#Soil was significant

lan20_soy_can = glmmTMB(crop_canopy~ soil*herb + (1|rep:site_crop_yr), data= (filter(Soybean1, site_crop_yr == "Lancaster 2020")), beta_family(link="logit"))

Anova(lan20_soy_can)
#nothing significant
```

```{r, echo= FALSE, include= FALSE}
lan19_soy_can_lsmeans= emmeans(lan19_soy_can ,~ soil, contr="pairwise", adjust="none", type="response")
lan19_soy_canCLD<- CLD(lan19_soy_can_lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

lan19_soy_canCLD <- as_tibble(lan19_soy_canCLD) %>%
  rename(crop_canopy = response) %>%
  mutate(soil= factor(soil, levels= order)) %>%
  arrange(soil)

lan19_soy_can_herblsmeans= emmeans(lan19_soy_can ,~ herb, contr="pairwise", adjust="none", type="response")
lan19_soy_canherbCLD<- CLD(lan19_soy_can_herblsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

lan19_soy_canherbCLD <- as_tibble(lan19_soy_canherbCLD) %>%
  rename(crop_canopy = response) %>%
  mutate(herb= factor(herb, levels= order2)) %>%
  arrange(herb)

lan20_soy_can_lsmeans= emmeans(lan20_soy_can ,~ soil, contr="pairwise", adjust="none", type="response")
lan20_soy_canCLD<- CLD(lan20_soy_can_lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

lan20_soy_canCLD <- as_tibble(lan20_soy_canCLD) %>%
  rename(crop_canopy = response) %>%
  mutate(soil= factor(soil, levels= order)) %>%
  arrange(soil)

lan20_soy_can_herblsmeans= emmeans(lan20_soy_can ,~ herb, contr="pairwise", adjust="none", type="response")
lan20_soy_canherbCLD<- CLD(lan20_soy_can_herblsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

lan20_soy_canherbCLD <- as_tibble(lan20_soy_canherbCLD) %>%
  rename(crop_canopy = response) %>%
  mutate(herb= factor(herb, levels= order2)) %>%
  arrange(herb)
```

```{r, include= FALSE, echo= FALSE}
hav19_soy_can = glmmTMB(crop_canopy~ soil*herb + (1|rep:site_crop_yr), data= (filter(Soybean1, site_crop_yr == "Havelock 2019")), beta_family(link="logit"))

Anova(hav19_soy_can)
#Soil and herbicide fixed effects were significant

hav20_soy_can = glmmTMB(crop_canopy~ soil*herb + (1|rep:site_crop_yr), data= (filter(Soybean1, site_crop_yr == "Havelock 2020")), beta_family(link="logit"))

Anova(hav20_soy_can)
#Soil and herbicide fixed effects were significant
```

```{r, echo= FALSE, include= FALSE}
hav19_soy_can_lsmeans= emmeans(hav19_soy_can ,~ soil, contr="pairwise", adjust="none", type="response")
hav19_soy_canCLD<- CLD(hav19_soy_can_lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

hav19_soy_canCLD <- as_tibble(hav19_soy_canCLD) %>%
 rename(crop_canopy = response) %>%
  mutate(soil= factor(soil, levels= order)) %>%
  arrange(soil)

hav19_soy_can_herblsmeans= emmeans(hav19_soy_can ,~ herb, contr="pairwise", adjust="none", type="response")
hav19_soy_canherbCLD<- CLD(hav19_soy_can_herblsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

hav19_soy_canherbCLD <- as_tibble(hav19_soy_canherbCLD) %>%
 rename(crop_canopy = response) %>%
  mutate(herb= factor(herb, levels= order2)) %>%
  arrange(herb)

hav20_soy_can_lsmeans= emmeans(hav20_soy_can ,~ soil, contr="pairwise", adjust="none", type="response")
hav20_soy_canCLD<- CLD(hav20_soy_can_lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

hav20_soy_canCLD <- as_tibble(hav20_soy_canCLD) %>%
 rename(crop_canopy = response) %>%
  mutate(soil= factor(soil, levels= order)) %>%
  arrange(soil)

hav20_soy_can_herblsmeans= emmeans(hav20_soy_can ,~ herb, contr="pairwise", adjust="none", type="response")
hav20_soy_canherbCLD<- CLD(hav20_soy_can_herblsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

hav20_soy_canherbCLD <- as_tibble(hav20_soy_canherbCLD) %>%
 rename(crop_canopy = response) %>%
  mutate(herb= factor(herb, levels= order2)) %>%
  arrange(herb)
```

```{r, echo= FALSE, include= FALSE}
sbcanopy= data.frame()

arl19_soy_canCLD$location= c("Arlington")
arl20_soy_canCLD$location= c("Arlington")
lan19_soy_canCLD$location= c("Lancaster")
lan20_soy_canCLD$location= c("Lancaster")
hav19_soy_canCLD$location= c("Havelock")
hav20_soy_canCLD$location= c("Havelock")

arl19_soy_canCLD$year= c("2019")
arl20_soy_canCLD$year= c("2020")
lan19_soy_canCLD$year= c("2019")
lan20_soy_canCLD$year= c("2020")
hav19_soy_canCLD$year= c("2019")
hav20_soy_canCLD$year= c("2020")

sbcanopy<- rbind(arl19_soy_canCLD, sbcanopy)
sbcanopy<- rbind(arl20_soy_canCLD, sbcanopy)
sbcanopy<- rbind(lan19_soy_canCLD, sbcanopy)
sbcanopy<- rbind(lan20_soy_canCLD, sbcanopy)
sbcanopy<- rbind(hav19_soy_canCLD, sbcanopy)
sbcanopy<- rbind(hav20_soy_canCLD, sbcanopy)
```

```{r, echo= FALSE, include= FALSE}
sbcanopy2= data.frame()

arl19_soy_canherbCLD$location= c("Arlington")
arl20_soy_canherbCLD$location= c("Arlington")
lan19_soy_canherbCLD$location= c("Lancaster")
lan20_soy_canherbCLD$location= c("Lancaster")
hav19_soy_canherbCLD$location= c("Havelock")
hav20_soy_canherbCLD$location= c("Havelock")

arl19_soy_canherbCLD$year= c("2019")
arl20_soy_canherbCLD$year= c("2020")
lan19_soy_canherbCLD$year= c("2019")
lan20_soy_canherbCLD$year= c("2020")
hav19_soy_canherbCLD$year= c("2019")
hav20_soy_canherbCLD$year= c("2020")

sbcanopy2<- rbind(arl19_soy_canherbCLD, sbcanopy2)
sbcanopy2<- rbind(arl20_soy_canherbCLD, sbcanopy2)
sbcanopy2<- rbind(lan19_soy_canherbCLD, sbcanopy2)
sbcanopy2<- rbind(lan20_soy_canherbCLD, sbcanopy2)
sbcanopy2<- rbind(hav19_soy_canherbCLD, sbcanopy2)
sbcanopy2<- rbind(hav20_soy_canherbCLD, sbcanopy2)
```

```{r, echo= FALSE, include= FALSE}
sb_canopy1 <- ggplot(sbcanopy, aes(x= soil, y= crop_canopy*100, color= soil)) +
  geom_point(size= 3) +
  geom_errorbar(aes(ymin= lower.CL*100, ymax= upper.CL*100), width= .3) +
  geom_text(aes(label = .group), nudge_y = 10, size= 6) +
  geom_jitter(data = Soybean1 ,mapping = aes(y = crop_canopy*100), alpha = 0.2) +
  #coord_flip() +
  facet_grid(year~location) +
  theme_bw() +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "none",
        plot.title = element_text(size = 25),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25),
        axis.text.x = element_text(size = 15)) +
  labs(title = "Soybean Canopy")


```



```{r, echo= FALSE, include= FALSE}
sb_canopyherb1 <- ggplot(sbcanopy2, aes(x= herb, y= crop_canopy*100, color= herb)) +
  geom_point(size= 3) +
  geom_errorbar(aes(ymin= lower.CL*100, ymax= upper.CL*100), width= .3) +
  geom_text(aes(label = .group), nudge_y = 10, size = 6) +
  geom_jitter(data = Soybean1 ,mapping = aes(y = crop_canopy*100), alpha = 0.2) +
  #coord_flip() +
  facet_grid(year~location) +
  theme_bw() +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "none",
        plot.title = element_text(size = 25),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25),
        axis.text.x = element_text(size = 15)) +
  labs(title = "Soybean Canopy")

### end of analysis with separate models for each site-year
```



## Comparison of Canopy Coverage Analysis Approaches

### Separated
```{r, echo= FALSE, include= FALSE}
separated<- {sb_canopy + sb_canopyherb1}+
  plot_annotation(title = 'Soybean Canopy Coverage',
  subtitle = "Separated Analysis",
  theme = theme(plot.title = element_text(size = 30),
                plot.subtitle= element_text(size = 20))) & 
  theme(plot.title = element_text("serif"))
```

```{r, fig.height= 6, fig.width= 12}
separated
```



### Condensed
```{r, echo= FALSE, include= FALSE}
composite1<- {sb_canopyherb + sb_canopy}+
  #plot_annotation(title = 'Soybean Canopy Coverage',
                  #subtitle = 'V3-V4 Crop Stage') &
  theme(plot.title = element_text(size = 30),
                plot.subtitle= element_text(size = 20))
```

```{r, echo= FALSE, fig.height= 6, fig.width= 12}
composite1
```


Results are very similar for all site-years with the soybean canopy models.  


In the analysis with separate models for every site-year it seems weird that the control trt has
a lower canopy coverage compared to some of the herbicide trts at Havelock 2019. 

Personally I prefer the simplified analyses with one model for all site-years 
to keep things straightforward, especially for the herbicide results.


# Soybean Yield Analysis

```{r}
Soybean1 %>%
  ggplot(aes(x = soil, y = yield, color = location)) +
  geom_boxplot() +
  facet_grid(crop ~ year)
```
It appears as if there are similar trends and yield within a location across years.
I will proceed with performing a separate analysis for each location.

```{r}
sb_yield= lmer(yield~ location*soil*herb + (1|rep:site_crop_yr), data= Soybean1)

qqnorm(resid(sb_yield))

plot(sb_yield)

anova(sb_yield)
#Location and soil:Location interaction significant
```


### Arlington Analysis

Only going to show the first model. Repeated for each location
```{r, fig.show= 'hide'}
arl_sb_yield= lmer(yield~ soil*herb + (1|rep/year), data= (filter(Soybean1, location == "Arlington")))

qqnorm(resid(arl_sb_yield))

plot(arl_sb_yield)
#assumptions look good
```

```{r,}
anova(arl_sb_yield)
#nothing significant
```


### Lancaster Analysis
```{r, include= FALSE, echo= FALSE}
lan_sb_yield= lmer(yield~ soil*herb + (1|rep/year), data= (filter(Soybean1, location == "Lancaster")))

qqnorm(resid(lan_sb_yield))

plot(lan_sb_yield)
#assumptions look good
```

```{r}
anova(lan_sb_yield)
#Soil management fixed effect significant
```

### Havelock Analysis
```{r, include= FALSE, echo= FALSE}
hav_sb_yield= lmer(yield~ soil*herb + (1|rep/year), data= (filter(Soybean1, location == "Lincoln")))

qqnorm(resid(hav_sb_yield))

plot(hav_sb_yield)
#assumptions look good
```

```{r}
anova(hav_sb_yield)
#nothing significant
```


```{r, echo= FALSE, include= FALSE}
arl_sb_lsmeans<- lsmeans(arl_sb_yield, ~ soil, contr="pairwise", adjust="none")

arl_sb_cld <- cld(arl_sb_lsmeans$lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

arl_sb_cld <- as_tibble(arl_sb_cld) %>% 
  rename(yield = lsmean) %>%
  mutate(soil= factor(soil, levels= order)) %>%
  arrange(soil)
```

```{r, echo= FALSE, include= FALSE}
lan_sb_lsmeans<- lsmeans(lan_sb_yield, ~ soil, contr="pairwise", adjust="none")

lan_sb_cld <- cld(lan_sb_lsmeans$lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

lan_sb_cld <- as_tibble(lan_sb_cld) %>% 
  rename(yield = lsmean) %>%
  mutate(soil= factor(soil, levels= order)) %>%
  arrange(soil)
```

```{r, echo= FALSE, include= FALSE}
hav_sb_lsmeans<- lsmeans(hav_sb_yield, ~ soil, contr="pairwise", adjust="none")

hav_sb_cld <- cld(hav_sb_lsmeans$lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

hav_sb_cld <- as_tibble(hav_sb_cld) %>% 
  rename(yield = lsmean) %>%
  mutate(soil= factor(soil, levels= order)) %>%
  arrange(soil)
```

```{r, echo= FALSE, include= FALSE}
sbyield= data.frame()

arl_sb_cld$location= c("Arlington")
lan_sb_cld$location= c('Lancaster')
hav_sb_cld$location= c('Lincoln')

sbyield<- rbind(arl_sb_cld, sbyield)
sbyield<- rbind(lan_sb_cld, sbyield)
sbyield<- rbind(hav_sb_cld, sbyield)
```

```{r, echo= FALSE, include= FALSE}
ylab.text1= expression("Yield kg ha"^"-1")

y2<- ggplot(sbyield, aes(x= soil, y= yield, color= soil)) +
  geom_point(size= 3) +
  geom_errorbar(aes(ymin= lower.CL, ymax= upper.CL), width= .3) +
  geom_text(aes(label = .group), nudge_y = 800, size= 6) +
  geom_jitter(data = Soybean1 ,mapping = aes(y = yield), alpha = 0.2) +
  #coord_flip() +
  facet_grid(~location) +
  theme_bw() +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "none",
        plot.title = element_text(size = 25),
        axis.title.x = element_text(size= 25),
        axis.title.y = element_text(size = 25),
        axis.text.x = element_text(size = 15),
        strip.text = element_text(size= 15)) +
  labs(title = "Soybean Yield", x= "Soil Management Practice" , y= bquote('Yield'~(kg~ha^-1)))
```

## Soybean Yield Figure
```{r, fig.width= 12, fig.height=6}
y2
```
Lancaster had reduced yield for the cover crop soil management treatment. If you look back to
the analysis we did on stand counts you can see that this treatment had the lowest mean
early-season stand count at Lancaster both years.

Now to see if the yield components explain more.


# Soybean Yield Components

Based on issues with storage and seed damage in 2019, the only data I feel comfortable with
using for 2019 yield components are the pod counts. Even the seed counts seem suspiciously 
low for 2019 so seeds per pod is probably useless for 2019.

```{r}
#SoybeanComp1 %>%
  #ggplot(aes(x = Soil, y = Total, color = Soil)) +
  #geom_boxplot() +
  #facet_grid(~ Site_Yr)

#Plot of Soybean seed density (Hundred weight)
SoybeanComp1 %>%
  ggplot(aes(x = Soil, y = Hundred, color = Soil)) +
  geom_boxplot() +
  facet_grid(~ Site_Yr)

#Plot of Soybean seed counts
SoybeanComp %>%
  ggplot(aes(x = Soil, y = Seeds, color = Soil)) +
  geom_boxplot() +
  facet_grid(~ Site_Yr)
#All of the 2019 data is lower, I am not confident in using it.

#Plot of Soybean pods per plant
SoybeanComp %>%
  ggplot(aes(x = Soil, y = Pods.Plant, color = Soil)) +
  geom_boxplot() +
  facet_grid(~ Site_Yr)

#Plot of Soybean seeds per pod
SoybeanComp1 %>%
  ggplot(aes(x = Soil, y = Seeds.Pod, color = Soil)) +
  geom_boxplot() +
  facet_grid(~ Site_Yr)
```

I decided to only analyse models for the lancaster location. I have done more analysis
on these yield components looking at site-year but they aren't really useful if we can't
compare them to yield so I am not sharing them with this document. Remember,
we are trying to explain reduced yield from soil management with a cover crop at Lancaster.
I made models for all locations for fun though.

## Pods per plant 

```{r, fig.show='hide'}
lan_SBPods_Plant= lmer(Pods.Plant~Soil*Herb+ (1|Rep/Year) , data= (filter(SoybeanComp, Location == "Lancaster")))
qqnorm(resid(lan_SBPods_Plant))
plot(lan_SBPods_Plant)
#assumptions met

```

```{r,}
anova(lan_SBPods_Plant)
#nothing significant, soil almost significant
```

```{r,}
lan_pod_plantlsmeans<- lsmeans(lan_SBPods_Plant, ~ Soil, contr="pairwise", adjust="none")

lan_pod_plantcld <- cld(lan_pod_plantlsmeans$lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

lan_pod_plantcld <- as_tibble(lan_pod_plantcld) %>% 
  rename(Pods.Plant = lsmean) %>%
  mutate(Soil= factor(Soil, levels= order)) %>%
  arrange(Soil)

lan_pod_plantcld
```


## Seed density 

```{r, fig.show='hide'}
lan_dens= lmer(Hundred~Soil*Herb+ (1|Rep) , data= (filter(SoybeanComp1, Location == "Lancaster")))
qqnorm(resid(lan_dens))
plot(lan_dens)
#assumptions met

```

```{r,}
anova(lan_dens)
#nothing significant
```

```{r,}
lan_denslsmeans<- lsmeans(lan_dens, ~ Soil, contr="pairwise", adjust="none")

lan_denscld <- cld(lan_denslsmeans$lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

lan_denscld <- as_tibble(lan_denscld) %>% 
  rename(Hundred = lsmean) %>%
  mutate(Soil= factor(Soil, levels= order)) %>%
  arrange(Soil)

lan_denscld 
```

## Seed Counts 

```{r, fig.show='hide'}
lan_seeds= lmer(Seeds~Soil*Herb+ (1|Rep) , data= (filter(SoybeanComp1, Location == "Lancaster")))
qqnorm(resid(lan_seeds))
plot(lan_seeds)
#assumptions met

```

```{r}
anova(lan_seeds)
#nothing significant
```

```{r}
lan_seedslsmeans<- lsmeans(lan_seeds, ~ Soil, contr="pairwise", adjust="none")

lan_seedscld <- cld(lan_seedslsmeans$lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

lan_seedscld <- as_tibble(lan_seedscld) %>% 
  rename(Seeds = lsmean) %>%
  mutate(Soil= factor(Soil, levels= order)) %>%
  arrange(Soil)

lan_seedscld
```

## Soybean Seeds/Pod
```{r, fig.show='hide'}
lan_seeds_pod= lmer(Seeds.Pod~Soil*Herb+ (1|Rep) , data= (filter(SoybeanComp1, Location == "Lancaster")))
qqnorm(resid(lan_seeds_pod))
plot(lan_seeds_pod)
#assumptions met

```

```{r}
anova(lan_seeds_pod)
#nothing significant
```

```{r}
lan_seeds_podlsmeans<- lsmeans(lan_seeds_pod, ~ Soil, contr="pairwise", adjust="none")

lan_seeds_podcld <- cld(lan_seeds_podlsmeans$lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

lan_seeds_podcld <- as_tibble(lan_seeds_podcld) %>% 
  rename(Seeds = lsmean) %>%
  mutate(Soil= factor(Soil, levels= order)) %>%
  arrange(Soil)

lan_seeds_podcld
```

## Summary of yield components analysis

There weren't any useful yield components to explain what we saw at Lancaster with reduced yield
for the cover crop soil management treatment. I think the lowest mean stand count early-season is
our best explanation. Only having 2020 data for most yield components didn't help this and if you
look at the boxplots generated prior to the yield analysis the reduced yield appears to have 
occurred in 2019. There was also reduced canopy cover for this treatment in 2019.


# Cover crop analysis

## Biomass analysis

```{r}
SBCC %>%
  ggplot(aes(x = herb, y = biomass_kg, color = site)) +
  geom_boxplot() +
  facet_grid(crop ~ year)
```
based on separations between locations and across years, I think it is appropriate to look for significant differences between site-years and separate means accordingly in CC biomass models.

### Soybean
```{r, fig.show='hide'}
sb_cc_bio= lmer(biomass_kg~ site_crop_yr * herb + (1|site_crop_yr:rep), data=SBCC)

qqnorm(resid(sb_cc_bio))
plot(sb_cc_bio)
#assumptions for equal variance not met

sb_cc_bio1= lmer(sqrt(biomass_kg)~ site_crop_yr * herb + (1|site_crop_yr:rep), data=SBCC)

qqnorm(resid(sb_cc_bio1))
plot(sb_cc_bio1)
#assumption improved.

sb_cc_bio2= lmer(log(biomass_kg)~ site_crop_yr * herb + (1|site_crop_yr:rep), data=SBCC)

qqnorm(resid(sb_cc_bio2))
plot(sb_cc_bio2)
```


```{r}
anova(sb_cc_bio1)
#Site-year significant
```

There is no evidence that herbicide carryover reduced CC biomass.

```{r, include= FALSE, echo= FALSE}
sb_cc_lsmeans= emmeans(sb_cc_bio1 ,~ site_crop_yr, contr="pairwise", adjust="none", type="response")
sb_ccCLD<- CLD(sb_cc_lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

sb_ccCLD <- as_tibble(sb_ccCLD) %>%
 rename(biomass_kg = response)
```

```{r, include= FALSE, echo= FALSE}
Soy_CCbio= ggplot(sb_ccCLD, aes(x= site_crop_yr, y= biomass_kg, color= site_crop_yr)) +
  geom_point(size= 3) +
  geom_errorbar(aes(ymin= lower.CL, ymax= upper.CL), width= .3) +
  geom_text(aes(label = .group), nudge_y = 150, size= 5) +
  geom_jitter(data = SBCC ,mapping = aes(y = biomass_kg), alpha = 0.2) +
  coord_flip() +
  #facet_grid(year~location) +
  theme_bw() +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "none") +
  labs(title = "Soybean")
```


## Cover Crop Canopy



```{r, include= FALSE, echo= FALSE}
SBCC1 <- SBCC %>%
  filter(!is.na(canopy))
```

```{r}
SBCC1 %>%
  ggplot(aes(x = herb, y = canopy, color = site)) +
  geom_boxplot() +
  facet_grid(crop ~ year)
```

Based on differences between location across growing seasons I thought it was best to proceed with testing site-year as a fixed effect.


```{r, fig.show='hide'}
sb_cc_can= glmmTMB(canopy~ site_crop_yr*herb + (1|site_crop_yr:rep), data=SBCC, beta_family(link="logit"))
```

```{r}
Anova(sb_cc_bio)

glmmTMB:::Anova.glmmTMB(sb_cc_bio)
#Site-Year significant
```
There is no evidence that herbicid carryover influenced cover crop canopy


```{r, include= FALSE, echo= FALSE}
sb_cc_can_lsmeans= emmeans(sb_cc_can ,~ site_crop_yr, contr="pairwise", adjust="none", type="response")
sb_cc_canCLD<- CLD(sb_cc_can_lsmeans, alpha=0.05, Letters=letters, adjust="none", sort=TRUE, reverse=TRUE)

sb_cc_canCLD <- as_tibble(sb_cc_canCLD) %>%
 rename(canopy = response)
```


```{r, include= FALSE, echo= FALSE}
Soy_CCcan= ggplot(sb_cc_canCLD, aes(x= site_crop_yr, y= canopy*100, color= site_crop_yr)) +
  geom_point(size= 3) +
  geom_errorbar(aes(ymin= lower.CL*100, ymax= upper.CL*100), width= .3) +
  geom_text(aes(label = .group), nudge_y = 5, size= 5) +
  geom_jitter(data = SBCC1 ,mapping = aes(y = canopy*100), alpha = 0.2) +
  coord_flip() +
  #facet_grid(year~location) +
  theme_bw() +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "none") +
  labs(title = "Soybean")
```