# Title: " Functions for the statistical models_ Model 1:Environmental factors"-----
# subtitle: "Data and R code used in: Plant geographic distribution influences chemical defenses in native and introduced Plantago lanceolata populations"
# author: "Pamela Medina-van Berkum (pberkum@ice.mpg.de)"
# date: "16-January-2024"
# doi: "xxx"
# email: "pberkum@ice.mpg.de"
# license: "CC BY 4.0"

#-----------------------------------------------------------------------------

## Model 1: Environmental factors --------------------------------------------

# The model include: (y~ Flowering +PC1* Range*Herbivory +[1|Range/Population] + [1|Harvest date])


# LMER- Linear mixed model:

Model1_lm<-function(data, dv){

  # Null model:
  Null <-lmer(paste(dv, "~ 1+(1|Range:Population)+(1|Day)"),
                REML = F, data = data)
  # Flowering effect:
  Flowering <-lmer(paste(dv, "~ Flowering+(1|Range:Population)+(1|Day)"),
             REML = F, data = data)

  # PCA1 effect:
  PCA1 <-lmer(paste(dv, "~Flowering+ PCA1+(1|Range:Population)+(1|Day)"),
               REML = F, data = data)

  # Range effect:
  Range<-lmer(paste(dv, "~ Flowering+PCA1+Range+(1|Range:Population)+(1|Day)"),
              REML = F, data = data)
  # Herbivory effect:
  Herbivory<-lmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+(1|Range:Population)+(1|Day)"),
              REML = F, data = data)

  # PCA and Range interaction effect:
  PR<-lmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+PCA1:Range+(1|Range:Population)+(1|Day)"),
               REML = F, data = data)

  # PCA and Herbivory interaction effect:
  PH<-lmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+PCA1:Range+PCA1:Herbivory+(1|Range:Population)+(1|Day)"),
               REML = F, data = data)

  # Range and Herbivory interaction effect:
  RH<-lmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+PCA1:Range+PCA1:Herbivory+Range:Herbivory+(1|Range:Population)+(1|Day)"),
               REML = F, data = data)


  # Interaction effect-Full model:
  PRH<- lmer(paste(dv, "~ Flowering+PCA1*Range*Herbivory+(1|Range:Population)+(1|Day)"),
                  REML = F, data = data)

  anova(Null,Flowering,PCA1,Range,Herbivory,PR,PH,RH,PRH)

}

# GLMER - generalized mixed model family gaussian(log):

Model1_glm<-function(data, dv){

  # Null model:
  Null <-glmer(paste(dv, "~ 1+(1|Range:Population)+(1|Day)"),
               nAGQ=0,
               control=glmerControl(optimizer="bobyqa",
                                    optCtrl=list(maxfun=2e5)),
               family = gaussian(link = "log"),data=data)

  # Flowering effect:
  Flowering <-glmer(paste(dv, "~ Flowering+(1|Range:Population)+(1|Day)"),
                    nAGQ=0,
                    control=glmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                    family = gaussian(link = "log"),data=data)

  # PCA1 effect:
  PCA1 <-glmer(paste(dv, "~Flowering+ PCA1+(1|Range:Population)+(1|Day)"),
               nAGQ=0,
               control=glmerControl(optimizer="bobyqa",
                                    optCtrl=list(maxfun=2e5)),
               family = gaussian(link = "log"),data=data)

  # Range effect:
  Range<-glmer(paste(dv, "~ Flowering+PCA1+Range+(1|Range:Population)+(1|Day)"),
               nAGQ=0,
               control=glmerControl(optimizer="bobyqa",
                                    optCtrl=list(maxfun=2e5)),
               family = gaussian(link = "log"),data=data)

  # Herbivory effect:
  Herbivory<-glmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+(1|Range:Population)+(1|Day)"),
                   nAGQ=0,
                   control=glmerControl(optimizer="bobyqa",
                                        optCtrl=list(maxfun=2e5)),
                   family = gaussian(link = "log"),data=data)

  # PCA and Range interaction effect:
  PR<-glmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+PCA1:Range+(1|Range:Population)+(1|Day)"),
            nAGQ=0,
            control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e5)),
            family = gaussian(link = "log"),data=data)

  # PCA and Herbivory interaction effect:
  PH<-glmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+PCA1:Range+PCA1:Herbivory+(1|Range:Population)+(1|Day)"),
            nAGQ=0,
            control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e5)),
            family = gaussian(link = "log"),data=data)

  # Range and Herbivory interaction effect:
  RH<-glmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+PCA1:Range+PCA1:Herbivory+Range:Herbivory+(1|Range:Population)+(1|Day)"),
            nAGQ=0,
            control=glmerControl(optimizer="bobyqa",
                                 optCtrl=list(maxfun=2e5)),
            family = gaussian(link = "log"),data=data)


  # Interaction effect-Full model:
  PRH<- glmer(paste(dv, "~ Flowering+PCA1*Range*Herbivory+(1|Range:Population)+(1|Day)"),
              nAGQ=0,
              control=glmerControl(optimizer="bobyqa",
                                   optCtrl=list(maxfun=2e5)),
              family = gaussian(link = "log"),data=data)

  anova(Null,Flowering,PCA1,Range,Herbivory,PR,PH,RH,PRH)

}

