if(!require(multiplestressR)){
  install.packages("multiplestressR")
  library(multiplestressR)
}

if(!require(dplyr)){
  install.packages("dplyr")
  library(dplyr)
}

if(!require(tidyverse)){
  install.packages("tidyverse")
  library(tidyverse)
}

if(!require(readxl)){
  install.packages("readxl")
  library(readxl)
}


if(!require(metafor)){
  install.packages("metafor")
  library(metafor)
}


if(!require(forestplot)){
  install.packages("forestplot")
  library(forestplot)
}

if(!require(ggplot2)){
  install.packages("ggplot2")
  library(ggplot2)
}


if(!require(MuMIn)){
  install.packages("MuMIn")
  library(MuMIn)
}


if(!require(tidyr)){
  install.packages("tidyr")
  library(tidyr)
}


if(!require(esc)){
  install.packages("esc")
  library(esc)
}

if(!require(clubSandwich)){
  install.packages("clubSandwich")
  library(clubSandwich)
}


if(!require(dmetar)){
  install.packages("dmetar")
  library(dmetar)
}

options(scipen=999)

setwd("") ### Your working directory

##### Load extracted means, sample sizes, and sd

Stressor_stats_df <- read_excel("AppendixA.xlsx", sheet=2)
Stressor_stats_df<-Stressor_stats_df[-75,]

###Load all accepted manuscripts (with and without required data, first 2 criteria)


### Run MultiplestressR analysis

Stressor_effect<- effect_size_additive(Control_N= Stressor_stats_df$Control_N,           
                                       Control_SD= Stressor_stats_df$Control_SD,    
                                       Control_Mean= Stressor_stats_df$Control_Mean,                  
                                       StressorA_N= Stressor_stats_df$StressorA_N,         
                                       StressorA_SD= Stressor_stats_df$StressorA_SD,  
                                       StressorA_Mean= Stressor_stats_df$StressorA_Mean,                
                                       StressorB_N= Stressor_stats_df$StressorB_N,         
                                       StressorB_SD= Stressor_stats_df$StressorB_SD,  
                                       StressorB_Mean= Stressor_stats_df$StressorB_Mean,                
                                       StressorsAB_N= Stressor_stats_df$StressorsAB_N,       
                                       StressorsAB_SD= Stressor_stats_df$StressorsAB_SD,
                                       StressorsAB_Mean = Stressor_stats_df$StressorsAB_Mean,
                                       Small_Sample_Correction= TRUE,
                                       Significance_Level      = 0.05)


###Create duplicate of stressor stats df just in case


Stressor_stats_meta <- Stressor_stats_df

colnames_merge <- colnames(Stressor_effect)


Stressor_stats_meta[,colnames_merge]<-effect_size_additive(Control_N= Stressor_stats_df$Control_N,           
                                                           Control_SD= Stressor_stats_df$Control_SD,    
                                                           Control_Mean= Stressor_stats_df$Control_Mean,                  
                                                           StressorA_N= Stressor_stats_df$StressorA_N,         
                                                           StressorA_SD= Stressor_stats_df$StressorA_SD,  
                                                           StressorA_Mean= Stressor_stats_df$StressorA_Mean,                
                                                           StressorB_N= Stressor_stats_df$StressorB_N,         
                                                           StressorB_SD= Stressor_stats_df$StressorB_SD,  
                                                           StressorB_Mean= Stressor_stats_df$StressorB_Mean,                
                                                           StressorsAB_N= Stressor_stats_df$StressorsAB_N,       
                                                           StressorsAB_SD= Stressor_stats_df$StressorsAB_SD,
                                                           StressorsAB_Mean = Stressor_stats_df$StressorsAB_Mean,
                                                           Small_Sample_Correction= TRUE,
                                                           Significance_Level      = 0.05)

Stressor_effect_class <- classify_interactions(effect_size_dataframe = Stressor_effect,
                                               assign_reversals = TRUE,
                                               remove_directionality = TRUE)

colnames_merge_2 <- colnames(Stressor_effect_class)

Stressor_stats_removed_directionality <- Stressor_stats_meta


### Without directionality

