MODELOS HIPSOMÉTRICOS GENERALIZADOS MISTOS NA PREDIÇÃO DA ALTURA DE Eucalyptus sp.

GENERALIZED MIXED HEIGHT-DIAMETER MODELS FOR Eucalyptus sp. HEIGHT PREDICTION

Resumos

Conduziu-se, este trabalho, com o objetivo de avaliar a teoria dos modelos generalizados mistos e a modelagem da heterogeneidade de variância no ajuste de dois modelos não lineares para a relação altura-diâmetro em Eucalyptus sp. Foram utilizados dados de árvores-amostras de diferentes materiais genéticos (clone e semente) e idades. De acordo com as análises, o modelo assintótico obteve os melhores resultados em relação ao modelo sigmoidal. A precisão dos modelos avaliados aumentou com a utilização da teoria dos efeitos mistos com adição de covariáveis e a modelagem da heterogeneidade da variância.

efeito misto; modelos estocásticos; heterocedasticidade; covariáveis

The purpose of this study was to evaluate the theory of the generalized mixed models and the modeling of the heteroscedasticity of variance in the fit of two nonlinear models for the height-diameter relationship in Eucalyptus sp. Data from sample trees of different genotypes (clone and seed) and ages were used. According to the analysis, the asymptotic model achieved the best results compared to the sigmoid model. The accuracy of the models increased using the theory of mixed effects with the addition of covariates and modeling of the heteroscedasticity of variance.

mixed effect; probability models; heteroscedasticity; covariates

INTRODUÇÃO

A medição de alturas em povoamentos florestais é uma atividade onerosa, uma vez que, em comparação com a obtenção do diâmetro, sua obtenção não é fácil. Com isso, nos inventários florestais é utilizada uma ferramenta que ajuda na economia de tempo e capital, ou seja, os modelos baseados em relações hipsométricas.

Alguns autores estudaram modelos para representarem a relação hipsométrica em povoamentos florestais, dentre eles: Barros et al. (2002), Bartoszeck et al. (2003), Batista et al. (2001), Couto e Bastos (1987), Soares et al. (2004) e Tomé et al. (2007). Só que, na maioria desses estudos, não foram utilizadas a teoria dos modelos generalizados de efeito misto, modelagem da heterogeneidade da variância e adição de covariáveis e, também, muitos destes utilizaram apenas modelos lineares. No Brasil, o uso de modelos mistos ou modelos generalizados mistos podem ser encontrados nos trabalhos de Calegario et al. (2005b, 2005c) sobre crescimento e/ou produção e Horle et al. (2010), estudando modelos de afilamento. Entre os trabalhos que utilizaram modelos mistos na relação hipsométrica podem ser citados: Jayaraman e Lappi (2001), em Tectona grandis L.f., na Índia; Sharma e Parton (2007), em oito espécies da floresta Boreal, no Canadá; e Adame et al. (2008), em Quercus pyrenaica Willd, na Espanha. A modelagem da heterogeneidade da variância foi pouco estudada em modelos biométricos florestais, podendo citar o trabalho de Calegario et al. (2005a). O uso de covariáveis pode ser encontrado em Guimarães et al. (2009) e Mendonça et al. (2011).

Visto que, a homogeneidade da variância é uma pressuposição da análise de regressão e que existe a necessidade de obter modelos de relação hipsométrica mais precisos, é necessário mais estudos nessa área. Nesse contexto, o objetivo do trabalho foi avaliar a teoria dos modelos generalizados mistos e a modelagem da heterogeneidade de variância no ajuste de modelos não lineares para a relação altura-diâmetro em Eucalyptus sp.

MATERIAL E MÉTODOS

Caracterização e análise dos dados

Os dados utilizados, coletados em plantios de Eucalyptus sp. da empresa International Paper do Brasil, localizada em Mogi das Cruzes, São Paulo. O clima da região é classificado por Köppen como Cwa. A temperatura média é inferior a 18 ºC e a precipitação média de 1300 mm. As árvores-amostra são provenientes de 12 (doze) materiais genéticos diferentes, sendo dez de plantios clonais (1, 3, 4, 5, 6, 7, 8, 9, 11 e 12) e dois de plantios por semente (2 e 10), com idade variando entre 2 a 7,45 anos. As árvores-amostra tiveram sua altura total (H) medida por hipsômetro e a circunferência a 1,30m do solo (CAP) por meio de uma fita métrica e foram divididas em 7 classes diamétricas (Tabela 1).