# GLMER Poisson-  - generalized mixed model family poisson:

Model1_glmP<-function(data, dv){

  # Null model:
  Null <-glmer(paste(dv, "~ 1+(1|Range:Population)+(1|Day)"),
                 family="poisson", data = data)
  # Flowering effect:
  Flowering <-glmer(paste(dv, "~ Flowering+(1|Range:Population)+(1|Day)"),
              family="poisson", data = data)

  # PCA1 effect:
  PCA1 <-glmer(paste(dv, "~Flowering+ PCA1+(1|Range:Population)+(1|Day)"),
                family="poisson", data = data)

  # Range effect:
  Range<-glmer(paste(dv, "~ Flowering+PCA1+Range+(1|Range:Population)+(1|Day)"),
               family="poisson", data = data)
  # Herbivory effect:
  Herbivory<-glmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+(1|Range:Population)+(1|Day)"),
               family="poisson", data = data)

  # PCA and Range interaction effect:
  PR<-glmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+PCA1:Range+(1|Range:Population)+(1|Day)"),
                family="poisson", data = data)

  # PCA and Herbivory interaction effect:
  PH<-glmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+PCA1:Range+PCA1:Herbivory+(1|Range:Population)+(1|Day)"),
                family="poisson", data = data)

  # Range and Herbivory interaction effect:
  RH<-glmer(paste(dv, "~ Flowering+PCA1+Range+Herbivory+PCA1:Range+PCA1:Herbivory+Range:Herbivory+(1|Range:Population)+(1|Day)"),
                family="poisson", data = data)


  # Interaction effect-Full model:
  PRH<- glmer(paste(dv, "~ Flowering+PCA1*Range*Herbivory+(1|Range:Population)+(1|Day)"),
                   family="poisson", data = data)

  anova(Null,Flowering,PCA1,Range,Herbivory,PR,PH,RH,PRH)

}




## Model 1a : Without flowering as co-variante-------------------------------
# Size related traits without flowering as co-variante (Total biomass, stems and flowers). Only undamaged controls
# The model include: (y~ PC1* Range +[1|Range/Population] + [1|Harvest date])

# LMER- Linear mixed model:

Model1a_lm<-function(data, dv){

  # Null model:
  Null <-lmer(paste(dv, "~ 1+(1|Range:Population)+(1|Day)"),
                REML = F, data = data)

  # PCA1 effect:
  PCA1 <- lmer(paste(dv, "~ PCA1+(1|Range:Population)+(1|Day)"),
       REML = F, data = data)

  # Range effect:
  Range<-lmer(paste(dv, "~ PCA1+Range+(1|Range:Population)+(1|Day)"),
       REML = F, data = data)

  # Interaction effect-Full model:
  PR <- lmer(paste(dv, "~ PCA1+Range+PCA1:Range+(1|Range:Population)+(1|Day)"),
              REML = F, data = data)

  anova(Null,PCA1,Range,PR)

}

# GLMER - generalized mixed model family Gaussian(log):

Model1a_glm<-function(data, dv){

  # Null model:
  Null <-glmer(paste(dv, "~ 1+(1|Range:Population)+(1|Day)"),
               family = gaussian(link = "log"),
               nAGQ=0,
               control=glmerControl(optimizer="bobyqa",
                                    optCtrl=list(maxfun=2e5)),
               data = data)

  # PCA1 effect:
  PCA1 <- glmer(paste(dv, "~ PCA1+(1|Range:Population)+(1|Day)"),
                family = gaussian(link = "log"),
                nAGQ=0,
                control=glmerControl(optimizer="bobyqa",
                                     optCtrl=list(maxfun=2e5)),
                data = data)

  # Range effect:
  Range<-glmer(paste(dv, "~ PCA1+Range+(1|Range:Population)+(1|Day)"),
               family = gaussian(link = "log"),
               nAGQ=0,
               control=glmerControl(optimizer="bobyqa",
                                    optCtrl=list(maxfun=2e5)),
               data = data)

  # Interaction effect-Full model:
  PR <- glmer(paste(dv, "~ PCA1+Range+PCA1:Range+(1|Range:Population)+(1|Day)"),
              family = gaussian(link = "log"),
              nAGQ=0,
              control=glmerControl(optimizer="bobyqa",
                                   optCtrl=list(maxfun=2e5)),
              data = data)

  anova(Null,PCA1,Range,PR)

}




