######################################################################
###
### Landscape diversity promotes landscape functioning in North
### America
###
### Sarah Mayor, Florian Altermatt, Tom W. Crowther, Iris Hordijk,
### Simon Landauer, Jacqueline Oehri, Merin Reji Chacko,
### Michael E. Schaepman, Bernhard Schmid, Pascal A. Niklaus (*)
###
### (*) correspondence: pascal.niklaus@uzh.ch
###
### Analysis of diversity effects
###

library(pascal) ## https://github.com/pascal-niklaus/pascal
library(asreml)

d <-read.csv("../dryad/data_CEC.csv", stringsAsFactors = T)

d <- transform(d, low_prod = S_frac + G_frac + U_frac)

######################################################################
###
### Analyses of EVI data

## ------------------- Growing-season EVI (EVIgs) --------------------

d.asr.gs <- asreml(EVIgs ~ rich + low_prod,
                   random = ~ block + comp,
                   data = d,
                   trace = F)
test.asreml(d.asr.gs)

## ---- Wald tests:
##             Df denDF  F.inc        Pr
## (Intercept)  1 327.8 2921.0 < 2.2e-16 ***
## rich         1  53.4    9.4  0.003403 **
## low_prod     1 958.4  314.3 < 2.2e-16 ***

## ---- Stratum variances:
##                  df   Variance     comp        block units!R
## comp       55.02883 0.33646538 258.0129   0.01462316       1
## block     287.95890 9.23085891   0.0000 128.09431483       1
## units!R 57934.01227 0.03293399   0.0000   0.00000000       1

## ---- Variance components:
##          component    std.error    z.ratio bound %ch
## comp    0.00117235 0.0002486155   4.715515     P   0
## block   0.07180588 0.0060056272  11.956433     P   0
## units!R 0.03293399 0.0001935052 170.196961     P   0

## ---- Dispersion:
## 0.1814773


## ---------------- Growing-season mean EVI (EVImean) ----------------

d.asr.mn <- asreml(EVImean ~ rich + low_prod,
                   random = ~ block + comp,
                   data = d,
                   trace = F)
test.asreml(d.asr.mn)

## ---- Wald tests:
##             Df denDF F.inc     Pr
## (Intercept)  1   313  1625 <2e-16 ***
## rich         1    54     8 0.0059 **
## low_prod     1  1102   272 <2e-16 ***

## ---- Stratum variances:
##              df Variance comp    block units!R
## comp       54.5    1.253  282 9.88e-03       1
## block     287.8   64.388    0 1.42e+02       1
## units!R 57934.6    0.105    0 0.00e+00       1

## ---- Variance components:
##         component std.error z.ratio bound %ch
## comp      0.00406  0.000852    4.77     P   0
## block     0.45376  0.037887   11.98     P   0
## units!R   0.10483  0.000616  170.20     P   0

## ---- Dispersion:
## 0.324

## --------------- Growing-season maximum EVI (EVImax) ---------------

d.asr.mx <- asreml(EVImax ~ rich + low_prod,
                   random = ~ block + comp,
                   data = d,
                   trace = F)
test.asreml(d.asr.mx)

## ---- Wald tests:
##             Df denDF F.inc     Pr
## (Intercept)  1   315  1362 <2e-16 ***
## rich         1    54     3  0.088 .
## low_prod     1  1312   238 <2e-16 ***

## ---- Stratum variances:
##              df Variance comp    block units!R
## comp       54.2    2.748  298 8.48e-03       1
## block     287.9  107.872    0 1.22e+02       1
## units!R 57934.9    0.197    0 0.00e+00       1

## ---- Variance components:
##         component std.error z.ratio bound %ch
## comp      0.00852   0.00177    4.82     P   0
## block     0.88429   0.07383   11.98     P   0
## units!R   0.19745   0.00116  170.20     P   0

## ---- Dispersion:
## 0.444

######################################################################
###
### Analyses of net diversity effects (NE) on EVI data

d.mix <- subset(d, rich > 1)
means_block <-aggregate(cbind(NE_EVIgs, NE_EVImean, NE_EVImax) ~ block + comp,
                    data = d.mix, FUN = mean)
