Supporting information


The document includes the following analysis:

Index

  1. Access and data description
  2. Data exploration
  3. Analysis
    1. Host tree effect
    2. Distance effect
  4. Conclusion
Clicking on each of the sections takes you directly to them. Likewise, at the end of each section, there is a link to return to the index.
Back to index

Access and data description

The data consists of records of vascular epiphytes using the SVERA methodology (Sampling Vascular Epiphytes Richness and Abundance). The sampling involves recording 35 randomly selected trees at the sampling site in a cloud forests in Colombia. Both the field phase and the herbarium phase (taxonomic identifications) were carried out by myself. The taxonomic identifications were validated by experts in each of the families. The data are loaded below.


Back to index 
data <- read.csv("DataCloudForest.csv", header = T)
data %>%
  head() %>%
  kbl() %>%
  kable_styling()
id HostID DAP.cm DAP.m Height.host Bark x.host y.host Class Family Genus Species Biomass Habit
1 1 23.68226 0.2368226 16 Smooth 24.4 28 Fern ASPLENIACEAE Asplenium Asplenium rutaceum 5.2842533 E
2 1 23.68226 0.2368226 16 Smooth 24.4 28 Fern HYMENOPHYLLACEAE Trichomanes Trichomanes sp 1.4492150 E
3 1 23.68226 0.2368226 16 Smooth 24.4 28 Fern ASPLENIACEAE Asplenium Asplenium rutaceum 0.4450922 E
4 1 23.68226 0.2368226 16 Smooth 24.4 28 Fern POLYPODIACEAE Terpsichore Terpsichore sp 0.9373825 E
5 1 23.68226 0.2368226 16 Smooth 24.4 28 Fern POLYPODIACEAE Terpsichore Terpsichore sp 0.9373825 E
6 1 23.68226 0.2368226 16 Smooth 24.4 28 Fern POLYPODIACEAE Terpsichore Terpsichore sp 0.4450922 E
The database includes:
  • id: consecutive number
  • HostID: host tree unique id
  • DAP.cm y DAP.m: diámetro a la altura del pecho (1.3 m de altura) en centímetros y en metros
  • Height.host: height in meters of the host tree
  • Bark: host tree bark type: furrowed (rugosa), plated (desprendible), ridged (fisurada), and smooth (suave)
  • x.host and y.host: cartesian coordinates of the tree in 1 hectare of forest
  • Class, Family, Genus, Species: taxonomic identity of each record
  • Biomass: estimated from height and lenght of the plant sensu Benavides et al. (2019). In grams
  • Habit: sensu Zotz (2016) are epiphyte (E), hemiepiphyte (HE) and Nomadic vines (NO).
    • Species list 
    data %>%
  group_by(Class, Family, Species) %>% 
  summarise(Individuos = n()) %>%
  as.data.frame() %>%
  kbl() %>%
  kable_styling()
    ## `summarise()` has grouped output by 'Class', 'Family'. You can override using