## Model 1b : With flowering as co-variante. Only undamaged controls or only damaged plants  ------------------------------------

# LMER- Linear mixed model:

Model1b_lm<-function(data, dv){

  # Null model:
  Null <-lmer(paste(dv, "~ 1+(1|Range:Population)+(1|Day)"),
                REML = F, data = data)
  # Flowering effect:
  Flowering <-lmer(paste(dv, "~ Flowering+(1|Range:Population)+(1|Day)"),
                REML = F, data = data)

  # PCA1 effect:
  PCA1 <- lmer(paste(dv, "~Flowering+ PCA1+(1|Range:Population)+(1|Day)"),
              REML = F, data = data)

  # Range effect:
  Range<-lmer(paste(dv, "~ Flowering+PCA1+Range+(1|Range:Population)+(1|Day)"),
            REML = F, data = data)

  # Interaction effect-Full model:
  PR <- lmer(paste(dv, "~ Flowering+PCA1+Range+PCA1:Range+(1|Range:Population)+(1|Day)"),
               REML = F, data = data)

  anova(Null,Flowering,PCA1,Range,PR)

}

# GLMER - generalized mixed model family Gaussian(log):

Model1b_glm<-function(data, dv){

  # Null model:
  Null <-glmer(paste(dv, "~ 1+(1|Range:Population)+(1|Day)"),
               family = gaussian(link = "log"),
               nAGQ=0,
               control=glmerControl(optimizer="bobyqa",
                                    optCtrl=list(maxfun=2e5)),
               data = data)

  # Flowering effect:
  Flowering <-glmer(paste(dv, "~ Flowering+(1|Range:Population)+(1|Day)"),
                    family = gaussian(link = "log"),
                    nAGQ=0,
                    control=glmerControl(optimizer="bobyqa",
                                         optCtrl=list(maxfun=2e5)),
                    data = data)

  # PCA1 effect:
  PCA1 <- glmer(paste(dv, "~Flowering+ PCA1+(1|Range:Population)+(1|Day)"),
                family = gaussian(link = "log"),
                nAGQ=0,
                control=glmerControl(optimizer="bobyqa",
                                     optCtrl=list(maxfun=2e5)),
                data = data)

  # Range effect:
  Range<-glmer(paste(dv, "~ Flowering+PCA1+Range+(1|Range:Population)+(1|Day)"),
               family = gaussian(link = "log"),
               nAGQ=0,
               control=glmerControl(optimizer="bobyqa",
                                    optCtrl=list(maxfun=2e5)),
               data = data)

  # Interaction effect-Full model:
  PR <- glmer(paste(dv, "~ Flowering+PCA1+Range+PCA1:Range+(1|Range:Population)+(1|Day)"),
              family = gaussian(link = "log"),
              nAGQ=0,
              control=glmerControl(optimizer="bobyqa",
                                   optCtrl=list(maxfun=2e5)),
              data = data)

  anova(Null,Flowering,PCA1,Range,PR)

}

# GLMER Poisson-  - generalized mixed model family poisson:

Model1b_glmP<-function(data, dv){

  # Null model:
  Null <-glmer(paste(dv, "~ 1+(1|Range:Population)+(1|Day)"),
                family = "poisson", data = data)
  # Flowering effect:
  Flowering <-glmer(paste(dv, "~ Flowering+(1|Range:Population)+(1|Day)"),
             family = "poisson", data = data)

  # PCA1 effect:
  PCA1 <- glmer(paste(dv, "~Flowering+ PCA1+(1|Range:Population)+(1|Day)"),
                family = "poisson", data = data)

  # Range effect:
  Range<-glmer(paste(dv, "~ Flowering+PCA1+Range+(1|Range:Population)+(1|Day)"),
              family = "poisson", data = data)

  # Interaction effect-Full model:
  PR <- glmer(paste(dv, "~ Flowering+PCA1+Range+PCA1:Range+(1|Range:Population)+(1|Day)"),
                 family = "poisson", data = data)

  anova(Null,Flowering,PCA1,Range,PR)

}