---
title: "germination"
output: html_document
date: "2025-01-03"
---

```{r setup, include=FALSE}
library(ggh4x)
library(forcats)
library(ggplot2)
library(stringr)
library(dplyr)
library(reshape)
library(pairwiseAdonis)
library(reshape2)
library(drc)
library(vegan)
library(lubridate)
library(plotly)
library(ggcorrplot)
library(patchwork)
library(data.table)
library(lme4)
library(MuMIn)
library(xtable)
library(lmerTest)
library(tidyverse)
library(ggpubr)
library(rstatix)
library(emmeans)
library(gridExtra) 
library(rstatix)
setwd("~/PhD/Germination project/publishing_data_and_code")
speciesinfo<-read.csv2("speciesinfo.csv", fileEncoding = "ISO-8859-1")
plotdata<-read.csv2("plotdata.csv", fileEncoding = "ISO-8859-1")
plotdata$plot<-paste(plotdata$Block,plotdata$Plot,sep="")
plotdata$Site<-tolower(plotdata$Site)
```

```{r filter on viability}
germ_indoor<-read.csv2("Germination_indoors.csv",sep=";")

germ_indoor<-germ_indoor%>%
  pivot_longer(cols = c(5:194), names_to = "variable", values_to = "value")%>%
  na.omit()%>%
  mutate(variable = gsub("^X", "", variable),
         variable=round(as.numeric(variable))) %>% 
  na.omit() %>%
  group_by(Code) %>%
  mutate(cumulative_sum = cumsum(value)) %>%
  mutate(viability=(max(cumulative_sum)/Seeds)*100)%>%
  dplyr::select(-variable,-value,-cumulative_sum)%>%
  distinct()%>%
  as.data.frame()

levels(germ_indoor$Temp) <- c("10°C", "15°C","20°C")

species_filter_viability<-germ_indoor%>%
  ungroup()%>%
  dplyr::select(-Code, -Temp, -Seeds)%>%
  group_by(Species)%>%
  mutate(viability=mean(viability))%>%
  distinct()%>%
  as.data.frame()%>%
  mutate(Germination_rate_over_20=ifelse(viability>20,"Yes","No"))

species_filter_viability_list<-species_filter_viability%>%filter(viability>20)%>%dplyr::select(Species)%>%as.vector()

### check if germination rate does not differ significantly between final group selection

t_test<-germ_indoor%>%
  inner_join(speciesinfo)%>%
  dplyr::select(Species,Speed,Temp,viability)%>%
  filter(Temp==10)%>%
  group_by(Species)%>%
  mutate(viability=mean(viability))%>%
  distinct(Speed,viability)%>%
  ungroup()

t_test%>%filter(Speed=="Slow")%>%dplyr::select(viability)%>%
  mutate(mean=mean(viability),
         sd=sd(viability))%>%
  distinct(mean,sd)

t_test%>%filter(Speed=="Fast")%>%dplyr::select(viability)%>%
  mutate(mean=mean(viability),
         sd=sd(viability))%>%
  distinct(mean,sd)

t.test(t_test%>%filter(Speed=="Slow")%>%dplyr::select(viability),t_test%>%filter(Speed=="Fast")%>%dplyr::select(viability))
```