Tabela 1:

Distribuição das classes diamétricas das árvores-amostra utilizadas no ajuste dos modelos.

Distribution of the diameter classes of the sample-trees used for fitting the models.

Modelos Avaliados

Foram avaliados os seguintes modelos não lineares:

a) Regressão fixa

a.1) Sigmoidal (1), em que: H_i = altura total da árvore (m); DAP_i = diâmetro a 1,30 m do solo (cm); e = base do logaritmo neperiano; = parâmetros do modelo; = erro estocástico. O parâmetro representa a assíntota horizontal superior (AHS), isto é, o valor máximo da resposta (H _i) quando o valor de DAP _i . O parâmetro representa o ponto de inflexão da curva, ou seja, o DAP _i onde a H _i atingir metade de . O parâmetro (escala) representa a diferença entre os DAP _s onde H _i atingir aproximadamente 73% de e o DAP _i no ponto de inflexão.

a.2) Assintótico (2), onde o parâmetro representa a assíntota horizontal superior (AHS). O parâmetro representa o valor da resposta (Hi) quando o valor de DAPi=0 e o parâmetro (escala) representa o logaritmo da taxa constante. Maiores detalhes sobre esses modelos podem ser encontrados em Pinheiro e Bates (2000).

b) Regressão mista com inclusão de covariáveis (3), em que m é o número de grupos (classes de DAP); n_i é número de observações do i-ésimo grupo; f é uma função geral, real e diferenciável de um grupo específico de vetores do parâmetro e de uma covariável , e é, normalmente, distribuído dentro do grupo.

A função f deve ser não linear em no mínimo um componente de um grupo específico de vetor de parâmetros , o qual é representado pela equação (4), onde β é um vetor p dimensional de efeitos fixos em que são inseridos os efeitos de covariáveis (idade e material genético) e b_i é um vetor q dimensional de efeitos aleatórios (classe de DAP) associados com o i-ésimo grupo com matriz de variância e covariância . Maiores detalhes sobre modelos mistos podem ser encontrados em Pinheiro e Bates (2000).

Modelagem da variância

Para a modelagem da variância foi utilizada uma função de variância. Funções de variância são utilizadas para modelar a estrutura de variância de erros dentro dos grupos, usando covariáveis com um parâmetro delta (δ) exponencial, estimado por processos iterativos. Foi escolhida a classe de função de variância denominada varPower. O modelo de variância (5) proposto por essa classe é expresso por Pinheiro e Bates (2000).

Correspondente à função de variância (6), onde é uma potência do valor absoluto da variância da covariável. O parâmetro δ é irrestrito, ou seja, pode assumir qualquer valor real e, assim, pode modelar casos onde a variância aumenta ou diminui com o valor absoluto da variância da covariável.

Avaliação dos modelos

Os modelos utilizados tiveram sua precisão testada por meio das seguintes estatísticas:

a) Erro padrão (7), onde S _yx = erro padrão da estimativa (m); S _yx (%) = erro padrão relativo (%); altura total (H _i) estimada pelo modelo (m); Y = altura total (H _i) observada (m); n = número de observações; p = número de parâmetros.

b) Critério de Informação de Akaike (8) (AIC) (SAKAMOTO et al., 1986), onde p = número de parâmetros do modelo; mv = valor da máxima verossimilhança. Quanto menor o valor de AIC, melhor o modelo avaliado.

c) Critério de Informação Bayesiano (9) (BIC) (SCHWARZ, 1978). Quanto menor o valor de BIC, melhor o modelo avaliado.

Os critérios AIC e BIC são de grande relevância, quando se trata da análise de regressão, pois penalizam modelos com grande número de parâmetros, escolhendo modelos mais parcimoniosos.

Além dessas estatísticas, foram feitas, também, a análise gráfica dos valores residuais (%) dados por (10).

Para análise dos dados, foi utilizado o pacote nlme do software estatístico R versão 2.10.1. (R DEVELOPMENT CORE TEAM, 2011).

RESULTADOS E DISCUSSÃO

Análise do ajuste dos modelos avaliados sem correção da heterocedasticidade

