Resumo

O padrão espacial de doenças é estudado para se conhecer o mecanismo de atuação da doença e para orientar estratégias de controle. A distribuição beta-binomial é uma generalização da distribuição binomial quando assume-se que a incidência não é constante na região em estudo. Ela pode ser escrita de forma que θ seja um parâmetro de agregação espacial da doença. Porém, não há solução analítica para a estimaço dos parâmetros. Dessa forma, devem ser utilizados procedimentos iterativos de maximização numérica da função de verossimilhança. É demonstrado o uso das funções optimize(), optim(), nlm() e nlminb() disponíveis em R para maximização numérica e das funções fitdistr() e mle() disponíveis para estimação por máxima verossimilhança. A metodologia foi aplicada a dados simulados, dados de Pinta Preta dos Citrus e de várias avaliações em tempos diferentes de Morte Súbita dos Citrus. Os resultados foram similares aos métodos tradicionais de análise por quadrats.

1 Introdução

Na análise do padrão espacial da incidência de doenças em plantas, costuma-se fazer a análise por quadrat counts, em que é estudado o número de plantas doentes por quadrat, Y . Os quadrats são áreas retangulares que abrangem n plantas em l linhas e c plantas na linha. Se o padrão espacial da doença no talhão é aleatório, pode-se assumir que Y ~ Bin(n,π), com n o número de plantas por quadrat e π a taxa de incidência da doença no talhão. Se n →∞ e π → 0, pode-se utilizar a distribuição de Poisson(θ), assumindo que θ ≈ nπ.

O teste da hipótese de aleatoriedade espacial nesse caso, consiste em verificar se é razoável assumir uma distribuição binomial ou de Poisson para o número plantas doentes. No contexto de incidência de doenças, Madden & Hughes (1995) sugerem que um índice de dispersão D seja calculado para medir a agregação e para testar superdispersão. D é a razão entre a variância observada e a variância esperada segundo um modelo. Se a incidência tem distribuição espacial aleatória, D será igual a 1 e pode-se assumir que a contagem segue uma distribuição binomial. Se a distribuição espacial for regular, D é menor que 1. Se o padrão for agregado, D é maior que 1.

A variância observada de Y é dada por

(1)

em que y_i é o número de plantas doentes no i-ésimo quadrat, n é o número de plantas por quadrats, p é a proporção de plantas doentes no talhão e N o número de quadrats avaliados. A variância esperada considerando um modelo Binomial é dada por

(2)

Nota-se que n é fixo, impondo uma restrição, ou seja, o número de plantas por quadrats devem ser iguais. Diante disso, pode-se utilizar a variância esperada segundo o modelo de Poisson, que é a média, podendo ser calculada por

(3)

em que n_i é o número de plantas no i-ésimo quadrat. Para o teste de significância para agregação utiliza-se a (N - 1)D, que tem distribuição χ_N-1².

Nos modelos binomial e de Poisson, π é constante em toda a região em estudo e o índice de dispersão é uma maneira simplista de testar isso. Em situações em que π varia, há evidências para afirmar que existe agregação espacial. Para considerar a variação de π, pode-se assumir uma distribuição Beta(α,β). A variação de π pode ser dada por α e β, ou então, pode ser feita uma reparametrização de forma que a variação de π seja expressa em um único parâmetro θ. No contexto de doenças de plantas é de grande importância, pois permite concluir se o padrão é agregado e fornece uma medida da agregação.

Na Seção 2 são implementadas as funções de densidade, de distribuição, de quantis e geradora de números aleatórios. Na Seção 3 são apresentadas as funções disponíveis em R para maximização numérica e para estimação por máxima verossimilhança. Na Seção 5 são feitas aplicações a dois conjuntos de dados de incidência de doenças em plantas.

2 Distribuição beta-binomial

A distribuição beta-binomial é obtida, segundo Skellam(1948) apud Mendonça, Bearzoti, Aquino &MFerreira (2004), compondo-se a distribuição binomial com a densidade beta para π. A função de densidade beta pode ser escrita de várias formas, uma delas é:

(4)

0 ≤ π ≤ 1, α e β são constantes reais positivas e Γ(.) é a função gama.

A função de densidade de probabilidade da distribuição beta-binomial pode ser expressa por:

(5)

com y = 0, 1,...,n; α ≥ 0, β ≥ 0 conhecidos e B a função beta definida por B(α,β) = .

