diff --git a/inst/doc/main_analyses.R b/inst/doc/main_analyses.R
index d322b95..ca81b93 100644
--- a/inst/doc/main_analyses.R
+++ b/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)
-
diff --git a/inst/doc/main_analyses.Rmd b/inst/doc/main_analyses.Rmd
index d948c38..e723243 100644
--- a/inst/doc/main_analyses.Rmd
+++ b/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
diff --git a/inst/doc/main_analyses.pdf b/inst/doc/main_analyses.pdf
index 0d1e5c5..c218347 100644
Binary files a/inst/doc/main_analyses.pdf and b/inst/doc/main_analyses.pdf differ
diff --git a/inst/doc/simdat.R b/inst/doc/simdat.R
index d5c331c..d82a212 100644
--- a/inst/doc/simdat.R
+++ b/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")
+
+
diff --git a/inst/doc/simdat.pdf b/inst/doc/simdat.pdf
index 93a6dc2..ef9e688 100644
Binary files a/inst/doc/simdat.pdf and b/inst/doc/simdat.pdf differ
diff --git a/inst/doc/workflow.R b/inst/doc/workflow.R
index d78fe39..cd86328 100644
--- a/inst/doc/workflow.R
+++ b/inst/doc/workflow.R
@@ -1,4 +1,109 @@
+## ---- echo=TRUE, message=FALSE, warning=FALSE, cache=TRUE----------------
+library(devtools)
+install_github("TobiasRoth/detectionfilter")
+library(detectionfilter)
+
+## ---- cache=TRUE---------------------------------------------------------
+set.seed(1234)
+dat <- simcom(mu.FTfilter.lpsi = -0.5, mean.psi = 0.8, 
+              mu.FTfilter.lp = 1, mean.p = 0.8,
+              nsite = 200, nspec = 100, nrep = 2)
+
+## ---- cache=TRUE---------------------------------------------------------
+# Calculate observed meta-community (i.e. merge the observations from the two visits)
+commat_obs <- apply(dat$y, c(1,2), max)
+
+# Show species with few observations
+sort(apply(commat_obs, 2, sum))[1:10]
+
+## ---- message=FALSE, warning=FALSE, cache=TRUE---------------------------
+# Package to fit hierarchical models of occurrence and abundance data
+library(unmarked)
+
+# Prepare detection-corrected meta-community matrix
+z <- array(NA, dim = c(dim(dat$y)[1], dim(dat$y)[2]))
+
+# Apply hierarchical model to all species seperate
+for(k in 1:dim(dat$y)[2]) {
+  d <- unmarked::unmarkedFrameOccu(y = dat$y[,k,], obsCovs = list(date = dat$date),
+                         siteCovs = data.frame(gradient = dat$gradient))
+  res <- unmarked::occu(~ date ~ gradient, data = d, se = TRUE)
+  z[,k] <- unmarked::bup(unmarked::ranef(res), stat = "mode")
+}
+
+## ----fig1, cache=TRUE, fig.height=5, fig.width= 4------------------------
+# Function to calculate mean trait expression of species in a community (CM)
+CM <- function(x) mean(dat$traitmat[x])
+
+# Calculate CMs from observed meta-community matrix, from detection-corrected 
+# meta-community matrix and from true mata-community matrix
+CM_obs <- apply(commat_obs==1, 1, CM)
+CM_cor <- apply(z>0.5, 1, CM)
+CM_true <- apply(dat$z_true==1, 1, CM)
+
+# Plot CM along gradient
+plot(dat$gradient, CM_obs, ylim = c(-1, 1), xlim = c(-5, 5), axes=F, xlab = "", ylab = "", pch = 16, cex = 0.7)
+points(dat$gradient, CM_true, cex = 0.7, col = "orange")
+points(dat$gradient, CM_cor, pch = "+", cex = 0.7, col = "red")
+a <- coef(lm(CM_obs ~ dat$gradient))
+curve(a[1] + a[2]*x, -5, 5, add = TRUE)
+a <- coef(lm(CM_true ~ dat$gradient))
+curve(a[1] + a[2]*x, -5, 5, add = TRUE, col = "orange")
+a <- coef(lm(CM_cor ~ dat$gradient))
+curve(a[1] + a[2]*x, -5, 5, add = TRUE, col = "red")
+polygon(c(-10,10,10,-10), c(-10, -10, -1, -1), col = "white", border = "white")
+axis(side=1, seq(-5, 5, 2.5), pos = -1, cex = 0.7)
+axis(side = 2, las = 1, pos = -5, cex = 0.7)
+mtext(text = "Gradient", side = 1, line = 2, cex = 1)
+mtext(text = "Community mean (CM)", side = 2, line = 3, cex = 1)
+
+## ---- cache=TRUE---------------------------------------------------------
+set.seed(1234)
+dat <- simcom(mu.FTfilter.lpsi = -0.5, mean.psi = 0.5, 
+              mu.FTfilter.lp = 2, mean.p = 0.5,
+              nsite = 200, nspec = 100, nrep = 2)
+
+# Observed community matrix
+commat_obs <- apply(dat$y, c(1,2), max)
+
+# Number of observed species
+ncol(commat_obs)
+
+## ---- cache=TRUE---------------------------------------------------------
+sum(apply(dat$z_true, 2, sum)>0)
+
 ## ------------------------------------------------------------------------
 # Proportion of  species observed at least once but in less than 10% of sites
 mean(apply(commat_obs, 2, sum) < 20)
 