## the `.groups` argument.
    Class Family Species Individuos
    Eudicot ARALIACEAE Schefflera trianae 2
    Eudicot ARECACEAE Wettinia kalbreyeri 2
    Eudicot BEGONIACEAE Begonia urticae 10
    Eudicot CAMPANULACEAE Burmeistera longifolia 3
    Eudicot CLUSIACEAE Clusia alata 4
    Eudicot CLUSIACEAE Clusia sp 9
    Eudicot CYCLANTHACEAE Sphaeradenia danielii 38
    Eudicot ERICACEAE Cavendishia sp2 17
    Eudicot ERICACEAE Diogenesia sp 15
    Eudicot ERICACEAE Disterigma sp1 19
    Eudicot ERICACEAE Disterigma sp3 16
    Eudicot ERICACEAE Psamisia sp1 22
    Eudicot ERICACEAE Psamisia sp2 25
    Eudicot ERICACEAE Psamisia sp3 5
    Eudicot ERICACEAE Sphyrospermum sp1 13
    Eudicot ERICACEAE Sphyrospermum sp2 1
    Eudicot ERICACEAE Sphyrospermum sp4 13
    Eudicot ERICACEAE Sphyrospermum sp5 9
    Eudicot GESNERIACEAE Besleria formosa 12
    Eudicot GESNERIACEAE Besleria sp2 3
    Eudicot GESNERIACEAE Columnea dimidiata 10
    Eudicot GESNERIACEAE Columnea fuscihirta 6
    Eudicot GESNERIACEAE Columnea sanguinea 45
    Eudicot GESNERIACEAE Columnea sp1 22
    Eudicot GESNERIACEAE Drymonia sp1 9
    Eudicot GESNERIACEAE Drymonia sp3 1
    Eudicot MELASTOMATACEAE Miconia caesia 2
    Eudicot MELASTOMATACEAE Miconia quintuplinervia 10
    Eudicot MELASTOMATACEAE Miconia sp 17
    Eudicot MELASTOMATACEAE Miconia theaezans 19
    Eudicot PIPERACEAE Peperomia acuminata 18
    Eudicot PIPERACEAE Peperomia sp1 4
    Eudicot PIPERACEAE Peperomia sp2 15
    Eudicot PIPERACEAE Peperomia tenella 14
    Eudicot PIPERACEAE Peperomia urocarpa 2
    Eudicot PIPERACEAE Piper sp1 1
    Eudicot PIPERACEAE Piper sp2 3
    Eudicot POACEAE Chusquea sp2 4
    Eudicot RUBIACEAE Nertera granadensis 2
    Eudicot URTICACEAE Bohemeria sp 2
    Eudicot URTICACEAE Pilea ceratocalyx 28
    Eudicot URTICACEAE Pilea sp1 36
    Eudicot URTICACEAE Urtica sp 6
    Eudicot VITACEAE Vitis sp 1
    Fern ASPLENIACEAE Asplenium harpeodes 27
    Fern ASPLENIACEAE Asplenium rutaceum 155
    Fern ATHYRIACEAE Diplazium hians 5
    Fern BLECHNACEAE Blechnum lherminieri 4
    Fern DRYOPTERIDACEAE Elaphoglossum erinaceum 80
    Fern DRYOPTERIDACEAE Elaphoglossum huacsaro 14
    Fern DRYOPTERIDACEAE Elaphoglossum obovatum 6
    Fern DRYOPTERIDACEAE Elaphoglossum sp1 1
    Fern DRYOPTERIDACEAE Elaphoglossum sp5 67
    Fern DRYOPTERIDACEAE Elaphoglossum sp6 15
    Fern DRYOPTERIDACEAE Megalastrum subincisum 16
    Fern DRYOPTERIDACEAE Polybotrya sp 2
    Fern HYMENOPHYLLACEAE Hymenophyllum fucoides 3
    Fern HYMENOPHYLLACEAE Hymenophyllum microcarpum 13
    Fern HYMENOPHYLLACEAE Trichomanes hymenophylloides 98
    Fern HYMENOPHYLLACEAE Trichomanes sp 84
    Fern POLYPODIACEAE Alansmia cultrata 23
    Fern POLYPODIACEAE Campyloneurum repens 34
    Fern POLYPODIACEAE Campyloneurum vulpinum 32
    Fern POLYPODIACEAE Enterosora trifurcata 2
    Fern POLYPODIACEAE Melpomene moniliformis 42
    Fern POLYPODIACEAE Pecluma camptophyllaria 28
    Fern POLYPODIACEAE Pecluma divaricata 7
    Fern POLYPODIACEAE Pleopeltis macrocarpa 1
    Fern POLYPODIACEAE Pleopeltis remota 4
    Fern POLYPODIACEAE Serpocaulon fraxinifolium 2
    Fern POLYPODIACEAE Serpocaulon levigatum 2
    Fern POLYPODIACEAE Terpsichore sp 43
    Fern VITTARIACEAE Radiovittaria remota 22
    Monocot ARACEAE Anthurium brachypodum 2
    Monocot ARACEAE Anthurium caucanum 6
    Monocot ARACEAE Anthurium cupreum 2
    Monocot ARACEAE Anthurium herthae 2
    Monocot ARACEAE Anthurium panduriforme 1
    Monocot ARACEAE Anthurium pulchellum 50
    Monocot ARACEAE Anthurium sp1 4
    Monocot ARACEAE Anthurium sp2 99
    Monocot ARACEAE Anthurium sp3 12
    Monocot ARACEAE Anthurium subcarinatum 28
    Monocot ARACEAE Anthurium tenuifolium 1
    Monocot ARACEAE Philodendron sp1 22
    Monocot ARACEAE Philodendron sp2 1
    Monocot BROMELIACEAE Guzmania triangularis 271
    Monocot BROMELIACEAE Pitcairnia auriculata 8
    Monocot BROMELIACEAE Racinaea adpressa 19
    Monocot BROMELIACEAE Racinaea michelii 2
    Monocot BROMELIACEAE Racinaea spiculosa 7
    Monocot BROMELIACEAE Racinaea tetrantha 5
    Monocot BROMELIACEAE Tillandsia sp1 17
    Monocot BROMELIACEAE Vriesea sp 10
    Monocot ORCHIDACEAE Cryptocentrum sp 16
    Monocot ORCHIDACEAE Cyrtochilum murinum 8
    Monocot ORCHIDACEAE Cyrtochilum ventilabrum 1
    Monocot ORCHIDACEAE Dichaea morrisii 2
    Monocot ORCHIDACEAE Elleanthus maculatus 2
    Monocot ORCHIDACEAE Elleanthus robustus 33
    Monocot ORCHIDACEAE Epidendrum cottoniflorum 6
    Monocot ORCHIDACEAE Epidendrum cylindrostachys 1
    Monocot ORCHIDACEAE Epidendrum geminiflorum 3
    Monocot ORCHIDACEAE Epidendrum maderoi 19
    Monocot ORCHIDACEAE Gomphicis sp 1
    Monocot ORCHIDACEAE Lepanthes discolor 7
    Monocot ORCHIDACEAE Lepanthes manabina 1
    Monocot ORCHIDACEAE Lepanthes meleagris 7
    Monocot ORCHIDACEAE Lepanthes ophelma 12
    Monocot ORCHIDACEAE Masdevallia sp 11
    Monocot ORCHIDACEAE Maxillaria acuminata 1
    Monocot ORCHIDACEAE Maxillaria aurea 13
    Monocot ORCHIDACEAE Maxillaria meridensis 18
    Monocot ORCHIDACEAE Maxillaria nubigena 11
    Monocot ORCHIDACEAE Oncidium anomalum 2
    Monocot ORCHIDACEAE Oncidium obryzatum 7
    Monocot ORCHIDACEAE Ornithosepalus sp 2
    Monocot ORCHIDACEAE Platystele sp2 13
    Monocot ORCHIDACEAE Pleurothallis canaligera 2
    Monocot ORCHIDACEAE Pleurothallis odobeniceps 18
    Monocot ORCHIDACEAE Pleurothallis pulvinaris 17
    Monocot ORCHIDACEAE Pleurothallis sp3 42
    Monocot ORCHIDACEAE Pleurothallis sp4 1
    Monocot ORCHIDACEAE Scaphyglottis sp3 12
    Monocot ORCHIDACEAE Stelis sp 9
    Monocot ORCHIDACEAE Stelis sp10 1
    Monocot ORCHIDACEAE Stelis sp4 1
    Monocot ORCHIDACEAE Stelis sp5 52
    Monocot ORCHIDACEAE Stelis sp6 8
    Monocot ORCHIDACEAE Stelis sp9 3

    Data exploration 

    data %>%
  summarise(AbFamilia = n_distinct(Family), # conteos
            AbGenero = n_distinct(Genus),
            AbEspecie = n_distinct(Species),
            Individuos = n()) %>%
  as.data.frame() %>%
  kbl() %>%
  kable_styling()
    AbFamilia AbGenero AbEspecie Individuos
    24 62 130 2279 
    data %>% 
  group_by(HostID) %>% 
  summarise(Riqueza = n_distinct(Species),
            Abundancia = n(),
            Biomasa = mean(Biomass),
            DAP = mean(DAP.cm)) %>%
  as.data.frame() -> my_data