Analisando os resultados da Tabela 2 das estatísticas para os modelos, relacionando a altura total (H _i) como variável resposta e o diâmetro a 1,3 m (DAP _i) como variável regressora (modelos simples), observa-se que os modelos analisados tiveram significância em seus parâmetros. Se for considerado somente esse critério, os dois modelos podem ser utilizados para estimativa da altura total das árvores de Eucalyptus sp. Mas, considerando o erro padrão relativo [Syx(%)], o modelo assintótico deve ser escolhido para o objetivo em questão.

Tabela 2:

Estatísticas do ajuste dos modelos na sua forma original.

Statistics for the fitted models in the original form.

• Em que: simples = modelo ajustado de altura total em função do DAP.

Após a análise dos modelos simples, os mesmos foram analisados considerando a classe de DAP como efeito aleatório. Na Tabela 3, estão as estatísticas de ajuste, para os modelos analisados, considerando a teoria dos modelos mistos. Na Tabela 4, estão apresentados os efeitos da classe de DAP sobre os parâmetros do modelo sigmoidal e o modelo assintótico da Tabela 3.

Tabela 3:

Estatísticas do ajuste dos modelos usando a teoria dos modelos mistos.

Statistics of the fitted models using the theory of mixed models.

• Em que: misto = modelo ajustado usando a teoria dos efeitos mistos.

Tabela 4:

Efeito das classes de DAP sobre os parâmetros dos modelos da Tabela 3.

DBH classes effects on the parameters of the models from Table 3.

Pela Tabela 3, nota-se que os parâmetros do modelo sigmoidal foram significativos e o parâmetro do modelo assintótico foi não significativo, a 5% de probabilidade. Isso significa que o uso do modelo assintótico para esse conjunto de dados não garante estimativas de altura próximas de zero (0), quando o DAP das árvores é igual a zero (0), como preconizado pela interpretação biológica do modelo. Com isso, o modelo sigmoidal é o mais adequado nessa situação, fato comprovado pelo menor valor de Syx (%) em relação ao modelo assintótico. Observa-se, também, um ganho em precisão (avaliado como a redução do Syx %), quando se aplica à classe diamétrica como efeito aleatório em relação ao modelo simples (Tabela 2) para o modelo sigmoidal (4,14%) e modelo assintótico (0,98%). Essa melhoria na precisão dos modelos com a adição de um efeito aleatório, além dos efeitos fixos, é similar a pesquisas realizadas por outros autores (ADAME et al., 2008; CALEGARIO et al., 2005b). Segundo Adame et al. (2008), tanto o viés quanto a raiz quadrada da média dos resíduos é menor no modelo com parâmetro aleatório, e diminuem com o número de árvores de amostra utilizados para fins de calibração do modelo avaliado para estimar a relação hipsométrica de Quercus pyrenaica. Já, Dorado et al. (2006), afirma que a calibração do modelo hipsométrico avaliado permitiu resultados precisos a serem obtidos com uma amostra pequena, fazendo com que essa abordagem altamente eficaz e útil.

Na Tabela 5, apresentam-se as estatísticas de ajuste, para os modelos analisados, considerando a teoria dos modelos mistos com inclusão das covariáveis material genético e idade (modelo completo). Na Tabela 6, apresentam-se as estimativas dos parâmetros, para cada classe de DAP, para o modelo sigmoidal e o modelo assintótico completo (Tabela 5).

Tabela 5:

Estatísticas de ajuste, para os modelos analisados, considerando a teoria dos modelos mistos e a inclusão dos efeitos das covariáveis material genético e idade.

Statistics of fit for the tested models, considering the theory of mixed models and the inclusion of the covariates effects of the genetic material and age.

• Em que: completo = modelo ajustado usando a teoria dos efeitos mistos e adição dos efeitos das covariáveis idade e material genético nos parâmetros.

Tabela 6:

Efeito das classes de DAP, idade e material genético (MG)sobre os parâmetros .dos modelos da Tabela 5.

Effect of DBH classes, age and genetic material (MG) on the parameters of the models from Table 5.

• Em que: MG = material genético.