```{r germination curves indoors }
germ_indoor<-read.csv2("Germination_indoors.csv", check.names = T,sep=";")

germ_indoor<-germ_indoor%>%
  pivot_longer(cols = c(5:ncol(germ_indoor)), names_to = "variable", values_to = "value")%>%
  na.omit()%>%
  mutate(variable = gsub("^X", "", variable),
         variable=round(as.numeric(variable)))%>%
  filter(Species%in%species_filter_viability_list$Species) # Filter out species with low viability

total_plate_seeds<- 
germ_indoor%>%
  distinct(Code,Species,Temp,Seeds)%>%
  group_by(Species,Temp)%>%
  mutate(Total_seeds=sum(Seeds))%>%
  ungroup()%>%
  dplyr::select(-Code,-Seeds)%>%
  distinct()

germ_indoor<-germ_indoor %>%
  dplyr::select(-Code, -Seeds) %>%
  inner_join(total_plate_seeds) %>%
  na.omit() %>%
  group_by(Species, Temp, variable, Total_seeds) %>%
  summarise(value = sum(as.numeric(value)), .groups = "drop") %>%
  group_by(Species, Temp, Total_seeds) %>%
  mutate(germinated = cumsum(value), Code = paste(Species, Temp)) %>%
  ungroup() %>%
  as.data.frame()

## Time-to-event models ##

T50_indoor <- data.frame()

Codes<-unique(germ_indoor%>%
                dplyr::select(Code)
              )

for (i in Codes$Code) {
  temp <- germ_indoor[germ_indoor$Code== i, ]
  
  timeBef <- c(0, temp$variable)
  timeAf <- c(temp$variable, Inf)
  Counts <- c(temp$value, unique(temp$Total_seeds) - sum(temp$value))
  
  df <- data.frame(timeBef, timeAf, Counts)
  
  # Use tryCatch to handle potential errors
  modTE <- tryCatch({
    drm(Counts ~ timeBef + timeAf, data = df, 
        fct = LN.3(), type = "event")
  }, error = function(e) {
    # If an error occurs, return NULL and skip to the next iteration
    return(NULL)
  })
  
  # Check if modTE is NULL (i.e., an error occurred)
  if (!is.null(modTE)) {
    # If no error, store the results
    time_range <- seq(0,1300,2)
    T50_indoor <- rbind(T50_indoor,
                 data.frame(
                   Code= i,
                   Species = unique(temp$Species),
                   Temp = unique(temp$Temp),
                   Speed_T50 = summary(modTE)[[3]][3, 1],
                   ES = summary(modTE)[[3]][3, 2],
                   Time = time_range,
                   Proportion = predict(modTE, newdata = data.frame(timeBef = time_range), type = "response"),
                   Proportion_T50 = predict(modTE, 
                                        newdata = data.frame(timeBef = summary(modTE)[[3]][3, 1]), 
                                        type = "response")
                 ))

  } else {
    # If an error occurred, add NA for that iteration
    T50_indoor <- rbind(T50_indoor,
                 data.frame(
                   Code= i,
                   Species = unique(temp$Species),
                   Temp = unique(temp$Temp),
                   Speed_T50 = NA,
                   ES = NA,
                   Time = NA,
                   Proportion = NA,
                   Proportion_T50 = NA
                 ))
  }
}

head(T50_indoor)
T50_indoor%>%
  filter(Temp=="10"|Species=="Plantago media"&Temp=="15")%>%
  dplyr::select(Species,Speed_T50)%>%
  distinct()%>%
  na.omit()%>%
  arrange(Speed_T50)%>%
  mutate(Speed_T50=round(Speed_T50, digits = 0))

T50_indoor%>% # Find mean and median of T50 values
  na.omit()%>%
  filter(Temp!=20)%>%
  mutate(mean=mean(Speed_T50))%>%
  mutate(median=median(Speed_T50))%>%
  distinct(mean,median)

T50_indoor<-T50_indoor%>%
  inner_join(speciesinfo)

lmer_indoor_germ<-T50_indoor%>%
  na.omit()%>%
  filter(Temp!=20)%>%
  group_by(Code)%>%
  mutate(Proportion=max(Proportion))%>%
  ungroup()%>%
  dplyr::select(-ES,-Code,-Time)%>%
  distinct()%>%
  as.data.frame()

model_speed<-lmer(data=lmer_indoor_germ,Speed_T50~Temp+Speed+(1|Species))
anova(model_speed)

model_speed<-lmer(data=lmer_indoor_germ,Proportion ~Temp+Speed+(1|Species))
anova(model_speed)

## Create model for all species combined ##

total_plate_seeds_all<-
  total_plate_seeds%>%
  inner_join(speciesinfo)%>%
  filter(Temp!=20)%>%
  group_by(Speed,Temp)%>%
  mutate(Total_seeds=sum(Total_seeds))%>%
  distinct(Total_seeds)

all_indoor<-germ_indoor%>%
  ungroup()%>%
  inner_join(speciesinfo)%>%
  filter(Temp!=20)%>%
  dplyr::select(Speed,Temp,variable,value)%>%
  group_by(Speed,variable)%>%
  mutate(value=sum(value))%>%
  ungroup()%>%
  distinct()%>%
  arrange(Speed,Temp,variable)%>%
  group_by(Speed,Temp)%>%
  mutate(germinated=cumsum(value))%>%
  as.data.frame()%>%
  inner_join(total_plate_seeds_all)

all_slow_10<-all_indoor%>%filter(Speed=="Slow"&Temp==10)

timeBef <- c(0, all_slow_10$variable)
timeAf <- c(all_slow_10$variable, Inf)
Counts <- c(all_slow_10$value, unique(all_slow_10$Total_seeds) - sum(all_slow_10$value))
  
df <- data.frame(timeBef, timeAf, Counts)
  
modTE_Slow_10 <-drm(Counts ~ timeBef + timeAf, data = df, 
        fct = LN.3(), type = "event")

time_range <- seq(0,1300,2)

predictions <- as.data.frame(predict(modTE_Slow_10, 
                       newdata = data.frame(timeBef = time_range), 
                       type = "response", se.fit = TRUE))

Slow_curve_10<-data.frame(Speed = "Slow",
          Temp=10,
          ES = predictions$SE,
          Time = time_range,
          Proportion = predictions$Prediction,
          Proportion_lower = predictions$Prediction - 1.96 * predictions$SE,
          Proportion_upper = predictions$Prediction + 1.96 * predictions$SE)

####

all_fast_10<-all_indoor%>%filter(Speed=="Fast"&Temp==10)

timeBef <- c(0, all_fast_10$variable)
timeAf <- c(all_fast_10$variable, Inf)
Counts <- c(all_fast_10$value, unique(all_fast_10$Total_seeds) - sum(all_fast_10$value))
  
df <- data.frame(timeBef, timeAf, Counts)
  
modTE_Fast_10 <-drm(Counts ~ timeBef + timeAf, data = df, 
        fct = LN.3(), type = "event")

time_range <- seq(0,1300,2)

predictions <- as.data.frame(predict(modTE_Fast_10, 
                       newdata = data.frame(timeBef = time_range), 
                       type = "response", se.fit = TRUE))

Fast_curve_10<-data.frame(Speed = "Fast",
          Temp=10,
          ES = predictions$SE,
          Time = time_range,
          Proportion = predictions$Prediction,
          Proportion_lower = predictions$Prediction - 1.96 * predictions$SE,
          Proportion_upper = predictions$Prediction + 1.96 * predictions$SE)


##
all_slow_15<-all_indoor%>%filter(Speed=="Slow"&Temp==15)

timeBef <- c(0, all_slow_15$variable)
timeAf <- c(all_slow_15$variable, Inf)
Counts <- c(all_slow_15$value, unique(all_slow_15$Total_seeds) - sum(all_slow_15$value))
  
df <- data.frame(timeBef, timeAf, Counts)
  
modTE_Slow_15 <-drm(Counts ~ timeBef + timeAf, data = df, 
        fct = LN.3(), type = "event")

time_range <- seq(0,1300,2)

predictions <- as.data.frame(predict(modTE_Slow_15, 
                       newdata = data.frame(timeBef = time_range), 
                       type = "response", se.fit = TRUE))

Slow_curve_15<-data.frame(Speed = "Slow",
          Temp=15,
          ES = predictions$SE,
          Time = time_range,
          Proportion = predictions$Prediction,
          Proportion_lower = predictions$Prediction - 1.96 * predictions$SE,
          Proportion_upper = predictions$Prediction + 1.96 * predictions$SE)


####

all_fast_15<-all_indoor%>%filter(Speed=="Fast"&Temp==15)

timeBef <- c(0, all_fast_15$variable)
timeAf <- c(all_fast_15$variable, Inf)
Counts <- c(all_fast_15$value, unique(all_fast_15$Total_seeds) - sum(all_fast_15$value))
  
df <- data.frame(timeBef, timeAf, Counts)
  
modTE_Fast_15 <-drm(Counts ~ timeBef + timeAf, data = df, 
        fct = LN.3(), type = "event")

time_range <- seq(0,1300,2)

predictions <- as.data.frame(predict(modTE_Fast_15, 
                       newdata = data.frame(timeBef = time_range), 
                       type = "response", se.fit = TRUE))

Fast_curve_15<-data.frame(Speed = "Fast",
          Temp=15,
          ES = predictions$SE,
          Time = time_range,
          Proportion = predictions$Prediction,
          Proportion_lower = predictions$Prediction - 1.96 * predictions$SE,
          Proportion_upper = predictions$Prediction + 1.96 * predictions$SE)


combined_curves<-Fast_curve_10%>%
  bind_rows(Fast_curve_15)%>%
  bind_rows(Slow_curve_10)%>%
  bind_rows(Slow_curve_15)
  
### Create plots ###
extra_rows <- germ_indoor %>%
  distinct(Species, Temp, Total_seeds, Code) %>%
  mutate(variable = 0, value = 0, germinated = 0)

# Bind the new rows to the original dataframe
germ_indoor_zeros <- bind_rows(germ_indoor, extra_rows) %>%
  arrange(Species, Temp, Total_seeds, variable)  # Arrange the data nicely

raw_data_indoors<-ggplot(germ_indoor_zeros,aes(variable,germinated/Total_seeds,col=as.factor(Temp)))+
  geom_line()+
  geom_point()+
  facet_wrap(~Species)+
  xlab("Hours")+
  ylab("Proportion of seed germinated")+ 
  theme( # Remove x-axis ticks
        panel.grid.major.x = element_blank(), # Remove major gridlines
        panel.grid.minor.x = element_blank(), # Remove minor gridlines
        panel.border = element_rect(color = "black", fill = NA, size = 1),
        strip.background = element_blank(),
        panel.background = element_blank(),
        legend.box.background = element_blank(),
        legend.position = "bottom")+
  labs(color = "Temperature")


time_to_event_invidivual_species_indoor <- ggplot(T50_indoor%>%
  filter(Temp!=20)) +
  geom_line(aes(x = Time, y = Proportion, col = as.factor(Temp))) +
  facet_wrap(~Speed + Species) +
  labs(x = "Time (hours)", 
       y = "Proportion germinated") +
  geom_vline(xintercept = 230.5912, linetype = "dashed") +
  geom_point(aes(x = Speed_T50, y = Proportion_T50, col = as.factor(Temp))) +
  xlim(-20, 1300)+
  #ylim(0,1)+
  theme( # Remove x-axis ticks
        panel.grid.major.x = element_blank(), # Remove major gridlines
        panel.grid.minor.x = element_blank(), # Remove minor gridlines
        panel.border = element_rect(color = "black", fill = NA, size = 1),
        strip.background = element_blank(),
        panel.background = element_blank(),
        legend.box.background = element_blank(),
        legend.position = "bottom")+
  labs(color = "Temperature")

time_to_event_groups_indoor <- ggplot() +
  labs(title = "Indoor", 
       x = "Time (hours)", 
       y = "Proportion germinated") +
   geom_line(data=combined_curves, aes(Time,Proportion,col=Speed, linetype = as.factor(Temp)),size=1)+
  geom_point(data=germ_indoor%>%inner_join(speciesinfo)%>%filter(Temp!=20),aes(variable,germinated/Total_seeds,col=Speed, shape=as.factor(Temp)),position = position_dodge(2),alpha=.3)+  
  scale_color_manual(values = c("Fast" = "#E69F00",
                               "Slow" = "#009E73"))+
  labs(color = "Speed\ngroup",
       shape = "Temperature",
       linetype = "Temperature")+
  xlim(-20, 1300)+ 
  theme( # Remove x-axis ticks
        panel.grid.major.x = element_blank(), # Remove major gridlines
        panel.grid.minor.x = element_blank(), # Remove minor gridlines
        panel.border = element_rect(color = "black", fill = NA, size = 1),
        strip.background = element_blank(),
        panel.background = element_blank(),
        legend.box.background = element_blank())


raw_data_indoors
time_to_event_invidivual_species_indoor
time_to_event_groups_indoor

#ggsave(file.path("C:/Users/tava0008/OneDrive - Umeå universitet/Documents/PhD/Germination project/latex/indoor_germination_raw.png"), plot = raw_data_indoors, width =9, height =6, units = "in", dpi = 300)

#ggsave(file.path("C:/Users/tava0008/OneDrive - Umeå universitet/Documents/PhD/Germination project/latex/indoor_germination_individual.png"), plot = time_to_event_invidivual_species_indoor, width = 7.5, height =6, units = "in", dpi = 300)
```