skim(my_data)
    Data summary
    Name my_data
    Number of rows 35
    Number of columns 5
    _______________________
    Column type frequency:
    numeric 5
    ________________________
    Group variables None

    Variable type: numeric

    skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
    HostID 0 1 18.00 10.25 1.00 9.50 18.00 26.50 35.00 ‡‡‡‡‡
    Riqueza 0 1 21.86 11.37 5.00 13.00 19.00 28.50 42.00 †‡…ƒ†
    Abundancia 0 1 65.11 46.47 12.00 34.00 50.00 96.50 175.00 ‡†‚‚‚
    Biomasa 0 1 11.01 5.00 2.57 7.19 12.03 13.32 25.65 …ƒ‡‚
    DAP 0 1 23.94 12.84 7.07 13.26 22.28 32.93 56.34 ‡†ƒƒ

    Here I can observe the extreme values (P100, Q4). For example, 75% of the richness data is below 29, but trees with 43 epiphyte species were recorded. Similarly, with abundance, 75% of my abundance data (individuals on the tree) is less than 97 individuals, but there are trees with 175 individuals.

    data %>% 
  group_by(Class) %>% 
  summarise(Riqueza = n_distinct(Species),
            Abundancia = n()) %>%
  as.data.frame()  %>%
  mutate(
    Total_Abundancia = sum(Abundancia),
    Percentage_Abundancia = (Abundancia / Total_Abundancia) * 100
  ) %>%
  as.data.frame()
    ##     Class Riqueza Abundancia Total_Abundancia Percentage_Abundancia