Pela Tabela 5, verifica-se que o modelo assintótico tem uma pequena vantagem, considerando o valor de Syx (%), em relação ao modelo sigmoidal. Nota-se, também, que os parâmetros dos modelos sigmoidal completo e assintótico completo não tiveram comportamentos semelhantes aos modelos sigmoidal misto e assintótico misto. Isso é verificado, principalmente, para o parâmetro . Isso ocorre, em decorrência, principalmente, do efeito da idade e material genético influenciar, diretamente na máxima resposta das árvores em altura total. Os efeitos dessas covariáveis (Tabela 6) são adicionados ao valor médio das estimativas dos parâmetros sem o efeito das mesmas (Tabela 5), para obter o valor das estimativas dos parâmetros para cada combinação idade, material genético e classe de diâmetro. Nota-se que, para o material genético 10 (Tabela 6), obtiveram-se estimativas dos efeitos desse material genético diferentemente do esperado para o modelo sigmoidal. Esperava-se que os valores de , e encontrados fossem próximos à média de todos os outros materiais genéticos. Isso pode ter ocorrido, pela quantidade de árvores-amostra terem sido reduzidas em relação aos demais materiais genéticos, o que gerou estimativas do efeito desse material genético para o modelo sigmoidal sem realismo biológico e, estatisticamente, não significativas. Esse fato significa que as estimativas das alturas das árvores provenientes desse material genético são de baixa precisão. Nota-se um ganho em precisão de 33,06% para o modelo sigmoidal e 36,50% para o modelo assintótico.

Na Figura 1, apresenta-se a distribuição gráfica dos resíduos para os modelos avaliados na estimativa da relação hipsométrica de Eucalyptus sp.

Figura 1:

Distribuição dos resíduos, em percentagem, em função da altura total estimada.

Residual distribution, in percent, in function of the estimated total height.

Analisando os resultados da Figura 1, nota-se uma melhoria na estimativa quando se aplica a teoria dos modelos generalizados mistos com adição de covariáveis. Isso pode ser explicado, quando se compara a distribuição das alturas estimadas e a distribuição das alturas reais. Considerando que a distribuição real da altura total varia de 5,7 metros a 35,2 metros, nota-se que os modelos sigmoidal e assintótico completos conseguem estimar alturas mais condizentes às alturas reais. Pela Figura 1, observa-se, também, que o modelo assintótico foi mais preciso, em relação ao modelo sigmoidal, independente do tipo de ajuste feito.

Na Tabela 7, apresentam-se as estatísticas de precisão para os modelos ajustados. Os resultados da Tabela 7 mostram que o modelo assintótico completo foi mais preciso, seguido pelo modelo sigmoidal completo, assintótico misto, sigmoidal misto, assintótico simples e sigmoidal simples, respectivamente. Os resultados apresentados na Tabela 7 estão de acordo com as estatísticas das Tabelas 2 e 5 e da análise gráfica dos resíduos (Figura 1). O ganho de precisão com adição de covariáveis relacionadas ao povoamento, também foi obtido por Guimarães et al. (2009) e Mendonça et al. (2011). Soares et al. (2004) obtiveram resultados precisos, quando utilizaram modelos que consideravam o efeito de variáveis do povoamento nos modelos hipsométricos genéricos.

Tabela 7:

Critério de Informação de Akaike (AIC), Critério de Informação Bayesiano (BIC) dos modelos ajustados.

Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC) of the adjusted models.

• Em que: simples = modelo ajustado de altura total em função do DAP; misto = modelo ajustado usando a teoria dos efeitos mistos; e completo = modelo ajustado usando a teoria dos efeitos mistos e adição dos efeitos das covariáveis idade e material genético nos parâmetros.

Comparação dos modelos com a modelagem da heterogeneidade da variância

Após as análises, foram escolhidos os modelos sigmoidal completo e assintótico completo e procedeu-se â correção da heterocedasticidade. Na Tabela 8, apresenta-se o erro padrão da estimativa e o erro padrão relativo aos modelos ajustados com a correção da heterocedasticidade.

Tabela 8:

Erro padrão da estimativa (Syx) e erro padrão relativo [Syx (%)] dos modelos ajustados com a correção da heterocesdaticidade.

Standard error of estimate (Syx) and relative standard error [Syx (%)] of the adjusted models with the correction of the heteroscedasticity.