```{r germination curves outdoors}
germ_tests<-fread("Germination_outdoors.csv")
germ_tests<-as.data.frame(germ_tests[,-c(1,3,4,5)])

total_plot_seeds<-
  germ_tests%>%
  group_by(Species,start)%>%
  mutate(Total_seeds=sum(sown))%>%
  ungroup()%>%
  distinct(Species,start,Total_seeds)%>%
  as.data.frame()
  
germ_outdoor<-germ_tests %>%
  mutate(across(4:12, as.integer))%>%
  pivot_longer(cols=c(4:12), names_to="variable", values_to="value")%>%
  mutate(sowing.moment=ifelse(start=="2021-07-01","late","early"),
         variable=as.numeric(variable),
         variable=ifelse(sowing.moment=="late",variable-21,variable))%>%
  filter(variable>0&variable<40) %>%
  na.omit()%>%
  inner_join(total_plot_seeds)%>%
  group_by(Total_seeds,Species,sowing.moment,variable)%>%
  mutate(value=sum(value))%>%
  ungroup()%>%
  distinct(Species,sowing.moment,Total_seeds,variable,value)%>%
  as.data.frame()%>%
  mutate(Code=paste(Species,sowing.moment))%>%
  group_by(Code)%>%
  mutate(germinated=cumsum(value))
  
## Time-to-event models ##

T50_outdoor <- data.frame()

Codes<-unique(germ_outdoor%>%
                dplyr::select(Code))

for (i in Codes$Code) {
  temp <- germ_outdoor[germ_outdoor$Code== i, ]
  
  timeBef <- c(0, temp$variable)
  timeAf <- c(temp$variable, Inf)
  Counts <- c(temp$value, unique(temp$Total_seeds) - sum(temp$value))
  
  df <- data.frame(timeBef, timeAf, Counts)
  
  # Use tryCatch to handle potential errors
  modTE <- tryCatch({
    drm(Counts ~ timeBef + timeAf, data = df, 
        fct = LN.3(), type = "event")
  }, error = function(e) {
    # If an error occurs, return NULL and skip to the next iteration
    return(NULL)
  })
  
  # Check if modTE is NULL (i.e., an error occurred)
  if (!is.null(modTE)) {
    # If no error, store the results
    time_range <- seq(0,60,1)
    T50_outdoor <- rbind(T50_outdoor,
                 data.frame(
                   Code= i,
                   Species = unique(temp$Species),
                   sowing.moment = unique(temp$sowing.moment),
                   Speed_T50 = summary(modTE)[[3]][3, 1],
                   ES = summary(modTE)[[3]][3, 2],
                   Time = time_range,
                   Proportion = predict(modTE, newdata = data.frame(timeBef = time_range), type = "response"),
                   Proportion_T50 = predict(modTE, 
                                        newdata = data.frame(timeBef = summary(modTE)[[3]][3, 1]), 
                                        type = "response")
                 ))

  } else {
    # If an error occurred, add NA for that iteration
    T50_outdoor <- rbind(T50_outdoor,
                 data.frame(
                   Code= i,
                   Species = unique(temp$Species),
                   sowing.moment = unique(temp$sowing.moment),
                   Speed_T50 = NA,
                   ES = NA,
                   Time = NA,
                   Proportion = NA,
                   Proportion_T50 = NA
                 ))
  }
}

head(T50_outdoor)

###

T50_outdoor%>% # Find mean and median of T50 values
  na.omit()%>%
  mutate(mean=mean(Speed_T50))%>%
  mutate(median=median(Speed_T50))%>%
  distinct(mean,median)

T50_outdoor<-T50_outdoor%>%
  inner_join(speciesinfo)

total_plot_seeds_all<-
  total_plot_seeds%>%
  inner_join(speciesinfo)%>%
  group_by(Speed,start)%>%
  mutate(Total_seeds=sum(Total_seeds))%>%
  distinct(start,Total_seeds)%>%
  mutate(sowing.moment=ifelse(start=="2021-07-01","late","early"))

all_outdoor<-germ_outdoor%>%
  ungroup()%>%
  inner_join(speciesinfo)%>%
  dplyr::select(Speed,sowing.moment,variable,value)%>%
  group_by(Speed,sowing.moment,variable)%>%
  mutate(value=sum(value))%>%
  ungroup()%>%
  distinct()%>%
  arrange(Speed,sowing.moment,variable)%>%
  group_by(Speed,sowing.moment)%>%
  mutate(germinated=cumsum(value))%>%
  as.data.frame()%>%
  inner_join(total_plot_seeds_all)

# Slow, early
all_slow_early<-all_outdoor%>%filter(Speed=="Slow"&sowing.moment=="early")

timeBef <- c(0, all_slow_early$variable)
timeAf <- c(all_slow_early$variable, Inf)
Counts <- c(all_slow_early$value, unique(all_slow_early$Total_seeds) - sum(all_slow_early$value))
  
df <- data.frame(timeBef, timeAf, Counts)
  
modTE_Slow_early <-drm(Counts ~ timeBef + timeAf, data = df, 
        fct = LN.3(), type = "event")

time_range <- seq(0,35,1)

predictions <- as.data.frame(predict(modTE_Slow_early, 
                       newdata = data.frame(timeBef = time_range), 
                       type = "response", se.fit = TRUE))

Slow_early_curve<-data.frame(Speed = "Slow",
          Speed_T50 = summary(modTE_Slow_early)[[3]][3, 1],
          sowing.moment="early",
          ES = predictions$SE,
          Time = time_range,
          Proportion = predictions$Prediction,
          Proportion_lower = predictions$Prediction - 1.96 * predictions$SE,
          Proportion_upper = predictions$Prediction + 1.96 * predictions$SE)

#Slow, late
all_slow_late<-all_outdoor%>%filter(Speed=="Slow"&sowing.moment=="late")

timeBef <- c(0, all_slow_late$variable)
timeAf <- c(all_slow_late$variable, Inf)
Counts <- c(all_slow_late$value, unique(all_slow_late$Total_seeds) - sum(all_slow_late$value))
  
df <- data.frame(timeBef, timeAf, Counts)
  
modTE_Slow_late <-drm(Counts ~ timeBef + timeAf, data = df, 
        fct = LN.3(), type = "event")

time_range <- seq(0,35,1)

predictions <- as.data.frame(predict(modTE_Slow_late, 
                       newdata = data.frame(timeBef = time_range), 
                       type = "response", se.fit = TRUE))

Slow_late_curve<-data.frame(Speed = "Slow",
          Speed_T50 = summary(modTE_Slow_late)[[3]][3, 1],
          sowing.moment="late",
          ES = predictions$SE,
          Time = time_range,
          Proportion = predictions$Prediction,
          Proportion_lower = predictions$Prediction - 1.96 * predictions$SE,
          Proportion_upper = predictions$Prediction + 1.96 * predictions$SE)

#Fast, early
all_fast_early<-all_outdoor%>%filter(Speed=="Fast"&sowing.moment=="early")

timeBef <- c(0, all_fast_early$variable)
timeAf <- c(all_fast_early$variable, Inf)
Counts <- c(all_fast_early$value, unique(all_fast_early$Total_seeds) - sum(all_fast_early$value))
  
df <- data.frame(timeBef, timeAf, Counts)
  
modTE_fast_early <-drm(Counts ~ timeBef + timeAf, data = df, 
        fct = LN.3(), type = "event")

time_range <- seq(0,35,1)
predictions <- as.data.frame(predict(modTE_fast_early, 
                       newdata = data.frame(timeBef = time_range), 
                       type = "response", se.fit = TRUE))

Fast_early_curve<-data.frame(Speed = "Fast",
          Speed_T50 = summary(modTE_fast_early)[[3]][3, 1],
          sowing.moment="early",
          ES = predictions$SE,
          Time = time_range,
          Proportion = predictions$Prediction,
          Proportion_lower = predictions$Prediction - 1.96 * predictions$SE,
          Proportion_upper = predictions$Prediction + 1.96 * predictions$SE)

#Fast, late
all_fast_late<-all_outdoor%>%filter(Speed=="Fast"&sowing.moment=="late")

timeBef <- c(0, all_fast_late$variable)
timeAf <- c(all_fast_late$variable, Inf)
Counts <- c(all_fast_late$value, unique(all_fast_late$Total_seeds) - sum(all_fast_late$value))
  
df <- data.frame(timeBef, timeAf, Counts)
  
modTE_fast_late <-drm(Counts ~ timeBef + timeAf, data = df, 
        fct = LN.3(), type = "event")

time_range <- seq(0,35,1)

predictions <- as.data.frame(predict(modTE_fast_late, 
                       newdata = data.frame(timeBef = time_range), 
                       type = "response", se.fit = TRUE))

Fast_late_curve<-data.frame(Speed = "Fast",
          Speed_T50 = summary(modTE_fast_late)[[3]][3, 1],
          sowing.moment="late",
          ES = predictions$SE,
          Time = time_range,
          Proportion = predictions$Prediction,
          Proportion_lower = predictions$Prediction - 1.96 * predictions$SE,
          Proportion_upper = predictions$Prediction + 1.96 * predictions$SE)

####

combined_curves<-Fast_early_curve%>%
  bind_rows(Fast_late_curve)%>%
  bind_rows(Slow_early_curve)%>%
  bind_rows(Slow_late_curve)

##
## Create plots ##
extra_rows <- germ_outdoor %>%
  distinct(Species, sowing.moment, Total_seeds, Code) %>%
  mutate(variable = 0, value = 0, germinated = 0)

# Bind the new rows to the original dataframe
germ_outdoor_zeros <- bind_rows(germ_outdoor, extra_rows) %>%
  arrange(Species, sowing.moment, Total_seeds, variable)  # Arrange the data nicely

raw_data_outdoor<-ggplot(germ_outdoor_zeros,aes(variable,germinated/Total_seeds,col=sowing.moment))+
  geom_line()+
  geom_point()+
  facet_wrap(~Species)+
  xlab("Days")+
  ylab("Proportion of seed germinated")+ 
  theme( # Remove x-axis ticks
        panel.grid.major.x = element_blank(), # Remove major gridlines
        panel.grid.minor.x = element_blank(), # Remove minor gridlines
        panel.border = element_rect(color = "black", fill = NA, size = 1),
        strip.background = element_blank(),
        panel.background = element_blank(),
        legend.box.background = element_blank(),
        legend.position = "bottom")+
  labs(color = "Sowing moment")

median<-T50_outdoor%>%
  dplyr::select(Species,Speed_T50)%>%
  distinct()%>%
  na.omit()%>%
  mutate(median=median(Speed_T50))
  
time_to_event_invidivual_species_outdoor <- ggplot(T50_outdoor) +
  geom_line(aes(x = Time, y = Proportion, col = sowing.moment)) +
  facet_wrap(~Speed + Species) +
  labs(x = "Time (days)", 
       y = "Proportion germinated") +
  geom_vline(xintercept = 13.48184, linetype = "dashed") +
  geom_point(aes(x = Speed_T50, y = Proportion_T50, col = sowing.moment)) +
  xlim(0, 60)+
  ylim(0,1)+
  theme( # Remove x-axis ticks
        panel.grid.major.x = element_blank(), # Remove major gridlines
        panel.grid.minor.x = element_blank(), # Remove minor gridlines
        panel.border = element_rect(color = "black", fill = NA, size = 1),
        strip.background = element_blank(),
        panel.background = element_blank(),
        legend.box.background = element_blank(),
        legend.position = "bottom")+
  labs(color = "Sowing moment")

time_to_event_groups_outdoor <- ggplot() +
  labs(title = "Field", 
       x = "Time (days)", 
       y = "Proportion germinated") +
  xlim(0, 35)+ 
  theme( # Remove x-axis ticks
        panel.grid.major.x = element_blank(), # Remove major gridlines
        panel.grid.minor.x = element_blank(), # Remove minor gridlines
        panel.border = element_rect(color = "black", fill = NA, size = 1),
        strip.background = element_blank(),
        panel.background = element_blank(),
        legend.box.background = element_blank())+
  labs(color = "Speed\ngroup",
       shape = "Sowing\nmoment",
       linetype = "Sowing\nmoment")+  
  scale_color_manual(values = c("Fast" = "#E69F00",
                               "Slow" = "#009E73"))+
  geom_line(data=combined_curves,aes(Time,Proportion,col=Speed,linetype = sowing.moment),size=1)+
  geom_point(data=germ_outdoor%>%inner_join(speciesinfo),aes(variable,germinated/Total_seeds,col=Speed, shape=sowing.moment),position = position_dodge(2),alpha=.3)



raw_data_outdoor
time_to_event_invidivual_species_outdoor 

combined_groups<-time_to_event_groups_indoor+time_to_event_groups_outdoor+
  plot_layout(ncol = 1)+
                 plot_annotation(tag_levels = 'a')
combined_groups

#ggsave(file.path("C:/Users/tava0008/OneDrive - Umeå universitet/Documents/PhD/Germination project/latex/outdoor_germination_raw.png"), plot = raw_data_outdoor, width =9, height =6, units = "in", dpi = 300)

#ggsave(file.path("C:/Users/tava0008/OneDrive - Umeå universitet/Documents/PhD/Germination project/latex/outdoor_germination_individual.png"), plot = time_to_event_invidivual_species_outdoor, width = 7.5, height =6, units = "in", dpi = 300)

#ggsave(file.path("C:/Users/tava0008/OneDrive - Umeå universitet/Documents/PhD/Germination project/latex/combined_groups.png"), plot = combined_groups, width = 5, height =6, units = "in", dpi = 300)

#svg("C:/Users/tava0008/OneDrive - Umeå universitet/Documents/PhD/Germination #project/latex/combined_groups.svg", width = 4, height = 5)  # specify the file name and dimensions
#print(combined_groups)
#dev.off()
```

