#The comments are "muted" by "#", which result that section is not executed in R.

##Data management (optional)
#Set up a working directory
#setwd("/home/zoltan/!myR/Rgraphics") # path for the folder is optional, the syntax for linux/unix  
#setwd("c:/!myR/Rgraphics") # the same as above for windows
#import external (and manage) data in R

## The packages below suggested to be loaded prior the running the codes.
library(nlme)
library(MuMIn)

## Elek-etal-aphids.xls - dataset for aphids (data = aphids) 
## Elek-etal-pathogens.xls - dataset for leaf spots (data = spots) 
## Elek-etal-pitfall.xls - dataset for carabids and spiders (data = pitfall) 
## Elek-etal-harvest.xls - dataset for wheat yield (data = harvest) 

## please convert *.xls files above  into *.csv (with comma as a separation mark) and then load datasets

aphids <- read.csv("Elek-etal-aphids.csv", sep=",", header=T)
spots <- read.csv("Elek-etal-pathogens.csv", sep=",", header=T)
pitfall <- read.csv("Elek-etal-pitfall.csv", sep=",", header=T)
harvest <- read.csv("Elek-etal-harvest.csv", sep=",", header=T)

## --------------------------- Aphids ---------------------------------------------------
## single-argument model set
aphid_m1 <- lme(sqrt(grain.aphid)~pH, random= ~1|farm, method="ML", data = aphids)
aphid_m2 <- lme(sqrt(grain.aphid)~nitrogen, random= ~1|farm, method="ML", data = aphids)
aphid_m3 <- lme(sqrt(grain.aphid)~SOC, random= ~1|farm, method="ML", data = aphids)
aphid_m4 <- lme(sqrt(grain.aphid)~config, random= ~1|farm, method="ML", data = aphids)
aphid_m5 <- lme(sqrt(grain.aphid)~landscape, random= ~1|farm, method="ML", data = aphids)

#selection of the best model
model.sel(aphid_m1, aphid_m2, aphid_m3, aphid_m4, aphid_m5)

#the best model
summary(aphid_m4)


## --------------------------- Leaf spots --------------------------------------------------
## single-argument model set
spots_m1 <- lme(sqrt(leaf.spots)~pH, random= ~1|farm, method="ML", data = spots)
spots_m2 <- lme(sqrt(leaf.spots)~nitrogen, random= ~1|farm, method="ML", data = spots)
spots_m3 <- lme(sqrt(leaf.spots)~SOC, random= ~1|farm, method="ML", data = spots)
spots_m4 <- lme(sqrt(leaf.spots)~config, random= ~1|farm, method="ML", data = spots)
spots_m5 <- lme(sqrt(leaf.spots)~landscape, random= ~1|farm, method="ML", data = spots)

#selection of the best model
model.sel(spots_m1, spots_m2, spots_m3, spots_m4, spots_m5)

#the best model
summary(spots_m5)


## --------------------------- Carabids ---------------------------------------------------
## single-argument model set
carabid_m1 <- lme(sqrt(carabids)~pH, random= ~1|farm, method="ML", data = pitfall)
carabid_m2 <- lme(sqrt(carabids)~nitrogen, random= ~1|farm, method="ML", data = pitfall)
carabid_m3 <- lme(sqrt(carabids)~SOC, random= ~1|farm, method="ML", data = pitfall)
carabid_m4 <- lme(sqrt(carabids)~config, random= ~1|farm, method="ML", data = pitfall)
carabid_m5 <- lme(sqrt(carabids)~landscape, random= ~1|farm, method="ML", data = pitfall)

#selection of the best model
model.sel(carabid_m1, carabid_m2, carabid_m3, carabid_m4, carabid_m5)

#model averaging for models with delta AICc < 2
carabid.mod.avg=dredge(lme(sqrt(carabids)~pH+nitrogen, random= ~1|farm, method="ML", data = pitfall))
summary(model.avg(carabid.mod.avg, subset = delta < 2))


## --------------------------- Spiders ---------------------------------------------------
## single-argument model set
spider_m1 <- lme(sqrt(spider)~pH, random= ~1|farm, method="ML", data = pitfall)
spider_m2 <- lme(sqrt(spider)~nitrogen, random= ~1|farm, method="ML", data = pitfall)
spider_m3 <- lme(sqrt(spider)~SOC, random= ~1|farm, method="ML", data = pitfall)
spider_m4 <- lme(sqrt(spider)~config, random= ~1|farm, method="ML", data = pitfall)
spider_m5 <- lme(sqrt(spider)~landscape, random= ~1|farm, method="ML", data = pitfall)

#selection of the best model
model.sel(spider_m1, spider_m2, spider_m3, spider_m4, spider_m5)

#model averaging for models with delta AICc < 2
spider.mod.avg=dredge(lme(sqrt(spider)~landscape+SOC+config+nitrogen, random= ~1|farm, method="ML", data = pitfall))
summary(model.avg(spider.mod.avg, subset = delta < 2))


## --------------------------- Yield - ear weight ---------------------------------------------------
## single-argument model set
ear_m1 <- lme(sqrt(ear.weight)~pH, random= ~1|farm, method="ML", data = harvest)
ear_m2 <- lme(sqrt(ear.weight)~nitrogen, random= ~1|farm, method="ML", data = harvest)
ear_m3 <- lme(sqrt(ear.weight)~SOC, random= ~1|farm, method="ML", data = harvest)
ear_m4 <- lme(sqrt(ear.weight)~config, random= ~1|farm, method="ML", data = harvest)
ear_m5 <- lme(sqrt(ear.weight)~landscape, random= ~1|farm, method="ML", data = harvest)

#selection of the best model
model.sel(ear_m1, ear_m2, ear_m3, ear_m4, ear_m5)

#the best model
summary(ear_m4)


## --------------------------- Yield - grain weight ---------------------------------------------------
## single-argument model set
grain_m1 <- lme(sqrt(grain.weight)~pH, random= ~1|farm, method="ML", data = harvest)
grain_m2 <- lme(sqrt(grain.weight)~nitrogen, random= ~1|farm, method="ML", data = harvest)
grain_m3 <- lme(sqrt(grain.weight)~SOC, random= ~1|farm, method="ML", data = harvest)
grain_m4 <- lme(sqrt(grain.weight)~config, random= ~1|farm, method="ML", data = harvest)
grain_m5 <- lme(sqrt(grain.weight)~landscape, random= ~1|farm, method="ML", data = harvest)

#selection of the best model
model.sel(grain_m1, grain_m2, grain_m3, grain_m4, grain_m5)

#the best model
summary(grain_m4)

