Calculate in how many communities functional diversity and compositio…

…n measures are relevantly affected by removing rare (< 4 observation) species.

inst/doc/main_analyses.R

@@ -125,53 +125,3 @@ used <- apply(commat, 1, sum)
 allobs <- apply(commat_obs, 1, sum)
 cor(used, allobs)
 
-## ------------------------------------------------------------------------
-# Specific leaf area (SLA)
-f <- function(x) mean(traitmat$sla[as.logical(x)])
-used <- apply(commat, 1, f)
-f <- function(x) mean(traitmat_all$sla[as.logical(x)])
-allobs <- apply(commat_obs, 1, f)
-cor(used, allobs)
-
-# Canopy height (CH)
-f <- function(x) mean(traitmat$ch[as.logical(x)])
-used <- apply(commat, 1, f)
-f <- function(x) mean(traitmat_all$ch[as.logical(x)])
-allobs <- apply(commat_obs, 1, f)
-cor(used, allobs)
-
-# Seed mass (SM)
-f <- function(x) mean(traitmat$sm[as.logical(x)])
-used <- apply(commat, 1, f)
-f <- function(x) mean(traitmat_all$sm[as.logical(x)])
-allobs <- apply(commat_obs, 1, f)
-cor(used, allobs)
-
-# Functional trait space (FRic)
-f <- function(x) convhulln(traitmat[as.logical(x),], "FA")$vol
-used <- apply(commat, 1, f)
-f <- function(x) convhulln(traitmat_all[as.logical(x),], "FA")$vol
-allobs <- apply(commat_obs, 1, f)
-cor(used, allobs)
-
-## Calculate distance matrix from complete trait matrix 
-dist_all <- as.matrix(dist(traitmat_all))
-
-# Mean nearest neighbout distance (mnnd)
-f <- function(x) {
-  sample.dis <- dist[as.logical(x),as.logical(x)]
-  mean(apply(sample.dis, 1, function(x) min(x[x>0])))
-}
-used <- apply(commat, 1, f)
-f <- function(x) {
-  sample.dis <- dist_all[as.logical(x),as.logical(x)]
-  mean(apply(sample.dis, 1, function(x) min(x[x>0])))
-}
-allobs <- apply(commat_obs, 1, f)
-cor(used, allobs)
-
-# Species richness SR
-used <- apply(commat, 1, sum)
-allobs <- apply(commat_obs, 1, sum)
-cor(used, allobs)
-

inst/doc/main_analyses.Rmd

@@ -339,6 +339,7 @@ In communities below 750m, the species that most often remain undetected were *B
 
 
 
+
 # Community composition and diversity along elevational gradient
 In this chapter we will calculate functional composition (i.e. single trait measures such as community mean)  and diversity (i.e. multi trait measures) from observed (`commat`) and detection-corrected meta-community (`z`). Differences between measures from observed and detection-corrected communities should be due to detection filtering.
 
