estatistica aplicada a ecologia, Notas de estudo de Engenharia Florestal

Baixe estatistica aplicada a ecologia e outras Notas de estudo em PDF para Engenharia Florestal, somente na Docsity! Universidade Estadual Paulista Programa de Pós-Graduação Biologia Animal Estatística aplicada à ecologia usando o R Professores responsáveis: Diogo Borges Provete (dbprovete@gmail.com) Fernando Rodrigues da Silva (bigosbio@yahoo.com.br) Thiago Gonçalves Souza (tgoncalves.souza@gmail.com) São José do Rio Preto, SP Abril, 2011     2  SUMÁRIO Objetivo do curso 4 O que você não encontrará nesta apostila 5 Introdução: integrando questões ecológicas e análises estatísticas 6 O melhor caminho para fazer a pergunta certa 8 Introdução ao ambiente de programação R 9 Baixando e instalando a versão base do R 10 Porque usar o R 10 O “workspace” do R e o Tinn-R 11 Os tipos de objeto: criação e manipulação 12 Operações aritméticas básicas 15 Entendendo o arquivo de ajuda 16 Instalando e carregando pacotes 17 Importação e exportação de dados 18 Criação e manipulação de gráficos no R 20 Distribuições estatísticas 18 Funções de probabilidade 23 Funções de distribuição acumulada 24 Distribuição binomial 24 Distribuição Poisson 28 Distribuição Normal 32 Modelos Lineares Generalizados 36 Curva de acumulação de espécies 65     5  Aprofundamento teórico, detalhes matemáticos, e explicação dos algoritmos são informações que infelizmente não serão abordadas neste curso. O foco do curso é a explicação de como cada teste funciona (teoria e procedimentos matemáticos básicos) e sua aplicação em testes ecológicos usando o programa R. Para tanto, o livro dos irmãos Pierre e Louis Legendre (Legendre & Legendre 1998) é uma leitura que permite o aprofundamento de cada uma das análises propostas aqui. Além disso, são de fundamental importância para o amadurecimento em análises ecológicas as seguintes leituras: Manly (1991), Pinheiro & Bates (2000), Scheiner & Gurevitch (2001), Quinn & Keough (2002), Venables & Ripley (2002), Magurran (2004) e Gotelli & Ellison (2004). Figura 1. Estrutura lógica para integrar teorias/questões ecológicas com análises estatísticas (e vice-versa). Lembre-se de que omitimos etapas importantes desta estrutura lógica, como o delineamento experimental, a coleta e organização dos dados, que estão além do objetivo desta apostila. P, R2, F, t, r, Z, AIC, AICc ... Observação Questões Hipóteses biológicas Hipóteses estatísticas Análises estatísticas DECISÃO Unidade amostral Variáveis Covariáveis Escala Predições Hipótese nula Hipótese alternativa TEORIA G en er al iz aç ão O QUE VOCÊ NÃO ENCONTRARÁ NESTA APOSTILA       6  Para a grande maioria dos estudantes [e professores] de biologia a palavra “estatística” traz certa vertigem e aversão. Em geral, alunos e professores consideram este passo um dos mais (se não o mais) problemáticos da pesquisa científica. Para ecologia e, especialmente, ecologia de comunidades, métodos analíticos complexos e que consomem muito tempo para serem realizados tornam a estatística uma tarefa ainda mais distante de ser alcançada (e compreendida). Infelizmente, a maioria opta por não cumprir esta tarefa. Em nossa opinião, muito dessa aversão à estatística se deve às disciplinas introdutórias do curso de graduação em Ciências Biológicas (a maioria, é claro) estarem baseados em um contexto puramente estatístico e com exemplos não-biológicos, sem um programa que integre a ferramenta analítica a um “problema de pesquisa”. De fato, entender exemplos estatísticos com uma lógica puramente estatística não parece uma tarefa trivial para alunos que buscam entender, por exemplo, como processos populacionais, de comunidades e ecossistêmicos determinam a distribuição das espécies. Uma alternativa que pode facilitar a compreensão das análises estatísticas para biólogos (e para todos os cientistas!) é a utilização da lógica do método científico tomando como fator de decisão os resultados estatísticos. Ao final do curso, ou da leitura desta apostila, gostaríamos de que você refletisse um pouco sobre as seguintes questões: (1) qual a principal teoria do meu trabalho? (2) Qual a principal pergunta do meu trabalho? (3) Qual é a unidade amostral, a variável dependente e independente do meu trabalho? A seguir, apresentamos a seqüência lógica que sugerimos que seja aplicada a todo e qualquer teste que utilize estatística frequentista (interpretação objetiva da probabilidade baseada no critério de falseamento de Karl R. Popper). Esta interpretação é, por sua vez, diferente da interpretação subjetiva da probabilidade utilizada no arcabouço da estatística Bayesiana e da Maxima Verossimilhança. É importante ressaltar ainda que a probabilidade (o fator de decisão dos frequentistas, i.e., o tão sonhado “p < 0,05”) representa uma classe de eventos (observados) comparados com uma série de repetições, e portanto o grau de incerteza relacionada a eventos. Todo este arcabouço dos testes de hipóteses estatísticas foi desenvolvido por Jerzy Neyman e Egon S. Pearson (Neyman & Pearson, 1933) adotando a visão Popperiana de que uma observação não fornece confirmação para uma teoria, devido ao problema da indução (para uma discussão mais detalhada veja os cap. 2 e 3 de Godfrey-Smith, 2003). Ao contrário, um teste deveria procurar refutar uma teoria, somente desta forma haveria ganhado conhecimento. Então, segundo o arcabouço de Neyman- Pearson, o teste estatístico procura rejeitar a hipótese nula, e não a confirmação da hipótese alternativa. Numa regressão, por exemplo, se o teste verificar que o coeficiente β é significativo, isto quer dizer que a inclinação da reta é diferente de zero, no entanto a interpretação biológica de uma relação linear entre as duas variáveis deve ser feita à luz das predições da teoria que se pretende testar. Por outro lado, os testes de modelos lineares generalizados em mistos utiliza a INTRODUÇÃO – INTEGRANDO QUESTÕES ECOLÓGICAS E ANÁLISES ESTATÍSTICAS       7  lógica da estatística Bayesiana e da Maxima Verossimilhança. Estes arcabouços utilizam a interpretaçãoo subjetiva da probabilidade. Como uma analogia, o arcabouço frequentista presume que a “verdade” ou todo o universo amostral está numa nuvem, distante e inalcançável, e que somente temos acesso a pequenas amostras de dados, que nesta metáfora, seriam um monte, com o qual chegaríamos o mais próximo possível da nuvem. Seguindo esta metáfora, a estatística Bayesiana e Maxima Verossimilhança assumem que j que a “nuvem” é algo inatingível não devemos considerá-la na análise e que a melhor estimativa que temos são os dados reais que coletamos. Portanto, neste contexto, devemos considerar nossos dados como o universo amostral total. Ao definir a questão de pesquisa é essencial conhecer como a teoria pode ser usada e como e porque ela pode explicar ou ser aplicada à sua questão (Ford 2000). Os modelos gerados pelas teorias podem ser aproveitados para criar suas hipóteses e predições. As hipóteses [científicas] são definidas como explicações potenciais que podem ser retiradas de observações do mundo externo (processo indutivo) ou de componentes de uma teoria (processo dedutivo). Uma hipótese científica, do ponto de vista de Popper, deve ser falseável. As predições são afirmações deduzidas de uma estrutura lógica ou causal de uma teoria, ou induzidas a partir de informações empíricas; em outras palavras, a predição é a conseqüência da hipótese, o resultado esperado se a hipótese for verdadeira. Uma hipótese bem articulada deve ser capaz de gerar predições. Um exercício fundamental para a criação de hipóteses e articulação de suas predições se faz a partir da construção de fluxogramas (Fig. 2). No fluxograma você pode separar cada variável e a relação esperada entre cada uma delas. As setas indicam a relação esperada entre as variáveis (os sinais acima das setas mostram a direção da relação). Setas com espessuras diferentes podem ser usadas como forma de demonstrar a importância relativa esperada para cada variável. Figura 2. Fluxograma representando as predições que foram articuladas a partir da hipótese “as florestas ripárias aumentam a riqueza de macro-invertebrados”.     10  tempo que os executa no seu computador, e não só os leia passivamente. Ainda, por motivo de tempo e espaço não abordaremos todas as questões relacionadas ao uso do R nesta apostila. Logo, aconselhamos que o leitor ao final das aulas você consulte o material sugerido para poder se aprofundar nas questões abordadas. Para começarmos a trabalhar com o R é necessário baixá-lo na página do R project da internet. Então, digite http://www.r-project.org na barra de endereços do seu navegador. Em seguida, clique no link download R embaixo da página, que o levará à pagina do CRAN (Comprehensive R Archive Network). Escolha qualquer página espelho do Brasil para baixar o programa. Escolha o sistema operacional do seu computador e clique em base. Reserve algum tempo posteriormente para explorar esta página do R-project. Existem vários livros (http://www.r-project.org/doc/bib/R-books.html) dedicados a diversos assuntos baseados no R, além disso, estão disponíveis manuais (http://cran.r-project.org/manuals.html) em diversas línguas (http://cran.r-project.org/other-docs.html) para serem baixados gratuitamente. Como o R é um software livre, não existe a possibilidade de o usuário entrar em contato com um serviço de suporte de usuários, muito comuns em softwares pagos. Ao invés disso, existem várias listas de correio eletrônico que fornecem suporte à comunidade de usuários (http://www.r-project.org/mail.html). Nós, particularmente, recomendamos o ingresso nas seguintes listas: R-help, R-sig-ecology, e R_BR (http://www.leg.ufpr.br/doku.php/software:rbr). Este último representa um grupo de usuários brasileiro do programa R. Ainda, existem vários blogs e páginas com arquivos de ajuda e planilhas com comandos, alguns deles podem ser baixados aqui: http://www.nceas.ucsb.edu/scicomp/software/r e http://devcheatsheet.com/tag/r/. Os criadores do R o chamam de uma linguagem e ambiente de programação estatística e gráfica. O R também é chamado de programa “orientado ao objeto” (object oriented programming), o que significa que utilizar o R envolve basicamente a criação e manipulação de objetos em uma tela branca em que o usuário tem de dizer exatamente o que deseja que o BAIXANDO E INSTALANDO A VERSÃO BASE DO R   PORQUE USAR O R?       11  programa execute ao invés de simplesmente pressionar um botão. E vem daí uma das grandes vantagens em se usar o R: o usuário tem total controle sobre o que está acontecendo e também tem de compreender totalmente o que deseja antes de executar uma análise. Na página pessoal do Prof. Nicolas J. Gotelli existem vários conselhos para um estudante iniciante de ecologia. Dentre esses conselhos, o Prof. Gotelli menciona que o domínio de uma linguagem de programação é uma das mais importantes, porque dá liberdade ao ecólogo para executar tarefas que vão além daquelas disponíveis em pacotes comerciais. Além disso, a maioria das novas análises propostas nos mais reconhecidos periódicos em ecologia normalmente são implementadas em linguagem R, e os autores incluem normalmente o código fonte no material suplementar dos artigos, tornando a análise acessível. A partir do momento que essas análises ficam disponíveis (seja por código fornecido pelo autor ou por implementação em pacotes pré-existentes), é mais simples entendermos a lógicas de análises complexas, especialmente as multivariadas, com nossos próprios dados realizando-as passo a passo. Sem a utilização do R, normalmente temos que contatar os autores que nem sempre são acessíveis. Uma última vantagem é que por ser um software livre, a citação do R em artigos é permitida e até aconselhável. Para saber como citar o R, digite citation()na linha de comando. Para citar um pacote específico, digite citation()com o nome do pacote entre aspas dentro dos parênteses. Neste ponto, esperamos ter convencido você leitor de que aprender a utilizar o R tem inúmeras vantágens, vai ser difícil no começo mas continue e perceberá que o investimento vai valer à pena no futuro. Com o R é possível manipular e analisar dados, visualizar gráficos e escrever desde pequenas linhas de comando até programas inteiros. O R é a versão em código aberto de uma linguagem de programação inventada nos anos 1980 no Bell Labs chamada de S. Essa linguagem tornou-se bastante popular e vários produtos comerciais que a usam estão disponíveis, como o S-PLUS, SPSS, STATA e SAS. Um aspecto digno de nota é que a linguagem R, ao contrário de outras linguagem como Fortran e C, é uma linguagem interpretada, o que a faz ser mais fácil de programar, pois processa linhas de comando e as transforma em linguagem de máquina (código binário que o computador efetivamente lê), mas isso diminui a velocidade de processamento. O “WORKSPACE” DO R E O TINN-R       12  Nas linhas de comandos do R haverá um sinal de >, que indica o prompt, representando que o R está pronto para receber comandos. Se uma linha de comando não está completa, aparecerá um sinal de +, indicando que você poderá continuar a digitar aquela linha. Para que o prompt apareça novamente, pressione Esc. Para que os comandos sejam executados, pressione Enter. Para criar objetos, podemos utilizar os símbolos -> ou = . Estes símbolos representam que queremos “guardar” a informação dentro do objeto. Neste curso iremos utilizar o R em conjunto com um editor, o Tinn-R. Existem vários editores para a linguagem R, como o RStudio, Eclipse etc. (veja uma lista não exaustiva em http://en.wikipedia.org/wiki/R_(programming_language)), mas preferimos o Tinn-R por ser de mais fácil utilização e por possibilitar o destaque das sintaxes de programação, diminuindo erros de digitação tão comuns. E ainda, é possível salvar os scripts para continuar a trabalhar neles posteriormente. Para baixá-lo, vá até http://www.sciviews.org/Tinn-R/ e faça o download do programa. Assim que o instalar, somente será necessário clicar no ícone do Tinn-R e o R abrirá automaticamente. Toda vez que terminar de escrever uma linha de comando, pressione Ctrl+Enter para enviá-la para o R. Para saber qual é o diretório de trabalho do R, ou seja, em qual pasta o programa salvará arquivos, digite: >get.wd() É possível mudar o diretório de trabalho do R de acordo com as necessidades do usuário. Então, como exercício para este curso, clique em Arquivo>mudar dir. e defina o diretório para uma pasta deste curso dentro de Meus documentos. Nós recomendamos mudar o diretório sempre que um novo conjunto de análises for feito como, por exemplo, quando for mudar das análises do primeiro capítulo da sua dissertação para o segundo, escolha a pasta onde estarão os dados deste capítulo como diretório de trabalho. Existem cinco classes de objetos na linguagem R: vetor, matriz, data frame, funções e lista. Vetor Existem três tipos de vetores: o vetor de caracteres, numérico e o lógico. OS TIPOS DE OBJETOS: CRIAÇÃO E MANIPULAÇÃO       15  Data frame O mesmo que uma matriz, mas aceita vetores de tipos diferentes. Este é o tipo mais comum de objeto que iremos usar ao longo deste curso. Um data frame permite incluir num mesmo objeto vetores numéricos e de caracteres, por exemplo: >comunidade<- data.frame(especies = c("D.nanus", "S.alter","I.guentheri", "A. callipygius"), habitat = factor(c("Folhiço", "Arbóreo", "Riacho", "Poça")), altura = c(1.1, 0.8, 0.9, 1), distancia = c(1, 1.7, 0.6, 0.2)) >class(comunidade) >xy=as.data.frame(xy)#converte (coerce) a matriz que criamos acima numa data frame >class(xy) #testa a conversão >str(comunidade) >fix(comunidade) >edit(comunidade) Lista Uma lista é um objeto que consiste de um conjunto de objetos ou componentes ordenados de forma hierárquica. Por exemplo, é possível construir uma lista com uma matriz, um vetor lógico, etc. > Lista.ex <- list(name="Toyoyo", wife="Rafaela", no.children=2, child.ages=c(2,6)) Muitos testes produzem objetos em formato de listas como resultado. Às vezes é útil extrair partes de uma lista para que possam ser utilizados posteriormente. >Lista.ex$name O R também pode ser utilizado como uma calculadora. Faça algumas operações aritméticas com os objetos que você acabou de criar, por exemplo: OPERAÇÕES ARITMÉTICAS BÁSICAS       16  >a*2 >b*3 #observe o que aconteceu? Como foi feita essa operação? >b[1]*3 #e agora? >b/4 >2+3 >3^3 >log(2)#observe o que aconteceu? Este é a função que calcula o logaritmo neperiano (ln). >log10(2) #compare o resultado anterior com este. São diferentes? >sqrt(3) >sum(a) >mean(b) >sum(b)/length(a) >pi >cor(a,b) >cor.test(a,b) ?cor.test Um importante passo para ter certa intimidade com a linguagem R é aprender a usar a ajuda de cada função. Além disso, existem uma função (RSiteSearch) e um pacote (sos) que também auxiliam o usuário a realizar uma análise quando não se sabe qual (e se) a mesma já foi implementada no R. Para utilizar o RSiteSearch, digite um tema ou o nome de uma análise entre aspas no argumento da função, como no exemplo abaixo: >RSiteSearch("analysis of variance") A função irá buscar na página do R na internet qual(is) função está(ão) disponível(is) para implementar aquela dada análise. Se o pacote sos estiver instalado e carregado, basta digitar: >???”analysis of variance” e o navegador de internet abrirá uma página mostrando qual(is) funções executam aquela análise. Também é necessário acesso à internet. Outra ferramenta de busca é a página ENTENDENDO O ARQUIVO DE AJUDA       17  http://www.rseek.org na qual é possível buscar por um termo não só nos pacotes do R, mas também em listas de emails, manuais, páginas na internet e livros sobre o programa. Vamos fazer um exercício para nos ambientarmos com a página de ajuda do R, digite: >?aov O arquivo de ajuda do R possui geralmente nove ou dez tópicos: Description - resumo da função Usage*- como utilizar a função e quais os seus argumentos Arguments* - detalha os argumentos e como os mesmos devem ser especifidados Details - detalhes importantes para se usar a função Value - mostra como interpretar a saída (output) da função (os resultados) Notes - notas gerais sobre a função Authors - autores da função References - referências bibliográficas para os métodos usados pra construir a função See also - funções relacionadas Examples* - exemplos do uso da função. Às vezes pode ser útil copiar esse trecho e colar no R para ver como funciona e como usar a função. O R é um ambiente de programação e existem atualmente mais de 3000 pacotes que desempenham funções específicas e que precisam ser instalados e carregados independentemente. Os pacotes stats e base já vêm instalados e carregados, são estes pacotes que possuem as funções para o cálculo de modelos lineares simples, como teste t, ANOVA, χ2, glm etc. A função que instala pacotes no R é a install.packages(). Ao longo deste curso utilizaremos vários pacotes, entre eles o vegan, para instalá-lo, utilize: >install.packages(“vegan”) para instalar vários pacotes ao mesmo tempo, utilize a função c()para criar um vetor: INSTALANDO E CARREGANDO PACOTES       20  O R é uma poderosa ferramenta para criação e manipulação de gráficos. Os pacotes graphics e grid, que já vêm instalados no R, possuem a função genérica plot(), além de outras como hist(). As funções par() e layout() permitem ainda plotar vários gráficos conjuntamente, formando uma única figura. Alguns pacotes foram desenvolvidos especialmente para manipulação de gráficos, como lattice, ggplot2, ggobi e rgl. Estes pacotes nos permitem fazer praticamente todos os tipos de gráficos, incluindo 3-D e mapas em relevo. Para visualizar uma parte das potencialidades dos pacotes, instale e carregue-os. Digite no prompt do R demo(lattice) e vá apertando Enter. Faça o mesmo com o ggplot2. Neste módulo iremos demonstrar algumas das potencialidades gráficas do R. Reiteramos que esses pacotes são um mundo em si só. Logo, convidamos o leitor a ler e explorar a literatura sugerida abaixo, consultar os quadros resumos, além de acessar as seguintes páginas da internet: http://research.stowers-institute.org/efg/R/ http://addictedtor.free.fr/graphiques/ http://www.gnuplot.info/ http://gnuplot.sourceforge.net/demo_4.2/ http://www.statmethods.net/advgraphs/parameters.html. As principais funções que possibilitam modificar gráficos no R são: plot()#Função genérica para plotar gráficos #utilize os argumentos xlab e ylab para adicionar legendas aos eixos, use aspas. # bty=”L” retira as molduras das partes direita e superior. # xlim e ylim determina os limites das escalas dos eixos. # cex modifica o tamanho dos pontos. # pch modifica o tipo do ponto # col modifica as cores dos pontos. Veja também a ajuda da função par(). hist()# plota um histograma barchart()# plota um gráfico de barras CRIAÇÃO E MANIPULAÇÃO DE GRÁFICOS NO R       21  locator()#localiza uma coordenda x-y no gráfico, utilize o argumento 1, 2 etc para definir quantos pontos quer localizar text()#adiciona um texto arrows()#adiciona uma seta mtext()adiciona um texto nas margens do gráfico box()#adiciona uma moldura segments()#adiciona uma linha legend()#adiciona legendas no alto e embaixo points()#adiciona pontos no gráfico lines()#adiciona linhas no gráfico par()#divide o layout e plota vários gráficos, utilize o argumento mfrow=c(2,2) para especificar o número de linhas e colunas. Neste caso a função par(mfrow=c(2,2)) cria uma janela para que quatro gráficos sejam visualizados (i.e., duas linhas e duas colunas) layout()#divide o layout e plota vários gráficos, utilize o argumento layout(matrix(1:4, ncol=2, nrow=2)) pra definir o número de colunas e linhas. O pacote lattice permite fazer gráficos univariados e multivariados de alto nível. Além disso, ele permite criar objetos da classe trellis que podem ser exportados e modificados. xyplot()#função do lattice para gráficos univariados bwplot()# plota um boxplotcoplot()#plota vários gráficos com estilos diferentes Exercícios 1) Carregue o pacote lattice e o conjunto de dados quakes, data(quakes), plote os dados utilizando a função xyplot(). 2) Carregue o conjunto de dados melanoma e utilizando a função plot() faça um gráfico com o tamanho dos pontos 24, legenda do eixo x “Frequência”, legenda do eixo y “Anos” e sem as molduras da direita e superior. 3) Crie dois conjuntos de dados quaisquer e combinando as funções abline() e lm() calcule uma regressão linear simples e ajuste uma reta que indique o modelo.     22  4) Crie um conjunto aleatório de números com distribuição normal e dê nome a este objeto. Utilize a função hist() para plotar um gráfico com as barras em cor cinza. a) Utilize a função points() para criar um ponto em formato de círculo no eixo x no lugar da média. b) Agora crie dois pontos verdes em formato de triângulo verde invertido no lugar dos 2 quantis. c) Crie uma legenda no canto superior esquerdo com os símbolos utilizados (triângulo e círculo), com os significado (média e quantil). d) Pinte de vermelho e verde os símbolos. Quais funções você aprendeu? Uma linguagem de programação é uma linguagem como qualquer outra, e sua aprendizagem exige domínio de vocabulário e sintaxe. O vocabulário da linguagem R são as funções e comandos. Então, sempre que um módulo acabar, lembre-se de tomar nota das funções e comandos, bem como para que serve cada uma delas. Utilize o marcador # em frente a uma função para explicar a sua utilidade. Você se lembra de todas que aprendeu hoje? Uma distribuição estatística é definida como uma função que define uma curva. A área sob essa curva determina a probabilidade de ocorrência de um dado evento. Variáveis aleatórias: A variável aleatória (X) é uma variável que tem um valor único (determinado aleatoriamente) para cada resultado de um experimento. A palavra aleatória indica que em geral só conhecemos aquele valor depois do experimento ser realizado. Exemplos de variáveis aleatórias: a. Número de presas capturadas em um determinado dia; b. Comprimento de um peixe adulto selecionado aleatoriamente. As variáveis aleatórias podem ser discretas ou contínuas. DISTRIBUIÇÕES ESTATÍSTICAS     25  Exemplo Há uma probabilidade de 0,30 de um girino, ao forragear em um corpo d’água, ser predado por uma larva de odonata. Determine as probabilidades de que, dentre seis girinos que estão forrageando no corpo d’água, 0, 1, 2, 3, 5 ou 6 sejam predados. Trace um histograma dessa distribuição de probabilidade. Solução Admitindo que a escolha seja aleatória, fazemos n = 6, q = 0,30 e, respectivamente, X = 0, 1, 2, 3, 4, 5 e 6 na fórmula da distribuição binomial: 𝑝 𝑋 = 𝑛 𝑋 𝑞!(1 − 𝑞)!!!   Figura 5. Histograma da distribuição binomial com n = 6 e q = 0,30. Número de girinos predados  ( ) ( ) ( ) 118,070,030,0 0 6 0 60 ≈      =p ( ) ( ) ( ) 060,070,030,04 6 4 24 ≈      =p ( ) ( ) ( ) 303,070,030,0 1 6 1 51 ≈      =p ( ) ( ) ( ) 010,070,030,05 6 5 15 ≈      =p ( ) ( ) ( ) 324,070,030,0 2 6 2 42 ≈      =p ( ) ( ) ( ) 001,070,030,0 6 6 6 06 ≈      =p ( ) ( ) ( ) 185,070,030,0 2 6 3 33 ≈      =p     26  REALIZANDO O MESMO EXERCÍCIO NO PROGRAMA R: Comandos Existem quatro funções que podem ser utilizadas para gerar os valores associados à distribuição binomial. Você pode obter uma lista completa das mesmas e as suas opções com o comando help: >help(Binomial) Quando o número de tentativas (size) e a probabilidade de sucesso são conhecidos para cada evento (prob) é possível utilizar o comando abaixo para descobrir a probabilidade para qualquer valor da variável x. >dbinom(x, size, prob) No caso do exemplo acima, para descobrirmos qual a probabilidade de dois girinos serem predados, precisamos digitar o seguinte comando: >dbinom (2, size = 6, prob = 0.3) 0.324135 A probabilidade de três girinos serem predados >dbinom (3, size = 6, prob = 0.3) 0.18522 Função de probabilidade acumulativa - Para descobrir a probabilidade de valores menores ou iguais a X utilizamos o comando: >pbinom(q, size, prob) Para descobrirmos qual a probabilidade de dois ou menos girinos (0, 1) serem predados, precisamos digitar o seguinte comando: >pbinom (2, size = 6, prob = 0.3) 0.74431 Para descobrirmos qual a probabilidade de que cinco ou menos girinos (0, 1, 2, 3, 4) sejam predados, precisamos digitar o seguinte comando:     27  >pbinom (5, size = 6, prob = 0.3) 0.999271 Inverso da função de probabilidade acumulativa - Um exemplo contrário ao comando anterior é utilizado quando um valor de probabilidade é fornecido e o programa retorna o valor de X associado a ele. Para isso utiliza-se o seguinte comando: >qbinom(p, size, prob) Qual o valor de X (número de girinos predados) associado à probabilidade de 0,74? >qbinom(0.74, size = 6, prob = 0.3) 2 Qual o valor de X (número de girinos predados) associado a probabilidade de 0,99? >qbinom(0.99, size = 6, prob = 0.3) 5 Finalmente, números aleatórios podem ser gerados de acordo com a distribuição binomial com o seguinte comando: >rbinom(n, size, prob) Por exemplo, para gerar dez números aleatórios de uma distribuição binomial com 20 tentativas e probabilidade 0,63. >rbinom(10, size = 20, prob = 0.63) Você pode plotar o gráfico da função massa de distribuição através do seguinte comando: >plot(dbinom(seq(0,6, by =1), size = 6, prob = 0.3), type ="h", xlab = "Número de girinos predados", ylab = "Probabilidade", main = "Função massa de probabilidade") O gráfico da função de probabilidade acumulada pode ser plotado com o seguinte comando: >plot(pbinom(seq(0,6, by =1), size = 6, prob = 0.3),type ="h", xlab = "Número de girinos predados", ylab = "Probabilidade", main = "Função de probabilidade acumulada")     30  >dpois (5, lambda = 10) 0.03783327 A probabilidade de que oito borboletas visitem uma flor é: >dpois (8, lambda = 10) 0.1125 Função de probabilidade acumulativa - Para descobrir a probabilidade de valores menores ou iguais a X utilizamos o comando: >ppois(x, lambda) Para descobrirmos qual a probabilidade de duas ou menos visitas (1) à flor, precisamos digitar o seguinte comando: >ppois (2, lambda = 10) 0.00276 A probabilidade de cinco ou menos visitas (1, 2, 3, 4) à flor é: >ppois (5, lambda = 10) 0.06708 Inverso da função de probabilidade acumulativa - Um exemplo contrário ao comando anterior é quando você fornece um valor de probabilidade e o programa retorna o valor de X associado a ele. Para isso usa-se o seguinte comando: >qpois (p, lambda) Qual o valor de X (número de visitas) associado à probabilidade de 0.8? >qpois (0.8, lambda = 10) 13 Qual o valor de X (número de visitas) associado a probabilidade de 0.1? >qpois (0.1, lambda = 10) 6     31  Finalmente números aleatórios podem ser gerados de acordo com a distribuição Poisson com o seguinte comando: >rpois (n, lambda) Por exemplo, para gerar dez números aleatórios de uma distribuição Poisson com média (λ ) 22. >rbinom(10, lambda = 22) Você pode plotar o gráfico da função massa de distribuição através do seguinte comando: >plot(dpois(seq(1,10, by =1), lambda = 10), type ="h",xlab = "Número de visitas", ylab = "Probabilidade", main = "Função massa de probabilidade") O gráfico da função de probabilidade acumulada pode ser plotado com o seguinte comando: >plot(ppois(seq(1,10, by =1), lambda = 10),type ="h", xlab = "Número visitas", ylab = "Probabilidade", main = "Função de probabilidade acumulada") Podemos usar a distribuição de Poisson como uma aproximação da distribuição Binomial quando n, o número de tentativas, for grande e p ou 1 – p for pequeno (eventos raros). Um bom princípio básico é usar a distribuição de Poisson quando n ≥ 30 e n.p ou n.(1- p) < 5%. Quando n for grande, pode consumir muito tempo em usar a distribuição binomial e tabelas para probabilidades binomiais, para valores muito pequenos de p podem não estar disponíveis. Se n(1-p) < 5, sucesso e fracasso deverão ser redefinidos de modo que Np < 5 para tornar a aproximação precisa. >plot(dbinom(seq(1,50, by =1), size =50, prob = 0.09), type ="h", ylab = "Probabilidade", main = "Distribuição Binomial") >plot(dpois(seq(1,50, by =1), lambda = 50*0.09), type ="h", ylab = "Probabilidade", main = "Distribuição Poisson")     32  A distribuição normal é uma das mais importantes distribuições com probabilidades contínuas. Conhecida também como Distribuição de Gauss ou Gaussiana. Esta distribuição é inteiramente descrita por parâmetros de média (µ) e desvio padrão (σ), ou seja, conhecendo-se estes parâmetros consegue-se determinar qualquer probabilidade em uma distribuição Normal. A importância da distribuição normal como um modelo de fenômenos quantitativos é devido em parte ao Teorema do Limite Central. O teorema afirma que "toda soma de variáveis aleatórias independentes de média finita e variância limitada é aproximadamente Normal, desde que o número de termos da soma seja suficientemente grande" (Fig. 7). Independentemente do tipo de distribuição da população, na medida em que o tamanho da amostra aumenta, a distribuição das médias amostrais tende a uma distribuição Normal. Figura 7. Gráficos demonstrando que mesmo com um grande número de variáveis aleatórias, as distribuições têm um padrão aproximadamente normal. A distribuição binomial B(n, p) é aproximadamente normal N(np, np(1 − p)) para grande n e para p não tão próximos de 0 ou 1. Enquanto que a distribuição Poisson Pois(λ) é aproximadamente Normal N(λ,  λ) para grandes valores de λ . A função de densidade de probabilidade da distribuição normal com média µ e variância σ2 (de forma equivalente, desvio padrão σ) é assim definida, Variáveis aleatórias com distribuição aproximadamente normal apresentam as seguintes propriedades: – Metade (50%) está acima (e abaixo) da média – Aproximadamente 68% está dentro de 1 desvio padrão da média DISTRIBUIÇÃO NORMAL  ! ( ) 22 2 22 1)( σ µ πσ − = − xexf     35  x = seq(70,130,length = 200) y = pnorm(x, mean=100, sd=10) plot(x, y, type="l", lwd=2, col="red", ylab = "Probabilidade",main ="Função de probabilidade acumulada") Exercícios 1) Uma aranha predadora que vive em flores polinizadas por pequenas mariposas consome em média cinco mariposas por hora. Qual a probabilidade da aranha predar duas mariposas em uma hora selecionada aleatoriamente? 2) Um pesquisador verificou que seis ovos de uma determinada ave são consumidos em média por hora em uma área de nidificação. a) Qual é a probabilidade de que três ovos sejam predados? b) Qual é a probabilidade de que três ou menos ovos sejam predados? 3) Um trabalho recente verificou que 1% dos fígados de cobaias submetidas ao tratamento com álcool apresentavam danos teciduais. Encontre a probabilidade de que mais de um fígado em uma amostra aleatória de 30 fígados apresente danos teciduais usando: a) Distribuição Binomial b) Distribuição Poisson 4) Uma nova técnica de amostragem registra dez indivíduos de lagartos por hora em uma área florestal. Encontre a probabilidade de que quatro ou menos indivíduos sejam registrados em uma hora aleatória. 5) Supondo que a probabilidade de um casal de ursos pandas ter filhotes albinos é de ¼. Se um casal produzir seis filhotes, qual é a probabilidade de que metade deles sejam albinos? 6) Se a probabilidade de um sapo capturar uma mosca em movimento é de 30%. Qual é a probabilidade de que em quatro tentativas ele capture no mínimo três moscas? 7) Um pesquisador extrai 15 amostras de DNA aleatoriamente de um banco de dados que produz 85% de amostras aceitáveis. Qual é a probabilidade de que dez amostras extraídas sejam aceitáveis?     36  8) Um população de crocodilos tem tamanho corporal médio de 400 cm e desvio padrão de 50 cm. Qual a probabilidade de capturarmos um crocodilo dessa população com tamanho entre 390 e 450 cm? 9) O comprimento do antebraço de uma espécie de morcego endêmica do Cerrado é de 4 cm com desvio padrão de 0,25 cm. A partir de qual comprimento os morcegos teriam os antebraços mais compridos nessa população? 10) Suponha que o tempo necessário para um leão consumir sua presa siga uma distribuição normal de média de 8 minutos e desvio padrão de 2 minutos. (a) Qual é a probabilidade de que um leão consuma sua presa em menos de 5 minutos? (b) E mais do que 9,5 minutos? (c) E entre 7 e 10 minutos? 11) A distribuição dos pesos de coelhos criados em uma granja pode muito bem ser representada por uma distribuição Normal, com média 5 kg e desvio padrão 0,9 kg. Um pesquisador comprará 5000 coelhos e pretende classificá-los de acordo com o peso do seguinte modo: 15% dos mais leves como pequenos, os 50% seguintes como médios, os 20% seguintes como grandes e os 15% mais pesados como extras. Quais os limites de peso para cada classificação? Classificação do pesquisador Seja, x1 o valor do peso que separa os 15% mais leves dos demais, x2 o valor do peso que separa os 65% mais leves dos demais, x3 o valor do peso que separa os 85% mais leves dos demais. Muitos métodos estatísticos populares são baseados em modelos matemáticos que assumem que os dados seguem uma distribuição Normal, dentre eles a análise de variância e a Generalized Linear Models (GLM) – Modelos Lineares Generalizados  15%  x1  x2  x3  50%  20%  15%      37  regressão múltipla. No entanto, em muitas situações a suposição de normalidade não é plausível. Conseqüentemente, o uso de métodos que assumem a normalidade pode ser insatisfatório e aumentam a probabilidade de cometermos erros inferenciais (erros do Tipo I e II). Nestes casos, outras alternativas que não pressupoem distribuição normal dos dados são atraentes e mais robustas. Podemos usar modelos lineares generalizados (GLM) quando a variância não é constante, e/ou quando os erros não são normalmente distribuídos. Muitos tipos de dados têm erros não normais. No passado, as únicas maneiras capazes de lidar com esse problema eram a transformação da variável resposta ou a adoção de métodos não paramétricos. Em GLM, assumimos que cada resultado da variável dependente Y seja gerado a partir de uma variedade de diferentes tipos de distribuições que lidam com esse problema: Poisson – úteis para dados de contagem Binomial – úteis para dados com proporções Gamma – úteis para dados mostrando um coeficiente constante de variância Exponencial – úteis com dados de análises de sobrevivência Existem muitas razões para usar GLMs, em vez de regressão linear. Dados de presença-ausência são (geralmente) codificados como 1 e 0, os dados proporcionais são sempre entre 0 e 100%, e os dados de contagem são sempre não-negativos. GLMs usados para 0-1 e dados proporcionais são normalmente baseados em distribuição binomial e para dados de contagem as distribuições de Poisson e binomial negativa são opções comuns. A média, µ, da distribuição depende das variáveis independentes, X, e é calculada através de: 𝐸 𝑌 =  𝜇 =  g!𝟏(𝑋𝛽) onde E (Y) é o valor esperado de Y; Xβ é o preditor linear, uma combinação linear de parâmetros desconhecidos, β; g é a função de ligação. GLM consiste em três etapas: 1. Uma hipótese sobre a distribuição da variável resposta Yi. Isso também define a média e a variância de Yi. (e.x., Distribuição Poisson, Binomial, Gamma). 2. Especificação da parte sistemática. Esta é uma função das variáveis explicativas. 𝑛! =  𝛼 +  𝛽!  𝑋!! +  𝛽!  𝑋!! +⋯+  𝛽!  𝑋!!       40  onde k é o número de parâmetros no modelo estatístico, e L é o valor maximizado da função likelihood para o modelo estimado. Dado um conjunto de modelos candidatos, o modelo preferido é aquele com o valor mínimo de AIC. O valor de AIC não só recompensa goodness- of-fit, mas inclui também uma penalização que é uma função crescente do número de parâmetros estimados. Esta penalidade desencoraja overfitting (aumentando o número de parâmetros livres no modelo melhora a qualidade do ajuste, independentemente do número de parâmetros livres no processo de geração de dados). AICC é AIC com uma correção para amostras finitas: 𝐴𝐼𝐶! = 𝐴𝐼𝐶 +   2𝐾 (𝐾 + 1) 𝑛  − 𝐾 − 1  onde k denota o número de parâmetros do modelo. Assim, AICC é AIC com uma maior penalização para os parâmetros extra. Burnham & Anderson (2002) recomendam o uso do AICC, ao invés de AIC, se n for pequeno ou k é grande. Uma vez que o valor de AICc converge para AIC quando n se torna grande, AICc geralmente devem ser empregados independentemente do tamanho da amostra. Usar AIC, em vez de AICC, quando n não é muitas vezes maior do k2 aumenta a probabilidade de seleção dos modelos que têm muitos parâmetros (overfitting). Uma outra comparação entre os modelos pode ser baseada no cálculo do Peso do Akaike (Akaike weigths - Buckland et al. 1997). Se existem M modelos candidatos, então o peso para o modelo i é: 𝑊𝑖 =   𝑒𝑥𝑝 (∆/2) exp ∆12 + exp ∆2 2 +⋯ exp ( ∆𝑚 2 )   onde Δ é a diferença entre o valor do AIC entre modelo i e os modelos restantes. Os pesos do Akaike calculados desta forma são usados para medir a força da evidência em favor de cada um dos modelos, com um grande peso indicando alta evidência. Dez orientações para Seleção de Modelo 1) Cada modelo deve representar uma hipótese (interessante) específica a ser testada. 2) Mantenha os sub-grupos de modelos candidatos curtos. É desaconselhável considerar tantos modelos quanto o número de dados que você tem.     41  3) Verificar a adequação do modelo: use o seu modelo global (modelo mais complexo) ou modelos subglobais para determinar se as hipóteses são válidas. Se nenhum dos modelos se ajustar aos dados, critérios de informação indicarão apenas o mais parcimonioso dos modelos mais pobres. 4) Evitar a dragagem de dados (e.g., procura de padrões após uma rodada inicial de análise). 5) Evite modelos overfitted. 6) Tenha cuidado com os valores faltantes (NA). Lembre-se de que valores faltantes somente para algumas variáveis alteram o tamanho do conjunto de dados e amostras dependendo de qual variável é incluída em um dado modelo. É sugirido remover casos omissos antes de iniciar a seleção de modelos. 7) Use a mesma variável resposta para todos os modelos candidatos. É inadequado executar alguns modelos com variável resposta transformados e outros com a variável não transformada. A solução é usar uma função de ligação diferente para alguns modelos (e.g., identity vs. log link). 8) Quando se trata de modelos com overdispersion, utilize o mesmo valor de c-hat para todos os modelos em um conjunto de modelos candidatos. Para modelos binomiais com trials > 1 ou com Poisson GLM, deve-se estimar o c-hat do modelo mais complexo (modelo global). Se c hat > 1, deve-se usar o mesmo valor para cada modelo do conjunto de modelos candidatos e inclui- lo na contagem dos parâmetros (K). Da mesma forma, para binomial negativa, você deve estimar o parâmetro de dispersão do modelo global e usar o mesmo valor em todos os modelos. 9) Burnham e Anderson (2002) recomendam evitar misturar a abordagem da teoria da informação e noções de significância (ou seja, os valores P). É melhor fornecer estimativas e uma medida de sua precisão (erro padrão, intervalos de confiança). 10) Determinar o ranking das modelos é apenas o primeiro passo. A soma do Peso Akaike é 1 para o modelo de todo o conjunto e pode ser interpretado como o peso das evidências em favor de um determinado modelo. Modelos com grandes valores do Peso Akaike têm forte apoio. Taxas de evidências, valores de importância, e intervalo de confianca para o melhor modelo são outras medidas que auxiliam na interpretação. Nos casos em que o melhor modelo do ranking tem um Peso Akaike > 0,9, pode-se inferir que este modelo é o mais parcimonioso. Quando muitos modelos são classificados por valores altos (ou seja, o delta (Q) AIC (c) < 2 ou 4), deve- se considerar a média dos parâmetors dos modelos de interesse que aparecem no topo. A média dos modelos consiste em fazer inferências com base no conjunto de modelos candidatos, em vez     42  de basear as conclusões em um único "melhor" modelo. É uma maneira elegante de fazer inferências com base nas informações contidas no conjunto inteiro de modelos. Exemplos A partir dos exemplos a seguir irei explicar os comandos básicos necessários para realizar as análises de GLM. É altamente recomendável que vocês recorram aos livros sugeridos no início desta apostila para um aprofundamento no assunto e para que possam realizar análises mais complexas. Carregando pacotes necessários para as análises >library(languageR) >library(nlme) >library(glmmML) >library(lme4) >library(AICcmodavg) >library(bestglm) >library(mgcv) >library(MuMIn) >library(pscl) >library(MASS) >library(bbmle) >library(lattice) >library(AED) ## Esse pacote tem deve ser baixado da página #http://www.highstat.com/book2.htm Primeiro Exemplo >data(RoadKills) ## Carregando dados - Os dados consistem do número de mortes de anfíbios em uma rodovia em 52 sítios em Portugal Teoria: Ecologia de Paisagem Variável dependente: Número de anfíbios mortos Questão: Quais variáveis da paisagem melhor explicam a mortalidade de anfíbios? >RK <- RoadKills ## Renomeando para facilitar     45  TESTE DE HIPÓTESES - Likelihood ratio test (LRT) DEVIANCE = RESIDUAL DEVIANCE = É 2 x a diferença entre o log likelihood do modelo que apresenta um ajuste perfeito (modelo saturado) e o modelo em questão. Quanto menor o residual deviance, melhor o modelo. >drop1(M1,test = "Chi") # A diferença entre as deviance dos modelos apresenta uma distribuição chi- square com p1 - p2 graus de liberdade >DM1 <- glm(TOT.N ~ OPEN.L + MONT.S + SQ.POLIC + SQ.SHRUB + SQ.WATRES + L.WAT.C + SQ.LPROAD + D.PARK, family = poisson, data = RK) >drop1(DM1, test = "Chi") Este resultado indica que podemos retirar a variável SQ.DWATCOUR, pois o modelo sem esta variável tem o mesmo poder de explicação do modelo com esta variável. Repita o processo até que nenhuma variável possa ser retirada do modelo. OVERDISPERSION Contudo a vida não é tão simples, antes de analisar os resultados e realizar as análises de seleção você precisa checar se os seus dados possuem overdispersion. A overdispersion significa que a variância é maior do que a média.     46  Como saber se os dados apresentam overdispersion? >M1 <- glm (TOT.N ~ OPEN.L + MONT.S + SQ.POLIC + SQ.SHRUB + SQ.WATRES + L.WAT.C + SQ.LPROAD + SQ.DWATCOUR + D.PARK, family = poisson, data=RK) >summary(M1) Veja que o resultado mostra que o parâmetro de dispersão para família Poisson tem que ser 1. Nesse caso o parâmetro de dispersão do seu modelo é 270,23/42 = 6,43. Desse modo, seu modelo apresenta overdispersion e você não pode continuar a análise considerando a família Poisson. Existem duas alternativas: corrigir o Poisson com Quasi-Poisson ou usar a distribuição Binomial Negativa. QUASI-POISSON >M4 <- glm(TOT.N ~ OPEN.L + MONT.S + SQ.POLIC+ SQ.SHRUB + SQ.WATRES + L.WAT.C + SQ.LPROAD+ SQ.DWATCOUR + D.PARK, family = quasipoisson, data = RK) >summary(M4)     47  Veja que o parâmetro de dispersão f é estimado em 5,93. Isto significa que todos os erros padrões foram multiplicados por 2,43 (a raiz quadrada de 5,93), e como resultado, a maioria dos parâmetros não são mais significativos. Não escreva na sua dissertação ou artigo que usou uma distribuição Quasi-Poisson. Quasi-Poisson não é uma distribuição. Basta dizer que você fez GLM com distribuição Poisson, detectou overdispersion, e corrigiu os erros padrões usando um modelo Quasi-GLM, onde a variância é dada por f × µ, onde µ é a média e f é o parâmetro de dispersão. Seleção modelos em Quasi-Poisson Quando inserirmos uma variável para a dispersão, os modelos não podem ser comparados por qui-quadrado. Eles são comparados por distribuição F. >drop1(M4, test = "F") Repita o procedimento até que nenhuma variável possa ser retirada do modelo. Modelo final selecionado >M12 <- glm (TOT.N ~ D.PARK, family = quasipoisson, data = RK) Gráfico com os dados ajustado para a curva Quasi-Poisson-Glm e intervalo de confiança de 95% (IC 95%). >G <- predict (M12, newdata = RK, type = "link", se = TRUE) >F <- exp(G$fit)     50  >NB <- glm.nb(TOT.N ~ OPEN.L + MONT.S + SQ.POLIC + SQ.SHRUB + SQ.WATRES + L.WAT.C + SQ.LPROAD + SQ.DWATCOUR + D.PARK, link="log", data=RK) >odTest(NB) O resultado mostra que a LRT entre Poisson e Binomial Negativa com uma diferença na deviance de 141.515 e com grau de liberdade 1 é p < 0.0000. Portanto, Binomial Negativa é melhor que Poisson. Modelos de Binomial Negativa: >NB1 <- glm.nb (TOT.N ~ OPEN.L + MONT.S + SQ.POLIC + SQ.SHRUB + SQ.WATRES + L.WAT.C + SQ.LPROAD + SQ.DWATCOUR + D.PARK, link="log", data=RK) >NB2 <- glm.nb (TOT.N ~ OPEN.L + MONT.S + SQ.POLIC + SQ.SHRUB + SQ.WATRES + L.WAT.C + SQ.LPROAD + D.PARK, link = "log", data = RK) >NB3 <- glm.nb (TOT.N ~ OPEN.L + MONT.S + SQ.SHRUB + SQ.WATRES + L.WAT.C + SQ.LPROAD + D.PARK, link = "log", data = RK) >NB4 <- glm.nb (TOT.N ~ OPEN.L + MONT.S + SQ.SHRUB + L.WAT.C + SQ.LPROAD + D.PARK, link = "log", data = RK) >NB5 <- glm.nb (TOT.N ~ OPEN.L + MONT.S + L.WAT.C + SQ.LPROAD + D.PARK, link = "log", data = RK) >NB6 <- glm.nb (TOT.N ~ OPEN.L + L.WAT.C + SQ.LPROAD + D.PARK, link = "log", data = RK) >NB7 <- glm.nb (TOT.N ~ OPEN.L + L.WAT.C + D.PARK, link = "log", data = RK) >NB8 <- glm.nb (TOT.N ~ OPEN.L + D.PARK, link = "log", data = RK)     51  Seleção automática por AIC: >AIC <- stepAIC(NB1) >AIC Seleção dos modelos por AICc: >AICc <- ICtab (NB1, NB2, NB3, NB4, NB5, NB6, NB7, NB8, type = c("AICc"), weights = TRUE, delta = TRUE, sort = TRUE, nobs = 52) >AICc Likelihood Ratio Test (LRT) >drop1(NB1,test="Chi") Repita o procedimento até que nenhuma variável retirada apresente efeito siginificativo na comparação. Para o modelo final, os autores justificaram a retirada de L.WAT.C porque seu valor estava muito próximo de 0.05. Modelo Final: >NB8 <- glm.nb(TOT.N ~ OPEN.L + D.PARK, link="log", data = RK) >summary(NB8) BINOMIAL NEGATIVA >plot (NB8) QUASI-POISSON >mu <- predict (M12, type = "response") >E <- RK$TOT.N - mu >EP2 <- E / sqrt (7.630148 * mu) >plot(x = mu, y = EP2, main = "Quasi-Poisson",     52  ylab = "residuos", xlab = "predito") abline(h = 0, v = 0) Comparando os resíduos do modelo final da Binomial Negativa e Quasi-Poisson vemos que os resíduos da Binomial não apresentam um padrão, enquanto a Quasi-Poisson apresenta. Então, Binomial é melhor. GLM BINOMIAL Agora mostraremos um exemplo bem simples com dados de presença e ausência. GLM com dados binários ou proporção são também chamados de regressão logística. >data(Boar) >head(Boar) Variável dependente: presença ou ausência de tuberculose. Variável independente: Comprimento do javali (cabeça-tronco). >B1 = glm ( Tb ~ LengthCT, family = binomial, data = Boar) >summary(B1) Likelihood Ratio Test: >drop1 (B1, test="Chi") Função para fazer o gráfico: >MyData <- data.frame (LengthCT = seq (from = 46.5, to = 165, by = 1)) >Pred <- predict (B1, newdata = MyData, type = "response")     55  Visualização dos resíduos: >EP = resid(Deer8,type = "pearson") >mu = predict(Deer8,type = "response") >E = Tbdeer$DeerPosProp - mu >plot(x = mu,y = EP, main="Pearson residuals") >plot(Deer8) Generalized Mixed Effects Models São usados para modelos mais complexos com design em blocos, medidas repetidas, split plot e dados aninhados. Aprensenta dois efeitos dentro da formúla do modelo: EFEITO FIXO - depende somente da média – as variáveis independentes de interesse. EFEITO ALEATÓRIO - depende somente da variância (não queremos medir o efeito, e.g. blocos). Exemplo 1 >data(RIKZ) Riqueza de animais marinhos bentônicos em nove praias, cada praia com cinco amostras. NAP = altura da estação de amostral em relação ao nível da maré PERGUNTA: Existe relação positiva entre a riqueza e a NAP? Transforma praia em fator: >RIKZ$fBeach <- factor(RIKZ$Beach) Modelo >Mlme1 <- lme (Richness ~ NAP, random = ~1 | fBeach, data=RIKZ) summary (Mlme1)     56  Utilizando praia como efeito aleatório permite que cada praia tenha um intercepto diferente. Se o StdDev do efeito aleatório for zero, todos os interceptos ficam na linha predita. Veja o gráfico abaixo. Função para fazer o gráfico: >F0 <- fitted(Mlme1,level=0) >F1 <- fitted(Mlme1,level=1) >I <- order(RIKZ$NAP) >NAPs <- sort(RIKZ$NAP) >plot(NAPs,F0[I],lwd=4,type="l",ylim=c(0,22), ylab="Riqueza de espécies",xlab="NAP") for (i in 1:9){ x1<-RIKZ$NAP[RIKZ$Beach==i] y1<-F1[RIKZ$Beach==i] K<-order(x1) lines(sort(x1),y1[K]) } >text(RIKZ$NAP,RIKZ$Richness,RIKZ$Beach,cex=0.9) Suponha que a relação entre riqueza de espécies e NAP é diferente em cada praia. Isto implica em que temos de incluir um interação entre NAP*Praia no modelo. Mas isso tem um custo muito alto elevando o modelo para 17 parâmetros. E não estamos interessados no efeito da praia. Contudo, não podemos ignorar uma possível variação entre praias e na interação NAP*Praias. Se fizermos isso, a variação sistemática vai aparecer nos resíduos, levando à inferências erradas. Podemos aplicar o Mixed Effects Model com intercepto e slope (inclinação) aleatórios. >Mlme2 <- lme (Richness ~ NAP, random = ~ 1 + NAP | fBeach, data = RIKZ) >summary(Mlme2)     57  O valor 3,54 é a quantidade de variação no intercepto da população. O valor 1,71 é a variação no slope (inclinação) nas nove praias. A correlação mostra que praias com interceptos mais altos também tem inclinação negativa mais alta. Veja o gráfico abaixo. Função para fazer o gráfico: >F0 <- fitted(Mlme2,level=0) >F1 <- fitted(Mlme2,level=1) >I <- order(RIKZ$NAP) >NAPs <- sort(RIKZ$NAP) >plot(NAPs,F0[I],lwd=4,type="l",ylim=c(0,22), ylab="Riqueza de espécies",xlab="NAP") for (i in 1:9){ x1<-RIKZ$NAP[RIKZ$Beach==i] y1<-F1[RIKZ$Beach==i] K<-order(x1) lines(sort(x1),y1[K]) } >text(RIKZ$NAP,RIKZ$Richness,RIKZ$Beach,cex=0.9) Likelihood em Mixed Models MAXIMUM LIKELIHOOD (ML) - escolhe os parâmetros tal que o valor de L é máximo. O problema é que ML ignora o fato que intercepto e slope são estimados no modelo. RESTRICTED MAXIMUM LIKELIHOOD (REML) - corrige o grau de liberdade incluindo o intercepto e o slope. Transformar algumas variáveis em fatores: >RIKZ$fExp <- RIKZ$Exposure >RIKZ$fExp[RIKZ$fExp==8]<- 10 >RIKZ$fExp <- factor(RIKZ$fExp,levels = c (10,11)) Modelos com ML e com REML:     60  PASSO 4 - Modelo Final com REML >B2 <- lme (Richness ~ NAP + fExp, data = RIKZ, random = ~1 | fBeach, method = "REML") >plot(B2) Exemplo Abelhas Os dados são aninhados com múltiplas observações na mesma colméia. No total são 24 colméias com três medidas por colméia. Mostrar comando VarIdent >data(Bees) Como variável dependente temos densidade de esporos medido em cada colméia. A variável independente Infection quantifica o grau de infecção, com valores 0, 1, 2 e 3. Embora mixed effects modelling podem lidar com um certo grau de dados desbalanceados, neste caso, é melhor converter a variável Infection em 0 (sem infecção) e 1 (infectado) porque existem poucas observações com valores 2 e 3. Transformar a variável Infection em presença e ausência: >Bees$Infection01 <- Bees$Infection >Bees$Infection01[Bees$Infection01 > 0] <- 1 >Bees$fInfection01 <- factor(Bees$Infection01) Transformar colméia em fator e logaritimizar esporos: >Bees$fHive <- factor(Bees$Hive) >Bees$LSpobee <- log10(Bees$Spobee + 1) Plotar os dados por colméia: >op <- par(mfrow = c(1, 2), mar = c(3, 4, 1, 1)) >dotchart(Bees$Spobee, groups = Bees$fHive) >dotchart(Bees$LSpobee, groups = Bees$fHive)     61  >par(op) Começaremos com uma regressão linear e plotaremos os resíduos por colmeia: >M1 <- lm (LSpobee ~ fInfection01 * BeesN, data = Bees) >E1 <- rstandard(M1) >plot (E1 ~ Bees$fHive, ylab = "Resíduos", xlab = "Colméias") >abline (0, 0) Veja que algumas colméias apresentam os três resíduos acima do esperado, enquanto outras possuem três resíduos abaixo do esperado. Temos a opção de colocar colméia como random effect. Vantagens (1) requer um parâmetro extra (variância do intercepto), comparado com regressão linear que requer 23 parâmetros extras. (2) Podemos fazer afirmações para colméias em geral não só para as 24 colméias do estudo. Selecionando random effect >M1 <- lme(LSpobee ~ fInfection01 * BeesN, random = ~ 1 | fHive, method = "REML", data = Bees) >M2 <- lme(LSpobee ~ fInfection01 * BeesN, random = ~ 1 + BeesN | fHive, method = "REML", data = Bees) >M3 <- lme (LSpobee ~ fInfection01 * BeesN, random = ~ 1 + fInfection01 | fHive, method = "REML", ` data = Bees) >anova(M1,M2) >anova(M1,M3) Verificando o modelo selecionado: >plot (M1, col = 1) plota por infecção: >boxplot (LSpobee ~ fInfection01, data = Bees, varwidth = TRUE)     62  Veja que há diferença na variação entre as categorias. Inserimos um comando para dizer que as variâncias para infecção são diferentes. varIdent = permite modelar diferentes variâncias para variáveis categóricas. >M1 <- lme (LSpobee ~ fInfection01 * BeesN, random = ~ 1 | fHive, method = "REML", data = Bees) >M4 <- lme (LSpobee ~ fInfection01 * BeesN, random = ~ 1 | fHive, method = "REML", data = Bees, weights = varIdent (form = ~ 1 | fInfection01)) >anova (M1,M4) Selecionando estrutura fixa: >M7full<- lme (LSpobee ~ fInfection01 * BeesN, random = ~ 1|fHive, weights = varIdent(form = ~ 1 | fInfection01), method = "ML", data = Bees) >M7sub <- update(M7full, .~. -fInfection01 : BeesN ) >anova (M7full,M7sub) >M8full <- lme (LSpobee ~ fInfection01 + BeesN, random = ~ 1|fHive, method = "ML", data = Bees, weights = varIdent(form =~ 1 | fInfection01)) >M8sub1 <- update (M8full, .~. -fInfection01 ) >M8sub2 <- update (M8full, .~. -BeesN ) >anova(M8full,M8sub1) >anova(M8full,M8sub2) >M9full<-lme(LSpobee ~ fInfection01, random = ~ 1|fHive, method="ML", data = Bees, weights = varIdent(form =~ 1 | fInfection01)) >M9sub1<-update(M9full, .~. -fInfection01 ) >anova(M9full,M9sub1) Modelo final: >Mfinal <- lme (LSpobee ~ fInfection01, random =~ 1|fHive, data = Bees, weights = varIdent (form = ~ 1 | fInfection01), method = "REML")     65  Responda: A probabilidade da presença de predadores está relacionada com o tamanho das poças d’água? EXERCÍCIO 3 - Carregue os dados da planilha Solea.csv Variável dependente = presença ou ausência do peixe Solea solea num estuário em Portugal. Variáveis independentes = 11 variáveis preditoras Unidade amostral = cada área de coleta ou ponto de coleta no estuário Teoria: Ecologia de Paisagem Responda: Quais variáveis melhor explicam a presença de Solea solea nos berçários de Portugal? Curvas de acumulação de espécies, algumas vezes chamadas de curva do coletor, são representações gráficas que demonstram o número acumulado de espécies registradas (S) em função do esforço amostral (n). O esforço amostral pode ser o número de indivíduos coletados, ou uma medida tal como o número de amostras (e.g., quadrados) ou tempo amostral (e.g., meses). Colwell & Coddington (1994) sugeriram um método que consiste em montar várias curvas adicionando-se as amostras em uma ordem aleatória. Após construir várias curvas com este método, pode-se calcular uma curva do coletor média (baseada na riqueza média para cada número de amostra) e expressar a variação possível em torno dessa média. É importante frisar que esta variação não corresponde ao conceito estatístico de intervalo de confiança, já que é calculada por repetições das mesmas unidades amostrais (Santos 2003). Se as curvas de acumulação de espécies atingem um ponto em que o aumento do esforço de coleta não implica num aumento no número de espécies, isto significa que aproximadamente toda a riqueza da área foi amostrada (Fig. 8). CURVA DE ACUMULAÇÃO DE ESPÉCIES      66  Figura 8. Exemplo de uma curva de acumulação de espécies.   RAREFAÇÃO Esse método nos permite comparar o número de espécies entre comunidades quando o tamanho da amostra ou o número de indivíduos (abundância) não são iguais. A rarefação calcula o número esperado de espécies em cada comunidade tendo como base comparativa um valor em que todas as amostras atinjam um tamanho padrão, ou comparações baseadas na menor amostra ou com menos indivíduos (dentre todas amostras possíveis). Se considerarmos n indivíduos (n < N) para cada comunidade, quantas espécies iríamos registrar?   𝐸(𝑆) = 1 −   (𝑁 −  𝑁!)/𝑛 𝑁/𝑛 Onde: E(S) = Número de espécies esperado N = Número total de indivíduos na amostra Ni = Número de indivíduos da iésima espécie n = tamanho da amostra padronizada (menor amostra) Gotelli & Collwel (2001) descrevem este método e discutem em detalhes as restrições sobre seu uso na ecologia: i) as amostras a serem comparados devem ser consistentes do ponto de vista taxonômico, ou seja, todos os indivíduos devem pertencer ao mesmo grupo taxonômico; ii) as comparações devem ser realizadas somente entre amostras com as mesmas técnicas de coleta; iii) os tipos de hábitat onde as amostras são obtidas devem ser semelhantes; e iv) é um método para estimar a riqueza de espécies em uma amostra menor – não pode ser usado para extrapolar e estimar riqueza.     67  Exemplo: Uma amostra de roedores tem quatro espécies e 42 indivíduos. A abundância de cada espécie foi 21, 16, 3, e 2 indivíduos. Desejamos calcular a riqueza de espécies esperada para amostras com 30 indivíduos. 𝐸 𝑠 = 1 − 42 − 21 /30 42/30   + 1 − 42 − 16 /30 42/30   + 1 − 42 − 3 /30 42/30   + 1 − 42 − 2 /30 42/30   E(30) = 1 + 1 +0.981 + 0.923 E(30) = 3.9 espécies REALIZANDO O MESMO EXERCÍCIO NO PROGRAMA R: Comandos Primeiramente carregue o pacote vegan: >library(vegan) O comando geral para realizar a análise de rarefação é: >rarefy(x, sample, se = FALSE, MARG = 1) Onde: x = comunidade para a qual se deseja estimar a riqueza de espécies sample = tamanho da sub-amostra (n) se = desvio padrão MARG = maneiras de visualizar o resultado – Utilizar número 2 Imagine que você tenha uma planilha aberta no R com o nome rare. Nesta planilha, existem três colunas referentes à três comunidades de roedores, e em cada linha a abundância de cada espécie (exemplo abaixo):     70  Tabela 5. Número de indivíduos registrados de cada espécie de anuros em 14 amostras no noroeste de São Paulo, Brasil. Será utilizado nos exemplos abaixo. Espécies AMOSTRAS Total 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Delian 0 0 6 15 2 2 0 0 0 1 0 5 5 2 38 Dmelan 0 0 0 0 1 0 0 1 0 0 0 0 0 0 2 Dminu 0 2 1 15 8 2 0 1 2 2 0 4 0 2 39 Dnanu 4 0 3 15 2 2 0 7 0 2 0 3 2 2 42 Dmulle 0 0 0 3 12 0 2 0 0 0 0 0 0 0 17 Ebic 0 0 0 0 1 0 0 1 0 0 0 0 1 0 3 Esp 0 0 2 0 0 0 0 0 0 1 0 0 0 0 3 Enat 0 4 1 0 17 0 2 0 1 0 4 0 0 1 30 Halb 5 0 0 0 0 1 0 9 0 1 0 0 4 0 20 Hfab 0 0 0 0 0 0 0 4 0 0 0 0 0 0 4 Hran 14 0 0 5 0 1 0 0 0 2 0 0 8 0 30 Lchaq 0 0 0 0 11 0 3 0 0 0 0 0 0 0 14 Lfus 8 3 2 5 4 2 1 6 1 3 1 2 3 6 47 Llab 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 Riqueza Total 4 3 6 6 9 6 4 7 3 8 2 4 6 5 CHAO 1 Estimador simples do número absoluto de espécies em uma comunidade. É baseado no número de espécies raras dentro de uma amostra. Esse método requer a abundância das espécies. 𝐶ℎ𝑎𝑜! = 𝑆!"# +   𝐹!! 2𝐹!  onde: Sobs = o número de espécies na comunidade F1 = número de espécies observadas com abundância de um indivíduo (espécies singleton) F2 = número de espécies observadas com abundância de dois indivíduos (espécies doubletons). O valor de Chao 1 é máximo quando todas as espécies menos uma são únicas (singleton). Neste caso, a riqueza estimada é aproximadamente o dobro da riqueza observada. Exemplo: Usando os dados da tabela 1 calcule o valor de Chao 1 para a comunidade: Chao 1 = 14 + [(12)/(2*1)] = 14 + (1/2) = 14 + 0,5 Chao 1 = 14,5     71  REALIZANDO O MESMO EXERCÍCIO NO PROGRAMA R: Comandos Carregue os pacotes Vegan e BiodiversityR >library(vegan) >library(BiodiversityR) Imagine que você tenha a mesma tabela acima salva no R com o nome “est”. Após carregar essa tabela você pode obter o valor de Chao 1 através do seguinte comando: >est <- read.table (estimadores, h = T) >Chao1 <-estaccumR (est, permutations = 100) >summary(Chao1, display = “chao”) Outra maneira de conseguir o mesmo valor: >est1 <- colSums(est)## soma abundância de cada linha = abundância total por espécie >Chao1 <- estimateR (est1) >Chao1 CHAO 2 De acordo com Anne Chao, o estimador Chao 1 pode ser modificado para uso com dados de presença/ausência levando em conta a distribuição das espécies entre amostras. Neste caso é necessário somente conhecer o número de espécies encontradas em somente uma amostra e o número de espécies encontradas exatamente em duas amostras. Essa variação ficou denominada Chao 2: 𝐶ℎ𝑎𝑜! = 𝑆!"# +   𝐿! 2𝑀  onde: L = número de espécies que ocorrem apenas em uma amostra (espécies uniques) M = número de espécies que ocorrem em exatamente duas amostras (espécies duplicates) O valor de Chao 2 é máximo quando todas as espécies menos uma são únicas (singletons). Neste caso, a riqueza estimada é aproximadamente o dobro da riqueza observada.     72  Collwel & Coddington (1994) encontraram que o valor de Chao 2 mostrou ser o estimador menos enviesado para amostras com tamanho pequeno. Exemplo: Usando os dados da tabela 1 calcule o valor de Chao 2 para a comunidade: Chao 2 = 14 + [(22)/(2*3)] = 14 + (4/6) = 14 + 0.66 Chao 2 = 14.66 REALIZANDO O MESMO EXERCÍCIO NO PROGRAMA R: Comandos A função poolaccum do pacote vegan apresenta resultados mais completos com valores de riqueza de espécie estimado para cada amostra >est <- read.table (estimadores, h = T) >Chao2 <- poolaccum (est, permutations = 100) >summary(Chao2, display = “chao”) Os comandos specpool e diversityresult são mais simples e diretos, pois apresentam somente o valor final estimado: >Chao2 <- specpool(est) >Chao2 >Chao2 <- diversityresult(est, index = “chao”) JACKKNIFE 1 Este estimador baseia-se no número de espécies que ocorrem em somente uma amostra (Q1).   𝐽𝑎𝑐𝑘! = 𝑆!"# + 𝑄!   𝑚 − 1 𝑚 Onde: m = número de amostras     75  𝛾!"#! = 𝑚𝑎𝑥 𝑆!"!# 𝐶!"# 𝑖(𝑖 − 1)𝐹!!"!!! (𝑁!"!#)(𝑁!"!# − 1) − 1   𝐶!"# = 1 + 𝐹! 𝑁!"!#   𝑁!"!# = 𝑖𝐹! !" !!! Não precisa fazer cara feia, é óbvio que iremos usar o programa para fazer esses cálculos. REALIZANDO O EXERCÍCIO NO PROGRAMA R: Comandos >est <- read.table(“estimadores.txt”, h = T) >ACE <- estaccumR(est, permutations = 100) >summary(ACE, display = “ace”) Outra maneira de conseguir o mesmo valor: >est1<-colSums(est) ## soma abundância de cada linha= abundância total por espécie >ACE <- estimateR(est1) >ACE   ICE (Incidence-based Coverage Estimator) Este método trabalha com o número de espécies infreqüentes (que ocorrem em poucas unidades amostrais). Esse método permite ao pesquisador determinar os limites para os quais uma espécie seja considerada infreqüente. Em geral, são consideradas como tal espécies com incidência entre 1 e 10 indivíduos (Chazdon et al. 1998) ou 1 a 20 (Walther & Morand 1998). A riqueza estimada pode variar conforme se aumente ou diminua o limiar de incidência, e     76  infelizmente não existem critérios biológicos definidos para a escolha do melhor intervalo (Santos 2003).   𝐼𝐶𝐸 = 𝑆!"#$ +   𝑆!"# ! 𝐶!"# + 𝑄! 𝐶!"#  𝛾!"#! onde:   𝛾!"#! = 𝑚𝑎𝑥 𝑆!"# ! 𝐶!"# 𝑚!"# ! (𝑚!"# !!!) 𝑗(𝑗 − 1)𝑄!!"!!! (𝑁!"# !)! − 1   𝐶!"# = 1 + 𝑄! 𝑁!"# !   𝑁!"# ! = 𝑗𝑄! !" !!!   REALIZANDO O EXERCÍCIO NO PROGRAMA R: Comandos >est <-read.table(estimadores, h = T) >ICE <- poolaccum(est, permutations = 100) >summary(ICE, display = “ice”) Outra maneira de conseguir o mesmo valor: >ICE <- specpool(est) >ICE   BOOTSTRAP Este método difere dos demais por utilizar dados de todas as espécies coletadas para estimar a riqueza total, não se restringindo às espécies raras. Ele requer somente dados de     77  incidência. A estimativa pelo bootstrap é calculada somando-se a riqueza observada à soma do inverso da proporção de amostras em que cada espécie ocorre. 𝐵𝑜𝑜𝑡 = 𝑆!"# +   (1 − 𝑃!)! !!"# !!! Onde: Pk = proporção do número de amostras em que cada espécie foi registrada m = número de amostras Exemplo: Usando os dados da tabela 1 calcule o valor de bootstrap para a comunidade: Bootstrap = 14 + [ (1- 8/14)14 +(1- 2/14)14 +(1- 10/14)14 +(1- 10/14)14 +(1- 3/14)14 +(1- 3/14)14 +(1- 2/14)14 + (1- 7/14)14 +(1- 5/14)14 +(1- 1/14)14 +(1- 5/14)14 +(1- 2/14)14 +(1- 14/14)14 +(1- 1/14)14] Bootstrap = 14 + 1 ,127 Boostrap = 15,127 REALIZANDO O MESMO EXERCÍCIO NO PROGRAMA R: Comandos >est <-read.table(estimadores, h = T) >BOOT <- poolaccum(est, permutations = 100) >summary(BOOT, display = “boot”) Outra maneira de conseguir o mesmo valor: >BOOT <- specpool (est) >BOOT >BOOT <- diversityresult (est, index = “boot”)     80  9 – Agora é necessário configurar o programa para realizar os testes; 10 – Clicar em “DIVERSITY”, como demonstrado na tela abaixo; 11 – Escolham a opção “Diversity Settings” 12 – Coloque 500 no lugar de 50 aleatorizações 13 – Depois de colocarem 500 cliquem na aba Estimators (destacado em amarelo) e depois em OK;     81  14 - Determine o número de espécies raras para o ACE e ICE. Esse número corresponde ao número de espécies que o programa irá considerar como espécies raras; 15 – Clicar em OK; 16 - Agora é só correr o teste. Clicar em Compute Diversity Stats; 17 – Aparecerá uma tela com os resultados do teste; 18 – Clicar em Export e salvar em algum lugar no seu computador, depois é só abrir com o Excel e fazer os gráficos no R;     82  Índices de diversidade Os índices de diversidade representam uma medida que combina a riqueza e abundância relativa (equitabilidade) das espécies de uma comunidade. O índice de Shannon (H’) é um dos mais utilizados na literatura para medir a diversidade de espécies. Este índice é derivado da teoria da informação e sua função foi derivada como: H’ = − 𝑝! ln 𝑝! Onde pi representa a proporção de indivíduos na i-nésima espécie em relação à abundância total na comunidade. Quanto maior o valor de H’, maior a diversidade da comunidade. Os valores de H’ raramente ultrapassam 4, sendo que para que H’ seja maior do que 5 a comunidade precisa ter mais de 105 espécies. Um dos problemas do índice de Shannon é que a diversidade é confundida pela riqueza de espécies e equitabilidade. Desse modo, tanto o número de espécies quanto o esforço amostral afetam o valor final do índice. Além disso, quando confrontamos valores de diversidade entre duas comunidades, por exemplo, H’ = 2,71 e H’ = 2,59, temos dificuldade para decidir se os valores são, de fato, diferentes. Outro índice de diversidade muito usado por ecólogos é o índice de Simpson (D). Este índice mede a probabilidade de dois indivíduos coletados ao acaso pertencerem à mesma espécie através da fórmula: D = 𝑝!! Onde pi representa a proporção de indivíduos na i-nésima espécie em relação às abundância total na comunidade. Quanto maior o valor de D, menor a diversidade da comunidade. Alguns autores expressam a fórmula do índice de Simpson como 1 – D ou 1 / D. Este índice é considerado uma das medidas de diversidade mais robustas. Apesar de existir um número impressionante de métricas para medir a diversidade biológica (Hulbert 1972, Magurran 2004), diversos autores desencorajam o uso dessas métricas para testar hipóteses ecológicas. Dentre os principais motivos destacamos: (1) ausência de uma base probabilística que nos permita assinalar valores de significância que, por sua vez, impede que façamos comparações biológicas entre duas comunidades; (2) todos os índices de diversidade são fortemente sensíveis ao número de indivíduos e de espécies; (3) problemas conceituais e de múltiplas definições que trazem pouco sentido biológico e dificultam a interpretação de padrões ecológicos. Dentre os autores que criticam a utilização de índices de diversidade na ecologia, se destacam pela clareza dos argumentos o trabalho marcante de ÍNDICES DE DIVERSIDADE E DIVERSIDADE BETA (β)       85  Onde N representa o desvio Normal e µ e σ são os coeficientes da fórmula. A abundância esperada (BSar) para a espécie na ordem (do inglês “rank”) r para o modelo Broken- Stick é: 𝐵𝑆𝑎! = ( 𝐽/𝑆 ) (1/𝑥 ) ! !!! Onde J representa o número total de indivíduos na comunidade e S o número total de espécies. Para o modelo Série Geométrica, a abundância esperada (GSar) para a espécie da ordem r é: 𝐺𝑆𝑎! = 𝐽𝛼(1 − 𝛼)!!! Onde J representa o número total de indivíduos na comunidade e o coeficiente α é uma estimativa da taxa de decréscimo da abundância por ordem r. Para o modelo Zipf, a abundância esperada (Zar) para a espécie da ordem r é: 𝑍𝑎! = 𝐽𝑝!𝑟! Onde J representa o número total de indivíduos na comunidade, p1 é a proporção ajustada da espécie mais abundante e γ é o coeficiente de decréscimo da abundância por ordem r. O modelo Zipf-Mandelbrot acrescenta um parâmetro na fórmula do Zipf para estimar a abundância (ZMar) da espécie da ordem r: 𝑍𝑀𝑎! = 𝐽𝑐(𝑟 + 𝛽)! Onde J representa o número total de indivíduos na comunidade, c e β são constantes de escala e γ é o coeficiente de decréscimo da abundância por ordem r (Wilson 1991). 0 20 40 60 80 100 120 140 160 1 2 3 4 5 10 20 40 60 N úm er o  de  e sp éc ie s Número de indivíduos 0 10 20 30 40 50 60 70 80 90 100 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Ab un dâ nc ia Ordem das espécies Comunidade A Comunidade B Comunidade C Figura 9. Duas representações comuns do padrão de distribuição da abundância das espécies. (A) Representação básica com o número de espécies com suas respectivas abundâncias organizadas em ordem decrescente. (B) Diagramas de abundância relativa (ou curvas de dominância) que podem ser utilizados para comparar o padrão de dominância entre diferentes comunidades. A) B)     86  Escolhendo o melhor modelo teórico no R > library(vegan) > rios=read.table("rios.txt", h=T) > rios > rad.rio1=radfit(rios[1,]) > rad.rio1 > plot(rad.rio1, xlab="Ordem das espécies", ylab="Abundância", pch=19) > rad.rio2=radfit(rios[2,]) > rad.rio2 > plot(rad.rio2, xlab="Ordem das espécies", ylab="Abundância", pch=19) > rad.rio3=radfit(rios[3,]) > rad.rio3 > plot(rad.rio3, xlab="Ordem das espécies", ylab="Abundância", pch=19) > par(mfrow=c(2, 2)) > plot(rad.rio1, main="Rio 1", xlab="Ordem das espécies", ylab="Abundância", pch=19) > plot(rad.rio2, main="Rio 2", xlab="Ordem das espécies", ylab="Abundância", pch=19) > plot(rad.rio3, main="Rio 3", xlab="Ordem das espécies", ylab="Abundância", pch=19) Praticando: Exercício 1: A bióloga responsável pela Secretaria de Meio Ambiente do Município de Florianópolis/SC precisa determinar a qualidade da água das seis praias mais movimentadas da cidade. Este trabalho surgiu após reclamações de banhistas e de pescadores de algumas dessas praias. A bióloga mediu os níveis de colifórmes fecais e coletou peixes em vários pontos de cada praia. Um estagiário derrubou o computador da bióloga e perdeu todos os dados dessa pesquisa. Por sorte, a bióloga havia anotado todos os dados referentes aos peixes coletados nas praias. Porém, os dados sobre os níveis de colifórmes fecais só foram anotados em arquivo digital. Com recursos limitados, a bióloga não pôde refazer as análises da qualidade da água e precisa realizar uma avaliação indireta a partir dos dados de riqueza e abundância de peixes. Teoria: Teoria do distúrbio + Distribuição da Abundância das Espécies (SADs)     87  Pergunta: Praias mais poluídas possuem padrão de distribuição da abundância da espécies mais equitativo? Unidade amostral: Pontos de amostragem em cada praia Variável dependente: Abundância relativa Variável independente: Praia Importe a planilha “peixes.floripa.txt” e indique a partir dos diagramas de abundância relativa qual a praia com melhor e pior qualidade da água. Informe os modelos teóricos que melhor explicam o padrão de distribuição de abundância de cada praia e faça um diagrama de abundância relativa para cada praia e uma figura contendo todos os diagramas na mesma janela. Diversidade beta Desde o início da ecologia, a identidade das espécies que constituem determinada comunidade (i.e., composição de espécies) tem gerado uma série de hipóteses importantes para o entendimento de como os organismos se distribuem no espaço e no tempo. Uma das principais perguntas sobre esse assunto é “O que torna comunidades de espécies mais ou menos similares em diferentes lugares e tempos?” (Vellend 2010). Após os influentes estudos do ecólogo Robert Whittaker (Whittaker 1960, 1972), o termo diversidade beta (i.e., variação na composição de espécies entre áreas) ganhou força na literatura ecológica. Nas duas últimas décadas, o número de trabalhos aumentou expressivamente com o desenvolvimento de novos métodos para medir a diversidade beta e de novos pacotes estatísticos. A grande quantidade de medidas, abordagens estatísticas, termos e interpretações para a diversidade beta aumentaram a confusão em relação às maneiras corretas de acessar e testar os padrões de modificação na composição de espécies (Tuomisto 2010a,b, Anderson et al. 2011). Nesta apostila utilizaremos um roteiro prático baseado em hipóteses sugerido recentemente por Anderson et al. (2011). Primeiro, é importante diferenciar dois tipos de conceito de diversidade beta, o conceito de substituição (turnover) e de variação. A substituição representa a modificação na composição de espécies de uma unidade amostral para a outra ao longo de um gradiente espacial, temporal ou ambiental. A substituição requer um gradiente que indique direção como, por exemplo, investigar a mudança na composição de espécies ao longo de um gradiente de profundidade em um lago (Fig. 10a). As principais questões testadas na análise de substituição são: (1) quantas novas espécies são encontradas ao longo de um gradiente e quantas delas foram inicialmente presentes e agora foram perdidas? (2) Qual a proporção de espécies encontradas em uma unidade amostral que não são compartilhadas com a próxima unidade do gradiente?     90  Onde y1j representa a abundância da espécie j na localidade x1 e y2j na localidade x2. Esse cálculo prossegue até a espécie p. Medidas multivariadas Uma medida de diversidade beta interessante para comparar N amostras é a dispersão em um espaço multivariado, com uma análise conhecida como teste de homogeneidade de dispersões multivariadas (Anderson 2006). Esta análise calcula o centróide (ou mediana especial) de um grupo específico (e.g., lagoa 1) e compara a dissimilaridade média das n observações individuais dentro desse grupo (e.g., abundância de cada espécie p na lagoa 1) utilizando uma medida apropriada de dissimilaridade (e.g., Bray-Curtis, Chao-Sørensen, Distância Euclideana, Jaccard, Sørensen). O cálculo do centróide para medidas que utilizam distância euclidiana é a média aritmética de cada variável. Porém, para calcular o centróide para índice de distância não-euclidianos (e.g., Jaccard) é necessário fazer uma análise de coordenadas principais (Anderson 2006). A hipótese nula desta análise é a de que a diversidade beta não é diferente entre as amostras de interesse. Para acessar a probabilidade de a hipótese nula ser verdadeira utiliza-se a estatística F de Levene comparando a distância média de cada observação ao centróide do seu grupo que, por sua vez, é definido por uma medida de dissimilaridade. Para gerar os valores do P são realizadas n permutações (e.g., 1000) (detalhes em Anderson 2006). Calculando os índices de diversidade no R 1. Calculando o índice clássico de Whittaker (βw): > salinidade=read.table("salinidade.txt", header=T) > salinidade > diversidade.beta=betadiver(salinidade, "w") > diversidade.beta 2. Calculando índices de Jaccard e Sørensen: > jaccard=betadiver(salinidade, "j") > sorensen=betadiver(salinidade, "sor") > scores(jaccard) > scores(sorensen)     91  3. Calculando os índices de Bray-Curtis e Morisita-Horn: > library(vegan) > data(mite) > bray=vegdist(mite, "bray") > bray > morisita.horn=vegdist(mite, "horn") > morisita.horn # Testando hipóteses com as matrizes de similaridade/dissimilaridade > library(vegan) > data(varespec) > data(varechem) > dist.species=vegdist(varespec, "bray") > dist.chemical=vegdist(scale(varechem), "euclidean") > associacao=mantel(dist.species, dist.chemical) > associacao 4. Calculando os índices de Chao-Jaccard e Chao-Sørensen: > CSoren.dist=ecol.dist(ilhas, chao.sorenson, type="dis") > CSoren.simi=ecol.dist(ilhas, chao.sorenson, type="sim") > CJaccar.dist=ecol.dist(ilhas, chao.jaccard, type="dis") > CJaccar.simi=ecol.dist(ilhas, chao.jaccard, type="sim") # se optar por calcular a similaridade entre duas localidades use a seguinte função: > IlhaA=ilhas[,1] > IlhaB=ilhas[,2] > CSoren.A.B=chao.sorenson(IlhaA, IlhaB) > CJaccar.A.B=chao.jaccard(IlhaA, IlhaB) > CSoren.A.B > CJaccar.A.B 5. Calculando outros índices de similaridade com o pacote fossil: > library (fossil) > Comunidade.A <- c(1,0,4,3,5,0,0,7) > Comunidade.B <- c(2,1,3,0,0,1,0,6) > bray.curtis(Comunidade.A, Comunidade.B)     92  > jaccard(Comunidade.A, Comunidade.B) > simpson(Comunidade.A, Comunidade.B) > sorenson(Comunidade.A, Comunidade.B) > morisita.horn(Comunidade.A, Comunidade.B) 6. Teste de homogeneidade de dispersões multivariadas: > library(vegan) > cafe=read.table("cafe.txt", header=T) > tipo.matriz=factor(c(rep(1,16), rep(2,8)), labels = c("com.mata","sem.mata")) > dissimilaridade=vegdist(cafe, "bray") > HDM=betadisper(dissimilaridade, tipo.matriz) > valor.P=permutest(HDM, pairwise = F) > plot(HDM) Praticando: Exercício 1: Baseado na teoria de que os organismos selecionam sua planta hospedeira considerando características fisiológicas e estruturais, um biólogo pretende testar se três clones (clones x1, x2, e x3) de uma planta X possuem composição de espécies de ácaros diferente. Ele coletou ácaros em 60 plantas (20 plantas de cada clone) em uma estação experimental que cultiva a planta X. Em cada planta, o biólogo coletou 10 folhas e identificou e quantificou todos os ácaros. Além disso, o biólogo mensurou o comprimento, largura e área foliar e a densidade de tricomas. Pergunta: O clone afeta a composição de espécies de ácaros? Teoria: Teoria do nicho (species sorting) Unidade amostral: Folha Variável dependente: Composição de espécies Variável independente: comprimento, largura e área foliar, e a densidade de tricomas. Importe a planilha “clone.col1.txt” e “clone.col2.txt”e verifique se os clones possuem composição semelhante ou diferente nas duas coletas hipotéticas. Após as análises, responda a pergunta do biólogo para cada coleta. Os resultados realmente permitem que a pergunta seja respondida? O que você pode interpretar com a coleta 1 e com a coleta 2?     95  Agrupamento Análise de agrupamento hirerárquico (cluster) A análise de agrupamento hierárquico é a mais utilizada em ecologia. No entanto, existem também outras análises não hierárquicas, como a K-means, que não serão abordadas neste curso. O objetivo da análise de agrupamento é agrupar objetos admitindo que haja um grau de similaridade entre eles. Esta análise pode ser utilizada ainda para classificar uma população em grupos homogêneos de acordo com uma característica de interesse. A grosso modo, uma análise de agrupamento tenta resumir uma grande quantidade de dados e apresentá- la de maneira fácil de visualizar e entender (em geral, na forma de um dendrograma). No entanto, os resultados da análise podem não refletir necessariamente toda a informação originalmente contida na matriz de dados. Para avaliar o quão bem uma análise de agrupamento representa os dados originais existe uma métrica — o coeficiente de correlação cofenético — o qual discutiremos em detalhes mais adiante. Apesar da sua versatilidade, deve-se ressaltar que nem todos os problemas em ecologia são problemas de agrupamento. Antes de considerar algum método de agrupamento, pense porque você esperaria que houvesse uma descontinuidade nos dados; ou ainda, considere se existe algum ganho prático em dividir uma nuvem de objetos contínuos em grupos. Além disso, existem algumas críticas que merecem atenção: mesmo para um conjunto de dados aleatórios é possível encontrar grupos; o padrão apresentado pelo dendograma depende do protocolo utilizado (método de agrupamento e índice de dissimilaridade); os grupos formados dependem do nível de corte escolhido. Normalmente, a análise de agrupamento tenta arranjar os objetos em grupos que são mutuamente excludentes, ou seja, o mesmo objeto não pode fazer parte de mais de um grupo. No entanto, existem algumas técnicas, chamadas de fuzzy clustering, que permitem uma gradação na classificação de objetos. Esta técnica não será abordada neste módulo, mas o leitor interessado é remetido à duas referências: Legendre & Legendre (1998) e Borcard et al. (2011). Os passos para a análise de agrupamento são os seguintes: 1) A matriz deve conter os objetos a serem agrupados (p.ex. espécies) nas linhas e as variáveis (p.ex., locais de coleta ou medidas morfológicas) nas colunas. Primeiramente, se os dados forem de abundância, é mais correto realizar a transformação de Hellinger (Legendre & Gallagher, 2001). Se a matriz original contiver muitos valores     96  discrepantes (p.ex., uma espécie muito mais ou muito menos abundante que outras) é necessário transformar os dados usando Log (x+1)1. Se as variáveis forem medidas tomadas em diferentes escalas (metros, graus celcius etc), é necessário padronizar cada variável utilizando a seguinte fórmula: desvio médiaobsZ −= Onde obs representa o valor da unidade amostral de interesse e os valores da média e do desvio padrão são calculados para cada variável. 2) Escolha do método de agrupamento A escolha do método de agrupamento é crítico para a escolha de um coeficiente de associação. É importante compreender completamente as propriedades dos métodos de agrupamento para interpretar corretamente a estrutura ecológica que eles evidenciam (Legendre & Legendre, 1998). De acordo com a classificação de Sneath & Sokal (1973) existem cinco tipos de métodos: 1) seqüenciais ou simultâneos; 2) aglomerativo ou divisivo; 3) monotéticos ou politéticos; 4) hierárquico ou não hierárquicos e 5) probabilístico. Por motivos de espaço e tempo discutiremos somente os métodos hierárquicos, que são os mais comumente encontrados na literatura ecológica. Métodos hierárquicos podem ser divididos naqueles que consideram o centróide ou a média aritmética entre os grupos. O principal método hierárquico que utiliza a média aritmética é o UPGMA (Agrupamento pelas médias aritméticas não ponderadas), e o principal método que utiliza centróides é a Distância mínima de Ward. O UPGMA funciona da seguinte forma: a maior similaridade (ou menor distância) identifica os próximos agrupamentos a serem formados. Após esse evento, o método calcula a média aritmética das similaridades ou distâncias entre um objeto e cada um dos membros do grupo ou, no caso de um grupo previamente formado, entre todos os membros dos dois grupos. Todos os objetos recebem pesos iguais no cálculo. A matriz de similaridade ou distância é atualizada e reduzida de tamanho em cada etapa do agrupamento, por isso não exige tanto do computador (Legendre & Legendre, 1998).                                                               1 O uso do 1 é obrigatório pois Log de zero na base 10 não existe.       97  O método de Ward é baseado no critério de quadrados mínimos dos modelos lineares. O objetivo é definir os grupos de maneira que a soma de quadrados (i.e. similar ao erro quadrado da ANOVA) dentro dos grupos seja minimizada (Borcard et al. 2011). 3) Escolha dos índices de similaridade (coeficientes de distância ou de associação, ou índices de dissimilaridade). Os índices de similaridade medem a distância entre dois objetos ou quantificam o quanto eles são parecidos. Lembre-se: as questões e hipóteses iniciais do estudo devem ser levadas em conta na escolha do índice (veja Anderson et al. 2011). Índices binários assimétricos Se os dados disponíveis foram de presença-ausência (binários), os índices recomendados são os de Jaccard e Sørensen. Os índices tradicionais de Jaccard e Sørensen são chamados de índices assimétricos, pois ao fazerem a comparação entre amostras não levam em conta duplas ausências. Essa característica é desejável ao analisar dados ecológicos porque o não encontro de duas espécies em duas localidades não é um indicativo de que duas localidades sejam similares, já que isto pode ter surgido por variação estocástica na amostragem, padrões de dispersão, etc. Além disso, as duplas-ausências não refletem necessariamente diferenças nas localidades (Legendre & Legendre, 1998; Anderson et al., 2011). Desta forma, somente serão considerados similares localidades que de fato compartilhem espécies. Compare as fórmulas dos coeficientes de Jaccard e Sørensen (Pag. 89): Como é possível perceber pelas fórmulas, o coeficiente de Sørensen dá um peso maior para as duplas presenças, pois elas são um indicativo mais forte de semelhança. No entanto, o índice de Sørensen é sensível à variações na riqueza entre as localidades. Como uma alternativa, o índice de Simpson para similaridade múltipla entre comunidades foi proposto recentemente por Baselga et al. (2007) como uma modificação do índice de diversidade de Simpson. Este índice tem a vantagem de ser independente da riqueza e assim, consegue distinguir entre a substituição verdadeira e a simples perda de espécies. Isto é importante porque, como visto anteriormente, a diversidade beta pode ser causada por dois distintos fenômenos: aninhamento e substituição de espécies (turnover) que, por sua vez, são causados por processos ecológicos diferentes. Além disso, este índice leva em consideração a similaridade em toda comunidade e não par-a-par, como outros índices tradicionais (Baselga et