Map 1 MAP (a) and MSI (b) for the region of study. Red line depicts the Cerrado border, the blue line corresponds to Amazon border, purple line to Caatinga border and the black dot is placed in the Parque Nacional da Chapada dos Veadeiros (PNCV). Monthly precipitation for South America from CHIRPS V2.0.
The visual inspection of the residuals distribution shows that for MAP (Figure 1a) it is normaly distributed, with the mean near zero. For the MSI, residuals presented a tailed distribution and to deal with heterocedasticity of residuals we used the function exp to transform the data (Figure 1b).
Figure 1: Residuals inspection.
Firstly, we used the gls function from the R package nlme (Pinheiro et al., 2019) to fit a linear mixed-effects model. This model considered MAP (and MSI) as the response variable, with the factors of interest being Biome (with two categories: Amazon, Caatinga and Cerrado) and Vegetation Type (with three categories: Savanna, Gallery Forest, and Dry Forest) (Table 1).
m1 <- nlme::gls(Rainfall ~ as.factor(Biome),
data = data,
weights = nlme::varIdent(form = ~ 1 | as.factor(Biome)),
method = "ML", na.action = na.omit)
summary(m1)$tTable %>% knitr::kable()
| Value | Std.Error | t-value | p-value | |
|---|---|---|---|---|
| (Intercept) | 2268.5796 | 1.053557 | 2153.2581 | 0 |
| as.factor(Biome)Caatinga | -1577.1555 | 1.445449 | -1091.1179 | 0 |
| as.factor(Biome)Cerrado | -832.1394 | 1.516494 | -548.7260 | 0 |
| as.factor(Biome)Dry Forest | -972.7419 | 19.925029 | -48.8201 | 0 |
| as.factor(Biome)Gallery Forest | -867.0943 | 2.885341 | -300.5171 | 0 |
| as.factor(Biome)Savanna | -859.5753 | 2.113434 | -406.7198 | 0 |
Next, we performed an analysis of variance (ANOVA) using the anova function in R to assess the statistical significance of the main effects of factors.
an_m <- anova(m1)
an_m %>% data.frame %>% knitr::kable()
| numDF | F.value | p.value | |
|---|---|---|---|
| (Intercept) | 1 | 6535985.7 | 0 |
| as.factor(Biome) | 5 | 238171.8 | 0 |
ANOVA revealed statistically significant differences (Table 2) and we proceeded with post hoc tests to conduct multiple comparisons between groups. For this purpose, we performed pairwise comparisons between factors with the emmeans function from the package emmeans (Lenth, 2022) to compute marginal means and confidence intervals for each combination of factors (Table 3).
means <- emmeans::emmeans(m1, specs = list("Biome"))
means %>% data.frame %>% knitr::kable()
| Biome | emmean | SE | df | lower.CL | upper.CL |
|---|---|---|---|---|---|
| Amazon | 2268.5796 | 1.0535568 | 138088.10590 | 2265.8000 | 2271.3592 |
| Caatinga | 691.4241 | 0.9896169 | 28495.18816 | 688.8131 | 694.0352 |
| Cerrado | 1436.4402 | 1.0907662 | 59247.67705 | 1433.5624 | 1439.3180 |
| Dry Forest | 1295.8377 | 19.8971552 | 12.38576 | 1233.4788 | 1358.1965 |
| Gallery Forest | 1401.4852 | 2.6861143 | 11.10997 | 1392.8858 | 1410.0847 |
| Savanna | 1409.0042 | 1.8321083 | 12.08495 | 1403.2359 | 1414.7726 |
Subsequently, we conducted multiple comparison tests using the contrast function, also from emmeans package, with the argument method = “pairwise”. This test allowed us to identify significant differences between groups (or not) by comparing all factor levels, such as biomes and vegetation types. Adjusted p-values were considered to control the family-wise error rate by Tukey method (see contrast-methods from the package emmeans).
contrasts <- emmeans::contrast(means, method = "pairwise")
contrasts %>%
data.frame %>%
rename(Pair = "emmeans.of.Biome.contrast") %>%
rename(estimate = "emmeans.of.Biome.estimate") %>%
rename(se = "emmeans.of.Biome.SE") %>%
rename(df = "emmeans.of.Biome.df") %>%
rename(`t ratio` = "emmeans.of.Biome.t.ratio") %>%
rename(`p-value` = "emmeans.of.Biome.p.value") %>% knitr::kable()
| Pair | estimate | se | df | t ratio | p-value |
|---|---|---|---|---|---|
| Amazon - Caatinga | 1577.155475 | 1.445449 | 102390.78352 | 1091.117868 | 0.0000000 |
| Amazon - Cerrado | 832.139410 | 1.516494 | 160995.80547 | 548.725977 | 0.0000000 |
| Amazon - Dry Forest | 972.741893 | 19.925029 | 12.45527 | 48.820100 | 0.0000000 |
| Amazon - Gallery Forest | 867.094348 | 2.885341 | 14.79007 | 300.517127 | 0.0000000 |
| Amazon - Savanna | 859.575345 | 2.113434 | 21.39806 | 406.719782 | 0.0000000 |
| Caatinga - Cerrado | -745.016065 | 1.472791 | 81755.45753 | -505.853311 | 0.0000000 |
| Caatinga - Dry Forest | -604.413582 | 19.921750 | 12.44713 | -30.339382 | 0.0000000 |
| Caatinga - Gallery Forest | -710.061127 | 2.862613 | 14.33063 | -248.046519 | 0.0000000 |
| Caatinga - Savanna | -717.580130 | 2.082297 | 20.16873 | -344.609824 | 0.0000000 |
| Cerrado - Dry Forest | 140.602483 | 19.927031 | 12.46039 | 7.055867 | 0.0001242 |
| Cerrado - Gallery Forest | 34.954938 | 2.899135 | 15.07896 | 12.057025 | 0.0000000 |
| Cerrado - Savanna | 27.435935 | 2.132227 | 22.16928 | 12.867268 | 0.0000000 |
| Dry Forest - Gallery Forest | -105.647545 | 20.077649 | 12.71369 | -5.261948 | 0.0017983 |
| Dry Forest - Savanna | -113.166548 | 19.981327 | 12.59476 | -5.663615 | 0.0009654 |
| Gallery Forest - Savanna | -7.519003 | 3.251435 | 18.00246 | -2.312518 | 0.2393606 |
We found differences in MAP between dry forests and both gallery forests and savannas (Table 4). Unexpectedly, the MAP of our samples are distinct from the MAP of both Biomes, including the Cerrado where they are located (Table 4).
Furthermore, we found that the MAP for more than half of all dry forests samples are outside the bistability zone (1300-2100 mm), i.e. the range of MAP in which forests and savannas may occur, defined by Ciemer et al. (2019) for the neotropics (Figure 2). However they are near the lower limit. Despite the unexpected differences in MAP for dry forests, most of our samples are inside the bistability zone for the neotropics (1300-2100mm; Ciemer et al., 2019) (Figure 2). By the other hand, Staver et al. (2011) defined a more wide range (1000-2500 mm), but using a global model. Thus, the limits defined by Staver et al. (2011) correctly includes dry forests into the bistability zone.
Figure 2: MAP distribution for Amazon, Caatinga and Cerrado Biomes, and for our sampled vegetation types. Red lines are the bistability zone defined by Staver et al. (2011).
To compare between vegetation types and Biomes regardind MSI (i.e. seasonality lenght) (Markham, 1970), we performed the same method used to test for differences in MAP (Tables 5-8).
m2 <- nlme::gls(exp(Rainfall) ~ as.factor(Biome),
data = data1,
weights = nlme::varIdent(form = ~ 1 | as.factor(Biome)),
method = "ML", na.action = na.omit)
summary(m2)$tTable %>% knitr::kable()
| Value | Std.Error | t-value | p-value | |
|---|---|---|---|---|
| (Intercept) | 1.4846488 | 0.0004642 | 3197.9982 | 0 |
| as.factor(Biome)Caatinga | 0.3037148 | 0.0014085 | 215.6307 | 0 |
| as.factor(Biome)Cerrado | 0.2894743 | 0.0005671 | 510.4890 | 0 |
| as.factor(Biome)Dry Forest | 0.3248452 | 0.0031265 | 103.8992 | 0 |
| as.factor(Biome)Gallery Forest | 0.3262108 | 0.0012231 | 266.7144 | 0 |
| as.factor(Biome)Savanna | 0.3295360 | 0.0010822 | 304.5066 | 0 |
The ANOVA revealed significative differences among groups (Table 6).
an_m <- anova(m2)
an_m %>% data.frame %>% knitr::kable()
| numDF | F.value | p.value | |
|---|---|---|---|
| (Intercept) | 1 | 47773529.45 | 0 |
| as.factor(Biome) | 5 | 59169.47 | 0 |
The range of values of MSI is very similar among groups. However, the range of MSI values for the Amazon biome clearly differs from the other groups and do not overlap with Cerrado and Caatinga nor with our samples (Table 7; Figure 3).
means <- emmeans::emmeans(m2, specs = list("Biome"))
means %>% data.frame %>% knitr::kable()
| Biome | emmean | SE | df | lower.CL | upper.CL |
|---|---|---|---|---|---|
| Amazon | 1.484649 | 0.0004642 | 1.336860e+05 | 1.483424 | 1.485874 |
| Caatinga | 1.788364 | 0.0013298 | 2.444528e+04 | 1.784855 | 1.791872 |
| Cerrado | 1.774123 | 0.0003256 | 5.886448e+04 | 1.773264 | 1.774982 |
| Dry Forest | 1.809494 | 0.0030919 | 8.823026e+00 | 1.799037 | 1.819951 |
| Gallery Forest | 1.810860 | 0.0011315 | 6.742534e+00 | 1.806689 | 1.815030 |
| Savanna | 1.814185 | 0.0009776 | 7.525629e+00 | 1.810718 | 1.817652 |
The comparisons among groups showed differences between biomes and between biomes and all sampled vegetation types (Table 8). Dry forests, gallery forests and savannas at the PNCV do not differ in terms of length of dry season.
contrasts <- emmeans::contrast(means, method = "pairwise")
contrasts %>%
data.frame %>%
rename(Pair = "emmeans.of.Biome.contrast") %>%
rename(estimate = "emmeans.of.Biome.estimate") %>%
rename(se = "emmeans.of.Biome.SE") %>%
rename(df = "emmeans.of.Biome.df") %>%
rename(`t ratio` = "emmeans.of.Biome.t.ratio") %>%
rename(`p-value` = "emmeans.of.Biome.p.value") %>% knitr::kable()
| Pair | estimate | se | df | t ratio | p-value |
|---|---|---|---|---|---|
| Amazon - Caatinga | -0.3037148 | 0.0014085 | 3.087989e+04 | -215.630689 | 0.0000000 |
| Amazon - Cerrado | -0.2894743 | 0.0005671 | 1.921352e+05 | -510.488963 | 0.0000000 |
| Amazon - Dry Forest | -0.3248452 | 0.0031265 | 9.225340e+00 | -103.899241 | 0.0000000 |
| Amazon - Gallery Forest | -0.3262108 | 0.0012231 | 9.202904e+00 | -266.714364 | 0.0000000 |
| Amazon - Savanna | -0.3295360 | 0.0010822 | 1.130368e+01 | -304.506613 | 0.0000000 |
| Caatinga - Cerrado | 0.0142405 | 0.0013691 | 2.744074e+04 | 10.401580 | 0.0000000 |
| Caatinga - Dry Forest | -0.0211304 | 0.0033657 | 1.238590e+01 | -6.278118 | 0.0003959 |
| Caatinga - Gallery Forest | -0.0224959 | 0.0017461 | 3.841572e+01 | -12.883842 | 0.0000000 |
| Caatinga - Savanna | -0.0258212 | 0.0016504 | 6.115541e+01 | -15.645018 | 0.0000000 |
| Cerrado - Dry Forest | -0.0353709 | 0.0031090 | 9.019871e+00 | -11.377016 | 0.0000123 |
| Cerrado - Gallery Forest | -0.0367365 | 0.0011775 | 7.904815e+00 | -31.199792 | 0.0000000 |
| Cerrado - Savanna | -0.0400618 | 0.0010304 | 9.287198e+00 | -38.881078 | 0.0000000 |
| Dry Forest - Gallery Forest | -0.0013656 | 0.0032924 | 1.203992e+01 | -0.414754 | 0.9980062 |
| Dry Forest - Savanna | -0.0046908 | 0.0032427 | 1.096499e+01 | -1.446566 | 0.7012329 |
| Gallery Forest - Savanna | -0.0033253 | 0.0014953 | 1.135516e+01 | -2.223782 | 0.2982834 |
MSI do not differ among our sample plots, but they differ biomes (Table 8). However, the lenght of the dry season (i.e MSI) of the sampled vegetation types is in the range of values found for the Cerrado Biome (Table 8; Fig. 3), as they are located within the biome.
Figure 3: MSI distribution for Amazon, Caatinga and Cerrado Biomes, and for our sampled vegetation types.
Gallery forests and savannas do not differ in MAP, but dry forests differ from them. Despite differences showed in the analysis of MAP for dry forests, the values are very close compared with the range of values for the Cerrado biome. Furthermore, all vegetation types are inside the bistability zone for MAP defined by Staver et al. (2011). Vegetation types do not differ in dry season length.
Ciemer, C., Boers, N., Hirota, M., Kurths, J., Müller-Hansen, F., Oliveira, R.S., & Winkelmann, R. (2019). Higher resilience to climatic disturbances in tropical vegetation exposed to more variable rainfall. Nat. Geosci. 12, 174–179. https://doi.org/10.1038/s41561-019-0312-z
Markham, C. G. (1970). Seasonality of precipitation in the United States. Annals of the Association of American Geographers, 60(3), 593–597. https://doi.org/10.1111/J.1467-8306.1970.TB00743.X
Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core Team (2019). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-140, <URL: https://CRAN.R-project.org/package=nlme>
Lenth, R.V. (2022). emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.7.4-1. https://CRAN.R-project.org/package=emmeans
Staver, A. C., Archibald, S., & Levin, S. A. (2011). The global extent and determinants of savanna and forest as alternative biome states. Science, 230–232(334). https://doi.org/10.1126/science.1210465