+## ----fig2, message=FALSE, warning=FALSE, cache=TRUE, fig.height=5, fig.width= 4----
+# Estimate detection-corrected meta-community matrix
+z <- array(NA, dim = c(dim(dat$y)[1], dim(dat$y)[2]))
+for(k in 1:dim(dat$y)[2]) {
+  d <- unmarked::unmarkedFrameOccu(y = dat$y[,k,], obsCovs = list(date = dat$date),
+                         siteCovs = data.frame(gradient = dat$gradient))
+  res <- unmarked::occu(~ date ~ gradient, data = d, se = TRUE)
+  z[,k] <- unmarked::bup(unmarked::ranef(res), stat = "mode")
+}
+
+# Calculate CMs
+CM_obs <- apply(commat_obs==1, 1, CM)
+CM_cor <- apply(z>0.5, 1, CM)
+CM_true <- apply(dat$z_true==1, 1, CM)
+
+# Plot CMs along gradient
+plot(dat$gradient, CM_obs, ylim = c(-1.5, 1.5), xlim = c(-5, 5), axes=F, xlab = "", ylab = "", pch = 16, cex = 0.7)
+points(dat$gradient, CM_true, cex = 0.7, col = "orange")
+points(dat$gradient, CM_cor, pch = "+", cex = 0.7, col = "red")
+a <- coef(lm(CM_obs ~ dat$gradient))
+curve(a[1] + a[2]*x, -5, 5, add = TRUE)
+a <- coef(lm(CM_true ~ dat$gradient))
+curve(a[1] + a[2]*x, -5, 5, add = TRUE, col = "orange")
+a <- coef(lm(CM_cor ~ dat$gradient))
+curve(a[1] + a[2]*x, -5, 5, add = TRUE, col = "red")
+polygon(c(-10,10,10,-10), c(-10, -10, -1.5, -1.5), col = "white", border = "white")
+axis(side=1, seq(-5, 5, 2.5), pos = -1.5, cex = 0.7)
+axis(side = 2, las = 1, pos = -5, cex = 0.7)
+mtext(text = "Gradient", side = 1, line = 2, cex = 1)
+mtext(text = "Community mean (CM)", side = 2, line = 3, cex = 1)
+
diff --git a/inst/doc/workflow.pdf b/inst/doc/workflow.pdf
index 0ec1fea..ff1da56 100644
Binary files a/inst/doc/workflow.pdf and b/inst/doc/workflow.pdf differ
diff --git a/vignettes/deteccor.RDATA b/vignettes/deteccor.RDATA
index 1cc4c4f..77d2ca6 100644
Binary files a/vignettes/deteccor.RDATA and b/vignettes/deteccor.RDATA differ
diff --git a/vignettes/main_analyses.Rmd b/vignettes/main_analyses.Rmd
index 5a3a310..e723243 100644
--- a/vignettes/main_analyses.Rmd
+++ b/vignettes/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 inferred 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
diff --git a/vignettes/main_analyses_files/figure-latex/fig1-1.pdf b/vignettes/main_analyses_files/figure-latex/fig1-1.pdf
index 341b208..bd131f8 100644
Binary files a/vignettes/main_analyses_files/figure-latex/fig1-1.pdf and b/vignettes/main_analyses_files/figure-latex/fig1-1.pdf differ
diff --git a/vignettes/main_analyses_files/figure-latex/fig2.1-1.pdf b/vignettes/main_analyses_files/figure-latex/fig2.1-1.pdf
index 5e989d4..432eddb 100644
Binary files a/vignettes/main_analyses_files/figure-latex/fig2.1-1.pdf and b/vignettes/main_analyses_files/figure-latex/fig2.1-1.pdf differ
diff --git a/vignettes/main_analyses_files/figure-latex/fig2.2-1.pdf b/vignettes/main_analyses_files/figure-latex/fig2.2-1.pdf
index 5ebbfb2..d299a81 100644
Binary files a/vignettes/main_analyses_files/figure-latex/fig2.2-1.pdf and b/vignettes/main_analyses_files/figure-latex/fig2.2-1.pdf differ
diff --git a/vignettes/main_analyses_files/figure-latex/fig3.1-1.pdf b/vignettes/main_analyses_files/figure-latex/fig3.1-1.pdf
index b5e059a..1d1326c 100644
Binary files a/vignettes/main_analyses_files/figure-latex/fig3.1-1.pdf and b/vignettes/main_analyses_files/figure-latex/fig3.1-1.pdf differ
diff --git a/vignettes/main_analyses_files/figure-latex/fig3.2-1.pdf b/vignettes/main_analyses_files/figure-latex/fig3.2-1.pdf
index 22e3919..a3f7dcf 100644
Binary files a/vignettes/main_analyses_files/figure-latex/fig3.2-1.pdf and b/vignettes/main_analyses_files/figure-latex/fig3.2-1.pdf differ
diff --git a/vignettes/main_analyses_files/figure-latex/fig3.3-1.pdf b/vignettes/main_analyses_files/figure-latex/fig3.3-1.pdf
index e5e9b92..8321617 100644
Binary files a/vignettes/main_analyses_files/figure-latex/fig3.3-1.pdf and b/vignettes/main_analyses_files/figure-latex/fig3.3-1.pdf differ
diff --git a/vignettes/main_analyses_files/figure-latex/fig4-1.pdf b/vignettes/main_analyses_files/figure-latex/fig4-1.pdf
new file mode 100644
index 0000000..15a71fa
Binary files /dev/null and b/vignettes/main_analyses_files/figure-latex/fig4-1.pdf differ
diff --git a/vignettes/simdat_files/figure-latex/unnamed-chunk-10-1.pdf b/vignettes/simdat_files/figure-latex/unnamed-chunk-10-1.pdf
index 6c22bb4..e32020d 100644
Binary files a/vignettes/simdat_files/figure-latex/unnamed-chunk-10-1.pdf and b/vignettes/simdat_files/figure-latex/unnamed-chunk-10-1.pdf differ
diff --git a/vignettes/simdat_files/figure-latex/unnamed-chunk-3-1.pdf b/vignettes/simdat_files/figure-latex/unnamed-chunk-3-1.pdf
index a62e3cf..b7c029b 100644
Binary files a/vignettes/simdat_files/figure-latex/unnamed-chunk-3-1.pdf and b/vignettes/simdat_files/figure-latex/unnamed-chunk-3-1.pdf differ
diff --git a/vignettes/simdat_files/figure-latex/unnamed-chunk-4-1.pdf b/vignettes/simdat_files/figure-latex/unnamed-chunk-4-1.pdf
index fb7968f..f20b184 100644
Binary files a/vignettes/simdat_files/figure-latex/unnamed-chunk-4-1.pdf and b/vignettes/simdat_files/figure-latex/unnamed-chunk-4-1.pdf differ
diff --git a/vignettes/simdat_files/figure-latex/unnamed-chunk-6-1.pdf b/vignettes/simdat_files/figure-latex/unnamed-chunk-6-1.pdf
index 8745309..6b7fd1e 100644
Binary files a/vignettes/simdat_files/figure-latex/unnamed-chunk-6-1.pdf and b/vignettes/simdat_files/figure-latex/unnamed-chunk-6-1.pdf differ
diff --git a/vignettes/simdat_files/figure-latex/unnamed-chunk-7-1.pdf b/vignettes/simdat_files/figure-latex/unnamed-chunk-7-1.pdf
index 3c19e6c..be70a73 100644
Binary files a/vignettes/simdat_files/figure-latex/unnamed-chunk-7-1.pdf and b/vignettes/simdat_files/figure-latex/unnamed-chunk-7-1.pdf differ
diff --git a/vignettes/simdat_files/figure-latex/unnamed-chunk-9-1.pdf b/vignettes/simdat_files/figure-latex/unnamed-chunk-9-1.pdf
index 85f67f8..1a7037b 100644
Binary files a/vignettes/simdat_files/figure-latex/unnamed-chunk-9-1.pdf and b/vignettes/simdat_files/figure-latex/unnamed-chunk-9-1.pdf differ
diff --git a/vignettes/workflow_files/figure-latex/fig1-1.pdf b/vignettes/workflow_files/figure-latex/fig1-1.pdf
index 16909ef..f7d01ba 100644
Binary files a/vignettes/workflow_files/figure-latex/fig1-1.pdf and b/vignettes/workflow_files/figure-latex/fig1-1.pdf differ
diff --git a/vignettes/workflow_files/figure-latex/fig2-1.pdf b/vignettes/workflow_files/figure-latex/fig2-1.pdf
index 8a4ef4c..b3e2732 100644
Binary files a/vignettes/workflow_files/figure-latex/fig2-1.pdf and b/vignettes/workflow_files/figure-latex/fig2-1.pdf differ