---
title: "The influence of climate warming on flowering phenology in relation to historical annual and seasonal temperatures and plant functional traits - Revised 2023"
author: "Cole Geissler, Allison Davidson, Richard A. Niesenbaum"
date: "January 3, 2023"
output: 
  word_document:
    fig_height: 3
    fig_width: 5
  pdf_document:
    fig_height: 3
    fig_width: 5
  html_document:
    fig_height: 3
    fig_width: 5
---

```{r, setup, include=FALSE}
require(mosaic)
require(readr)
require(tidyr)
require(dplyr)
require(lme4)
library(ggpubr)
library(rstatix)
knitr::opts_chunk$set(
  tidy=FALSE,     # display code as typed
  size="small")   # slightly smaller font for code
```

#Influence of average yearly temperature on Collection Date | Analysis

```{r}
#Load phenology data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Load climate data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join data and set new name
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

lmer<-lmer(JulianCollectionDay ~ YearAvgTAVG + (1 + YearAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))
summary(lmer)

#Model
#Function to extract the slopes from the lmer( ) model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ YearAvgTAVG + (1 + YearAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta2, data=boot.samples, p=0.95)
```

#Influence of average yearly temperature on Collection Date | Plot (Figure 1)

```{r}
#Load phenology data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Load climate data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join data and set new name
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

gf_point(JulianCollectionDay ~ YearAvgTAVG, alpha=0.3, data=subset(new.data, CollectionYear > 1884)) %>%
  gf_abline(slope = -2.2693, intercept = 206.9257) +
  theme_classic() +
  xlab("Year Average Temperature (degrees C)") +
  ylab("Collection Date") +
  theme(legend.position="none")
```

#Influence of average yearly temperature on  Colection Date - Exclude species with >= 5 month flowering period | Analysis

```{r}
#Load phenology data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Load climate data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join data and set new name
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

lmer<-lmer(JulianCollectionDay ~ YearAvgTAVG + (1 + YearAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))
summary(lmer)

#Model
#Function to extract the slopes from the lmer( ) model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ YearAvgTAVG + (1 + YearAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod <= 4), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta2, data=boot.samples, p=0.95)
```

#Influence of winter average temperature on Collection Date | Analysis

```{r}
#Load phenology data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Load climate data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join data and set new name
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

lmer<-lmer(JulianCollectionDay ~ DJFAvgTAVG + (1 + DJFAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))
summary(lmer)

#Model
#Function to extract the slopes from the lmer( ) model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ DJFAvgTAVG + (1 + DJFAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta2, data=boot.samples, p=0.95)
```

#Influence of winter average temperature on Collection Date | Plot (Figure 2)

```{r}
#Load phenology data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Load climate data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join data and set new name
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

gf_point(JulianCollectionDay ~ DJFAvgTAVG, alpha=0.3, data=subset(new.data, CollectionYear > 1884)) %>%
  gf_abline(slope = -.1618, intercept = 181.3843) +
  theme_classic() +
  xlab("Winter Average Temperature (degrees C)") +
  ylab("Collection Date") +
  theme(legend.position = "none")
```

#Influence of winter average temperature on Collection Date - Exclude species with >= 5 month flowering period | Analysis

```{r}
#Load phenology data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Load climate data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join data and set new name
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

lmer<-lmer(JulianCollectionDay ~ DJFAvgTAVG + (1 + DJFAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))
summary(lmer)

#Model
#Function to extract the slopes from the lmer( ) model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ DJFAvgTAVG + (1 + DJFAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta2, data=boot.samples, p=0.95)
```

#Influence of average spring onset temperature on Collection Date | Analysis

```{r}
#Load phenology data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Load climate data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join data and set new name
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

lmer<-lmer(JulianCollectionDay ~ MAMAvgTAVG + (1 + MAMAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))
summary(lmer)

#Model
#Function to extract the slopes from the lmer( ) model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ MAMAvgTAVG + (1 + MAMAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta2, data=boot.samples, p=0.95)
```

#Influence of average spring onset temperature on Collection Date | Plot (Figure 3)

```{r}
#Load phenology data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Load climate data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join data and set new name
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

gf_point(JulianCollectionDay ~ MAMAvgTAVG, alpha=0.3, data=subset(new.data, CollectionYear > 1884)) %>%
  gf_abline(slope = -2.9283, intercept = 211.4333) +
  theme_classic() +
  xlab("Spring Onset Average Temperature (degrees C)") +
  ylab("Collection Date") +
  theme(legend.position="none")
```

#Influence of average spring onset temperature on Collection Date - Exclude species with >= 5 month flowering period | Analysis

```{r}
#Load phenology data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Load climate data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join data and set new name
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

lmer<-lmer(JulianCollectionDay ~ MAMAvgTAVG + (1 + MAMAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))
summary(lmer)

#Model
#Function to extract the slopes from the lmer( ) model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ MAMAvgTAVG + (1 + MAMAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta2, data=boot.samples, p=0.95)
```