Griffiths (1973) sugere uma reparametrização implicando, segundo Mendonça et al. (2004), em simplicidade e interpretabilidade dos parâmetros. Os novos parâmetros são p = , o valor esperado do parâmetro π da distribuição binomial e θ = , uma medida da variação em π.

O valor esperado de Y e sua variância são dados por: E(Y ) = np e V ar(Y ) = . Assumindo-se que p varia (θ > 0), pode-se observar que a variância da distribuição beta-binomial é maior que da binomial. O parâmetro θ, no contexto dos dados de incidência de doenças, pode ser interpretado como um índice de agregação espacial da doença. Se θ > 0, o padrão da doença é agregado.

2.1 Definindo a função de densidade

Pode-se calcular a densidade de probabilidade definindo uma função em R. A forma mais fácil de implementar é a Equação 5. Para manter a interpretabilidade dos parâmetros no contexto de doenças em plantas, a função terá p e θ como argumentos e define internamente α = e β = .

Um detalhe que se deve acrescentar é que para valores de θ muito próximos de zero, os parâmetros α e β são muito grandes. Isso faz com que o valor da função beta seja muito próximo de zero. Exemplo:

Os valores de α e β não muito grandes, fazem a função se aproximar muito de zero. A escala logarítmica ameniza este problema:

A escala logarítmica é a utilizada no cálculo da densidade:

A partir desta função pode-se estudar a densidade de probabilidade em gráficos de linhas, plotando o valor da densidade variando os parâmetros da distribuição. Na Figura 1, visualiza-se a densidade para algumas variações dos parâmetros.

Observa-se na Figura 1 que quanto mais próximo θ está de zero, a densidade da distribuição beta-binomial fica mais próxima da densidade da distribuição binomial. Além disso, valores de θ, próximos de zero comprometem a precisão no cálculo da densidade:

A precisão pode ser melhorada considerando que a distribuição beta-binomial se aproxima da distribuição binomial quando θ ≈ 0. Assim, pode-se definir a função de densidade de forma que ser o valor de θ for muito próximo de zero, seja tomada a densidade da distribuição binomial. No R, .Machine é uma lista com as características numéricas do computador utilizado. .Machine[[’double.eps’]] é o menor número positivo tal que o computador em uso seja capaz de reconhecer como verdadeiro que 1 + .Machine[[′double.eps′]] > 1:

A função de densidade pode ser re-definida, usando a raiz de .Machine[[’double.eps’]] como valor mínimo para θ:

2.2 Definindo as funções de distribuição, de quantis e de números aleatórios

A partir da função de densidade, podem ser definidas também a função de probabilidade acumulada, sua inversa e a geradora de números aleatórios.

A função de distribuição, P(Y ≤ y) é obtida somando-se as densidades avaliadas em {0,..,y}.

No cálculo do inverso da função de distribuição acumulada, ou quantis, a partir de {p (0, 1)}, encontra-se o menor valor de y = {0,...,n} tal que p > P(Y ≤ y).

A função geradora de números aleatórios para qualquer densidade de probabilidade pode ser genericamente obtida a partir de um gerador de números aleatórios da distribuição uniforme no intervalo (0, 1) e da função que calcula quantis. Dessa forma, pode-se gerar amostras da distribuição U(0, 1) e calcular quantis tomando essas amostras como probabilidade acumulada.

No caso de distribuições cujo número de valores possíveis é finito, o uso da função acumulada inversa pode ser evitado. É possível utilizar uma função para sortear valores a partir de um conjunto finito de valores, com probabilidade dada pela função de densidade. No R, pode-se utilizar a função sample() para isso.

O pacote rmutil, de J. K. Lindsey, contém essas funções. Porém, há um inconveniente apontado na Seção 3 que impede utilizá-las. Esse pacote, não está disponível no cran, mas está disponível em http://alpha.luc.ac.be/~jlindsey/rcode.html. Na parametrização adotada por Lindsey, é utilizado ϕ = 1∕θ. Nesse caso ϕ é um parâmetro de dispersão de π.

3 Funções para maximização numérica

