# Title: " Functions for the statistical models_ Model 2:Native vs Introduced populations"-----
# 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 2: Native vs Introduced populations --------------------------------------------

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


# LMER- Linear mixed model:

Model2_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)


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

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


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


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

  anova(Null,Flowering,Range,Herbivory,FH,RH)

}

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

Model2_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)


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

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


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

  anova(Null,Flowering,Range,Herbivory,FH,RH)

}


# GLMER Poisson-  - generalized mixed model family poisson:

Model2_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)


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

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


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


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

  anova(Null,Flowering,Range,Herbivory,FH,RH)

}




## Model 2a : 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~ Range +[1|Range/Population] + [1|Harvest date])

# LMER- Linear mixed model:

Model2a_lm<-function(data, dv){

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

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


  anova(Null,Range)

}

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

Model2a_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)

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


  anova(Null,Range)

}



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

# LMER- Linear mixed model:

Model2b_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)

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

  anova(Null,Flowering,Range)

}

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

Model2b_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)

  # Range effect:
  Range<-glmer(paste(dv, "~ Flowering+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,Range)

}

# GLMER Poisson-  - generalized mixed model family poisson:

Model2b_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)


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


  anova(Null,Flowering,Range)

}