```{r LMER outdoor germination }
LMER_data_outdoor_germination<-germ_outdoor%>%
  group_by(Species,sowing.moment)%>%
  mutate(max_germination=max(germinated)/Total_seeds)%>%
  ungroup()%>%
  distinct(Species,sowing.moment,max_germination)%>%
  inner_join((T50_outdoor%>%
  distinct(Species,sowing.moment,Speed,Speed_T50)))%>%
  as.data.frame()

max_germination_outdoor_plot<-LMER_data_outdoor_germination%>%
  ggplot(aes(Speed,max_germination,col=sowing.moment))+
  geom_boxplot()+ 
  theme( # Remove x-axis ticks
        panel.grid.major.x = element_blank(), # Remove major gridlines
        panel.grid.minor.x = element_blank(), # Remove minor gridlines
        panel.border = element_rect(color = "black", fill = NA, size = 1),
        strip.background = element_blank(),
        panel.background = element_blank(),
        legend.box.background = element_blank(),
        legend.position = "bottom") +
  labs(x = "Speed group", 
       y = "Germination rate",
       color = "Sowing moment") 

max_germination_outdoor_plot
#ggsave(file.path("C:/Users/tava0008/OneDrive - Umeå universitet/Documents/PhD/Germination project/latex/max_germination_outdoor_plot.png"), plot = max_germination_outdoor_plot, width = 3, height =2.5, units = "in", dpi = 300)

T50_outdoor_plot<-LMER_data_outdoor_germination%>%
  ggplot(aes(Speed,Speed_T50,col=sowing.moment))+
  geom_boxplot()+ 
  theme( # Remove x-axis ticks
        panel.grid.major.x = element_blank(), # Remove major gridlines
        panel.grid.minor.x = element_blank(), # Remove minor gridlines
        panel.border = element_rect(color = "black", fill = NA, size = 1),
        strip.background = element_blank(),
        panel.background = element_blank(),
        legend.box.background = element_blank(),
        legend.position = "bottom") +
  labs(title = "Germination speed (T50)", 
       x = "Speed group", 
       y = "T50",
       color = "Sowing moment") 
T50_outdoor_plot

#ggsave(file.path("C:/Users/tava0008/OneDrive - Umeå universitet/Documents/PhD/Germination project/latex/T50_outdoor_plot.png"), plot = T50_outdoor_plot, width = 3, height =3, units = "in", dpi = 300)

anova(lmer(max_germination~Speed*sowing.moment+(1|Species),data = LMER_data_outdoor_germination))

anova(lmer(Speed_T50~Speed*sowing.moment+(1|Species),data = LMER_data_outdoor_germination))

LMER_data_outdoor_germination %>%
  group_by(Speed)%>%
  summarize(
    mean_max_germination = mean(max_germination, na.rm = TRUE),
    sd_max_germination = sd(max_germination, na.rm = TRUE),
    mean_Speed_T50 = mean(Speed_T50, na.rm = TRUE),
    sd_Speed_T50 = sd(Speed_T50, na.rm = TRUE)
  )

LMER_data_outdoor_germination %>%
  group_by(sowing.moment)%>%
  summarize(
    mean_max_germination = mean(max_germination, na.rm = TRUE),
    sd_max_germination = sd(max_germination, na.rm = TRUE),
    mean_Speed_T50 = mean(Speed_T50, na.rm = TRUE),
    sd_Speed_T50 = sd(Speed_T50, na.rm = TRUE)
  )
```