No procedimento de estimação por máxima verossimilhança é necessário encontrar as derivadas em relação aos parâmetros. No caso da estimação dos parâmetros da distribuição beta-binomial, não há expressão analítica fechada para os estimadores de máxima verossimilhança. Nessa situação, a solução que maximiza a verossimilhança deve ser encontrada iterativamente, via algoritmo de maximização numérica. Devido à estabilidade numérica, é conveniente maximizar o logaritmo da função de verossimilhança.

Para a maximização numérica, será utilizada a função dbetabinom() definida anteriormente, embora já esteja implementada no pacote rmutil. A função dbetabinom() do pacote rmutil apresenta um inconveniente que impossibilita utilizá-la na função de maximização numérica: quando os valores dos parâmetros estão fora do espaço paramétrico, a função é interrompida e não retorna valor algum. Na função dbetabinom() definida anteriormente, ao ser informado um valor fora do espaço paramétrico, a expressão é avaliada e é retornado zero ou NA, dependendo do caso. Isto permite utilizá-la no procedimento de maximização numérica.

No R, estão disponíveis funções para minimização ou maximização numérica a partir de uma função objetivo definida. A menos que seja informado em argumentos adicionais dessas funções, estas buscarão o valor que minimiza a função objetivo. Neste caso, para simplificar, as mesmas são utilizadas para encontrar a solução que minimiza o negativo do logaritmo da função de verossimilhança, (R Development Core Team 2005).

3.1 Função optimize()

A função optimize() ou seu pseudônimo optimise(), pode ser utilizada para encontrar o mínimo de uma função em um intervalo. Esta função usa uma tradução em C do código Fortran (Netlib) baseada no procedimento do ’Algol 60 localmin’ dado na referência.

Os argumentos obrigatórios da função optimize() são: a função e o intervalo. Os demais argumentos são alternativos. O argumento maximum indica se o objetivo é maximizar ou minimizar.

Para exemplificar, seja uma amostra da distribuição binomial com n = 5 e π = 0.5.

Suponha desconhecida a expressão analítica do estimador de máxima verossimilhança de π. Pode-se maximizar numericamente a função de verossimilhança em relação a π e, para isto, deve-se defini-la de forma que o parâmetro a ser estimado seja o primeiro argumento:

A função vero.binom() retorna a soma do logaritmo das densidades avaliadas em um valor de x informado, considerando os parâmetros π e n também informados.

O valor estimado para π foi 0.4899824 e o correspondente valor do logaritmo da verossimilhança de -146.2198. Conhecendo a expressão analítica do estimador, pode-se confirmar:

3.2 Função nlm()

A função nlm() pode ser utilizada para minimização ou maximização numérica de funções não-lineares. O procedimento iterativo é um algoritmo do tipo de Newton. Maiores detalhes podem ser vistos na documentação.

A função a ser minimizada deve ser informada no argumento f. O gradiente e o hessiano da função podem ser informados como atributos da função. Quando estes argumentos não são informados, são utilizadas derivadas numéricas da função para a minimização.

A expressão do negativo do logaritmo da função de verossimilhança, considerando a distribuição beta-binomial pode ser implementada por:

Para simular 100 amostras com distribuição beta-binomial com parâmetros π = 0.3, n = 5 e θ = 0.1 utiliza-se a função rbetabinom() definida anteriormente:

Os códigos 1 e 2 no elemento code indicam que, provavelmente, chegou-se à solução.

A matriz de informação de Fischer pode ser aproximada pelo hessiano numérico obtido numericamente, neste caso.

3.3 Função optim()

A função optim() permite obter os valores de parâmetros de uma função que minimizam essa função. Os métodos de optimização disponíveis (R versão 2.1.1) são: Nelder-Mead, Quasi-Newton, Gradiente-Conjugado, Confinamentos de Caixa e Simulated Annealing. Na documentação há detalhes desses métodos.

Nesta função, os argumentos obrigatórios são: o valor inicial dos parâmetros e a expressão da função. Os argumentos opcionais são: a função gradiente, o método de minimização (na documentação há detalhes sobre cada método disponível), os limites inferior e superior do espaço paramétrico se for utilizado o método ’L-BFGS-B’, uma lista de parâmetros de controle e demais argumentos da função a ser minimizada ou da função gradiente.

Utiliza-se a função optim() para a estimação de π e θ. Neste caso o valor verdadeiro de θ é zero. Pode-se utilizar o método de otimização default:

