library("drc")
library("readxl")
library("caroline")
library("ggplot2")
library("dplyr")
library("multcompView")
library("ggthemes")
library("car")
library("tidyverse")
library("ggpubr")
library("writexl")


results$RGR[results$RGR < 0] <- 0
results$RGR_OECD[results$RGR < 0] <- 0

results <- results %>% 
  rename(
    conc_Me1 = conc_As,
    conc_Me2 = conc_Zn,
    conc_Me3 = conc_Pb,
    conc_Me4 = conc_Ag,
    conc_Me5 = conc_Cu
  )

Me1 <- subset(results,metal=="As")
Me2 <- subset(results,metal=="Zn")
Me3 <- subset(results,metal=="Pb")
Me4 <- subset(results,metal=="Ag")
Me5 <- subset(results,metal=="Cu")
B  <- subset(results,metal=="B")
T  <- subset(results,metal=="T")
Q  <- subset(results,metal=="Q")
R  <- subset(results,metal=="R")
X  <- subset(results,metal=="X")

#to remove outliers (#below the LOQ)
Q <- Q[Q$No !=1 , ]  
X <- X[X$No !=1 , ]  
R <- R[R$No !=1 , ] 
R <- R[R$No !=2 , ]
R <- R[R$No !=3 , ]
#N.B. don't forget to remove the outliers also from the CA calculation!!

#remove the No columns
Me1 <- Me1[ , -which(names(Me1) %in% c("No"))]
Me2 <- Me2[ , -which(names(Me2) %in% c("No"))]
Me3 <- Me3[ , -which(names(Me3) %in% c("No"))]
Me4 <- Me4[ , -which(names(Me4) %in% c("No"))]
Me5 <- Me5[ , -which(names(Me5) %in% c("No"))]

B <- B[ , -which(names(B) %in% c("No"))]
T <- T[ , -which(names(T) %in% c("No"))]
Q <- Q[ , -which(names(Q) %in% c("No"))]
R <- R[ , -which(names(R) %in% c("No"))]
X <- X[ , -which(names(X) %in% c("No"))]


#############################################################
#SUMMARY TABLE 
#############################################################

Me1_DRM      <- drm(Me1$RGR~Me1$conc_Me1, data = Me1, fct = LL.2())
Me1_DRM_OECD <- drm(Me1$RGR_OECD~Me1$conc_Me1, data = Me1, fct = LL.2())

Me2_DRM      <- drm(Me2$RGR~Me2$conc_Me2, data = Me2, fct = LL.2())
Me2_DRM_OECD <- drm(Me2$RGR_OECD~Me2$conc_Me2, data = Me2, fct = LL.2())

Me3_DRM      <- drm(Me3$RGR~Me3$conc_Me3, data = Me3, fct = LL.2())
Me3_DRM_OECD <- drm(Me3$RGR_OECD~Me3$conc_Me3, data = Me3, fct = LL.2())

Me4_DRM      <- drm(Me4$RGR~Me4$conc_Me4, data = Me4, fct = LL.2())
Me4_DRM_OECD <- drm(Me4$RGR_OECD~Me4$conc_Me4, data = Me4, fct = LL.2())

Me5_DRM      <- drm(Me5$RGR~Me5$conc_Me5, data = Me5, fct = LL.2())
Me5_DRM_OECD <- drm(Me5$RGR_OECD~Me5$conc_Me5, data = Me5, fct = LL.2())

#EC AND SLOPE TABLE --> SE
ECMe1 <- ED(Me1_DRM,c(10,20,50), interval="delta")
ECMe2 <- ED(Me2_DRM,c(10,20,50), interval="delta")
ECMe3 <- ED(Me3_DRM,c(10,20,50), interval="delta")
ECMe4 <- ED(Me4_DRM,c(10,20,50), interval="delta")
ECMe5 <- ED(Me5_DRM,c(10,20,50), interval="delta")


EC10Me1 <- ECMe1[1:1]
EC20Me1 <- ECMe1["e:1:20" , "Estimate"]
EC50Me1 <- ECMe1["e:1:50" , "Estimate"]
SE10Me1 <- ECMe1["e:1:10" , "Std. Error"]
SE20Me1 <- ECMe1["e:1:20" , "Std. Error"]
SE50Me1 <- ECMe1["e:1:50" , "Std. Error"]
U10Me1 <-  ECMe1["e:1:10" , "Upper"]
L10Me1 <-  ECMe1["e:1:10" , "Lower"]
U20Me1 <-  ECMe1["e:1:20" , "Upper"]
L20Me1 <-  ECMe1["e:1:20" , "Lower"]
U50Me1 <-  ECMe1["e:1:50" , "Upper"]
L50Me1 <-  ECMe1["e:1:50" , "Lower"]


EC10Me2 <- ECMe2[1:1]
EC20Me2 <- ECMe2["e:1:20" , "Estimate"]
EC50Me2 <- ECMe2["e:1:50" , "Estimate"]
SE10Me2 <- ECMe2["e:1:10" , "Std. Error"]
SE20Me2 <- ECMe2["e:1:20" , "Std. Error"]
SE50Me2 <- ECMe2["e:1:50" , "Std. Error"]
U10Me2 <-  ECMe2["e:1:10" , "Upper"]
L10Me2 <-  ECMe2["e:1:10" , "Lower"]
U20Me2 <-  ECMe2["e:1:20" , "Upper"]
L20Me2 <-  ECMe2["e:1:20" , "Lower"]
U50Me2 <-  ECMe2["e:1:50" , "Upper"]
L50Me2 <-  ECMe2["e:1:50" , "Lower"]


EC10Me3 <- ECMe3[1:1]
EC20Me3 <- ECMe3["e:1:20" , "Estimate"]
EC50Me3 <- ECMe3["e:1:50" , "Estimate"]
SE10Me3 <- ECMe3["e:1:10" , "Std. Error"]
SE20Me3 <- ECMe3["e:1:20" , "Std. Error"]
SE50Me3 <- ECMe3["e:1:50" , "Std. Error"]
U10Me3 <-  ECMe3["e:1:10" , "Upper"]
L10Me3 <-  ECMe3["e:1:10" , "Lower"]
U20Me3 <-  ECMe3["e:1:20" , "Upper"]
L20Me3 <-  ECMe3["e:1:20" , "Lower"]
U50Me3 <-  ECMe3["e:1:50" , "Upper"]
L50Me3 <-  ECMe3["e:1:50" , "Lower"]


EC10Me4 <- ECMe4[1:1]
EC20Me4 <- ECMe4["e:1:20" , "Estimate"]
EC50Me4 <- ECMe4["e:1:50" , "Estimate"]
SE10Me4 <- ECMe4["e:1:10" , "Std. Error"]
SE20Me4 <- ECMe4["e:1:20" , "Std. Error"]
SE50Me4 <- ECMe4["e:1:50" , "Std. Error"]
U10Me4 <-  ECMe4["e:1:10" , "Upper"]
L10Me4 <-  ECMe4["e:1:10" , "Lower"]
U20Me4 <-  ECMe4["e:1:20" , "Upper"]
L20Me4 <-  ECMe4["e:1:20" , "Lower"]
U50Me4 <-  ECMe4["e:1:50" , "Upper"]
L50Me4 <-  ECMe4["e:1:50" , "Lower"]


EC10Me5 <- ECMe5[1:1]
EC20Me5 <- ECMe5["e:1:20" , "Estimate"]
EC50Me5 <- ECMe5["e:1:50" , "Estimate"]
SE10Me5 <- ECMe5["e:1:10" , "Std. Error"]
SE20Me5 <- ECMe5["e:1:20" , "Std. Error"]
SE50Me5 <- ECMe5["e:1:50" , "Std. Error"]
U10Me5 <-  ECMe5["e:1:10" , "Upper"]
L10Me5 <-  ECMe5["e:1:10" , "Lower"]
U20Me5 <-  ECMe5["e:1:20" , "Upper"]
L20Me5 <-  ECMe5["e:1:20" , "Lower"]
U50Me5 <-  ECMe5["e:1:50" , "Upper"]
L50Me5 <-  ECMe5["e:1:50" , "Lower"]


Me1_DRM_df <-  summary(Me1_DRM)
Me1_DRM_df <- Me1_DRM_df$coefficients
Me1_DRM_df <- as.matrix(Me1_DRM_df, header=TRUE)
BMe1   <- Me1_DRM_df[1:1]
SEBMe1 <- Me1_DRM_df["b:(Intercept)", "Std. Error"]


Me2_DRM_df <-  summary(Me2_DRM)
Me2_DRM_df <- Me2_DRM_df$coefficients
Me2_DRM_df <- as.matrix(Me2_DRM_df, header=TRUE)
BMe2   <- Me2_DRM_df[1:1]
SEBMe2 <- Me2_DRM_df["b:(Intercept)", "Std. Error"]


Me3_DRM_df <-  summary(Me3_DRM)
Me3_DRM_df <- Me3_DRM_df$coefficients
Me3_DRM_df <- as.matrix(Me3_DRM_df, header=TRUE)
BMe3   <- Me3_DRM_df[1:1]
SEBMe3 <- Me3_DRM_df["b:(Intercept)", "Std. Error"]

Me4_DRM_df <-  summary(Me4_DRM)
Me4_DRM_df <- Me4_DRM_df$coefficients
Me4_DRM_df <- as.matrix(Me4_DRM_df, header=TRUE)
BMe4   <- Me4_DRM_df[1:1]
SEBMe4 <- Me4_DRM_df["b:(Intercept)", "Std. Error"]

Me5_DRM_df <-  summary(Me5_DRM)
Me5_DRM_df <- Me5_DRM_df$coefficients
Me5_DRM_df <- as.matrix(Me5_DRM_df, header=TRUE)
BMe5   <- Me5_DRM_df[1:1]
SEBMe5 <- Me5_DRM_df["b:(Intercept)", "Std. Error"]

Results_summary <- data.frame(Me1 = numeric(14), Me2 = numeric(14), Me3 = numeric(14), 
                              Me4 = numeric(14), Me5 =numeric(14),
                              row.names = c("EC10", "EC20", "EC50", "slope", "SE_slope", "SE_EC10", "SE_EC20","SE_EC50", 
                                            "Upper_EC10", "Lower_EC10", "Upper_EC20", "Lower_EC20", "Upper_EC50", "Lower_EC50"))

Results_summary[1, "Me1"] <- as.numeric(EC10Me1)
Results_summary[1, "Me2"] <- as.numeric(EC10Me2)
Results_summary[1, "Me3"] <- as.numeric(EC10Me3)
Results_summary[1, "Me4"] <- as.numeric(EC10Me4)
Results_summary[1, "Me5"] <- as.numeric(EC10Me5)

Results_summary[2, "Me1"] <- as.numeric(EC20Me1)
Results_summary[2, "Me2"] <- as.numeric(EC20Me2)
Results_summary[2, "Me3"] <- as.numeric(EC20Me3)
Results_summary[2, "Me4"] <- as.numeric(EC20Me4)
Results_summary[2, "Me5"] <- as.numeric(EC20Me5)

Results_summary[3, "Me1"] <- as.numeric(EC50Me1)
Results_summary[3, "Me2"] <- as.numeric(EC50Me2)
Results_summary[3, "Me3"] <- as.numeric(EC50Me3)
Results_summary[3, "Me4"] <- as.numeric(EC50Me4)
Results_summary[3, "Me5"] <- as.numeric(EC50Me5)

Results_summary[4, "Me1"] <- as.numeric(BMe1)
Results_summary[4, "Me2"] <- as.numeric(BMe2)
Results_summary[4, "Me3"] <- as.numeric(BMe3)
Results_summary[4, "Me4"] <- as.numeric(BMe4)
Results_summary[4, "Me5"] <- as.numeric(BMe5)

Results_summary[5, "Me1"] <- as.numeric(SEBMe1)
Results_summary[5, "Me2"] <- as.numeric(SEBMe2)
Results_summary[5, "Me3"] <- as.numeric(SEBMe3)
Results_summary[5, "Me4"] <- as.numeric(SEBMe4)
Results_summary[5, "Me5"] <- as.numeric(SEBMe5)

Results_summary[6, "Me1"] <- as.numeric(SE10Me1)
Results_summary[6, "Me2"] <- as.numeric(SE10Me2)
Results_summary[6, "Me3"] <- as.numeric(SE10Me3)
Results_summary[6, "Me4"] <- as.numeric(SE10Me4)
Results_summary[6, "Me5"] <- as.numeric(SE10Me5)

Results_summary[7, "Me1"] <- as.numeric(SE20Me1)
Results_summary[7, "Me2"] <- as.numeric(SE20Me2)
Results_summary[7, "Me3"] <- as.numeric(SE20Me3)
Results_summary[7, "Me4"] <- as.numeric(SE20Me4)
Results_summary[7, "Me5"] <- as.numeric(SE20Me5)

Results_summary[8, "Me1"] <- as.numeric(SE50Me1)
Results_summary[8, "Me2"] <- as.numeric(SE50Me2)
Results_summary[8, "Me3"] <- as.numeric(SE50Me3)
Results_summary[8, "Me4"] <- as.numeric(SE50Me4)
Results_summary[8, "Me5"] <- as.numeric(SE50Me5)

Results_summary[9, "Me1"] <- as.numeric(U10Me1)
Results_summary[9, "Me2"] <- as.numeric(U10Me2)
Results_summary[9, "Me3"] <- as.numeric(U10Me3)
Results_summary[9, "Me4"] <- as.numeric(U10Me4)
Results_summary[9, "Me5"] <- as.numeric(U10Me5)

Results_summary[10, "Me1"] <- as.numeric(L10Me1)
Results_summary[10, "Me2"] <- as.numeric(L10Me2)
Results_summary[10, "Me3"] <- as.numeric(L10Me3)
Results_summary[10, "Me4"] <- as.numeric(L10Me4)
Results_summary[10, "Me5"] <- as.numeric(L10Me5)

Results_summary[11, "Me1"] <- as.numeric(U20Me1)
Results_summary[11, "Me2"] <- as.numeric(U20Me2)
Results_summary[11, "Me3"] <- as.numeric(U20Me3)
Results_summary[11, "Me4"] <- as.numeric(U20Me4)
Results_summary[11, "Me5"] <- as.numeric(U20Me5)

Results_summary[12, "Me1"] <- as.numeric(L20Me1)
Results_summary[12, "Me2"] <- as.numeric(L20Me2)
Results_summary[12, "Me3"] <- as.numeric(L20Me3)
Results_summary[12, "Me4"] <- as.numeric(L20Me4)
Results_summary[12, "Me5"] <- as.numeric(L20Me5)

Results_summary[13, "Me1"] <- as.numeric(U50Me1)
Results_summary[13, "Me2"] <- as.numeric(U50Me2)
Results_summary[13, "Me3"] <- as.numeric(U50Me3)
Results_summary[13, "Me4"] <- as.numeric(U50Me4)
Results_summary[13, "Me5"] <- as.numeric(U50Me5)

Results_summary[14, "Me1"] <- as.numeric(L50Me1)
Results_summary[14, "Me2"] <- as.numeric(L50Me2)
Results_summary[14, "Me3"] <- as.numeric(L50Me3)
Results_summary[14, "Me4"] <- as.numeric(L50Me4)
Results_summary[14, "Me5"] <- as.numeric(L50Me5)

write.table(Results_summary, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Mixture_interaction_analysis//input_df//Results_summary_SE.txt", 
            sep = "\t", quote = FALSE, row.names = TRUE)

Results_summary <- Results_summary %>% 
  rename(
    As = Me1,
    Zn = Me2,
    Pb = Me3,
    Ag = Me4,
    Cu = Me5
  )

write.table(Results_summary, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Summaries//Results_summary_SE.txt", 
            sep = "\t", quote = FALSE, row.names = TRUE)

#writexl::write_xlsx(Results_summary, "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Summaries//Summary_results.xlsx")


#################################################################################

#TU calculation
Me1$TU10 <- as.numeric((Me1$conc_Me1) / EC10Me1)
Me1$TU20 <- as.numeric((Me1$conc_Me1) / EC20Me1)
Me1$TU50 <- as.numeric((Me1$conc_Me1) / EC50Me1)

Me2$TU10 <- as.numeric((Me2$conc_Me2) / EC10Me2)
Me2$TU20 <- as.numeric((Me2$conc_Me2) / EC20Me2)
Me2$TU50 <- as.numeric((Me2$conc_Me2) / EC50Me2)

Me3$TU10 <- as.numeric((Me3$conc_Me3) / EC10Me3)
Me3$TU20 <- as.numeric((Me3$conc_Me3) / EC20Me3)
Me3$TU50 <- as.numeric((Me3$conc_Me3) / EC50Me3)

Me4$TU10 <- as.numeric((Me4$conc_Me4) / EC10Me4)
Me4$TU20 <- as.numeric((Me4$conc_Me4) / EC20Me4)
Me4$TU50 <- as.numeric((Me4$conc_Me4) / EC50Me4)

Me5$TU10 <- as.numeric((Me5$conc_Me5) / EC10Me5)
Me5$TU20 <- as.numeric((Me5$conc_Me5) / EC20Me5)
Me5$TU50 <- as.numeric((Me5$conc_Me5) / EC50Me5)


Me1$IA_Me1 <- 1/(1+(Me1$TU50)^BMe1)   
Me2$IA_Me2 <- 1/(1+(Me2$TU50)^BMe2)
Me3$IA_Me3 <- 1/(1+(Me3$TU50)^BMe3)
Me4$IA_Me4 <- 1/(1+(Me4$TU50)^BMe4)
Me5$IA_Me5 <- 1/(1+(Me5$TU50)^BMe5)

#-------------------------------------------
#BINARY - TU SUM

#10
B$TU10_Me1 <- as.numeric((B$conc_Me1) / EC10Me1)
B$TU10_Me2 <- as.numeric((B$conc_Me2) / EC10Me2)

B$B_sum_TU10 <- as.numeric(B$TU10_Me1 + B$TU10_Me2)

#20
B$TU20_Me1 <- as.numeric((B$conc_Me1) / EC20Me1)
B$TU20_Me2 <- as.numeric((B$conc_Me2) / EC20Me2)

B$B_sum_TU20 <- as.numeric(B$TU20_Me1 + B$TU20_Me2)

#50
B$TU_Me1 <- as.numeric((B$conc_Me1) / EC50Me1)
B$TU_Me2 <- as.numeric((B$conc_Me2) / EC50Me2)

B$B_sum_TU <- as.numeric(B$TU_Me1 + B$TU_Me2)

B$IA <- (1/(1+(B$conc_Me1/EC50Me1)^BMe1))*(1/(1+(B$conc_Me2/EC50Me2)^BMe2))

#---------------------------------------------
#TERNARY - TU SUM

#10
T$TU10_Me1 <- as.numeric((T$conc_Me1) / EC10Me1)
T$TU10_Me2 <- as.numeric((T$conc_Me2) / EC10Me2)
T$TU10_Me3 <- as.numeric((T$conc_Me3) / EC10Me3)

T$T_sum_TU10 <- as.numeric(T$TU10_Me1 + T$TU10_Me2 + T$TU10_Me3)

#20
T$TU20_Me1 <- as.numeric((T$conc_Me1) / EC20Me1)
T$TU20_Me2 <- as.numeric((T$conc_Me2) / EC20Me2)
T$TU20_Me3 <- as.numeric((T$conc_Me3) / EC20Me3)

T$T_sum_TU20 <- as.numeric(T$TU20_Me1 + T$TU20_Me2 + T$TU20_Me3)

#50
T$TU_Me1 <- as.numeric((T$conc_Me1) / EC50Me1)
T$TU_Me2 <- as.numeric((T$conc_Me2) / EC50Me2)
T$TU_Me3 <- as.numeric((T$conc_Me3) / EC50Me3)

T$T_sum_TU <- as.numeric(T$TU_Me1 + T$TU_Me2 + T$TU_Me3)

T$IA <- (1/(1+(T$conc_Me1/EC50Me1)^BMe1))*(1/(1+(T$conc_Me2/EC50Me2)^BMe2))*(1/(1+(T$conc_Me3/EC50Me3)^BMe3))

#-----------------------------------------------------
#QUATERNARY - TU SUM

#10
Q$TU10_Me1 <- as.numeric((Q$conc_Me1) / EC10Me1)
Q$TU10_Me2 <- as.numeric((Q$conc_Me2) / EC10Me2)
Q$TU10_Me3 <- as.numeric((Q$conc_Me3) / EC10Me3)
Q$TU10_Me4 <- as.numeric((Q$conc_Me4) / EC10Me4)

Q$Q_sum_TU10 <- as.numeric(Q$TU10_Me1 + Q$TU10_Me2 + Q$TU10_Me3 + Q$TU10_Me4)

#20
Q$TU20_Me1 <- as.numeric((Q$conc_Me1) / EC20Me1)
Q$TU20_Me2 <- as.numeric((Q$conc_Me2) / EC20Me2)
Q$TU20_Me3 <- as.numeric((Q$conc_Me3) / EC20Me3)
Q$TU20_Me4 <- as.numeric((Q$conc_Me4) / EC20Me4)

Q$Q_sum_TU20 <- as.numeric(Q$TU20_Me1 + Q$TU20_Me2 + Q$TU20_Me3 + Q$TU20_Me4)

#50
Q$TU_Me1 <- as.numeric((Q$conc_Me1) / EC50Me1)
Q$TU_Me2 <- as.numeric((Q$conc_Me2) / EC50Me2)
Q$TU_Me3 <- as.numeric((Q$conc_Me3) / EC50Me3)
Q$TU_Me4 <- as.numeric((Q$conc_Me4) / EC50Me4)

Q$Q_sum_TU <- as.numeric(Q$TU_Me1 + Q$TU_Me2 + Q$TU_Me3 + Q$TU_Me4)

Q$IA <- (1/(1+(Q$conc_Me1/EC50Me1)^BMe1))*(1/(1+(Q$conc_Me2/EC50Me2)^BMe2))*(1/(1+(Q$conc_Me3/EC50Me3)^BMe3))*(1/(1+(Q$conc_Me4/EC50Me4)^BMe4))

#--------------------------------------------------------
#QUINARY - environmentally relevant ray - TU SUM

#10
R$TU10_Me1 <- as.numeric((R$conc_Me1) / EC10Me1)
R$TU10_Me2 <- as.numeric((R$conc_Me2) / EC10Me2)
R$TU10_Me3 <- as.numeric((R$conc_Me3) / EC10Me3)
R$TU10_Me4 <- as.numeric((R$conc_Me4) / EC10Me4)
R$TU10_Me5 <- as.numeric((R$conc_Me5) / EC10Me5)

R$R_sum_TU10 <- as.numeric(R$TU10_Me1 + R$TU10_Me2 + R$TU10_Me3 + R$TU10_Me4 + R$TU10_Me5)

#20
R$TU20_Me1 <- as.numeric((R$conc_Me1) / EC20Me1)
R$TU20_Me2 <- as.numeric((R$conc_Me2) / EC20Me2)
R$TU20_Me3 <- as.numeric((R$conc_Me3) / EC20Me3)
R$TU20_Me4 <- as.numeric((R$conc_Me4) / EC20Me4)
R$TU20_Me5 <- as.numeric((R$conc_Me5) / EC20Me5)

R$R_sum_TU20 <- as.numeric(R$TU20_Me1 + R$TU20_Me2 + R$TU20_Me3 + R$TU20_Me4 + R$TU20_Me5)