```{r indoor outdoor correlations}
germ_rate<-germ_indoor%>%
  inner_join(speciesinfo)%>%
  group_by(Species,Temp)%>%
  filter(Temp!=20)%>%
  mutate(germinated=max(germinated)/Total_seeds)%>%
  ungroup()%>%
  distinct(Species,Temp,germinated)%>%
  rename(condition=Temp)%>%
  mutate(condition=as.factor(condition))%>%
  bind_rows(

germ_outdoor%>%
  inner_join(speciesinfo)%>%
  group_by(Species,sowing.moment)%>%
  mutate(germinated=max(germinated)/Total_seeds)%>%
  ungroup()%>%
  distinct(Species,sowing.moment,germinated)%>%
  rename(condition=sowing.moment)
)

germ_rate_correlations<-germ_rate%>%
  pivot_wider(names_from = condition, values_from = germinated)%>%
  column_to_rownames(var = "Species")%>%
  na.omit()
  
T50_correlations<-T50_indoor%>%
  distinct(Species,Temp,Speed_T50)%>%
  filter(Temp!=20)%>%
  rename(condition=Temp)%>%
  mutate(condition=as.factor(condition))%>%
  bind_rows(
T50_outdoor%>%
  rename(condition=sowing.moment)%>%
  distinct(Species,condition,Speed_T50))

T50_correlations<-T50_correlations%>%
  pivot_wider(names_from = condition, values_from = Speed_T50)%>%
  column_to_rownames(var = "Species")%>%
  na.omit()
  
# Compute the correlation matrix
cor_matrix <- cor(germ_rate_correlations, use = "complete.obs")

# Create the correlation plot
rate_corr<-ggcorrplot(cor_matrix, 
           method = "square",     # Use "number", "circle", "square"
           type = "lower",        # "full", "lower", or "upper" triangle
           lab = TRUE,            # Show correlation values
           lab_size = 3,          # Size of correlation text
           colors = c("blue", "white", "red"),  # Color gradient
           title = "Germination rate",
           outline.col = "black")+
  theme(legend.position = "none")

cor_matrix <- cor(T50_correlations, use = "complete.obs")

# Create the correlation plot
rate_speed<-ggcorrplot(cor_matrix, 
           method = "square",     # Use "number", "circle", "square"
           type = "lower",        # "full", "lower", or "upper" triangle
           lab = TRUE,            # Show correlation values
           lab_size = 3,          # Size of correlation text
           colors = c("blue", "white", "red"),  # Color gradient
           title = "Germination speed",
           outline.col = "black")

corr_plot<-rate_speed|rate_corr
corr_plot
#ggsave(file.path("C:/Users/tava0008/OneDrive - Umeå universitet/Documents/PhD/Germination project/latex/germ_corr.png"), plot = corr_plot, width = 6, height =3, units = "in", dpi = 300)

```