Os valores retornados foram: os valores dos parâmetros com os quais encontrou-se o mínimo da função; o correspondente valor da função avaliado nos parâmetros encontrados; o número de valores avaliados na função e no gradiente, neste caso não foi avaliado o gradiente; o código indicando a convergência, zero indica que convergiu; e, alguma mensagem adicional quando houver.

Apresenta-se a seguir as estimativas de máxima verossimilhança da amostra 2. Obtém-se, também, o hessiano numérico para aproximar a matriz de informação de Fisher observada colocando o argumento hessian=TRUE.

A matriz de informação de Fisher aproximada é obtida fazendo:

3.4 Função nlminb()

A função nlminb() faz minimização numérica usando rotinas PORT. Também realiza minimização sob restrições no espaço paramétrico. Exemplo, com restrição no espaço paramétrico:

4 Funções para estimação por máxima verossimilhança

A função fitdistr() e mle() são funções especificamente implementadas para estimação por máxima verossimilhança.

4.1 Função fitdistr()

A função fitdistr() do pacote MASS é específica para estimação por máxima verossimilhança (Venables & Ripley 2002).

Para utilizá-la, basta entrar com os dados, a função de densidade e os valores iniciais em formato de lista nomeada. Esta função utiliza a função optim() e demais argumentos podem ser passados para função optim() e também para a função de densidade nomeando-os, por exemplo, o argumento size. Essa função retorna também as estimativas dos erros padrão, obtidos através do hessiano numérico.

Na amostra 1 não foi possível utilizar a função fitdistr(). Sabe-se que o valor do parâmetro θ é zero e é o limite inferior do espaço paramétrico. A dificuldade aqui é encontrar o hessiano numérico, pois para isso é preciso avaliar a função na região próxima ao valor estimado na minimização numérica e para qualquer valor de θ menor que zero, a função de densidade não é definida. Neste caso, visualiza-se apenas a função desvio (ou deviance), Figura 2.

Aplicando às outras duas amostras, podemos utilizar o método “L-BFGS-B” que obrigatoriamente necessita a definição do espaço paramétrio. Para isto, informa-se os limites de cada um dos parâmetros nos argumentos lower e upper.

4.2 Função mle()

A função mle() do pacote stats4 também é específica para a estimação por máxima verossimilhança, oferecendo a possibilidade de estudar as verossimilhanças perfilhadas.

Deve ser informada a expressão do negativo do logaritmo da verossimilhança e os valores iniciais para os parâmetros em forma de lista nomeada. Esta função também utiliza internamente a função optim().

O resultado da função mle() é um objeto da classe mle, a partir do qual pode-se extrair a matriz de informação de Fisher, obtida através do hessiano numérico e as verossimilhanças perfilhadas para os parâmetros e com isso, construir intervalos de confiança assintóticos e visualizar as verossimilhanças perfilhadas. Esses resultados são extraídos utilizando-se os métodos implementados para o objeto da classe mle, tais como summary(), vcov(), profile(). Ao resultado da função profile() podem ser aplicados os métodos plot() e confint().

Para a amostra 2, simulada com p = 0.3 e θ = 0.1, tem-se:

A função vcov() extrai a matriz aproximada de variância e covariância das estimativas dos parâmetros, ou seja, o inverso da informação de Fisher observada aproximada pelo hessiano numérico.

As verossimilhanças perfilhadas são obtidas pela função profile. A partir de seu resultado, pode-se calcular o intervalo de confiança e plotar esses resultados. A verossimilhança perfilhada de um parâmetro é a análise do comportamento da função de verossimilhança, avaliada em valores próximos ao valor estimado desse parâmetro, mantendo os outros fixos nos respectivos valores estimados. O intervalo de confiança é obtido a partir das verossimilhanças perfilhadas para cada um dos parâmetros estimados.

O resultado da função profile também pode ser vizualizado em gráfico. O método plot() para essa classe produz um gráfico para cada parâmetro, tendo-se no eixo horizontal os valores do parâmetro e no eixo vertical, o valor da verossimilhança padronizado para z ~ N(0, 1). Isto facilita a vizualização da precisão da estimativa, sendo traçadas linhas de referência para intervalos de confiança. Fazendo z = 1.96 na linha vertical, obtém-se o valor correspondente ao limite do intervalo de confiança de 95% na linha horizontal. Na Figura 3, visualiza-se esse gráfico para os parâmetros estimados na amostra 2.

