### R script for multiple comparison procedures of species-specific browsing rates: 

require(betareg)
require(multcomp)
require(lmtest)
require(AICcmodavg)
data<-read.table("C:\\data\\data.txt",header=T)
attach(data)
names(data)
summary(data)
Sp<-factor(Sp)
Stand<-factor(Stand)
Year<-factor(Year)
Tr<-factor(Tr)

### comparison of null model (overall browsing probability equals on all species) vs species model (msp, browsing probability varies among species);
### error structure considers the sampling design, with stand and year random effects; 

### A1: multiple pair-wise comparison of proportion of browsed plants, generalized linear mixed models with lmer4: 
require(lme4)
require(arm)
null <- glmer(BrY ~ 1 + (1 | Stand/Year/Tr), data = data, family = binomial)
msp <- glmer(BrY ~ -1 + Sp + (1 | Stand/Year/Tr), data = data, family = binomial)
AICc(null)
AICc(msp)
lrtest(null, msp)
summary(msp)
mctukey <- glht(msp, linfct = mcp(Sp = "Tukey"))
summary(mctukey)
### palatable (Qi and Pl) and unpalatable species (Au, Cm, Cs, and Ea) significantly differ; Qi and Pl do not significantly differ in the proportion of browsed plants.
### extract information and plot:
mctukey.cld <- cld(mctukey)
### use sufficiently large upper margin
old.par <- par( mai=c(1,1,1.25,1))
### plot
plot(mctukey.cld)
### figure 1 shows that there are two groups of species, palatable (Pl and Qi) and unpalatable (Au, Cm, Cs, and Ea); 
### get confidence intervals
ci.glht <- confint(mctukey)
### plot them
plot(ci.glht)


### A2: multiple pair-wise comparison of logit transformed proportion of browsed shoots, linear mixed models with lmer4: 
require(lme4)
require(arm)
null <- lmer(logit(Br) ~ 1 + (1 | Stand/Year/Tr), data = data)
msp <- lmer(logit(Br) ~ -1 + Sp + (1 | Stand/Year/Tr), data = data)
AICc(null)
AICc(msp)
lrtest(null, msp)
summary(msp)
mctukey <- glht(msp, linfct = mcp(Sp = "Tukey"))
summary(mctukey)
### Pl and Qi significantly differ in the proportion of browsed shoots.
### extract information and plot:
mctukey.cld <- cld(mctukey)
### use sufficiently large upper margin
old.par <- par( mai=c(1,1,1.25,1))
### plot
plot(mctukey.cld)
### figure 2 confirms that there are two groups of species, palatable (Pl and Qi) and unpalatable (Au, Cm, Cs, and Ea);
### Additionally, the two palatable species significantly differs in the proportion of browsed shoots. 
### get confidence intervals
ci.glht <- confint(mctukey)
### plot them
plot(ci.glht)


### A2b: multiple comparison of proportion of browsed shoots, untransformed proportion data (0<y<1) analysed with beta regression:  
null<-betareg(Br~1|-1+Stand+Year,data=data)
msp<-betareg(Br~-1+Sp|-1+Stand+Year,data=data)
AICc(null)
AICc(msp)
lrtest(null, msp)
summary(msp)
mctukey <- glht(msp, linfct = c("SpCm - SpAu == 0", "SpCs - SpAu == 0", "SpEa - SpAu == 0", "SpPl - SpAu == 0", "SpQi - SpAu == 0", "SpCs - SpCm == 0", "SpEa - SpCm == 0", "SpPl - SpCm == 0", "SpQi - SpCm == 0", "SpEa - SpCs == 0", "SpPl - SpCs == 0", "SpQi - SpCs == 0", "SpPl - SpEa == 0", "SpQi - SpEa == 0", "SpQi - SpPl == 0")) 
summary(mctukey)
### all results of A2 confirmed, Pl and Qi significantly differ in the proportion of browsed shoots.

### A2c: we also performed simple ANOVA (no random effects) and Tukey test for multiple pair-wise comparison of logit transformed data of the proportion of browsed shoots: 
maov <- aov(logit(Br) ~ Sp, data = data)
mctukey <- glht(maov, linfct = mcp(Sp = "Tukey"))
summary(maov)
summary(mctukey)
### results of A2 and A2b confirmed, Pl and Qi significantly differ in the proportion of browsed shoots per plant.