## 1 Eudicot      44        515             2279              22.59763
## 2    Fern      29        832             2279              36.50724
## 3 Monocot      57        932             2279              40.89513

    Análisis (pruebas estadísticas) 

    sp.site <- as.data.frame.array(with(data,table(HostID,Species)))

diversityresult(sp.site, y=NULL, index="chao", method = "pooled")
    ##             chao
## pooled 151.49925
    s <- data %>%
  group_by(HostID) %>%
  summarise(SpeciesCounts = n_distinct(Species)) 
sc <- iNEXT(s$SpeciesCounts, q = 0)
sc$DataInfo
    ##   Assemblage   n S.obs SC f1 f2 f3 f4 f5 f6 f7 f8 f9 f10
## 1     site.1 765    35  1  0  0  0  0  1  0  0  1  5   1

    Normality test

    perform_shapiro <- function(column) {
  test <- shapiro.test(column)
  data.frame(
    Statistic = test$statistic,
    P_Value = test$p.value
  )
}

my_data %>%
  select(Riqueza, Abundancia, Biomasa, DAP) %>%
  map_df(perform_shapiro, .id = "Variable") %>%
  mutate(Significant = ifelse(P_Value < 0.05, "Yes", "No"))
    ##         Variable  Statistic      P_Value Significant
## W...1    Riqueza 0.92241321 0.0167422274         Yes
## W...2 Abundancia 0.88189072 0.0013349247         Yes
## W...3    Biomasa 0.94892905 0.1050487675          No
## W...4        DAP 0.94161599 0.0626808920          No

    We tested homocedasticity

    lmMod <- lm(Riqueza ~ DAP, data=my_data) # initial model
par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(lmMod)

    See top-left and bottom-left plots. The top-left plot is the residuals versus fitted values plot, while in the bottom-left, there are standardized residuals on the Y-axis. If theres absolutely no heteroscedasticity, you should see an evenly distributed and completely random scatter of points across the entire range of the X-axis and a flat red line. So, the inference here is that there is NO heteroscedasticity.

    To establish the presence or absence of heteroscedasticity, I used the NCV test (Score test for non-constant error variance. It calculates a score test for the hypothesis of constant error variance against the alternative that error variance changes with the level of the response (fitted values) or with a linear combination of predictors).

    That is, we want the variance to be constant (homoscedastic data), and its a desirable property of simple regression models because we seek to have errors with constant variance.

    library(car)
