Mestranda em Ecologia, Instituto de Biociências, USP.

Título do meu projeto: “Padrões de diversidade de estratégias ecológicas das espécies arbóreas da Restinga Alta, Ilha do Cardoso/SP”.

exec

Diversidade funcional corresponde ao tipo, amplitude e abundância de atributos (características) presentes em uma comunidade que influenciam o funcionamento dos ecossistemas (1).

A minha proposta é criar uma função para medir a diversidade funcional de uma comunidade. Essa função calcularia os índices de riqueza, equitatividade e divergência (2) e o índice FD (functional diversity)(3). A função retornaria uma tabela com esses índices e a significância de cada um, calculado através de permutação.

Referências: 1)Díaz, S., Lavorel, S., Chapin III, F.S., Tecco, P.A., Gurvich,D.E. & Grigulis, K. 2007. Functional Diversity – at the Crossroads between Ecosystem Functioning and Environmental Filters. In: Canadell, J.G., Pataki, D., Pitelka, L. (eds). Terrestrial Ecosystems in a Changing World. The IGBP Series, Springer-Verlag, Berlin Heidelber. 2)Mason, N. W. H., Mouillot, D., Lee, W. G. & Wilson, J. B. 2005. Functional richness, functional evenness and functional divergence: the primary components of functional diversity. Oikos 111: 112-118. 3)Petchey, O.L. & Gaston,K.J. 2002. Functional diversity (FD), species richness, and community composition. Ecol. Lett. 5:402-11.

A proposta alternativa seria criar uma função de decomposição. A partir dos dados brutos, essa função calcularia a(s) taxa(s) de decomposição (k) e aplicaria um teste de estatística (teste-t ou Anova ou Monte Carlo - o usuário poderia escolher). A função também retornaria um gráfico dos ks.

As duas ideias são boas e bem dimensionadas para o trabalho final. Mas é preciso definir claramente entrada e saída da sua função (seções 'arguments' e 'value' da página de ajuda). Só comece depois de fazer isto.

Para o plano A, explique também o que o teste pretende comparar, qual a hipótese nula, e como a randomização irá reproduzir este cenário nulo. Para o plano B, não está claro o que é testado e o é o gráfico.

Os índices propostos por Mason et al. (2005) são baseados em dados univariados. Ao contrário dos índices propostos por Villéger et al. (2008), que são baseados em dados multivariados. Dessa forma, obtei pelo uso de dados multivariados que permitem uma melhor compreensão sobre as estrutura das comunidades.

diversidade.funcional                           package:world.domination                              R Documentation

                             Calcula três componentes da diversidade funcional de uma comunidade

Description
A função calcula três índices de diversidade funcional para uma comunidade: riqueza, equitatividade e divergência. E a significância dos índices de equitatividade e divergência funcional. 

Usage
diversidade.funcional (x, y, nsim=1000, salvar.graficos=FALSE)

Arguments
x                  matriz numérica contínua de espécies por atributos 
y                  matriz numérica de abundância das espécies na comunidade.
nsim               número de simulações para a criar o cenário nulo
salvar.graficos    salva os gráficos gerados em um arquivo pdf

Details
Os valores dos atributos da amostra (x) são padronizados (média=0 e variância=1) e convertidos para uma matriz de distâncias. 
Para o índice de riqueza funcional (FD), calcula-se o comprimento total dos braços ligando todas as espécies em um dendrograma, excluindo os segmentos desnecessários de raiz da árvore. 
Para o índice de equitatividade funcional (FEve), calcula-se as distâncias entre pares de espécies dividida pela soma da abundância ponderada dos pares de espécies(EWi). Cada EWi é dividido pela somatória total dos valores de EWi e são comparados com os valores de 1/(nº de espécies - 1) para a obtenção do mínimo. Os mínimos são subtraídos e divididos pelo termo 1/(nº de espécies - 1) e 1 - 1/(nº de espécies - 1), respectivamente. Esses valores são somados para gerar o índice em questão.
Para o índice de divergência funcional (FDiv), calcula-se a distância Euclidiana e média entre as espécies e o centro de gravidade de um polígono formado no espaço funcional (convex hull). Os valores da distância Euclidiana são subtraídos da distância média em módulo e em não-módulo (desvio modular e desvio não-modular, respectivamente). Cada desvio é multiplicado pela abundância ponderada das espécies da comunidade.  Em seguida, calcula-se a somatória dos desvios não modulares dividido pela somatória dos desvios modulares. 