mixed_NE <- aggregate(cbind(NE_EVIgs, NE_EVImean, NE_EVImax) ~ comp,
                      data = means_block, FUN = mean)

## ----------------- Overyielding of mixtures: EVIgs -----------------

t.test(mixed_NE$NE_EVIgs, alternative = "greater", mu = 0)
## t = 9, df = 49, p-value = 1e-12

## ---------------- Overyielding of mixtures: EVImean ----------------

t.test(mixed_NE$NE_EVImean, alternative = "greater", mu = 0)
## t = 9, df = 49, p-value = 7e-13

## ---------------- Overyielding of mixtures: EVImax ----------------

t.test(mixed_NE$NE_EVImax, alternative = "greater", mu = 0)
## t = 8, df = 49, p-value = 4e-10

##  ---------------- Dependency of NE(EVIgs) on LCR  -----------------

d.mix.asr.gs <- asreml(NE_EVIgs ~ rich,
                       random = ~ block + comp,
                       data = d.mix, trace = F)
test.asreml(d.mix.asr.gs)

## ---- Wald tests:
##             Df denDF F.inc      Pr
## (Intercept)  1 144.2  59.5 1.9e-12 ***
## rich         1  51.4   8.7  0.0048 **

## ---- Variance components:
##         component std.error z.ratio bound %ch
## comp     0.000353  8.41e-05     4.2     P   0
## block    0.001868  1.80e-04    10.4     P   0
## units!R  0.021080  1.51e-04   139.4     P   0

## ---- Dispersion:
## 0.145

##  --------------- Dependency of NE(EVImean) on LCR  ----------------

d.mix.asr.mn <- asreml(NE_EVImean ~ rich,
                       random = ~ block + comp,
                       data = d.mix,
                       trace = F)
test.asreml(d.mix.asr.mn)

## ---- Wald tests:
##             Df denDF F.inc      Pr
## (Intercept)  1 148.6  66.6 1.3e-13 ***
## rich         1  51.5   8.2  0.0061 **

## ---- Stratum variances:
##              df Variance comp block units!R
## comp       45.9   0.6782  534  2.59       1
## block     295.6   0.4574    0 63.65       1
## units!R 38836.5   0.0655    0  0.00       1

## ---- Variance components:
##         component std.error z.ratio bound %ch
## comp      0.00112  0.000265    4.21     P   0
## block     0.00616  0.000591   10.41     P   0
## units!R   0.06553  0.000470  139.35     P   0

## ---- Dispersion:
## 0.256

##  ---------------- Dependency of NE(EVImax) on LCR  ----------------

d.mix.asr.mx <- asreml(NE_EVImax ~ rich,
                       random = ~ block + comp,
                       data = d.mix,
                       trace = F)
test.asreml(d.mix.asr.mx)

## ---- Wald tests:
##             Df denDF F.inc      Pr
## (Intercept)  1 186.4  42.9 5.4e-10 ***
## rich         1  53.1   5.1   0.027 *

## ---- Stratum variances:
##              df Variance comp block units!R
## comp       44.2    1.108  657  2.01       1
## block     299.3    0.801    0 57.20       1
## units!R 38834.5    0.117    0  0.00       1

## ---- Variance components:
##         component std.error z.ratio bound %ch
## comp      0.00147  0.000359     4.1     P   0
## block     0.01196  0.001145    10.4     P   0
## units!R   0.11675  0.000838   139.3     P   0

## ---- Dispersion:
## 0.342

######################################################################
###
### Analysis of temporal stability of EVIgs

dicv <- subset(d, !is.na(invCV_EVIgs))
dicv <- transform(dicv, low_var = A_frac + G_frac + S_frac)

## -------------------------- CV^-1{EVIgs} ---------------------------

dicv.asr.gs <- asreml(invCV_EVIgs ~ rich + low_var,
                      random = ~ block + comp,
                      data = dicv,
                      trace = F)
test.asreml(dicv.asr.gs)

