### R script for Qi data (YfC=T): 

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

### multivariate analysis with bacward model selection based on AIC: 
full<-betareg(Br~Du+H+YfC+WB+Dp|-1+Stand+Year+Tr,data=data)
M1<-betareg(Br~Du+H+YfC+WB|-1+Stand+Year+Tr,data=data)
M2<-betareg(Br~Du+H+YfC+Dp|-1+Stand+Year+Tr,data=data)
M3<-betareg(Br~Du+H+WB+Dp|-1+Stand+Year+Tr,data=data)
M4<-betareg(Br~Du+YfC+WB+Dp|-1+Stand+Year+Tr,data=data)
M5<-betareg(Br~H+YfC+WB+Dp|-1+Stand+Year+Tr,data=data)
summary(full)
summary(M1)
summary(M2)
summary(M3)
summary(M4)
summary(M5)
AICc(M1)
AICc(M2)
AICc(M3)
AICc(M4)
AICc(M5)
AICc(full)
lrtest(full, M2)
summary(M2)

### results and parameter estimates:
Call:
betareg(formula = Br ~ Du + H + YfC + Dp | -1 + Stand + Year + Tr, data = data)

Standardized weighted residuals 2:
    Min      1Q  Median      3Q     Max 
-3.2506 -0.4704  0.1499  0.6177  2.8216 

Coefficients (mean model with logit link):
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.255824   0.370274   8.793  < 2e-16 ***
Du          -0.037369   0.010042  -3.721 0.000198 ***
H           -0.008176   0.001483  -5.511 3.56e-08 ***
YfC         -0.419508   0.081923  -5.121 3.04e-07 ***
Dp          -0.075332   0.011724  -6.426 1.31e-10 ***

Phi coefficients (precision model with log link):
       Estimate Std. Error z value Pr(>|z|)    
StandA  0.02996    0.23985   0.125   0.9006    
StandB -0.23313    0.22506  -1.036   0.3003    
StandC -0.46408    0.26567  -1.747   0.0807 .  
YearB   1.09216    0.20045   5.449 5.08e-08 ***
YearC   0.31804    0.20810   1.528   0.1264    
YearD   1.06107    0.22284   4.762 1.92e-06 ***
Tr     -0.06662    0.05323  -1.252   0.2107    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 306.2 on 12 Df
Pseudo R-squared: 0.3341
Number of iterations: 24 (BFGS) + 3 (Fisher scoring)