Preparação

Instalar os pacotes:

install.packages("vegan")

Carrega os pacotes:

library(vegan)

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.5-2

Dados de comunidades

Vamos utilizar os dados de comunidades de plantas:

data(varespec)

Nas linhas temos os locais, nas colunas as espécies.

# Mostrar apenas as 5 primeiras linhas e colunas
varespec[1:5, 1:5]

Riqueza de espécies

A medida mais simples para medir a diversidade, é a riqueza de espécies. Podemos fazer isso somando o número de espécies de uma comunidade que apresentaram valores de abundância (biomassa, etc.) ou abundânica relativa maiores que zero.

abun_total <- colSums(varespec) # adunância total por espécie
riqueza <- sum(abun_total > 0)

Temos exatamente 44 espécies.

Algumas medidas buscam compensar o efeito da área, calculando a densidade de espécies. Por exemplo, se supomos que as áreas desse conjuto de dados possuem 55m2, então a densidade de espécies por m2 é:

area <- 55
densidade <- riqueza/area

Temos exatamente 0.8 espécies por m2.

Outros índices buscam compensar o efeito da amostragem, por exemplo, o de margalef.

amostras <- nrow(varespec)
margalef <- (riqueza - 1) / amostras
margalef

## [1] 1.791667

Outro exemplo, Índice de Menhinick.

menhinick <- riqueza / sqrt(amostras)
menhinick

## [1] 8.981462

Mas esses índices não controlam de fato o efeito da amostragem. Por exemplo Menhinick:

menhinick.seq <- numeric(amostras)
for (i in 1:amostras) {
  riq.seq <- sum(colSums(varespec[1:i, ]) > 1)
  menhinick.seq[i] <- riq.seq / sqrt(i)
}
plot(y = menhinick.seq, x = 1:amostras, type = "l", ylab = "menhinick", xlab = "sample size")

Uma solução é utilizar uma curva de acumulação de espécies.

curacum <- function(comunidade) {
  comunidade <- comunidade[, colSums(comunidade) != 0]
  primeira <- apply(comunidade, 2, function(x){min(which(x > 0))})
  frequencias <- table(factor(primeira, levels = 0:nrow(comunidade)))
  curva <- cumsum(frequencias)
  return(curva)
}

resultado <- curacum(varespec)
plot(resultado, ylab = "Species richness", xlab = "Samples", log = "x")
lines(resultado)

A posição das amostras altera a o formato do resultado.

resultado <- curacum(varespec[sample(1:amostras), ])
plot(resultado, ylab = "Species richness", xlab = "Samples", xlog = TRUE)
lines(resultado)

Por isso é importante variar a posição das amostras, já que o formato muda de acordo com isso.

rarefaction <- function(comunidade) {
  rep <- 100
  amostras_rare <- nrow(comunidade)
  resultado.rare <- matrix(ncol = amostras_rare + 1, nrow = rep)
  for (i in 1:rep) {
    resultado.rare[i, ] <- curacum(varespec[sample(1:amostras_rare), ])
  }
  return(resultado.rare)
}

resultado.rare <- rarefaction(varespec)
plot_prep <- function(rare_res) {
  media <- apply(rare_res, 2, mean)
  sd <- apply(rare_res, 2, sd)
  max_media <- media + sd
  min_media <- media - sd
  amostras_prep <- ncol(rare_res) - 1
  pol <- cbind(c(0:amostras_prep, amostras_prep:0), c(max_media, rev(min_media)))
  return(list(media, pol, max_media))
}
preps <- plot_prep(resultado.rare)
plot(0:amostras, preps[[1]], ylab = "Species richness", xlab = "Samples", ylim = c(0, max(preps[[3]])), type = "l", col = "darkblue", lwd = 2)
polygon(preps[[2]], col = rgb(0, 0, .2, .5), border = "transparent")

Comparar duas:

resultado.rare1 <- rarefaction(varespec[1:12, ])
preps <- plot_prep(resultado.rare1)
plot(0:12, preps[[1]], ylab = "Species richness", xlab = "Samples", 
     ylim = c(0, max(preps[[3]])), type = "l", col = "darkblue",
     lwd = 2)
polygon(preps[[2]], col = rgb(0, 0, .2, .5), border = "transparent")

resultado.rare2 <- rarefaction(varespec[13:24, ])
preps2 <- plot_prep(resultado.rare2)
polygon(preps2[[2]], col = rgb(.5, 0, 0, 0.5), border = "transparent")
lines(0:12, preps2[[1]], col = "red")
lines(0:12, preps[[1]], col = "darkblue")

Podemos também tentar estimar a riqueza total em uma amostra. Chao 2

# Vamos criar umas espécies raras
comunidade_rara <- varespec
comunidade_rara$rara1 <- c(rep(0, amostras - 1), 1)
comunidade_rara$rara2 <- c(rep(0, amostras - 2), 1, 1)

freqs <- apply(comunidade_rara, 2, function(x){sum(x > 0)})
unicas <- sum(freqs == 1)
duplas <- sum(freqs == 2)

# Sem espécies que ocorrem em 1 ou 2 amostras, não se pode estimar a riqueza.
chao2 <- riqueza + ((unicas ^ 2) / (2 * duplas))

Jacknife 1 ordem:

jack <- riqueza + unicas * ((amostras - 1) / amostras)

Apresentando o levantamento:

resultado.rare <- rarefaction(comunidade_rara)
preps <- plot_prep(resultado.rare)
plot(0:amostras, preps[[1]], ylab = "Species richness", xlab = "Samples",
     ylim = c(0, max(c(preps[[3]], jack))), type = "l", col = "darkblue",
     lwd = 2, bty = "l")
polygon(preps[[2]], col = rgb(0, 0, .2, .5), border = "transparent")
points(x = 24, y = jack, col = "red", pch = 18, cex = 1.5)

Índices de diversidade

Índice de Simpson (sensível a espécies dominantes):

simpson <- function(amostra) {
  pi <- amostra/sum(amostra)
  pi <- pi[pi > 0]
  D <- sum(pi ^ 2)
  return(D)
}

# Para um
simpson(varespec[1, ])

## [1] 0.1782885

# Para todos
resul_simpson <- numeric(amostras)
for (i in 1:amostras) {
  resul_simpson[i] <- simpson(varespec[i, ])
}
resul_simpson

##  [1] 0.1782885 0.2372446 0.2189888 0.2558636 0.1589201 0.1818067 0.1969030
##  [8] 0.1752328 0.4400378 0.1817193 0.1700607 0.1538473 0.1600909 0.2988451
## [15] 0.4385110 0.2611191 0.3581880 0.2173892 0.4498912 0.5038560 0.3243201
## [22] 0.4973852 0.1953731 0.1410379

O índice de Shannon (sensível em amostras pequenas):

shannon <- function(amostra) {
  pi <- amostra/sum(amostra)
  pi <- pi[pi > 0]
  H <- -sum(pi * log(pi))
  return(H)
}

# Para um
shannon(varespec[1, ])

## [1] 2.017763

# Para todos
resul_shannon <- numeric(amostras)
for (i in 1:amostras) {
  resul_shannon[i] <- shannon(varespec[i, ])
}
resul_shannon

##  [1] 2.017763 1.837101 1.834652 1.871294 2.146992 1.980967 2.014880
##  [8] 2.041299 1.238434 1.966392 2.084649 2.111599 2.121102 1.527634
## [15] 1.210139 1.634510 1.387198 1.813530 1.096279 1.189561 1.564995
## [22] 1.192426 1.921902 2.283356

Comparar o resultado do simpson com o do shannon:

plot(resul_simpson, resul_shannon, xlab = "Simpson (D)", ylab = "Shannon (H)")

Equabilidade de Shannon (ou uniformidade, ou heterogeneidade…)