## ---- Wald tests:
##             Df denDF F.inc     Pr
## (Intercept)  1   178  3050 <2e-16 ***
## rich         1    54    10 0.0031 **
## low_var      1  2084   616 <2e-16 ***

## ---- Stratum variances:
##              df Variance comp   block units!R
## comp       54.9    362.9  224   0.197       1
## block     257.7    956.1    0 126.753       1
## units!R 52784.4     22.9    0   0.000       1

## ---- Variance components:
##         component std.error z.ratio bound %ch
## comp         1.51     0.309    4.89     P   0
## block        7.36     0.664   11.08     P   0
## units!R     22.87     0.141  162.46     P   0

## ---- Dispersion:
## 4.78

## ------------------------- CV^-1{EVImean} --------------------------

dicv.asr.mn <- asreml(invCV_EVImean ~ rich + low_var,
                      random = ~ block + comp,
                      data = dicv,
                      trace = F)
test.asreml(dicv.asr.mn)

## ---- Wald tests:
##             Df denDF F.inc     Pr
## (Intercept)  1   178  3050 <2e-16 ***
## rich         1    54    10 0.0031 **
## low_var      1  2084   616 <2e-16 ***

## ---- Stratum variances:
##              df Variance comp   block units!R
## comp       54.9    362.9  224   0.197       1
## block     257.7    956.1    0 126.753       1
## units!R 52784.4     22.9    0   0.000       1

## ---- Variance components:
##         component std.error z.ratio bound %ch
## comp         1.51     0.309    4.89     P   0
## block        7.36     0.664   11.08     P   0
## units!R     22.87     0.141  162.46     P   0

## ---- Dispersion:
## 4.78

## -------------------------- CV^-1{EVImax} --------------------------

dicv.asr.mx <- asreml(invCV_EVImax ~ rich + low_var,
                      random = ~ block + comp,
                      data = dicv,
                      trace = F)
test.asreml(dicv.asr.mx)

## ---- Wald tests:
##             Df denDF F.inc     Pr
## (Intercept)  1   202  3093 <2e-16 ***
## rich         1    54     5  0.034 *
## low_var      1  2389   547 <2e-16 ***

## ---- Stratum variances:
##              df Variance comp   block units!R
## comp       54.4    271.4  280   0.375       1
## block     256.7    796.3    0 146.100       1
## units!R 52785.9     12.8    0   0.000       1

## ---- Variance components:
##         component std.error z.ratio bound %ch
## comp        0.916    0.1857    4.93     P   0
## block       5.362    0.4811   11.15     P   0
## units!R    12.815    0.0789  162.46     P   0

## ---- Dispersion:
## 3.58

######################################################################
###
### Alternative analysis of EVIgs: Meta-analysis

r <- lapply(levels(d$block), function(b) {
    tmp <- subset(d, block == b)
    tmp2 <- aggregate(cbind(EVIgs, low_prod) ~ rich + comp, data=tmp, FUN=mean)
    tmp2.lm <- if(nrow(tmp2)>3) summary(lm(EVIgs ~ rich + low_prod, data=tmp2))
               else summary(lm(EVIgs~rich, data=tmp2))
    tmp3 <- as.data.frame(tmp2.lm$coefficients)["rich",1:2]
    data.frame(block=b, rich=tmp3[1,1], se.rich=tmp3[1,2])
})

r <- do.call("rbind", r)
r <- transform(r, block = as.factor(block), w = se.rich^-2)

r.asr <- asreml(rich ~ 1, weights= w, random=~block,
                family=asr_gaussian(dispersion=1), data=r)
test.asreml(r.asr)

## -------------------- meta-analytic effect size --------------------

summary(r.asr, coef=TRUE)$coef.fixed
##             solution std error z.ratio
## (Intercept)   0.0127   0.00224    5.64


## ---------- for reference: linear mixed model effect size ----------

d.asr.gs <- asreml(EVIgs ~ rich + low_prod,
                   random = ~ block + comp,
                   data = d, trace = F)
summary(d.asr.gs, coef=TRUE)$coef.fixed["rich", ]
##             solution std error z.ratio
## rich          0.0155   0.00503    3.08