#50
R$TU_Me1 <- as.numeric((R$conc_Me1) / EC50Me1)
R$TU_Me2 <- as.numeric((R$conc_Me2) / EC50Me2)
R$TU_Me3 <- as.numeric((R$conc_Me3) / EC50Me3)
R$TU_Me4 <- as.numeric((R$conc_Me4) / EC50Me4)
R$TU_Me5 <- as.numeric((R$conc_Me5) / EC50Me5)

R$R_sum_TU <- as.numeric(R$TU_Me1 + R$TU_Me2 + R$TU_Me3 + R$TU_Me4 + R$TU_Me5)

R$IA <- (1/(1+(R$conc_Me1/EC50Me1)^BMe1))*(1/(1+(R$conc_Me2/EC50Me2)^BMe2))*(1/(1+(R$conc_Me3/EC50Me3)^BMe3))*(1/(1+(R$conc_Me4/EC50Me4)^BMe4))*(1/(1+(R$conc_Me5/EC50Me5)^BMe5))

#--------------------------------------------------------------
#QUINARY - equitoxic ray - TU SUM

#10
X$TU10_Me1 <- as.numeric((X$conc_Me1) / EC10Me1)
X$TU10_Me2 <- as.numeric((X$conc_Me2) / EC10Me2)
X$TU10_Me3 <- as.numeric((X$conc_Me3) / EC10Me3)
X$TU10_Me4 <- as.numeric((X$conc_Me4) / EC10Me4)
X$TU10_Me5 <- as.numeric((X$conc_Me5) / EC10Me5)

X$X_sum_TU10 <- as.numeric(X$TU10_Me1 + X$TU10_Me2 + X$TU10_Me3 + X$TU10_Me4 + X$TU10_Me5)

#20
X$TU20_Me1 <- as.numeric((X$conc_Me1) / EC20Me1)
X$TU20_Me2 <- as.numeric((X$conc_Me2) / EC20Me2)
X$TU20_Me3 <- as.numeric((X$conc_Me3) / EC20Me3)
X$TU20_Me4 <- as.numeric((X$conc_Me4) / EC20Me4)
X$TU20_Me5 <- as.numeric((X$conc_Me5) / EC20Me5)

X$X_sum_TU20 <- as.numeric(X$TU20_Me1 + X$TU20_Me2 + X$TU20_Me3 + X$TU20_Me4 + X$TU20_Me5)

#50
X$TU_Me1 <- as.numeric((X$conc_Me1) / EC50Me1)
X$TU_Me2 <- as.numeric((X$conc_Me2) / EC50Me2)
X$TU_Me3 <- as.numeric((X$conc_Me3) / EC50Me3)
X$TU_Me4 <- as.numeric((X$conc_Me4) / EC50Me4)
X$TU_Me5 <- as.numeric((X$conc_Me5) / EC50Me5)

X$X_sum_TU <- as.numeric(X$TU_Me1 + X$TU_Me2 + X$TU_Me3 + X$TU_Me4 + X$TU_Me5)

X$IA <- (1/(1+(X$conc_Me1/EC50Me1)^BMe1))*(1/(1+(X$conc_Me2/EC50Me2)^BMe2))*(1/(1+(X$conc_Me3/EC50Me3)^BMe3))*(1/(1+(X$conc_Me4/EC50Me4)^BMe4))*(1/(1+(X$conc_Me5/EC50Me5)^BMe5))

#################################################################################################################################################################

smaller_Me1 <- Me1 %>%
  group_by(conc_Me1) %>%
  summarise(across(where(is.numeric), mean))

smaller_Me2 <- Me2 %>%
  group_by(conc_Me2) %>%
  summarise(across(where(is.numeric), mean))

smaller_Me3 <- Me3 %>%
  group_by(conc_Me3) %>%
  summarise(across(where(is.numeric), mean))

smaller_Me4 <- Me4 %>%
  group_by(conc_Me4) %>%
  summarise(across(where(is.numeric), mean))

smaller_Me5 <- Me5 %>%
  group_by(conc_Me5) %>%
  summarise(across(where(is.numeric), mean))

smaller_B <- B %>%
  group_by(conc_Me1) %>%
  summarise(across(where(is.numeric), mean))

smaller_T <- T %>%
  group_by(conc_Me1) %>%
  summarise(across(where(is.numeric), mean))

smaller_Q <- Q %>%
  group_by(conc_Me1) %>%
  summarise(across(where(is.numeric), mean))

smaller_R <- R %>%
  group_by(conc_Me1) %>%
  summarise(across(where(is.numeric), mean))

smaller_X <- X %>%
  group_by(conc_Me1) %>%
  summarise(across(where(is.numeric), mean))

#####################

smaller_Me1 <- smaller_Me1[-1 ,]
smaller_Me2 <- smaller_Me2[-1 ,]
smaller_Me3 <- smaller_Me3[-1 ,]
smaller_Me4 <- smaller_Me4[-1 ,]
smaller_Me5 <- smaller_Me5[-1 ,]
smaller_B <- smaller_B[-1 ,]
smaller_T <- smaller_T[-1 ,]
smaller_Q <- smaller_Q[-1 ,]
smaller_R <- smaller_R[-1 ,]
smaller_X <- smaller_X[-1 ,]

B_stat <- rbind(smaller_Me1[, c("conc_Me1", "conc_Me2", "RGR")], 
                smaller_Me2[, c("conc_Me1", "conc_Me2", "RGR")], 
                smaller_B[, c("conc_Me1", "conc_Me2", "RGR")])

T_stat <- rbind(smaller_Me1[, c("conc_Me1", "conc_Me2", "conc_Me3", "RGR")], 
                smaller_Me2[, c("conc_Me1", "conc_Me2", "conc_Me3", "RGR")], 
                smaller_Me3[, c("conc_Me1", "conc_Me2", "conc_Me3", "RGR")], 
                smaller_T[, c("conc_Me1", "conc_Me2", "conc_Me3", "RGR")])

Q_stat <- rbind(smaller_Me1[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "RGR")], 
                smaller_Me2[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "RGR")], 
                smaller_Me3[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "RGR")],
                smaller_Me4[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "RGR")], 
                smaller_Q[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "RGR")])