Para o calculo da significância dos índices de equitatividade e divergência funcional, comparou-se os valores observados com os valores dos índices obtidos a partir de um cenário nulo.  A hipótese nula testada é de que a comunidade não tem um padrão de estruturação e os valores observados são ao acaso. 
Para o calculo do cenário nulo, as espécies da matriz x (linhas) foram aleatorizadas com reposição, mantendo-se constante a identidade dos atributos. Sendo calculados novamente os índices para cada nova comunidade gerada.

Value
Um “data frame” é exibido na tela com os seguintes vetores:
Índices: valores da riqueza (FD), equitatividade (FEve) e divergência funcional (FDiv). São valores positivos, sendo que 
         o FD pode ir até infinito. Já os índices FEve e FDiv estão restritos entre 0 e 1;
p-valor: valor da significância dos índices FEve e FDiv;
IC 2.5 e 97.5%: percentis do cenário nulo correspondentes a cada índice (FEve e FDiv).
E cinco gráficos são gerados em um mesmo dispositivo gráfico.

Warnings
As matrizes x e y devem ter a primeira dimensão de comprimentos iguais. E não é permitido NAs em nenhuma das matrizes.
O processamento da função pode levar vários minutos e é aconselhável usar o programa do R sem nenhum editor.
Para matrizes (x) com mais de oito colunas e/ou um número de simulações maior que 1000, é necessário um aporte computacional superior a um processador Intel Core i3 de 2.27GHz e memória RAM de 4Gb. 

Author
Juliana L. Vendrami juliana.lv@gmail.com 

References
Petchey, O.L. & Gaston, K.J. (2002) Functional diversity (FD), species richness and community composition. Ecology Letters 5: 402-411.
Villéger, S., Mason, N.W.H. & Mouillot, D. (2008) New multimensional functional diversity indices for a multifaceted framework in functional ecology. Ecology 89(8):2290-2301.
http://cran.r-project.org/web/packages/ape/ape.pdf
http://127.0.0.1:17239/library/vegan/doc/diversity-vegan.pdf 
http://www.qhull.org

Note
Os pacotes "vegan","ape", "geometry" e "plotrix" serão instalados durante a execução da função.

See Also
"hclust" do pacote base para a construção do dendograma.
"treeheight" do pacote vegan para o calculo da soma dos braços do dendograma.
"mst" do pacote ape para a construção do gráfico de distância entre as espécies.
"convhulln" do pacote geometry para a construção do polígono no espaço funcional (convex hull).

Examples
## 1º exemplo
x<-matrix(rnorm(90), ncol=3, nrow=30)
rownames(x)<-paste(letters, 1:nrow(x), sep="")
colnames(x)<-paste("trait", 1:ncol(x), sep="")

y<-matrix(sample(30:200,30),ncol=1)
rownames(y)<-paste(letters, 1:nrow(y), sep="")

diversidade.funcional(x,y, salvar.graficos=FALSE)

## 2º exemplo
a<-matrix(sample(0:100000, 180), ncol=3, nrow=60)
rownames(a)<-paste(letters, 1:nrow(a), sep="")
colnames(a)<-paste("trait", 1:ncol(a), sep="")

b<-matrix(sample(150:500,60),ncol=1)
rownames(b)<-paste(letters, 1:nrow(b), sep="")