Stressor_stats_removed_directionality[,colnames_merge_2] <- classify_interactions(effect_size_dataframe = Stressor_stats_meta[,colnames_merge],
                                                                                  assign_reversals = TRUE, 
                                                                                  remove_directionality = TRUE)

Stressor_stats_with_directionality <- Stressor_stats_meta


###With directionality

Stressor_stats_with_directionality[,colnames_merge_2] <-classify_interactions(effect_size_dataframe = Stressor_stats_meta[,colnames_merge],
                                                                              assign_reversals = TRUE, 
                                                                              remove_directionality = FALSE)


Stressor_stats_removed_directionality <- Stressor_stats_removed_directionality%>% 
  relocate(any_of(c("Code_Manuscript", "Subextraction", "ID_Code", "Year","First_author")))

Stressor_stats_with_directionality <- Stressor_stats_with_directionality%>% 
  relocate(any_of(c("Code_Manuscript", "Subextraction", "ID_Code", "Year","First_author")))


###Add scaled columns for continuous variables

Stressor_stats_removed_directionality$Duration_sc <- scale(as.numeric(Stressor_stats_removed_directionality$Duration))
#Stressor_stats_removed_directionality$Colonization_sc <- scale(as.numeric(Stressor_stats_removed_directionality$Time.of.colonization))
#Stressor_stats_removed_directionality$acclimation_sc <- scale(as.numeric(Stressor_stats_removed_directionality$Time.of.acclimation.of.MZB..max.))
Stressor_stats_removed_directionality$litter_sc <- scale(as.numeric(Stressor_stats_removed_directionality$Litter_mass))
#Stressor_stats_removed_directionality$volume_sc <- scale(as.numeric(Stressor_stats_removed_directionality$Volume_ml))



Stressor_stats_removed_directionality$ID_Code <- 1:nrow(Stressor_stats_removed_directionality)
### Visual inspection

dfa_plots <- summary_plots(effect_size_dataframe = Stressor_stats_removed_directionality[,colnames_merge_2],
                           Small_Sample_Correction = TRUE,
                           Significance_Level = 0.05)
dfa_plots


dat<-Stressor_stats_removed_directionality

###Keep only datapoints of longest duration for each manuscript and treatment

datred <- dat %>%
  arrange(Code_manuscript, desc(Duration)) %>%
  group_by(Code_manuscript) %>%
  filter(Duration == max(Duration)) %>%
  ungroup()%>%
  mutate(Code_manuscript= as.numeric(sub('.*_', '',Code_manuscript)))%>%
  mutate(pooled_n= StressorA_N+StressorB_N+StressorsAB_N, Control_N)


#### Now fit models

### In order to run a multivariate meta-analysis we require a var-cov matrix

###We calculate the first one with the ClubSandwich package

V <- impute_covariance_matrix(datred$Interaction_Variance,
                              cluster=datred$Code_manuscript,
                              r=0.6, ### assume a 0.6 correlation
                              smooth_vi = FALSE,
                              check_PD = TRUE,
                              return_list = FALSE
)


rand.mv <- rma.mv(yi=Interaction_Effect_Size, 
                  V,
                  data = datred,
                  random= c(~1|Code_manuscript/ID_Code),
                  test = "t", 
                  method = "REML"
                  #struct= "CAR",
                  # rho=0.
)





randrob <- robust(rand.mv, # Note that this reports model-based (not robust) standard errors
                  cluster=Code_manuscript, clubSandwich = TRUE) 

randrob ### These are robust

mlm.variance.distribution(rand.mv)

I2rand <-mlm.variance.distribution(rand.mv)

Macroinvertebrates <- rma.mv(yi=Interaction_Effect_Size, 
                             V,
                             #W=W,
                             data = datred,
                             mods= ~Macroinvertebrates-1,
                             random= c(~1|Code_manuscript/ID_Code),
                             test = "t", 
                             method = "REML"
)

mlm.variance.distribution(Macroinvertebrates) ###I squared values


MIrob <- robust(Macroinvertebrates, # Note that this reports model-based (not robust) standard errors
                cluster=Code_manuscript, clubSandwich = TRUE) 