Avaliando a Tabela 8, nota-se um ganho de 52,35% para o modelo sigmoidal e 53,89% para o modelo assintótico com a modelagem da heterocedasticidade em relação ao modelo homocedástico (Tabela 5). Pela Tabela 8, verifica-se, também, que o modelo assintótico foi o mais preciso em relação ao modelo sigmoidal, como nas análises sem a correção da heterocedasticidade (Tabelas 2 e 5 e Figura 1).

Na Figura 2, apresenta-se a distribuição dos resíduos para os modelos avaliados na estimativa da relação hipsométrica de árvores Eucalyptus sp., após a correção da heterocedasticidade.

Figura 2:

Distribuição dos resíduos, em percentagem, em função da altura total estimada para os modelos sigmoidal e assintótico heterocedásticos.

Residual distribution, in percentage, in function of the estimated total height for the sigmoidal and asymptotic heterocedastic models.

Pela Figura 2, nota-se uma distribuição residual semelhante. Esses modelos tiveram a distribuição de alturas estimadas muito próximas à distribuição real da altura total (5,7 metros a 35,2 metros). A análise da Figura 2 está de acordo com a análise da Tabela 8.

De acordo com esses resultados, nota-se a grande importância de se estudar a aplicação de covariáveis, a aplicação dos modelos mistos e a correção da heterogeneidade de variância na modelagem da relação hipsométrica de árvores de Eucalyptus sp. Percebe-se, também, que esses estudos podem englobar a adição de outras covariáveis que expliquem o crescimento em altura de árvores de espécies florestais, como, por exemplo, variáveis ambientais. Isso melhora a precisão e, também, diminui o número de equações utilizadas quando se tem uma base de dados com diferentes materiais genéticos, idades, densidades de plantio e capacidade produtiva. Outro fato interessante a se ressaltar é que o uso das técnicas utilizadas, neste trabalho, acarreta na diminuição do número de árvores-amostras necessária para a calibração do modelo com precisão adequada para ser utilizado, fato esse já citado por Dorado et al. (2006). Com isso, diminui-se o tempo e, consequentemente, os custos dos inventários florestais. O uso de modelos com explicação biológica, como os utilizados neste trabalho, é de grande valia, pois explica a relação altura-diâmetro dos dados reais e possuem interpretação biológica dos parâmetros. Esse fato é importante, uma vez que o uso de modelos não lineares se limita pela necessidade do conhecimento de estimativas iniciais dos parâmetros para começar o processo iterativo.

CONCLUSÕES

De acordo com os resultados obtidos, para as condições em que foi desenvolvido este estudo, conclui-se que: a) O modelo assintótico foi o mais preciso na estimativa da relação hipsométrica de Eucalyptus sp.; b) A utilização da teoria dos modelos de regressão mista e a modelagem da heterogeneidade de variância possuem grande potencial de utilização em modelagem florestal, uma vez que essas técnicas, juntamente, forneceram um grande aumento na precisão dos modelos analisados e seu uso, muitas vezes, restringe-se na medicina e no melhoramento genético.

AGRADECIMENTOS

À Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), pela concessão da bolsa de estudos.