@@ -570,7 +571,6 @@ f.plotFD(d, "(c) Taxonomic richness", "Number of species",
 **Fig. 3.2:** Changes in (a) functional richness (convex hull volume of the three functional dimensions’ specific leaf area, canopy height and seed mass) (b) functional packing (mean nearest neighbour distance) and (c) taxonomic diversity (number of species) along the elevational gradient. The lines represent the predictions from the generalized additive model (GAM) applied to the observed communities (dotted lines) and to the detection-corrected communities (solid line); the shaded area give the region of the bootstraped 95% confidence intervals of predictions from the GAM.
 
 
-
 ```{r, cache=TRUE}
 nsim <- 10
 
@@ -666,61 +666,73 @@ cor(used, allobs)
 ```
 We found that measures of functional composition and diversity calculated from the observations of all `r ncol(commat_obs)` species were very strongly correlated with measures calculated from the `r ncol(commat)` species with at least four observations (all r > 0.995). We therefore confident that removing species with less than four observations did not strongly bias our results.
 
+To further infer the effect of removing rare species, we compared functional composition and diversity calculated from all observations and only from species with at least four observations. Similar to how we estimated detection effects on communities, we calculated the proportion of communities removing rare species hat a relevant effect.
 
-# Effect of remooving rare species
-For all analyses we removed species that were observed at less than 4 sites. This is because it was not able to apply the hierarchical models to some of these species. To infer whether this may have impacted our results we calculated all measures of functional composition and diversity with and without missing species and infered how strong the measures were correlated. 
-
-```{r}
+```{r, cache=TRUE}
 # Specific leaf area (SLA)
-f <- function(x) mean(traitmat$sla[as.logical(x)])
-used <- apply(commat, 1, f)
-f <- function(x) mean(traitmat_all$sla[as.logical(x)])
-allobs <- apply(commat_obs, 1, f)
-cor(used, allobs)
+f.sla <- function(x) mean(traitmat$sla[as.logical(x)])
+CM.sla.obs <- apply(commat, 1, f.sla)
+f.sla <- function(x) mean(traitmat_all$sla[as.logical(x)])
+CM.sla.cor <- apply(commat_obs, 1, f.sla)
+rel <- 100 * abs(lm(CM.sla.cor ~ plantsBDM$elevation)$coef[2])
+dif.sla <- abs(CM.sla.obs - CM.sla.cor) > rel
 
 # Canopy height (CH)
-f <- function(x) mean(traitmat$ch[as.logical(x)])
-used <- apply(commat, 1, f)
-f <- function(x) mean(traitmat_all$ch[as.logical(x)])
-allobs <- apply(commat_obs, 1, f)
-cor(used, allobs)
+f.ch <- function(x) mean(traitmat$ch[as.logical(x)])
+CM.ch.obs <- apply(commat, 1, f.ch)
+f.ch <- function(x) mean(traitmat_all$ch[as.logical(x)])
+CM.ch.cor <- apply(commat_obs, 1, f.ch)
+rel <- 100 * abs(lm(CM.ch.cor ~ plantsBDM$elevation)$coef[2])
+dif.ch <- abs(CM.ch.obs - CM.ch.cor) > rel
 
 # Seed mass (SM)
-f <- function(x) mean(traitmat$sm[as.logical(x)])
-used <- apply(commat, 1, f)
-f <- function(x) mean(traitmat_all$sm[as.logical(x)])
-allobs <- apply(commat_obs, 1, f)
-cor(used, allobs)
+f.sm <- function(x) mean(traitmat$sm[as.logical(x)])
+CM.sm.obs <- apply(commat, 1, f.sm)
+f.sm <- function(x) mean(traitmat_all$sm[as.logical(x)])
+CM.sm.cor <- apply(commat_obs, 1, f.sm)
+rel <- 100 * abs(lm(CM.sm.cor ~ plantsBDM$elevation)$coef[2])
+dif.sm <- abs(CM.sm.obs - CM.sm.cor) > rel
+```
+
+Bias due to removing rare species was relevant in `r round(100*mean(dif.sla),1)`% of communities for SLA, in `r round(100*mean(dif.ch),1)`% of communities for canopy height and in `r round(100*mean(dif.sm),1)`% of communities for seed mass. 
+
 
+```{r, cache=TRUE}
 # Functional trait space (FRic)
-f <- function(x) convhulln(traitmat[as.logical(x),], "FA")$vol
-used <- apply(commat, 1, f)
-f <- function(x) convhulln(traitmat_all[as.logical(x),], "FA")$vol
-allobs <- apply(commat_obs, 1, f)
-cor(used, allobs)
+f.FRic <- function(x) convhulln(traitmat[as.logical(x),], "FA")$vol
+FRic.obs <- apply(commat, 1, f.FRic)
+f.FRic <- function(x) convhulln(traitmat_all[as.logical(x),], "FA")$vol
+FRic.cor <- apply(commat_obs, 1, f.FRic)
+rel <- 100 * abs(lm(FRic.obs ~ plantsBDM$elevation)$coef[2])
+FRic.dif <- abs(FRic.obs - FRic.cor) > rel
 
-## Calculate distance matrix from complete trait matrix 
-dist_all <- as.matrix(dist(traitmat_all))
+## Calculate distance matrix from trait matrix 
+dist <- as.matrix(dist(traitmat))
 
 # Mean nearest neighbout distance (mnnd)
-f <- function(x) {
+dist <- as.matrix(dist(traitmat))
+f.mnnd <- function(x) {
   sample.dis <- dist[as.logical(x),as.logical(x)]
   mean(apply(sample.dis, 1, function(x) min(x[x>0])))
 }
-used <- apply(commat, 1, f)
-f <- function(x) {
-  sample.dis <- dist_all[as.logical(x),as.logical(x)]
+mnnd.obs <- apply(commat, 1, f.mnnd)
+dist <- as.matrix(dist(traitmat_all))
+f.mnnd <- function(x) {
+  sample.dis <- dist[as.logical(x),as.logical(x)]
   mean(apply(sample.dis, 1, function(x) min(x[x>0])))
 }
-allobs <- apply(commat_obs, 1, f)
-cor(used, allobs)
+mnnd.cor <- apply(commat_obs, 1, f.mnnd)
+rel <- 100 * abs(lm(mnnd.obs ~ plantsBDM$elevation)$coef[2])
+mnnd.dif <- abs(mnnd.obs - mnnd.cor) > rel
 
 # Species richness SR
-used <- apply(commat, 1, sum)
-allobs <- apply(commat_obs, 1, sum)
-cor(used, allobs)
+SR.obs <- apply(commat, 1, sum)
+SR.cor <- apply(z, 1, sum)
+rel <- 100 * abs(lm(SR.obs ~ plantsBDM$elevation)$coef[2])
+SR.dif <- abs(SR.obs - SR.cor) > rel
 ```
-We found that measures of functional composition and diversity calculated from the observations of all `r ncol(commat_obs)` species were very strongly correlated with measures calculated from the `r ncol(commat)` species with at least four observations (all r > 0.995). We therefore confident that removing species with less than four observations did not strongly bias our results.
+
+Bias due to removing rare species was relevant in `r round(100*mean(FRic.dif),1)`% of communities for functional richness, and in `r round(100*mean(mnnd.dif),1)`% of communities for functional packing.  For all measures bias due to removing rare species was less strong than bias due to detection filtering. Again functional packing is the measure that reacts most sensitive to undersampling. 
 
 
 # References

inst/doc/simdat.R

@@ -3,6 +3,78 @@ library(devtools)
 install_github("TobiasRoth/detectionfilter")
 library(detectionfilter)
 
+## ---- cache=TRUE---------------------------------------------------------
+# Calculate observed community matrix
+z_obs <- apply(plantsBDM$y, c(1,2), max)
+
+# Calculate observed occupancies for all species
+psi <- apply(z_obs, 2, mean)
+
+# Mean and sd at logit-scale
+(mean.lpsi <- mean(qlogis(psi)))
+(sig.lpsi <- sd(qlogis(psi)))
+
+## ---- fig.height=4, cache=TRUE-------------------------------------------
+par(mfrow = c(1,2))
+hist(psi, breaks = 100, ylim = c(0,500), main = "", 
+     ylab = "Number of species", xlab = "Occupancy (observed)")
+hist(plogis(rnorm(1733, mean.lpsi, sig.lpsi)), 
+     breaks = 100, ylim = c(0,500), main = "",
+     ylab = "Number of species", xlab = "Occupancy (simulated)")
+
+## ---- fig.height=4, cache=TRUE-------------------------------------------
+# Choose how many species should be simulated
+nspec <- 50  
+
+# Choose the values for the four parameters to describe species' occupancies along gradient
+mean.lpsi <- 1.0        # Mean intercept of species' occupancies
+sigma.lpsi <- 1.0       # SD of species differences in intercept
+mu.beta.lpsi <- -0.3    # Mean slope of species' occupancies along gradient
+sig.beta.lpsi  <- 0.5   # SD of species differences in slopes
+
+# Simulate intercept (b0) and slope (b1) for each species
+b0 <- rnorm(nspec, mean.lpsi, sigma.lpsi) 
+b1 <- rnorm(nspec, mu.beta.lpsi, sig.beta.lpsi)
+
+# Plot species occupancies along gradient
+plot(NA, xlim = c(-3, 3), ylim = c(0, 1), 
+     xlab = "gradient (e.g. elevation)", ylab = "Occupancy")
+gradient <- seq(-3, 3, 0.1)
+for(k in 1:nspec) {
+  lpsi <- b0[k] + b1[k] * gradient 
+  points(gradient, plogis(lpsi), lwd = 0.5, ty = "l", col = "grey")
+}
+
+## ---- cache=TRUE---------------------------------------------------------
+# Number of sites in the meta-community 
+nsite <- 100
+
+# Distribution of sites along the gradient
+# Note that the choice of the distribution is not of particular importance, and 
+# we may also use a uniform distiruntion
+gradient <- rnorm(nsite, 0, 1.5) 
+
+# Simulation of the realized meta-community
+z <- matrix(NA, nsite, nspec)
+for(k in 1:nspec) {
+  lpsi <- b0[k] + b1[k] * gradient 
+  z[,k] <- rbinom(nspec, 1, plogis(lpsi))
+}
+
+## ---- cache=TRUE---------------------------------------------------------
+# Distribution of functional trait expression across species
+# (we assume that we standardized the trait values and the range 
+# of trait values covers about 6 standard deviations)
+trait <- rnorm(nspec, 0, 1.5)
+
+# Calculate community mean
+CM <- apply(z, 1, function(x) mean(trait[x==1]))
+
+# Plot CMs along gradient
+plot(gradient, CM, pch = 16, cex = 0.7, ylim = c(-1.5, 1.5))
+abline(lm(CM ~ gradient))
+abline(h = mean(trait), lty = 2)
+
 ## ------------------------------------------------------------------------
 # Effect of functional trait on species response to gradient
 mu.FTfilter.lpsi <- -1
@@ -26,6 +98,27 @@ plot(gradient, CM, pch = 16, cex = 0.7, ylim = c(-1.5, 1.5),
 abline(lm(CM~gradient))
 abline(h = mean(trait), lty = 2)
 
+## ---- fig.height=4, cache=TRUE-------------------------------------------
+# Choose the values for the four parameters to describe species' detection probability
+mean.lp <- 0.8        # Mean intercept of species' detection probability
+sigma.lp <- 1.0       # SD of species differences in intercept
+mu.beta.lp <- 1       # Mean species' phenological effect on detection probability
+sig.beta.lp  <- 1     # SD of species differences in phenological effect
+
+# Simulate intercept (b0) and slope (b1) for each species
+a0 <- rnorm(nspec, mean.lp, sigma.lp) 
+a1 <- rnorm(nspec, mu.beta.lp, sig.beta.lp)
+
+# Distribution of data of surveys
+date <- round(rnorm(nspec, 0, 10), 0)
+
+# Simulation of observed community
+y <- matrix(NA, nsite, nspec)
+for(k in 1:nspec) {
+  lp <- a0[k] + a1[k] * date 
+  y[,k] <- rbinom(nspec, 1, z[,k]*plogis(lpsi))
+}
+
 ## ---- fig.height=4-------------------------------------------------------
 # Calculate species richness
 SR_true <- apply(z, 1, sum)
@@ -47,3 +140,30 @@ points(gradient, CM_true, cex = 0.7, col = "orange")
 abline(lm(CM_true ~ gradient), col = "orange")
 
 
+## ---- fig.height=4, cache=TRUE-------------------------------------------
+# Effect of functional trait on the average detection probability of a species
+mu.FTfilter.lp <- 2
+
+# Average detection probability (intercept) also depends on functional trait 
+a0 <- rnorm(nspec, mean.lp + mu.FTfilter.lp * trait,  sigma.lp) 
+
+# Simulation of observed community (same code as above)
+y <- matrix(NA, nsite, nspec)
+for(k in 1:nspec) {
+  lp <- a0[k] + a1[k] * date 
+  y[,k] <- rbinom(nspec, 1, z[,k]*plogis(lp))
+}
+
+# Calculate community mean of trait expression
+CM_obs <- apply(y, 1, function(x) mean(trait[x==1]))
+
+# Make plots
+par(mfrow = c(1,2))
+plot(trait, plogis(a0), pch = 16, cex = 0.7, ylim = c(0,1),
+     ylab = "Avg.  detection probability", xlab = "Trait (e.g. plant height)")
+plot(gradient, CM_obs, pch = 16, cex = 0.7, ylab = "CM", ylim = c(-2,2))
+abline(lm(CM_obs ~ gradient))
+points(gradient, CM_true, cex = 0.7, col = "orange")
+abline(lm(CM_true ~ gradient), col = "orange")
+
+