diversidade.funcional(a,b, salvar.graficos=TRUE)

diversidade.funcional<- function(x, y, nsim=1000, salvar.graficos=FALSE)
{
# classe do objeto
if(class(x)!="matrix") stop("\t", "Erro: os dados devem ser da classe matriz!\n")
if(class(y)!="matrix") stop("\t", "Erro: os dados devem ser da classe matriz!\n")

# tamanho dos objetos
if(dim(x)[1]!=dim(y)[1]) stop("\t", "Erro: número de espécies diferentes nas matrizes!\n")

# premissas
if(nrow(x)<3) stop("\t", "Erro: número mínimo de espécies para o calculo é três!\n")
if(nrow(x)<ncol(x)) stop ("\t", "Erro: número de espécies menor que o número de atributos!\n")

# presença de NA
if(any(is.na(x)==TRUE)) stop("\t", "Erro: presença de NAs na matriz de atributos!\n")
if(any(is.na(y)==TRUE)) stop("\t", "Erro: presença de NAs na matriz de abundância!\n")

# alterando o nome dos objetos
dados<-x
d.abundancia<-y

#salvar os gráficos
if(salvar.graficos==TRUE)
{pdf("Divfuncional%02d.pdf")}

#janela para mostrar os gráficos
par(mfrow=c(2,3), family="serif")

# padronizar os dados de atributos para média=0 e desvio-padrão=1
padronizacao<-scale(dados, center=TRUE, scale=TRUE)


##################### Indice de riqueza funcional - FD (Petchey & Gaston, 2002)

  # calculo da distancia entre as espécies (matriz de distancia)
distancia.sp<-dist(padronizacao, method="euclidean")

  # dendograma
dendograma<-hclust(distancia.sp, method="complete")

  # construção do gráfico
par(las=1, family="serif", cex=0.5)
plclust(dendograma, xlab="", ylab="Distância entre as espécies", sub="", main="Dendograma")

  # verificando se o pacote vegan já está instalado. Em caso negativo, instalar o pacote.
is.installed <- function(m.pct)
{is.element(m.pct, installed.packages()[,1])}

if(is.installed("vegan")==FALSE)
{install.packages("vegan")}
if(is.installed("vegan")==TRUE)
{require("vegan", quietly=TRUE)}

  # calculo do comprimentos dos braços do dendograma  (indice FD)
FD<-treeheight(dendograma)


##################### Indice de uniformidade - FEve Villéger et al 2008

   # instalando o pacote ape
if(is.installed("ape")==FALSE)
{install.packages("ape")}
if(is.installed("ape")==TRUE)
{require("ape", quietly=TRUE)}

   # MST - minimum spanning tree - soma do comprimento mínimo conectando todos os pontos
mst.arvore <- mst(distancia.sp)
mst.arvore.dist<-as.dist(mst.arvore)

   # construção do gráfico
par(las=1, family="serif")
plot(mst.arvore, graph = "nsca",  bty="l", main="Árvore de distância\n entre as espécies", las=1)

   # calculo da abundancia relativa de cada espécie
ab.0<-d.abundancia[d.abundancia>0]
abtotal<-sum(ab.0)
ab.relativa<-as.matrix(ab.0/abtotal)

   # calculo da abundancia relativa entre os pares de espécies
s<-length(d.abundancia)

ab.par.sp<-matrix(0,nrow=s,ncol=s)
for (i in 1:s)
for (j in 1:s)
{ab.par.sp[i,j]<-ab.relativa[i]+ab.relativa[j]}

ab.par.sp.dist<-as.dist(ab.par.sp)  # mostra somente os dados uma única vez

   # calculo da uniformidade ponderada pela abundancia (EW)
EW<-rep(0,(s-1))
int<-1
for (a in 1:length(distancia.sp))
{if(mst.arvore.dist[a]!=0)
{EW[int]<- distancia.sp[a]/ab.par.sp.dist[a]
int<-int+1
}}

  # calculo do indice de uniformidade FEve
PEW<-EW/sum(EW)

PEW.min<-rep(0,s-1)
for (i in 1:(s-1))
{PEW.min[i]<-min(PEW[i] , (1/(s-1)))}

FEve<-round((sum(PEW.min) - (1/(s-1))) / (1-(1/(s-1))),6)

##################### Indice de divergencia (FDiv) - Villéger et al 2008

  # verificando se o pacote geometry está instalado
if(is.installed("geometry")==FALSE)
{install.packages("geometry")}
if(is.installed("geometry")==TRUE)
{require("geometry", quietly=TRUE)}

  # calculo das espécies que ocupam o vertice de um convexo
vert<-convhulln(padronizacao, "Fx") # Fx lista os vertices
vert.num<-as.numeric(vert)
vertice<-sort(unique(vert.num))

sp.vertices<-padronizacao[vertice,]

  # calculo do centro de gravidade
gk<-apply(sp.vertices,2,mean)

  # calculo da distancia Euclidiana ao centro
dist.eucl<-rep(0,s)
for(i in 1:s)
{dist.eucl[i]<-sqrt(sum((padronizacao[i,]-gk)^2))}

  # calculo da distancia media das espécies ao centro de gravidade
media.dist<- (1/s)*sum(dist.eucl)

  # construção do gráfico
par(las=1)
plot(sp.vertices, bty="l", xlab="Atributos1", ylab="Atributos2", main="Divergência com relação\n ao centro de gravidade")
points(x=gk[1], y=gk[2], pch=19)

 ## instalação do pacote plotrix
if(is.installed("plotrix")==FALSE)
{install.packages("plotrix")}
if(is.installed("plotrix")==TRUE)
{require("plotrix", quietly=TRUE)}

draw.circle(x=gk[1], y=gk[2],radius=media.dist)

coordenadas.grafico.x<-padronizacao[,1]
coordenadas.grafico.y<-padronizacao[,2]
for(i in 1:length(dist.eucl))
{
arrows(gk[1],gk[2],coordenadas.grafico.x[i],coordenadas.grafico.y[i],col="blue", angle=0)
}

  # soma dos desvios da abundancia ponderada
ab.pond<-sum( (ab.relativa*(dist.eucl-media.dist)))

  # soma dos desvios absolutos da abundancia ponderada
ab.pond.abs<-sum( (ab.relativa*(abs(dist.eucl-media.dist))))

# calculo do indice de divergencia
FDiv<-(ab.pond+media.dist)/(ab.pond.abs+media.dist)


############################### Modelo nulo ###############################
# FEve
par(las=1)
plot(0:49,0:49, type="n", xlim=c(0,1), ylim=c(0,(0.04*nsim)), xlab=" FEve ", ylab="Frequência", main="Simulação")

require("ape", quietly=TRUE)
div.sample.t <- rep(NA,nsim)
for (q in 1:nsim)
{
amostragem.t<- dados[sample(1:nrow(dados),nrow(dados), replace=T),]
rownames(amostragem.t)<-NULL
padronizacao.t<-scale(amostragem.t, center=TRUE, scale=TRUE)
distancia.sp.t<-dist(padronizacao.t, method="euclidean")
mst.arvore.t <- mst(distancia.sp.t)
mst.arvore.dist.t<-as.dist(mst.arvore.t)
ab.relativa.t<-d.abundancia/sum(d.abundancia[d.abundancia>0])
rownames(ab.relativa.t)<-NULL
s.t<-length(padronizacao.t[,1])
ab.par.sp.t<-matrix(0,nrow=s.t,ncol=s.t)
for (i in 1:s.t)
for (j in 1:s.t)
{ab.par.sp.t[i,j]<-ab.relativa.t[i]+ab.relativa.t[j]}
ab.par.sp.dist.t<-as.dist(ab.par.sp.t)
EW.t<-rep(0,(s.t-1))
int<-1
for (a in 1:length(distancia.sp.t))
{if(mst.arvore.dist.t[a]!=0)
{EW.t[int]<- distancia.sp.t[a]/(ab.par.sp.dist.t[a])
int<-int+1
}}
PEW.min.t<-rep(0,s.t-1)
for (i in 1:(s.t-1))
{PEW.min.t[i]<-min((EW.t[i]/sum(EW.t)) , (1/(s.t-1)))}
div.sample.t[q]<-round((sum(PEW.min.t) - (1/(s.t-1))) / (1-(1/(s.t-1))),3)

stripchart(div.sample.t, method="stack", add=T, pch=15, col="gray")
abline(v=FEve, col="black", lty=2)
}

maior.obs.t<-sum(div.sample.t>=FEve)
p.value.t<-maior.obs.t/length(div.sample.t)
ic.25.t<-round(quantile(div.sample.t, probs=0.025),3)
ic.97.t<-round(quantile(div.sample.t, probs=0.975),3)

# FDiv
par(las=1)
plot(0:49,0:49, type="n", xlim=c(0,1), ylim=c(0,(0.04*nsim)), xlab=" FDiv ", ylab="Frequência", main="Simulação")

require("geometry", quietly=TRUE)
require("ape", quietly=TRUE)
div.sample.2 <- rep(NA,nsim)
for (p in 1:nsim)
{
amostragem.2<- dados[sample(1:nrow(dados),nrow(dados), replace=T),]
rownames(amostragem.2)<-NULL
padronizacao.2<-scale(amostragem.2, center=TRUE, scale=TRUE)
vert.2<-convhulln(padronizacao.2,"Fx")
vert.num.2<-as.numeric(vert.2)
vertice.2<-sort(unique(vert.num.2))
sp.vertices.2<-padronizacao.2[vertice.2,]
gk.2<-apply(sp.vertices.2,2,mean)
s.2<-length(padronizacao.2[,1])
dist.eucl.2<-rep(0,s.2)
for(i in 1:s.2)
{dist.eucl.2[i]<-sqrt(sum((padronizacao.2[i,]-gk.2)^2))}
media.dist.2<- (1/s.2)*sum(dist.eucl.2)
ab.relativa.2<-d.abundancia/sum(d.abundancia[d.abundancia>0])
rownames(ab.relativa.2)<-NULL
ab.pond.2<-sum((ab.relativa.2*(dist.eucl.2-media.dist.2)))
ab.pond.abs.2<-sum((ab.relativa.2*(abs(dist.eucl.2-media.dist.2))))
div.sample.2[p]<-round((ab.pond.2+media.dist.2)/(ab.pond.abs.2+media.dist.2),3)

stripchart(div.sample.2, method="stack", add=T, pch=15, col="gray")
abline(v=FDiv, col="black", lty=2)
}

maior.obs.2<-sum(div.sample.2>=FDiv)
p.value.2<-maior.obs.2/length(div.sample.2)
ic.25.2<-round(quantile(div.sample.2, probs=0.025),3)
ic.97.2<-round(quantile(div.sample.2, probs=0.975),3)

## restaurando o padrao do dispositivo gráfico
if(salvar.graficos==TRUE)
{dev.off()}
par(mfrow=c(1,1))

########## Resultado
Índices<-round(c(FD, FEve, FDiv),3)
p.value<-c("-",p.value.t,p.value.2)
IC25<- c("-",ic.25.t,ic.25.2)
IC97<-c("-",ic.97.t, ic.97.2)
resultado<-data.frame(Índices,p.value,IC25, IC97)
rownames(resultado)<-c("FD", "FEve", "FDiv")
colnames(resultado)<- c("Índices","p-value","IC 2.5%", "97.5%")

return(resultado)
}

diversidade_funcional.r