• ADAME, P.; DEL-RÍO, M.; CAÑELLAS, I. A mixed nonlinear height-diameter model for pyrenean oak (Quercus pyrenaica Willd.). Forest Ecology and Management, Amsterdam, v. 256, n. 1/2, p. 88-98, 2008.
• BARROS, D. A.; MACHADO, S. A.; ACERBI JUNIOR, F. W.; SCOLFORO, J. R. S. Comportamento de modelos hipsométricos tradicionais e genéricos para plantação de Pinus oocarpa em diferentes tratamentos. Boletim de Pesquisa Florestal, Colombo, n. 45, p. 3-28, 2002.
• BARTOSZECK, A. C. P. S.; MACHADO, S. A.; FIGUEIREDO FILHO, A.; OLIVEIRA, E. B. Modelagem da relação hipsométrica para bracatingais da região metropolitana de Curitiba, PR. Floresta, Curitiba, v. 32, n. 2, p. 189-204, 2003.
• BATISTA, J. L. F.; COUTO, H. T. Z.; MARQUESINI, M. Desempenho de modelos de relação hipsométrica: estudo em três tipos de floresta. Scientia Forestalis, Piracicaba, n. 60, p. 149-163, 2001.
• CALEGARIO, N.; CALEGARIO, C. L. L.; MAESTRI, R.; DANIELS, R. F. Melhoria da qualidade de ajuste de modelos biométricos florestais pelo emprego da teoria dos modelos não lineares generalizados. Scientia Forestalis, Piracicaba, n. 69, p. 38-50, 2005a.
• CALEGARIO, N.; DANIELS, R. F.; MAESTRI, R.; NEIVA, R. Modeling dominant height growth based on nonlinear mixed-effects model: a clonal Eucalyptus plantation case study. Forest Ecology and Management, Amsterdam, v. 204, n. 1, p. 11-20, 2005b.
• CALEGARIO, N.; MAESTRI, R.; LEAL, C. L.; DANIELS, R. F. Estimativa de crescimento de povoamentos de Eucalyptus baseada na teoria dos modelos não-lineares em multinível de efeito misto. Ciência Florestal, Santa Maria, v. 15, n. 3, p. 285-292, 2005c.
• COUTO, H. T. Z.; BASTOS, N. L. M. Modelos de equações de volume e relações hipsométricas para plantações de Eucalyptus no estado de São Paulo. IPEF, Piracicaba, n. 37, p. 33-44, 1987.
• DORADO, F. C.; DIÉGUEZ-ARANDA, U.; ANTA, M. B.; RODRIGUEZ, M. S.; VON GADOW, K. A generalized height-diameter model including random components for radiata pine plantations in northwestern Spain. Forest Ecology and Management, Amsterdam, v. 229, p. 202-213, 2006.
• GUIMARÃES, M. A. M.; CALEGARIO, N.; CARVALHO, L. M. T.; TRUGILHO, P. F. Height-diameter models in forestry with inclusion of covariates. Cerne, Lavras, v. 15, n. 3, p. 313-321, 2009.
• HORLE, D. C.; MENDONÇA, A. R.; CARVALHO, S. P. C.; CALEGARIO, N. Modelagem não linear do perfil longitudinal de fustes de Pinus oocarpa. Cerne, Lavras, v. 16, n. 2, p. 177-184, 2010.
• JAYARAMAN, K.; LAPPI, J. Estimation of height-diameter curves through multilevel models with special reference to even-aged teak stands. Forest Ecology and Management, Amsterdam, v. 142, n. 1/3, p. 155-162, 2001.
• MENDONÇA, A. R.; CALEGARIO, N.; SILVA, G. F.; BORGES, L. A. C.; CARVALHO, S. P. C. Modelos hipsométricos e de crescimento em altura das árvores dominantes e codominantes para Pinus caribaea var. hondurensis. Scientia Forestalis, Piracicaba, v. 39, n. 90, p. 151-160, 2011.
• PINHEIRO, J. C.; BATES, D. M. Mixed-effects models in S and S-PLUS. New York: Springer-Verlag, 2000. 528 p.
• R DEVELOPMENT CORE TEAM. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing, 2011. Disponível em: <http://www.R-project.org/>. Acesso em: 10 mar. 2012.
  » http://www.R-project.org/
• SAKAMOTO, Y.; ISHIGURO, M.; KITAGAWA, G. Akaike information criterion statistics. Dordrecht: Kluwer Academic, 1986.
• SCHWARZ, G. Estimating the dimension of a model. Annals of Statistics, Philadelphia, v. 6, p. 461-464, 1978.
• SHARMA, M.; PARTON, J. Height-diameter equations for boreal tree species in Ontario using a mixed-effects modeling approach. Forest Ecology and Management, Amsterdam, v. 249, n. 3, p. 187-198, 2007.
• SOARES, T. S.; SCOLFORO, J. R. S.; FERREIRA, S. O.; MELLO, J. M. Uso de diferentes alternativas para viabilizar a relação hipsométrica no povoamento florestal. Revista Árvore, Viçosa, v. 28, n. 6, p. 845-854, nov./dez. 2004.
• TOMÉ, M.; RIBEIRO, F.; FAIAS, S. Relação hipsométrica geral para Eucalyptus globulus Labill. em Portugal. Silva Lusitânica, Lisboa, v. 15, n. 1, p. 41-55, 2007.

Datas de Publicação

• Publicação nesta coleção
  Jan-Mar 2015

Histórico

• Recebido
  18 Abr 2011
• Aceito
  31 Ago 2014