ncvTest(lmMod)  # Non-constant Variance Score Test
    ## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.031594382, Df = 1, p = 0.858921

    The test has a p-value greater than a significance level of 0.05, so we can accept the null hypothesis that the variance of the residuals is constant and infer that heteroscedasticity is not present in the errors, confirming my graphical inference.

    Since the data meet the assumptions of normality and homoscedasticity (homogeneity of variances), we carry out the tests.

    Host tree effect

    Since there is a relationship between the variables, I will perform a linear regression.

    ### TREE SIZE ####
perform_lm <- function(data, response) {
  lmMod <- lm(response ~ DAP, data = data)
  tidy(lmMod) %>%
    mutate(Significant = ifelse(p.value < 0.05, "Yes", "No"))
}

lm_results <- my_data %>%
  summarise(
    Riqueza = list(perform_lm(my_data, Riqueza)),
    Abundancia = list(perform_lm(my_data, Abundancia)),
    Biomasa = list(perform_lm(my_data, Biomasa))
  )

lm_results %>%
  pivot_longer(cols = everything(), names_to = "Response", values_to = "Results") %>%
  unnest(cols = "Results")
    ## # A tibble: 6 × 7
##   Response   term        estimate std.error statistic     p.value Significant
##   <chr>      <chr>          <dbl>     <dbl>     <dbl>       <dbl> <chr>      
## 1 Riqueza    (Intercept)    5.77     2.72       2.12  0.0415      Yes        
## 2 Riqueza    DAP            0.672    0.100      6.69  0.000000129 Yes        
## 3 Abundancia (Intercept)    4.27    12.2        0.352 0.727       No         
## 4 Abundancia DAP            2.54     0.449      5.66  0.00000260  Yes        
## 5 Biomasa    (Intercept)    6.69     1.63       4.11  0.000243    Yes        
## 6 Biomasa    DAP            0.180    0.0601     3.00  0.00506     Yes

    plotting

    # Plot for richness
a <- ggplot(my_data, aes(x = DAP, y = Riqueza)) +
  geom_smooth(method = 'lm') +
  labs(x = "DBH", y = "Epiphyte richness") +
  theme_classic()

# Plot for abundance
b <- ggplot(my_data, aes(x = DAP, y = Abundancia)) +
  geom_smooth(method = 'lm') +
  labs(x = "DBH", y = "Abundance") +
  theme_classic()

# Plot for biomass
c <- ggplot(my_data, aes(x = DAP, y = Biomasa)) +
  geom_smooth(method = 'lm') +
  labs(x = "DBH", y = "Biomass") +
  theme_classic()

plot <- plot_grid(a, b, c, labels = "AUTO")
plot

    # plot <- plot_grid(a, b, c, nrow = 1, labels = "AUTO")
# ggsave("plot.png", plot, width = 16, height = 5, units = "in", dpi = 300)

    We had the size of trees as DAP (independent variable), richness of epiphytes (dependent variable), and the type of epiphyte habit (categorical predictor variable). We use linear regression analysis to assess the relationship between the size of trees and the richness of epiphytes, while controlling for the effect of epiphyte habit.

    epiphyte_data <- data %>%
  filter(Habit != "ACC") %>%
  group_by(HostID, DAP.cm, Habit) %>%
  summarize(ri_per_ha = n_distinct(Species), .groups = 'drop')

model <- lm(ri_per_ha ~ DAP.cm + Habit, data = epiphyte_data)
summary(model)
    ## 
## Call:
## lm(formula = ri_per_ha ~ DAP.cm + Habit, data = epiphyte_data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -15.47473  -4.12504  -0.32734   2.88749  12.28746 
## 
## Coefficients:
##               Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)  10.710022   1.570971   6.8175 2.267e-09 ***
## DAP.cm        0.285365   0.051814   5.5075 5.171e-07 ***
## HabitHE     -19.191989   2.113331  -9.0814 1.311e-13 ***
## HabitNO     -14.685887   1.395547 -10.5234 2.771e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.7015 on 73 degrees of freedom
## Multiple R-squared:  0.68707,    Adjusted R-squared:  0.67421 
## F-statistic: 53.426 on 3 and 73 DF,  p-value: < 2.22e-16
    Anova(model, type="III") # Test for overall significance of the model
    ## Anova Table (Type III tests)