No gráfico da função padronizada de verossimilhança perfilhada de θ, observa-se que não é razoável construir um intervalo de confiança simétrico para θ a partir o erro padrão obtido numericamente. Pode-se obter os valores numéricos do intervalo de confiança a partir da verossimilhança perfilhada, fazendo:

O teste de hipótese pode ser feito a partir da razão de verossimilhanças ou da diferença de log-verossimilhanças. No caso particular da distribuição beta-binomial, para testar a hipótese de que θ é igual a zero, basta comparar a verossimilhança da hipótese com a verossimilhança segundo a distribuição binomial.

5 Aplicação em dados reais

A metodologia foi aplicada a dois conjuntos de dados. O primeiro contém dados de ppc em um talhão do município de Itajobi - MG e o segundo, dados de msc em um talhão localizado no município de Comendador Gomes - MG.

Na aplicação a dados de doenças de plantas, fez-se a contagem do número de plantas doentes por quadrat e aplicou-se a metodologia na contagem. O Rcitrus é um pacote desenvolvido para análise de dados de doenças em Citrus e contém funções para análise por quadrats, incluindo a metodologia da seção anterior.

5.1 Dados de Pinta Preta dos Citrus

O patógeno causador da ppc pode se reproduzir por ascósporos ou por pcnidósporos. Da primeira forma, o esporo está na fase sexual do patógeno e a disseminação ocorre pelo vento causando um padrão espacial aleatório. Da segunda forma, relacionado à fase assexual, a disseminação ocorre por escorrimento da água da chuva ou irrigação e a contaminação de uma árvore para outra se dá pelo contato com essa água.

A metodologia de estimação dos parâmetros da distribuição beta-binomial foi aplicada em um conjunto de dados de incidência de ppc. Estes dados são de uma fazenda localizada em Itajobi, Estado de Minas Gerais. Esses dados estão disponíveis no pacote Rcitrus.

O mapa da incidência pode ser visto na Figura 4. Pode-se observar que existe irregularidade no número de plantas nas linhas.

É preciso avaliar se é possível assumir que a incidência é constante em toda a área em estudo. A seguir utiliza-se a função quadrat.count() implementada no pacote Rcitrus para obter a contagem do número de plantas doentes por quadrats.

Observa-se quadrats incompletos, com menos de 100 plantas. Com esses dados explora-se as verossimilhanças perfilhadas fazendo a estimação que internamente utiliza a função mle(). Para isso utiliza-se a função fit.betabinom() disponível no pacote Rcitrus.

O intervalo de confiança para o parâmetro de agregação não abrange o zero, ou seja, o modelo binomial é insuficiente para explicar a variabilidade da incidência da doença. Com isso, pode-se concluir que o padrão espacial da incidência de Pinta Preta dos Citrus nesse talhão é agregado. O gráfico das verossimilhanças perfilhadas padronizadas para essa análise estão na Figura 5.

No pacote Rcitrus, está implementado a função disp.quadrats() para o teste de agregação baseado no índice de dispersão binomial e Poisson. Aplicando aos dados, fazendo a contagem por quadrats de tamanho 10 × 10, tem-se:

A variância observada é significativamente maior que a variância esperada segundo o modelo binomial e conclui-se que o padrão é agregado. Pode-se propor o modelo beta-binomial e testar se o parâmetro de agregação da distribuição beta-binomial θ é igual a zero. No argumento model, pode-se escolher a opção de utilizar a distribuição beta-binomial para testar a hipótese de aleatoriedade espacial. O p-valor retornado é o resultado de um teste de razão de verossimilhanças. O teste avalia a verossimilhança da distribuição beta-binomial com parâmetro θ igual ao estimado com os dados contra a verossimilhança da distribuição binomial.

Foi avaliado um número maior de quadrats que na análise do índice de dispersão binomial.

5.2 Dados de Morte Súbita dos Citrus

Segundo Bassanezi, Fernandes &MYammamoto (2003) a msc é uma nova doença dos Citrus que provoca rápido definhamento e morte de variedades de laranjas enchertadas, praticamente, apenas em limoeiro Cravo. O Estado de São Paulo responde por 80% da produção citrícola nacional e 85% de seus pomares apresentam laranjeiras doces enxertadas sobre limoeiro cravo. Portanto, há uma elevada vulnerabilidade da cultura à ocorrência de novas epidemias e particularmente para a msc.