MIrob ### These are robust



Dispersal <- rma.mv(yi=Interaction_Effect_Size, 
                    V=V,
                    # W=W,
                    data = datred,
                    mods= ~Dispersal_allowed-1,
                    random= c(~1|Code_manuscript/ID_Code),
                    test = "t", 
                    method = "REML"
)


Drob <-robust(Dispersal, # Note that this reports model-based (not robust) standard errors
              cluster = Code_manuscript, clubSandwich = TRUE) 

Drob ### These are robust

mlm.variance.distribution(Dispersal) 

Mycorrhizae <- rma.mv(yi=Interaction_Effect_Size, 
                      V,
                      #W=W,
                      data = datred,
                      mods= ~Mycorrhizae-1,
                      random= c(~1|Code_manuscript/ID_Code),
                      test = "t", 
                      method = "REML"
)


Myrob <- robust(Mycorrhizae, # Note that this reports model-based (not robust) standard errors
                cluster = Code_manuscript) 
Myrob ### These are robust

mlm.variance.distribution(Mycorrhizae)

Habitat <- rma.mv(yi=Interaction_Effect_Size, 
                  V,
                  #W=W,
                  data = datred,
                  mods= ~Habitat-1,
                  random= c(~1|Code_manuscript/ID_Code),
                  test = "t", 
                  method = "REML"
)


Hrob <- robust(Habitat, # Note that this reports model-based (not robust) standard errors
               cluster = Code_manuscript, clubSandwich = TRUE) 
Hrob ### These are robust



mlm.variance.distribution(Habitat)



Duration <-rma.mv(yi=Interaction_Effect_Size, 
                  V,
                  #W=W,
                  data = datred,
                  mods= ~Duration_sc,
                  random= c(~1|Code_manuscript/ID_Code),
                  test = "t", 
                  method = "REML"
)



Duration


mlm.variance.distribution(Duration)

Litter <-rma.mv(yi=Interaction_Effect_Size, 
                V,
                #W=W,
                data = datred,
                mods= ~litter_sc,
                random= c( ~1|Code_manuscript/ID_Code),
                test = "t", 
                method = "REML"
)

Litter


mlm.variance.distribution(Litter) ###I squared values


###This is how the weight is assigned in a multivariate analysis

wi <- weights(rand.mv, type="rowsum")
data.frame(k = c(table(datred$Code_manuscript)),
           weight = tapply(wi, datred$Code_manuscript, sum))



Eggtest <- rma.mv(yi=Interaction_Effect_Size, 
                  V=V,
                  data = datred,
                  mods=Interaction_Variance,
                  random= c(~1|Code_manuscript/ID_Code),
                  test = "t", 
                  method = "REML"
)


Eggtest
coef_test(Eggtest,vcov = "CR2", p_values = T)
#### No significant correlation between interaction variance and interaction effect size, which can be interpreted as no bias

Eggtest2 <- rma.mv(yi=Interaction_Effect_Size, 
                   V=V,
                   data = datred,
                   mods=I(1/pooled_n),
                   random= c(~1|Code_manuscript/ID_Code),
                   test = "t", 
                   method = "REML"
)


Eggtest2
coef_test(Eggtest2,vcov = "CR2", p_values = T)

#### No significant correlation between sample size and interaction effect size, which can be interpreted as no bias


#
#
#
#
#



full_model <- rma.mv(Interaction_Effect_Size, V, 
                     mods = ~ Macroinvertebrates + Dispersal_allowed + Mycorrhizae+ Habitat +Duration_sc,
                     random=c( ~1|Code_manuscript/ID_Code),
                     data= datred, method="ML")
eval(metafor:::.MuMIn)
rankings <- dredge(full_model)
subset(rankings, delta <= 4, recalc.weights=FALSE)
rank_df <- as.data.frame(subset(rankings, delta <= 4, recalc.weights=FALSE))

colnames(rank_df) <- c("Intercept", "Dispersal potential", "Duration","Habitat type", 
                       "Macroinvertebrate presence", "Mycorrhizal type",
                       "df", "logLik","AICc", "ΔAIC", "weight")