#Native vs. Non-Native | Year Average Temperature (Table 2)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ YearAvgTAVG * NativeStatus + (1 + YearAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ YearAvgTAVG * NativeStatus + (1 + YearAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Native vs. Non-Native | Year Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ YearAvgTAVG * NativeStatus + (1 + YearAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ YearAvgTAVG * NativeStatus + (1 + YearAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Native vs. Non-Native | Winter Average Temperature (Table 2)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not being used in publication

#Regression
lmer(JulianCollectionDay ~ DJFAvgTAVG * NativeStatus + (1 + DJFAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ DJFAvgTAVG * NativeStatus + (1 + DJFAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Native vs. Non-Native | Winter Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not being used in publication

#Regression
lmer(JulianCollectionDay ~ DJFAvgTAVG * NativeStatus + (1 + DJFAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ DJFAvgTAVG * NativeStatus + (1 + DJFAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Native vs. Non-Native | Spring Onset Average Temperature (Table 2)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not being used in publication

#Regression
lmer(JulianCollectionDay ~ MAMAvgTAVG * NativeStatus + (1 + MAMAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ MAMAvgTAVG * NativeStatus + (1 + MAMAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Native vs. Non-Native | Spring Onset Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not being used in publication

#Regression
lmer(JulianCollectionDay ~ MAMAvgTAVG * NativeStatus + (1 + MAMAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ MAMAvgTAVG * NativeStatus + (1 + MAMAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Woody vs. Herbaceous | Year Average Temperature (Table 2 and Figure 4)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot
gf_point(JulianCollectionDay ~ YearAvgTAVG, alpha=0.3, shape=~GrowthHabit, data=subset(new.data, CollectionYear > 1884)) %>%
  gf_abline(slope = -1.681, intercept = 205.980, linewidth=0.6) %>% #Herbaceous line
  gf_abline(slope = -4.533, intercept = 215.321, linewidth=0.6, linetype = "dashed") + #Woody line
  theme_classic() +
  xlab("Year Average Temperature (degrees C)") +
  ylab("Collection Date") +
  theme(legend.position="none") +
  scale_shape_manual(values=c(16, 17))

#Regression
lmer(JulianCollectionDay ~ YearAvgTAVG * GrowthHabit + (1 + YearAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ YearAvgTAVG * GrowthHabit + (1 + YearAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Woody vs. Herbaceous | Year Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot
gf_point(JulianCollectionDay ~ YearAvgTAVG, alpha=0.3, shape~GrowthHabit, data=subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5)) %>%
  gf_abline(slope = -1.681, intercept = 210.410, linewidth=0.6) %>% #Herbaceous line
  gf_abline(slope = -4.515, intercept = 215.135, linewidth=0.6, linetype = "dashed") + #Woody line
  theme_classic() +
  xlab("Year Average Temperature (degrees C)") +
  ylab("Collection Date") +
  theme(legend.position="none") +
  scale_shape_manual(values=c(16, 17))

#Regression
lmer(JulianCollectionDay ~ YearAvgTAVG * GrowthHabit + (1 + YearAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ YearAvgTAVG * GrowthHabit + (1 + YearAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Woody vs. Herbaceous | Winter Average Temperature (Table 2)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ DJFAvgTAVG * GrowthHabit + (1 + DJFAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ DJFAvgTAVG * GrowthHabit + (1 + DJFAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Woody vs. Herbaceous | Winter Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ DJFAvgTAVG * GrowthHabit + (1 + DJFAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ DJFAvgTAVG * GrowthHabit + (1 + DJFAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Woody vs. Herbaceous | Spring Onset Average Temperature (Table 2)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ MAMAvgTAVG * GrowthHabit + (1 + MAMAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ MAMAvgTAVG * GrowthHabit + (1 + MAMAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Woody vs. Herbaceous | Spring Onset Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ MAMAvgTAVG * GrowthHabit + (1 + MAMAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ MAMAvgTAVG * GrowthHabit + (1 + MAMAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Dry vs. Fleshy Fruit | Year Average Temperature (Table 2)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ YearAvgTAVG * FruitType + (1 + YearAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ YearAvgTAVG * FruitType + (1 + YearAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Dry vs. Fleshy Fruit | Year Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ YearAvgTAVG * FruitType + (1 + YearAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ YearAvgTAVG * FruitType + (1 + YearAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Dry vs. Fleshy Fruit | Winter Average Temperature (Table 2)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot
gf_point(JulianCollectionDay ~ DJFAvgTAVG, alpha=0.3, shape~FruitType, data=subset(new.data, CollectionYear > 1884)) %>%
  gf_abline(slope = -0.4553, intercept = 185.7856, linewidth=0.6) %>% #Dry line
  gf_abline(slope = 0.712, intercept = 165.4542, linewidth=0.6, linetype = "dashed") + #Fleshy line
  theme_classic() +
  xlab("Winter Average Temperature (degrees C)") +
  ylab("Collection Date") +
  theme(legend.position="none") +
  scale_shape_manual(values=c(16, 17))

#Regression
lmer(JulianCollectionDay ~ DJFAvgTAVG * FruitType + (1 + DJFAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ DJFAvgTAVG * FruitType + (1 + DJFAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Dry vs. Fleshy Fruit | Winter Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot
gf_point(JulianCollectionDay ~ DJFAvgTAVG, alpha=0.3, shape~FruitType, data=subset(new.data, CollectionYear > 1884)) %>%
  gf_abline(slope = -0.5048, intercept = 186.1987, linewidth=0.6) %>% #Dry line
  gf_abline(slope = 0.5772, intercept = 165.3987, linewidth=0.6, linetype = "dashed") + #Fleshy line
  theme_classic() +
  xlab("Winter Average Temperature (degrees C)") +
  ylab("Collection Date") +
  theme(legend.position="none") +
  scale_shape_manual(values=c(16, 17))

#Regression
lmer(JulianCollectionDay ~ DJFAvgTAVG * FruitType + (1 + DJFAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ DJFAvgTAVG * FruitType + (1 + DJFAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Dry vs. Fleshy Fruit | Spring Onset Average Temperature (Table 2)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ MAMAvgTAVG * FruitType + (1 + MAMAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ MAMAvgTAVG * FruitType + (1 + MAMAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Dry vs. Fleshy Fruit | Spring Onset Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ MAMAvgTAVG * FruitType + (1 + MAMAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ MAMAvgTAVG * FruitType + (1 + MAMAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Spring- vs. Summer-Blooming | Year Average Temperature (Table 2 and Figure 5)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot
gf_point(JulianCollectionDay ~ YearAvgTAVG, alpha=0.3, shape=~Seasonality, data=subset(new.data, CollectionYear > 1884)) %>%
  gf_abline(slope = -3.284, intercept = 198.472, linewidth=0.6) %>% #Spring line
  gf_abline(slope = -1.014, intercept = 217.668, linewidth=0.6, linetype = "dashed") + #Summer line
  theme_classic() +
  xlab("Year Average Temperature (degrees C)") +
  ylab("Collection Date") +
  theme(legend.position="none") +
  scale_shape_manual(values=c(16, 17))

#Regression
lmer(JulianCollectionDay ~ YearAvgTAVG * Seasonality + (1 + YearAvgTAVG|Species), data = subset(new.data, CollectionYear >
1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ YearAvgTAVG * Seasonality + (1 + YearAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Spring- vs. Summer-Blooming | Year Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot
gf_point(JulianCollectionDay ~ YearAvgTAVG, alpha=0.3, shape~Seasonality, data=subset(new.data, CollectionYear > 1884)) %>%
  gf_abline(slope = -3.953, intercept = 203.235, linewidth=0.6) %>% #Spring line
  gf_abline(slope = -1.295, intercept = 219.989, linewidth=0.6, linetype = "dashed") + #Summer line
  theme_classic() +
  xlab("Year Average Temperature (degrees C)") +
  ylab("Collection Date") +
  theme(legend.position="none") +
  scale_shape_manual(values=c(16, 17))

#Regression
lmer(JulianCollectionDay ~ YearAvgTAVG * Seasonality + (1 + YearAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ YearAvgTAVG * Seasonality + (1 + YearAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Spring- vs. Summer-Blooming | Winter Average Temperature (Table 2)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ DJFAvgTAVG * Seasonality + (1 + DJFAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ DJFAvgTAVG * Seasonality + (1 + DJFAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Spring- vs. Summer-Blooming | Winter Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ DJFAvgTAVG * Seasonality + (1 + DJFAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ DJFAvgTAVG * Seasonality + (1 + DJFAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Spring- vs. Summer-Blooming | Spring Onset Average Temperature (Table 2)

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ MAMAvgTAVG * Seasonality + (1 + MAMAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ MAMAvgTAVG * Seasonality + (1 + MAMAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```

#Spring- vs. Summer-Blooming | Spring Onset Average Temperature - Exclude species with >= 5 month flowering period

```{r}
#Enter Phenology Data
library(readr)
PhenologyData <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/PhenologyData.csv")

#Enter Climate Data
library(readr)
ClimateDataRevised <- read_csv("C:/Users/geiss/Desktop/HerbariumResearch/DataAndAnalysisFiles/ClimateDataRevised.csv")

#Join and Rename Data
new.data=(dplyr::left_join(PhenologyData, ClimateDataRevised, by = "UniqueYear"))

#Plot not used in publication

#Regression
lmer(JulianCollectionDay ~ MAMAvgTAVG * Seasonality + (1 + MAMAvgTAVG|Species), data = subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5))

#BELOW CODE FOR CI
#Function to extract the slopes from the lmer() model
mySlopes <- function(.) {
  c(beta = getME(., "beta"))
}

#Initialize i and form a space to store all of our slopes
i=1
boot.samples = NULL

#Run a for loop to do the same thing as do(1000) command does
for(i in 1:1000) {
  boot.samples = rbind(boot.samples, mySlopes(lmer(JulianCollectionDay ~ MAMAvgTAVG * Seasonality + (1 + MAMAvgTAVG|Species), data = resample(subset(new.data, CollectionYear > 1884 & FloweringPeriod < 5), groups = Species))))
}

head(boot.samples)
boot.samples = data.frame(boot.samples)
cdata(~beta4, data=boot.samples, p=0.95)
```