pseudoreplication issue case study elevation richness.R

# Diversity along elevation gradient - pseudoreplication issue

# Elevation 200-1000 m asl, five elevation zones
# Richness (number of species) between 1-100 species

seed <- 56654

par (mfrow = c(2,2))

# Scenario 1: 5 replicate at each zone, each replicate at different mountain (random diversity)
set.seed (seed)
elev_1 <- rep (c (200, 400, 600, 800, 1000), each = 5)
rich_1 <- runif (n = 5*5, min = 0, max = 100)

plot (rich_1 ~ elev_1, main = 'Scenario 1\n (n = 5*5, zones = 5, mountains = 5)', ylim = c(0, 110), col = rep (1:5, times = 5), pch = rep (1:5, times = 5))
cor.test (elev_1, rich_1)
summary (lm (rich_1 ~ elev_1))

# Scenario 2: 5 replicates at each zone, each replicate at THE SAME mountain (similar diversity)
set.seed (seed)
elev_2 <- rep (c (200, 400, 600, 800, 1000), each = 5)
rich_2 <- unlist (lapply (runif (n = 5, min = 0, max = 100), FUN = function (x) x + rnorm (n = 5, mean = 0, sd = 10)))

plot (rich_2 ~ elev_2, main = 'Scenario 2\n (n = 5*5, zones = 5, mountains = 1)', ylim = c(0, 110))
cor.test (elev_2, rich_2)
summary (lm (rich_2 ~ elev_2))

# Scenario 4: as scenario 3, but 100 replicates
set.seed (seed)
elev_3 <- rep (c (200, 400, 600, 800, 1000), each = 100)
rich_3 <- unlist (lapply (runif (n = 5, min = 0, max = 100), FUN = function (x) x + rnorm (n = 100, mean = 0, sd = 10)))

boxplot (rich_3 ~ elev_3, main = 'Scenario 3\n (n = 5*100, zones = 5, mountains = 1)', ylim = c(0, 110))
cor.test (elev_3, rich_3)
summary (lm (rich_3 ~ elev_3))

# Repeat 100 times each scenario and count significant P-values ----
# Scenario 1: 5 replicate at each zone, each replicate at different mountain (random diversity)
res_1 <- replicate (n = 100, simplify = FALSE, expr = {
  elev_1 <- rep (c (200, 400, 600, 800, 1000), each = 5)
  rich_1 <- runif (n = 5*5, min = 0, max = 100)
  summary (lm (rich_1 ~ elev_1))
})
P_1 <- sum (sapply (res_1, FUN = function (x) x$coefficients[2,4]) < 0.05)/100

# Scenario 2: 5 replicates at each zone, each replicate at THE SAME mountain (similar diversity)
res_2 <- replicate (n = 100, simplify = FALSE, expr = {
  elev_2 <- rep (c (200, 400, 600, 800, 1000), each = 5)
  rich_2 <- unlist (lapply (runif (n = 5, min = 0, max = 100), FUN = function (x) x + rnorm (n = 5, mean = 0, sd = 10)))
  summary (lm (rich_2 ~ elev_2))
})
P_2 <- sum (sapply (res_2, FUN = function (x) x$coefficients[2,4]) < 0.05)/100

# Scenario 3: as scenario 3, but 100 replicates
res_3 <- replicate (n = 100, simplify = FALSE, expr = {
  elev_3 <- rep (c (200, 400, 600, 800, 1000), each = 100)
  rich_3 <- unlist (lapply (runif (n = 5, min = 0, max = 100), FUN = function (x) x + rnorm (n = 100, mean = 0, sd = 10)))
  summary (lm (rich_3 ~ elev_3))
})
P_3 <- sum (sapply (res_3, FUN = function (x) x$coefficients[2,4]) < 0.05)/100

barplot (c(P_1, P_2, P_3), names.arg = c('#1', '#2', '#3'), xlab = 'Scenario', ylab = 'Proportion P < 0.05', las = 1,
         main = 'Proportion of significant results\n (out of 100 replications)', ylim = c(0,1))
abline (h = 0.05, lty = 2)