## 
## Response: ri_per_ha
##              Sum Sq Df F value     Pr(>F)    
## (Intercept) 1510.84  1 46.4777 2.2670e-09 ***
## DAP.cm       986.01  1 30.3327 5.1711e-07 ***
## Habit       4765.50  2 73.3003 < 2.22e-16 ***
## Residuals   2372.99 73                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    plot(model)

    Distance effect

    Overall analysis

    richness_data <- data %>%
  group_by(HostID) %>%
  summarize(overall_richness = n_distinct(Species),
            bark = first(Bark),
            dbh = first(DAP.m),
            height = first(Height.host),
            latitude = first(y.host),
            longitude = first(x.host))

comp <- vegdist(decostand(sp.site, "hell"), "euclidean")
tree.d <- vegdist(richness_data %>% 
                      select(latitude, longitude), 
                      method = "euclidian")
biotic.d  <-  dist(richness_data %>% 
                select(bark, dbh, height))
    ## Warning in dist(richness_data %>% select(bark, dbh, height)): NAs introduced by
## coercion
    ctree <- mantel(comp ~ biotic.d); ctree
    ##     mantelr       pval1       pval2       pval3   llim.2.5%  ulim.97.5% 
## 0.097587640 0.068000000 0.933000000 0.153000000 0.047977088 0.198314835
    dtree <- mantel(comp ~ tree.d); dtree
    ##    mantelr      pval1      pval2      pval3  llim.2.5% ulim.97.5% 
## 0.21810408 0.01000000 0.99100000 0.01300000 0.13532489 0.29377922

    Effect of the distance per habit

    sp.siteE <- data %>% filter(Habit == "HO")
tree.dE <- sp.siteE %>%
  group_by(HostID) %>%
  summarize(latitude = first(y.host),
            longitude = first(x.host))
tree.dE <- vegdist(tree.dE %>% 
            select(latitude, longitude),
            method = "euclidian")
sp.siteE <- as.data.frame.array(with(sp.siteE,table(HostID,Species)))
compE <- vegdist(decostand(sp.siteE, "hell"), "euclidean")
mantel(compE ~ tree.dE)
    ##    mantelr      pval1      pval2      pval3  llim.2.5% ulim.97.5% 
##         NA          1          1          1         NA         NA
    sp.siteHE <- data %>% filter(Habit == "HE")
tree.dHE <- sp.siteHE %>%
  group_by(HostID) %>%
  summarize(latitude = first(y.host),
            longitude = first(x.host))
tree.dHE <- vegdist(tree.dHE %>% 
            select(latitude, longitude),
            method = "euclidian")
sp.siteHE <- as.data.frame.array(with(sp.siteHE,table(HostID,Species)))
compHE <- vegdist(decostand(sp.siteHE, "hell"), "euclidean")
mantel(compHE ~ tree.dHE)
    ##     mantelr       pval1       pval2       pval3   llim.2.5%  ulim.97.5% 
## -0.13325754  0.72400000  0.27700000  0.51600000 -0.26663207  0.10964575
    sp.siteNO <- data %>% filter(Habit == "NO")
tree.dNO <- sp.siteNO %>%
  group_by(HostID) %>%
  summarize(latitude = first(y.host),
            longitude = first(x.host))
tree.dNO <- vegdist(tree.dNO %>% 
            select(latitude, longitude),
            method = "euclidian")
sp.siteNO <- as.data.frame.array(with(sp.siteNO,table(HostID,Species)))
compNO <- vegdist(decostand(sp.siteNO, "hell"), "euclidean")
mantel(compNO ~ tree.dNO)
    ##      mantelr        pval1        pval2        pval3    llim.2.5%   ulim.97.5% 
##  0.111429307  0.102000000  0.899000000  0.209000000 -0.037614555  0.216634031

    Conclusion

    Tree size and distance has an effect of the epiphytic distribution but it depends on the habit.