O primeiro registro oficial da doença foi realizado em fevereiro de 2001 no município de Comendador Gomes, Estado de Minas Gerais. Em 2002 a msc atingiu o Estado de São Paulo e vários trabalhos tem sido conduzidos para a descoberta de informações a respeito. Esses trabalhos incluem a coleta de dados

O Citrus Sudden Death Virus - csdv, novo vírus da família Tymoviridae, tem sido associado à msc em pesquisas feitas pela empresa de biotecnologia Alellyx. Em abril de 2004 a mesma empresa anunciou a descoberta da presença do patógeno da MSC. O CSDV foi encontrado em três insetos, dois deles sendo capazes de transmitir o vírus para as plantas, (Alellyx 2004).

No talhão 303 da fazenda Vale Verde, localizada no município de Comendador Gomes - MG, foram feitas 25 avaliações. Esses dados estão disponíveis no pacote Rcitrus.

Na Figura 6, estão os mapas da evolução da msc ao long das primeiras 24 das 25 avaliações feitas. Não se apresenta o mapa da última avaliação por restrição de espaço, mas este pode ser gerado a partir dos dados disponíveis no Rcitrus.

Os códigos atribuídos à cada planta em cada avaliação são: 0 para plantas sadias, 1 para plantas com sintomas de msc, 2 para plantas em estado avançado de msc e 3 para plantas mortas. Considere o número de plantas em cada status, em cada avaliação:

O espaçamento nesse talhão é de 7.5 metros entre linhas e de 4 metros entre plantas na linha. Na análise, adota-se quadrats 2 × 4 considerando a incidência como sendo plantas identificadas pelos códigos 1, 2 e 3:

Observa-se que o padrão é aleatório nas avaliações 2 e 3 e agregado na primeira avaliação e nas avaliações 4 a 25. Esses resultados não mudam muito fazendo a análise com quadrats 3 × 5 ou 3 × 7 ou utilizando o índice de dispersão binomial.

A visualização da função desvio para θ em cada avaliação pode ser feita a partir das estimativas obtidas e da contagem por quadrats. Na Figura 7 estão plotadas as funções de deviance para θ estimado nas 25 avaliações utilizando quadrats 2 × 4.

6 Discussão e Resultados

A idéia de estimação por máxima verossimilhança foi apresentada em detalhes. Todos os passos envolvidos foram apresentados: problema da obtenção analítica dos estimadores, definição da função de densidade, definição da função de verossimilhança, maximização da função de verossimilhança via algoritmos de minimização numérica disponíveis em R, estimativa dos erros dos estimadores a partir do hessiano numérico, intervalos de confiança e teste de razão de verossimilhanças.

Foram apresentadas diferentes funções para maximização numérica, demonstrado o seu uso e, apontadas as particularidades de cada uma. Os procedimentos para estimação, obtenção dos erros padrão aproximados numericamente, obtenção da verossimilhança perfilhada e da função deviance foram exemplificados.

A metodologia foi exemplificada na estimação por máxima verossimilhança dos parâmetros da distribuição beta-binomial. Foi utilizada uma parametrização de forma que houvesse uma interpretação prática no contexto de análise de padrões espaciais de incidência de doenças. O teste da hipótese de aleatoriedade espacial deve ser realizado a partir do teste de razão de verossimilhanças, pois o estimador de θ não tem comportamento simétrico, quando a estimativa está próxima da borda do espaço paramétrico.

A análise do padrão espacial de doenças foi aplicada a dados de ppc e dados de msc. Para os dados de msc, foram analisados os dados de 25 avaliações feitas, em que se observou incidências desde 2.43% até 97.14%. A partir do teste de aleatoriedade espacial utilizando o parâmetro de agregação da distribuição beta-binomial, verificou-se que o padrão espacial da ppc é agregado no talhão analisado. Na análise dos dados de msc, verificou-se que as evidências de padrão agregado são menores nas avaliações com incidência baixa da msc.

Agradecimentos

Ao Laboratório de Estatística e Geoinformação (leg), pelo ambiente computacional. Ao Fundo de Defesa da Citricultura (fundecitrus) pelo financiamento do projeto em convênio com o Departamento de Estatística da Universidade Federal do Paraná, em cujo escopo estas atividades foram desenvolvidas. À Ana Beatriz Tozzo Martins pela revisão e comentários.

Referências