not_all_na <- function(rank_df) any(!is.na(rank_df))
rank_df<- cbind(Ranking= 1:nrow(rank_df), rank_df)%>% 
  mutate_if(is.numeric, round, digits = 2)%>%select(where(not_all_na))%>%
  mutate_all(funs(ifelse(is.na(.), "", .)))%>%select(-c(Intercept, logLik, AICc))

rank_df[,-c(1,7,8,9)] <- ifelse(rank_df[,-c(1,7,8,9)]=="","","+")

rank_df<- rank_df%>%
  select(Ranking, `Macroinvertebrate presence`, `Dispersal potential`,
         `Mycorrhizal type`, `Habitat type`, Duration, df, ΔAIC, weight)
sw(rankings)

sw_r <- as.data.frame(sw(rankings))
sw_r$Variables <- c("Mycorrhizal type", "Experiment duration", "Dispersal potential",
                    "Macroinvertebrate presence", "Habitat type"             
)
colnames(sw_r)[1]<- "Value"


##### Plots, ignore :)


ggplot(data=sw_r, aes(x=reorder(Variables,-Value), y=Value))+
  geom_bar(stat = "identity", width= 0.7)+
  xlab("\n\nModerator variables")+ ylab("Relative importance\n")+
  theme_classic()+
  theme(axis.title.x = element_text(face="bold", colour="black", size=16),
        axis.text.x  = element_text(vjust=0.1, size=15),
        axis.title.y = element_text(face="bold", colour="black", size=16),
        axis.text.y  = element_text(vjust=0.1, size=15),
        axis.line = element_line(colour = 'black', size = 1),
        plot.margin = unit(c(5, 5, 10, 10), "pt")) +
  scale_x_discrete(labels=function(x){sub("\\s", "\n", x)})



cols <- c("red", "gray", "cornflowerblue", "dark blue")[match(datred$Interaction_Classification, c("Synergistic", "Null", "Antagonistic", "Reversal"))]


###Remove again all previous plots (little broom)
par(mar=c(2.5,0.5,1,1), mgp=c(2,0,0), tcl=0.15)
forest(
  x=datred$Interaction_Effect_Size,
  ci.ub = datred$Interaction_CI_Upper,
  ci.lb = datred$Interaction_CI_Lower,
  order= "obs",
  slab = paste0(
    gsub(",.*$", "", 
         datred$First_author),
    " (",
    datred$Year,")_",
    datred$Subextraction),
  xlim = c(-15, 10), alim= c(-11,12), at=c(-0.5,-1,-2,-5,-10,0,0.5,1,2,5,10),
  cex = 0.65, xlab = "Interaction effect size", 
  cex.lab = 1.2, col = cols, lwd= 3,
  header= "Publication")
addpoly(rand.mv, row=-1, mlab = "", annotate = FALSE, cex= 0.8, col= "cornflowerblue", border = "cornflowerblue")
text(-15, -1, pos=4, cex=0.8, "Overall (RE model)", col= "cornflowerblue")
text(8, -1, pos=4, cex=0.7, "-0.49 [-0.87, -0.12]", col= "cornflowerblue")



par(mar= c(5,5,5,5))
funnel(rand.mv, 
       main= "Standard Error",
       xlab = "Interaction Effect Size", xlim= c(-20,10),
       cex.lab= 1.5, refline = 0,
       back = adjustcolor("cornflowerblue", alpha.f = 0.1), 
       col= "cornflowerblue", lwd= 3 
)


par(mar = c(5, 5, 5, 5))
plot(pooled_n~Interaction_Effect_Size,
     data = datred,
     xlab="Interaction effect Size",
     cex.lab=2,
     xlim=c(-10,10),
     ylab = "Pooled sample size",
     ylim=c(5,30),
     pch=21,
     col= c("cornflowerblue"),
     bg= c("lightskyblue"),
     lwd=2.5,
     cex=2,
     axes=FALSE)
abline(v = 0, lty = 2.5, lwd= 3, col= c("black",alpha= 0.2))
axis(1,lwd=3, tick=TRUE)
axis(2,lwd=3, tick=TRUE)