R_stat <- rbind(smaller_Me1[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")], 
                smaller_Me2[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")], 
                smaller_Me3[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")],
                smaller_Me4[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")],
                smaller_Me5[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")], 
                smaller_R[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")])

X_stat <- rbind(smaller_Me1[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")], 
                smaller_Me2[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")], 
                smaller_Me3[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")],
                smaller_Me4[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")],
                smaller_Me5[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")], 
                smaller_X[, c("conc_Me1", "conc_Me2", "conc_Me3", "conc_Me4", "conc_Me5", "RGR")])


B_stat[is.na(B_stat)] <- 0
T_stat[is.na(T_stat)] <- 0
Q_stat[is.na(Q_stat)] <- 0
R_stat[is.na(R_stat)] <- 0
X_stat[is.na(X_stat)] <- 0

write.table(B_stat, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Mixture_interaction_analysis//input_df//B_stat.txt", sep = "\t", quote = FALSE, row.names = TRUE)
write.table(T_stat, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Mixture_interaction_analysis//input_df//T_stat.txt", sep = "\t", quote = FALSE, row.names = TRUE)
write.table(Q_stat, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Mixture_interaction_analysis//input_df//Q_stat.txt", sep = "\t", quote = FALSE, row.names = TRUE)
write.table(R_stat, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Mixture_interaction_analysis//input_df//R_stat.txt", sep = "\t", quote = FALSE, row.names = TRUE)
write.table(X_stat, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Mixture_interaction_analysis//input_df//X_stat.txt", sep = "\t", quote = FALSE, row.names = TRUE)


#########################################################################################################

#remove all NA columns
Me1 <- Me1[ , colSums(is.na(Me1))==0]
Me2 <- Me2[ , colSums(is.na(Me2))==0]
Me3 <- Me3[ , colSums(is.na(Me3))==0]
Me4 <- Me4[ , colSums(is.na(Me4))==0]
Me5 <- Me5[ , colSums(is.na(Me5))==0]

B <- B[ , colSums(is.na(B))==0]
T <- T[ , colSums(is.na(T))==0]
Q <- Q[ , colSums(is.na(Q))==0]
R <- R[ , colSums(is.na(R))==0]
X <- X[ , colSums(is.na(X))==0]


#the first TU to not be 0 (to avoid problems with log graph)
Me1$TU10[1] <- 0.00001
Me1$TU20[1] <- 0.00001
Me1$TU50[1] <- 0.00001

Me2$TU10[1] <- 0.00001
Me2$TU20[1] <- 0.00001
Me2$TU50[1] <- 0.00001

Me3$TU10[1] <- 0.00001
Me3$TU20[1] <- 0.00001
Me3$TU50[1] <- 0.00001

Me4$TU10[1] <- 0.00001
Me4$TU20[1] <- 0.00001
Me4$TU50[1] <- 0.00001

Me5$TU10[1] <- 0.00001
Me5$TU20[1] <- 0.00001
Me5$TU50[1] <- 0.00001

B <- within(B, {TU_Me1[1] <- 0.00001; 
                TU_Me2[1] <- 0.00001; 
                B_sum_TU[1] <- 0.00001})

T <- within(T, {TU_Me1[1] <- 0.00001; 
                TU_Me2[1] <- 0.00001; 
                TU_Me3[1] <- 0.00001;
                T_sum_TU[1] <- 0.00001})

Q <- within(Q, {TU_Me1[1] <- 0.00001; 
                TU_Me2[1] <- 0.00001; 
                TU_Me3[1] <- 0.00001;
                TU_Me4[1] <- 0.00001;
                Q_sum_TU[1] <- 0.00001})

R <- within(R, {TU_Me1[1] <- 0.00001; 
                TU_Me2[1] <- 0.00001; 
                TU_Me3[1] <- 0.00001;
                TU_Me4[1] <- 0.00001;
                TU_Me5[1] <- 0.00001;
                R_sum_TU[1] <- 0.00001})

X <- within(X, {TU_Me1[1] <- 0.00001; 
                TU_Me2[1] <- 0.00001; 
                TU_Me3[1] <- 0.00001;
                TU_Me4[1] <- 0.00001;
                TU_Me5[1] <- 0.00001;
                X_sum_TU[1] <- 0.00001})


# drm for the mixture curves, i.e., x= TU, TU EC50
Me1_DRM_TU <- drm(Me1$RGR~Me1$TU50, data = Me1, fct = LL.2())
Me2_DRM_TU <- drm(Me2$RGR~Me2$TU50, data = Me2, fct = LL.2())
Me3_DRM_TU <- drm(Me3$RGR~Me3$TU50, data = Me3, fct = LL.2())
Me4_DRM_TU <- drm(Me4$RGR~Me4$TU50, data = Me4, fct = LL.2())
Me5_DRM_TU <- drm(Me5$RGR~Me5$TU50, data = Me5, fct = LL.2())

B_DRM    <- drm(B$RGR~B$B_sum_TU, data = B, fct = LL.2())
B_IA_DRM <- drm(B$IA~B$B_sum_TU, data = B, fct = LL.2())

T_DRM    <- drm(T$RGR~T$T_sum_TU, data = T, fct = LL.2())
T_IA_DRM <- drm(T$IA~T$T_sum_TU, data = T, fct = LL.2())

Q_DRM    <- drm(Q$RGR~Q$Q_sum_TU, data = Q, fct = LL.2())
Q_IA_DRM <- drm(Q$IA~Q$Q_sum_TU, data = Q, fct = LL.2())

R_DRM    <- drm(R$RGR~R$R_sum_TU, data = R, fct = LL.2())
R_IA_DRM <- drm(R$IA~R$R_sum_TU, data = R, fct = LL.2())

X_DRM    <- drm(X$RGR~X$X_sum_TU, data = X, fct = LL.2())
X_IA_DRM <- drm(X$IA~X$X_sum_TU, data = X, fct = LL.2())

#######################################################################################

#EC AND SLOPE TABLE --> (EC50)
EC_Me1_TU50 <- ED(Me1_DRM_TU, c(10,20,50), interval="delta")
EC_Me2_TU50 <- ED(Me2_DRM_TU, c(10,20,50), interval="delta")
EC_Me3_TU50 <- ED(Me3_DRM_TU, c(10,20,50), interval="delta")
EC_Me4_TU50 <- ED(Me4_DRM_TU, c(10,20,50), interval="delta")
EC_Me5_TU50 <- ED(Me5_DRM_TU, c(10,20,50), interval="delta")

EC_B_TU50 <- ED(B_DRM, c(10,20,50), interval="delta")
EC_T_TU50 <- ED(T_DRM, c(10,20,50), interval="delta")
EC_Q_TU50 <- ED(Q_DRM, c(10,20,50), interval="delta")
EC_R_TU50 <- ED(R_DRM, c(10,20,50), interval="delta")
EC_X_TU50 <- ED(X_DRM, c(10,20,50), interval="delta")

EC10_Me1_TU50 <- EC_Me1_TU50[1:1]
SE10_Me1_TU50 <- EC_Me1_TU50["e:1:10" , "Std. Error"]
EC20_Me1_TU50 <- EC_Me1_TU50["e:1:20" , "Estimate"]
SE20_Me1_TU50 <- EC_Me1_TU50["e:1:20" , "Std. Error"]
EC50_Me1_TU50 <- EC_Me1_TU50["e:1:50" , "Estimate"]
SE50_Me1_TU50 <- EC_Me1_TU50["e:1:50" , "Std. Error"]


EC10_Me2_TU50 <- EC_Me2_TU50[1:1]
SE10_Me2_TU50 <- EC_Me2_TU50["e:1:10" , "Std. Error"]
EC20_Me2_TU50 <- EC_Me2_TU50["e:1:20" , "Estimate"]
SE20_Me2_TU50 <- EC_Me2_TU50["e:1:20" , "Std. Error"]
EC50_Me2_TU50 <- EC_Me2_TU50["e:1:50" , "Estimate"]
SE50_Me2_TU50 <- EC_Me2_TU50["e:1:50" , "Std. Error"]


EC10_Me3_TU50 <- EC_Me3_TU50[1:1]
SE10_Me3_TU50 <- EC_Me3_TU50["e:1:10" , "Std. Error"]
EC20_Me3_TU50 <- EC_Me3_TU50["e:1:20" , "Estimate"]
SE20_Me3_TU50 <- EC_Me3_TU50["e:1:20" , "Std. Error"]
EC50_Me3_TU50 <- EC_Me3_TU50["e:1:50" , "Estimate"]
SE50_Me3_TU50 <- EC_Me3_TU50["e:1:50" , "Std. Error"]


EC10_Me4_TU50 <- EC_Me4_TU50[1:1]
SE10_Me4_TU50 <- EC_Me4_TU50["e:1:10" , "Std. Error"]
EC20_Me4_TU50 <- EC_Me4_TU50["e:1:20" , "Estimate"]
SE20_Me4_TU50 <- EC_Me4_TU50["e:1:20" , "Std. Error"]
EC50_Me4_TU50 <- EC_Me4_TU50["e:1:50" , "Estimate"]
SE50_Me4_TU50 <- EC_Me4_TU50["e:1:50" , "Std. Error"]


EC10_Me5_TU50 <- EC_Me5_TU50[1:1]
SE10_Me5_TU50 <- EC_Me5_TU50["e:1:10" , "Std. Error"]
EC20_Me5_TU50 <- EC_Me5_TU50["e:1:20" , "Estimate"]
SE20_Me5_TU50 <- EC_Me5_TU50["e:1:20" , "Std. Error"]
EC50_Me5_TU50 <- EC_Me5_TU50["e:1:50" , "Estimate"]
SE50_Me5_TU50 <- EC_Me5_TU50["e:1:50" , "Std. Error"]


EC10_B_TU50 <- EC_B_TU50[1:1]
SE10_B_TU50 <- EC_B_TU50["e:1:10" , "Std. Error"]
EC20_B_TU50 <- EC_B_TU50["e:1:20" , "Estimate"]
SE20_B_TU50 <- EC_B_TU50["e:1:20" , "Std. Error"]
EC50_B_TU50 <- EC_B_TU50["e:1:50" , "Estimate"]
SE50_B_TU50 <- EC_B_TU50["e:1:50" , "Std. Error"]
U10_B_TU50 <-  EC_B_TU50["e:1:10" , "Upper"]
L10_B_TU50 <-  EC_B_TU50["e:1:10" , "Lower"]
U20_B_TU50 <-  EC_B_TU50["e:1:20" , "Upper"]
L20_B_TU50 <-  EC_B_TU50["e:1:20" , "Lower"]
U50_B_TU50 <-  EC_B_TU50["e:1:50" , "Upper"]
L50_B_TU50 <-  EC_B_TU50["e:1:50" , "Lower"]


EC10_T_TU50 <- EC_T_TU50[1:1]
SE10_T_TU50 <- EC_T_TU50["e:1:10" , "Std. Error"]
EC20_T_TU50 <- EC_T_TU50["e:1:20" , "Estimate"]
SE20_T_TU50 <- EC_T_TU50["e:1:20" , "Std. Error"]
EC50_T_TU50 <- EC_T_TU50["e:1:50" , "Estimate"]
SE50_T_TU50 <- EC_T_TU50["e:1:50" , "Std. Error"]
U10_T_TU50 <-  EC_T_TU50["e:1:10" , "Upper"]
L10_T_TU50 <-  EC_T_TU50["e:1:10" , "Lower"]
U20_T_TU50 <-  EC_T_TU50["e:1:20" , "Upper"]
L20_T_TU50 <-  EC_T_TU50["e:1:20" , "Lower"]
U50_T_TU50 <-  EC_T_TU50["e:1:50" , "Upper"]
L50_T_TU50 <-  EC_T_TU50["e:1:50" , "Lower"]


EC10_Q_TU50 <- EC_Q_TU50[1:1]
SE10_Q_TU50 <- EC_Q_TU50["e:1:10" , "Std. Error"]
EC20_Q_TU50 <- EC_Q_TU50["e:1:20" , "Estimate"]
SE20_Q_TU50 <- EC_Q_TU50["e:1:20" , "Std. Error"]
EC50_Q_TU50 <- EC_Q_TU50["e:1:50" , "Estimate"]
SE50_Q_TU50 <- EC_Q_TU50["e:1:50" , "Std. Error"]
U10_Q_TU50 <-  EC_Q_TU50["e:1:10" , "Upper"]
L10_Q_TU50 <-  EC_Q_TU50["e:1:10" , "Lower"]
U20_Q_TU50 <-  EC_Q_TU50["e:1:20" , "Upper"]
L20_Q_TU50 <-  EC_Q_TU50["e:1:20" , "Lower"]
U50_Q_TU50 <-  EC_Q_TU50["e:1:50" , "Upper"]
L50_Q_TU50 <-  EC_Q_TU50["e:1:50" , "Lower"]


EC10_R_TU50 <- EC_R_TU50[1:1]
SE10_R_TU50 <- EC_R_TU50["e:1:10" , "Std. Error"]
EC20_R_TU50 <- EC_R_TU50["e:1:20" , "Estimate"]
SE20_R_TU50 <- EC_R_TU50["e:1:20" , "Std. Error"]
EC50_R_TU50 <- EC_R_TU50["e:1:50" , "Estimate"]
SE50_R_TU50 <- EC_R_TU50["e:1:50" , "Std. Error"]
U10_R_TU50 <-  EC_R_TU50["e:1:10" , "Upper"]
L10_R_TU50 <-  EC_R_TU50["e:1:10" , "Lower"]
U20_R_TU50 <-  EC_R_TU50["e:1:20" , "Upper"]
L20_R_TU50 <-  EC_R_TU50["e:1:20" , "Lower"]
U50_R_TU50 <-  EC_R_TU50["e:1:50" , "Upper"]
L50_R_TU50 <-  EC_R_TU50["e:1:50" , "Lower"]


EC10_X_TU50 <- EC_X_TU50[1:1]
SE10_X_TU50 <- EC_X_TU50["e:1:10" , "Std. Error"]
EC20_X_TU50 <- EC_X_TU50["e:1:20" , "Estimate"]
SE20_X_TU50 <- EC_X_TU50["e:1:20" , "Std. Error"]
EC50_X_TU50 <- EC_X_TU50["e:1:50" , "Estimate"]
SE50_X_TU50 <- EC_X_TU50["e:1:50" , "Std. Error"]
U10_X_TU50 <-  EC_X_TU50["e:1:10" , "Upper"]
L10_X_TU50 <-  EC_X_TU50["e:1:10" , "Lower"]
U20_X_TU50 <-  EC_X_TU50["e:1:20" , "Upper"]
L20_X_TU50 <-  EC_X_TU50["e:1:20" , "Lower"]
U50_X_TU50 <-  EC_X_TU50["e:1:50" , "Upper"]
L50_X_TU50 <-  EC_X_TU50["e:1:50" , "Lower"]


Me1_DRM_df_TU50 <-  summary(Me1_DRM_TU)
Me1_DRM_df_TU50 <- Me1_DRM_df_TU50$coefficients
Me1_DRM_df_TU50 <- as.matrix(Me1_DRM_df_TU50, header=TRUE)
Me1_B_TU50      <- Me1_DRM_df_TU50[1:1]
Me1_SEB_TU50    <- Me1_DRM_df_TU50["b:(Intercept)", "Std. Error"]

Me2_DRM_df_TU50 <-  summary(Me2_DRM_TU)
Me2_DRM_df_TU50 <- Me2_DRM_df_TU50$coefficients
Me2_DRM_df_TU50 <- as.matrix(Me2_DRM_df_TU50, header=TRUE)
Me2_B_TU50      <- Me2_DRM_df_TU50[1:1]
Me2_SEB_TU50    <- Me2_DRM_df_TU50["b:(Intercept)", "Std. Error"]

Me3_DRM_df_TU50 <-  summary(Me3_DRM_TU)
Me3_DRM_df_TU50 <- Me3_DRM_df_TU50$coefficients
Me3_DRM_df_TU50 <- as.matrix(Me3_DRM_df_TU50, header=TRUE)
Me3_B_TU50      <- Me3_DRM_df_TU50[1:1]
Me3_SEB_TU50    <- Me3_DRM_df_TU50["b:(Intercept)", "Std. Error"]

Me4_DRM_df_TU50 <-  summary(Me4_DRM_TU)
Me4_DRM_df_TU50 <- Me4_DRM_df_TU50$coefficients
Me4_DRM_df_TU50 <- as.matrix(Me4_DRM_df_TU50, header=TRUE)
Me4_B_TU50      <- Me4_DRM_df_TU50[1:1]
Me4_SEB_TU50    <- Me4_DRM_df_TU50["b:(Intercept)", "Std. Error"]

Me5_DRM_df_TU50 <-  summary(Me5_DRM_TU)
Me5_DRM_df_TU50 <- Me5_DRM_df_TU50$coefficients
Me5_DRM_df_TU50 <- as.matrix(Me5_DRM_df_TU50, header=TRUE)
Me5_B_TU50      <- Me5_DRM_df_TU50[1:1]
Me5_SEB_TU50    <- Me5_DRM_df_TU50["b:(Intercept)", "Std. Error"]

B_DRM_df_TU50 <-  summary(B_DRM)
B_DRM_df_TU50 <- B_DRM_df_TU50$coefficients
B_DRM_df_TU50 <- as.matrix(B_DRM_df_TU50, header=TRUE)
B_B_TU50   <- B_DRM_df_TU50[1:1]
B_SEB_TU50 <- B_DRM_df_TU50["b:(Intercept)", "Std. Error"]

T_DRM_df_TU50 <-  summary(T_DRM)
T_DRM_df_TU50 <- T_DRM_df_TU50$coefficients
T_DRM_df_TU50 <- as.matrix(T_DRM_df_TU50, header=TRUE)
T_B_TU50   <- T_DRM_df_TU50[1:1]
T_SEB_TU50 <- T_DRM_df_TU50["b:(Intercept)", "Std. Error"]

Q_DRM_df_TU50 <-  summary(Q_DRM)
Q_DRM_df_TU50 <- Q_DRM_df_TU50$coefficients
Q_DRM_df_TU50 <- as.matrix(Q_DRM_df_TU50, header=TRUE)
Q_B_TU50   <- Q_DRM_df_TU50[1:1]
Q_SEB_TU50 <- Q_DRM_df_TU50["b:(Intercept)", "Std. Error"]

R_DRM_df_TU50 <-  summary(R_DRM)
R_DRM_df_TU50 <- R_DRM_df_TU50$coefficients
R_DRM_df_TU50 <- as.matrix(R_DRM_df_TU50, header=TRUE)
R_B_TU50   <- R_DRM_df_TU50[1:1]
R_SEB_TU50 <- R_DRM_df_TU50["b:(Intercept)", "Std. Error"]

X_DRM_df_TU50 <-  summary(X_DRM)
X_DRM_df_TU50 <- X_DRM_df_TU50$coefficients
X_DRM_df_TU50 <- as.matrix(X_DRM_df_TU50, header=TRUE)
X_B_TU50   <- X_DRM_df_TU50[1:1]
X_SEB_TU50 <- X_DRM_df_TU50["b:(Intercept)", "Std. Error"]


Results_summary_TU50 <- data.frame(Me1 = numeric(14), Me2 = numeric(14), Me3 = numeric(14), Me4 = numeric(14), Me5 = numeric(14),
                                   B = numeric(14), T = numeric(14), Q = numeric(14), R = numeric(14), X = numeric(14),
                                   row.names = c("EC10", "EC20", "EC50", "slope", "SE_slope", "SE_EC10", "SE_EC20","SE_EC50", 
                                                 "Upper_EC10", "Lower_EC10", "Upper_EC20", "Lower_EC20", "Upper_EC50", "Lower_EC50"))

Results_summary_TU50[1, "Me1"] <- as.numeric(EC10_Me1_TU50)
Results_summary_TU50[1, "Me2"] <- as.numeric(EC10_Me2_TU50)
Results_summary_TU50[1, "Me3"] <- as.numeric(EC10_Me3_TU50)
Results_summary_TU50[1, "Me4"] <- as.numeric(EC10_Me4_TU50)
Results_summary_TU50[1, "Me5"] <- as.numeric(EC10_Me5_TU50)
Results_summary_TU50[1, "B"]  <- as.numeric(EC10_B_TU50)
Results_summary_TU50[1, "T"]  <- as.numeric(EC10_T_TU50)
Results_summary_TU50[1, "Q"]  <- as.numeric(EC10_Q_TU50)
Results_summary_TU50[1, "R"]  <- as.numeric(EC10_R_TU50)
Results_summary_TU50[1, "X"]  <- as.numeric(EC10_X_TU50)

Results_summary_TU50[2, "Me1"] <- as.numeric(EC20_Me1_TU50)
Results_summary_TU50[2, "Me2"] <- as.numeric(EC20_Me2_TU50)
Results_summary_TU50[2, "Me3"] <- as.numeric(EC20_Me3_TU50)
Results_summary_TU50[2, "Me4"] <- as.numeric(EC20_Me4_TU50)
Results_summary_TU50[2, "Me5"] <- as.numeric(EC20_Me5_TU50)
Results_summary_TU50[2, "B"]  <- as.numeric(EC20_B_TU50)
Results_summary_TU50[2, "T"]  <- as.numeric(EC20_T_TU50)
Results_summary_TU50[2, "Q"]  <- as.numeric(EC20_Q_TU50)
Results_summary_TU50[2, "R"]  <- as.numeric(EC20_R_TU50)
Results_summary_TU50[2, "X"]  <- as.numeric(EC20_X_TU50)

Results_summary_TU50[3, "Me1"] <- as.numeric(EC50_Me1_TU50)
Results_summary_TU50[3, "Me2"] <- as.numeric(EC50_Me2_TU50)
Results_summary_TU50[3, "Me3"] <- as.numeric(EC50_Me3_TU50)
Results_summary_TU50[3, "Me4"] <- as.numeric(EC50_Me4_TU50)
Results_summary_TU50[3, "Me5"] <- as.numeric(EC50_Me5_TU50)
Results_summary_TU50[3, "B"]  <- as.numeric(EC50_B_TU50)
Results_summary_TU50[3, "T"]  <- as.numeric(EC50_T_TU50)
Results_summary_TU50[3, "Q"]  <- as.numeric(EC50_Q_TU50)
Results_summary_TU50[3, "R"]  <- as.numeric(EC50_R_TU50)
Results_summary_TU50[3, "X"]  <- as.numeric(EC50_X_TU50)

Results_summary_TU50[4, "Me1"] <- as.numeric(Me1_B_TU50)
Results_summary_TU50[4, "Me2"] <- as.numeric(Me2_B_TU50)
Results_summary_TU50[4, "Me3"] <- as.numeric(Me3_B_TU50)
Results_summary_TU50[4, "Me4"] <- as.numeric(Me4_B_TU50)
Results_summary_TU50[4, "Me5"] <- as.numeric(Me5_B_TU50)
Results_summary_TU50[4, "B"]  <- as.numeric(B_B_TU50)
Results_summary_TU50[4, "T"]  <- as.numeric(T_B_TU50)
Results_summary_TU50[4, "Q"]  <- as.numeric(Q_B_TU50)
Results_summary_TU50[4, "R"]  <- as.numeric(R_B_TU50)
Results_summary_TU50[4, "X"]  <- as.numeric(X_B_TU50)

Results_summary_TU50[5, "Me1"] <- as.numeric(Me1_SEB_TU50)
Results_summary_TU50[5, "Me2"] <- as.numeric(Me2_SEB_TU50)
Results_summary_TU50[5, "Me3"] <- as.numeric(Me3_SEB_TU50)
Results_summary_TU50[5, "Me4"] <- as.numeric(Me4_SEB_TU50)
Results_summary_TU50[5, "Me5"] <- as.numeric(Me5_SEB_TU50)
Results_summary_TU50[5, "B"]  <- as.numeric(B_SEB_TU50)
Results_summary_TU50[5, "T"]  <- as.numeric(T_SEB_TU50)
Results_summary_TU50[5, "Q"]  <- as.numeric(Q_SEB_TU50)
Results_summary_TU50[5, "R"]  <- as.numeric(R_SEB_TU50)
Results_summary_TU50[5, "X"]  <- as.numeric(X_SEB_TU50)

Results_summary_TU50[6, "Me1"] <- as.numeric(SE10_Me1_TU50)
Results_summary_TU50[6, "Me2"] <- as.numeric(SE10_Me2_TU50)
Results_summary_TU50[6, "Me3"] <- as.numeric(SE10_Me3_TU50)
Results_summary_TU50[6, "Me4"] <- as.numeric(SE10_Me4_TU50)
Results_summary_TU50[6, "Me5"] <- as.numeric(SE10_Me5_TU50)
Results_summary_TU50[6, "B"]  <- as.numeric(SE10_B_TU50)
Results_summary_TU50[6, "T"]  <- as.numeric(SE10_T_TU50)
Results_summary_TU50[6, "Q"]  <- as.numeric(SE10_Q_TU50)
Results_summary_TU50[6, "R"]  <- as.numeric(SE10_R_TU50)
Results_summary_TU50[6, "X"]  <- as.numeric(SE10_X_TU50)

Results_summary_TU50[7, "Me1"] <- as.numeric(SE20_Me1_TU50)
Results_summary_TU50[7, "Me2"] <- as.numeric(SE20_Me2_TU50)
Results_summary_TU50[7, "Me3"] <- as.numeric(SE20_Me3_TU50)
Results_summary_TU50[7, "Me4"] <- as.numeric(SE20_Me4_TU50)
Results_summary_TU50[7, "Me5"] <- as.numeric(SE20_Me5_TU50)
Results_summary_TU50[7, "B"]  <- as.numeric(SE20_B_TU50)
Results_summary_TU50[7, "T"]  <- as.numeric(SE20_T_TU50)
Results_summary_TU50[7, "Q"]  <- as.numeric(SE20_Q_TU50)
Results_summary_TU50[7, "R"]  <- as.numeric(SE20_R_TU50)
Results_summary_TU50[7, "X"]  <- as.numeric(SE20_X_TU50)

Results_summary_TU50[8, "Me1"] <- as.numeric(SE50_Me1_TU50)
Results_summary_TU50[8, "Me2"] <- as.numeric(SE50_Me2_TU50)
Results_summary_TU50[8, "Me3"] <- as.numeric(SE50_Me3_TU50)
Results_summary_TU50[8, "Me4"] <- as.numeric(SE50_Me4_TU50)
Results_summary_TU50[8, "Me5"] <- as.numeric(SE50_Me5_TU50)
Results_summary_TU50[8, "B"]  <- as.numeric(SE50_B_TU50)
Results_summary_TU50[8, "T"]  <- as.numeric(SE50_T_TU50)
Results_summary_TU50[8, "Q"]  <- as.numeric(SE50_Q_TU50)
Results_summary_TU50[8, "R"]  <- as.numeric(SE50_R_TU50)
Results_summary_TU50[8, "X"]  <- as.numeric(SE50_X_TU50)

Results_summary_TU50[9, "B"]  <- as.numeric(U10_B_TU50)
Results_summary_TU50[9, "T"]  <- as.numeric(U10_T_TU50)
Results_summary_TU50[9, "Q"]  <- as.numeric(U10_Q_TU50)
Results_summary_TU50[9, "R"]  <- as.numeric(U10_R_TU50)
Results_summary_TU50[9, "X"]  <- as.numeric(U10_X_TU50)

Results_summary_TU50[10, "B"]  <- as.numeric(L10_B_TU50)
Results_summary_TU50[10, "T"]  <- as.numeric(L10_T_TU50)
Results_summary_TU50[10, "Q"]  <- as.numeric(L10_Q_TU50)
Results_summary_TU50[10, "R"]  <- as.numeric(L10_R_TU50)
Results_summary_TU50[10, "X"]  <- as.numeric(L10_X_TU50)

Results_summary_TU50[11, "B"]  <- as.numeric(U20_B_TU50)
Results_summary_TU50[11, "T"]  <- as.numeric(U20_T_TU50)
Results_summary_TU50[11, "Q"]  <- as.numeric(U20_Q_TU50)
Results_summary_TU50[11, "R"]  <- as.numeric(U20_R_TU50)
Results_summary_TU50[11, "X"]  <- as.numeric(U20_X_TU50)

Results_summary_TU50[12, "B"]  <- as.numeric(L20_B_TU50)
Results_summary_TU50[12, "T"]  <- as.numeric(L20_T_TU50)
Results_summary_TU50[12, "Q"]  <- as.numeric(L20_Q_TU50)
Results_summary_TU50[12, "R"]  <- as.numeric(L20_R_TU50)
Results_summary_TU50[12, "X"]  <- as.numeric(L20_X_TU50)

Results_summary_TU50[13, "B"]  <- as.numeric(U50_B_TU50)
Results_summary_TU50[13, "T"]  <- as.numeric(U50_T_TU50)
Results_summary_TU50[13, "Q"]  <- as.numeric(U50_Q_TU50)
Results_summary_TU50[13, "R"]  <- as.numeric(U50_R_TU50)
Results_summary_TU50[13, "X"]  <- as.numeric(U50_X_TU50)

Results_summary_TU50[14, "B"]  <- as.numeric(L50_B_TU50)
Results_summary_TU50[14, "T"]  <- as.numeric(L50_T_TU50)
Results_summary_TU50[14, "Q"]  <- as.numeric(L50_Q_TU50)
Results_summary_TU50[14, "R"]  <- as.numeric(L50_R_TU50)
Results_summary_TU50[14, "X"]  <- as.numeric(L50_X_TU50)


Results_summary_TU50 <- Results_summary_TU50 %>% 
  rename(
    As = Me1,
    Zn = Me2,
    Pb = Me3,
    Ag = Me4,
    Cu = Me5
  )


write.table(Results_summary_TU50, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Summaries//Results_summary_TU50.txt", 
            sep = "\t", quote = FALSE, row.names = TRUE)


#--------------------------------------------------------------------------------------------------------
# CA CALCULATION
#--------------------------------------------------------------------------------------------------------
EC50_B <- c(EC50Me1, EC50Me2)
EC50_T <- c(EC50Me1, EC50Me2, EC50Me3)
EC50_Q <- c(EC50Me1, EC50Me2, EC50Me3, EC50Me4)
EC50_R <- c(EC50Me1, EC50Me2, EC50Me3, EC50Me4, EC50Me5)
EC50_X <- c(EC50Me1, EC50Me2, EC50Me3, EC50Me4, EC50Me5)

B_B <- c(BMe1, BMe2)
B_T <- c(BMe1, BMe2, BMe3)
B_Q <- c(BMe1, BMe2, BMe3, BMe4)
B_R <- c(BMe1, BMe2, BMe3, BMe4, BMe5)
B_X <- c(BMe1, BMe2, BMe3, BMe4, BMe5)

Clist_B <- as.list(as.data.frame(t(smaller_B[, 1:2])))  # Define mixture concentrations for metal n the same order everytime
Clist_T <- as.list(as.data.frame(t(smaller_T[, 1:3])))
Clist_Q <- as.list(as.data.frame(t(smaller_Q[, 1:4])))
Clist_R <- as.list(as.data.frame(t(smaller_R[, 1:5])))
Clist_X <- as.list(as.data.frame(t(smaller_X[, 1:5])))

-#CA BINARY#--------------------------------------------------------------------------------------------------------

# Define CA function
CAcalc_B <-function(x,C_B,EC50_B,B_B) {  # Function will return the sum TU for a given x
  Sum <- 0
  for (metal in 1:length(C_B))  # Loop allows to run over all metals for which the information is provided
  {
    Sum <- Sum + C_B[metal]/(x^(1/B_B[metal])*EC50_B[metal])  # CA equation
  }
  return(abs(1-Sum))  # Return absolute value of 1-Sum because optimize searches will search for the minimum value
}


output_B <- data.frame(matrix(ncol=2,nrow=length(Clist_B))) # Make empty list to fill with output per mixture
colnames(output_B) <- c("CA","x")
for (mixture in 1:length(Clist_B)) {
  C_B <- Clist_B[[mixture]]
  # Look for the x that leads to Sum = 1 i.e. optimize the function 
  xcalc <- optimize(CAcalc_B,interval=c(1e-6,1e4),C=C_B,EC50=EC50_B,B=B_B)
  xmin <- xcalc$minimum
  output_B[mixture,] <- c(100/(xmin+1),xmin)
}


output_B <- output_B[rep(1:nrow(output_B), each = 3), ]
output_B <-  output_B/100
output_B <-  rbind(c(1, 0), output_B)
B$CA <- output_B[, 1]

-#CA TERNARY#--------------------------------------------------------------------------------------------------------

# Define CA function
CAcalc_T <-function(x,C_T,EC50_T,B_T) {  # Function will return the sum TU for a given x
  Sum <- 0
  for (metal in 1:length(C_T))  # Loop allows to run over all metals for which the information is provided
  {
    Sum <- Sum + C_T[metal]/(x^(1/B_T[metal])*EC50_T[metal])  # CA equation
  }
  return(abs(1-Sum))  # Return absolute value of 1-Sum because optimize searches will search for the minimum value
}


output_T <- data.frame(matrix(ncol=2,nrow=length(Clist_T))) # Make empty list to fill with output per mixture
colnames(output_T) <- c("CA","x")
for (mixture in 1:length(Clist_T)) {
  C_T <- Clist_T[[mixture]]
  # Look for the x that leads to Sum = 1 i.e. optimize the function 
  xcalc <- optimize(CAcalc_T,interval=c(1e-6,1e4),C=C_T,EC50=EC50_T,B=B_T)
  xmin <- xcalc$minimum
  output_T[mixture,] <- c(100/(xmin+1),xmin)
}

output_T <- output_T[rep(1:nrow(output_T), each = 3), ]
output_T <-  output_T/100
output_T <-  rbind(c(1, 0), output_T)
T$CA <- output_T[, 1]

-#CA QUATERNARY#--------------------------------------------------------------------------------------------------------

# Define CA function
CAcalc_Q <-function(x,C_Q,EC50_Q,B_Q) {  # Function will return the sum TU for a given x
  Sum <- 0
  for (metal in 1:length(C_Q))  # Loop allows to run over all metals for which the information is provided
  {
    Sum <- Sum + C_Q[metal]/(x^(1/B_Q[metal])*EC50_Q[metal])  # CA equation
  }
  return(abs(1-Sum))  # Return absolute value of 1-Sum because optimize searches will search for the minimum value
}


output_Q <- data.frame(matrix(ncol=2,nrow=length(Clist_Q))) # Make empty list to fill with output per mixture
colnames(output_Q) <- c("CA","x")
for (mixture in 1:length(Clist_Q)) {
  C_Q <- Clist_Q[[mixture]]
  # Look for the x that leads to Sum = 1 i.e. optimize the function 
  xcalc <- optimize(CAcalc_Q,interval=c(1e-6,1e4),C=C_Q,EC50=EC50_Q,B=B_Q)
  xmin <- xcalc$minimum
  output_Q[mixture,] <- c(100/(xmin+1),xmin)
}


output_Q <- output_Q[rep(1:nrow(output_Q), each = 3), ]
output_Q <-  output_Q/100
output_Q <-  rbind(c(1, 0), output_Q)
Q$CA <- output_Q[, 1]

-#CA ENV REL#--------------------------------------------------------------------------------------------------------

# Define CA function
CAcalc_R <-function(x,C_R,EC50_R,B_R) {  # Function will return the sum TU for a given x
  Sum <- 0
  for (metal in 1:length(C_R))  # Loop allows to run over all metals for which the information is provided
  {
    Sum <- Sum + C_R[metal]/(x^(1/B_R[metal])*EC50_R[metal])  # CA equation
  }
  return(abs(1-Sum))  # Return absolute value of 1-Sum because optimize searches will search for the minimum value
}


output_R <- data.frame(matrix(ncol=2,nrow=length(Clist_R))) # Make empty list to fill with output per mixture
colnames(output_R) <- c("CA","x")
for (mixture in 1:length(Clist_R)) {
  C_R <- Clist_R[[mixture]]
  # Look for the x that leads to Sum = 1 i.e. optimize the function 
  xcalc <- optimize(CAcalc_R,interval=c(1e-6,1e4),C=C_R,EC50=EC50_R,B=B_R)
  xmin <- xcalc$minimum
  output_R[mixture,] <- c(100/(xmin+1),xmin)
}


output_R <- output_R[rep(1:nrow(output_R), each = 3), ]
output_R <-  output_R/100
output_R <-  rbind(c(1, 0), output_R)
R$CA <- output_R[, 1]

-#CA EQUITOXIC#--------------------------------------------------------------------------------------------------------

# Define CA function
CAcalc_X <-function(x,C_X,EC50_X,B_X) {  # Function will return the sum TU for a given x
  Sum <- 0
  for (metal in 1:length(C_X))  # Loop allows to run over all metals for which the information is provided
  {
    Sum <- Sum + C_X[metal]/(x^(1/B_X[metal])*EC50_X[metal])  # CA equation
  }
  return(abs(1-Sum))  # Return absolute value of 1-Sum because optimize searches will search for the minimum value
}


output_X <- data.frame(matrix(ncol=2,nrow=length(Clist_X))) # Make empty list to fill with output per mixture
colnames(output_X) <- c("CA","x")
for (mixture in 1:length(Clist_X)) {
  C_X <- Clist_X[[mixture]]
  # Look for the x that leads to Sum = 1 i.e. optimize the function 
  xcalc <- optimize(CAcalc_X,interval=c(1e-6,1e4),C=C_X,EC50=EC50_X,B=B_X)
  xmin <- xcalc$minimum
  output_X[mixture,] <- c(100/(xmin+1),xmin)
}


output_X <- output_X[rep(1:nrow(output_X), each = 3), ]
output_X <-  output_X/100
output_X <-  rbind(c(1, 0), output_X)
X$CA <- output_X[, 1]

# DRM -----------------------------------------------------------------
B_CA_DRM <- drm(B$CA~B$B_sum_TU, data = B, fct = LL.2())
T_CA_DRM <- drm(T$CA~T$T_sum_TU, data = T, fct = LL.2())
Q_CA_DRM <- drm(Q$CA~Q$Q_sum_TU, data = Q, fct = LL.2())
R_CA_DRM <- drm(R$CA~R$R_sum_TU, data = R, fct = LL.2())
X_CA_DRM <- drm(X$CA~X$X_sum_TU, data = X, fct = LL.2())

#############################################################################
#EC10 
#############################################################################

Me1_DRM_TU10 <- drm(Me1$RGR~Me1$TU10, data = Me1, fct = LL.2())
Me2_DRM_TU10 <- drm(Me2$RGR~Me2$TU10, data = Me2, fct = LL.2())
Me3_DRM_TU10 <- drm(Me3$RGR~Me3$TU10, data = Me3, fct = LL.2())
Me4_DRM_TU10 <- drm(Me4$RGR~Me4$TU10, data = Me4, fct = LL.2())
Me5_DRM_TU10 <- drm(Me5$RGR~Me5$TU10, data = Me5, fct = LL.2())

B_DRM_TU10    <- drm(B$RGR~B$B_sum_TU10, data = B, fct = LL.2())
B_IA_DRM_TU10 <- drm(B$IA~B$B_sum_TU10,  data = B, fct = LL.2())
B_CA_DRM_TU10 <- drm(B$CA~B$B_sum_TU10,  data = B, fct = LL.2())

T_DRM_TU10    <- drm(T$RGR~T$T_sum_TU10, data = T, fct = LL.2())
T_IA_DRM_TU10 <- drm(T$IA~T$T_sum_TU10,  data = T, fct = LL.2())
T_CA_DRM_TU10 <- drm(T$CA~T$T_sum_TU10,  data = T, fct = LL.2())

Q_DRM_TU10    <- drm(Q$RGR~Q$Q_sum_TU10, data = Q, fct = LL.2())
Q_IA_DRM_TU10 <- drm(Q$IA~Q$Q_sum_TU10,  data = Q, fct = LL.2())
Q_CA_DRM_TU10 <- drm(Q$CA~Q$Q_sum_TU10,  data = Q, fct = LL.2())

R_DRM_TU10    <- drm(R$RGR~R$R_sum_TU10, data = R, fct = LL.2())
R_IA_DRM_TU10 <- drm(R$IA~R$R_sum_TU10,  data = R, fct = LL.2())
R_CA_DRM_TU10 <- drm(R$CA~R$R_sum_TU10,  data = R, fct = LL.2())

X_DRM_TU10    <- drm(X$RGR~X$X_sum_TU10, data = X, fct = LL.2())
X_IA_DRM_TU10 <- drm(X$IA~X$X_sum_TU10,  data = X, fct = LL.2())
X_CA_DRM_TU10 <- drm(X$CA~X$X_sum_TU10,  data = X, fct = LL.2())


#EC AND SLOPE TABLE --> (EC10)
EC_B_TU10 <- ED(B_DRM_TU10, c(10,20,50), interval="delta")
EC_T_TU10 <- ED(T_DRM_TU10, c(10,20,50), interval="delta")
EC_Q_TU10 <- ED(Q_DRM_TU10, c(10,20,50), interval="delta")
EC_R_TU10 <- ED(R_DRM_TU10, c(10,20,50), interval="delta")
EC_X_TU10 <- ED(X_DRM_TU10, c(10,20,50), interval="delta")


EC10_B_TU10 <- EC_B_TU10[1:1]
SE10_B_TU10 <- EC_B_TU10["e:1:10" , "Std. Error"]
EC20_B_TU10 <- EC_B_TU10["e:1:20" , "Estimate"]
SE20_B_TU10 <- EC_B_TU10["e:1:20" , "Std. Error"]
EC50_B_TU10 <- EC_B_TU10["e:1:50" , "Estimate"]
SE50_B_TU10 <- EC_B_TU10["e:1:50" , "Std. Error"]
U10_B_TU10 <-  EC_B_TU10["e:1:10" , "Upper"]
L10_B_TU10 <-  EC_B_TU10["e:1:10" , "Lower"]
U20_B_TU10 <-  EC_B_TU10["e:1:20" , "Upper"]
L20_B_TU10 <-  EC_B_TU10["e:1:20" , "Lower"]
U50_B_TU10 <-  EC_B_TU10["e:1:50" , "Upper"]
L50_B_TU10 <-  EC_B_TU10["e:1:50" , "Lower"]


EC10_T_TU10 <- EC_T_TU10[1:1]
SE10_T_TU10 <- EC_T_TU10["e:1:10" , "Std. Error"]
EC20_T_TU10 <- EC_T_TU10["e:1:20" , "Estimate"]
SE20_T_TU10 <- EC_T_TU10["e:1:20" , "Std. Error"]
EC50_T_TU10 <- EC_T_TU10["e:1:50" , "Estimate"]
SE50_T_TU10 <- EC_T_TU10["e:1:50" , "Std. Error"]
U10_T_TU10 <-  EC_T_TU10["e:1:10" , "Upper"]
L10_T_TU10 <-  EC_T_TU10["e:1:10" , "Lower"]
U20_T_TU10 <-  EC_T_TU10["e:1:20" , "Upper"]
L20_T_TU10 <-  EC_T_TU10["e:1:20" , "Lower"]
U50_T_TU10 <-  EC_T_TU10["e:1:50" , "Upper"]
L50_T_TU10 <-  EC_T_TU10["e:1:50" , "Lower"]


EC10_Q_TU10 <- EC_Q_TU10[1:1]
SE10_Q_TU10 <- EC_Q_TU10["e:1:10" , "Std. Error"]
EC20_Q_TU10 <- EC_Q_TU10["e:1:20" , "Estimate"]
SE20_Q_TU10 <- EC_Q_TU10["e:1:20" , "Std. Error"]
EC50_Q_TU10 <- EC_Q_TU10["e:1:50" , "Estimate"]
SE50_Q_TU10 <- EC_Q_TU10["e:1:50" , "Std. Error"]
U10_Q_TU10 <-  EC_Q_TU10["e:1:10" , "Upper"]
L10_Q_TU10 <-  EC_Q_TU10["e:1:10" , "Lower"]
U20_Q_TU10 <-  EC_Q_TU10["e:1:20" , "Upper"]
L20_Q_TU10 <-  EC_Q_TU10["e:1:20" , "Lower"]
U50_Q_TU10 <-  EC_Q_TU10["e:1:50" , "Upper"]
L50_Q_TU10 <-  EC_Q_TU10["e:1:50" , "Lower"]


EC10_R_TU10 <- EC_R_TU10[1:1]
SE10_R_TU10 <- EC_R_TU10["e:1:10" , "Std. Error"]
EC20_R_TU10 <- EC_R_TU10["e:1:20" , "Estimate"]
SE20_R_TU10 <- EC_R_TU10["e:1:20" , "Std. Error"]
EC50_R_TU10 <- EC_R_TU10["e:1:50" , "Estimate"]
SE50_R_TU10 <- EC_R_TU10["e:1:50" , "Std. Error"]
U10_R_TU10 <-  EC_R_TU10["e:1:10" , "Upper"]
L10_R_TU10 <-  EC_R_TU10["e:1:10" , "Lower"]
U20_R_TU10 <-  EC_R_TU10["e:1:20" , "Upper"]
L20_R_TU10 <-  EC_R_TU10["e:1:20" , "Lower"]
U50_R_TU10 <-  EC_R_TU10["e:1:50" , "Upper"]
L50_R_TU10 <-  EC_R_TU10["e:1:50" , "Lower"]


EC10_X_TU10 <- EC_X_TU10[1:1]
SE10_X_TU10 <- EC_X_TU10["e:1:10" , "Std. Error"]
EC20_X_TU10 <- EC_X_TU10["e:1:20" , "Estimate"]
SE20_X_TU10 <- EC_X_TU10["e:1:20" , "Std. Error"]
EC50_X_TU10 <- EC_X_TU10["e:1:50" , "Estimate"]
SE50_X_TU10 <- EC_X_TU10["e:1:50" , "Std. Error"]
U10_X_TU10 <-  EC_X_TU10["e:1:10" , "Upper"]
L10_X_TU10 <-  EC_X_TU10["e:1:10" , "Lower"]
U20_X_TU10 <-  EC_X_TU10["e:1:20" , "Upper"]
L20_X_TU10 <-  EC_X_TU10["e:1:20" , "Lower"]
U50_X_TU10 <-  EC_X_TU10["e:1:50" , "Upper"]
L50_X_TU10 <-  EC_X_TU10["e:1:50" , "Lower"]


B_DRM_df_TU10 <-  summary(B_DRM_TU10)
B_DRM_df_TU10 <- B_DRM_df_TU10$coefficients
B_DRM_df_TU10 <- as.matrix(B_DRM_df_TU10, header=TRUE)
B_B_TU10   <- B_DRM_df_TU10[1:1]
B_SEB_TU10 <- B_DRM_df_TU10["b:(Intercept)", "Std. Error"]

T_DRM_df_TU10 <-  summary(T_DRM_TU10)
T_DRM_df_TU10 <- T_DRM_df_TU10$coefficients
T_DRM_df_TU10 <- as.matrix(T_DRM_df_TU10, header=TRUE)
T_B_TU10   <- T_DRM_df_TU10[1:1]
T_SEB_TU10 <- T_DRM_df_TU10["b:(Intercept)", "Std. Error"]

Q_DRM_df_TU10 <-  summary(Q_DRM_TU10)
Q_DRM_df_TU10 <- Q_DRM_df_TU10$coefficients
Q_DRM_df_TU10 <- as.matrix(Q_DRM_df_TU10, header=TRUE)
Q_B_TU10   <- Q_DRM_df_TU10[1:1]
Q_SEB_TU10 <- Q_DRM_df_TU10["b:(Intercept)", "Std. Error"]

R_DRM_df_TU10 <-  summary(R_DRM_TU10)
R_DRM_df_TU10 <- R_DRM_df_TU10$coefficients
R_DRM_df_TU10 <- as.matrix(R_DRM_df_TU10, header=TRUE)
R_B_TU10   <- R_DRM_df_TU10[1:1]
R_SEB_TU10 <- R_DRM_df_TU10["b:(Intercept)", "Std. Error"]

X_DRM_df_TU10 <-  summary(X_DRM_TU10)
X_DRM_df_TU10 <- X_DRM_df_TU10$coefficients
X_DRM_df_TU10 <- as.matrix(X_DRM_df_TU10, header=TRUE)
X_B_TU10   <- X_DRM_df_TU10[1:1]
X_SEB_TU10 <- X_DRM_df_TU10["b:(Intercept)", "Std. Error"]


Results_summary_TU10 <- data.frame(B = numeric(14), T = numeric(14), Q = numeric(14), R = numeric(14), X = numeric(14),
                                   row.names = c("EC10", "EC20", "EC50", "slope", "SE_slope", "SE_EC10", "SE_EC20","SE_EC50", 
                                                 "Upper_EC10", "Lower_EC10", "Upper_EC20", "Lower_EC20", "Upper_EC50", "Lower_EC50"))

Results_summary_TU10[1, "B"]  <- as.numeric(EC10_B_TU10)
Results_summary_TU10[1, "T"]  <- as.numeric(EC10_T_TU10)
Results_summary_TU10[1, "Q"]  <- as.numeric(EC10_Q_TU10)
Results_summary_TU10[1, "R"]  <- as.numeric(EC10_R_TU10)
Results_summary_TU10[1, "X"]  <- as.numeric(EC10_X_TU10)

Results_summary_TU10[2, "B"]  <- as.numeric(EC20_B_TU10)
Results_summary_TU10[2, "T"]  <- as.numeric(EC20_T_TU10)
Results_summary_TU10[2, "Q"]  <- as.numeric(EC20_Q_TU10)
Results_summary_TU10[2, "R"]  <- as.numeric(EC20_R_TU10)
Results_summary_TU10[2, "X"]  <- as.numeric(EC20_X_TU10)

Results_summary_TU10[3, "B"]  <- as.numeric(EC50_B_TU10)
Results_summary_TU10[3, "T"]  <- as.numeric(EC50_T_TU10)
Results_summary_TU10[3, "Q"]  <- as.numeric(EC50_Q_TU10)
Results_summary_TU10[3, "R"]  <- as.numeric(EC50_R_TU10)
Results_summary_TU10[3, "X"]  <- as.numeric(EC50_X_TU10)

Results_summary_TU10[4, "B"]  <- as.numeric(B_B_TU10)
Results_summary_TU10[4, "T"]  <- as.numeric(T_B_TU10)
Results_summary_TU10[4, "Q"]  <- as.numeric(Q_B_TU10)
Results_summary_TU10[4, "R"]  <- as.numeric(R_B_TU10)
Results_summary_TU10[4, "X"]  <- as.numeric(X_B_TU10)

Results_summary_TU10[5, "B"]  <- as.numeric(B_SEB_TU10)
Results_summary_TU10[5, "T"]  <- as.numeric(T_SEB_TU10)
Results_summary_TU10[5, "Q"]  <- as.numeric(Q_SEB_TU10)
Results_summary_TU10[5, "R"]  <- as.numeric(R_SEB_TU10)
Results_summary_TU10[5, "X"]  <- as.numeric(X_SEB_TU10)

Results_summary_TU10[6, "B"]  <- as.numeric(SE10_B_TU10)
Results_summary_TU10[6, "T"]  <- as.numeric(SE10_T_TU10)
Results_summary_TU10[6, "Q"]  <- as.numeric(SE10_Q_TU10)
Results_summary_TU10[6, "R"]  <- as.numeric(SE10_R_TU10)
Results_summary_TU10[6, "X"]  <- as.numeric(SE10_X_TU10)

Results_summary_TU10[7, "B"]  <- as.numeric(SE20_B_TU10)
Results_summary_TU10[7, "T"]  <- as.numeric(SE20_T_TU10)
Results_summary_TU10[7, "Q"]  <- as.numeric(SE20_Q_TU10)
Results_summary_TU10[7, "R"]  <- as.numeric(SE20_R_TU10)
Results_summary_TU10[7, "X"]  <- as.numeric(SE20_X_TU10)

Results_summary_TU10[8, "B"]  <- as.numeric(SE50_B_TU10)
Results_summary_TU10[8, "T"]  <- as.numeric(SE50_T_TU10)
Results_summary_TU10[8, "Q"]  <- as.numeric(SE50_Q_TU10)
Results_summary_TU10[8, "R"]  <- as.numeric(SE50_R_TU10)
Results_summary_TU10[8, "X"]  <- as.numeric(SE50_X_TU10)

Results_summary_TU10[9, "B"]  <- as.numeric(U10_B_TU10)
Results_summary_TU10[9, "T"]  <- as.numeric(U10_T_TU10)
Results_summary_TU10[9, "Q"]  <- as.numeric(U10_Q_TU10)
Results_summary_TU10[9, "R"]  <- as.numeric(U10_R_TU10)
Results_summary_TU10[9, "X"]  <- as.numeric(U10_X_TU10)

Results_summary_TU10[10, "B"]  <- as.numeric(L10_B_TU10)
Results_summary_TU10[10, "T"]  <- as.numeric(L10_T_TU10)
Results_summary_TU10[10, "Q"]  <- as.numeric(L10_Q_TU10)
Results_summary_TU10[10, "R"]  <- as.numeric(L10_R_TU10)
Results_summary_TU10[10, "X"]  <- as.numeric(L10_X_TU10)

Results_summary_TU10[11, "B"]  <- as.numeric(U20_B_TU10)
Results_summary_TU10[11, "T"]  <- as.numeric(U20_T_TU10)
Results_summary_TU10[11, "Q"]  <- as.numeric(U20_Q_TU10)
Results_summary_TU10[11, "R"]  <- as.numeric(U20_R_TU10)
Results_summary_TU10[11, "X"]  <- as.numeric(U20_X_TU10)

Results_summary_TU10[12, "B"]  <- as.numeric(L20_B_TU10)
Results_summary_TU10[12, "T"]  <- as.numeric(L20_T_TU10)
Results_summary_TU10[12, "Q"]  <- as.numeric(L20_Q_TU10)
Results_summary_TU10[12, "R"]  <- as.numeric(L20_R_TU10)
Results_summary_TU10[12, "X"]  <- as.numeric(L20_X_TU10)

Results_summary_TU10[13, "B"]  <- as.numeric(U50_B_TU10)
Results_summary_TU10[13, "T"]  <- as.numeric(U50_T_TU10)
Results_summary_TU10[13, "Q"]  <- as.numeric(U50_Q_TU10)
Results_summary_TU10[13, "R"]  <- as.numeric(U50_R_TU10)
Results_summary_TU10[13, "X"]  <- as.numeric(U50_X_TU10)

Results_summary_TU10[14, "B"]  <- as.numeric(L50_B_TU10)
Results_summary_TU10[14, "T"]  <- as.numeric(L50_T_TU10)
Results_summary_TU10[14, "Q"]  <- as.numeric(L50_Q_TU10)
Results_summary_TU10[14, "R"]  <- as.numeric(L50_R_TU10)
Results_summary_TU10[14, "X"]  <- as.numeric(L50_X_TU10)


write.table(Results_summary_TU10, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Summaries//Results_summary_TU10.txt", 
            sep = "\t", quote = FALSE, row.names = TRUE)

#############################################################################
#EC20 
#############################################################################

Me1_DRM_TU20 <- drm(Me1$RGR~Me1$TU20, data = Me1, fct = LL.2())
Me2_DRM_TU20 <- drm(Me2$RGR~Me2$TU20, data = Me2, fct = LL.2())
Me3_DRM_TU20 <- drm(Me3$RGR~Me3$TU20, data = Me3, fct = LL.2())
Me4_DRM_TU20 <- drm(Me4$RGR~Me4$TU20, data = Me4, fct = LL.2())
Me5_DRM_TU20 <- drm(Me5$RGR~Me5$TU20, data = Me5, fct = LL.2())

B_DRM_TU20    <- drm(B$RGR~B$B_sum_TU20, data = B, fct = LL.2())
B_IA_DRM_TU20 <- drm(B$IA~B$B_sum_TU20,  data = B, fct = LL.2())
B_CA_DRM_TU20 <- drm(B$CA~B$B_sum_TU20,  data = B, fct = LL.2())

T_DRM_TU20   <- drm(T$RGR~T$T_sum_TU20, data = T, fct = LL.2())
T_IA_DRM_TU20 <- drm(T$IA~T$T_sum_TU20, data = T, fct = LL.2())
T_CA_DRM_TU20 <- drm(T$CA~T$T_sum_TU20, data = T, fct = LL.2())

Q_DRM_TU20    <- drm(Q$RGR~Q$Q_sum_TU20, data = Q, fct = LL.2())
Q_IA_DRM_TU20 <- drm(Q$IA~Q$Q_sum_TU20,  data = Q, fct = LL.2())
Q_CA_DRM_TU20 <- drm(Q$CA~Q$Q_sum_TU20,  data = Q, fct = LL.2())

R_DRM_TU20    <- drm(R$RGR~R$R_sum_TU20, data = R, fct = LL.2())
R_IA_DRM_TU20 <- drm(R$IA~R$R_sum_TU20,  data = R, fct = LL.2())
R_CA_DRM_TU20 <- drm(R$CA~R$R_sum_TU20,  data = R, fct = LL.2())

X_DRM_TU20    <- drm(X$RGR~X$X_sum_TU20, data = X, fct = LL.2())
X_IA_DRM_TU20 <- drm(X$IA~X$X_sum_TU20,  data = X, fct = LL.2())
X_CA_DRM_TU20 <- drm(X$CA~X$X_sum_TU20,  data = X, fct = LL.2())


#EC AND SLOPE TABLE --> (EC20)
EC_B_TU20 <- ED(B_DRM_TU20, c(10,20,50), interval="delta")
EC_T_TU20 <- ED(T_DRM_TU20, c(10,20,50), interval="delta")
EC_Q_TU20 <- ED(Q_DRM_TU20, c(10,20,50), interval="delta")
EC_R_TU20 <- ED(R_DRM_TU20, c(10,20,50), interval="delta")
EC_X_TU20 <- ED(X_DRM_TU20, c(10,20,50), interval="delta")


EC10_B_TU20 <- EC_B_TU20[1:1]
SE10_B_TU20 <- EC_B_TU20["e:1:10" , "Std. Error"]
EC20_B_TU20 <- EC_B_TU20["e:1:20" , "Estimate"]
SE20_B_TU20 <- EC_B_TU20["e:1:20" , "Std. Error"]
EC50_B_TU20 <- EC_B_TU20["e:1:50" , "Estimate"]
SE50_B_TU20 <- EC_B_TU20["e:1:50" , "Std. Error"]
U10_B_TU20 <-  EC_B_TU20["e:1:10" , "Upper"]
L10_B_TU20 <-  EC_B_TU20["e:1:10" , "Lower"]
U20_B_TU20 <-  EC_B_TU20["e:1:20" , "Upper"]
L20_B_TU20 <-  EC_B_TU20["e:1:20" , "Lower"]
U50_B_TU20 <-  EC_B_TU20["e:1:50" , "Upper"]
L50_B_TU20 <-  EC_B_TU20["e:1:50" , "Lower"]


EC10_T_TU20 <- EC_T_TU20[1:1]
SE10_T_TU20 <- EC_T_TU20["e:1:10" , "Std. Error"]
EC20_T_TU20 <- EC_T_TU20["e:1:20" , "Estimate"]
SE20_T_TU20 <- EC_T_TU20["e:1:20" , "Std. Error"]
EC50_T_TU20 <- EC_T_TU20["e:1:50" , "Estimate"]
SE50_T_TU20 <- EC_T_TU20["e:1:50" , "Std. Error"]
U10_T_TU20 <-  EC_T_TU20["e:1:10" , "Upper"]
L10_T_TU20 <-  EC_T_TU20["e:1:10" , "Lower"]
U20_T_TU20 <-  EC_T_TU20["e:1:20" , "Upper"]
L20_T_TU20 <-  EC_T_TU20["e:1:20" , "Lower"]
U50_T_TU20 <-  EC_T_TU20["e:1:50" , "Upper"]
L50_T_TU20 <-  EC_T_TU20["e:1:50" , "Lower"]


EC10_Q_TU20 <- EC_Q_TU20[1:1]
SE10_Q_TU20 <- EC_Q_TU20["e:1:10" , "Std. Error"]
EC20_Q_TU20 <- EC_Q_TU20["e:1:20" , "Estimate"]
SE20_Q_TU20 <- EC_Q_TU20["e:1:20" , "Std. Error"]
EC50_Q_TU20 <- EC_Q_TU20["e:1:50" , "Estimate"]
SE50_Q_TU20 <- EC_Q_TU20["e:1:50" , "Std. Error"]
U10_Q_TU20 <-  EC_Q_TU20["e:1:10" , "Upper"]
L10_Q_TU20 <-  EC_Q_TU20["e:1:10" , "Lower"]
U20_Q_TU20 <-  EC_Q_TU20["e:1:20" , "Upper"]
L20_Q_TU20 <-  EC_Q_TU20["e:1:20" , "Lower"]
U50_Q_TU20 <-  EC_Q_TU20["e:1:50" , "Upper"]
L50_Q_TU20 <-  EC_Q_TU20["e:1:50" , "Lower"]


EC10_R_TU20 <- EC_R_TU20[1:1]
SE10_R_TU20 <- EC_R_TU20["e:1:10" , "Std. Error"]
EC20_R_TU20 <- EC_R_TU20["e:1:20" , "Estimate"]
SE20_R_TU20 <- EC_R_TU20["e:1:20" , "Std. Error"]
EC50_R_TU20 <- EC_R_TU20["e:1:50" , "Estimate"]
SE50_R_TU20 <- EC_R_TU20["e:1:50" , "Std. Error"]
U10_R_TU20 <-  EC_R_TU20["e:1:10" , "Upper"]
L10_R_TU20 <-  EC_R_TU20["e:1:10" , "Lower"]
U20_R_TU20 <-  EC_R_TU20["e:1:20" , "Upper"]
L20_R_TU20 <-  EC_R_TU20["e:1:20" , "Lower"]
U50_R_TU20 <-  EC_R_TU20["e:1:50" , "Upper"]
L50_R_TU20 <-  EC_R_TU20["e:1:50" , "Lower"]


EC10_X_TU20 <- EC_X_TU20[1:1]
SE10_X_TU20 <- EC_X_TU20["e:1:10" , "Std. Error"]
EC20_X_TU20 <- EC_X_TU20["e:1:20" , "Estimate"]
SE20_X_TU20 <- EC_X_TU20["e:1:20" , "Std. Error"]
EC50_X_TU20 <- EC_X_TU20["e:1:50" , "Estimate"]
SE50_X_TU20 <- EC_X_TU20["e:1:50" , "Std. Error"]
U10_X_TU20 <-  EC_X_TU20["e:1:10" , "Upper"]
L10_X_TU20 <-  EC_X_TU20["e:1:10" , "Lower"]
U20_X_TU20 <-  EC_X_TU20["e:1:20" , "Upper"]
L20_X_TU20 <-  EC_X_TU20["e:1:20" , "Lower"]
U50_X_TU20 <-  EC_X_TU20["e:1:50" , "Upper"]
L50_X_TU20 <-  EC_X_TU20["e:1:50" , "Lower"]


B_DRM_df_TU20 <-  summary(B_DRM_TU20)
B_DRM_df_TU20 <- B_DRM_df_TU20$coefficients
B_DRM_df_TU20 <- as.matrix(B_DRM_df_TU20, header=TRUE)
B_B_TU20   <- B_DRM_df_TU20[1:1]
B_SEB_TU20 <- B_DRM_df_TU20["b:(Intercept)", "Std. Error"]

T_DRM_df_TU20 <-  summary(T_DRM_TU20)
T_DRM_df_TU20 <- T_DRM_df_TU20$coefficients
T_DRM_df_TU20 <- as.matrix(T_DRM_df_TU20, header=TRUE)
T_B_TU20   <- T_DRM_df_TU20[1:1]
T_SEB_TU20 <- T_DRM_df_TU20["b:(Intercept)", "Std. Error"]

Q_DRM_df_TU20 <-  summary(Q_DRM_TU20)
Q_DRM_df_TU20 <- Q_DRM_df_TU20$coefficients
Q_DRM_df_TU20 <- as.matrix(Q_DRM_df_TU20, header=TRUE)
Q_B_TU20   <- Q_DRM_df_TU20[1:1]
Q_SEB_TU20 <- Q_DRM_df_TU20["b:(Intercept)", "Std. Error"]

R_DRM_df_TU20 <-  summary(R_DRM_TU20)
R_DRM_df_TU20 <- R_DRM_df_TU20$coefficients
R_DRM_df_TU20 <- as.matrix(R_DRM_df_TU20, header=TRUE)
R_B_TU20   <- R_DRM_df_TU20[1:1]
R_SEB_TU20 <- R_DRM_df_TU20["b:(Intercept)", "Std. Error"]

X_DRM_df_TU20 <-  summary(X_DRM_TU20)
X_DRM_df_TU20 <- X_DRM_df_TU20$coefficients
X_DRM_df_TU20 <- as.matrix(X_DRM_df_TU20, header=TRUE)
X_B_TU20   <- X_DRM_df_TU20[1:1]
X_SEB_TU20 <- X_DRM_df_TU20["b:(Intercept)", "Std. Error"]


Results_summary_TU20 <- data.frame(B = numeric(14), T = numeric(14), Q = numeric(14), R = numeric(14), X = numeric(14),
                                   row.names = c("EC10", "EC20", "EC50", "slope", "SE_slope", "SE_EC10", "SE_EC20","SE_EC50", 
                                                 "Upper_EC10", "Lower_EC10", "Upper_EC20", "Lower_EC20", "Upper_EC50", "Lower_EC50"))

Results_summary_TU20[1, "B"]  <- as.numeric(EC10_B_TU20)
Results_summary_TU20[1, "T"]  <- as.numeric(EC10_T_TU20)
Results_summary_TU20[1, "Q"]  <- as.numeric(EC10_Q_TU20)
Results_summary_TU20[1, "R"]  <- as.numeric(EC10_R_TU20)
Results_summary_TU20[1, "X"]  <- as.numeric(EC10_X_TU20)

Results_summary_TU20[2, "B"]  <- as.numeric(EC20_B_TU20)
Results_summary_TU20[2, "T"]  <- as.numeric(EC20_T_TU20)
Results_summary_TU20[2, "Q"]  <- as.numeric(EC20_Q_TU20)
Results_summary_TU20[2, "R"]  <- as.numeric(EC20_R_TU20)
Results_summary_TU20[2, "X"]  <- as.numeric(EC20_X_TU20)

Results_summary_TU20[3, "B"]  <- as.numeric(EC50_B_TU20)
Results_summary_TU20[3, "T"]  <- as.numeric(EC50_T_TU20)
Results_summary_TU20[3, "Q"]  <- as.numeric(EC50_Q_TU20)
Results_summary_TU20[3, "R"]  <- as.numeric(EC50_R_TU20)
Results_summary_TU20[3, "X"]  <- as.numeric(EC50_X_TU20)

Results_summary_TU20[4, "B"]  <- as.numeric(B_B_TU20)
Results_summary_TU20[4, "T"]  <- as.numeric(T_B_TU20)
Results_summary_TU20[4, "Q"]  <- as.numeric(Q_B_TU20)
Results_summary_TU20[4, "R"]  <- as.numeric(R_B_TU20)
Results_summary_TU20[4, "X"]  <- as.numeric(X_B_TU20)

Results_summary_TU20[5, "B"]  <- as.numeric(B_SEB_TU20)
Results_summary_TU20[5, "T"]  <- as.numeric(T_SEB_TU20)
Results_summary_TU20[5, "Q"]  <- as.numeric(Q_SEB_TU20)
Results_summary_TU20[5, "R"]  <- as.numeric(R_SEB_TU20)
Results_summary_TU20[5, "X"]  <- as.numeric(X_SEB_TU20)

Results_summary_TU20[6, "B"]  <- as.numeric(SE10_B_TU20)
Results_summary_TU20[6, "T"]  <- as.numeric(SE10_T_TU20)
Results_summary_TU20[6, "Q"]  <- as.numeric(SE10_Q_TU20)
Results_summary_TU20[6, "R"]  <- as.numeric(SE10_R_TU20)
Results_summary_TU20[6, "X"]  <- as.numeric(SE10_X_TU20)

Results_summary_TU20[7, "B"]  <- as.numeric(SE20_B_TU20)
Results_summary_TU20[7, "T"]  <- as.numeric(SE20_T_TU20)
Results_summary_TU20[7, "Q"]  <- as.numeric(SE20_Q_TU20)
Results_summary_TU20[7, "R"]  <- as.numeric(SE20_R_TU20)
Results_summary_TU20[7, "X"]  <- as.numeric(SE20_X_TU20)

Results_summary_TU20[8, "B"]  <- as.numeric(SE50_B_TU20)
Results_summary_TU20[8, "T"]  <- as.numeric(SE50_T_TU20)
Results_summary_TU20[8, "Q"]  <- as.numeric(SE50_Q_TU20)
Results_summary_TU20[8, "R"]  <- as.numeric(SE50_R_TU20)
Results_summary_TU20[8, "X"]  <- as.numeric(SE50_X_TU20)

Results_summary_TU20[9, "B"]  <- as.numeric(U10_B_TU20)
Results_summary_TU20[9, "T"]  <- as.numeric(U10_T_TU20)
Results_summary_TU20[9, "Q"]  <- as.numeric(U10_Q_TU20)
Results_summary_TU20[9, "R"]  <- as.numeric(U10_R_TU20)
Results_summary_TU20[9, "X"]  <- as.numeric(U10_X_TU20)

Results_summary_TU20[10, "B"]  <- as.numeric(L10_B_TU20)
Results_summary_TU20[10, "T"]  <- as.numeric(L10_T_TU20)
Results_summary_TU20[10, "Q"]  <- as.numeric(L10_Q_TU20)
Results_summary_TU20[10, "R"]  <- as.numeric(L10_R_TU20)
Results_summary_TU20[10, "X"]  <- as.numeric(L10_X_TU20)

Results_summary_TU20[11, "B"]  <- as.numeric(U20_B_TU20)
Results_summary_TU20[11, "T"]  <- as.numeric(U20_T_TU20)
Results_summary_TU20[11, "Q"]  <- as.numeric(U20_Q_TU20)
Results_summary_TU20[11, "R"]  <- as.numeric(U20_R_TU20)
Results_summary_TU20[11, "X"]  <- as.numeric(U20_X_TU20)

Results_summary_TU20[12, "B"]  <- as.numeric(L20_B_TU20)
Results_summary_TU20[12, "T"]  <- as.numeric(L20_T_TU20)
Results_summary_TU20[12, "Q"]  <- as.numeric(L20_Q_TU20)
Results_summary_TU20[12, "R"]  <- as.numeric(L20_R_TU20)
Results_summary_TU20[12, "X"]  <- as.numeric(L20_X_TU20)

Results_summary_TU20[13, "B"]  <- as.numeric(U50_B_TU20)
Results_summary_TU20[13, "T"]  <- as.numeric(U50_T_TU20)
Results_summary_TU20[13, "Q"]  <- as.numeric(U50_Q_TU20)
Results_summary_TU20[13, "R"]  <- as.numeric(U50_R_TU20)
Results_summary_TU20[13, "X"]  <- as.numeric(U50_X_TU20)

Results_summary_TU20[14, "B"]  <- as.numeric(L50_B_TU20)
Results_summary_TU20[14, "T"]  <- as.numeric(L50_T_TU20)
Results_summary_TU20[14, "Q"]  <- as.numeric(L50_Q_TU20)
Results_summary_TU20[14, "R"]  <- as.numeric(L50_R_TU20)
Results_summary_TU20[14, "X"]  <- as.numeric(L50_X_TU20)


write.table(Results_summary_TU20, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Summaries//Results_summary_TU20.txt", 
            sep = "\t", quote = FALSE, row.names = TRUE)

#################################################################################
# MIF (EC10)
#################################################################################

MIF_df <- data.frame(n = numeric(5), MIF = numeric(5), MIF_SE = numeric(5), group = numeric(5),
                     row.names = c("B", "T", "Q", "X", "R"))

MIF_df["B", 1] <- 2
MIF_df["T", 1] <- 3
MIF_df["Q", 1] <- 4
MIF_df["X", 1] <- 4.9999999
MIF_df["R", 1] <- 5.0000001

MIF_df["B", 2] <- as.numeric(EC10_B_TU10)
MIF_df["T", 2] <- as.numeric(EC10_T_TU10)
MIF_df["Q", 2] <- as.numeric(EC10_Q_TU10)
MIF_df["X", 2] <- as.numeric(EC10_X_TU10)
MIF_df["R", 2] <- as.numeric(EC10_R_TU10)

MIF_df["B", 3] <- as.numeric(SE10_B_TU10)
MIF_df["T", 3] <- as.numeric(SE10_T_TU10)
MIF_df["Q", 3] <- as.numeric(SE10_Q_TU10)
MIF_df["X", 3] <- as.numeric(SE10_X_TU10)
MIF_df["R", 3] <- as.numeric(SE10_R_TU10)

write.table(MIF_df, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Summaries//MIF.txt", 
            sep = "\t", quote = FALSE, row.names = TRUE)


MIF_df$group <- c("Equitoxic ray", 
                  "Equitoxic ray",
                  "Equitoxic ray",
                  "Equitoxic ray",
                  "Environmentally relevant ray")

MIF_df$RGR <- c(0.9,0.9,0.9,0.9,0.9)    
MIF_df$name <- c("B", "T", "Q", "X", "R")

#################################################################################
# MIF (EC20)
#################################################################################

MIF20_df <- data.frame(n = numeric(5), MIF20 = numeric(5), MIF20_SE = numeric(5), group = numeric(5),
                       row.names = c("B", "T", "Q", "X", "R"))

MIF20_df["B", 1] <- 2
MIF20_df["T", 1] <- 3
MIF20_df["Q", 1] <- 4
MIF20_df["X", 1] <- 4.9999999
MIF20_df["R", 1] <- 5.0000001

MIF20_df["B", 2] <- as.numeric(EC20_B_TU20)
MIF20_df["T", 2] <- as.numeric(EC20_T_TU20)
MIF20_df["Q", 2] <- as.numeric(EC20_Q_TU20)
MIF20_df["X", 2] <- as.numeric(EC20_X_TU20)
MIF20_df["R", 2] <- as.numeric(EC20_R_TU20)

MIF20_df["B", 3] <- as.numeric(SE20_B_TU20)
MIF20_df["T", 3] <- as.numeric(SE20_T_TU20)
MIF20_df["Q", 3] <- as.numeric(SE20_Q_TU20)
MIF20_df["X", 3] <- as.numeric(SE20_X_TU20)
MIF20_df["R", 3] <- as.numeric(SE20_R_TU20)

write.table(MIF20_df, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Summaries//MIF20.txt", 
            sep = "\t", quote = FALSE, row.names = TRUE)


MIF20_df$group <- c("Equitoxic ray", 
                    "Equitoxic ray",
                    "Equitoxic ray",
                    "Equitoxic ray",
                    "Environmentally relevant ray")

MIF20_df$RGR <- c(0.8,0.8,0.8,0.8,0.8)    
MIF20_df$name <- c("B", "T", "Q", "X", "R")


#################################################################################
# MIF_IA (EC10)
#################################################################################

MIF_IA_df <- data.frame(n = numeric(5), MIF_IA = numeric(5), MIF_IA_SE = numeric(5), group = numeric(5),
                        row.names = c("B", "T", "Q", "X", "R"))

IA_STU10_B <- ED(B_IA_DRM_TU10, 10) 
IA_STU10_T <- ED(T_IA_DRM_TU10, 10)
IA_STU10_Q <- ED(Q_IA_DRM_TU10, 10) 
IA_STU10_X <- ED(X_IA_DRM_TU10, 10)
IA_STU10_R <- ED(R_IA_DRM_TU10, 10) 

MIF_IA_df["B", 1] <- 2
MIF_IA_df["T", 1] <- 3
MIF_IA_df["Q", 1] <- 4
MIF_IA_df["X", 1] <- 4.9999999
MIF_IA_df["R", 1] <- 5.0000001

MIF_IA_df["B", 2] <- EC10_B_TU10/IA_STU10_B[1,1]
MIF_IA_df["T", 2] <- EC10_T_TU10/IA_STU10_T[1,1]
MIF_IA_df["Q", 2] <- EC10_Q_TU10/IA_STU10_Q[1,1]
MIF_IA_df["X", 2] <- EC10_X_TU10/IA_STU10_X[1,1]
MIF_IA_df["R", 2] <- EC10_R_TU10/IA_STU10_R[1,1]

#calculation of the standard error of a division 
MIF_IA_df["B", 3] <- (EC10_B_TU10/IA_STU10_B[1,1])*sqrt(((EC_B_TU10[1,2]/EC10_B_TU10)^2)+((IA_STU10_B[1,2]/IA_STU10_B[1,1])^2))
MIF_IA_df["T", 3] <- (EC10_T_TU10/IA_STU10_T[1,1])*sqrt(((EC_T_TU10[1,2]/EC10_T_TU10)^2)+((IA_STU10_T[1,2]/IA_STU10_T[1,1])^2))
MIF_IA_df["Q", 3] <- (EC10_Q_TU10/IA_STU10_Q[1,1])*sqrt(((EC_Q_TU10[1,2]/EC10_Q_TU10)^2)+((IA_STU10_Q[1,2]/IA_STU10_Q[1,1])^2))
MIF_IA_df["X", 3] <- (EC10_X_TU10/IA_STU10_X[1,1])*sqrt(((EC_X_TU10[1,2]/EC10_X_TU10)^2)+((IA_STU10_X[1,2]/IA_STU10_X[1,1])^2))
MIF_IA_df["R", 3] <- (EC10_R_TU10/IA_STU10_R[1,1])*sqrt(((EC_R_TU10[1,2]/EC10_R_TU10)^2)+((IA_STU10_R[1,2]/IA_STU10_R[1,1])^2))

write.table(MIF_IA_df, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Summaries//MIF_IA.txt", 
            sep = "\t", quote = FALSE, row.names = TRUE)


MIF_IA_df$group <- c("Equitoxic ray", 
                     "Equitoxic ray",
                     "Equitoxic ray",
                     "Equitoxic ray",
                     "Environmentally relevant ray")

MIF_IA_df$RGR <- c(0.9,0.9,0.9,0.9,0.9)    
MIF_IA_df$name <- c("B", "T", "Q", "X", "R")


########################################################################################################
#CHECK IF CA IS CORRECT
########################################################################################################

check_CA_STU10_B <- ED(B_CA_DRM_TU10, 10) 
check_CA_STU10_T <- ED(T_CA_DRM_TU10, 10)
check_CA_STU10_Q <- ED(Q_CA_DRM_TU10, 10) 
check_CA_STU10_R <- ED(R_CA_DRM_TU10, 10)
check_CA_STU10_X <- ED(X_CA_DRM_TU10, 10) 

check_CA_STU20_B <- ED(B_CA_DRM_TU20, 20) 
check_CA_STU20_T <- ED(T_CA_DRM_TU20, 20)
check_CA_STU20_Q <- ED(Q_CA_DRM_TU20, 20)
check_CA_STU20_R <- ED(R_CA_DRM_TU20, 20)
check_CA_STU20_X <- ED(X_CA_DRM_TU20, 20)

check_CA_STU50_B <- ED(B_CA_DRM, 50) 
check_CA_STU50_T <- ED(T_CA_DRM, 50)
check_CA_STU50_Q <- ED(Q_CA_DRM, 50)
check_CA_STU50_R <- ED(R_CA_DRM, 50)
check_CA_STU50_X <- ED(X_CA_DRM, 50)

check_CA <- data.frame(CA_10 = numeric(5), CA_20 = numeric(5), CA_50 = numeric(5),
                       row.names = c("B", "T", "Q", "R", "X"))

check_CA[1, 1]  <- round(as.numeric(check_CA_STU10_B[1, 1]), 2)
check_CA[2, 1]  <- round(as.numeric(check_CA_STU10_T[1, 1]), 2)
check_CA[3, 1]  <- round(as.numeric(check_CA_STU10_Q[1, 1]), 2)
check_CA[4, 1]  <- round(as.numeric(check_CA_STU10_R[1, 1]), 2)
check_CA[5, 1]  <- round(as.numeric(check_CA_STU10_X[1, 1]), 2)

check_CA[1, 2]  <- round(as.numeric(check_CA_STU20_B[1, 1]), 2)
check_CA[2, 2]  <- round(as.numeric(check_CA_STU20_T[1, 1]), 2)
check_CA[3, 2]  <- round(as.numeric(check_CA_STU20_Q[1, 1]), 2)
check_CA[4, 2]  <- round(as.numeric(check_CA_STU20_R[1, 1]), 2)
check_CA[5, 2]  <- round(as.numeric(check_CA_STU20_X[1, 1]), 2)

check_CA[1, 3]  <- round(as.numeric(check_CA_STU50_B[1, 1]), 2)
check_CA[2, 3]  <- round(as.numeric(check_CA_STU50_T[1, 1]), 2)
check_CA[3, 3]  <- round(as.numeric(check_CA_STU50_Q[1, 1]), 2)
check_CA[4, 3]  <- round(as.numeric(check_CA_STU50_R[1, 1]), 2)
check_CA[5, 3]  <- round(as.numeric(check_CA_STU50_X[1, 1]), 2)

write.table(check_CA, file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Summaries//check_CA.txt", 
            sep = "\t", quote = FALSE, row.names = TRUE)


#########################################################
# Dose-Response Models and plots per Metal
#########################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//all_SE")

plot(Me1_DRM, broken = TRUE, type = "none", ylab = "Relative growth rate", cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3,
     xlab = "Concentration [?g/L]",  main = "Single exposure curves", lty=1, lwd=2,
     ylim = c(0,1), xlim=c(0.01,100000), col = "purple", bg= "purple") 

plot(Me2_DRM, broken = TRUE, type = "none", ylab = "Relative growth rate", lty=1, lwd=2,
     xlab = "Concentration [?g/L]", main = "Single exposure curves", 
     ylim = c(0,1), col = "#d11141", bg= "#d11141", add = TRUE) 

plot(Me3_DRM, broken = TRUE, type = "none", ylab = "Relative growth rate", lty=1, lwd=2,
     xlab = "Concentration [?g/L]", main = "Single exposure curves", 
     ylim = c(0,1), col = "#999999", bg=  "#999999", add = TRUE) 

plot(Me4_DRM, broken = TRUE, type = "none", ylab = "Relative growth rate", lty=1, lwd=2,
     xlab = "Concentration [?g/L]", main = "Single exposure curves", 
     ylim = c(0,1), col = "#00ccff", bg= "#00ccff", add = TRUE) #col="#FFFF00"

plot(Me5_DRM, broken = TRUE, type = "none", ylab = "Relative growth rate", lty=1, lwd=2,
     xlab = "Concentration [?g/L]", main = "Single exposure curves", 
     ylim = c(0,1), col = "black", bg= "black", add = TRUE) 

#if you want without point put type = "none", and remove pch = 1, lty=2, lwd=2, and bg= "color"

legend(0.8, 0.27, cex=1.2, legend=c("As", "Zn", "Pb", "Ag", "Cu"), 
       col = c("purple","#d11141", "#999999", "#00ccff", "black"),
       lty=c(1, 1, 1, 1, 1), lwd=c(2, 2, 2, 2, 2), box.col = "transparent")

dev.off()

################################
#       Arsenic          
################################

pdf(file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//Arsenic")

plot(Me1_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Concentration [?g/L]", main = "Arsenic", pch = 5,
     cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3, lwd=2.5,
     ylim = c(0,1.05), xlim=c(0,10000), col = "purple", bg="purple")
plot(Me1_DRM, type = "confidence", col = "purple", add = TRUE)


legend("bottomleft", cex=1.2, legend=c("Arsenic"), 
       col = c("purple"), pt.bg = c("purple"),
       lwd=2.5, pch=5, box.col = "transparent")

dev.off()

################
# 
# pdf(file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//Arsenic_OECD")
# 
# plot(Me1_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Concentration [?g/L]", main = "Arsenic", pch = 5, ylim = c(0,1), add = TRUE)
# plot(Me1_DRM, type = "confidence", col = "grey", add = TRUE)
# 
# plot(Me1_DRM_OECD, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Concentration [?g/L]", main = "Arsenic", pch = 5, ylim = c(0,1), add = TRUE)
# plot(Me1_DRM_OECD, type = "confidence", add = TRUE, col = "blue")
# legend("bottomleft", cex=1, legend=c("Arsenic", "Arsenic OECD"), fill = c("grey","blue"))
# 
# dev.off()


################################
#           Zinc        
################################

pdf(file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//Zinc")

plot(Me2_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Concentration [?g/L]", main = "Zinc", pch = 4, 
     cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3, lwd=2.5,
     ylim = c(0,1.05), xlim=c(0,300), col = "black", bg="black")
plot(Me2_DRM, type = "confidence", col = "black", add = TRUE)

legend("bottomleft", cex=1.2, legend=c("Zinc"), 
       col = c("black"), pt.bg = c("black"),
       lwd=2.5 ,pch=4, box.col = "transparent")

dev.off()


###############

# pdf(file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//Zinc_OECD")
# 
# plot(Cu_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Concentration [?g/L]", main = "Copper", pch = 5,ylim = c(0,1), add = TRUE)
# plot(Cu_DRM, type = "confidence", col = "grey", add = TRUE)
# 
# plot(Cu_DRM_OECD, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Concentration [?g/L]", main = "Copper OECD", pch = 5, ylim = c(0,1), add = TRUE)
# plot(Cu_OECD, type = "confidence", col = "blue", add = TRUE)
# legend("bottomleft", cex=1, legend=c("Copper", "Copper OECD"), fill = c("grey","blue"))
# 
# dev.off()

################################
#           Lead            
################################

pdf(file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//Lead")

plot(Me3_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Concentration [?g/L]", main = "Lead", pch = 19, 
     cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3, lwd=2.5,
     ylim = c(0,1.05), xlim=c(0,300), col = "#d11141", bg="#d11141")
plot(Me3_DRM, type = "confidence", col = "#d11141", add = TRUE)

legend("bottomleft", cex=1.2, legend=c("Lead"), 
       col = c("#d11141"), pt.bg = c("#d11141"),
       lwd=2.5, pch=19, box.col = "transparent")

dev.off()

########################

#pdf(file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//Lead_OECD")
# 
# plot(Cu_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Concentration [?g/L]", main = "Copper", pch = 5, add = TRUE)
# plot(Cu_DRM, type = "confidence", col = "grey", add = TRUE)
# 
# plot(Cu_DRM_OECD, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Concentration [?g/L]", main = "Copper OECD", pch = 5, add = TRUE)
# plot(Cu_DRM_OECD, type = "confidence", col = "blue", add = TRUE)
# legend("bottomleft", cex=1, legend=c("Copper", "Copper OECD"), fill = c("grey","blue"))
# 
# dev.off()
# 
################################
#           Silver             
################################

pdf(file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//Silver")

plot(Me4_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Concentration [?g/L]", main = "Silver", pch = 23,
     cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3, lwd=2.5,
     ylim = c(0,1.05), xlim=c(0,10), col = "#00ccff", bg="#00ccff")
plot(Me4_DRM, type = "confidence", col = "#00ccff", add = TRUE)

legend("bottomleft", cex=1.2, legend=c("Silver"), 
       col = c("#00ccff"), pt.bg = c("#00ccff"),
       lwd=2.5, pch=23, box.col = "transparent")

dev.off()


############################

#pdf(file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//Silver_OECD")
# 
# plot(Ni_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Concentration [?g/L]", main = "Nickel", pch = 5, add = TRUE)
# plot(Ni_DRM, type = "confidence", col = "grey", add = TRUE)
# 
# plot(Ni_DRM_OECD, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Concentration [?g/L]", main = "Nickel", pch = 5, add = TRUE)
# plot(Ni_DRM_OECD, type = "confidence", col = "blue", add = TRUE)
# legend("bottomleft", cex=1, legend=c("Nickel", "Nickel OECD"), fill = c("grey","blue"))
# 
# dev.off()

################################
#           Copper            
################################

pdf(file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//Copper")

plot(Me5_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Concentration [?g/L]", main = "Copper", pch = 22,
     cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3, lwd=2.5,
     ylim = c(0,1.05), xlim=c(0,200), col = "#999999", bg="#999999")
plot(Me5_DRM, type = "confidence", col = "#999999", add = TRUE)

legend("bottomleft", cex=1.2, legend=c("Copper"), 
       col = c("#999999"), pt.bg = c("#999999"),
       lwd=2.5, pch=22, box.col = "transparent")

dev.off()


###############
# 
# pdf(file = "C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//SE//Copper_OECD")
# 
# plot(Pb_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Concentration [?g/L]", main = "Lead", pch = 5, add = TRUE)
# plot(Pb_DRM, type = "confidence", col = "grey", add = TRUE)
# 
# plot(Pb_DRM_OECD, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Concentration [?g/L]", main = "Lead", pch = 5, add = TRUE)
# plot(Pb_DRM_OECD, type = "confidence", col = "blue", add = TRUE)
# legend("bottomleft", cex=1, legend=c("Lead", "Lead OECD"), fill = c("grey","blue"))
# 
# dev.off()


################################
#           Binary            
################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC50_graphs//B vs TU (EC50)")

plot(Me1_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units",  main = "Binary", cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3, 
     xlim=c(0.01,20), col = "purple", lwd=2.5, ylim = c(0,1)) 

plot(Me2_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Binary", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#d11141", lwd=2.5, add = TRUE) 

plot(B_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Binary", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.01,20), col = "#00cc00", bg="#00cc00", add = TRUE) 

plot(B_IA_DRM, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Binary", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.01,20), col = "orange", bg="orange", add = TRUE) 

plot(B_CA_DRM, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units", main = "Binary", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.01,20), col = "blue", bg="blue", add = TRUE)

legend(0.01, 0.27, cex=1.2, legend=c("As", "Zn", "Observed As-Zn", "IA prediction", "CA prediction"), 
       col = c("purple","#d11141", "#00cc00", "orange", "blue"), pt.bg = c("purple","#d11141", "#00cc00", "orange", "blue"),
       lty=c(1, 1, 2, NA, NA),lwd=c(2.5, 2.5, 2, NA, NA),pch=c(NA,NA, 1, 22, 23), box.col = "transparent")

dev.off()

################################
#           Ternary            
################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC50_graphs//T vs TU (EC50)")

plot(Me1_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units",  main = "Ternary", cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.01,20), col = "purple", lwd=2.5) 

plot(Me2_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Ternary",  
     ylim = c(0,1), xlim=c(0.01,20), col = "#d11141", lwd=2.5, add = TRUE) 

plot(Me3_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Ternary", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#999999", lwd=2.5, add = TRUE) 

plot(T_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Ternary", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.01,20), col = "#00cc00", bg="#00cc00", add = TRUE) 

plot(T_IA_DRM, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Ternary", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.01,20), col = "orange", bg="orange", add = TRUE)

plot(T_CA_DRM, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units", main = "Ternary", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.01,20), col = "blue", bg="blue", add = TRUE)

legend(0.01, 0.32, cex=1.2, legend=c("Me2", "Me2", "Me3", "Observed As-Me2-Me3", "IA prediction", "CA prediction"), 
       col = c("purple","#d11141", "#999999", "#00cc00", "orange", "blue"), pt.bg = c("purple","#d11141", "#999999", "#00cc00", "orange", "blue"),
       lty=c(1, 1, 1, 2, NA, NA),lwd=c(2.5, 2.5, 2.5, 2, NA, NA),pch=c(NA,NA, NA, 1, 22, 23), box.col = "transparent")

dev.off()

################################
#           Quaternary            
################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC50_graphs//Q vs TU (EC50)")

plot(Me1_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units",  main = "Quaternary", cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3, 
     ylim = c(0,1), xlim=c(0.01,20), col = "purple", lwd=2.5) 

plot(Me2_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quaternary", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#d11141", lwd=2.5, add = TRUE) 

plot(Me3_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quaternary", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#999999", lwd=2.5, add = TRUE) 

plot(Me4_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quaternary", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#00ccff",lwd=2.5, add = TRUE) 

plot(Q_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quaternary", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.01,20), col = "#00cc00", bg="#00cc00", add = TRUE) 

plot(Q_IA_DRM, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quaternary", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.01,20), col = "orange", bg="orange", add = TRUE) 

plot(Q_CA_DRM, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units", main = "Quaternary", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.01,20), col = "blue", bg="blue", add = TRUE)

legend(0.01, 0.37, cex=1.2, legend=c("Me2", "Me2", "Me3", "Me4", "Observed Me1-Me2-Me3-Ni", "IA prediction", "CA prediction"), 
       col = c("purple","#d11141", "#999999", "#00ccff", "#00cc00", "orange", "blue"), pt.bg = c("purple","#d11141", "g#999999rey", "#00ccff", "#00cc00", "orange", "blue"),
       lty=c(1, 1, 1, 1, 2, NA, NA),lwd=c(2.5, 2.5, 2.5, 2.5, 2, NA, NA),pch=c(NA,NA, NA,NA, 1, 22, 23), box.col = "transparent")

dev.off()


########################################
#           Quinary (equitoxic ray)            
########################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC50_graphs//X vs TU (EC50)")

plot(Me1_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units",  main = "Quinary (equitoxic ray)", cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.01,20), col = "purple", lwd=2.5) 

plot(Me2_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (equitoxic ray)", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#d11141", lwd=2.5, add = TRUE)

plot(Me3_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (equitoxic ray)", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#999999", lwd=2.5, add = TRUE) 

plot(Me4_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (equitoxic ray)", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#00ccff", lwd=2.5, add = TRUE) 

plot(Me5_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (equitoxic ray)", 
     ylim = c(0,1), xlim=c(0.01, 20), col = "black", lwd=2.5, add = TRUE) 

plot(X_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (equitoxic ray)", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.01,20), col = "#00cc00", bg="#00cc00", add = TRUE) 

plot(X_IA_DRM, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (equitoxic ray)", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.01,20), col = "orange", bg="orange", add = TRUE) 

plot(X_CA_DRM, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units", main = "Quinary (equitoxic ray)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.01,20), col = "blue", bg="blue", add = TRUE)

legend(0.01, 0.42, cex=1.2, legend=c("Me2", "Me2", "Me3", "Me4", "Me5", "Observed Me1-Me2-Me3-Me4-Me5", "IA prediction", "CA prediction"), 
       col = c("purple","#d11141", "#999999", "#00ccff", "black", "#00cc00", "orange", "blue"), pt.bg = c("purple","#d11141", "#999999", "#00ccff", "black", "#00cc00", "orange", "blue"),
       lty=c(1, 1, 1, 1, 1, 2, NA, NA),lwd=c(2.5, 2.5, 2.5, 2.5, 2.5, 2, NA, NA),pch=c(NA, NA, NA, NA, NA, 1, 22, 23), box.col = "transparent")

dev.off()

################################################################
#           Quinary (environmentally relevant ray)            
################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC50_graphs//R vs TU (EC50)")

plot(Me1_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units",  main = "Quinary (environmentally relevant ray)", cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.01,20), col = "purple", lwd=2.5) 

plot(Me2_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (environmentally relevant ray)", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#d11141", lwd=2.5, add = TRUE) 

plot(Me3_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (environmentally relevant ray)", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#999999", lwd=2.5, add = TRUE) 

plot(Me4_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (environmentally relevant ray)", 
     ylim = c(0,1), xlim=c(0.01,20), col = "#00ccff",  lwd=2.5, add = TRUE) #col="#FFFF00"

plot(Me5_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (environmentally relevant ray)", 
     ylim = c(0,1), xlim=c(0.01, 20), col = "black", lwd=2.5, add = TRUE) 

plot(R_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (environmentally relevant ray)", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.01,20), col = "#33cc33", bg="#33cc33", add = TRUE) 

plot(R_IA_DRM, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units", main = "Quinary (environmentally relevant ray)", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.01,20), col = "orange", bg="orange", add = TRUE) 

plot(R_CA_DRM, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units", main = "Quinary (environmentally relevant ray)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.01,20), col = "blue", bg="blue", add = TRUE)

legend(0.01, 0.42, cex=1.2, legend=c("Me2", "Me2", "Me3", "Me4", "Me5", "Observed Me1-Me2-Me3-Me4-Me5", "IA prediction", "CA prediction"), 
       col = c("purple","#d11141", "#999999", "#00ccff", "black", "#00cc00", "orange", "blue"), pt.bg = c("purple","#d11141", "#999999", "#00ccff", "black", "#00cc00", "orange", "blue"),
       lty=c(1, 1, 1, 1, 1, 2, NA, NA),lwd=c(2.5, 2.5, 2.5, 2.5, 2.5, 2, NA, NA),pch=c(NA, NA, NA, NA, NA, 1, 22, 23), box.col = "transparent")

dev.off()

################################################################
#           Quinary (together)            
################################################################

# pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC50_graphs//Quinary together (EC50)")
# 
# plot(Me1_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
#      xlab = "Sum of toxic units",  main = "Quinary", cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3,
#      ylim = c(0,1), xlim=c(0.01,20), col = "#6A1B9A", lwd=2.5) 
# 
# plot(Me2_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
#      xlab = "Sum of toxic units", main = "Quinary", 
#      ylim = c(0,1), xlim=c(0.01,20), col = "#d11141", lwd=2.5, add = TRUE)
# 
# plot(Me3_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
#      xlab = "Sum of toxic units", main = "Quinary", 
#      ylim = c(0,1), xlim=c(0.01,20), col = "#757575", lwd=2.5, add = TRUE) 
# 
# plot(Me4_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
#      xlab = "Sum of toxic units", main = "Quinary", 
#      ylim = c(0,1), xlim=c(0.01,20), col = "#fbc02d", lwd=2.5, add = TRUE) 
# 
# plot(Me5_DRM_TU, broken = TRUE, type = "none", ylab = "Relative growth rate", 
#      xlab = "Sum of toxic units", main = "Quinary", 
#      ylim = c(0,1), xlim=c(0.01, 20), col = "#212121", lwd=2.5, add = TRUE) 
# 
# plot(X_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Sum of toxic units", main = "Quinary", pch = 1, lty=2, lwd=2,
#      ylim = c(0,1), xlim=c(0.01,20), col = "#43A047", bg="#4caf50", add = TRUE) 
# 
# plot(R_DRM, broken = TRUE, type = "all", ylab = "Relative growth rate", 
#      xlab = "Sum of toxic units", main = "Quinary", pch = 1, lty=2, lwd=2,
#      ylim = c(0,1), xlim=c(0.01,20), col = "#D84315", bg="#D84315", add = TRUE) 
# 
# legend(0.01, 0.42, cex=1.2, legend=c("Me2", "Me2", "Me3", "Me4", "Me5", "Quinary (ENR)", "Quinary (EQT)"), #"IA prediction", "CA prediction"
#        col = c("purple","#d11141", "grey", "#40ECD0", "black", "#00cc00", "#D84315"), pt.bg = c("purple","#d11141", "grey", "black", "#40ECD0", "#00cc00", "#D84315"),
#        lty=c(1, 1, 1, 1, 1, 2, 2),lwd=c(2.5, 2.5, 2.5, 2.5, 2.5, 2, 2),pch=c(NA, NA, NA, NA, NA, 1, 1), box.col = "transparent")
# 
# dev.off()
# 
# 
# # col = c("purple","#d11141", "grey", "#40ECD0", "black", "#00cc00", "orange", "blue"), pt.bg = c("purple","#d11141", "grey", "black", "#40ECD0", "#00cc00", "orange", "blue"),
# # lty=c(1, 1, 1, 1, 1, 2, NA, NA),lwd=c(2.5, 2.5, 2.5, 2.5, 2.5, 2, NA, NA),pch=c(NA, NA, NA, NA, NA, 1, 22, 23), box.col = "transparent")


#############################################################################################
#############################################################################################
#############################################################################################

################################
#     Binary (EC10)            
################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC10_graphs//B vs TU (EC10)")

plot(Me1_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "purple", lwd=2.5) 
title("Binary", cex.main=1.6, line = 0.5)

plot(Me2_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)",
     ylim = c(0,1), xlim=c(0.05,100), col = "#d11141", lwd=2.5, add = TRUE)

plot(B_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.05,100), col = "#00cc00", bg="#00cc00", add = TRUE) 

plot(B_IA_DRM_TU10, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.05,100), col = "orange", bg="orange", add = TRUE)

plot(B_CA_DRM_TU10, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "blue", bg="blue", add = TRUE)

points(MIF_df$MIF[MIF_df$name == B][1], MIF_df$RGR[MIF_df$name == B][1],
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.2)

# text(x=70, y=0.9, "CA/IA prediction:\n < observed     antagonism\n> observed    synergism", 
#      cex=1.5, adj=0.5)

# points(IA_STU10_B[1,1], MIF_df$RGR[MIF_df$name == B][1],
#        pch = 4, col = "red", bg="red", lwd=4, cex=1.2)


legend(0.05, 0.32, cex=1.2, legend=c("As", "Zn", "Observed As-Zn", "IA prediction", "CA prediction", "MIF"), 
       col = c("purple","#d11141", "#00cc00", "orange", "blue", "#00ccff"), 
       pt.bg = c("purple","#d11141", "#00cc00", "orange", "blue", "#00ccff"),
       lty=c(1, 1, 2, NA, NA, NA),lwd=c(2.5, 2.5, 2, NA, NA, 4),pch=c(NA,NA, 1, 22, 23, 4), box.col = "transparent")

dev.off()

################################
#           Ternary (EC10)            
################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC10_graphs//T vs TU (EC10)")

plot(Me1_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "purple", lwd=2.5)  #, ylim = c(0,1), xlim=c(0.01,20)
title("Ternary", cex.main=1.6, line = 0.5)

plot(Me2_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#d11141", lwd=2.5, add = TRUE)

plot(Me3_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "grey", lwd=2.5, add = TRUE) 

plot(T_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.05,100), col = "#00cc00", bg="#00cc00", add = TRUE) 

plot(T_IA_DRM_TU10, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.05,100), col = "orange", bg="orange", add = TRUE) 

plot(T_CA_DRM_TU10, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "blue", bg="blue", add = TRUE)

points(MIF_df$MIF[MIF_df$name == T][1], MIF_df$RGR[MIF_df$name == T][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.5)

# points(IA_STU10_T[1,1], MIF_df$RGR[MIF_df$name == T][1],
#        pch = 4, col = "red", bg="red", lwd=4, cex=1.2)

legend(0.05, 0.37, cex=1.2, legend=c("As", "Zn", "Pb", "Observed As-Zn-Pb", "IA prediction", "CA prediction", "MIF"), 
       col = c("purple","#d11141", "grey", "#00cc00", "orange", "blue", "#00ccff"), 
       pt.bg = c("purple","#d11141", "grey", "#00cc00", "orange", "blue", "#00ccff"),
       lty=c(1, 1, 1, 2, NA, NA, NA),lwd=c(2.5, 2.5, 2.5, 2, NA, NA, 4),pch=c(NA,NA, NA, 1, 22, 23, 4), box.col = "transparent")

dev.off()

################################
#           Quaternary (EC10)            
################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC10_graphs//Q vs TU (EC10)")

plot(Me1_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "purple", lwd=2.5) 
title("Quaternary", cex.main=1.6, line = 0.5)

plot(Me2_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#d11141", lwd=2.5, add = TRUE)

plot(Me3_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "grey", lwd=2.5, add = TRUE) 

plot(Me4_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#40ECD0", lwd=2.5, add = TRUE) 

plot(Q_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)",  pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.05,100), col = "#00cc00", bg="#00cc00", add = TRUE) 

plot(Q_IA_DRM_TU10, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)",  pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.1,100), col = "orange", bg="orange", add = TRUE) 

plot(Q_CA_DRM_TU10, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "blue", bg="blue", add = TRUE)

points(MIF_df$MIF[MIF_df$name == Q][1], MIF_df$RGR[MIF_df$name == Q][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.5)

# points(IA_STU10_Q[1,1], MIF_df$RGR[MIF_df$name == Q][1],
#        pch = 4, col = "red", bg="red", lwd=4, cex=1.2)

legend(0.05, 0.42, cex=1.2, legend=c("As", "Zn", "Pb", "Ag", "Observed As-Zn-Pb-Ag", "IA prediction", "CA prediction", "MIF"), 
       col = c("purple","#d11141", "grey", "#40ECD0", "#00cc00", "orange", "blue", "#00ccff"), 
       pt.bg = c("purple","#d11141", "grey", "#40ECD0", "#00cc00", "orange", "blue", "#00ccff"),
       lty=c(1, 1, 1, 1, 2, NA, NA, NA),lwd=c(2.5, 2.5, 2.5, 2.5, 2, NA, NA, 4),pch=c(NA,NA, NA,NA, 1, 22, 23, 4), box.col = "transparent")

dev.off()


########################################
#   Quinary (equitoxic ray) (EC10)            
########################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC10_graphs//X vs TU (EC10)")

plot(Me1_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)",   cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "purple", lwd=2.5) 
title("Quinary (equitoxic ray)", cex.main=1.6, line = 0.5)

plot(Me2_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#d11141", lwd=2.5, add = TRUE)

plot(Me3_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "grey", lwd=2.5, add = TRUE) 

plot(Me4_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#40ECD0", lwd=2.5, add = TRUE) 

plot(Me5_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", 
     ylim = c(0,1), xlim=c(0.05, 100), col = "black", lwd=2.5, add = TRUE) 

plot(X_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.05,100), col = "#00cc00", bg="#00cc00", add = TRUE)  

plot(X_IA_DRM_TU10, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC10)", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.05,100), col = "orange", bg="orange", add = TRUE) 

plot(X_CA_DRM_TU10, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "blue", bg="blue", add = TRUE)

points(MIF_df$MIF[MIF_df$name == X][1], MIF_df$RGR[MIF_df$name == X][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.5)

# points(IA_STU10_X[1,1], MIF_df$RGR[MIF_df$name == X][1],
#        pch = 4, col = "red", bg="red", lwd=4, cex=1.2)

legend(0.05, 0.47, cex=1.2, legend=c("As", "Zn", "Pb", "Ag", "Cu", "Observed As-Zn-Pb-Ag-Cu", "IA prediction", "CA prediction", "MIF"), 
       col = c("purple","#d11141", "grey", "#40ECD0", "black", "#00cc00", "orange", "blue", "#00ccff"), 
       pt.bg = c("purple","#d11141", "grey", "#40ECD0", "black", "#00cc00", "orange", "blue", "#00ccff"),
       lty=c(1, 1, 1, 1, 1, 2, NA, NA, NA),lwd=c(2.5, 2.5, 2.5, 2.5, 2.5, 2, NA, NA, 4),pch=c(NA, NA, NA, NA, NA, 1, 22, 23, 4), box.col = "transparent")

dev.off()

################################################################
#      Quinary (environmentally relevant ray) (EC10)            
################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC10_graphs//R vs TU (EC10)")

plot(Me1_DRM_TU10, broken = TRUE, type = "none",
     ylab = "Relative growth rate", xlab = "Sum of toxic units (EC10)",
     cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "purple", lwd=2.5)
title("Quinary (environmentally relevant ray)", cex.main=1.6, line = 0.5)


plot(Me2_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)",
     ylim = c(0,1), xlim=c(0.05,100), col = "#d11141", lwd=2.5, add = TRUE)

plot(Me3_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)",
     ylim = c(0,1), xlim=c(0.05,100), col = "grey", lwd=2.5, add = TRUE)

plot(Me4_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)",
     ylim = c(0,1), xlim=c(0.05,100), col = "#40ECD0", lwd=2.5, add = TRUE)

plot(Me5_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)",
     ylim = c(0,1), xlim=c(0.05, 100), col = "black", lwd=2.5, add = TRUE)

plot(R_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.05,100), col = "#00cc00", bg="#00cc00", add = TRUE)

plot(R_IA_DRM_TU10, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.05,20), col = "orange", bg="orange", add = TRUE)

plot(R_CA_DRM_TU10, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "blue", bg="blue", add = TRUE)

points(MIF_df$MIF[MIF_df$name == R][1], MIF_df$RGR[MIF_df$name == R][1],
       pch = 7, col = "#f48fb1", bg="#f48fb1", lwd=1.7, cex=1.7)

# points(IA_STU10_R[1,1], MIF_df$RGR[MIF_df$name == R][1],
#        pch = 4, col = "red", bg="red", lwd=4, cex=1.2)

legend(0.05, 0.47, cex=1.2, legend=c("As", "Zn", "Pb", "Ag", "Cu", "Observed As-Zn-Pb-Ag-Cu", "IA prediction", "CA prediction", "MIF"),
       col = c("purple","#d11141", "grey", "#40ECD0", "black", "#00cc00", "orange", "blue", "#f48fb1"),
       pt.bg = c("purple","#d11141", "grey",  "#40ECD0","black", "#00cc00", "orange", "blue", "#f48fb1"),
       lty=c(1, 1, 1, 1, 1, 2, NA, NA, NA),lwd=c(2.5, 2.5, 2.5, 2.5, 2.5, 2, NA, NA, 1.7),pch=c(NA, NA, NA, NA, NA, 1, 22, 23, 7), box.col = "transparent")

dev.off()

################################################################
#      Quinary (together) (EC10)            
################################################################

# pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC10_graphs//Quinary together (EC10)")
# 
# plot(Me1_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate",
#      xlab = "Sum of toxic units (EC10)",cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
#      ylim = c(0,1), xlim=c(0.1,100), col = "purple", lwd=2.5)
# title("Quinary", cex.main=1.6, line = 0.5)
# 
# plot(Me2_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate",
#      xlab = "Sum of toxic units (EC10)", 
#      ylim = c(0,1), xlim=c(0.1,100), col = "#d11141", lwd=2.5, add = TRUE)
# 
# plot(Me3_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate",
#      xlab = "Sum of toxic units (EC10)", 
#      ylim = c(0,1), xlim=c(0.1,100), col = "grey", lwd=2.5, add = TRUE)
# 
# plot(Me4_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate",
#      xlab = "Sum of toxic units (EC10)", 
#      ylim = c(0,1), xlim=c(0.1,100), col = "#40ECD0", lwd=2.5, add = TRUE)
# 
# plot(Me5_DRM_TU10, broken = TRUE, type = "none", ylab = "Relative growth rate",
#      xlab = "Sum of toxic units (EC10)", 
#      ylim = c(0,1), xlim=c(0.1,100), col = "black", lwd=2.5, add = TRUE)
# 
# plot(R_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
#      xlab = "Sum of toxic units (EC10)", pch = 1, lty=2, lwd=2,
#      ylim = c(0,1), xlim=c(0.1,100), col = "#00cc00", bg="#00cc00", add = TRUE)
# 
# plot(X_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
#      xlab = "Sum of toxic units (EC10)", pch = 1, lty=2, lwd=2,
#      ylim = c(0,1), xlim=c(0.1,100), col = "#f48fb1", bg="#f48fb1", add = TRUE)
# 
# legend(0.1, 0.42, cex=1.2, legend=c("Me2", "Me2", "Me3", "Me4", "Me5", "Quinary (ENR)", "Quinary (EQT)"), #"IA prediction", "CA prediction"
#        col = c("purple","#d11141", "grey", "#40ECD0", "black", "#00cc00", "#f48fb1"), 
#        pt.bg = c("purple","#d11141", "grey", "black", "#40ECD0", "#00cc00", "#f48fb1"),
#        lty=c(1, 1, 1, 1, 1, 2, 2),lwd=c(2.5, 2.5, 2.5, 2.5, 2.5, 2, 2),
#        pch=c(NA, NA, NA, NA, NA, 1, 1), box.col = "transparent")
# 
# dev.off()
# 

################################################################################
#Antagonism vs synergism #EC50
################################################################################

#BINARY -> 50
pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Antagonism vs synergism//50_Antagonism vs synergism")
par(mfrow= c(2,3), par(pty="s"), cex = 0.5)

plot(B$RGR, B$CA, xlab = "Observed relative growth rate", ylab = "Predicted relative growth rate", 
     cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1.1), xlim=c(0,1.1), pch = 23, col = "blue", bg="blue")  
title("Binary", cex.main=1.5, line = 0.5)

points(B$RGR, B$IA, xlab = "Observed relative growth rate", 
       ylab = "Predicted relative growth rate", cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
       ylim = c(0,1), xlim=c(0,1),pch = 22, col = "orange", bg="orange", add = TRUE) 
abline(coef = c(0, 1),col = "black", lwd = 1, lty=2)
text(0.2, 1.05, "Synergism", cex=1.2, lwd = 4)
text(0.9, 0.05, "Antagonism", cex=1.2, lwd = 4)
legend(0.85, 0.45, cex=1, legend=c("CA", "IA"), 
       col = c("blue","orange"), pt.bg = c("blue","orange"),pch=c(23, 22), box.col = "transparent")


#TERNARY -> 50
plot(T$RGR, T$CA, xlab = "Observed relative growth rate", 
     ylab = "Predicted relative growth rate", 
     cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1.1), xlim=c(0,1.1), pch = 23, col = "blue", bg="blue") 
title("Ternary", cex.main=1.5, line = 0.5)

points(T$RGR, T$IA, xlab = "Observed relative growth rate", 
       ylab = "Predicted relative growth rate", ylim = c(0,1.1), xlim=c(0,1.1), 
       pch = 22, col = "orange", bg="orange", add = TRUE) 
abline(coef = c(0, 1),col = "black", lwd = 1, lty=2)
text(0.2, 1.05, "Synergism", cex=1.2, lwd = 4)
text(0.9, 0.05, "Antagonism", cex=1.2, lwd = 4)
legend(0.85, 0.45, cex=1, legend=c("CA", "IA"), 
       col = c("blue","orange"), pt.bg = c("blue","orange"),pch=c(23, 22), box.col = "transparent")

#QUATERNARY -> 50
plot(Q$RGR, Q$CA, xlab = "Observed relative growth rate", 
     ylab = "Predicted relative growth rate", 
     cex.lab=1.3, cex.axis=1.2, cex.sub=1.3, 
     ylim = c(0,1.1), xlim=c(0,1.1), pch = 23, col = "blue", bg="blue") 
title("Quaternary", cex.main=1.5, line = 0.5)

points(Q$RGR, Q$IA, xlab = "Observed relative growth rate", 
       ylab = "Predicted relative growth rate", ylim = c(0,1.1), xlim=c(0,1.1), 
       pch = 22, col = "orange", bg="orange", add = TRUE) 
abline(coef = c(0, 1),col = "black", lwd = 1, lty=2)
text(0.2, 1.05, "Synergism", cex=1.2, lwd = 4)
text(0.9, 0.05, "Antagonism", cex=1.2, lwd = 4)
legend(0.85, 0.45, cex=1, legend=c("CA", "IA"), 
       col = c("blue","orange"), pt.bg = c("blue","orange"),pch=c(23, 22), box.col = "transparent")


#Quinary (equitoxic ray) -> 50 
plot(X$RGR, X$CA, xlab = "Observed relative growth rate", cex.lab=1.3, cex.axis=1.2, cex.main=1.5, cex.sub=1.3,
     ylab = "Predicted relative growth rate", 
     ylim = c(0,1.1), xlim=c(0,1.1), 
     pch = 23, col = "blue", bg="blue") 
title("Quinary (EQT)", cex.main=1.5, line = 0.5)

points(X$RGR, X$IA, xlab = "Observed relative growth rate", 
       ylab = "Predicted relative growth rate", ylim = c(0,1.1), xlim=c(0,1.1), 
       pch = 22, col = "orange", bg="orange", add = TRUE) 
abline(coef = c(0, 1),col = "black", lwd = 1, lty=2)
text(0.2, 1.05, "Synergism", cex=1.2, lwd = 4)
text(0.9, 0.05, "Antagonism", cex=1.2, lwd = 4)
legend(0.85, 0.45, cex=1, legend=c("CA", "IA"), 
       col = c("blue","orange"), pt.bg = c("blue","orange"),pch=c(23, 22), box.col = "transparent")


#Quinary (environmentally relevant ray) -> 50
# plot(R$RGR, R$CA, xlab = "Observed relative growth rate", 
#      cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
#      ylab = "Predicted relative growth rate", 
#      ylim = c(0,1.1), xlim=c(0,1.1), pch = 23, col = "blue", bg="blue") 
# title("Quinary (ENR)", cex.main=1.5, line = 0.5)
# 
# points(R$RGR, R$IA, xlab = "Observed relative growth rate", 
#        ylab = "Predicted relative growth rate", ylim = c(0,1.1), xlim=c(0,1.1), 
#        pch = 22, col = "orange", bg="orange", add = TRUE) 
# abline(coef = c(0, 1),col = "black", lwd = 1, lty=2)
# text(0.2, 1.05, "Synergism", cex=1.2, lwd = 4)
# text(0.9, 0.05, "Antagonism", cex=1.2, lwd = 4)
# legend(0.85, 0.45, cex=1, legend=c("CA", "IA"), 
#        col = c("blue","orange"), pt.bg = c("blue","orange"),pch=c(23, 22), box.col = "transparent")

dev.off()

#####################################################################################################################
#CA VS IA --> EC50
#####################################################################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Antagonism vs synergism//50_Antagonism vs synergism_CAvsIA")
par(mfrow= c(1,2), par(pty="s"), cex = 0.5) #1 column 2 rows

#CA -> 50
plot(B$RGR, B$CA,
     xlab = "Observed relative growth rate", ylab = "Predicted relative growth rate", 
     cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0,1.1), pch = 21, col = "#d11141", bg="#d11141") 
title("Concentration addition", cex.main=1.6, line = 0.5)

points(T$RGR, T$CA, pch = 24, col = "#00ccff", bg="#00ccff") 

points(Q$RGR, Q$CA, pch = 4, col = "#00cc00", bg="#00cc00") 

points(X$RGR, X$CA,pch = 23, col = "orange", bg="orange") 

points(R$RGR, R$CA, pch = 10, col = "#f48fb1", bg="#f48fb1") 

abline(coef = c(0, 1),col = "black", lwd = 1, lty=2)

text(0.21, 0.9, "Synergism", cex=1.2, lwd = 4, font=2)
text(0.87, 0.1, "Antagonism", cex=1.2, lwd = 4, font=2)

legend(0.83, 0.45, cex=1, legend=c("Binary", "Ternary", "Quaternary", "Quinary (EQT)", "Quinary (ENR)"),
       col = c("#d11141","#00ccff", "#00cc00", "orange", "#f48fb1"), 
       pt.bg = c("#d11141","#00ccff", "#00cc00", "orange", "#f48fb1"),
       pch=c(21, 24, 4, 23, 10), box.col = "transparent", y.intersp= 0.9)

#IA -> 50
plot(B$RGR, B$IA, add = TRUE, 
     xlab = "Observed relative growth rate", ylab = "Predicted relative growth rate", 
     cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0,1.1), pch = 21, col = "#d11141", bg="#d11141") 
title("Independent action", cex.main=1.6, line = 0.5)

points(T$RGR, T$IA, pch = 24, col = "#00ccff", bg="#00ccff") 

points(Q$RGR, Q$IA, pch = 4, col = "#00cc00", bg="#00cc00") 

points(X$RGR, X$IA,pch = 23, col = "orange", bg="orange") 

points(R$RGR, R$IA, pch = 10, col = "#f48fb1", bg="#f48fb1") 

abline(coef = c(0, 1),col = "black", lwd = 1, lty=2)

text(0.21, 0.9, "Synergism", cex=1.2, lwd = 4, font=2)
text(0.87, 0.1, "Antagonism", cex=1.2, lwd = 4, font=2)

legend(0.834, 0.45, cex=1, legend=c("Binary", "Ternary", "Quaternary", "Quinary (EQT)", "Quinary (ENR)"),
       col = c("#d11141","#00ccff", "#00cc00", "orange", "#f48fb1"), 
       pt.bg = c("#d11141","#00ccff", "#00cc00", "orange", "#f48fb1"),
       pch=c(21, 24, 4, 23, 10), box.col = "transparent", y.intersp= 0.9)

dev.off()


##########################################################################################################
#TOGETHER --> EC50
##########################################################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Antagonism vs synergism//50_Antagonism vs synergism_allCAvsIA")

plot(B$RGR, B$CA, xlab = "Observed relative growth rate", 
     ylab = "Predicted relative growth rate", 
     cex.lab=1.3, cex.axis=1.2, cex.main=1.6, cex.sub=1.3,
     ylim = c(0,1.1), xlim=c(0,1.1), pch = 23, col = "blue", bg="blue")  

points(B$RGR, B$IA, pch = 22, col = "orange", bg="orange") 

points(T$RGR, T$CA, pch = 23, col = "blue", bg="blue") 

points(T$RGR, T$IA, pch = 22, col = "orange", bg="orange") 

points(Q$RGR, Q$CA, pch = 23, col = "blue", bg="blue") 

points(Q$RGR, Q$IA, pch = 22, col = "orange", bg="orange") 

points(X$RGR, X$CA,pch = 23, col = "blue", bg="blue") 

points(X$RGR, X$IA, pch = 22, col = "orange", bg="orange") 

# points(R$RGR, R$CA, pch = 23, col = "blue", bg="blue") 
# 
# points(R$RGR, R$IA, pch = 22, col = "orange", bg="orange") 

abline(coef = c(0, 1),col = "black", lwd = 1, lty=2)

text(0.21, 1, "Synergism", cex=1.2, lwd = 4)
text(0.87, 0.1, "Antagonism", cex=1.2, lwd = 4)
legend(0, 0.8, cex=1.3, legend=c("CA", "IA"),
       col = c("blue","orange"), pt.bg = c("blue","orange"),pch=c(23, 22), box.col = "transparent")

dev.off()

##########################################################################################################
#together but different shapes 
###########################################################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Antagonism vs synergism//50_Antagonism vs synergism_CAvsIA_different shapes")

plot(B$RGR, B$CA,
     xlab = "Observed relative growth rate", ylab = "Predicted relative growth rate", 
     cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1.1), xlim=c(0,1.1), pch = 21, col = "blue", bg="blue") 
title("Experiment 3", cex.main=1.6, line = 0.5)

points(T$RGR, T$CA, pch = 22, col = "blue", bg="blue") 

points(Q$RGR, Q$CA, pch = 23, col = "blue", bg="blue") 

points(X$RGR, X$CA,pch = 24, col = "blue", bg="blue") 

# points(R$RGR, R$CA, pch = 25, col = "blue", bg="blue") 


points(B$RGR, B$IA, pch = 21, col = "orange", bg="orange") 

points(T$RGR, T$IA, pch = 22, col = "orange", bg="orange") 

points(Q$RGR, Q$IA, pch = 23, col = "orange", bg="orange") 

points(X$RGR, X$IA,pch = 24, col = "orange", bg="orange") 

# points(R$RGR, R$IA, pch = 25, col = "orange", bg="orange") 

abline(coef = c(0, 1),col = "black", lwd = 1, lty=2)

text(0.21, 0.9, "Synergism", cex=1.2, lwd = 4, font=2)
text(0.87, 0.1, "Antagonism", cex=1.2, lwd = 4, font=2)

text(0.21, 0.9, "Synergism", cex=1.2, lwd = 4, font=2)
text(0.87, 0.1, "Antagonism", cex=1.2, lwd = 4, font=2)

# legend(0.7, 0.45, cex=1, legend=c("Binary", "Ternary", "Quaternary", "Quinary (EQT)", "Quinary (ENR)"),
#        col = c("orange","orange", "orange", "orange", "orange"), 
#        pt.bg = c("orange","orange", "orange", "orange", "orange"),      
#        pch=c(21, 22, 23, 24, 24), box.col = "transparent", y.intersp= 0.9)

# legend(0.7, 0.45, cex=1, legend=c("Binary", "Ternary", "Quaternary", "Quinary (EQT)"),
#        col = c("orange","orange", "orange",  "orange"), 
#        pt.bg = c("orange","orange", "orange", "orange"),      
#        pch=c(21, 22, 23, 24), box.col = "transparent", y.intersp= 0.9)


dev.off()

# col = c("blue","blue", "blue", "blue", "blue"), 
# pt.bg = c("blue","blue", "blue", "blue", "blue"),

#        col = c("orange","orange", "orange", "orange", "orange"), 
#        pt.bg = c("orange","orange", "orange", "orange", "orange"),


#################################################################################
# MIF WITH TEXT
#################################################################################



#################################################################################
# slope comparison 
#################################################################################

# Results_summary_TU50_df <- as.data.frame(t(Results_summary_TU50))
# Results_summary_TU50_df$metal <- as.numeric(c(1:10))
# Results_summary_TU50_df$group <- c("SE", "SE", "SE", "SE", "SE", "M", "M", "M", "M", "M")
# 
# pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//Slope comparison", width = 7, height = 4.7)
# 
# ggplot(data = Results_summary_TU50_df, aes(x= metal, y = slope, color = group, shape = group)) +
#   ylim(0, 10) +
#   geom_errorbar(aes(ymin=slope-SE_slope, ymax=slope+SE_slope), 
#                 colour="black", alpha=0.7, width=0.05) + 
#   theme(legend.position = "none")+
#   theme(panel.background = element_blank(),strip.background = element_rect(colour=NA, fill=NA),
#         panel.border = element_rect(fill = NA, color = "black"), axis.title.x = element_blank()) +
#   scale_x_continuous(breaks=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
#                      labels=c("1" = "As", "2" = "Zn","3" = "Pb","4" = "Ag", "5" = "Cu", "6" = "B", "7" = "T", "8" = "Q", "9" = "R", "10" = "X")) +
#   geom_point(size = 3, stroke = 1) +
#   scale_color_manual(values = c("#00ccff","#d11141")) +
#   scale_shape_manual(values = c(21, 24),
#                      labels = c("SE", "M")) +
#   labs(y = "Slope")
# 
# dev.off()


########################################################################################
#STATISTIC (http://www.sthda.com/english/wiki/comparing-means-in-r)

##################################################################

########################################################################################
#for slides 

##################################################
#           Binary (EC10) for slides            
##################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//for_slides//2 metals (EC10)")

plot(B_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     ylim = c(0,1.02), xlim=c(0, 100), 
     xlab = "Sum of toxic units (EC10)", 
     lty=2, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     pch = 1, col = "#00cc00", bg="#00cc00") 
title("Binary", cex.main=1.6, line = 0.5)

axis(side = 2, at = 0.9, labels = expression("0.9"), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9))) #, col.axis = "darkred"
axis(side = 2, at = 0.5, labels = expression("0.5"), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9))) #, col.axis = "darkred"
axis(side = 1, at = MIF_df$MIF[MIF_df$name == B][1], labels = round(MIF_df$MIF[MIF_df$name == B][1], digits = 1), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9)))


plot(B_CA_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", 
     lty=1, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     ylim = c(0,1.02), xlim=c(0,100), 
     pch = 23, col = "blue", bg="blue", add = TRUE)

plot(B_IA_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", 
     lty=1, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     ylim = c(0,1.02), xlim=c(0,100), 
     pch = 22, col = "orange", bg="orange", add = TRUE)

points(MIF_df$MIF[MIF_df$name == B][1], MIF_df$RGR[MIF_df$name == B][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.8)

legend(0.095, 0.25, cex=1.25, legend=c("MIF", "IA prediction", "CA prediction","Observed Mn-Zn-Ba-Cr-Cd"), 
       col = c("#00ccff","orange", "blue", "#00cc00"), 
       pt.bg = c("#00ccff","orange", "blue", "#00cc00"),
       lty=c(NA, 1,1,2),lwd=c(4, 2,2,2),pch=c(4, 22, 23, 1), box.col = "transparent")

dev.off()

##################################################
#           Ternary (EC10) for slides            
##################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//for_slides//3 metals (EC10)")

plot(T_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     ylim = c(0,1.02), xlim=c(0 ,100), 
     xlab = "Sum of toxic units (EC10)", 
     lty=2, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     pch = 1, col = "#00cc00", bg="#00cc00") 
title("Ternary", cex.main=1.6, line = 0.5)

axis(side = 2, at = 0.9, labels = expression("0.9"), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9))) #, col.axis = "darkred"
axis(side = 2, at = 0.5, labels = expression("0.5"), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9))) #, col.axis = "darkred"
axis(side = 1, at = MIF_df$MIF[MIF_df$name == T][1], labels = round(MIF_df$MIF[MIF_df$name == T][1], digits = 1), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9)))


plot(T_CA_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", 
     lty=1, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     ylim = c(0,1.02), xlim=c(0,100), 
     pch = 23, col = "blue", bg="blue", add = TRUE)

plot(T_IA_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", 
     lty=1, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     ylim = c(0,1.02), xlim=c(0,100), 
     pch = 22, col = "orange", bg="orange", add = TRUE)

points(MIF_df$MIF[MIF_df$name == T][1], MIF_df$RGR[MIF_df$name == T][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.8)

legend(0.095, 0.25, cex=1.25, legend=c("MIF", "IA prediction", "CA prediction","Observed Mn-Zn-Ba-Cr-Cd"), 
       col = c("#00ccff","orange", "blue", "#00cc00"), 
       pt.bg = c("#00ccff","orange", "blue", "#00cc00"),
       lty=c(NA, 1,1,2),lwd=c(4, 2,2,2),pch=c(4, 22, 23, 1), box.col = "transparent")

dev.off()


##################################################
#           Quaternary (EC10) for slides            
##################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//for_slides//4 metals (EC10)")

plot(Q_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     ylim = c(0,1.02), xlim=c(0 , 100), 
     xlab = "Sum of toxic units (EC10)", 
     lty=2, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     pch = 1, col = "#00cc00", bg="#00cc00") 
title("Quaternary", cex.main=1.6, line = 0.5)

axis(side = 2, at = 0.9, labels = expression("0.9"), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9))) #, col.axis = "darkred"
axis(side = 2, at = 0.5, labels = expression("0.5"), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9))) #, col.axis = "darkred"
axis(side = 1, at = MIF_df$MIF[MIF_df$name == Q][1], labels = round(MIF_df$MIF[MIF_df$name == Q][1], digits = 1), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9)))


plot(Q_CA_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", 
     lty=1, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     ylim = c(0,1.02), xlim=c(0,100), 
     pch = 23, col = "blue", bg="blue", add = TRUE)

plot(Q_IA_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", 
     lty=1, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     ylim = c(0,1.02), xlim=c(0,100), 
     pch = 22, col = "orange", bg="orange", add = TRUE)

points(MIF_df$MIF[MIF_df$name == Q][1], MIF_df$RGR[MIF_df$name == Q][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.8)

legend(0.095, 0.25, cex=1.25, legend=c("MIF", "IA prediction", "CA prediction","Observed Mn-Zn-Ba-Cr-Cd"), 
       col = c("#00ccff","orange", "blue", "#00cc00"), 
       pt.bg = c("#00ccff","orange", "blue", "#00cc00"),
       lty=c(NA, 1,1,2),lwd=c(4, 2,2,2),pch=c(4, 22, 23, 1), box.col = "transparent")

dev.off()


##################################################
#           Quinary R (EC10) for slides            
##################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//for_slides//5 metals R (EC10)")

plot(R_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     ylim = c(0,1.02), xlim=c(0 , 100), 
     xlab = "Sum of toxic units (EC10)", 
     lty=2, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     pch = 1, col = "#00cc00", bg="#00cc00") 
title("Quinary (ENR)", cex.main=1.6, line = 0.5)

axis(side = 2, at = 0.9, labels = expression("0.9"), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9))) #, col.axis = "darkred"
axis(side = 2, at = 0.5, labels = expression("0.5"), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9))) #, col.axis = "darkred"
axis(side = 1, at = MIF_df$MIF[MIF_df$name == R][1], labels = round(MIF_df$MIF[MIF_df$name == R][1], digits = 1), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9)))


plot(R_CA_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", 
     lty=1, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     ylim = c(0,1.02), xlim=c(0,100), 
     pch = 23, col = "blue", bg="blue", add = TRUE)

plot(R_IA_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", 
     lty=1, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     ylim = c(0,1.02), xlim=c(0,100), 
     pch = 22, col = "orange", bg="orange", add = TRUE)

points(MIF_df$MIF[MIF_df$name == R][1], MIF_df$RGR[MIF_df$name == R][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.8)

legend(0.095, 0.25, cex=1.25, legend=c("MIF", "IA prediction", "CA prediction","Observed Mn-Zn-Ba-Cr-Cd"), 
       col = c("#00ccff","orange", "blue", "#00cc00"), 
       pt.bg = c("#00ccff","orange", "blue", "#00cc00"),
       lty=c(NA, 1,1,2),lwd=c(4, 2,2,2),pch=c(4, 22, 23, 1), box.col = "transparent")

dev.off()

##################################################
#           Quinary X (EC10) for slides            
##################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//for_slides//5 metals X (EC10)")

plot(X_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     ylim = c(0,1.02), xlim=c(0 , 100), 
     xlab = "Sum of toxic units (EC10)", 
     lty=2, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     pch = 1, col = "#00cc00", bg="#00cc00") 
title("Quaternary (EQT)", cex.main=1.6, line = 0.5)

axis(side = 2, at = 0.9, labels = expression("0.9"), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9))) #, col.axis = "darkred"
axis(side = 2, at = 0.5, labels = expression("0.5"), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9))) #, col.axis = "darkred"
axis(side = 1, at = MIF_df$MIF[MIF_df$name == B][1], labels = round(MIF_df$MIF[MIF_df$name == B][1], digits = 1), tick = TRUE, 
     col.ticks = "black", las =1, cex.axis = rep(1.4, length(0.9)))


plot(X_CA_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", 
     lty=1, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     ylim = c(0,1.02), xlim=c(0,100), 
     pch = 23, col = "blue", bg="blue", add = TRUE)

plot(X_IA_DRM_TU10, broken = TRUE, type = "all", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC10)", 
     lty=1, lwd=2, cex.lab=1.5, cex.axis=1.4, cex.sub=1.4,
     ylim = c(0,1.02), xlim=c(0,100), 
     pch = 22, col = "orange", bg="orange", add = TRUE)

points(MIF_df$MIF[MIF_df$name == X][1], MIF_df$RGR[MIF_df$name == X][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.8)

legend(0.095, 0.25, cex=1.25, legend=c("MIF", "IA prediction", "CA prediction","Observed Mn-Zn-Ba-Cr-Cd"), 
       col = c("#00ccff","orange", "blue", "#00cc00"), 
       pt.bg = c("#00ccff","orange", "blue", "#00cc00"),
       lty=c(NA, 1,1,2),lwd=c(4, 2,2,2),pch=c(4, 22, 23, 1), box.col = "transparent")

dev.off()


################################
#     Binary (EC20)            
################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC20_graphs//B vs TU (EC20)")

plot(Me1_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "purple", lwd=2.5) 
title("Binary", cex.main=1.6, line = 0.5)

plot(Me2_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)",
     ylim = c(0,1), xlim=c(0.05,100), col = "#d11141", lwd=2.5, add = TRUE)

plot(B_DRM_TU20, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.05,100), col = "#00cc00", bg="#00cc00", add = TRUE) 

plot(B_IA_DRM_TU20, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.05,100), col = "orange", bg="orange", add = TRUE)