J <- riqueza_amo <- numeric(amostras)
for (i in 1:amostras) {
  riqueza_amo[i] <- sum(varespec[i, ] > 0)
  J[i] <- shannon(varespec[i, ]) / riqueza_amo[i]
}
# Equabilidade 
J

##  [1] 0.06957805 0.07065772 0.07976748 0.07485176 0.08257662 0.07074881
##  [7] 0.07462517 0.07560366 0.04586792 0.06343199 0.08338597 0.07820736
## [13] 0.07070338 0.05875515 0.05261475 0.06810459 0.05548791 0.06975114
## [19] 0.05769891 0.04405783 0.08694418 0.04968444 0.08356096 0.08154842

Comparar os resultados com equabilidade e riquesa:

par(mar = c(4, 4, 1, 4))
pos <- order(resul_shannon)
plot(J[pos], riqueza_amo[pos], ylab = "Riqueza", xlab = "J")
par(new = TRUE)
plot(resul_shannon[pos], axes=F, xlab=NA, ylab=NA, type = "l", col = "red")
axis(4)
mtext("Shannon", 4, line = 2)

Série de Hill unifica todas as medidas.

hill <- function(amostra, a) {
  if (a == 1) {
    a <- 0.99999999999
  }
  pi <- amostra/sum(amostra)
  pi <- pi[pi > 0] ^ a
  Na <- sum(pi) ^ (1/(1-a))
  return(Na)
}

hill(varespec[1, ], 0) # igual a riqueza

## [1] 29

sum(varespec[1, ] > 0)

## [1] 29

hill(varespec[1, ], 1) # igual a exp(shannon)

## [1] 7.521485

exp(shannon(varespec[1, ]))

## [1] 7.521483

hill(varespec[1, ], 2) # igual a 1/D

## [1] 5.608887

1 / simpson(varespec[1, ])

## [1] 5.608887

Como apresentar:

serie_hill <- function(amostra, a) {
  alphas <- length(a)
  hill_res <- numeric(alphas)
  for (i in 1:alphas) {
    hill_res[i] <- hill(amostra, a[i])
  }
  return(hill_res)
}

a <- seq(0, 5, .1)

plot(a, serie_hill(varespec[1, ], a), xlab = "alpha", ylab = "Na", type = "l")

Comparar as 5 primeiras amostras:

plot(a, serie_hill(varespec[1, ], a), xlab = "alpha", ylab = "Na", type = "l",
     ylim = c(0, max(riqueza_amo)), lwd = 2)

for (i in 2:5) {
  lines(a, serie_hill(varespec[i, ], a), col = rainbow(5)[i], lwd = 2)
}

Outra forma, distrbuição-abundância:

plot(sort(as.numeric(varespec[1, ]), decreasing = TRUE), type = "l",
     ylab = "Abundance", xlab = "Species position", ylim = c(0, 40),
     lwd = 2)
# Comparar com as três primeiras
for (i in 1:3) {
  lines(sort(as.numeric(varespec[i, ]), decreasing = TRUE),
     ylab = "Abundance", xlab = "Species position", col = rainbow(3, alpha = .5)[i],
     lwd = 2)
}

Abundância relativa

relative <- function(x){x/sum(x)}
plot(as.numeric(sort(relative(varespec[1, ]), decreasing = TRUE)), type = "l",
     ylab = "Relative abundance", xlab = "Species position", ylim = c(0, 1),
     lwd = 2)
# Comparar com as três primeiras
for (i in 1:3) {
  lines(as.numeric(sort(relative(varespec[i, ]), decreasing = TRUE)),
     ylab = "Abundance", xlab = "Species position", col = rainbow(3, alpha = .5)[i],
     lwd = 2)
}

Exercício

Refazer as análises com outros dados. Aqui um exemplo: https://www.davidzeleny.net/anadat-r/doku.php/en:data_import_examples

Medindo a diversidade - aula 1

Bruno Vilela

9/21/2018

Preparação

Dados de comunidades

Riqueza de espécies

Índices de diversidade

Exercício