plot(B_CA_DRM_TU20, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC20)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "blue", bg="blue", add = TRUE)

points(MIF20_df$MIF20[MIF20_df$name == B][1], MIF20_df$RGR[MIF20_df$name == B][1],
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.2)

# text(x=70, y=0.9, "CA/IA prediction:\n < observed     antagonism\n> observed    synergism", 
#      cex=1.5, adj=0.5)

legend(0.05, 0.32, cex=1.2, legend=c("Arsenic V", "Copper", "Observed As-Cu", "IA prediction", "CA prediction", "MIF"), 
       col = c("purple","#d11141", "#00cc00", "orange", "blue", "#00ccff"), 
       pt.bg = c("purple","#d11141", "#00cc00", "orange", "blue", "#00ccff"),
       lty=c(1, 1, 2, NA, NA, NA),lwd=c(2.5, 2.5, 2, NA, NA, 4),pch=c(NA,NA, 1, 22, 23, 4), box.col = "transparent")

dev.off()

################################
#           Ternary (EC20)            
################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC20_graphs//T vs TU (EC20)")

plot(Me1_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "purple", lwd=2.5)  #, ylim = c(0,1), xlim=c(0.01,20)
title("Ternary", cex.main=1.6, line = 0.5)

plot(Me2_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#d11141", lwd=2.5, add = TRUE)

plot(Me3_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "grey", lwd=2.5, add = TRUE) 

plot(T_DRM_TU20, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.05,100), col = "#00cc00", bg="#00cc00", add = TRUE) 

plot(T_IA_DRM_TU20, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.05,100), col = "orange", bg="orange", add = TRUE) 

plot(T_CA_DRM_TU20, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC20)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "blue", bg="blue", add = TRUE)

points(MIF20_df$MIF20[MIF20_df$name == T][1], MIF20_df$RGR[MIF20_df$name == T][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.5)


legend(0.05, 0.37, cex=1.2, legend=c("Arsenic V", "Copper", "Lead", "Observed As-Cu-Pb", "IA prediction", "CA prediction", "MIF"), 
       col = c("purple","#d11141", "grey", "#00cc00", "orange", "blue", "#00ccff"), 
       pt.bg = c("purple","#d11141", "grey", "#00cc00", "orange", "blue", "#00ccff"),
       lty=c(1, 1, 1, 2, NA, NA, NA),lwd=c(2.5, 2.5, 2.5, 2, NA, NA, 4),pch=c(NA,NA, NA, 1, 22, 23, 4), box.col = "transparent")

dev.off()

################################
#           Quaternary (EC20)            
################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC20_graphs//Q vs TU (EC20)")

plot(Me1_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "purple", lwd=2.5) 
title("Quaternary", cex.main=1.6, line = 0.5)

plot(Me2_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#d11141", lwd=2.5, add = TRUE)

plot(Me3_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "grey", lwd=2.5, add = TRUE) 

plot(Me4_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#40ECD0", lwd=2.5, add = TRUE) 

plot(Q_DRM_TU20, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)",  pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.05,100), col = "#00cc00", bg="#00cc00", add = TRUE) 

plot(Q_IA_DRM_TU20, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)",  pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.1,100), col = "orange", bg="orange", add = TRUE) 

plot(Q_CA_DRM_TU20, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC20)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "blue", bg="blue", add = TRUE)

points(MIF20_df$MIF20[MIF20_df$name == Q][1], MIF20_df$RGR[MIF20_df$name == Q][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.5)


legend(0.05, 0.42, cex=1.2, legend=c("Arsenic V", "Copper", "Lead", "Nickel", "Observed As-Cu-Pb-Ni", "IA prediction", "CA prediction", "MIF"), 
       col = c("purple","#d11141", "grey", "#40ECD0", "#00cc00", "orange", "blue", "#00ccff"), 
       pt.bg = c("purple","#d11141", "grey", "#40ECD0", "#00cc00", "orange", "blue", "#00ccff"),
       lty=c(1, 1, 1, 1, 2, NA, NA, NA),lwd=c(2.5, 2.5, 2.5, 2.5, 2, NA, NA, 4),pch=c(NA,NA, NA,NA, 1, 22, 23, 4), box.col = "transparent")

dev.off()


########################################
#   Quinary (equitoxic ray) (EC20)            
########################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC20_graphs//X vs TU (EC20)")

plot(Me1_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)",   cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "purple", lwd=2.5) 
title("Quinary (equitoxic ray)", cex.main=1.6, line = 0.5)

plot(Me2_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#d11141", lwd=2.5, add = TRUE)

plot(Me3_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "grey", lwd=2.5, add = TRUE) 

plot(Me4_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#40ECD0", lwd=2.5, add = TRUE) 

plot(Me5_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05, 100), col = "black", lwd=2.5, add = TRUE) 

plot(X_DRM_TU20, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.05,100), col = "#00cc00", bg="#00cc00", add = TRUE)  

plot(X_IA_DRM_TU20, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.05,100), col = "orange", bg="orange", add = TRUE) 

plot(X_CA_DRM_TU20, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC20)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "blue", bg="blue", add = TRUE)

points(MIF20_df$MIF20[MIF20_df$name == X][1], MIF20_df$RGR[MIF20_df$name == X][1], 
       pch = 4, col = "#00ccff", bg="#00ccff", lwd=4, cex=1.5)


legend(0.05, 0.47, cex=1.2, legend=c("Arsenic V", "Copper", "Lead", "Nickel", "Cadmium", "Observed As-Cu-Pb-Ni-Cd", "IA prediction", "CA prediction", "MIF"), 
       col = c("purple","#d11141", "grey", "#40ECD0", "black", "#00cc00", "orange", "blue", "#00ccff"), 
       pt.bg = c("purple","#d11141", "grey",  "#40ECD0","black", "#00cc00", "orange", "blue", "#00ccff"),
       lty=c(1, 1, 1, 1, 1, 2, NA, NA, NA),lwd=c(2.5, 2.5, 2.5, 2.5, 2.5, 2, NA, NA, 4),pch=c(NA, NA, NA, NA, NA, 1, 22, 23, 4), box.col = "transparent")

dev.off()

################################################################
#      Quinary (environmentally relevant ray) (EC20)            
################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//EC20_graphs//R vs TU (EC20)")

plot(Me1_DRM_TU20, broken = TRUE, type = "none", 
     ylab = "Relative growth rate", xlab = "Sum of toxic units (EC20)",  
     cex.lab=1.3, cex.axis=1.2, cex.sub=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "purple", lwd=2.5) 
title("Quinary (environmentally relevant ray)", cex.main=1.6, line = 0.5)


plot(Me2_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#d11141", lwd=2.5, add = TRUE)

plot(Me3_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "grey", lwd=2.5, add = TRUE) 

plot(Me4_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", 
     ylim = c(0,1), xlim=c(0.05,100), col = "#40ECD0", lwd=2.5, add = TRUE) 

plot(Me5_DRM_TU20, broken = TRUE, type = "none", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)",
     ylim = c(0,1), xlim=c(0.05, 100), col = "black", lwd=2.5, add = TRUE) 

plot(R_DRM_TU20, broken = TRUE, type = "all", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", pch = 1, lty=2, lwd=2,
     ylim = c(0,1), xlim=c(0.05,100), col = "#00cc00", bg="#00cc00", add = TRUE)

plot(R_IA_DRM_TU20, broken = TRUE, type = "obs", ylab = "Relative growth rate", 
     xlab = "Sum of toxic units (EC20)", pch = 22, cex=1.5,
     ylim = c(0,1), xlim=c(0.05,20), col = "orange", bg="orange", add = TRUE) 

plot(R_CA_DRM_TU20, broken = TRUE, type = "obs", ylab = "Relative growth rate",
     xlab = "Sum of toxic units (EC20)", pch = 23, cex=1.3,
     ylim = c(0,1), xlim=c(0.05,100), col = "blue", bg="blue", add = TRUE)

points(MIF20_df$MIF20[MIF20_df$name == R][1], MIF20_df$RGR[MIF20_df$name == R][1], 
       pch = 7, col = "#f48fb1", bg="#f48fb1", lwd=1.7, cex=1.7)  

legend(0.05, 0.47, cex=1.2, legend=c("Arsenic V", "Copper", "Lead", "Nickel", "Cadmium", "Observed As-Cu-Pb-Ni-Cd", "IA prediction", "CA prediction", "MIF"), 
       col = c("purple","#d11141", "grey", "#40ECD0", "black", "#00cc00", "orange", "blue", "#f48fb1"), 
       pt.bg = c("purple","#d11141", "grey",  "#40ECD0","black", "#00cc00", "orange", "blue", "#f48fb1"),
       lty=c(1, 1, 1, 1, 1, 2, NA, NA, NA),lwd=c(2.5, 2.5, 2.5, 2.5, 2.5, 2, NA, NA, 1.7),pch=c(NA, NA, NA, NA, NA, 1, 22, 23, 7), box.col = "transparent")

dev.off()

#################################################################################
# MIF WITH TEXT
#################################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//MIF_withtext", width = 3, height = 2.7)

ggplot(data = MIF_df, aes(x = n, y = MIF, color = group, shape = group, size = group)) +
  ylim(-0.1, 3.7) + xlim(1.95,5.56) +
  geom_errorbar(aes(x=n, ymin=MIF-MIF_SE, ymax=MIF+MIF_SE), 
                colour="black", alpha=0.7, width=0.06, size=0.34) +  
  geom_point(size = 2, stroke = 0.7) +
  scale_color_manual(values = c("#f48fb1", "#00ccff")) +
  scale_size_manual(values= c(2, 10)) +
  scale_shape_manual(values = c(7, 4),
                     labels = c("Environmentally relevant ray", "Equitoxic ray")) +
  labs(x = "No. of metals", y = "MIF (EC10)") +
  theme(panel.background = element_blank(),strip.background = element_rect(colour=NA, fill=NA),
        panel.border = element_rect(fill = NA, color = "black")) +
  guides(color = guide_legend(override.aes = list(size = 2.7, fill = "white"))) +
  theme(legend.position = c(0.04, 0.21),
        legend.key = element_rect(fill = "white"),
        legend.justification = c(0, 1),
        legend.box.background = element_rect(fill = "white", color = NA),
        legend.margin = margin(0,0,0,0),
        legend.title = element_blank(),
        legend.key.size = unit(0.4, "cm"),
        axis.title = element_text(size = 9)) +
  annotate("text", x = 5.45, y = 1.2, label = "antagonistic", angle = 90, hjust = 0, size = 2.3)  + # Annotate rotated text label
  annotate("text", x = 5.45, y = 0.8, label = "synergistic", angle = 90, hjust = 1, size = 2.3)  + 
  annotate("text", x = 5.42, y = 1.01, label = "additive", angle = 0, hjust = 0.5, size = 2.3, vjust = "center")  + # Annotate rotated text label
  # Annotate rotated text label
  # geom_label(aes(x = 5.4, y = 1.02, label = "additivity"),  angle = 0, hjust = 0.5,, 
  #            size = 2.3, fill = NA, color ="black", label.size = NA) +
  geom_segment(aes(x=5.56, y=1.2, xend=5.56, yend=2.2), arrow = arrow(length=unit(.10, 'cm')), color='black', lwd=0.2) + 
  geom_segment(aes(x=5.56, y=0.8, xend=5.56, yend=-0.1), arrow = arrow(length=unit(.10, 'cm')), color='black', lwd=0.2) + 
  geom_segment(x=1.8, y=1, xend=5.05, yend=1, col = "black", size = 0.4, linetype = "dashed")  #size =0.4,

dev.off()

############################################################################################################
#MIF WITHOUT TEXT
########################################################################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//for_slides//MIF_withouttext", width = 3, height = 2.7)

ggplot(data = MIF_df, aes(x = n, y = MIF, color = group, shape = group, size = group)) +
  ylim(-0.1, 3.7) + xlim(1.95,5.05) +
  geom_errorbar(aes(x=n, ymin=MIF-MIF_SE, ymax=MIF+MIF_SE), 
                colour="black", alpha=0.7, width=0.06, size=0.34) + 
  geom_point(size = 2, stroke = 0.7) +
  scale_color_manual(values = c("#f48fb1", "#00ccff")) +
  scale_size_manual(values= c(2, 10)) +
  scale_shape_manual(values = c(7, 4),
                     labels = c("Environmentally relevant ray", "Equitoxic ray")) +
  labs(x = "No. of metals", y = "MIF (EC10)") +
  theme(panel.background = element_blank(),strip.background = element_rect(colour=NA, fill=NA),
        panel.border = element_rect(fill = NA, color = "black", size=0.15),
        axis.text = element_text(size = 9), 
        axis.line = element_line(size = 0.10),
        axis.ticks = element_line(size = 0.15),
        axis.title = element_text(size = 8)) +
  guides(color = guide_legend(override.aes = list(size = 2.7, fill = "white"))) +
  #theme(legend.position = "none") +
  theme(legend.position = c(0.04, 0.21),
        legend.key = element_rect(fill = "white"),
        legend.justification = c(0, 1),
        legend.box.background = element_rect(fill = "white", color = NA),
        legend.margin = margin(0,0,0,0),
        legend.title = element_blank(),
        legend.key.size = unit(0.4, "cm"),
        legend.text = element_text(size = 6.5)) +
  geom_hline(yintercept=1, linetype="dashed", color = "black", alpha =0.5, size = 0.22)

dev.off()


############################################################################################################
#MIF20 WITHOUT TEXT
########################################################################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//for_slides//MIF20_withouttext", width = 3, height = 2.7)

ggplot(data = MIF20_df, aes(x = n, y = MIF20, color = group, shape = group, size = group)) +
  ylim(-0.1, 3.7) + xlim(1.95,5.05) +
  geom_errorbar(aes(x=n, ymin=MIF20-MIF20_SE, ymax=MIF20+MIF20_SE), 
                colour="black", alpha=0.7, width=0.06, size=0.34) + 
  geom_point(size = 2, stroke = 0.7) +
  scale_color_manual(values = c("#f48fb1", "#00ccff")) +
  scale_size_manual(values= c(2, 10)) +
  scale_shape_manual(values = c(7, 4),
                     labels = c("Environmentally relevant ray", "Equitoxic ray")) +
  labs(x = "No. of metals", y = "MIF (EC20)") +
  theme(panel.background = element_blank(),strip.background = element_rect(colour=NA, fill=NA),
        panel.border = element_rect(fill = NA, color = "black", size=0.15),
        axis.text = element_text(size = 9), 
        axis.line = element_line(size = 0.10),
        axis.ticks = element_line(size = 0.15),
        axis.title = element_text(size = 8)) +
  guides(color = guide_legend(override.aes = list(size = 2.7, fill = "white"))) +
  #theme(legend.position = "none") +
  theme(legend.position = c(0.04, 0.21),
        legend.key = element_rect(fill = "white"),
        legend.justification = c(0, 1),
        legend.box.background = element_rect(fill = "white", color = NA),
        legend.margin = margin(0,0,0,0),
        legend.title = element_blank(),
        legend.key.size = unit(0.4, "cm"),
        legend.text = element_text(size = 6.5)) +
  geom_hline(yintercept=1, linetype="dashed", color = "black", alpha =0.5, size = 0.22)

dev.off()

#################################################################################
#MIF20 WITH TEXT
#################################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//MIF20_withtext", width = 3, height = 2.7)

ggplot(data = MIF20_df, aes(x = n, y = MIF20, color = group, shape = group, size = group)) +
  ylim(-0.1, 3.7) + xlim(1.95,5.56) +
  geom_errorbar(aes(x=n, ymin=MIF20-MIF20_SE, ymax=MIF20+MIF20_SE), 
                colour="black", alpha=0.7, width=0.06, size=0.34) +  
  geom_point(size = 2, stroke = 0.7) +
  scale_color_manual(values = c("#f48fb1", "#00ccff")) +
  scale_size_manual(values= c(2, 10)) +
  scale_shape_manual(values = c(7, 4),
                     labels = c("Environmentally relevant ray", "Equitoxic ray")) +
  labs(x = "No. of metals", y = "MIF (EC20)") +
  theme(panel.background = element_blank(),strip.background = element_rect(colour=NA, fill=NA),
        panel.border = element_rect(fill = NA, color = "black")) +
  guides(color = guide_legend(override.aes = list(size = 2.7, fill = "white"))) +
  theme(legend.position = c(0.04, 0.21),
        legend.key = element_rect(fill = "white"),
        legend.justification = c(0, 1),
        legend.box.background = element_rect(fill = "white", color = NA),
        legend.margin = margin(0,0,0,0),
        legend.title = element_blank(),
        legend.key.size = unit(0.4, "cm"),
        axis.title = element_text(size = 9)) +
  annotate("text", x = 5.45, y = 1.2, label = "antagonistic", angle = 90, hjust = 0, size = 2.3)  + # Annotate rotated text label
  annotate("text", x = 5.45, y = 0.8, label = "synergistic", angle = 90, hjust = 1, size = 2.3)  + 
  annotate("text", x = 5.42, y = 1.01, label = "additive", angle = 0, hjust = 0.5, size = 2.3, vjust = "center")  + # Annotate rotated text label
  # Annotate rotated text label
  # geom_label(aes(x = 5.4, y = 1.02, label = "additivity"),  angle = 0, hjust = 0.5,, 
  #            size = 2.3, fill = NA, color ="black", label.size = NA) +
  geom_segment(aes(x=5.56, y=1.2, xend=5.56, yend=2.2), arrow = arrow(length=unit(.10, 'cm')), color='black', lwd=0.2) + 
  geom_segment(aes(x=5.56, y=0.8, xend=5.56, yend=-0.1), arrow = arrow(length=unit(.10, 'cm')), color='black', lwd=0.2) + 
  geom_segment(x=1.8, y=1, xend=5.05, yend=1, col = "black", size = 0.4, linetype = "dashed")  #size =0.4,

dev.off()


#################################################################################
# MIF IA WITH TEXT
#################################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//MIF_IA_withtext", width = 3, height = 2.7)

ggplot(data = MIF_IA_df, aes(x = n, y = MIF_IA, color = group, shape = group, size = group)) +
  ylim(-0.1, 3.7) + xlim(1.95,5.56) +
  geom_errorbar(aes(x=n, ymin=MIF_IA-MIF_IA_SE, ymax=MIF_IA+MIF_IA_SE), 
                colour="black", alpha=0.7, width=0.06, size=0.34) +  
  geom_point(size = 2, stroke = 0.7) +
  scale_color_manual(values = c("#f48fb1", "#00ccff")) +
  scale_size_manual(values= c(2, 10)) +
  scale_shape_manual(values = c(7, 4),
                     labels = c("Environmentally relevant ray", "Equitoxic ray")) +
  labs(x = "No. of metals", y = "MIF (IA, EC10)") +
  theme(panel.background = element_blank(),strip.background = element_rect(colour=NA, fill=NA),
        panel.border = element_rect(fill = NA, color = "black")) +
  guides(color = guide_legend(override.aes = list(size = 2.7, fill = "white"))) +
  theme(legend.position = "none") +
  #if you want back the legend hashtag the previous line and unhashtag the following:
  # theme(legend.position = c(0.04, 0.21),
  #     legend.key = element_rect(fill = "white"),
  #     legend.justification = c(0, 1),
  #     legend.box.background = element_rect(fill = "white", color = NA),
  #     legend.margin = margin(0,0,0,0),
  #     legend.title = element_blank(),
  #     legend.key.size = unit(0.4, "cm"),
  #     axis.title = element_text(size = 9)) +
  annotate("text", x = 5.45, y = 1.2, label = "antagonistic", angle = 90, hjust = 0, size = 2.3)  + # Annotate rotated text label
  annotate("text", x = 5.45, y = 0.8, label = "synergistic", angle = 90, hjust = 1, size = 2.3)  + 
  annotate("text", x = 5.42, y = 1.01, label = "additive", angle = 0, hjust = 0.5, size = 2.3, vjust = "center")  + # Annotate rotated text label
  # Annotate rotated text label
  # geom_label(aes(x = 5.4, y = 1.02, label = "additivity"),  angle = 0, hjust = 0.5,, 
  #            size = 2.3, fill = NA, color ="black", label.size = NA) +
  geom_segment(aes(x=5.56, y=1.2, xend=5.56, yend=2.2), arrow = arrow(length=unit(.10, 'cm')), color='black', lwd=0.2) + 
  geom_segment(aes(x=5.56, y=0.8, xend=5.56, yend=-0.1), arrow = arrow(length=unit(.10, 'cm')), color='black', lwd=0.2) + 
  geom_segment(x=1.8, y=1, xend=5.05, yend=1, col = "black", size = 0.4, linetype = "dashed")  #size =0.4,

dev.off()

############################################################################################################
#MIF_IA IA WITHOUT TEXT
########################################################################################################################

pdf("C://Users//dedol//OneDrive//Desktop//MEED-mix project//Experiments//Results experiments//2023-06-05_As,Zn,Pb,Ag,Cu+M+SE+B+T+Q1+Q2//R_050623//Results//for_slides//MIF_IA_withouttext", width = 3, height = 2.7)

ggplot(data = MIF_IA_df, aes(x = n, y = MIF_IA, color = group, shape = group, size = group)) +
  ylim(-0.1, 3.7) + xlim(1.95,5.05) +
  geom_errorbar(aes(x=n, ymin=MIF_IA-MIF_IA_SE, ymax=MIF_IA+MIF_IA_SE), 
                colour="black", alpha=0.7, width=0.06, size=0.34) + 
  geom_point(size = 2, stroke = 0.7) +
  scale_color_manual(values = c("#f48fb1", "#00ccff")) +
  scale_size_manual(values= c(2, 10)) +
  scale_shape_manual(values = c(7, 4),
                     labels = c("Environmentally relevant ray", "Equitoxic ray")) +
  labs(x = "No. of metals", y = "MIF (IA, EC10)") +
  theme(panel.background = element_blank(),strip.background = element_rect(colour=NA, fill=NA),
        panel.border = element_rect(fill = NA, color = "black", size=0.15),
        axis.text = element_text(size = 9), 
        axis.line = element_line(size = 0.10),
        axis.ticks = element_line(size = 0.15),
        axis.title = element_text(size = 8)) +
  guides(color = guide_legend(override.aes = list(size = 2.7, fill = "white"))) +
  theme(legend.position = "none") +
  #if you want back the legend hashtag the previous line and unhashtag the following:
  # theme(legend.position = c(0.04, 0.21),
  #       legend.key = element_rect(fill = "white"),
  #       legend.justification = c(0, 1),
  #       legend.box.background = element_rect(fill = "white", color = NA),
  #       legend.margin = margin(0,0,0,0),
  #       legend.title = element_blank(),
  #       legend.key.size = unit(0.4, "cm"),
  #       legend.text = element_text(size = 6.5)) +
  geom_hline(yintercept=1, linetype="dashed", color = "black", alpha =0.5, size = 0.22)

dev.off()