diff --git a/docs/01-tcc.Rnw b/docs/01-tcc.Rnw index a50f3e6fb82aad9dc4cd91753dc9e7c81cbd7f6c..c6070c01b9baf72dcf81e8ecaa1c09311f3c3b92 100644 --- a/docs/01-tcc.Rnw +++ b/docs/01-tcc.Rnw @@ -1,24 +1,3 @@ -%% abtex2-modelo-trabalho-academico.tex, v<VERSION> laurocesar -%% Copyright 2012-<COPYRIGHT_YEAR> by abnTeX2 group at http://www.abntex.net.br/ -%% -%% This work may be distributed and/or modified under the -%% conditions of the LaTeX Project Public License, either version 1.3 -%% of this license or (at your option) any later version. -%% The latest version of this license is in -%% http://www.latex-project.org/lppl.txt -%% and version 1.3 or later is part of all distributions of LaTeX -%% version 2005/12/01 or later. -%% -%% This work has the LPPL maintenance status `maintained'. -%% -%% The Current Maintainer of this work is the abnTeX2 team, led -%% by Lauro César Araujo. Further information are available on -%% http://www.abntex.net.br/ -%% -%% This work consists of the files abntex2-modelo-trabalho-academico.tex, -%% abntex2-modelo-include-comandos and abntex2-modelo-references.bib -%% - % ------------------------------------------------------------------------ % ------------------------------------------------------------------------ % abnTeX2: Modelo de Trabalho Academico (tese de doutorado, dissertacao de @@ -45,6 +24,21 @@ svgnames ]{abntex2} +% --- +% Novo list of (listings) para QUADROS +% --- + +\newcommand{\quadroname}{Quadro} +\newcommand{\listofquadrosname}{Lista de quadros} + +\newfloat[chapter]{quadro}{loq}{\quadroname} +\newlistof{listofquadros}{loq}{\listofquadrosname} +\newlistentry{quadro}{loq}{0} + +% configurações para atender à s regras da ABNT +\counterwithout{quadro}{chapter} +\renewcommand{\cftquadroname}{\quadroname\space} +\renewcommand*{\cftquadroaftersnum}{\hfill--\hfill} % --- % PACOTES @@ -97,6 +91,7 @@ \usepackage{tikz} \usepackage{pdflscape} % para ambiente landscape \usepackage{pgfgantt} % cronograma estilo gráfico de gantt +\usepackage{multicol} \usetikzlibrary{backgrounds} % --- @@ -142,7 +137,7 @@ % --- % Informações de dados para CAPA e FOLHA DE ROSTO % --- -\titulo{Extensões e Aplicações Modelo de Regressão +\titulo{Extensões e Aplicações do Modelo de Regressão Conway-Maxwell-Poisson para Modelagem de Dados de Contagem} \vspace{2cm} \autor{Eduardo Elias Ribeiro Junior} @@ -227,7 +222,7 @@ <<setup, include=FALSE, cache=FALSE>>= source("_setup.R") -library(tccPackage) +library(cmpreg) @ @@ -314,29 +309,29 @@ library(tccPackage) representam o número de ocorrências de um evento em um domÃnio discreto ou contÃnuo. Para análise estatÃstica dessas variáveis, o modelo de Poisson é amplamente utilizado. Porém, não são raras as - situações de sub ou superdispersão, que inviabilizam o emprego deste + situações de sub ou superdispersão, que inviabilizam o emprego desse modelo. Uma alternativa paramétrica é o modelo COM-Poisson que, com a adição de um parâmetro, contempla diferentes nÃveis de dispersão. - Outras caracterÃsticas bastantes frequentes em dados de contagem são - frequência excessiva de valores zeros e estrutura de correlação entre - observações, muitas vezes induzida pelo processo de casualização ou - amostragem. Nesses casos os modelos adotados devem ser - adaptados. Neste trabalho são exploradas as caracterÃsticas da - distribuição COM-Poisson e apresentados os modelos de regressão - COM-Poisson de efeitos fixos, com modelagem para excesso de zeros e - incluindo efeitos aleatórios. O emprego dos modelos COM-Poisson e suas - extensões é ilustrado com aplicações onde seus resultados são - comparados com as abordagens Poisson, Quasi-Poisson e Binomial - Negativa (para casos de superdispersão) via nÃveis descritivos de - testes de razão de verossimilhanças, critério de informação de Akaike - e predições pontuais e intervalares. O ajuste dos modelos é feito via - maximização da verossimilhança. Os resultados mostram que o modelo - Poisson é de fato restritivo, com ajustes inadequados na maioria das - aplicações. O modelo COM-Poisson, por sua vez, mostrou-se bastante - flexÃvel com resultados similares aos obtidos via abordagem - semi-paramétrica Quasi-Poisson. As extensões propostas para o modelo - COM-Poisson apresentaram resultados satisfatórios, sendo equivalentes - as abordagens já consolidadas na literatura. + Outras caracterÃsticas frequentes em dados de contagem são excesso de + contagens nulas e estrutura de correlação entre observações, muitas + vezes induzida pelo processo de casualização ou amostragem. Nesses + casos os modelos adotados devem ser adaptados. Neste trabalho são + exploradas as caracterÃsticas da distribuição COM-Poisson e + apresentados os modelos de regressão COM-Poisson de efeitos fixos, com + modelagem para excesso de zeros e incluindo efeitos aleatórios. O + emprego dos modelos COM-Poisson e suas extensões é ilustrado com + aplicações e seus resultados são comparados com as abordagens Poisson, + Quasi-Poisson e Binomial Negativa (para casos de superdispersão) via + nÃveis descritivos de testes de razão de verossimilhanças, critério de + informação de Akaike e predições pontuais e intervalares. O ajuste dos + modelos é feito via maximização da verossimilhança. Os resultados + mostram que o modelo Poisson é de fato restritivo, com ajustes + inadequados na maioria das aplicações. O modelo COM-Poisson, por sua + vez, mostrou-se bastante flexÃvel apresentando resultados similares + aos obtidos via abordagem semi-paramétrica Quasi-Poisson. As extensões + propostas para o modelo COM-Poisson apresentaram resultados + satisfatórios, sendo equivalentes à s abordagens já consolidadas na + literatura. \textbf{Palavras-chave}: COM-Poisson; dados de contagem; subdispersão; superdispersão; excesso @@ -359,6 +354,14 @@ library(tccPackage) \cleardoublepage % --- +% --- +% inserir lista de quadros +% --- +\pdfbookmark[0]{\listofquadrosname}{loq} +\listofquadros* +\cleardoublepage +% --- + % --- % inserir lista de abreviaturas e siglas % --- diff --git a/docs/01-tcc.pdf b/docs/01-tcc.pdf index 9ca8f38c8e826a459c9a1b8d1cab359cb059cc19..3c40772a8dbb3374674a78d31e2b50e756f2a21b 100644 Binary files a/docs/01-tcc.pdf and b/docs/01-tcc.pdf differ diff --git a/docs/cap01_introducao.Rnw b/docs/cap01_introducao.Rnw index 84c60fa620a657e47cfb802942494a6582150bf4..838213ce9014a8d6ee0dde2b64e057a9284c6f05 100644 --- a/docs/cap01_introducao.Rnw +++ b/docs/cap01_introducao.Rnw @@ -6,20 +6,20 @@ Em diversas áreas do conhecimento é comum o interesse em i) compreender o relacionamento entre variáveis de interesse e caracterÃsticas de uma amostra e ii) realizar predições por meio de modelos estatÃsticos ajustados por dados de uma amostra. A teoria de modelos de regressão -sustentam muitas das pesquisas na área de EstatÃstica aplicada. +sustenta muitas das pesquisas na área de EstatÃstica aplicada. Os modelos de regressão, na sua forma univariada e usual, consistem no estabelecimento de uma equação matemática que relaciona a média de uma variável aleatória de interesse (variável resposta) com as demais -variáveis observadas (covariáveis). Nesta metodologia considera-se uma -distribuição de probabilidades para a variável resposta condicionada as -covariáveis cuja média está associada a uma preditor que acomoda os -efeitos das covariáveis. +variáveis observadas (covariáveis). Nessa metodologia considera-se uma +distribuição de probabilidades para a variável resposta condicionada à s +covariáveis cuja média está associada a um preditor que acomoda os +efeitos dessas covariáveis. Pode-se destacar o modelo linear normal como o de uso predominante dentre os disponÃveis para análises estatÃsticas aplicadas. Esse modelo -estabelece que a variável resposta condicional as covariáveis tem -distribuição Normal de média descrita por um preditor linear das +estabelece que a variável resposta, condicional à s covariáveis, tem +distribuição Normal, de média descrita por um preditor linear das covariáveis. Todavia, não são raras as situações em que a variável resposta é uma contagem, assumindo valores inteiros não negativos. Variáveis aleatórias de contagem, de forma geral, representam @@ -28,7 +28,7 @@ ser contÃnuo, como um intervalo de tempo ou espaço, ou discreto, como indivÃduos ou grupos. A análise de dados de contagem pelo modelo linear normal produz -estimativas que contêm erros padrões inconsistentes e podem produzir +estimativas que contêm erros padrões inconsistentes e pode produzir predições negativas para o número de eventos \cite{King1989}. Uma alternativa adotada durante muitos anos, e ainda aplicada, é encontrar alguma forma de transformação da variável resposta a fim de atender aos @@ -40,57 +40,59 @@ a relação média e variância, caracterÃstica de dados de contagem e iv) o uso da transformação logarÃtmica é problemática quando há contagens nulas. -Diante do problema diferentes abordagens foram propostas, contudo -destaca-se o trabalho apresentado por \citeonline{Nelder1972} que -introduz a teoria dos modelos lineares generalizados (MLG's). Esta nova -classe de modelos flexibilizou a distribuição condicional permitindo que -outras distribuições pertencentes à famÃlia exponencial fossem -consideradas para a distribuição da variável resposta. Tal famÃlia -contempla as distribuições Poisson, Binomial, Gama entre outras bem -conhecidas na literatura, além da própria distribuição Normal. - -Com os MLG's a modelagem de dados passou a ser mais fiel a natureza da +Diante dos problemas relatados na aplicação de modelos normais para +análise de dados de contagem, diferentes abordagens foram +propostas. Destaca-se o trabalho apresentado por \citeonline{Nelder1972} +que introduz a teoria dos modelos lineares generalizados (MLG's). Essa +nova classe de modelos flexibilizou a distribuição condicional +permitindo que outras distribuições pertencentes à famÃlia exponencial +fossem consideradas para a distribuição da variável resposta. Tal +famÃlia contempla as distribuições Poisson, Binomial, Gama entre outras +bem conhecidas na literatura, além da própria distribuição Normal. + +Com os MLG's a modelagem de dados passou a ser mais fiel à natureza da variável resposta, principalmente no que diz respeito ao seu suporte. Nesse contexto, a análise de variáveis aleatórias de contagem, que têm suporte nos conjunto dos números naturais, foi enriquecida expressivamente. Para análise estatÃstica dessas variáveis, o modelo probabilÃstico de -Poisson, já consolidado na literatura é amplamente utilizado. Este +Poisson, já consolidado na literatura, é amplamente utilizado. Esse modelo possui apenas um parâmetro, denotado por $\lambda$, que representa a média e também a variância, o que implica em uma relação identidade ($\lambda = E(Y) = V(Y)$). Essa propriedade, chamada de equidispersão, é uma particularidade do modelo Poisson que pode não ser -adequada a diversas situações. Quando aplicado sob negligência desta +adequada a diversas situações. Quando aplicado sob negligência dessa suposição, o modelo Poisson apresenta erros padrões inconsistentes para -as estimativas dos parâmetros e por consequência, para toda função +as estimativas dos parâmetros e, por consequência, para toda função desses parâmetros \cite{Winkelmann1995, Winkelmann1994}. O caso de superdispersão, quando a variância é maior que a média, é o mais comum e existe uma variedade de métodos para análise de dados -assim. A superdispersão pode ocorrer pela ausência de covariáveis -importantes, excesso de zeros, diferentes amplitudes de domÃnio -(\textit{offset}) não consideradas, heterogeneidade de unidades +superdispersos. A superdispersão pode ocorrer pela ausência de +covariáveis importantes, excesso de zeros, diferentes amplitudes de +domÃnio (\textit{offset}) não consideradas, heterogeneidade de unidades amostrais, entre outros \cite{RibeiroJr2012}. Para tais casos, uma -abordagem é a adoção de modelos com efeitos aleatórios que capturam a -variabilidade extra com a adoção de um ou mais termos de efeito +abordagem é a adoção de modelos com efeitos aleatórios, que capturam a +variabilidade extra, com a adoção de um ou mais termos de efeito aleatório. Um caso particular do modelo Poisson de efeitos aleatórios, -muito adotado no campo aplicado da EstatÃstica, ocorre quando -considera-se a distribuição Gama para os efeitos aleatórios, nessa -situação há expressão fechada para a função de probabilidade -marginal, que assume a forma Binomial Negativa. +muito adotado no campo aplicado da EstatÃstica, ocorre quando a +distribuição Gama é assumida para os efeitos aleatórios. Nessa situação +há expressão fechada para a função de probabilidade marginal, que assume +a forma Binomial Negativa. Outra manifestação de fuga da suposição de equidispersão é a -subdispersão, situação menos comum na literatura. Os processos que -reduzem a variabilidade das contagens, abaixo do estabelecido pela -Poisson, não são tão conhecidos quanto os que produzem variabilidade -extra. Pela mesma razão, são poucas as abordagens descritas na -literatura capazes de tratar subdispersão, uma vez que efeitos -aleatórios só capturam a variabilidade extra. Cita-se os modelos -de quasi-verossimilhança como a abordagem mais utilizada. Todavia não é -possÃvel descrever uma distribuição de probabilidades para a variável -resposta nessa abordagem, pois a modelagem é baseada apenas nos dois -primeiros momentos da distribuição condicional \cite{Paula2013}. +subdispersão, situação menos comum na prática e menos relatada na +literatura. Os processos que reduzem a variabilidade das contagens, +abaixo do estabelecido pela Poisson, não são tão conhecidos quanto os +que produzem variabilidade extra. Pela mesma razão, são poucas as +abordagens descritas na literatura capazes de tratar subdispersão, uma +vez que efeitos aleatórios só capturam a variabilidade extra. Cita-se os +modelos de quasi-verossimilhança como a abordagem mais +utilizada. Todavia não é possÃvel descrever uma distribuição de +probabilidades para a variável resposta nessa abordagem, pois a +modelagem é baseada apenas nos dois primeiros momentos da distribuição +condicional \cite{Paula2013}. <<processo-pontual, fig.cap="Ilustração de diferentes tipos de processos pontuais. Da direita para esquerda têm-se processos sob padrões aleatório, aglomerado e uniforme.", fig.height=3, fig.width=7>>= @@ -132,21 +134,21 @@ xyplot(y ~ x | caso, data = da, @ -A figura \ref{fig:processo-pontual} ilustra, em duas dimensões, a -ocorrência de equi, super e subdispersão respectivamente. Nesta figura +A \autoref{fig:processo-pontual} ilustra, em duas dimensões, a +ocorrência de equi, super e subdispersão respectivamente. Nessa figura cada ponto representa a ocorrência de um evento e cada parcela, delimitada pelas linhas pontilhadas, representa a unidade (ou domÃnio) -na qual tem-se o número de eventos (como variável aleatória). O painel -da esquerda representa a situação de dados de contagem equidispersos, -nesse cenário as ocorrências da variável aleatória se dispõem -aleatoriamente. No painel central o padrão já se altera, tem-se a -representação do caso de superdispersão. Nesse cenário formam-se -aglomerados que deixam parcelas com contagens muito elevadas e parcelas -com contagens baixas. Uma possÃvel causa deste padrão se dá pelo -processo de contágio (e.g. contagem de casos de uma doença contagiosa, -contagem de frutos apodrecidos). No terceiro e último painel ilustra-se -o caso de subdispersão, em que as ocorrências se dispõem uniformemente -no espaço. Agora as contagens de ocorrências nas parcelas variam bem +na qual conta-se o número de eventos (como variável aleatória). O painel +da esquerda representa a situação de dados de contagem equidispersos. +Nesse cenário as ocorrências dos eventos se dispõem aleatoriamente. No +painel central o padrão já se altera, tem-se a representação do caso de +superdispersão. Nesse cenário formam-se aglomerados que deixam parcelas +com contagens muito elevadas e parcelas com contagens baixas. Uma +possÃvel causa desse padrão se dá pelo processo de contágio +(e.g. contagem de casos de uma doença contagiosa, contagem de frutos +apodrecidos). No terceiro e último painel ilustra-se o caso de +subdispersão, em que as ocorrências se dispõem uniformemente no +espaço. Agora as contagens de ocorrências nas parcelas variam bem pouco. Ao contrário do caso superdisperso uma causa provável seria o oposto de contágio, a repulsa, ou seja, uma ocorrência causa a repulsa de outras ocorrências em seu redor (e.g. contagem de árvores, contagem @@ -154,7 +156,7 @@ de animais territoriais ou que disputam por território). Uma alterativa paramétrica que contempla os casos de equi, super e subdispersão é a adoção de uma distribuição mais flexÃvel para a -variável resposta condicional as covariáveis. \citeonline{Conway1962}, +variável resposta condicional à s covariáveis. \citeonline{Conway1962}, antes da formalização dos MLG's, propuseram uma distribuição denominada COM-Poisson (nome em em homenagem aos seus autores Richard W. Conway, William L. Maxwell, \textbf{Co}nway-\textbf{M}axwell-Poisson) que @@ -165,10 +167,10 @@ contemplando os casos de sub e superdispersão \cite{Shmueli2005}. Uma caracterÃstica bastante relevante é que a COM-Poisson possui como casos particulares as distribuições Poisson, Geométrica e Binomial. Portanto, empregando a COM-Poisson como distribuição -condicional de um modelo de regressão, a imposição de equidispersão não +condicional em um modelo de regressão, a imposição de equidispersão não precisa ser satisfeita. Tal flexibilidade, considerando o amplo uso do modelo Poisson, significa que a COM-Poisson pode ser aplicada nessas -situações e será especialmente importante naquelas onde há fuga da +situações e será especialmente importante naquelas em que há fuga da equidispersão. Assim como no modelo COM-Poisson vários aspectos do COM-Poisson podem @@ -177,42 +179,43 @@ experimento sugere uma estrutura de covariância entre observações induzidas por um processo hierárquico de casualização ou amostragem. São casos assim os experimentos em parcelas subdivididas e experimentos com medidas repetidas ou longitudinais. Tais estruturas estabelecem modelos -com efeitos não observáveis que agem em grupos experimentais e isso pode -ser incorporado no modelo de regressão COM-Poisson com a inclusão de -efeitos aleatórios. Da mesma forma, excesso de zeros pode ser -introduzido a essa distribuição da mesma maneira que ocorre para o -modelo Poisson, através de truncamento (modelos Hurdle) ou inflação -(modelos de mistura) \cite{Sellers2016}. Estas extensões do modelo -COM-Poisson ainda não são bem consolidadas na literatura e são escassas -suas aplicações. Uma constatação do fato é que não há implementações -destas extensões nos principais softwares estatÃsticos. +com efeitos não observáveis e isso pode ser incorporado no modelo de +regressão COM-Poisson com a inclusão de efeitos aleatórios a nÃvel de +grupos experimentais. Da mesma forma, excesso de zeros pode ser +introduzido a essa distribuição como ocorre para o modelo Poisson, +através de truncamento (modelos Hurdle) ou inflação (modelos de mistura) +\cite{Sellers2016}. Estas extensões do modelo COM-Poisson ainda não são +bem consolidadas na literatura e são escassas suas aplicações. Uma +constatação do fato é que não há implementações destas extensões nos +principais softwares estatÃsticos. Na literatura brasileira, aplicações do modelo COM-Poisson são -escassas. Foram encontradas apenas aplicações na área de Análise de +raras. Foram encontradas apenas aplicações na área de Análise de Sobrevivência, mais especificamente em modelos com fração de cura \cite{Ribeiro2012, Borges2012}. Portanto, o presente trabalho visa colaborar com a literatura estatÃstica brasileira i) apresentando e -explorando o modelo de regressão COM-Poisson para dados de contagem, ii) +explorando o modelo de regressão COM-Poisson para dados de contagem; ii) estendendo as aplicações desse modelo para situações especÃficas como -inclusão de efeitos aleatórios e modelagem de excesso de zeros, iii) -discutindo os aspectos inferenciais por meio de análise de dados reais e -iv) disponibilizando os recursos computacionais, em formato de pacote R, -para ajuste dos modelos apresentados. Nas aplicações optou-se também -pela análise via modelos já disponÃveis para as situações estudadas. - -O trabalho é organizado em cinco capÃtulos. Esse primeiro capÃtulo visa +inclusão de efeitos aleatórios e modelagem de excesso de zeros; iii) +discutindo os aspectos inferenciais por meio de análise de dados reais; +e iv) disponibilizando os recursos computacionais, em formato de pacote +R, para ajuste dos modelos apresentados. Nas aplicações optou-se também +pela análise via modelos Poisson, Quasi-Poisson e Binomial Negativa para +comparação de resultados. + +O trabalho é organizado em cinco capÃtulos. O primeiro capÃtulo visa enfatizar as caracterÃsticas das variáveis aleatórias de contagem e suas lacunas que podem ser complementadas na análise estatÃstica dessas -variáveis. O capÃtulo \ref{cap:modelos-para-dados-de-contagem} é -dedicado a revisão bibliográfica dos modelos estatÃsticos empregados a -análise de dados de contagem. Nesse capÃtulo os modelos Poisson, -Binomial Negativo, COM-Poisson, as abordagens para excesso de zeros e a -estrutura dos modelos de efeitos aleatórios são apresentados. No -capÃtulo \ref{cap:material-e-metodos} são apresentos os conjuntos de -dados a serem analisados e os métodos para ajuste e comparação dos -modelos. O capÃtulo \ref{cap:resultados-e-discussao} traz os os -principais resultados da aplicação e comparação dos modelos estatÃsticos -com ênfase nas discussões sob aspectos inferenciais -empÃricos. Finalmente no capÃtulo \ref{cap:consideracoes-finais} são -apresentadas as considerações finais obtidas desse trabalho e listados -algumas possÃveis linhas de pesquisa para estudos futuros. +variáveis. O \autoref{cap:modelos-para-dados-de-contagem} é dedicado a +revisão bibliográfica dos modelos estatÃsticos empregados à análise de +dados de contagem. Nesse capÃtulo os modelos Poisson, Binomial Negativo, +COM-Poisson, as abordagens para excesso de zeros e a estrutura dos +modelos de efeitos aleatórios são apresentados. No +\autoref{cap:material-e-metodos} são apresentados os conjuntos de dados +a serem analisados e os métodos para ajuste e comparação dos modelos. O +\autoref{cap:resultados-e-discussao} traz os principais resultados da +aplicação e comparação dos modelos estatÃsticos, com ênfase nas +discussões sob aspectos inferenciais empÃricos. Finalmente, no +\autoref{cap:consideracoes-finais} são apresentadas as considerações +finais obtidas desse trabalho e listadas algumas possÃveis linhas de +pesquisa para estudos futuros. diff --git a/docs/cap02_revisao-de-literatura.Rnw b/docs/cap02_revisao-de-literatura.Rnw index 4d647741b534ce16a55d963d3ed794e913939fdd..70cd424db1942a4d6d8bb1dc5b2a01262ca34c14 100644 --- a/docs/cap02_revisao-de-literatura.Rnw +++ b/docs/cap02_revisao-de-literatura.Rnw @@ -7,40 +7,40 @@ quantidade disponÃvel para dados contÃnuos. Destaca-se o modelo log-linear Poisson como o modelo mais utilizado quando se trata de dados de contagem. Porém, não raramente os dados de contagens apresentam variância superior ou inferior à sua média. Esses são os casos de super -ou subdispersão já enunciados no capÃtulo \ref{cap:introducao}, que -quando ocorrem inviabilizam o uso da distribuição Poisson. +ou subdispersão já enunciados no \autoref{cap:introducao} que, +quando ocorrem, inviabilizam o uso da distribuição Poisson. -Nos casos de fuga da equidispersão algumas abordagens não paramétricas +Nos casos de fuga da equidispersão algumas abordagens semi-paramétricas são empregadas. Nesse contexto, são alternativas os métodos de estimação via quase-verossimilhança, estimação robusta dos erros padrões (estimador ``sanduÃche'') e estimação dos erros padrões via reamostragem (``\textit{bootstrap}'') \cite{Hilbe2014}. Desses métodos detalha-se, brevemente, somente o método de estimação via função de -quase-verossimilhança na seção -\ref{cap02:estimacao-via-quase-verossimilhanca}. +quase-verossimilhança na +\autoref{cap02:estimacao-via-quase-verossimilhanca}. No contexto paramétrico, pesquisas recentes trazem modelos bastante flexÃveis à fuga de equidispersão no campo da EstatÃstica aplicada, veja -\citeonline{Sellers2010, Zeviani2014, Lord2010}. Na tabela -\ref{tab:distribuicoes} são listadas as distribuições de probabilidades -consideradas por \citeonline{Winkelmann2008} e +\citeonline{Sellers2010, Zeviani2014, Lord2010}. No +\autoref{quad:distribuicoes} são listadas as distribuições de +probabilidades consideradas por \citeonline{Winkelmann2008} e \citeonline{Kokonendji2014} e as caracterÃsticas de dados de contagem -que são contempladas. Nota-se que a Poisson na verdade é um caso -particular, pois é a única das distribuições listadas que contempla -somente a caracterÃstica de equidispersão, ainda observa-se um -conjunto maior de distribuições para os casos de superdispersão com -relação os casos de subdispersão. Embora este grande número de -distribuições exista para lidar com os casos de fuga de equidispersão, -são poucos os pacotes estatÃsticos que as disponibilizam como -alternativas para ajuste de modelos de regressão para dados de contagem. +que são contempladas. Nota-se que a Poisson na verdade é a única das +distribuições listadas que contempla somente a caracterÃstica de +equidispersão. Observa-se um conjunto maior de distribuições para +os casos de superdispersão com relação aos casos de subdispersão. Embora +esse grande número de distribuições exista para lidar com os casos de +fuga de equidispersão, são poucos os pacotes estatÃsticos que as +disponibilizam como alternativas para ajuste de modelos de regressão +para dados de contagem. %%---------------------------------------------------------------------- %% Tabela das distribuições para dados de contagem -\begin{table} +\begin{quadro} \centering \caption{Distribuições de probabilidades para dados de contagem com indicação das caracterÃsticas contempladas} -\label{tab:distribuicoes} +\label{quad:distribuicoes} \begin{tabular}{lccc} \toprule \multirow{2}{*}{Distribuição} & \multicolumn{3}{c}{Contempla a caracterÃstica de} \\ @@ -64,17 +64,17 @@ Katz & \checkmark & \checkmark & \chec \small \item Fonte: Elaborado pelo autor. \end{tablenotes} -\end{table} +\end{quadro} %%---------------------------------------------------------------------- Dos modelos paramétricos, o Binomial Negativo aparece em destaque com implementações já consolidadas nos principais \textit{softwares} estatÃsticos e frequentes aplicações nos casos de superdispersão. Na -seção \ref{cap02:binomneg} detalhes da construção desses modelos são -apresentados. Dos demais modelos derivados das distribuições listadas na -tabela \ref{tab:distribuicoes} este trabalho abordará somente o -modelo COM-Poisson, que é apresentado com detalhes na seção -\ref{cap02:compoisson}. +\autoref{cap02:binomneg} detalhes da construção desses modelos são +apresentados. Dos demais modelos derivados das distribuições listadas no +\autoref{quad:distribuicoes} este trabalho abordará somente o modelo +COM-Poisson, que é apresentado com detalhes na +\autoref{cap02:compoisson}. Um outro fenômeno que é frequente em dados de contagem é a ocorrência excessiva de zeros. Esse fenômeno sugere a modelagem de dois processos @@ -82,44 +82,42 @@ geradores de dados, o gerador de zeros extra e o gerador das contagens. Existem ao menos duas abordagens pertinentes para estes casos que são os modelos de mistura e os modelos condicionais. Na abordagem por modelos de mistura a variável resposta é modelada como uma mistura -de duas distribuições, no trabalho de \citeonline{Lambert1992}, +de duas distribuições. \citeonline{Lambert1992} apresenta uma mistura da distribuição Bernoulli com uma distribuição de Poisson ou Binomial Negativa. Considerando os modelos condicionais, também chamados de modelos de barreira \cite{Ridout1998}, tem-se que a modelagem da variável resposta é realizada em duas etapas. A primeira refere-se ao processo gerador de contagens nulas e a segunda ao gerador de contagens não nulas. Nesse trabalho a modelagem de excesso de zeros se dará -somente via modelos de barreira. A seção \ref{cap02:zeros} é destinada a +somente via modelos de barreira. A \autoref{cap02:zeros} é destinada a um breve detalhamento desta abordagem. Neste capÃtulo também é abordada a situação da inclusão de efeitos -aleatórios na seção \ref{cap02:aleatorio}. Em análise de dados de -contagem a inclusão desses efeitos permitem acomodar variabilidade -extra e incorporar a estrutura amostral do problema como em experimentos -com medidas repetidas ou longitudinais e experimentos em parcelas +aleatórios na \autoref{cap02:aleatorio}. Em análise de dados de contagem +a inclusão desses efeitos permitem acomodar variabilidade extra e +incorporar a estrutura amostral do problema, como em experimentos com +medidas repetidas ou longitudinais e experimentos em parcelas subdivididas. \section{Modelo Poisson} \label{cap02:poisson} -A Poisson é uma das principais distribuição de probabilidades +A Poisson é uma das principais distribuições de probabilidades discretas. Com suporte nos inteiros não negativos, uma variável aleatória segue um modelo Poisson se sua função massa de probabilidade for - \begin{equation} \label{eqn:pmf-poisson} - \Pr(Y = y \mid \lambda) = \frac{\lambda^ye^{-\lambda}}{y!} - \qquad y = 0, 1, 2, \ldots + \Pr(Y = y \mid \lambda) = \frac{\lambda^ye^{-\lambda}}{y!}, + \qquad y = 0, 1, 2, \ldots \end{equation} - -\noindent -em que $\lambda > 0$ representa a taxa de ocorrência do evento. Uma -particularidade já destacada desta distribuição é que $E(X) = V(X) = -\lambda$. Isso torna a distribuição Poisson bastante restritiva. Na -figura \ref{fig:distr-poisson} são apresentadas as distribuições Poisson -para diferentes parâmetros, note que devido a propriedade $E(X) = V(X)$ -contagens maiores também são mais dispersas. +em que $\lambda > 0$ representa a taxa de ocorrência do +evento. Uma particularidade já destacada desta distribuição é que $E(X) += V(X) = \lambda$. Isso torna a distribuição Poisson bastante +restritiva. Na \autoref{fig:distr-poisson} são apresentadas as +distribuições Poisson para diferentes parâmetros. Note que, devido a +propriedade $E(X) = V(X)$, contagens de médias maiores também tem +probabilidades mais dispersas. <<distr-poisson, fig.height=3.3, fig.width=6.7, fig.cap="Probabilidades pela distribuição Poisson para diferentes parâmetros.">>= @@ -147,39 +145,37 @@ distribuição Exponencial. Essa relação estabelece que se os tempos entre a ocorrência de eventos se distribuem conforme modelo Exponencial de parâmetro $\lambda$ a contagem de eventos em um intervalo de tempo $t$ tem distribuição Poisson com média $\lambda t$. A distribuição -\textit{Gamma-Count}, citada na tabela \ref{tab:distribuicoes}, estende -esta propriedade do processo adotando a distribuição Gama para os tempos -entre eventos tornando a distribuição da contagem decorrente mais +\textit{Gamma-Count}, citada no \autoref{quad:distribuicoes}, estende +essa propriedade do processo adotando a distribuição Gama para os tempos +entre eventos, tornando a distribuição da contagem decorrente mais flexÃvel \cite{Winkelmann1995, Zeviani2014}. Outra propriedade que decorre da construção do modelo Poisson é sobre a razão entre probabilidades sucessivas, $\frac{\Pr(Y=y-1)}{\Pr(Y=y)} = \frac{y}{\lambda}$. Essa razão é linear em $y$ e tem sua taxa de -crescimento ou decrescimento como $\frac{1}{\lambda}$. Os modelos Katz e +variação instantânea igual a $\frac{1}{\lambda}$. Os modelos Katz e COM-Poisson se baseiam na generalização dessa razão de probabilidades a fim de flexibilizar a distribuição de probabilidades. A utilização do modelo Poisson na análise de dados se dá por meio do -modelo de regressão Poisson. Seja $Y_i$ variáveis aleatórias -condicionalmente independentes, dados as covariáveis $X_i$, -$i=1,2,\ldots,n$. O modelo de regressão log-linear Poisson, sob a teoria -dos MLG's é definido como - +modelo de regressão Poisson. Sejam $Y_1, Y_2, \ldots, Y_n$ variáveis +aleatórias condicionalmente independentes, dado o vetor de covariáveis +$\underline{x}_i^t=(x_{i1},x_{i2},\ldots,x_{ip})$. O modelo de regressão +log-linear Poisson, sob a teoria dos MLG's, é definido como \begin{equation} \label{eqn:reg-poisson} \begin{split} - Y_i \mid & X_i \sim \textrm{Poisson}(\mu_i) \\ - &\log(\mu_i) = X_i\beta + Y_i \mid & \underline{x}_i \sim \textrm{Poisson}(\mu_i) \\ + &\log(\mu_i) = \underline{x}_i^t\beta \end{split} \end{equation} - -\noindent -em que $\mu_i > 0$ é a média da variável aleatória $Y_i \mid X_i$ que é -calculada a partir do vetor $\beta \in \mathbb{R}^p$. +em que $\mu_i > 0$ é a média da variável aleatória $Y_i$ condicionada ao +vetor de covariáveis $\underline{x}_i^t$, que é calculada a partir do +vetor $\beta \in \mathbb{R}^p$. O processo de estimação do vetor $\beta$ é baseado na maximização da -verossimilhança, que nas distribuições pertencentes à famÃlia -exponencial, os MLG's, é realizada via algoritmo de mÃnimos quadrados +função de verossimilhança que, nas distribuições pertencentes à famÃlia +exponencial, é realizada via algoritmo de mÃnimos quadrados ponderados iterativamente, ou, do inglês \textit{Iteractive Weighted Least Squares - IWLS} \cite{Nelder1972}. @@ -188,50 +184,49 @@ ponderados iterativamente, ou, do inglês \textit{Iteractive Weighted \citeonline{Wedderburn1974} propôs uma forma de estimação a partir de uma função biparamétrica, denominada quase-verossimilhança. Suponha -$y_i$ observações independentes com esperanças $\mu_i$ e variâncias -$V(\mu_i)$, em que $V$ é uma função positiva e conhecida. A função de -quase-verossimilhança é expressa como - +$Y_1, Y_2, \ldots, Y_n$ variáveis aleatórias independentes com $E(Y_i) = +\mu_i$ e variâncias $V(\mu_i)$, em que $V$ é uma função +positiva e conhecida. A função de quase-verossimilhança é expressa como \begin{equation} \label{eqn:quase-verossimilhanca} - Q(\mu_i \mid y_i) = \int_{y_i}^{\mu_i} \frac{y_i - t}{\sigma^2 V(\mu_i)}dt + Q(\mu_i \mid y_i) = \int_{y_i}^{\mu_i} + \frac{y_i - \mu_i^t}{\sigma^2 V(\mu_i)}d\mu_i^t \end{equation} -Na expressão \ref{eqn:quase-verossimilhanca} a função de -quase-verossimilhança é definida a partir da especificação de $\mu_i$, -$V(\mu_i)$ e $\sigma^2$. O processo de estimação via maximização dessa -função compartilha, do método baseado na maximazação da verossimilhança, -as mesmas estimativas para $\mu_i$, porém a dispersão de $y_i$, $V(y_i) -= \theta V(\mu_i)$ é corrigida pelo parâmetro adicional $\sigma^2$. - -Assim os problemas com a fuga da suposição de equidispersão podem ser -superados quando a estimação por máxima quase-verossimilhança é -adotada. Porém um resultado dessa abordagem é que - +Na \autoref{eqn:quase-verossimilhanca} a função de quase-verossimilhança +é definida a partir da especificação de $\mu_i$, $V(\mu_i)$ e +$\sigma^2$. O processo de estimação via maximização dessa função +compartilha, do método baseado na maximazação da função de +verossimilhança, as mesmas estimativas para $\mu_i$, porém a dispersão +de $y_i$ é corrigida pelo parâmetro adicional $\sigma^2$, $V(y_i) = +\sigma^2 V(\mu_i)$. + +Com a adição desse parâmetro de dispersão $\sigma^2$, os problemas com a +fuga da suposição de equidispersão são superados. Porém um resultado +dessa abordagem é que \begin{equation} \label{eqn:quasi-informacao} -E\left ( \frac{\partial^2 Q(\mu \mid y)}{\partial \mu^2} \right) \leq -E\left ( \frac{\partial^2 \ell(\mu \mid y)}{\partial \mu^2} \right) \end{equation} - -\noindent -ou seja a informação a respeito de $\mu$ quando se conhece apenas +ou seja, a informação a respeito de $\mu$ quando se descreve apenas $\sigma^2$ e $V(\mu)$, a relação entre média e variância, é menor do que -a informação quando se conhece a distribuição da variável resposta, dada -pela log-verossimilhança $\ell(\mu \mid y)$. Além disso ressalta-se que, -de forma geral, não é possÃvel descrever uma distribuição de +a informação quando se descreve a distribuição da variável resposta, +dada pela log-verossimilhança $\ell(\mu \mid y)$. Além disso ressalta-se +que, de forma geral, não é possÃvel descrever uma distribuição de probabilides para $Y$ somente com as especificações de $\sigma^2$ e $V(\mu)$. -Em modelos de regressão, $g(\mu_i) = X\beta$ e $V(\mu_i)$ definem a -função de quase-verossimilhança. Nessa abordagem são estimados os -parâmetros $\beta$ e $\sigma^2$. A estimativa do vetor $\beta$ pode ser -obtidas pelo algoritmo \textit{IWLS}. Usando as funções quase-escore e -matriz de quase-informação chega-se ao mesmo algoritmo de estimação dado -no caso Poisson, que não depende de $\sigma^2$. O parâmetro $\sigma^2$ é -estimado separadamente, pós estimação dos $\beta$'s. Um estimador usual -é o baseado na estatÃstica $\chi^2$ de Pearson. - +Em modelos de regressão, $g(\mu_i) = \underline{x}_i^t\beta$ e +$V(\mu_i)$ definem a função de quase-verossimilhança. Nessa abordagem +são estimados os parâmetros $\beta$ e $\sigma^2$. A estimação do vetor +$\beta$ pode ser realizada pelo algoritmo \textit{IWLS}. Usando as +funções quase-escore, derivadas de primeira ordem da função $Q(\mu_i +\mid y_i)$ em relaçao aos $\beta$'s, e matriz de quase-informação, +derivadas de segunda ordem, chega-se ao mesmo algoritmo de estimação +dado no caso Poisson, que não depende de $\sigma^2$. O parâmetro +$\sigma^2$ é estimado separadamente, pós estimação dos $\beta$'s. Um +estimador usual é o baseado na estatÃstica $\chi^2$ de Pearson \begin{equation} \label{eqn:estimador-theta} \hat{\sigma^2} = \frac{1}{n-p} \sum_{i=1}^n @@ -241,11 +236,10 @@ estimado separadamente, pós estimação dos $\beta$'s. Um estimador usual \section{Modelo Binomial Negativo} \label{cap02:binomneg} -Uma das principais alternativas paramétricas para dados de contagem -superdispersos é a distribuição Binomial Negativa. A função massa de -probabilidade da distribuição Binomial Negativa pode ser deduzida de um -processo hierárquico de efeitos aleatórios onde se assume que - +Uma das principais distribuições paramétricas para dados de contagem +superdispersos é a Binomial Negativa. A função massa de probabilidade da +distribuição Binomial Negativa pode ser deduzida de um processo +hierárquico de efeitos aleatórios em que se assume \begin{equation} \label{eqn:proc-binomneg} \begin{split} @@ -254,44 +248,6 @@ processo hierárquico de efeitos aleatórios onde se assume que \end{split} \end{equation} -\noindent -A função massa de probabilidade decorrente da estrutura descrita em -\ref{eqn:proc-binomneg} é deduzida integrando os efeitos aleatórios. -Considere $f(y \mid b)$ como a função massa de probabilidade da -distribuição Poisson (vide expressão em \ref{eqn:pmf-poisson}) e $g(b -\mid \mu, \phi)$ a função densidade da distribuição Gama \footnote{O - desenvolvimento detalhado da integral pode ser visto em - \citeonline[pág. 303-305]{Paula2013}. Obs.: A função densidade do - modelo Gama está parametrizada para que $\mu$ represente a média da - distribuição.} - -\begin{equation} - \label{eqn:proc-binomneg} - \begin{split} - \Pr(Y = y \mid \mu,\theta) &= \int_0^\infty f(y \mid b) - g(b \mid \mu,\theta) db\\ - &= \frac{\theta^\theta}{y!\mu^\theta\Gamma(\theta)} - \int_0^\infty e^{-b(1 + \theta/\mu)} b^{y+\theta-1}db \\ - &= \frac{\Gamma(\theta + y)}{\Gamma(y+1)\Gamma(\theta)} - \left ( \frac{\mu}{\mu + \theta} \right )^y - \left ( \frac{\theta}{\mu + \theta} \right )^\theta - \qquad y = 0, 1, 2, \cdots - \end{split} -\end{equation} - -\noindent -com $\mu >0$ e $\theta > 0$. Esse é um caso particular de um modelo de -efeito aleatório cuja integral tem solução analÃtica e por consequência -o modelo marginal tem forma fechada. Outro caso que se baseia no mesmo -princÃpio é o modelo \textit{Inverse Gaussian Poisson}, que como o nome -sugere adota a distribuição Inversa Gaussiana para os efeitos -aleatórios. Na figura \ref{fig:distr-binomneg} são apresentadas as -distribuições Binomial Negativa para diferentes parâmetros $\theta$ em -comparação com a distribuição Poisson equivalente em locação. Note que -quanto menor o parâmetro $\theta$, maior a dispersão da -distribuição. Isso introduz uma propriedade importante desse modelo, -para $\theta \to \infty$ a distribuição reduz-se a Poisson. - <<distr-binomneg, fig.height=3.4, fig.width=6.7, fig.cap="Probabilidades pela distribuição Binomial Negativa para diferentes nÃveis de dispersão, fixando a média em 5.">>= ##------------------------------------------- @@ -343,13 +299,13 @@ xyplot(values ~ c(y - 0.15) | ind, data = da.po, type = "h", col = cols[2])) for(i in 1:3){ trellis.focus("panel", i, 1, highlight=FALSE) - grid::grid.text(label = sprintf("E[Y]: %.1f\nV[Y]: %.1f", + grid::grid.text(label = sprintf("E(Y): %.1f\nV(Y): %.1f", mu, mu), x = .62, y = 0.03, default.units = "npc", gp = grid::gpar(col = cols[1]), just = c("left", "bottom")) - grid::grid.text(label = sprintf("E[Y]: %.1f\nV[Y]: %.1f", + grid::grid.text(label = sprintf("E(Y): %.1f\nV(Y): %.1f", mu, vars[i]), x = .08, y = 0.03, default.units = "npc", @@ -360,14 +316,39 @@ trellis.unfocus() @ -Os momentos média e variância da distribuição Binomial Negativa são -expressos como $E(Y) = \mu$ e $V(Y) = \mu + \mu^2/\sigma^2$. Pelas -expressões fica evidente a caracterÃstica da Binomial Negativa de -acomodar somente superdispersão, pois $E(Y)$ é menor que $V(Y)$ para -qualquer $\sigma^2$. Percebe-se também que quanto maior o parâmetro -$\sigma^2$ mais $E(Y)$ se aproxima de $V(Y)$, e no limite, quando -$\sigma^2 \rightarrow \infty$, $E(Y) = V(Y)$ fazendo com que a -distribuição Binomial Negativa se reduza a Poisson. +A função massa de probabilidade de $Y$, decorrente da estrutura descrita +na \autoref{eqn:proc-binomneg} é deduzida integrando os efeitos +aleatórios. Considere $f(y \mid b)$ como a função massa de +probabilidade da distribuição Poisson (vide \autoref{eqn:pmf-poisson}) e +$g(b \mid \mu, \phi)$ a função densidade da distribuição Gama +\footnote{O desenvolvimento detalhado da integral pode ser visto em + \citeonline[pág. 303-305]{Paula2013}. Obs.: A função densidade do + modelo Gama está parametrizada para que $\mu$ represente a média da + distribuição.} +\begin{equation} + \label{eqn:proc-binomneg} + \begin{split} + \Pr(Y = y \mid \mu,\theta) &= \int_0^\infty f(y \mid b) + g(b \mid \mu,\theta) db\\ + &= \frac{\theta^\theta}{y!\mu^\theta\Gamma(\theta)} + \int_0^\infty e^{-b(1 + \theta/\mu)} b^{y+\theta-1}db \\ + &= \frac{\Gamma(\theta + y)}{\Gamma(y+1)\Gamma(\theta)} + \left ( \frac{\mu}{\mu + \theta} \right )^y + \left ( \frac{\theta}{\mu + \theta} \right )^\theta, + \qquad y = 0, 1, 2, \cdots + \end{split} +\end{equation} +com $\mu >0$ e $\theta > 0$. Esse é um caso particular de um modelo de +efeito aleatório, cuja integral tem solução analÃtica e, por +consequência, o modelo marginal tem forma fechada. Outro caso que se +baseia no mesmo princÃpio é o modelo \textit{Inverse Gaussian Poisson}, +que, como o nome sugere, adota a distribuição Inversa Gaussiana para os +efeitos aleatórios. Na \autoref{fig:distr-binomneg} são apresentadas +as distribuições Binomial Negativa para diferentes parâmetros $\theta$ +em comparação com a distribuição Poisson, equivalente em locação. Note +que quanto menor o parâmetro $\theta$, maior a dispersão da +distribuição. Isso introduz uma propriedade importante desse modelo, +para $\theta \to \infty$ a distribuição reduz-se a Poisson. <<mv-binomneg, fig.height=4, fig.width=4, fig.cap="Relação Média e Variância na distribuição Binomial Negativa.">>= @@ -409,51 +390,55 @@ fonte("Fonte: Elaborado pelo autor.", cex = 0.95) @ -A relação funcional entre média e variância é ilustrada na figura -\ref{fig:mv-binomneg} onde são apresentadas as médias e variâncias para -$\mu$ entre 0 e 10 e $\theta$ entre 0 e 50. O comportamento dessa -relação proporciona um maior flexibilidade à distribuição em acomodar -superdispersão, uma caracterÃstica importante exibida nesta figura é que -para a Binomial Negativa se aproximar a Poisson em contagens altas o -$\theta$ deve ser extremamente grande. +Os momentos média e variância da distribuição Binomial Negativa são +dados por $E(Y) = \mu$ e $V(Y) = \mu + \mu^2/\sigma^2$. Pelas expressões +fica evidente a caracterÃstica da Binomial Negativa de acomodar somente +superdispersão, pois $E(Y)$ é menor que $V(Y)$ para qualquer +$\sigma^2$. Percebe-se também que quanto maior o parâmetro $\sigma^2$ +mais $E(Y)$ se aproxima de $V(Y)$, e no limite, quando $\sigma^2 +\rightarrow \infty$, $E(Y) = V(Y)$ fazendo com que a distribuição +Binomial Negativa se reduza à Poisson. + +A relação funcional entre média e variância é ilustrada na +\autoref{fig:mv-binomneg} em que são apresentadas as médias e variâncias +para $\mu$, entre 0 e 10, e $\theta$, entre 0 e 50. O comportamento +dessa relação proporciona uma maior flexibilidade à distribuição em +acomodar superdispersão. Uma caracterÃstica importante exibida nessa +figura é que para a Binomial Negativa se aproximar da Poisson em médias +altas o $\theta$ deve ser extremamente grande. O emprego do modelo Binomial Negativo em problemas se regressão ocorre de maneira similar aos MLG's, com exceção de que a distribuição só -pertence a famÃlia exponencial de distribuições se o parâmetro $\theta$ -for conhecido e assim o processo sofre algumas -alterações. Primeiramente, assim como na Poisson, defini-se $g(\mu_i) = -X\beta$, comumente utiliza-se a função $g(\mu_i) = -\log(\mu_i)$. Desenvolvendo a log-verossimilhança e suas funções -derivadas, função escore e matriz de informação de Fisher, mostra-se que -matriz de informação é bloco diagonal caracterizando a ortogonalidade -dos parâmetros $\beta$ de locação e $\theta$ de dispersão. Deste fato -decorre que a estimação dos parâmetros pode ser realizada em paralelo, -ou seja, estima-se o vetor $\beta$ pelo método de \textit{IWLS} e -posteriormente o parâmetro $\theta$ pelo método de Newton-Raphson, -faz-se os dois procedimentos simultaneamente até a convergência das -estimativas. +pertence à famÃlia exponencial de distribuições se o parâmetro $\theta$ +for fixado e assim o processo sofre algumas alterações. Primeiramente, +assim como na Poisson, define-se $g(\mu_i) = \underline{x}_i^t\beta$, +comumente utiliza-se a função $g(\mu_i) = \log(\mu_i)$. A partir da +log-verossimilhança e suas funções derivadas, função escore e matriz de +informação de Fisher, mostra-se que matriz de informação é bloco +diagonal caracterizando a ortogonalidade dos parâmetros $\beta$ de +locação e $\theta$ de dispersão. Desse fato decorre que a estimação dos +parâmetros pode ser realizada em paralelo, ou seja, estima-se o vetor +$\beta$ pelo algoritmo \textit{IWLS} e posteriormente o parâmetro +$\theta$ pelo método de Newton-Raphson. Os dois procedimentos são +realizados simultaneamente até a convergência das estimativas. \section{Modelo COM-Poisson} \label{cap02:compoisson} A distribuição de probabilidades COM-Poisson foi proposta por -\citeonline{Conway1962}, em um contexto de filas e generaliza a Poisson +\citeonline{Conway1962}, em um contexto de filas, e generaliza a Poisson em termos da razão de probabilidades sucessivas, como será visto -adiante. Seja $Y$ uma variável aleatória COM-Poisson, então sua função +adiante. Seja $Y$ uma variável aleatória COM-Poisson então sua função massa de probabilidade é - \begin{equation} \label{eqn:pmf-compoisson} - \Pr(Y=y \mid \lambda, \nu) = \frac{\lambda^y}{(y!)^\nu Z(\lambda, \nu)} + \Pr(Y=y \mid \lambda, \nu) = \frac{\lambda^y}{(y!)^\nu Z(\lambda, \nu)}, \qquad y = 0, 1, 2, \ldots \end{equation} - -\noindent em que $\lambda > 0$, $\nu \geq 0$ e $Z(\lambda, \nu)$ é uma constante -de normalização, calculada para que de fato \ref{eqn:pmf-compoisson} -seja uma função massa de probabilidade ($\sum_{i=1}^\infty \Pr(Y = y) = -1$). $Z(\lambda, \nu)$ é definida como se segue - +de normalização, calculada para que de fato a \autoref{eqn:pmf-compoisson} +seja uma função massa de probabilidade ($\sum_{i=0}^\infty \Pr(Y = i) = +1$). A função $Z(\lambda, \nu)$ é definida como se segue \begin{equation} \label{eqn:constante-z} Z(\lambda, \nu) = \sum_{j=0}^\infty \frac{\lambda^j}{(j!)^\nu} @@ -461,26 +446,23 @@ seja uma função massa de probabilidade ($\sum_{i=1}^\infty \Pr(Y = y) = O fato que torna a distribuição COM-Poisson mais flexÃvel é a razão entre probabilidades sucessivas - \begin{equation} \label{eqn:prob-ratio} \frac{\Pr(Y=y-1)}{\Pr(Y=y)} = \frac{y^\nu}{\lambda} \end{equation} - -\noindent -que se caracteriza não necessariamente linear em $y$, diferentemente da -Poisson, o que permite caudas mais pesadas ou mais leves à distribuição -\cite{Sellers2010}. Na figura \ref{fig:distr-compoisson} são +que se caracteriza não, necessariamente, linear em $y$, diferentemente +da Poisson, o que permite caudas mais pesadas ou mais leves à +distribuição \cite{Sellers2010}. Na \autoref{fig:distr-compoisson} são apresentadas as distribuições COM-Poisson para diferentes valores de $\lambda$ e $\nu$, em contraste com as equivalentes, em locação, distribuições Poisson. Nessa figura pode-se ver a flexibilidade desse modelo, pois i) contempla o caso de subdispersão mesmo em contagens baixas ($E(Y)=3$, painel a esquerda), a distribuição permite caudas -pesadas e consequentemente uma dispersão extra Poisson, ii) contempla +pesadas e consequentemente uma dispersão extra Poisson; ii) contempla subdispersão mesmo em contagens altas, onde na Poisson tem-se variabilidade na mesma magnitude, na COM-Poisson pode-se ter caudas mais leves concentrando as probabilidades em torno da média (painel a -direita) e iii) tem como caso particular a Poisson quando o parâmetro +direita); e iii) tem como caso particular a Poisson quando o parâmetro $\nu = 1$ (painel central). <<distr-compoisson, fig.height=3.4, fig.width=6.7, fig.cap="Probabilidades pela distribuição COM-Poisson para diferentes parâmetros.">>= @@ -542,13 +524,13 @@ xyplot(values ~ c(y - 0.15) | ind, data = da.po, type = "h", col = cols[2])) for(i in 1:3){ trellis.focus("panel", i, 1, highlight=FALSE) - grid::grid.text(label = sprintf("E[Y]: %.1f\nV[Y]: %.1f", + grid::grid.text(label = sprintf("E(Y): %.1f\nV(Y): %.1f", mus[i], mus[i]), x = .57, y = 0.03, default.units = "npc", gp = grid::gpar(col = cols[1]), just = c("left", "bottom")) -grid::grid.text(label = sprintf("E[Y]: %.1f\nV[Y]: %.1f", +grid::grid.text(label = sprintf("E(Y): %.1f\nV(Y): %.1f", mus[i], vars[i]), x = .05, y = 0.03, default.units = "npc", @@ -561,19 +543,19 @@ trellis.unfocus() Uma das vantagens do modelo COM-Poisson é que possui, além da Poisson quando $\nu = 1$, outras distribuições bem conhecidas como casos -particulares. Esses casos particulares ocorrem essencialmente devido a -forma assumida pela série infinita $Z(\lambda, \nu)$. Quando $\lambda = -1$, $Z(\lambda, \nu = 1) = e^\lambda$ e substituindo na expressão -\ref{eqn:pmf-compoisson}, tem-se a distribuição Poisson +particulares. Esses casos particulares ocorrem essencialmente devido à +forma assumida pela série infinita $Z(\lambda, \nu)$. Quando $\nu = 1$, +$Z(\lambda, \nu = 1) = e^\lambda$ e substituindo na +\autoref{eqn:pmf-compoisson}, tem-se a distribuição Poisson resultante. Quando $\nu \rightarrow \infty,\, Z(\lambda, \nu) \rightarrow 1+\lambda$ e a distribuição COM-Poisson se aproxima de uma distribuição Bernoulli com $P(Y=1)=\frac{\lambda}{1+\lambda}$. E quando $\nu = 0$ e $\lambda < 1$ $Z(\lambda, \nu)$ é uma soma geométrica que -resulta em $(1-\lambda)^{-1}$ e a expressão \ref{eqn:pmf-compoisson} se -resume a uma distribuição Geométrica com $P(Y=0)=(1-\lambda)$ -\cite{Shmueli2005}. Os três respectivos casos particulares citados são -ilustrados na figura \ref{fig:casos-particulares}, onde os parâmetros -foram escolhidos conforme restrições para redução da distribuição. +resulta em $(1-\lambda)^{-1}$ e a \autoref{eqn:pmf-compoisson} se resume +a uma distribuição Geométrica com $P(Y=0)=(1-\lambda)$ +\cite{Shmueli2005}. Os três casos particulares citados são ilustrados na +\autoref{fig:casos-particulares}, onde os parâmetros foram escolhidos +conforme restrições para redução da distribuição. <<casos-particulares, fig.height=3, fig.width=7, fig.cap="Exemplos de casos particulares da distribuição COM-Poisson.">>= @@ -617,7 +599,7 @@ xyplot(values ~ y | ind, data = da, layout = c(NA, 1), par.strip = list(lines = 2, col = "transparent"), sub = "Fonte: Elaborado pelo autor.") -distr <- c("Poisson", "Bernoulli", "Geométrica") +distr <- expression("Poisson", ""%~~%"Bernoulli", "Geométrica") ##------------------------------------------- ## http://stackoverflow.com/questions/33632344/strip-with-two-lines-title-r-lattice-plot for (i in 1:3) { @@ -636,7 +618,6 @@ for (i in 1:3) { Um inconveniente desse modelo é que os momentos média e variância não tem forma fechada. Sendo assim, devem ser calculados a partir da definição - $$ E(Y) = \sum_{y = 0}^{\infty} y \cdot p(y) \qquad \textrm{e} \qquad V(Y) = \sum_{y = 0}^{\infty} y^2 \cdot p(y) - E^2(Y) @@ -644,27 +625,23 @@ $$ \citeonline{Shmueli2005}, a partir de uma aproximação para $Z(\lambda, \nu)$, apresenta uma forma aproximada para os momentos da distribuição - \begin{equation} \label{eqn:cmp-mean-aprox} E(Y) \approx \lambda^{1/\nu} - \frac{\nu - 1}{2\nu} \qquad \textrm{e} \qquad V(Y) \approx \frac{\lambda^{1/\nu}}{\nu} \end{equation} - -\noindent os autores ressaltam que essa aproximação é satisfatória para $\nu \leq -1$ ou $\lambda > 10^\nu$. Na figura \ref{fig:mv-compoisson} é -representada a relação média e variância aproximada pelas expressões em +1$ ou $\lambda > 10^\nu$. Na \autoref{fig:mv-compoisson} é representada +a relação média e variância aproximada pelas expressões em \ref{eqn:cmp-mean-aprox}. Percebe-se que a relação é praticamente linear entre média e variância, \citeonline{Sellers2010} descrevem que essa pode ser relação pode, ainda, ser aproximada por -$\frac{1}{\nu}E(Y)$. Dessas aproximações, bem como das visualizações em -\ref{fig:distr-compoisson}, \ref{fig:casos-particulares} e -\ref{fig:mv-compoisson} deduz-se que o parâmetros $\nu$, ou -$\frac{1}{\nu}$, controla a precisão da distribuição, sendo ela -equidispersa quando $\nu = 1$, superdispersa quando $\nu < 1$ e -subdispersa quando $\nu > 1$. +$\frac{1}{\nu}E(Y)$. Dessas aproximações, bem como das visualizações na +\autoref{fig:distr-compoisson}, \autoref{fig:casos-particulares} e +\autoref{fig:mv-compoisson}, deduz-se que o parâmetro $\nu$, controla a +precisão da distribuição, sendo ela equidispersa quando $\nu = 1$, +superdispersa quando $\nu < 1$ e subdispersa quando $\nu > 1$. <<mv-compoisson, fig.height=4, fig.width=4, fig.cap="Relação Média e Variância na distribuição COM-Poisson.">>= @@ -693,7 +670,7 @@ title(xlab = expression(E(X) == lambda^{1/nu} - frac(nu-1, 2*nu)), line = 4) grid() -## Curvas da relação média e variância da Binomial Negativa +## Curvas da relação média e variância da COM-Poisson for (a in seq_along(nu)) { curve((1/nu[a])*(mu + (nu[a] - 1)/(2*nu[a])), add = TRUE, xname = "mu", col = col[a], lwd = 2) @@ -705,46 +682,44 @@ plotrix::color.legend( rect.col = col) mtext(text = expression(nu), side = 3, cex = 1.3, line = -3.5, at = 11.5) +abline(0, 1) fonte("Fonte: Elaborado pelo autor.", cex = 0.95) -wrapfigure() - @ -Embora o modelo COM-Poisson não tenha expressão fechada para a média da -distribuição pode-se utilizá-lo como modelo associado a distribuição -condicional da variável resposta de contagem. Isso é feito incorporando -um preditor linear em $\lambda$, que mesmo não representando a média, +Embora a distribuição COM-Poisson não tenha expressão fechada para a +média, pode-se utilizá-la como distribuição condicional da variável +resposta de contagem em modelos de regressão. Isso é feito incorporando +um preditor linear em $\lambda$ que, mesmo não representando a média, está associado com a locação da distribuição, ou seja, modela-se a média indiretamente nessa abordagem. O modelo de regressão é definido com as variáveis aleatórias condicionalmente independentes $Y_1, Y_2, \ldots, -Y_n$, dado o vetor de covariáveis $X_i = (x_{i1}, x_{i2}, \ldots, -x_{ip})$ seguindo um modelo COM-Poisson de parâmetros $\lambda_i = -e^{X_i\beta}$, $i = 1, 2, \ldots, n$ e $\nu$ comum a todas as -observações. Na expressão \ref{eqn:reg-poisson} o modelo é devidamente -formulado, conforme a notação de MLG's +Y_n$, dado o vetor de covariáveis $\underline{x}_i = (x_{i1}, x_{i2}, +\ldots, x_{ip})$ seguindo um modelo COM-Poisson de parâmetros $\lambda_i += e^{\underline{x}_i^t\beta}$, $i = 1, 2, \ldots, n$ e $\nu$ comum a +todas as observações. Na \autoref{eqn:reg-compoisson} o modelo é +devidamente formulado, conforme a notação de MLG's \begin{equation} \label{eqn:reg-compoisson} \begin{split} - Y_i \mid & X_i \sim \textrm{COM-Poisson}(\lambda_i, \nu) \\ - &\eta(E(Y_i \mid X_i)) = \log(\lambda_i) = X_i\beta + Y_i \mid & \underline{x}_i \sim \textrm{COM-Poisson}(\lambda_i, + \nu) \\ + &g(E(Y_i \mid \underline{x}_i)) = + \log(\lambda_i) = \underline{x}_i^t\beta \end{split} \end{equation} O algoritmo para estimação do conjunto de parâmetros $\Theta = (\nu, -\beta)$ do modelo é baseado na maximização da log-verossimilhança, que -decorrente da especificação em \ref{eqn:reg-compoisson} é - +\beta)$ do modelo é baseado na maximização da log-verossimilhança que, +decorrente da especificação em \ref{eqn:reg-compoisson}, é dada por \begin{equation} \label{eqn:loglik-compoisson} - \ell(\nu, \beta \mid \underline{y}) = \sum_i^n y_i \log(\lambda_i) - - \nu \sum_i^n \log(y!) - \sum_i^n \log(Z(\lambda_i, \nu)) + \ell(\nu, \beta \mid \underline{y}) = \sum_{i=1}^n y_i + \log(\lambda_i) - \nu \sum_{i=1}^n \log(y!) - \sum_{i=1}^n + \log(Z(\lambda_i, \nu)) \end{equation} - -\noindent e então as estimativas de máxima verossimilhança são - $$ \hat{\Theta} = (\hat{\nu}, \hat{\beta}) = \underset{(\nu,\,\beta)}{\textrm{arg max }} \ell(\nu, \beta \mid @@ -753,21 +728,6 @@ $$ <<constante-z, fig.height=3, fig.width=6.7, fig.cap="Convergência da constante de normalização da COM-Poisson para diferentes conjuntos de parâmetros.">>= -##------------------------------------------- -## Calcula Z para um c(lambda, phi) -funZ <- function(lambda, nu, maxit = 500, tol = 1e-5) { - z <- rep(NA, maxit) - j = 1 - ## - z[j] <- exp(j * log(lambda) - nu * lfactorial(j)) - ## - while (abs(z[j] - 0) > tol && j <= maxit) { - j = j + 1 - z[j] <- exp(j * log(lambda) - nu * lfactorial(j)) - } - return(cbind("j" = 0:j, "z" = c(1, z[!is.na(z)]))) -} - ##------------------------------------------- ## Parametros da distribuição pars <- list("p1" = c(1.362, 0.4), "p2" = c(8, 1), "p3" = c(915, 2.5)) @@ -801,22 +761,22 @@ xyplot(z ~ j | .id, data = da, @ -Para avaliação da log-verossimilhança em \ref{eqn:loglik-compoisson} a -constante de normalização $Z(\lambda, \nu)$, conforme definida em -\ref{eqn:constante-z} é calcula para cada observação o que +Para avaliação da log-verossimilhança, \autoref{eqn:loglik-compoisson}, +a constante de normalização $Z(\lambda, \nu)$, conforme definida em +\ref{eqn:constante-z}, é calculada para cada observação, o que potencialmente torna o processo de estimação lento. Uma ilustração do número de incrementos considerados para cálculo da constante $Z(\lambda, -\nu)$ é apresentada na figura \ref{fig:constante-z}. Nesta ilustração -foram utilizados os mesmos parâmetros definidos em -\ref{fig:distr-compoisson} e o número de incrementos -considerados para convergência \footnote{Adotou-se como critério de - convergência a iteração $j$ tal que $\lambda^j/(j!)^\nu < - 0,00001$}. de $Z(\lambda, \nu)$ foram \Sexpr{c(table(da[, ".id"]))} -nos primeiro, segundo e terceiro painéis respectivamente. +\nu)$ é apresentada na \autoref{fig:constante-z}. Nesta ilustração +foram utilizados os mesmos parâmetros das distribuições na +\autoref{fig:distr-compoisson}. O número de incrementos necessários para +convergência\footnote{Adotou-se como critério de convergência a + iteração $j$ tal que $\lambda^j/(j!)^\nu < 0,00001$} de $Z(\lambda, +\nu)$ foram \Sexpr{c(table(da[, ".id"]))} nos primeiro, segundo e +terceiro painéis respectivamente. Detalhes computacionais do algoritmo de maximização e manipulações algébricas para eficiência na avaliação da log-verossimilhança no modelo -COM-Poisson são discutidos na seção \ref{cap03:metodos}. +COM-Poisson são discutidos na \autoref{cap03:metodos}. \section{Modelos para excesso de zeros} \label{cap02:zeros} @@ -824,31 +784,33 @@ COM-Poisson são discutidos na seção \ref{cap03:metodos}. Problemas com excesso de zeros são comuns em dados de contagem. Caracteriza-se como excesso de zeros casos em que a quantidade observada de contagens nulas supera substancialmente aquela esperada -pelo modelo de contagem adotado. No caso do modelo Poisson a fração de -zeros é $e^{-\lambda}$. +pelo modelo de contagem adotado. -As contagens nulas que geram o excesso de zeros podem ser explicadas de +As contagens nulas em dados com excesso de zeros podem ser explicadas de duas formas distintas. A primeira denomina-se de zeros estruturais, quando a ocorrência de zero se dá pela ausência de determinada -caracterÃstica na população e a segunda, que zeros amostrais que -ocorrem segundo um processo gerador de dados de contagem (e.g processo +caracterÃstica na população e a segunda de zeros amostrais, que ocorrem +segundo um processo gerador de dados de contagem (e.g processo Poisson). Por exemplo, considerando o número de dias que uma famÃlia consome um determinado produto, tem-se aquelas famÃlias que não consomem o produto (zeros estruturais) e as demais famÃlias que consomem o produto, porém não o consumiram no intervalo de tempo considerado no -estudo (zeros amostrais). Assim, de forma geral são dois processos +estudo (zeros amostrais). Assim, de forma geral, são dois processos geradores de dados em uma variável aleatória de contagem com excessivos zeros. -Em geral, quando dados de contagem apresentam excessos de valores zero -também apresentarão superdispersão. Todavia, essa dispersão pode ser -exclusivamente devido ao excesso de zeros e assim os modelos -alternativos já apresentados não terão um bom desempenho. Uma ilustração -deste fato é apresentada na figura \ref{fig:ilustra-zeros}, em que foram -simulados dados com excesso de zeros e ajustado um modelo -COM-Poisson. Em ambos os casos o modelo não se ajustou adequadamente, -indicando que os excessos de zeros devem ser abordados de forma -diferente. +Em geral, quando dados de contagem apresentam excesso de zeros também +apresentarão superdispersão. Todavia, essa dispersão pode ser +exclusivamente devido ao excesso de zeros, e os modelos alternativos já +apresentados não terão um bom desempenho. Uma ilustração desse fato é +apresentada na \autoref{fig:ilustra-zeros}, em que foram simulados dados +com excesso de zeros. A simulação foi realizada de forma hierárquica, +simulando valores $y_z$ de uma variável aleatória Bernoulli de parâmetro +$\pi$ e, para $y_z=0$ armazenou-se o zero e para $y_z=1$ simulou-se de +uma distribuição Poisson de paramêtro $\lambda$. Ajustando um modelo +COM-Poisson para as duas simulações com diferentes parâmetros $\pi$ e +$\lambda$, observa-se que o modelo não se mostra adequado, indicando que +os excessos de zeros devem ser abordados de forma diferente. <<ilustra-zeros, fig.cap="Ilustração de dados de contagem com excesso de zeros.", fig.height=3, fig.width=5>>= @@ -857,14 +819,14 @@ diferente. set.seed(20124689) n <- 1000 -lambda <- 2; pi <- 0.9 -y1 <- sapply(rbinom(n, 1, pi), function(x) { - ifelse(x == 0, 0, rpois(1, lambda)) +lambda1 <- 2; pi1 <- 0.9 +y1 <- sapply(rbinom(n, 1, pi1), function(x) { + ifelse(x == 0, 0, rpois(1, lambda1)) }) -lambda <- 5; pi <- 0.85 -y2 <- sapply(rbinom(n, 1, pi), function(x) { - ifelse(x == 0, 0, rpois(1, lambda)) +lambda2 <- 5; pi2 <- 0.85 +y2 <- sapply(rbinom(n, 1, pi2), function(x) { + ifelse(x == 0, 0, rpois(1, lambda2)) }) ##------------------------------------------- @@ -896,6 +858,11 @@ key <- list( columns = 2, lines = list(lty = 1, col = cols), text = list(c("Observado", "COM-Poisson"))) +fl <- substitute( + expression(pi == p1~","~lambda == l1, + pi == p2~","~lambda == l2), + list(p1 = pi1, p2 = pi2, l1 = lambda1, l2 = lambda2) +) ##------------------------------------------- ## Gráfico @@ -905,7 +872,9 @@ xyplot(py_real ~ c(yu - 0.15) | caso, data = da, ylab = expression(Pr(Y==y)), ylim = ylim, key = key, - strip = strip.custom(factor.levels = paste("Simulação", 1:2)), + strip = strip.custom( + factor.levels = fl + ), sub = "Fonte: Elaborado pelo autor.") + as.layer(xyplot( py_dcmp ~ c(yu + 0.15) | caso, data = da, @@ -923,57 +892,52 @@ as duas principais abordagens são i) os modelos de mistura \textit{Hurdle Models}. Neste trabalho somente a abordagem via modelos condicionais será considerada. A função massa de probabilidade do modelo Hurdle é - \begin{equation} \label{eqn:pmf-hurdle} \Pr(Y = y \mid \pi, \Theta_c) = \begin{dcases*} - \pi & \text{se } y = 0,\\ + \pi &, \text{se } y = 0\,;\\ (1 - \pi) \frac{\Pr(Z = z \mid \Theta_c)}{1 - \Pr(Z = 0 \mid - \Theta_c)} & \text{se } y = 1, 2, \dots + \Theta_c)} &, \text{se } y = 1, 2, \dots \end{dcases*} \end{equation} - -\noindent em que $0<\pi<1$, representa a probabilidade de ocorrência de zeros e $\Pr(Z = z \mid \Theta_c)$ a função massa de probabilidade de uma variável aleatória de contagem $Z$, como a Poisson ou a Binomial Negativa. -Da especificação em \ref{eqn:pmf-hurdle}, os momentos média e variância -são obtidos usando as definições $E(Y) = \sum_{y=1}^\infty y \cdot -\Pr(Y=y)$ e $V(Y) = \sum_{y=1}^\infty y^2 \cdot \Pr(Y=y) - E^2(Y)$ +Da especificação em \ref{eqn:pmf-hurdle}, a média e a variância +são obtidas usando as definições $E(Y) = \sum_{y=1}^\infty y \cdot +\Pr(Y=y)$ e $V(Y) = \sum_{y=1}^\infty y^2 \cdot \Pr(Y=y) - E^2(Y)$. $$ E(Y) = \frac{E(Z)(1-\pi)}{1-\Pr(Z = 0)} \quad \textrm{e} \quad -V(Y) = \frac{1-\pi}{1-Pr(Z = 0)} \left [ E(Z) \frac{(1-\pi)}{1-\Pr(Z = +V(Y) = \frac{1-\pi}{1-\Pr(Z = 0)} \left [ E(Z) \frac{(1-\pi)}{1-\Pr(Z = 0)} \right ] $$ Para a inclusão de covariáveis, caracterizando um problema de regressão, dado que o modelo tem dois processos modela-se ambos como se segue - \begin{equation} \label{eqn:reg-hurdle} - \log \left (\frac{\pi_i}{1-\pi_i} \right ) = G_i\gamma \qquad e \qquad + \log \left (\frac{\pi_i}{1-\pi_i} \right ) = + \underline{\textsf{z}}_i^t\gamma \qquad \textrm{e} \qquad \begin{matrix} Z_i \sim D(\mu_i, \phi) \\ - g(\mu_i) = X_i\beta + g(\mu_i) = \underline{x}_i^t\beta \end{matrix} \end{equation} - -\noindent -com $i = 1, 2, \ldots, n$, $G_i$ e $X_i$ as covariáveis da i-ésima -observação consideradas para explicação da contagens nulas e não nulas -respectivamente, $D(\mu_i, \phi)$ uma distribuição de probabilidades -considerada para as contagens não nulas que pode conter ou não um -parâmetro $\phi$ adicional, se Poisson, $D(\mu_i, \phi)$ se resume a -$\textrm{Poisson}(\mu_i)$ e $g(\mu_i)$ uma função de ligação, nos casos -Poisson e Binomial Negativa considera-se $\log(\mu_i)$. O que está -implÃcito na formulação \ref{eqn:reg-hurdle} é que para a componente que -explica a geração de zeros está sendo considerada a distribuição -Bernoulli de parâmetro $\pi_i$, contudo pode-se utilizar distribuições -censuradas a direita no ponto $y=1$ para estimação desta probabilidade, +com $i = 1, 2, \ldots, n$, $\underline{\textsf{z}}_i$ e +$\underline{x}_i$ as covariáveis da i-ésima observação consideradas para +explicação da contagens nulas e não nulas respectivamente, $D(\mu_i, +\phi)$ uma distribuição de probabilidades considerada para as contagens +não nulas que pode conter ou não um parâmetro $\phi$ adicional e +$g(\mu_i)$ uma função de ligação. Nos casos Poisson e Binomial +Negativa, em geral, considera-se $g(\mu_i) = \log(\mu_i)$. O que está +implÃcito na formulação em \ref{eqn:reg-hurdle} é que para a componente +que explica a geração de zeros está sendo considerada a distribuição +Bernoulli de parâmetro $\pi_i$. Contudo pode-se utilizar distribuições +censuradas à direita no ponto $y=1$ para estimação dessa probabilidade, como explicam \citeonline{Zeileis2007}. \section{Modelos de efeitos aleatórios} @@ -981,7 +945,7 @@ como explicam \citeonline{Zeileis2007}. Nas seções anteriores os modelos que flexibilizam algumas suposições do modelo Poisson, basicamente permitindo casos não equidispersos e -modelando conjuntamente um processo gerador de zeros extra foram +modelando conjuntamente um processo gerador de zeros extra, foram explorados. Contudo, uma suposição dos modelos de regressão para dados de contagem vistos até aqui é que as variáveis aleatórias $Y_1, Y_2, \ldots, Y_n$ são condicionalmente independentes, dado o vetor de @@ -989,50 +953,49 @@ covariáveis. Porém não são raras as situações em que essa suposição não se mostra adequada. \citeonline{Ribeiro2012} cita alguns exemplos: \begin{itemize} - \item as observações podem ser correlacionadas no espaço, - \item as observações podem ser correlacionadas no tempo, + \item as observações podem ser correlacionadas no espaço; + \item as observações podem ser correlacionadas no tempo; \item interações complexas podem ser necessárias para modelar o efeito - conjunto de algumas covariáveis, + conjunto de algumas covariáveis; \item heterogeneidade entre indivÃduos ou unidades podem não ser suficientemente descrita por covariáveis. \end{itemize} Nessas situações pode-se estender a classe de modelos de regressão com a adição de efeitos aleatórios que incorporam termos baseados em variáveis -não observáveis (latentes) ao modelo, permitindo assim acomodar uma +não observáveis (latentes) ao modelo, permitindo acomodar uma fonte de variabilidade, que pode ser ou não estruturada, não prescrita pelo -modelo. De forma geral a especificação dos modelos de efeitos aleatórios -segue uma especificação hierárquica - +modelo. De forma geral os modelos de efeitos aleatórios seguem uma +especificação hierárquica \begin{equation} \label{eqn:reg-misto} \begin{split} - Y_{ij} \mid b_{i},& X_{ij} \sim \textrm{D}(\mu_{ij}, \phi) \\ - g(&\mu_{ij}) = X_{ij}\beta + Z_ib_i \\ - & b \sim \textrm{K}(\Theta_b) + Y_{ij} \mid b_{i},\,& \underline{x}_{ij} \sim + \textrm{D}(\mu_{ij}, \phi) \\ + g(&\mu_{ij}) = \underline{x}_{ij}^t\beta + \underline{z}_i^tb_i \\ + & b_i \sim \textrm{K}(\Theta_b) \end{split} \end{equation} - -\noindent para $i = 1, 2, \ldots, m$ (grupos com efeitos aleatórios comuns) e $j = 1, 2, \ldots, n$ (observações) com D$(\mu_{ij}, \phi)$, uma distribuição considerada para as variáveis resposta condicionalmente independentes, $g(\mu_{ij})$ uma função de ligação conforme definida na teoria dos -MLG's, $X_{ij}$ e $Z_{i}$ os vetores conhecidos que representam os -efeitos das covariáveis de interesse e os termos que definem os grupos -considerados como aleatórios, $b_i$ uma quantidade aleatória provida de -uma distribuição K$(\Theta_b)$. Nesses modelos uma quantidade aleatória -é somada ao preditor linear, diferentemente dos modelos de efeitos fixos -e a partir desta quantidade é possÃvel induzir uma estrutura de -dependência entre as observações. - -Como são duas quantidades aleatórias no modelo, $Y \mid X$ e $b$, a -verossimilhança para um modelo de efeito aleatório é dada integrando-se -os efeitos aleatórios +MLG's, $\underline{x}_{ij}$ e $\underline{z}_{i}$ os vetores conhecidos +que representam os efeitos das covariáveis de interesse e os termos que +definem os grupos considerados como aleatórios, $b_i$ uma quantidade +aleatória provida de uma distribuição K$(\Theta_b)$. Nesses modelos um +termo aleatório é somado ao preditor linear, diferentemente dos +modelos de efeitos fixos, e a partir deste termo é possÃvel induzir +uma estrutura de dependência entre as observações. + +Como são dois termos aleatórios no modelo, $Y_{ij}$ condicional ao vetor +de covariáveis e $b_i$, a verossimilhança para um modelo de efeito +aleatório é dada integrando-se os efeitos aleatórios \begin{equation} \label{eqn:loglik-misto} - \Ell(\beta, \phi, \Theta_b \mid \underline{y}) = \prod_{i=1}^m \int_{\R^q} + \Ell(\beta, \phi, \Theta_b \mid \underline{y}, \underline{b}) = + \prod_{i=1}^m \int_{\R^q} \left ( \prod_{j=1}^{n_i} f_D(y_{ij}, \mu, b_i)\right ) \cdot f_K(b \mid \Theta_b) db_i \end{equation} @@ -1040,30 +1003,26 @@ os efeitos aleatórios Na avaliação da verossimilhança é necessário o cálculo de $m$ integrais de dimensão $q$. Para muitos casos essa integral não tem forma analÃtica sendo necessários métodos numéricos de intergração, que são discutidos -na seção \ref{cap03:metodos}. As estimativas de máxima verossimilhança +na \autoref{cap03:metodos}. As estimativas de máxima verossimilhança são - $$ \hat{\Theta} = (\hat{\beta}, \hat{\Theta_b}) = \underset{(\beta,\,\Theta_b)}{\textrm{arg max }} \log(\Ell(\beta, \phi, -\Theta_b \mid \underline{y})) +\Theta_b \mid \underline{y}, \underline{b})) $$ -\noindent Em modelos de efeitos mistos é comum adotar como distribuição para os efeitos aleatórios uma Normal $q$-variada com média 0 e matriz de -variância e covariâncias $\Sigma$, ou seja, na especificação -\ref{eqn:reg-misto} K$(\Theta_b) = NMV_q(0, \Sigma)$. Para estes casos -os principais métodos de aproximação da integral tem desempenhos -melhores \cite{Bates2015}. +variâncias e covariâncias $\Sigma$, ou seja, na especificação +\ref{eqn:reg-misto} K$(\Theta_b) = NMV_q(0, \Sigma)$. Como mencionado anteriormente modelos de efeitos aleatórios são -candidatos a modelagem de dados superdispersos. Quando não há uma +candidatos à modelagem de dados superdispersos. Quando não há uma estrutura de delineamento experimental ou observacional pode-se incluir efeitos aleatórios ao nÃvel de observação (e então $m=n$, ou seja, os -vetores $Y$ e $b$ tem mesma dimensão). Casos particulares de modelos de -efeitos aleatórios, onde o efeito aleatório é adicionado ao nÃvel de -observação são o modelo Binomial Negativo e o \textit{Inverse Gaussian - Model}, em ambos os casos a integral, definida em -\ref{eqn:loglik-misto} tem solução analÃtica e consequentemente a -marginal em $Y$ forma fechada. +vetores $Y$ e $\underline{b}$ tem mesma dimensão). Casos particulares de +modelos de efeitos aleatórios, onde o efeito aleatório é adicionado ao +nÃvel de observação são o modelo Binomial Negativo e o \textit{Inverse + Gaussian Model}. Em ambos os casos a integral, definida a +\autoref{eqn:loglik-misto}, tem solução analÃtica e, consequentemente, a +marginal em $Y$, forma fechada. diff --git a/docs/cap03_materiais-e-metodos.Rnw b/docs/cap03_materiais-e-metodos.Rnw index 9b49b2d225a36a9b3001a408e02f6156e6ed6a63..c6c632665c74cda35f9d2b5c1cfa9913eaf1be7e 100644 --- a/docs/cap03_materiais-e-metodos.Rnw +++ b/docs/cap03_materiais-e-metodos.Rnw @@ -2,14 +2,14 @@ % CAPÃTULO 3 - MATERIAIS E MÉTODOS % ------------------------------------------------------------------------ -Essa seção é destinada a apresentação dos conjuntos de dados analisados -no trabalho e descrição dos recursos computacionais e métodos utilizados -na análise. Na seção \ref{cap03:materiais-dados} seis conjuntos de dados -com diferentes caracterÃsticas são apresentados. Os recursos -computacionais utilizados são descritos na seção -\ref{cap03:materiais-recursos}. Na última seção desse capÃtulo, -\ref{cap03:metodos}, são apresentados os métodos para ajuste, avaliação -e comparação dos modelos propostos. +Esse capÃtulo é destinado à apresentação dos conjuntos de dados +analisados no trabalho, descrição dos recursos computacionais e métodos +utilizados na análise. Na \autoref{cap03:materiais-dados} seis conjuntos +de dados com diferentes caracterÃsticas são apresentados. Os recursos +computacionais utilizados são descritos na +\autoref{cap03:materiais-recursos}. Na última seção desse capÃtulo, +\autoref{cap03:metodos}, são apresentados os métodos para ajuste, +avaliação e comparação dos modelos propostos. \section{Materias} \label{cap03:materiais} @@ -21,12 +21,12 @@ A seguir são apresentados os seis conjuntos de dados utilizados para avaliar o desempenho dos modelos COM-Poisson. Os dados em estudo são, quase em sua totalidade, resultantes de experimentos agronômicos com delineamentos balanceados, o que é uma caracterÃstica vantajosa para -avaliação do desempenho do modelo COM-Poisson quando empregado a análise +avaliação do desempenho do modelo COM-Poisson quando empregado à análise desses dados. A apresentação dos conjuntos segue a ordem de 1) descrição do experimento ou estudo em destaque, 2) definição das variáveis e suas -unidades de medidas e 3) descrição de suas caracterÃsticas, +unidades de medidas e 3) descrição das caracterÃsticas dos dados, potencialmente contempladas por modelos alternativos ao Poisson. \subsubsection{Capulhos de algodão sob efeito de desfolha artificial} @@ -42,30 +42,31 @@ niveis.est <- paste(unique(cottonBolls$est), collapse = ", ") @ -Experimento conduzido sob delineamento inteiramente casualizado com -cinco repetições em casa de vegetação com plantas de algodão -\emph{Gossypium hirsutum} submetidas a diferentes nÃveis de desfolha -artificial de remoção foliar (\Sexpr{niveis.des}), em combinação com o -estágio fenológico no qual a desfolha foi aplicada -(\Sexpr{niveis.est}). A unidade experimental foi um vaso com duas -plantas onde avaliou-se o número de capulhos produzidos ao final da -ciclo cultura \cite{Silva2012}. O experimento contou com -\Sexpr{nrow(cottonBolls)} observações das quais têm-se as informações -das variáveis número de capulhos de algodão produzidos (\texttt{ncap}), -nÃvel de desfolha de remoção foliar (\texttt{des}) e estágio fenológico -das planta na unidade experimental (\texttt{est}). - -Esse conjunto de dados já fora publicado sob a motivação da +Experimento com plantas de algodão \emph{Gossypium hirsutum} submetidas +à diferentes nÃveis de desfolha artificial de remoção foliar, (0, 25, +50, 75 e 100\%), em combinação com o estágio fenológico no qual a +desfolha foi aplicada, (vegetativo, botão floral, florescimento, maça e +capulho). Esse experimento foi conduzido sob delineamento interamente +casualizado com cinco repetições, em casa de vegetação. A unidade +experimental foi um vaso com duas plantas, onde avaliou-se o número de +capulhos produzidos ao final da ciclo cultura \cite{Silva2012}. O +experimento contou com \Sexpr{nrow(cottonBolls)} observações das quais +têm-se as informações das variáveis número de capulhos de algodão +produzidos (\texttt{ncap}), nÃvel de desfolha de remoção foliar +(\texttt{des}) e estágio fenológico das plantas na unidade experimental +(\texttt{est}). + +Esse conjunto de dados já fora analisado e publicado sob a motivação da caracterÃstica de subdispersão, utilizando o modelo \textit{Gamma-Count} -\cite{Zeviani2014}. Na figura \ref{fig:descr-cottonBolls}, são -apresentados os dados do experimento. À esquerda apresenta-se a -disposição das cinco observações em cada tratamento (combinação de nÃvel -de desfolha e estágio fenológico do algodão) e à direita um gráfico -descritivo cruzando médias e variâncias amostrais calculadas em cada -tratamento, onde a linha pontilhada representa a caracterÃstica de -equidispersão, média igual a variância. Em todos os tratamentos -obteve-se a média menor que a variância apontando evidência de -subdispersão. +\cite{Zeviani2014}. Na \autoref{fig:descr-cottonBolls}, são apresentados +os dados do experimento. À esquerda apresenta-se a disposição das cinco +observações em cada tratamento (combinação de nÃvel de desfolha e +estágio fenológico do algodão) e à direita um gráfico descritivo +cruzando médias e variâncias amostrais calculadas em cada tratamento, +onde a linha pontilhada representa a caracterÃstica de equidispersão, +média igual a variância, e a contÃnua a reta de um ajuste de regressão +linear simples. Em todos os tratamentos obteve-se a média menor que a +variância apontando evidência de subdispersão. <<descr-cottonBolls, fig.height=4.2, fig.width=7, fig.cap="Número de capulhos produzidos para cada nÃvel de desfolha e estágio fenológico (esquerda) e médias e variâncias das cinco repetições em cada combinação de nÃvel de desfolha e estágio fenológico (direita).">>= @@ -124,15 +125,15 @@ Experimento conduzido na Universidade Federal da Grande Dourados (UFGD) em 2007, cujo objetivo foi avaliar os impactos da exposição de plantas à alta infestação de Mosca-branca \emph{Bemisia tabaci} em componentes de produção do algodão \cite{Martelli2008}. No experimento, plantas de -algodão foram expostas à alta infestação da praga por diferentes -perÃodos, \Sexpr{niveis.dexp} onde avaliou-se o número de capulhos +algodão foram expostas a alta infestação da praga por diferentes +perÃodos, 0, 1, 2, 3, 4, e 5 dias. Avaliou-se o número de capulhos produzidos (\texttt{ncapu}), o número de estruturas reprodutivas (\texttt{nerep}) e o número de nós (\texttt{nnos}), como variáveis de interesse que representam a produtividade do cultivo de algodão. A condução do estudo deu-se via delineamento inteiramente casualizado com cinco vasos contendo duas plantas, para cada perÃodo de exposição. -<<descr-cottonBolls2, fig.height=3.5, fig.width=7, fig.cap="Disposição das variáveis de contagem nº de estruturas reprodutivas, nº de capulhos produzidos e nº de nós da planta observadas sob diferentes dias de exposição à infestação de Mosca-branca.">>= +<<descr-cottonBolls2, fig.height=3.5, fig.width=7.2, fig.cap="Disposição das variáveis de contagem nº de estruturas reprodutivas, nº de capulhos produzidos e nº de nós da planta observadas sob diferentes dias de exposição à infestação de Mosca-branca.">>= vars <- c("dexp", "vaso", "planta", "nerep", "ncapu", "nnos") cottonBolls2 <- cottonBolls2[, vars] @@ -141,6 +142,7 @@ da <- reshape2::melt(cottonBolls2, id = c("dexp", "vaso", "planta"), variable.name = "va", value.name = "count") ## da <- aggregate(count ~ vaso + dexp + va, data = da, FUN = sum) +da$va <- relevel(da$va, "ncapu") xyplot(count ~ dexp | va, data = da, type = c("p", "g", "smooth"), layout = c(NA, 1), @@ -152,8 +154,8 @@ xyplot(count ~ dexp | va, data = da, column = 2, title = "Planta", cex.title = 1, lines = TRUE), strip = strip.custom( - factor.levels = c("Estruturas reprodutivas ", - "Capulhos produzidos", + factor.levels = c("Capulhos produzidos ", + "Estruturas reprodutivas", "Nós da planta") ), spread = 0.15, @@ -164,22 +166,34 @@ fonte.xy("Fonte: Elaborado pelo autor.") @ -Na figura \ref{fig:descr-cottonBolls2} a disposição de cada uma das -variáveis aleatórias de contagem número de estruturas reprodutivas, -número de capulhos produzidos e número de nós da planta para os -diferentes perÃodos em que as plantas estiveram sob alta infestação de -Mosca-branca é apresentada. Para todas as contagens parece haver um -comportamento subdisperso. A indicação de subdispersão também se observa -na tabela \ref{tab:mv-cottonBolls2}, onde as médias e variâncias -amostrais calculadas com as dez observações nos seis perÃodos de -exposição à infestação de Mosca-branca são exibidas. Em todos os casos -observa-se as variâncias amostrais substancialmente menores que -respectivas médias, ainda a manifestação de subdispersão é mais -expressiva na variável número de nós da planta. Portanto, nesse -experimento modelos alternativos ao Poisson devem ser empregados, pois a -suposição de equidispersão é violada. - -<<mv-cottonBolls2, include=FALSE>>= +Na \autoref{fig:descr-cottonBolls2} a disposição de cada uma das +variáveis aleatórias de contagem, \texttt{ncapu}, \texttt{nerep} e +\texttt{nnos}, para os diferentes perÃodos em que as plantas estiveram +sob alta infestação de Mosca-branca é apresentada. Para todas as +variáveis parece haver um comportamento subdisperso, são observadas +muitas contagens sobrepostas e dispostas em um intervalo pequeno de +valores. A indicação de subdispersão também se observa na +\autoref{tab:mv-cottonBolls2}, em que as médias e variâncias amostrais, +calculadas com as dez observações nos seis perÃodos de exposição à +infestação de Mosca-branca, são exibidas. Em todos os casos observa-se +as variâncias amostrais substancialmente menores que respectivas médias, +ainda a manifestação de subdispersão é mais expressiva na variável +número de nós da planta. Portanto, nesse experimento modelos +alternativos ao Poisson devem ser empregados, pois a suposição de +equidispersão é violada. + +\begin{table}[ht] +\centering +\caption{Médias e variâncias amostras das contagens avaliadas no + experimento de capulhos de algodão sob efeito de Mosca-Branca} +\label{tab:mv-cottonBolls2} +\begin{tabular}{>{\centering\arraybackslash} p{2cm}*{6}{c}} + \toprule + \multirow{2}{\linewidth}{Dias de Exposição} & \multicolumn{2}{c}{N. Capulhos} & \multicolumn{2}{c}{N. Estruturas} & \multicolumn{2}{c}{N. Nós} \\ + \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} + & média & variância & média & variância & média & variância \\ + \midrule +<<mv-cottonBolls2, results="asis">>= ##------------------------------------------- ## Calcula as médias e variâncias @@ -190,31 +204,15 @@ mvr <- mvr[order(mvr$va, mvr$dexp), ] ## Organiza lado a lado para apresenta em formato de tabela mv <- mvr[, -(1:2)] -mv <- cbind(mvr[1:6, 1], mv[1:6, ], mv[7:12, ], mv[13:18, ]) +mv <- cbind(mvr[1:6, 1], mv[7:12, ], mv[1:6, ], mv[13:18, ]) -## ## Formata o resultado no ambiente latex e não R (esperamos que eu não -## ## tenha feita nada errada, pois esse procedimento não é reproduzivel) -## xtable(mv) +## Resultados em formato de tabela Latex +print(xtable(mv), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) @ - -\begin{table}[ht] -\centering -\caption{Médias e variâncias amostras das contagens avaliadas no - experimento de capulhos de algodão sob efeito de Mosca-Branca} -\label{tab:mv-cottonBolls2} -\begin{tabular}{>{\centering\arraybackslash} p{2cm}*{6}{c}} - \toprule - \multirow{2}{\linewidth}{Dias de Exposição} & \multicolumn{2}{c}{N. Estruturas} & \multicolumn{2}{c}{N. Capulhos} & \multicolumn{2}{c}{N. Nós} \\ - \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} - & média & variância & média & variância & média & variância \\ - \midrule - 0 & 4,50 & 0,50 & 4,40 & 0,93 & 13,60 & 2,27 \\ - 1 & 4,20 & 1,29 & 3,90 & 1,43 & 16,30 & 0,90 \\ - 2 & 3,90 & 1,21 & 3,40 & 1,60 & 16,10 & 4,54 \\ - 3 & 3,50 & 1,17 & 3,40 & 1,16 & 15,40 & 3,38 \\ - 4 & 3,80 & 1,07 & 3,70 & 1,34 & 15,80 & 2,62 \\ - 5 & 3,80 & 1,07 & 3,80 & 1,07 & 15,70 & 2,68 \\ \bottomrule \end{tabular} \begin{tablenotes} @@ -239,20 +237,21 @@ niveis.K <- paste0( niveis.umid <- paste0( paste(unique(soyaBeans$umid), - collapse = ", "), " \\% do volume total dos poros") + collapse = ", "), "\\% do volume total dos poros") @ -Experimento fatorial 5 $\times$ 3 que estudou diferentes nÃveis de -adubação potássica aplicada ao solo, \Sexpr{niveis.K} e diferentes -nÃveis de umidade do solo, \Sexpr{niveis.umid}, que representam pouca -água, água em quantidade ideal e água em abundância, nos componentes de -produção da soja \cite{Serafim2012}. O experimento foi instalado em casa -de vegetação no delineamento de blocos casualizados completos e a -unidade experimental foi um vaso com duas plantas de soja. No -experimento foram medidas várias variáveis respostas (que representam a -produtividade), sendo que o número de vagens viáveis por vaso e o número -de grãos por vaso foram as variáveis de contagem. +Nesse experimento estudou-se os componentes de produção da soja com +relação à diferentes nÃveis de adubação potássica aplicada ao solo (0, +30, 60, 120 e 180 mg dm$^{-3}$) e diferentes nÃveis de umidade do solo +(37.5, 50, 62.5\%, que representam pouca água, água em quantidade ideal +e água em abundância respectivamente), caracterizando um experimento +fatorial 5 $\times$ 3 \cite{Serafim2012}. O experimento foi instalado em +casa de vegetação no delineamento de blocos casualizados completos e a +unidade experimental foi um vaso com duas plantas de soja. Foram medidas +várias variáveis respostas (que representam a produtividade), sendo que +o número de vagens viáveis por vaso e o número de grãos por vaso foram +as variáveis de contagem. <<descr-soyaBeans, fig.height=4, fig.width=7, fig.cap="Disposição das variáveis número de grãos e número de vagens nos diferentes nÃveis de adubação potássica e umidade do solo.">>= @@ -280,17 +279,19 @@ fonte.xy("Fonte: Elaborado pelo autor.") @ -Na figura \ref{fig:descr-soyaBeans} é apresentada a dispersão das -contagens nas combinações das covariáveis umidade do solo e adubação -potássica. As duas variáveis de contagem avaliadas no experimento -apresentam nÃveis de dispersão distintos, essa caracterÃstica fica -explÃcita na figura \ref{fig:mv-soyaBeans}, em que é exibida a dispersão -entre médias e variâncias amostrais para cada uma das variáveis. Para o -número de grãos por parcela, com contagens mais elevadas, as variâncias -amostrais são, quase em sua totalidade, superiores as médias +Na \autoref{fig:descr-soyaBeans} é apresentada a dispersão das contagens +nas combinações das covariáveis umidade do solo e adubação potássica. As +duas variáveis de contagem avaliadas no experimento apresentam nÃveis de +dispersão distintos. Essa caracterÃstica fica explÃcita na +\autoref{fig:mv-soyaBeans}, em que são exibidas as dispersões entre +médias e variâncias amostrais para cada uma das variáveis, com a linha +pontilhada representando a igualdade entre média e variância +(equidispersão) e a contÃnua um ajuste de regressão linear simples. Para +o número de grãos por parcela, com contagens mais elevadas, as +variâncias amostrais são, quase em sua totalidade, superiores à s médias, caracterizando uma evidência de superdispersão. Já para o número de -vagens por parcela as médias e variâncias são, em média, próximas, o que -indica que a suposição de equidispersão é razoável. +vagens por parcela, as médias e variâncias são próximas, o que indica +que a suposição de equidispersão é razoável. <<mv-soyaBeans, fig.height=3.5, fig.width=7, fig.cap="Médias e variâncias amostrais das contagens de grão e vagens, avaliadas no experimento com soja sob efeito umidade e adubação potássica.">>= @@ -373,19 +374,21 @@ niveis.data <- paste0(format(unique(whiteFly$data), "%d/%m/%y"), Nesse experimento também envolvendo a cultura de soja e a praga Mosca-branca, foram avaliadas plantas de quatro diferentes cultivares de -soja, \Sexpr{niveis.cult}, contabilizando o número de ninfas de -mosca-branca nos folÃolos dos terços superior, médio e inferior das -plantas em seis datas, \Sexpr{niveis.data} dentre os 38 dias de estudo -. O experimento foi conduzido em casa de vegetação sob o delineamento de +soja (BRS 245 RR, BRS 243 RR, BRS 246 RR e BRS 239), contabilizando o +número de ninfas de mosca-branca nos folÃolos dos terços superior, médio +e inferior das plantas em seis datas (11/12/09, 19/12/09, 24/12/09, +02/01/10, 11/01/10 e 18/01/10) dentre os 38 dias de estudo. O +experimento foi conduzido em casa de vegetação sob o delineamento de blocos casualizados para controle de variação local \cite{Suekane2011}. As contagens da praga para cada cultivar em cada uma das datas de avaliação, representadas pelos dias decorridos após a primeira -avaliação, em 11/12/09 são mostradas na figura \ref{fig:descr-ninfas} a -esquerda. As contagens são muito altas e dispersas, principalmente nas -quatro primeiras avaliações. A direita uma descrição no nÃvel de -dispersão da variável de contagem é apresentada. Esse é um conjunto de -dados extremamente superdisperso. Os pontos, que representam as médias e +avaliação, em 11/12/09, são apresentadas à esquerda na +\autoref{fig:descr-ninfas}. As contagens são muito elevadas e dispersas, +principalmente nas quatro primeiras avaliações. À direita da +\autoref{fig:descr-ninfas} uma descrição do nÃvel de dispersão da +variável de contagem é apresentada. Esse é um conjunto de dados +extremamente superdisperso. Os pontos, que representam as médias e variâncias em cada combinação de cultivares de soja e dias após a primeira avaliação, estão todos acima da reta identidade (de equidispersão) com variâncias em torno de 1.000 vezes maiores que as @@ -430,7 +433,7 @@ fonte.xy("Fonte: Elaborado pelo autor") @ -\subsubsection{Peixes Capturados por Pescadores em um Parque Estadual} +\subsubsection{Peixes Capturados por Visitantes em um Parque Estadual} \label{sec:fish} <<data-fish, include=FALSE, echo=FALSE>>= @@ -442,11 +445,11 @@ data(fish, package = "cmpreg") Diferentemente dos demais, esse é um estudo observacional feito por biólogos com interesse em modelar o número de peixes capturados por grupos de pescadores visitantes em um Parque Estadual \cite{UCLA}. Nesse -estudo tem-se como informações a respeito dos grupos de visitantes, o +estudo tem-se como informações referentes ao grupo de visitantes, o número de pessoas e de crianças e se há ou não a presença de -campista. Um fato interessante deste dado é que nem todos os grupos de -visitantes praticaram pescaria, portanto, nesses grupos o número de -peixes capturado será zero. +campista. Um fato interessante nesse estudo é que nem todos os grupos de +visitantes praticaram pescaria, portanto, para esses grupos o número de +peixes capturados será zero. <<descr-fish, fig.height=3.5, fig.width=7.2, fig.cap="LogarÃtmo neperiano do número de peixes capturados acrescido de 0,5 para as diferentes composições dos grupos (esquerda). Histograma do número de peixes capturados por grupo (direita).">>= @@ -485,17 +488,18 @@ fonte.xy("Fonte: Elaborado pelo autor.") @ -Na figura \ref{fig:descr-fish} é evidente o excesso de contagens -zero. No gráfico à esquerda tem-se a disposição das contagens, -transformadas por $\log(y_i|x_i + 0,5)$. É caracterÃstica marcante no -gráfico a grande quantidade de pontos dispostos no primeiro valor do -eixo $y$, \Sexpr{log(0.5)} = $\log(0.5)$. Embora seja um gráfico -marginal, não considerando as covariáveis de cada contagem, a direita um -histograma da variável resposta é realizado e percebe-se novamente a -grande quantidade de valores nulos, ao todo \Sexpr{with(fish, - sum(npeixes == 0)/length(npeixes))*100}\% dos dados são contagens -nulas. Portanto nesse problema, claramente modelos alternativos que -acomodem excesso de zeros se fazem necessários. +Nos gráficos apresentados na \autoref{fig:descr-fish} é evidente o +excesso de contagens zero. No gráfico à esquerda tem-se a disposição das +contagens, transformadas por $\log(y_i + 0,5)$. É caracterÃstica +marcante no gráfico a grande quantidade de pontos dispostos no primeiro +valor do eixo $y$, \Sexpr{log(0.5)} = $\log(0.5)$. À direita da +\autoref{fig:descr-fish} um histograma da variável resposta é +apresentado e, embora seja uma representação da distribuição marginal do +número de peixes capturados (não considera as covariáveis de cada +contagem), percebe-se novamente a grande quantidade de valores nulos, ao +todo \Sexpr{with(fish, sum(npeixes == 0)/length(npeixes))*100}\% dos +dados são contagens nulas. Portanto, nesse problema, modelos +alternativos que acomodem excesso de zeros se fazem necessários. \subsubsection{Número de nematoides em raizes de feijoeiro} \label{sec:nematodes} @@ -513,20 +517,21 @@ key <- list( ## corner = c(0.1, 0.9), type = "b", divide = 1, lines = list(pch = c(NA, 15), lty = c(2, 0), col = cols), - text = list(c("Média de nematoides por linhagem", - "Média de nematoides"))) + text = list(c("Média do nº de nematoides", + "Média do nº de nematoides por linhagem"))) -xyplot(nema ~ cult, data = nematodes, - type = c("p", "g"), +xyplot(nema/off ~ cult, data = nematodes, + type = c("p"), key = key, xlab = "Linhagem de feijoeiro", ylab = "Contagem de nematoides", panel = function(x, y, ...) { + panel.grid() means <- aggregate(y, list(x), mean) - panel.xyplot(x, y, ...) panel.abline(h = mean(y), lty = 2, col = cols[1]) panel.points(x = means[, 1], y = means[, 2], - pch = 15, col = cols[2]) + pch = 15, col = cols[2], cex = 1.1) + panel.xyplot(x, y, ...) }, par.settings = ps.sub) @@ -535,26 +540,26 @@ fonte.xy("Fonte: Elaborado pelo autor.") @ -Esse último conjunto de dados explorado no trabalho, é resultado de um -experimento em casa de vegetação que estudou a reprodução de nematoides -em cultivares de feijoeiro. No experimento, o solo de vasos com duas -plantas de feijão foi inicialmente contaminado com nematoides e as -raizes das duas plantas por vaso foram, ao final do experimento, -lavadas, trituradas, peneiradas e diluÃdas e, a partir de alÃquotas -dessa solução, contou-se o número de nematoides. Como denominador da -contagem tem-se a razão entre a massa fresca de raizes (em gramas) por -parcela e o volume de água (em milÃmetros) utilizado para diluir essa -quantidade \footnote{Cedido para fins acadêmicos por Andressa Cristina - Zamboni Machado, pesquisadora do Instituto Agronômico do Paraná - (IAPAR), e pelo técnico agrÃcola do IAPAR Santino Aleandro da Silva}. +Esse último conjunto de dados explorado no trabalho é resultado de um +experimento em casa de vegetação cujo intersse foi a reprodução de +nematoides em linhagens de feijoeiro. No experimento, o solo de vasos +com duas plantas de feijão foi inicialmente contaminado com nematoides e +as raizes das plantas por vaso foram, ao final do experimento, lavadas, +trituradas, peneiradas e diluÃdas e, a partir de alÃquotas dessa +solução, contou-se o número de nematoides. Como denominador da contagem +tem-se a razão entre a massa fresca de raizes (em gramas) por parcela e +o volume de água (em milÃmetros) utilizado para diluir essa quantidade +\footnote{Cedido para fins acadêmicos por Andressa Cristina Zamboni + Machado, pesquisadora do Instituto Agronômico do Paraná (IAPAR), e + pelo técnico agrÃcola do IAPAR, Santino Aleandro da Silva}. -Na figura \ref{fig:descr-nematodes} a dispersão das contagens de -nematoides em aliquotas da solução de uma grama de massa fresca de raiz +Na \autoref{fig:descr-nematodes} a dispersão das contagens de +nematoides em alÃquotas da solução de uma grama de massa fresca de raiz por um milÃmetro de água para cada linhagem é exibida. As contagens para cada uma das linhagens se distribuem em torno do perfil médio (linha pontilhada). Um detalhe interesse desse conjunto de dados é que o efeito -das linhagens pode ser considerado aleatório em certas fazes do programa +das linhagens pode ser considerado aleatório em certas fases do programa de melhoramento genético. Portanto, pode-se interpretar as linhagens escolhidas como um sorteio aleatório dentre uma população de linhagens de feijoeiro. Assim, modelos com efeitos aleatórios a nÃvel de linhagem @@ -565,25 +570,23 @@ por meio de uma distribuição de probabilidades. \label{cap03:materiais-recursos} O \textit{software} R, versão \Sexpr{with(R.version, paste0(major, ".", - minor))}, é utilizado tanto para a preparação e apresentação dos dados -quanto para ajuste dos modelos e apresentação de resultados. Pacotes -auxiliares utilizados no trabalho são: \texttt{MASS} + minor))}, foi utilizado tanto para a preparação e apresentação dos +dados quanto para ajuste dos modelos e apresentação de +resultados. Pacotes auxiliares utilizados no trabalho são: \texttt{MASS} (\Sexpr{packageVersion("MASS")}) para ajuste e inferências dos modelos -Binomial Negativo, \texttt{bbmle} (versão -\Sexpr{packageVersion("bbmle")}) para estimação via máxima -verossimilhança das funções implementadas para o modelo COM-Poisson , -\texttt{pscl} (\Sexpr{packageVersion("pscl")}) para ajuste dos -modelos Poisson e Binomial Negativo com componente de barreira para -modelagem de excesso de zeros e \texttt{lme4} (versão +Binomial Negativo, \texttt{bbmle} (\Sexpr{packageVersion("bbmle")}) para +estimação via máxima verossimilhança das funções implementadas para o +modelo COM-Poisson , \texttt{pscl} (\Sexpr{packageVersion("pscl")}) para +ajuste dos modelos Poisson e Binomial Negativo com componente de +barreira para modelagem de excesso de zeros e \texttt{lme4} (versão \Sexpr{packageVersion("lme4")}) para ajuste dos modelos Poisson com efeitos aleatórios normais. Para apresentação gráfica dos resultados os pacotes \texttt{lattice} (\Sexpr{packageVersion("lattice")}), \texttt{latticeExtra} (\Sexpr{packageVersion("latticeExtra")}) e \texttt{corrplot} (\Sexpr{packageVersion("corrplot")}) são exaustivamente utilizados. Finalmente, para elaboração do relatório, -mesclando códigos em R e escrita na linguagem de marcação \LaTeX{}, o -pacote \texttt{knitr} (\Sexpr{packageVersion("knitr")}) é -requerido. +mesclando códigos em R e escrita na linguagem de marcação \LaTeX{}, +utilizou-se o pacote \texttt{knitr} (\Sexpr{packageVersion("knitr")}). Destaca-se nesse trabalho que todas as funções implementadas para ajuste e inferência dos modelos de regressão COM-Poisson estão disponÃveis, em @@ -597,107 +600,96 @@ dados exibido no trabalho é ilustrado com códigos R. A estimação dos parâmetros do modelo de regressão COM-Poisson de efeitos fixos é realizada maximizando uma forma reparametrizada da -log-verossimilhança, definida na expressão \ref{eqn:loglik-compoisson}, -via algoritmo numérico de otimização \textit{BFGS}. O parâmetro extra da -COM-Poisson, $\nu$ tem suporte nos reais positivos, restringindo o -espaço paramétrico de busca do otimizador, o que é numericamente -indesejável. Para deixar o domÃnio de busca nos reais reparametrou-se o -modelo com o parâmetro $\phi = \log(\nu)$, como $0 < \nu < \infty$ então -$-\infty < \phi < \infty$. Sob a reparametrização a função a ser -maximizada é - +log-verossimilhança, definida na \autoref{eqn:loglik-compoisson}, via +algoritmo numérico de otimização \textit{BFGS} \cite{Nocedal1995}. O +parâmetro extra da COM-Poisson ($\nu$) tem suporte nos reais positivos, +restringindo o espaço paramétrico de busca do otimizador, o que é +numericamente indesejável. Para deixar o domÃnio de busca nos reais +reparametrizou-se o modelo com $\phi = \log(\nu)$. Como $0 < \nu < +\infty$, então $-\infty < \phi < \infty$. Sob a reparametrização a função +a ser maximizada é dada por \begin{equation} \label{loglik-compoissonr} - \ell(\phi, \beta \mid \underline{y}) = \sum_i^n y_i \log(\lambda_i) - - e^\phi \sum_i^n \log(y!) - \sum_i^n \log(Z(\lambda_i, \phi)) + \ell(\phi,\, \beta \mid \underline{y}) = \sum_{i=1}^n y_i + \log(\lambda_i) - e^\phi \sum_{i=1}^n \log(y!) - \sum_{i=1}^n + \log(Z(\lambda_i, \phi)) \end{equation} - -\noindent -em que $\lambda_i = e^{X_i\beta}$, com $X_i$ o vetor $(x_{i1}, x_{i2}, -\ldots x_{ip})$ de covariáveis da i-ésima observação, e $(\beta, \phi) -\in \R^{p+1}$. - -O ajuste do modelo é realizado sob $\phi$. Portanto as inferências -decorrentes são sobre esse parâmetro. Todavia pode-se retornar para -parametrização original utilizando a função inversa em valores pontuais -ou método delta para funções de $\phi$. Nesse trabalho as inferências -são realizadas sob o parâmetro $\phi$. Para esse parâmetro as -interpretações são como se segue - +em que $\lambda_i = e^{\underline{x}_i^t\beta}$, com $\underline{x}_i$ o +vetor $(x_{i1}, x_{i2}, \ldots x_{ip})$ de covariáveis da i-ésima +observação, e $(\beta, \phi) \in \R^{p+1}$. + +O ajuste do modelo é realizado sob $\phi$. + +As inferências com relação à dispersão, decorrentes do modelo +reparametrizado, são sobre o parâmetro $\phi$. Todavia pode-se retornar +para parametrização original utilizando a função inversa em valores +pontuais ou método delta para funções de $\phi$. Nesse trabalho as +inferências são realizadas sob o parâmetro $\phi$. Para esse parâmetro +as interpretações são como se segue $$ -\phi < 0 \Rightarrow \textrm{Superdispersão} \quad -\phi = 0 \Rightarrow \textrm{Equidispersão} \quad +\phi < 0 \Rightarrow \textrm{Superdispersão}; \quad +\phi = 0 \Rightarrow \textrm{Equidispersão}; \textrm{ e} \quad \phi > 0 \Rightarrow \textrm{Subdispersão} $$ - -\noindent -ou seja, possui a interpretação de um parâmetro de precisão. +ou seja, $\phi$ possui a interpretação de um parâmetro de precisão. A partir dessa reparametrização a condução de testes de hipóteses é -facilitada. Uma vez que $\phi = 0$, representa o caso particular em que +facilitada. Uma vez que $\phi = 0$ representa o caso particular em que a COM-Poisson se reduz a Poisson, a estatÃstica - \begin{equation*} - TRV = 2 \cdot \left ( \ell_{CMP} - \ell_{P} \right ) \sim \rchi^2_{1} + TRV = 2 \left ( \ell_{CMP} - \ell_{P} \right ) \sim \rchi^2_{1} \end{equation*} - -\noindent sendo $\ell_{CMP}$ e $\ell_{P}$ as log-verossimilhanças maximizadas dos modelos COM-Poisson e Poisson com mesmo preditor linear respectivamente, -se refere ao teste de razão de verossimilhanças para $H_0: \phi = 0$, ou -de forma mais apelativa, ao teste sobre a equivalência dos modelos -COM-Poisson e Poisson. +se refere ao teste de razão de verossimilhanças para $H_0: \phi = 0$, +equivalência dos modelos COM-Poisson e Poisson. -A partir da definição em \ref{eqn:pmf-hurdle}, para incluir um +A partir da \autoref{eqn:pmf-hurdle}, para incluir um componente de barreira no modelo COM-Poisson, acomodando excesso de -zeros, adota-se para $\Pr(Z = z \mid \Theta_c)$ a distribuição -COM-Poisson (\ref{eqn:pmf-compoisson}) resultando em - +zeros, adota-se, para $\Pr(Z = z \mid \Theta_c)$, a distribuição +COM-Poisson (\autoref{eqn:pmf-compoisson}), resultando em \begin{equation} \label{eqn:pmf-hurdlecmp} \Pr(Y = y \mid \pi, \phi, \lambda) = \begin{dcases*} - \pi & \text{se } y = 0,\\ + \pi &, \text{se } y = 0\,;\\ (1 - \pi) \frac{\lambda^y}{(y!)^{e^\phi}Z(\lambda, - \phi)}\left (1 - \frac{1}{Z(\lambda, \phi)} \right )^{-1} & + \phi)}\left (1 - \frac{1}{Z(\lambda, \phi)} \right )^{-1} &, \text{se } y = 1, 2, \dots \end{dcases*} \end{equation} - Para modelos de regressão com componente de barreira, são incorporados preditores lineares em $\pi$, -$\underline{\pi}=\frac{\exp(G\gamma)}{1+\exp(G\gamma)}$ e $\lambda$, +$\underline{\pi}=\frac{\exp(Z\gamma)}{1+\exp(Z\gamma)}$ e $\lambda$, $\underline{\lambda}=\exp(X\beta)$ e a verossimilhança desse modelo toma a forma - \begin{equation} \label{eqn:loglik-hurdlecmp} \Ell(\phi, \beta, \gamma \mid \underline{y}) = - \ind [\underline{\pi}] \cdot (1-\ind) \left [ - (1-\underline{\pi})\left ( - \frac{\underline{\lambda}^y}{(y!)^{e^\phi} - Z(\underline{\lambda}, \phi)} - \right ) \left ( - 1-\frac{1}{Z(\underline{\lambda}, \phi)} - \right ) \right ] + \prod_{i \in \Omega_0} \left [ \pi_i \right ] + \prod_{i \in \Omega_+} \left [ + (1-\pi_i) \left ( + \frac{\lambda_i^{y_i}}{(y_i!)^{e^\phi} + Z(\lambda_i, \phi)} + \right ) \left ( 1-\frac{1}{Z(\lambda_i, \phi)} \right ) + \right ] \end{equation} - -\noindent -em que $\ind$ é uma função indicadora para $y = 0$. Os argumentos +sendo $\Omega_0 = \{i \mid y_i = 0\}$ o conjunto de observações que +apresentam contagens 0 e $\Omega_+ = \{i \mid y_i > 0\}$ o conjunto de +observações que apresentam contagens não nulas. Os argumentos $\hat{\phi}$, $\hat{\beta}$ e $\hat{\gamma}$, que maximizam o logaritmo -neperiano da função \ref{eqn:loglik-hurdlecmp} serão as estimativas de +neperiano da \autoref{eqn:loglik-hurdlecmp} serão as estimativas de máxima verossimilhança do modelo COM-Poisson com componente de barreira. Uma outra extensão proposta para o modelo COM-Poisson é a inclusão de efeitos aleatórios a fim de modelar a estrutura experimental ou -observacional de um conjunto de dados. Este trabalho restringe-se a -inclusão de efeitos aleatórios Normais, ou seja, $b \sim +observacional de um conjunto de dados. Neste trabalho restringe-se à +inclusão de efeitos aleatórios Normais, ou seja, $b_j \sim \textrm{Normal}(0, \Sigma)$, que são incorporados sob a forma $\underline{\lambda} = X\beta + Z b$ conforme especificação em \ref{eqn:reg-misto}. Assim, considerando a distribuição COM-Poisson para -a variável resposta condicionada as covariáveis e os efeitos aleatórios, +a variável resposta condicionada à s covariáveis e aos efeitos aleatórios, a verossimilhança pode ser escrita como - \begin{equation} \label{eqn:loglik-mixedcmp} \Ell(\phi, \Sigma, \beta \mid \underline{y}) = @@ -709,30 +701,27 @@ a verossimilhança pode ser escrita como -\frac{1}{2}b^t \Sigma^{-1} b \right ) db_i \end{equation} - -\noindent sendo $m$ o número de grupos que compartilham do mesmo efeito aleatório, $q$ o número de efeitos aleatórios (intercepto aleatório, inclinação e intercepto aleatórios, etc.) e $n_i$ o número de observações no i-ésimo -grupo. A integração em \ref{eqn:loglik-mixedcmp}, necessária para a -avaliação da verossimilhança não tem forma analÃtica. Utiliza-se a +grupo. A integração na \autoref{eqn:loglik-mixedcmp}, necessária para a +avaliação da verossimilhança, não tem forma analÃtica. Utiliza-se a aproximação de Laplace da forma como apresentada em \citeonline[pág. 141]{RibeiroJr2012} para aproximação dessa integral. A -estimação dos parâmetros é realizada via maximização da $\log(\Ell(\phi, -\Sigma, \beta \mid \underline{y}))$ com métodos numéricos de -otimização. Ressalta-se que esse é um procedimento computacionalmente -intensivo, pois a cada iteração do algoritmo de maximização, $m$ -aproximações de Laplace para integrais de dimensão $q$ são -realizadas. Ainda, quando considerada a distribuição COM-Poisson para a -variável resposta condicionalmente independente, tem-se também o cálculo -de $n_m$ constantes normalizadoras $Z(\lambda, \phi)$ -(\ref{eqn:constante-z}) para cada $m$ grupo em cada iteração do -algoritmo de otimização. Com toda essa estrutura hierárquica, -procedimentos computacionais realizados a cada estágio são +estimação dos parâmetros é realizada via maximização da +log-verossimilhança, com métodos numéricos de otimização. Ressalta-se +que esse é um procedimento computacionalmente intensivo, pois a cada +iteração do algoritmo de maximização, $m$ aproximações de Laplace para +integrais de dimensão $q$ são realizadas. Ainda, quando considerada a +distribuição COM-Poisson para a variável resposta condicionalmente +independente, tem-se também o cálculo de $n_m$ constantes normalizadoras +$Z(\lambda, \phi)$ (\autoref{eqn:constante-z}) para cada um dos $m$ grupos +em cada iteração do algoritmo de otimização. Com toda essa estrutura +hierárquica, procedimentos computacionais realizados a cada estágio são potencialmente instáveis numericamente. Para comparação entre os modelos COM-Poisson e demais modelos -listados no capÃtulo \ref{cap:modelos-para-dados-de-contagem} utiliza-se +listados no \autoref{cap:modelos-para-dados-de-contagem} utiliza-se essencialmente o valor maximizado da log-verossimilhança e o critério de informação de Akaike (AIC) definido como @@ -743,16 +732,16 @@ de informação de Akaike (AIC) definido como \noindent sendo $k$ o número de parâmetros e $\ell(\Theta_k, \underline{y})$ a -log\-verossimilhança maximizada do modelo definido pelo conjunto -$\Theta_k$ de parâmetros. Nas análises compara-se também, os nÃveis -descritivos nos testes de razão de verossimilhanças entre modelos +log-verossimilhança maximizada do modelo definido pelo conjunto +$\Theta_k$ de parâmetros. Nas análises compara-se também os nÃveis +descritivos dos testes de razão de verossimilhanças entre modelos encaixados. Nos modelos de regressão de efeitos fixos os valores preditos pelos modelos COM-Poisson e demais alternativas pertinentes são exibidos graficamente com bandas de confiança. -Para maximização numérica das log\-verossimilhanças dos modelos de -regressão COM-Poisson e suas extensões utiliza-se um método de -otimização quasi-Newton bastante popular, denominado \textit{BFGS} -\cite{Nocedal1995}. As informações do vetor gradiente (derivadas de -primeira e matriz hessiana (derivadas de segunda ordem) são obtidos -numericamente via aproximação de diferenças finitas. +Para maximização numérica das log-verossimilhanças dos modelos de +regressão COM-Poisson e suas extensões utiliza-se o método de otimização +quasi-Newton, denominado \textit{BFGS}. O vetor gradiente (derivadas de +primeira ordem) e matriz hessiana (derivadas de segunda ordem) são +obtidos numericamente via aproximação de diferenças finitas +\cite{Nocedal1995}. diff --git a/docs/cap04_resultados-e-discussao.Rnw b/docs/cap04_resultados-e-discussao.Rnw index 7e9bcc834faec20496b049ab3ab76c1b20b19e99..152d68d1df243f47551831b9cb2364b071f0e233 100644 --- a/docs/cap04_resultados-e-discussao.Rnw +++ b/docs/cap04_resultados-e-discussao.Rnw @@ -2,14 +2,14 @@ % CAPÃTULO 4 - RESULTADOS E DISCUSSÃO % ------------------------------------------------------------------------ -Nesse capÃtulo são apresentados os resultados e discussões da aplicação -modelos de regressão COM-Poisson ajustados aos dados apresentados na -seção \ref{cap03:materiais-dados} comparando-as com abordagens já -utilizadas na EstatÃstica aplicada. As primeiras seis seções são -destinadas a apresentação das análises estatÃsticas de cada conjunto de -dados citado. Na seção \ref{cap04:discussao} discussões gerais sobre os -resultados dos modelos COM-Poisson empregados nas análises são -realizadas. +Neste capÃtulo são apresentados os resultados e discussões da aplicação +dos modelos de regressão COM-Poisson ajustados aos dados apresentados na +\autoref{cap03:materiais-dados}. Os resultados são comparados com +abordagens já utilizadas na EstatÃstica aplicada. As primeiras seis +seções são destinadas à apresentação das análises estatÃsticas de cada +conjunto de dados citado. Na \autoref{cap04:discussao} discussões +gerais sobre os resultados dos modelos COM-Poisson empregados nas +análises são realizadas. \section{Análise de dados de capulhos de algodão sob efeito de desfolha} \label{sec:analise-cottonBolls} @@ -54,11 +54,12 @@ prof.cottonBolls <- profile(m5C, which = "phi") @ -Diante da estrutura do experimento apresentada na seção -\ref{sec:cottonBolls} foram propostos, por \citeonline{Zeviani2014}, +Diante da estrutura do experimento apresentada na +\autoref{sec:cottonBolls} foram propostos, por \citeonline{Zeviani2014}, cinco preditores crescentes em complexidade que testam aspectos -interesses sobre os fatores experimentais. Abaixo os cinco -preditores considerados são descritos. +interesses sobre os fatores experimentais. Abaixo os cinco preditores +considerados são descritos, sendo \texttt{def} a covariável que +representa o nÃvel de desfolha artificial (0, 25, 50, 75 e 100\%). \noindent Preditor 1: $g(\mu) = \beta_0$ \\ @@ -73,27 +74,27 @@ Preditor 5: $g(\mu) = \beta_0 + \beta_{1j} \textrm{def} + \beta_{2j} \noindent onde $j$ varia nos nÃveis de estágio fenológico da planta (1: vegetativo, 2: botão floral, 3: florescimento, 4: maça, 5: capulho) e -$g(\mu)$ uma função de ligação. A proposta desses preditores foi -realizada de forma aninhada a fim de facilitar a condução de testes de -hipóteses. O modelo 1 contêm somente o intercepto, e é ajustado apenas -como ponto de partida para verificar como modelos mais estruturados -melhoram o ajuste. O modelo 2 apresenta apenas o efeito de desfolha de -forma linear, o modelo 3 é o modelo 2 somado um efeito de segunda -ordem. O modelo 4, apresenta o efeito de desfolha linear mudando de -acordo com o estágio de crescimento (interação entre o efeito linear de -desfolha e estágio), e por fim o modelo 5 não somente o efeito de -primeira ordem muda com o estágio de crescimento, mais também o efeito +$g(\mu)$ a função de ligação considerada no modelo. A proposta desses +preditores foi realizada de forma aninhada a fim de facilitar a condução +de testes de hipóteses. O modelo 1 contêm somente o intercepto, e é +ajustado apenas como ponto de partida para verificar como modelos mais +estruturados melhoram o ajuste. O modelo 2 apresenta apenas o efeito de +desfolha de forma linear. O modelo 3 é o modelo 2 somado um efeito de +segunda ordem. O modelo 4, apresenta o efeito de desfolha linear mudando +de acordo com o estágio de crescimento (interação entre o efeito linear +de desfolha e estágio). E por fim, no modelo 5 não somente o efeito de +primeira ordem muda com o estágio de crescimento, mas também o efeito de segunda ordem (interação entre o efeito de primeira e segunda ordem de desfolha e estágio). -A seguir são ajustados os modelos Poisson e COM-Poisson como -alternativas paramétricas à análise de dados e como alternativa -semi-paramétrica a estimação via quasi-verossimilhança Poisson. Na -tabela \ref{tab:ajuste-cottonBolls} os resultados dos três modelos -ajustados aos cinco preditores são apresentados. O modelo COM-Poisson -apresentou melhor ajuste dentre todos os preditores considerados quando -comparado ao Poisson, indicado pelas maiores log-verossimilhanças e -menores AIC's. +Na sequência da análise, foram ajustados os modelos Poisson e +COM-Poisson como alternativas paramétricas à análise de dados e, como +alternativa semi-paramétrica, a estimação via quasi-verossimilhança +Poisson. Na \autoref{tab:ajuste-cottonBolls} os resultados dos três +modelos ajustados aos cinco preditores são apresentados. O modelo +COM-Poisson apresentou melhor ajuste dentre todos os preditores +considerados quando comparado ao Poisson, indicado pelas maiores +log-verossimilhanças e menores AIC's. <<loglik-cottonBolls, include=FALSE>>= @@ -118,10 +119,9 @@ tab.ajuste <- rbind(tabP, tabC, tabQ) rownames(tab.ajuste) <- NULL tab.ajuste <- data.frame(etas, tab.ajuste) -## ##---------------------------------------------------------------------- -## ## Copiar e colar o corpo do resultado na customização latex abaixo -## digits <- c(1, 0, 0, 2, 2, 2, 0, -2, 3, -2) -## xtable(tab.ajuste, digits = digits) +## Para notação cientÃfica nas tabelas Latex +digits <- c(1, 0, 0, 2, 2, 2, 0, -2, 3, -2) + @ \begin{table}[ht] @@ -134,31 +134,45 @@ tab.ajuste <- data.frame(etas, tab.ajuste) \toprule Poisson & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & & \\ \midrule - Preditor 1 & 1 & -279,93 & 561,87 & & & & & \\ - Preditor 2 & 2 & -272,00 & 548,00 & 15,86 & 1 & 6,81E-05 & & \\ - Preditor 3 & 3 & -271,35 & 548,71 & 1,29 & 1 & 2,56E-01 & & \\ - Preditor 4 & 7 & -258,67 & 531,35 & 25,36 & 4 & 4,26E-05 & & \\ - Preditor 5 & 11 & -255,80 & 533,61 & 5,74 & 4 & 2,19E-01 & & \\[0.3cm] - COM-Poisson & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & $\hat{\phi}$ & P($>\rchi^2$) \\ +<<tab-cottonBolls1, results="asis">>= + +## Resultados em formato de tabela Latex +print(xtable(tab.ajuste[1:5, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento + COM-Poisson & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & $\hat{\phi}$ & P($>\rchi^2$) \\ \midrule - Preditor 1 & 2 & -272,48 & 548,96 & & & & 0,551 & 1,13E-04 \\ - Preditor 2 & 3 & -257,46 & 520,93 & 30,03 & 1 & 4,25E-08 & 0,794 & 6,97E-08 \\ - Preditor 3 & 4 & -256,09 & 520,18 & 2,75 & 1 & 9,73E-02 & 0,816 & 3,29E-08 \\ - Preditor 4 & 8 & -220,20 & 456,40 & 71,78 & 4 & 9,54E-15 & 1,392 & 1,75E-18 \\ - Preditor 5 & 12 & -208,25 & 440,50 & 23,90 & 4 & 8,38E-05 & 1,585 & 1,80E-22 \\[0.3cm] +<<tab-cottonBolls2, results="asis">>= + +## Resultados em formato de tabela Latex +print(xtable(tab.ajuste[6:10, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento Quase-Poisson & np & deviance & AIC & F & diff np & P($>F$) & $\hat{\sigma}^2$ & P($>\rchi^2$) \\ \midrule - Preditor 1 & 1 & 75,51 & & & & & 0,567 & 3,66E-04 \\ - Preditor 2 & 2 & 59,65 & & 34,21 & 1 & 4,17E-08 & 0,464 & 5,13E-07 \\ - Preditor 3 & 3 & 58,36 & & 2,81 & 1 & 9,62E-02 & 0,460 & 3,66E-07 \\ - Preditor 4 & 7 & 33,00 & & 22,77 & 4 & 5,89E-14 & 0,278 & 9,15E-16 \\ - Preditor 5 & 11 & 27,25 & & 5,96 & 4 & 2,18E-04 & 0,241 & 3,57E-18 \\ +<<tab-cottonBolls3, results="asis">>= + +## Resultados em formato de tabela Latex +print(xtable(tab.ajuste[11:15, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ \bottomrule \end{tabular} \begin{tablenotes} \footnotesize -\item np, número de parâmetros, diff $\ell$, diferença entre - log-verossimilhanças, F, estatÃstica F baseada nas quasi-deviances, +\item np, número de parâmetros; diff $\ell$, diferença entre + log-verossimilhanças; F, estatÃstica F baseada nas quasi-deviances; diff np, diferença entre o np. \\[0.1cm] \item Fonte: Elaborado pelo autor. \end{tablenotes} @@ -166,56 +180,20 @@ tab.ajuste <- data.frame(etas, tab.ajuste) As estimativas dos parâmetros extras $\phi$ e $\sigma^2$ dos modelos COM-Poisson e Quasi-Poisson respectivamente, também são apresentadas na -tabela \ref{tab:ajuste-cottonBolls} e indicam subdispersão ($\phi>0$ e -$\sigma^2<1$). Note que, mesmo não considerando covariáveis, preditor 1, -a hipótese de equidispersão foi rejeitada pelo modelos COM-Poisson e +\autoref{tab:ajuste-cottonBolls} e indicam subdispersão ($\phi>0$ e +$\sigma^2<1$). Note que, mesmo não considerando covariáveis (preditor 1) +a hipótese de equidispersão foi rejeitada pelos modelos COM-Poisson e Quasi-Poisson. Isso se reflete nos nÃveis descritivos dos testes de razão de verossimilhanças realizados, em que o modelo Poisson, em discordância com os demais, não indicou significância do efeito -quadrático por nÃvel de desfolha, preditor 5, pois superestima a -variabilidade do processo. Esses resultados estão de acordos com os -apresentados por \citeonline{Zeviani2014}, onde um modelo -\textit{Gamma-Count} foi ajustado, destaca-se a similaridade entre as +quadrático do nÃvel de desfolha por estágio fenológico (preditor 5), +pois superestima a variabilidade do processo. Esses resultados estão de +acordos com os apresentados por \citeonline{Zeviani2014}, onde um modelo +\textit{Gamma-Count} foi ajustado. Destaca-se a similaridade entre as medidas de ajuste dos modelos COM-Poisson e \textit{Gamma-Count}. Os valores das log-verossimilhanças maximizadas nos dois modelos difere somente nas casas decimais, para todos os preditores. -<<prof-cottonBolls, fig.height=4, fig.width=4.5, out.width="0.6\\linewidth", fig.cap="Perfil de log-verossimilhança para o parâmetro extra da COM-Poisson, estimado no modelo com o quinto preditor.">>= - -myprof(prof.cottonBolls, namestrip = expression("Perfil para"~phi), - par.settings = ps.sub) -fonte.xy("Fonte: Elaborado pelo autor.") - -@ - -Na figura \ref{fig:prof-cottonBolls} a avaliação do parâmetro $\phi$ do -modelo COM-Poisson com efeito de desfolha artificial de primeira e -segunda ordem para cada estágio fenológico, via verossimilhança -perfilhada é apresentada. O valor zero, que representa a não necessidade -de um modelo COM-Poisson está dentro dos limites de confiança de 90, 95 -e até 99\%. A simetria do perfil de verossimilhança também é algo para -se destacar, pois neste caso intervalos do tipo Wald (computacionalmente -mais fáceis), via aproximação quadrática da verossimilhança, podem ser -construÃdos, muito embora os construÃdos via perfil de -log-verossimilhança sejam preferÃveis. Em concordância com a figura, o -teste de hipóteses via razão de verossimilhanças para $H_0: \phi = 0$, -rejeitou a hipótese nula com um nÃvel de significância muito próximo a -zero, tabela \ref{tab:ajuste-cottonBolls}. - -<<coef-cottonBolls, include=FALSE>>= - -pnames <- c("phi", "beta0", paste0("beta1", 1:5), - paste0("beta2", 1:5)) - -## Tabela com os coeficientes -tab.coef <- coeftab(m5P, m5Q, m5C, - rownames = paste0("$++", pnames, "$")) - -## ## Código latex para tabela com os coeficientes -## print.xtable(xtable(tab.coef), include.rownames = TRUE) - -@ - \begin{table}[ht] \centering \caption{Estimativas dos parâmetros e razões entre as estimativa e erro @@ -227,18 +205,24 @@ tab.coef <- coeftab(m5P, m5Q, m5C, \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} Parâmetro & Estimativa & Est/EP & Estimativa & Est/EP & Estimativa & Est/EP \\ \midrule - $\sigma^2,\,\phi$ & -- & -- & 0,24 & -- & 1,58 & 12,42 \\ - $\beta_{0}$ & 2,19 & 34,57 & 2,19 & 70,42 & 10,90 & 7,76 \\ - $\beta_{11}$ & 0,44 & 0,85 & 0,44 & 1,73 & 2,02 & 1,77 \\ - $\beta_{12}$ & 0,29 & 0,57 & 0,29 & 1,16 & 1,34 & 1,21 \\ - $\beta_{13}$ & -1,24 & -2,06 & -1,24 & -4,19 & -5,75 & -3,89 \\ - $\beta_{14}$ & 0,36 & 0,64 & 0,36 & 1,31 & 1,60 & 1,30 \\ - $\beta_{15}$ & 0,01 & 0,02 & 0,01 & 0,04 & 0,04 & 0,03 \\ - $\beta_{21}$ & -0,81 & -1,38 & -0,81 & -2,81 & -3,72 & -2,78 \\ - $\beta_{22}$ & -0,49 & -0,86 & -0,49 & -1,75 & -2,26 & -1,80 \\ - $\beta_{23}$ & 0,67 & 0,99 & 0,67 & 2,01 & 3,13 & 2,08 \\ - $\beta_{24}$ & -1,31 & -1,95 & -1,31 & -3,97 & -5,89 & -3,66 \\ - $\beta_{25}$ & -0,02 & -0,04 & -0,02 & -0,07 & -0,09 & -0,08 \\ +<<coef-cottonBolls, results="asis">>= + +pnames <- c("\\sigma^2,\\,\\phi ", "\\beta_0", + paste0("\\beta_{1", 1:5, "}"), + paste0("\\beta_{2", 1:5, "}")) + +## Tabela com os coeficientes +tab.coef <- coeftab(m5P, m5Q, m5C, + rownames = paste0("$", pnames, "$")) + +print.xtable(xtable(tab.coef), + include.rownames = TRUE, + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE, + sanitize.text.function = function(x) x) + +@ \bottomrule \end{tabular} \begin{tablenotes} @@ -247,20 +231,27 @@ tab.coef <- coeftab(m5P, m5Q, m5C, \end{tablenotes} \end{table} -As estimativas dos efeitos lineares e quadráticos de desfolha -artificial, conforme notação do preditor 5, são apresentadas na tabela -\ref{tab:coef-cottonBolls} para os modelos Poisson, Quasi-Poisson e -COM-Poisson. Para os modelos Poisson e Quasi-Poisson as estimativas são -idênticas, por construção \ref{cap02:poisson}, o que difere são as -magnitudes dessas estimativas em comparação com seu erro padrão, que no -caso Quasi-Poisson é corrigido pelo parâmetro $\sigma^2$. Considerando o -modelo COM-Poisson as estimativas são notavelmente diferentes, pois o -preditor linear é construÃdo em $\lambda$, da expressão -\ref{eqn:pmf-compoisson}, e este parâmetro não descreve, diretamente, a -média da distribuição. Sendo assim as estimativas do COM-Poisson não -podem ser comparadas com as demais estimativas. Contudo, a magnitude -desses efeitos com relação ao efeito padrão sim. E neste caso os modelos -Quasi-Poisson e COM-Poisson levam as mesmas conclusões. +<<prof-cottonBolls, fig.height=3.5, fig.width=4, out.width="0.6\\linewidth", fig.cap="Perfil de log-verossimilhança para o parâmetro extra da COM-Poisson, estimado no modelo com o efeito quadrático do nÃvel de desfolha por cada estágio fenológico.">>= + +myprof(prof.cottonBolls, namestrip = expression("Perfil para"~phi), + par.settings = ps.sub) +fonte.xy("Fonte: Elaborado pelo autor.") + +@ + +Na \autoref{fig:prof-cottonBolls} a avaliação do parâmetro $\phi$ do +modelo COM-Poisson com efeito de desfolha artificial de primeira e +segunda ordem para cada estágio fenológico, via verossimilhança +perfilhada, é apresentada. O valor zero, que representa a não +necessidade de um modelo COM-Poisson, não está dentro dos limites de +confiança de 99, 95 e até 90\%. A simetria do perfil de verossimilhança +também é algo para se destacar, pois neste caso intervalos do tipo Wald +(computacionalmente mais fáceis), via aproximação quadrática da +verossimilhança, podem ser construÃdos, muito embora os construÃdos via +perfil de log-verossimilhança sejam preferÃveis. Em concordância com a +figura, o teste de hipóteses via razão de verossimilhanças para $H_0: +\phi = 0$ (última coluna da \autoref{tab:ajuste-cottonBolls}), rejeitou +a hipótese nula com um nÃvel de significância muito próximo a zero. <<corr-cottonBolls, fig.width=7, fig.height=7, fig.cap="Imagem da matriz de correlação entre os parâmetros do modelo COM-Poisson.">>= @@ -275,26 +266,38 @@ fonte("Fonte: Elaborado pelo autor.") @ -As covariâncias entre as estimativas dos parâmetros do modelo -COM-Poisson são apresentadas, na escala da correlação, na figura -\ref{fig:corr-cottonBolls}. Destaca-se nessa figura a forte correlação -do parâmetro de precisão $\phi$ com os $\beta$'s da regressão. Embora -seja uma representação empÃrica, observada a esse particular conjunto de -dados, nota-se a não ortogonalidade na matriz de informação observada, o -que implica que inferências sobre os $\beta$'s são condicionais a -$\phi$. Esse comportamento dos modelos COM-Poisson é recorrente, como -será visto nos demais conjuntos de dados. +As estimativas dos efeitos lineares e quadráticos de desfolha +artificial, conforme notação do preditor 5, são apresentadas na +\autoref{tab:coef-cottonBolls} para os modelos Poisson, Quasi-Poisson e +COM-Poisson. Para os modelos Poisson e Quasi-Poisson as estimativas são +idênticas, por construção (veja \autoref{cap02:poisson}), o que difere +são as magnitudes dessas estimativas em comparação com seu erro padrão, +que no caso Quasi-Poisson é corrigido pelo parâmetro +$\sigma^2$. Considerando o modelo COM-Poisson as estimativas são +notavelmente diferentes, pois o preditor linear é construÃdo em +$\lambda$, da \autoref{eqn:pmf-compoisson}, e esse parâmetro não +descreve, diretamente, a média da distribuição. Sendo assim as +estimativas do COM-Poisson não podem ser comparadas com as demais +estimativas. Contudo a magnitude desses efeitos, com relação ao seu erro +padrão, sim. E nesse caso, os modelos Quasi-Poisson e COM-Poisson levam +as mesmas conclusões. + +Devido ao modelo COM-Poisson não ser construÃdo diretamente para a +média, as estimativas dos parâmetros não refletem efeitos +multiplicativos, como ocorre nos casos Poisson e Quasi-Poisson. Com +isso, a interpretação dos efeitos nesse modelo é somente com relação ao +sinal da estimativa, quando positivo indica um aumento na média da +variável de interesse, e quando negativo uma diminuição. -Essa caracterÃstica de não ortogonalidade da matriz de informação -observada teve de ser levada em consideração para cálculo dos valores -preditos, uma vez que a informação sobre a incerteza das estimativas -contida na matriz de variâncias e covariâncias não pôde ser -marginalizada para os $\beta$'s, que efetivamente são utilizados para -cálculo de $\hat{\lambda}_i$ e consequentemente $\hat{\mu}_i$. Portanto, -para cálculo dos valores preditos utiliza-se a matriz de variâncias e -covariâncias condicionada a $\phi$, conforme \citeonline[teorema 3.6, -pág. 123]{Ferreira2011}. Essa é uma prática tomada também para cálculo -dos valores preditos nos demais conjunto de dados. +As covariâncias entre as estimativas dos parâmetros do modelo +COM-Poisson são apresentadas, na escala da correlação, na +\autoref{fig:corr-cottonBolls}. Destaca-se nessa figura a forte +correlação do parâmetro de precisão $\phi$ com os $\beta$'s da +regressão. Embora seja uma representação empÃrica, observada a esse +particular conjunto de dados, nota-se a não ortogonalidade na matriz de +informação observada, o que implica que inferências sobre os $\beta$'s +são condicionais a $\phi$. Esse comportamento dos modelos COM-Poisson é +recorrente, como será visto nos demais conjuntos de dados. <<pred-cottonBolls, fig.height=3.4, fig.width=6.7, fig.cap="Curva dos valores preditos com intervalo de confiança de (95\\%) como função do nÃvel de desfolha e do estágio fenológico da planta.">>= @@ -368,14 +371,30 @@ fonte.xy("Fonte: Elaborado pelo autor.") @ +Essa caracterÃstica de não ortogonalidade da matriz de informação +observada teve de ser levada em consideração para cálculo dos valores +preditos, uma vez que a informação sobre a incerteza das estimativas +contida na matriz de variâncias e covariâncias não pôde ser +marginalizada para os $\beta$'s, que efetivamente são utilizados para +cálculo de $\hat{\lambda}_i$ e consequentemente $\hat{\mu}_i$. Portanto, +no cálculo dos valores preditos utiliza-se a matriz de variâncias e +covariâncias condicional a $\hat{\phi}$, conforme teorema 3.6 +\citeonline[pág. 123]{Ferreira2011}. Para computação dos intervalos de +confiança utiliza-se o método delta \cite{Ribeiro2012}. A utilização da +matriz de variâncias e covariâncias condicional e o método delta para +computação dos valores preditos, são práticas tomadas também na análise +dos demais conjuntos de dados. + As médias com intervalos de confiança calculadas com os modelos -COM-Poisson e Quasi-Poisson são praticamente idênticas, conforme pode -ser visto na figura \ref{fig:pred-cottonBolls}. Contudo, destaca-se que -o modelo COM-Poisson é totalmente paramétrico permitindo representar uma +COM-Poisson e Quasi-Poisson são idênticas, conforme pode ser visto na +\autoref{fig:pred-cottonBolls}. Isso se deve ao fato da relação média e +variância ser aproximada de forma satisfatória por $\frac{1}{\nu}E(Y)$ +nos casos de subdispersão, no modelo COM-Poisson (vide +\autoref{fig:mv-compoisson}). Contudo, destaca-se que o modelo +COM-Poisson é totalmente paramétrico permitindo representar uma distribuição, calculando probabilidades, o que não é possÃvel com a -formulação Quasi-Poisson. Ainda nota-se claramente que o modelo Poisson -é inadequado a esse conjunto de dados e que inferências a partir deste -seriam incorretas. +formulação Quasi-Poisson. Como visto o modelo Poisson é inadequado a +esse conjunto de dados e inferências a partir deste são incorretas. \section{Análise de dados de capulhos de algodão sob efeito de Mosca-Branca} \label{sec:analise-cottonBolls2} @@ -448,8 +467,8 @@ prof.nnos <- profile(m3C.nnos, which = "phi") Nesse conjunto de dados também há indÃcios de subdispersão para as três variáveis de interesse mensuradas no estudo, conforme apresentado -na seção \ref{sec:cottonBolls2}. Para cada contagem procedeu-se com o -ajuste dos modelos Poisson, Quasi-Poisson e COM-Poisson adotando os +na \autoref{sec:cottonBolls2}. Para cada contagem procedeu-se com o +ajuste dos modelos Poisson, Quasi-Poisson e COM-Poisson com os preditores: \noindent @@ -501,8 +520,7 @@ phis <- cmptest(m3C.ncapu, m2C.nerep, m3C.nnos) ## ##---------------------------------------------------------------------- ## ## Copiar e colar o corpo do resultado na customização latex abaixo -## digits <- c(1, 1, 0, 2, 2, -2, 2, 2, -2, 2, -2) -## print(xtable(tab.ajuste, digits = digits)) +digits <- c(1, 1, 0, 2, 2, 4, 2, 2, 4, 2, 4) @ @@ -525,19 +543,36 @@ linear e quadrático dos dias de exposição, respectivamente. \midrule \multicolumn{10}{l}{{\bfseries \footnotesize Número de capulhos produzidos}} \\[0.0cm] \cline{1-5} \\[-0.2cm] - & 1 & -105,27 & 212,55 & & -92,05 & 188,09 & & 20,80 & \\ - & 2 & -105,03 & 214,05 & 4,83E-01 & -91,31 & 188,62 & 2,25E-01 & 20,31 & 2,23E-01 \\ - & 3 & -104,44 & 214,88 & 2,78E-01 & -89,47 & 186,95 & 5,52E-02 & 19,13 & 6,16E-02 \\[0.2cm] +<<tab-conttonBolls21, results="asis">>= + +print(xtable(tab.ajuste[1:3, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ + \specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento \multicolumn{10}{l}{{\bfseries \footnotesize Número de estruturas reprodutivas}} \\[0.0cm] \cline{1-5} \\[-0.2cm] - & 1 & -104,74 & 211,49 & & -86,41 & 176,82 & & 16,23 & \\ - & 2 & -104,27 & 212,54 & 3,32E-01 & -84,59 & 175,18 & 5,66E-02 & 15,29 & 6,19E-02 \\ - & 3 & -104,06 & 214,12 & 5,16E-01 & -83,73 & 175,47 & 1,90E-01 & 14,87 & 2,07E-01 \\[0.2cm] +<<tab-conttonBolls22, results="asis">>= + +print(xtable(tab.ajuste[4:6, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ + \specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento \multicolumn{10}{l}{{\bfseries \footnotesize Número de nós da planta}} \\[0.0cm] \cline{1-5} \\[-0.2cm] - & 1 & -143,79 & 289,59 & & -120,58 & 245,16 & & 12,69 & \\ - & 2 & -143,48 & 290,95 & 4,25E-01 & -119,03 & 244,06 & 7,87E-02 & 12,05 & 7,39E-02 \\ - & 3 & -142,95 & 291,89 & 3,04E-01 & -116,27 & 240,54 & 1,88E-02 & 11,00 & 2,23E-02 \\ +<<tab-conttonBolls23, results="asis">>= + +print(xtable(tab.ajuste[7:9, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ \bottomrule \end{tabular} \begin{tablenotes} @@ -547,33 +582,32 @@ linear e quadrático dos dias de exposição, respectivamente. \end{tablenotes} \end{table} -Na tabela \ref{tab:ajuste-cottonBolls2} são exibidas as medidas de -ajuste dos modelos para as três variáveis resposta. Em todos os casos o -modelo COM-Poisson apresentou maiores log-verossimilhanças indicando um -melhor ajuste, quando comparado ao Poisson, também indicado pelos os -valores de AIC que ponderam a log-verossimilhança pelo número de -parâmetros considerados no modelo. Para questões inferenciais novamente, -há um desacordo entre os modelos paramétricos. Pelos modelos Poisson -não há evidências para manutenção de nenhum efeito da variável número -de dias sob infestação, em todos os casos, ao passo que no modelo -COM-Poisson tem-se evidências do efeito quadrático quando considerado o -modelo para o número de nós da planta (nÃvel descritivo de -\Sexpr{1-anC.nnos[3,6]}) e o número de capulhos produzidos (nÃvel -descritivo de \Sexpr{1-anC.ncapu[3,6]}, na borda da região de -significância, mas com uma diminuição do AIC em favor do efeito -quadrático). Quando modelado o número de estruturas reprodutivas o -modelo COM-Poisson também não indicou efeito quadrático, contudo o -efeito linear de \texttt{dexp} pode ser discutido uma vez que a -significância do TRV foi de \Sexpr{anC.ncapu[3,6]} e o AIC apresentou um -pequeno aumento com relação ao modelo nulo. Considera-se nas demais -inferências os preditores com efeitos linear, para o número de -estruturas reprodutivas e quadrático, para o número de capulhos -produzidos e número de nós da planta. - -A especificação do modelo via Quasi-Verossimilhança Poisson obteve +Na \autoref{tab:ajuste-cottonBolls2} são exibidas as medidas de ajuste +dos modelos para as três variáveis resposta. Em todos os casos o modelo +COM-Poisson apresentou maiores log-verossimilhanças indicando um melhor +ajuste, quando comparado ao Poisson, também indicado pelos valores de +AIC que ponderam a log-verossimilhança pelo número de parâmetros +considerados no modelo. Para questões inferenciais, novamente, há um +desacordo entre os modelos paramétricos. Pelo modelo Poisson não há +evidências para manutenção de nenhum efeito da variável número de dias +sob infestação, em todos os casos, ao passo que no modelo COM-Poisson +tem-se evidências do efeito quadrático quando considerado o modelo para +o número de nós da planta (nÃvel descritivo de \Sexpr{anC.nnos[3,6]}) e +o número de capulhos produzidos (nÃvel descritivo de +\Sexpr{anC.ncapu[3,6]}, na borda da região de significância, mas com uma +diminuição do AIC em favor do efeito quadrático). Quando modelado o +número de estruturas reprodutivas, o modelo COM-Poisson também não +indicou efeito quadrático, contudo o efeito linear de \texttt{dexp} pode +ser discutido uma vez que a significância do TRV foi de +\Sexpr{anC.ncapu[3,6]} e o AIC apresentou um pequeno aumento com relação +ao modelo nulo. Considera-se nas demais inferências os preditores com +efeito linear, para o número de estruturas reprodutivas e quadrático, +para o número de capulhos produzidos e número de nós da planta. + +Na estimação dos parâmetros via quasi-Verossimilhança Poisson obteve-se nÃveis descritivos mais conservadores para a rejeição da hipótese nula -que o modelo COM-Poisson. Contudo, para escolha de preditores as mesmas -tendências apontadas pelo COM-Poisson foram seguidas. +que no modelo COM-Poisson. Contudo, para escolha de preditores os +resultados se mostram equivalentes. <<prof-cottonBolls2, fig.height=3, fig.width=7, fig.cap="Perfis de log-verossimilhança para o parâmetro extra da COM-Poisson nos modelos para número de capulhos produzidos (esquerda), número de estruturas reprodutivas (central) e número de nós (direira).">>= @@ -634,22 +668,22 @@ xyplot(abs(z) ~ focal | param, data = da, @ Para avaliação do parâmetro $\phi$ da COM-Poisson nos três modelos -considerados, intervalos de confiança construÃdos sob -perfilhamento da verossimilhança são exibidos na figura -\ref{fig:prof-cottonBolls2}. Para nenhum dos modelos o valor de $\phi = -0$ esteve dentro dos limites de confiança de 90, 95 e 99\%. Os valores -estimados dos parâmetros nos modelos para número de capulhos, número de -estruturas reprodutivas e número de nós da planta foram de \Sexpr{phis[, - 1]} respectivamente, indicando subdispersão em todos os casos. - -Na figura \ref{fig:corr-cottonBolls2} são representadas as matrizes de +considerados, intervalos de confiança construÃdos sob perfilhamento da +verossimilhança são exibidos na \autoref{fig:prof-cottonBolls2}. Nenhum +dos intervalos, de 99, 95 e 90\% de confiança, compreende o valor zero +para $\phi$. Os valores estimados dos parâmetros nos modelos para número +de capulhos, número de estruturas reprodutivas e número de nós da planta +foram de \Sexpr{phis[, 1]} respectivamente, indicando subdispersão em +todos os casos. + +Na \autoref{fig:corr-cottonBolls2} são representadas as matrizes de covariâncias (via correlações) entre as estimativas dos modelos para -número de capulhos, à esquerda, número de estruturas reprodutivas, ao -centro e número de nós da plantas, à direita. A forte correlação entre o -parâmetro de precisão $\phi$ e $\beta_0$ (principalmente) também foi +número de capulhos (à esquerda), número de estruturas reprodutivas (ao +centro) e número de nós da plantas (à direita). A forte correlação entre +o parâmetro de precisão $\phi$ e $\beta_0$ (principalmente) também foi observada no ajuste do modelo para esses conjuntos de dados. -<<corr-cottonBolls2, fig.height=1.7, fig.width=5, out.width="1\\linewidth", fig.cap="Imagem da matriz de correlação entre os parâmetros do modelo COM-Poisson. (Esquerda) Modelo para o número de capulhos por parcela, (centro) para o número de estruturas reprodutivas e (direita) para o número de nós por parcela.">>= +<<corr-cottonBolls2, fig.height=1.7, fig.width=5, out.width="1\\linewidth", fig.cap="Imagem da matriz de correlação entre os parâmetros do modelo COM-Poisson. (esquerda) Modelo para o número de capulhos por parcela, (centro) para o número de estruturas reprodutivas e (direita) para o número de nós por parcela.">>= pnames <- c("phi", "beta0", "beta1", "beta2") @@ -667,9 +701,12 @@ Corr.nerep <- cov2cor(Vcov) dimnames(Corr.nerep) <- list(pnames, pnames) par(mfrow = c(1, 3)) -mycorrplot(Corr.ncapu, mar = c(1.5, 0, 0, 0)) -mycorrplot(Corr.nerep, mar = c(1.5, 0, 0, 0)) -mycorrplot(Corr.nnos, mar = c(1.5, 0, 0, 0)) +mycorrplot(Corr.ncapu, mar = c(2.5, 0, 0, 0)) +mtext(text = "Capulhos produzidos", side = 1, cex = 0.7, adj = 0.65) +mycorrplot(Corr.nerep, mar = c(2.5, 0, 0, 0)) +mtext(text = "Estruturas reprodutivas", side = 1, cex = 0.7, adj = 0.75) +mycorrplot(Corr.nnos, mar = c(2.5, 0, 0, 0)) +mtext(text = "Número de nós", side = 1, cex = 0.7, adj = 0.6) fonte("Fonte: Elaborado pelo autor.", cex = 0.7) @@ -686,71 +723,71 @@ pred <- data.frame(dexp = with(cottonBolls2, ##====================================================================== ## Considerando a Poisson ##------------------------------------------- -## Para nerep -aux <- predict(m2P.nerep, newdata = pred, se.fit = TRUE) -aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*"))) -aux <- data.frame(modelo = "Poisson", aux) -predP.nerep <- cbind(var = "nerep", pred, aux) -##------------------------------------------- ## Para ncapu aux <- predict(m3P.ncapu, newdata = pred, se.fit = TRUE) aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*"))) aux <- data.frame(modelo = "Poisson", aux) predP.ncapu <- cbind(var = "ncapu", pred, aux) ##------------------------------------------- +## Para nerep +aux <- predict(m2P.nerep, newdata = pred, se.fit = TRUE) +aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*"))) +aux <- data.frame(modelo = "Poisson", aux) +predP.nerep <- cbind(var = "nerep", pred, aux) +##------------------------------------------- ## Para nnos aux <- predict(m3P.nnos, newdata = pred, se.fit = TRUE) aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*"))) aux <- data.frame(modelo = "Poisson", aux) predP.nnos <- cbind(var = "nnos", pred, aux) ## -predP <- rbind(predP.nerep, predP.ncapu, predP.nnos) +predP <- rbind(predP.ncapu, predP.nerep, predP.nnos) ##====================================================================== ## Considerando a COM-Poisson ##------------------------------------------- -## Para nerep -aux <- predict(m2C.nerep, newdata = pred, type = "response", - interval = "confidence") -aux <- data.frame(modelo = "COM-Poisson", aux) -predC.nerep <- cbind(var = "nerep", pred, aux) -##------------------------------------------- ## Para ncapu aux <- predict(m3C.ncapu, newdata = pred, type = "response", interval = "confidence") aux <- data.frame(modelo = "COM-Poisson", aux) predC.ncapu <- cbind(var = "ncapu", pred, aux) ##------------------------------------------- +## Para nerep +aux <- predict(m2C.nerep, newdata = pred, type = "response", + interval = "confidence") +aux <- data.frame(modelo = "COM-Poisson", aux) +predC.nerep <- cbind(var = "nerep", pred, aux) +##------------------------------------------- ## Para nnos aux <- predict(m3C.nnos, newdata = pred, type = "response", interval = "confidence") aux <- data.frame(modelo = "COM-Poisson", aux) predC.nnos <- cbind(var = "nnos", pred, aux) ## -predC <- rbind(predC.nerep, predC.ncapu, predC.nnos) +predC <- rbind(predC.ncapu, predC.nerep, predC.nnos) ##====================================================================== ## Considerando a Quasi-Poisson ##------------------------------------------- -## Para nerep -aux <- predict(m2Q.nerep, newdata = pred, se.fit = TRUE) -aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*"))) -aux <- data.frame(modelo = "Quasi-Poisson", aux) -predQ.nerep <- cbind(var = "nerep", pred, aux) -##------------------------------------------- ## Para ncapu aux <- predict(m3Q.ncapu, newdata = pred, se.fit = TRUE) aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*"))) aux <- data.frame(modelo = "Quasi-Poisson", aux) predQ.ncapu <- cbind(var = "ncapu", pred, aux) ##------------------------------------------- +## Para nerep +aux <- predict(m2Q.nerep, newdata = pred, se.fit = TRUE) +aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*"))) +aux <- data.frame(modelo = "Quasi-Poisson", aux) +predQ.nerep <- cbind(var = "nerep", pred, aux) +##------------------------------------------- ## Para nnos aux <- predict(m3Q.nnos, newdata = pred, se.fit = TRUE) aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*"))) aux <- data.frame(modelo = "Quasi-Poisson", aux) predQ.nnos <- cbind(var = "nnos", pred, aux) ## -predQ <- rbind(predQ.nerep, predQ.ncapu, predQ.nnos) +predQ <- rbind(predQ.ncapu, predQ.nerep, predQ.nnos) ##====================================================================== ##------------------------------------------- @@ -767,6 +804,7 @@ key <- list( lines = list(lty = c(3, 1, 2), lwd = 1, col = cols), text = list(c("Poisson", "COM-Poisson", "Quasi-Poisson"))) +da$va <- relevel(da$va, "ncapu") xyplot(count ~ dexp | va, data = da, key = key, type = c("p", "g"), @@ -776,9 +814,10 @@ xyplot(count ~ dexp | va, data = da, scales = list( y = list(relation = "free", rot = 0)), strip = strip.custom( - factor.levels = c("Estruturas reprodutivas ", - "Capulhos produzidos", - "Nós da planta")), + factor.levels = c( + "Capulhos produzidos", + "Estruturas reprodutivas", + "Nós da planta")), alpha = 0.3, spread = 0.15, panel = panel.beeswarm, @@ -788,7 +827,7 @@ xyplot(count ~ dexp | va, data = da, layout = c(NA, 1), scales = list( y = list(relation = "free", rot = 0)), - type = "l", col = cols[1], lty = c(1, 3, 3), lwd = 1) + type = "l", col = cols[1], lty = c(1, 4, 4), lwd = 1) ) + as.layer( xyplot(fit + lwr + upr ~ dexp | var, data = predC, @@ -809,13 +848,12 @@ fonte.xy("Fonte: Elaborado pelo autor.") @ -Finalmente a representação gráfica na figura \ref{fig:pred-cottonBolls2} +Finalmente a representação gráfica na \autoref{fig:pred-cottonBolls2} mostra os valores preditos pelos modelos Poisson, COM-Poisson e -Quasi-Poisson com intervalos de confiança para média com 95\% de -confiança. Assim como na análise realizada na seção -\ref{sec:analise-cottonBolls}, os valores preditos com bandas de -confiança obtidos dos modelos COM-Poisson e Quasi-Poisson, são -praticamente idênticos levando as mesmas interpretações. +Quasi-Poisson com intervalos de confiança de 95\% para média. Assim como +na análise realizada na \autoref{sec:analise-cottonBolls}, os valores +preditos com bandas de confiança obtidos dos modelos COM-Poisson e +Quasi-Poisson, são idênticos, levando à s mesmas interpretações. Com esse segundo exemplo de subdispersão, em que três contagens foram realizados em um único experimento. A flexibilidade do modelo @@ -883,47 +921,47 @@ prof.ng <- profile(m2C.ng, which = "phi") @ -Nesse experimento apresentado em \ref{sec:soyaBeans}, mais de uma -variável de interesse em forma de contagem é mensurada e pela descrição -dos dados caracterÃsticas relacionadas a dispersão da contagem são -distintas em ambas (equidispersão e superdispersão). Dos modelos -apresentados no capÃtulo \ref{cap:modelos-para-dados-de-contagem}, o -Poisson, COM-Poisson, Binomial-Negativo são as alternativas paramétricas -avaliadas e o Quasi-Poisson é tomado como a alternativa -semi-paramétrica. As variáveis de interesse números de grãos de soja e -de vagens viáveis foram contabilizados por unidade experimental (vaso -com duas plantas) e estão sob o efeito, controlado, de duas covariáveis, -nÃveis de adubação potássica (\Sexpr{niveis.K}) e nÃveis de umidade do -solo (\Sexpr{niveis.umid}), que foram considerados na análise como -fatores com 5 e 3 nÃveis respectivamente. Ainda têm-se, pela condução do +Nesse experimento, mais de uma variável de interesse em forma de +contagem é mensurada. Pela descrição dos dados, realizada na +\autoref{sec:soyaBeans}, caracterÃsticas relacionadas a dispersão da +contagem são distintas em ambas as variáveis (equidispersão e +superdispersão). Dos modelos apresentados no +\autoref{cap:modelos-para-dados-de-contagem}, o Poisson, COM-Poisson, +Binomial-Negativo são as alternativas paramétricas a serem consideradas +e o Quasi-Poisson é tomado como a alternativa semi-paramétrica. As +variáveis de interesse, números de grãos de soja e de vagens viáveis, +foram contabilizados por unidade experimental (vaso com duas plantas) e +estão sob o efeito, controlado dos nÃveis de adubação potássica (0, 30, +60, 120 e 180 mg dm$^{-3}$) e dos nÃveis de umidade do solo +(37.5, 50 e 62.5\%), que foram considerados na análise como fatores +com 5 e 3 nÃveis respectivamente. Ainda têm-se, pela condução do experimento, o efeito relacionado a blocagem realizada, foram cinco blocos utilizados para controle de variação local. Os preditores considerados são \noindent Preditor 1: $\eta_1 = g(\mu_{ijk}) = \beta_0 + \tau_i + \gamma_j + \delta_k $ \\ -Preditor 1: $\eta_2 = g(\mu_{ijk}) = \beta_0 + \tau_i + \gamma_j + \delta_k + +Preditor 2: $\eta_2 = g(\mu_{ijk}) = \beta_0 + \tau_i + \gamma_j + \delta_k + \alpha_{jk}$ \noindent -em que $\tau_i$ é o efeito do i-ésimo bloco, $i=$1: bloco II, 2: bloco -III, 3: bloco IV e 4: V; $\gamma_j$ o efeito do j-ésimo nÃvel de umidade -aplicado, $j=$1: 50\% e 2: 62,5\%; $\delta_k$ o efeito do k-ésimo nÃvel -de adubação potássica, $k=$ 1: 30, 2: 60, 3: 20 e 4: 180 mg dm$^{-3}$ e -$\alpha_{jk}$ o efeito da interação entre o j-ésimo nÃvel de umidade do -solo e o k-ésimo nÃvel de adubação potássica. Assim no modelo mais -completo, com interação, são 19 parâmetros de locação a serem -estimados. - -Para ajuste dos modelos COM-Poisson nesse exemplo o tempo computacional -foi ligeiramente mais demorado (em torno de 10s para os quatro modelos -considerando as duas contagens e os dois preditores). Isso se deve ao -fato das contagens serem altas (variando entre +em que $\tau_i$ é o efeito do i-ésimo bloco ($i=$1: bloco II, 2: bloco +III, 3: bloco IV e 4: V), $\gamma_j$ o efeito do j-ésimo nÃvel de +umidade aplicado ($j=$1: 50\% e 2: 62,5\%), $\delta_k$ o efeito do +k-ésimo nÃvel de adubação potássica ($k=$ 1: 30, 2: 60, 3: 20 e 4: 180 +mg dm$^{-3}$) e $\alpha_{jk}$ o efeito da interação entre o j-ésimo +nÃvel de umidade do solo e o k-ésimo nÃvel de adubação potássica. No +modelo mais completo, com interação, são 19 parâmetros de locação a +serem estimados. + +Na abordagem via modelos COM-Poisson nesse exemplo, o tempo para ajuste +foi ligeiramente maior com relação aos exemplos anteriores. Isso se deve +ao fato das contagens serem elevadas (variando entre \Sexpr{paste(range(soyaBeans[ , "ngra"]), collapse=" e ")} para o número de grãos e \Sexpr{paste(range(soyaBeans[ , "nvag"]), collapse=" e ")} para o número de vagens) e superdispersas ($\phi<0$). Nesse cenário os incrementos da constante normalizadora $Z(\lambda_i, \nu = \exp(\phi))$, -expressão \ref{eqn:constante-z}, convergem para 0 mais lentamente. +\autoref{eqn:constante-z}, convergem para 0 mais lentamente. <<convergez-soyaBeans, fig.height=3, fig.width=6.7, fig.cap="Convergência das constantes de normalização para cada indivÃduo no modelo para o número de vagens viáveis (esquerda) e para o número de grãos produzidos (direita)">>= @@ -977,19 +1015,19 @@ lambda.ng <- lam.ng[as.numeric(obs.ng)] @ -Na figura \ref{fig:convergez-soyaBeans} são exibidos os termos dessa +Na \autoref{fig:convergez-soyaBeans} são exibidos os termos dessa constante para cada observação nos modelos mais complexos (com -interação) para o número de vagens e para o número de grãos. O critério -de convergência adotado foi de $\lambda^j/(j!)^\nu < 1 \times -10^{-3}$. No modelo para número de vagens o maior valor para a constante -foi de \Sexpr{soma.nv}, soma de \Sexpr{nterm.nv} termos, calculados para a -observação \Sexpr{obs.nv}, que teve o maior valor estimado para o -parâmetro $\lambda$, $\hat{\lambda} = $\Sexpr{lambda.nv}. Nesse o modelo -o parâmetro $\phi$ foi estimado em \Sexpr{m2C.nv@coef[1]}. Já no modelo -para o número de grãos foram necessários \Sexpr{nterm.ng} termos que -somados resultaram em \Sexpr{soma.ng}, maior constante calculada. Isso -também se deu na observação \Sexpr{obs.ng} que para este modelo, com -$\hat{\phi} = $\Sexpr{m2C.ng@coef[1]}, teve um parâmetro $\lambda$ +interação), para o número de vagens e para o número de grãos. O critério +de convergência adotado foi $\lambda^j/(j!)^\nu < 1 \times 10^{-3}$. No +modelo para número de vagens o maior valor para a constante foi de +\Sexpr{soma.nv}, soma de \Sexpr{nterm.nv} termos, calculados para a +observação \Sexpr{obs.nv}, cujo valor estimado de $\lambda$, +$\hat{\lambda} = $\Sexpr{lambda.nv}, foi o maior. Nesse o modelo o +parâmetro $\phi$ foi estimado em \Sexpr{m2C.nv@coef[1]}. Já no modelo +para o número de grãos foram necessários \Sexpr{nterm.ng} termos que, +somados, resultaram em \Sexpr{soma.ng}, maior constante calculada. Isso +também se deu para observação \Sexpr{obs.ng}, que, para este modelo com +$\hat{\phi} = $\Sexpr{m2C.ng@coef[1]}, estimou-se o parâmetro $\lambda$ estimado em \Sexpr{lambda.ng}. <<loglik-soyaBeans, include=FALSE>>= @@ -1034,78 +1072,32 @@ dispersions.ng <- c( ## Adicionando os parametros de dispersão à tabela tab.nv <- cbind(tab.nv, dispersions.nv) tab.ng <- cbind(tab.ng, dispersions.ng) +rownames(tab.ng) <- rownames(tab.nv) <- NULL ## Juntando as tabelas -## tab.ajuste <- data.frame(pred = rep(paste0("$\\eta_", 1:2, "$"), 4)) -## tab.ajuste <- cbind(tab.ajuste, as.data.frame(cbind(tab.nv, tab.ng))) -## tab.ajuste <- as.data.frame(cbind(tab.nv, tab.ng)) - -## ##---------------------------------------------------------------------- -## ## Copiar e colar o corpo do resultado na customização latex abaixo -## digits <- c(1, 0, 2, 2, -2, -2, 2, 2, -2, -2) -## print(xtable(tab.ajuste, digits = digits)) +tab.ajuste <- data.frame(pred = rep(paste0("$\\eta_", 1:2, "$"), 4)) +tab.ajuste <- cbind(tab.ajuste, as.data.frame(cbind(tab.nv, tab.ng))) +pvals <- paste(round(cmptest(m1C.nv, m2C.nv)[, 2], 3), collapse = "e") @ Medidas de qualidade de ajuste calculadas sob os modelos Poisson, COM-Poisson, Binomial Negativo e Quasi-Poisson são apresentadas na -tabela \ref{tab:ajuste-soyaBeans}. Considerando a variável resposta -número de vagens viáveis, não há indÃcios de afastamento da -equidispersão indicados i) pelos parâmetros extras dos modelos -alternativos ao Poisson, em que estimativas $\hat{\phi}$, $\hat{\theta}$ -e $\hat{\sigma^2}$ estão próximas dos valores 0, $\infty$ e 1, que -compreendem o caso particular Poisson nos modelos COM-Poisson, Binomial -Negativo e Quasi-Poisson respectivamente, ii) pelas log-verossimilhanças -dos modelos paramétricos que resultaram em valores muito próximos, iii) +\autoref{tab:ajuste-soyaBeans}. Considerando a variável resposta número +de vagens viáveis, não há indÃcios de afastamento da equidispersão +indicados i) pelos parâmetros extras dos modelos alternativos ao +Poisson, em que as estimativas $\hat{\phi}$ e $\hat{\sigma^2}$ estão +próximas dos valores 0 e 1, que compreendem o caso particular Poisson +nos modelos COM-Poisson e Quasi-Poisson respectivamente e $\hat{\theta}$ +é um valor bastante elevado (lembre-se que a Binomial Negativa se reduz +à Poisson quando $\theta \to \infty$); e ii) pelas log-verossimilhanças +dos modelos paramétricos que resultaram em valores muito próximos; iii) pelos valores de AIC que foram menores nos modelos Poisson, mostrando -que não há ganho expressivo quando estimados os parâmetros extra dos -modelos alternativos. Os \textit{p-valores} associados ao TRV entre os -modelos COM-Poisson e Poisson com preditores 1 e 2 foram de -\Sexpr{cmptest(m1C.nv, m2C.nv)[, 2]}, evidenciando a não fuga de -equidispersão dos dados. Na figura \ref{fig:prof-soyaBeans} à esquerda -são apresentados os intervalos de confiança baseados no perfil de -verossimilhança para $\phi$, no modelo COM-Poisson com efeito de -interação, como esses intervalos contém o valor da hipótese nula 0, o -modelo COM-Poisson pode ser reduzido ao Poisson. Para avaliação dos -preditores, novamente tem-se um caso de valores na borda de -significância. Nas análises que a seguir o modelo mais completo com a -interação entre adubação e umidade é considerado. - -\begin{table}[ht] -\centering -\small -\caption{Medidas de ajuste para avaliação e comparação entre preditores - e modelos ajustados ao número de vagens e ao número de grão por parcela} -\label{tab:ajuste-soyaBeans} -\begin{tabular}{lcccrcccrc} - \toprule - & & \multicolumn{4}{c}{{\bfseries Número de vagens}} & \multicolumn{4}{c}{{\bfseries Número de grãos}} \\ - \cmidrule{3-6} \cmidrule{7-10} -{PO} & np & $\ell$ & AIC & P($>\rchi^2$) & & $\ell$ & AIC & P($>\rchi^2$) & \\ - \midrule - $\eta_1$ & 11 & -266,69 & 555,38 & & & -343,16 & 708,33 & & \\ - $\eta_2$ & 19 & -259,62 & 557,23 & 7,79E-02 & & -321,67 & 681,34 & 8,83E-07 & \\[0.3cm] -{CP} & np & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\phi}$ & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\phi}$ \\ - \midrule - $\eta_1$ & 12 & -266,60 & 557,20 & & -6,75E-02 & -326,61 & 677,21 & & -8,17E-01 \\ - $\eta_2$ & 20 & -259,33 & 558,65 & 6,85E-02 & 1,29E-01 & -315,64 & 671,29 & 5,06E-03 & -5,18E-01 \\[0.3cm] -{BN} & np & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\theta}$ & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\theta}$ \\ - \midrule - $\eta_1$ & 12 & -266,69 & 557,37 & & 4,59E+03 & -326,54 & 677,07 & & 1,42E+02 \\ - $\eta_2$ & 20 & -259,62 & 559,23 & 7,82E-02 & 1,03E+06 & -315,39 & 670,77 & 4,39E-03 & 2,61E+02 \\[0.3cm] -{QP} & np & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\sigma^2}$ & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\sigma^2}$ \\ - \midrule - $\eta_1$ & 11 & 79,43 & & & 1,28E+00 & 167,71 & & & 2,71E+00 \\ - $\eta_2$ & 19 & 65,28 & & 1,87E-01 & 1,20E+00 & 124,72 & & 3,00E-02 & 2,29E+00 \\ - \bottomrule -\end{tabular} -\begin{tablenotes} - \footnotesize -\item np, número de parâmetros, PO, Poisson, CP, COM-Poisson, BN, - Binomial Negativo, QP, Quasi-Poisson. -\item Fonte: Elaborado pelo autor. -\end{tablenotes} -\end{table} +que não há ganho expressivo quando estimados os parâmetros de +dispersão/precisão nos modelos alternativos. Os \textit{p-valores} +associados ao TRV entre os modelos COM-Poisson e Poisson com preditores +1 e 2 foram \Sexpr{pvals}, evidenciando a não fuga de equidispersão dos +dados. <<prof-soyaBeans, fig.height=2.7, fig.width=5.5, out.width="0.8\\linewidth", fig.cap="Perfis de log-verossimilhança para o parâmetro de precisão da COM-Poisson nos modelos para número de vagens viáveis por parcela (esquerda) e número grãos de soja por parcela (direira).">>= @@ -1161,23 +1153,104 @@ xyplot(abs(z) ~ focal | param, data = da.soya, @ -No fragmento direito da tabela \ref{tab:ajuste-soyaBeans} são -apresentados os resultados para os modelos que ajustam os efeitos para o + +Na figura \autoref{fig:prof-soyaBeans} (à esquerda) são +apresentados os intervalos de confiança baseados no perfil de +verossimilhança para $\phi$, no modelo COM-Poisson com efeito de +interação. Como esses intervalos contém o valor 0, da hipótese nula, o +modelo COM-Poisson pode ser reduzido ao Poisson. Para avaliação dos +preditores, novamente tem-se um caso de valores próximos ao nÃvel de +significância nominal de 0,05. Nas análises a seguir o modelo mais +completo, com a interação entre adubação e umidade, é considerado. + +\begin{table}[ht] +\centering +\small +\caption{Medidas de ajuste para avaliação e comparação entre preditores + e modelos ajustados ao número de vagens e ao número de grão por parcela} +\label{tab:ajuste-soyaBeans} +\begin{tabular}{lccccccccc} + \toprule + & & \multicolumn{4}{c}{{\bfseries Número de vagens}} & \multicolumn{4}{c}{{\bfseries Número de grãos}} \\ + \cmidrule(lr){3-6} \cmidrule(lr){7-10} +{PO} & np & $\ell$ & AIC & P($>\rchi^2$) & & $\ell$ & AIC & P($>\rchi^2$) & \\ + \midrule +<<tab-soyaBeans1, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 4, 3, 2, 2, -2, 3) +print.xtable(xtable(tab.ajuste[1:2, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE, + sanitize.text.function = function(x) x) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento +{CP} & np & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\phi}$ & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\phi}$ \\ + \midrule +<<tab-soyaBeans2, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 4, 3, 2, 2, 4, 4) +print.xtable(xtable(tab.ajuste[3:4, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE, + sanitize.text.function = function(x) x) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento +{BN} & np & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\theta}$ & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\theta}$ \\ + \midrule +<<tab-soyaBeans3, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 4, -1, 2, 2, 4, -1) +print.xtable(xtable(tab.ajuste[5:6, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE, + sanitize.text.function = function(x) x) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento +{QP} & np & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\sigma^2}$ & $\ell$ & AIC & P($>\rchi^2$) & $\hat{\sigma^2}$ \\ + \midrule +<<tab-soyaBeans4, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 4, 3, 2, 2, 4, 3) +print.xtable(xtable(tab.ajuste[7:8, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE, + sanitize.text.function = function(x) x) + +@ + \bottomrule +\end{tabular} +\begin{tablenotes} + \footnotesize +\item np, número de parâmetros; PO, Poisson; CP, COM-Poisson; BN, + Binomial Negativo; e QP, Quasi-Poisson. +\item Fonte: Elaborado pelo autor. +\end{tablenotes} +\end{table} + +Na tabela \autoref{tab:ajuste-soyaBeans} (resultados à direita) são +apresentados as medidas de ajuste para os modelos considerados para o número de grãos por parcela. Neste caso há evidências de superdispersão, pois as estimativas dos parâmetros $\phi$ e $\sigma^2$ foram menor que -zero e maior que 1 respectivamente. Os valores de AIC se apresentam -menores e as avaliações da log-verossimilhança no ponto máximo maiores -para os modelos paramétricos alternativos ao Poisson. Ainda a evidência -sobre o efeito de interação para essa variável resposta é maior. Na -figura \ref{fig:prof-soyaBeans} à direita, a verossimilhança perfilhada -em $\phi$ é apresentada com indicação dos intervalos de confiança e -estes não contém o valor zero. +zero e maior que 1, respectivamente. Os valores de AIC foram menores e +as avaliações da log-verossimilhança maiores, nos modelos paramétricos +alternativos ao Poisson, quando comparados ao Poisson. Na +\autoref{fig:prof-soyaBeans} à direita, a verossimilhança perfilhada em +$\phi$ é apresentada com indicação dos intervalos de confiança e estes +não contém o valor zero. A visualização das covariâncias entre as estimativas dos parâmetros no -modelo COM-Poisson para o número de vagens por parcela é feita na figura -\ref{fig:corr-soyaBeansa} e para o número de grãos por parcela na figura -\ref{fig:corr-soyaBeansb}. Em ambos os casos a correlação entre os -parâmetros de locação ($\beta$'s) e dispersão ($\phi$) ganha destaque. +modelo COM-Poisson para o número de vagens por parcela é feita na +\autoref{fig:corr-soyaBeansa} e, para o número de grãos por parcela na +\autoref{fig:corr-soyaBeansb}. Em ambos os casos a correlação entre os +parâmetros de locação ($\beta$'s) e dispersão ($\phi$) ganha destaque, +pois há uma forte correlação, principalmente entre $\phi$ e $\beta_0$. <<corr-soyaBeansa, fig.height=7, fig.width=7, fig.cap="Imagem da matriz de correlação entre os parâmetros do modelo COM-Poisson ajustados ao número de vagens por parcela.">>= @@ -1213,14 +1286,14 @@ fonte("Fonte: Elaborado pelo autor.") @ -Na figura \ref{fig:pred-soyaBeans} são apresentadas as médias calculadas +Na \autoref{fig:pred-soyaBeans} são apresentadas as médias calculadas com intervalos de confiança 95\% sob os modelos Poisson, COM-Poisson, Binomial-Negativo e Quasi-Poisson, considerando efeito de interação entre os nÃveis de umidade do solo e adubação potássica. Tomou-se o -efeito médio de bloco, uma vez que esse efeito aditivo não é de interesse -prático. +efeito médio de bloco, uma vez que esse efeito aditivo não é de +interesse prático. -<<pred-soyaBeans, fig.height=5, fig.width=7.5, fig.cap="Valores preditos com intervalos de confiança (95\\%) como função do nÃvel de adubação com potássio e do percentual de umidade do solo para cada variável de interesse mensurada (número de vagens e número de grãos por parcela).">>= +<<pred-soyaBeans, fig.height=4.5, fig.width=7.2, fig.cap="Valores preditos com intervalos de confiança (95\\%) como função do nÃvel de adubação com potássio e do percentual de umidade do solo para cada variável de interesse mensurada (número de vagens e número de grãos por parcela).">>= library(multcomp) @@ -1331,6 +1404,7 @@ da <- reshape2::melt( variable.name = "var", value.name = "count") key <- list(type = "o", divide = 1, + columns = 2, lines = list(pch = 1:nlevels(pred.all$model) + 4, cex = 0.8), text = list(c("Poisson", "COM-Poisson", @@ -1370,20 +1444,21 @@ fonte.xy("Fonte: Elaborado pelo autor.") @ -Para a contagem do número de vagens, observa-se os intervalos com +Para a contagem do número de vagens, observa-se intervalos com comprimento muito parecidos, ligeiramente menores para o caso -COM-Poisson e Binomial Negativo. Para a contagem do número de grão por +COM-Poisson e Binomial Negativo. Para a contagem do número de grãos por parcela, um caso superdisperso, percebe-se que o modelo Poisson nos leva -a uma falsa precisão, uma vez que os intervalos são menores não por se -ajustar melhor aos dados, mas sim por subestimar a variabilidade do -processo. Para as formulações alternativas, obteve-se intervalos de -confiança para média menores nos modelos paramétricos quando comparados -com o semi-paramétrico Quasi-Poisson, isso é razoável, pois nos -Quasi-Poisson somente a especificação de dois momentos é feita, enquanto -que nos paramétricos especifica-se a distribuição completa, ganhando -informação (ver equação \ref{eqn:quasi-informacao}). De forma geral os -intervalos sob os modelos COM-Poisson e Binomial Negativa são maiores, -porém fiéis a variabilidade inerente ao processo. +a uma falsa precisão, uma vez que os intervalos são menores não pelo +modelo se ajustar melhor aos dados, mas sim por subestimar a +variabilidade do processo. Para as formulações alternativas, obteve-se +intervalos de confiança menores nos modelos paramétricos quando +comparados com os intervalos obtidos da abordagem semi-paramétrico +Quasi-Poisson. Isso é razoável, pois nos modelos Quasi-Poisson somente a +especificação de dois momentos é feita, enquanto que nos paramétricos +especifica-se a distribuição completa, ganhando informação (ver +\autoref{eqn:quasi-informacao}). De forma geral os intervalos sob os +modelos COM-Poisson e Binomial Negativa são fiéis a variabilidade +inerente ao processo. \section{Análise de ninfas de mosca-branca em lavoura de soja} \label{sec:analise-whiteFly} @@ -1421,34 +1496,34 @@ prof.ntot <- profile(m1C.ntot, which = "phi") @ -Neste experimento também há fortes indÃcios de superdispersão, conforme -visto na seção \ref{sec:whiteFly}. Assim os modelos Poisson, -COM-Poisson, Binomial Negativo e Quasi-Poisson serão aplicados. A +Nesse experimento também há fortes indÃcios de superdispersão, conforme +visto na \autoref{sec:whiteFly}. Assim os modelos Poisson, +COM-Poisson, Binomial Negativo e Quasi-Poisson foram aplicados. A variável em estudo é a contagem da quantidade de ninfas de Mosca-branca -nos folÃolos de plantas de soja, ao longo dos dias nas diferentes +nos folÃolos de plantas de soja ao longo dos dias em diferentes cultivares. Como o experimento foi conduzido sob delineamento de blocos casualizados, os efeitos de bloco são considerados no modelo. As -covariáveis serão tratadas como fator, assim como na aplicação anterior, +covariáveis foram tratadas como fator, assim como na aplicação anterior, com seis nÃveis para o número de dias decorridos a partir da primeira avaliação e quatro nÃveis para o fator cultivar de soja. Os preditores em comparação são: \noindent Preditor 1: $\eta_1 = g(\mu_{ijk}) = \beta_0 + \tau_i + \gamma_j + \delta_k $ \\ -Preditor 1: $\eta_2 = g(\mu_{ijk}) = \beta_0 + \tau_i + \gamma_j + \delta_k + +Preditor 2: $\eta_2 = g(\mu_{ijk}) = \beta_0 + \tau_i + \gamma_j + \delta_k + \alpha_{jk}$ \noindent -em que $\tau_i$ é o efeito do i-ésimo bloco, $i=$1: bloco II, 2: bloco -III, 3: bloco IV e 4: V; $\gamma_j$ o efeito da j-ésima cultivar , -$j=$1: BRS 243 RR, 2: BRS 245 RR e 3: BRS 246 RR; $\delta_k$ o efeito do -k-ésimo nÃvel do número de dias após o inÃcio do experimento, $k=$ 8, -13, 22, 31 e 38 dias e $\alpha_{jk}$ o efeito da interação entre a -j-ésima cultivar e o k-ésimo nÃvel do número de dias após o inÃcio do -experimento. A avaliação do efeito de interação é de interesse prático, -pois informa se há um padrão distinto na quantidade de ninfas ao longo -do tempo entre as cultivares. No modelo com interação, 27 -parâmetros de locação a devem ser estimados. +em que $\tau_i$ é o efeito do i-ésimo bloco ($i=$1: bloco II, 2: bloco +III, 3: bloco IV e 4: V), $\gamma_j$ o efeito da j-ésima cultivar +($j=$1: BRS 243 RR, 2: BRS 245 RR e 3: BRS 246 RR), $\delta_k$ o efeito +do k-ésimo nÃvel do número de dias após o inÃcio do experimento ($k=$1: +8, 2: 13, 3: 22, 4: 31 e 5: 38 dias) e $\alpha_{jk}$ o efeito da +interação entre a j-ésima cultivar e o k-ésimo nÃvel do número de dias +após o inÃcio do experimento. A avaliação do efeito de interação é de +interesse prático, pois informa se há um padrão distinto na quantidade +de ninfas ao longo do tempo entre as cultivares. No modelo com +interação, 27 parâmetros de locação devem ser estimados. <<convergez-prof-whiteFly, fig.height=3, fig.width=6.5, fig.cap="Convergência das constantes de normalização para cada indivÃduo (direita) e perfil de log-verossimilhança para o parâmetro extra da COM-Poisson (esquerda) no modelo para o número de ninfas de Mosca-branca.">>= @@ -1478,12 +1553,14 @@ xy1 <- xyplot(z ~ j | var, data = const.ntot, expression(frac(lambda[i]^j, "(j!)"^nu)), rot = 0), strip = strip.custom( - factor.levels = c("Número de vagens")), + factor.levels = expression("Incrementos "~Z[j])), par.settings = ps.sub) ##------------------------------------------- ## Perfil de log-verossimilhanca para phi -xy2 <- myprof(prof.ntot, subset = 4, par.settings = ps.sub) +xy2 <- myprof(prof.ntot, + namestrip = expression("Perfil para "~phi), + subset = 4, par.settings = ps.sub) print(xy1, split = c(1, 1, 2, 1), more = TRUE) print(xy2, split = c(2, 1, 2, 1), more = FALSE) @@ -1491,25 +1568,25 @@ fonte.xy("Fonte: Elaborado pelo autor.") @ -Assim como na aplicação superdispersa apresentada na seção -\ref{sec:analise-soyaBeans}, nesse exemplo tem-se um cenário com -contagens altas (variando entre \Sexpr{paste(range(soyaBeans[ , +Assim como na aplicação superdispersa apresentada na +\autoref{sec:analise-soyaBeans}, nesse exemplo tem-se um cenário com +contagens elevadas (variando entre \Sexpr{paste(range(soyaBeans[ , "ngra"]), collapse=" e ")}) e ainda superdispersas (parâmetros $\phi$ estimados próximos à -3). Isso torna a convergência da função $Z(\lambda_i, \nu = \exp(\phi))$ demorada e o valor dessa constante, que normaliza a densidade, é altÃssimo para a maioria das observações. Considerando o modelo com interação, pode-se visualizar os termos, que somados compõem a constante $Z$, para cada observação, à -direira da figura \ref{fig:convergez-prof-whiteFly}. Para a observação +esquerda da \autoref{fig:convergez-prof-whiteFly}. Para a observação \Sexpr{obs.ntot} tem-se o maior valor calculado da constante $Z$, -\Sexpr{soma.ntot}. Para obtenção deste valor \Sexpr{nterm.ntot} termos +\Sexpr{soma.ntot}. Para obtenção desse valor \Sexpr{nterm.ntot} termos foram necessários, conforme exibido no eixo $x$ do gráfico. Em problemas com contagens altas e comportamento muito superdisperso a obtenção da constante Z pode se tornar proibitiva computacionalmente, -devido à \textit{overflow} (valores que ultrapassam o limite de -capacidade de cálculo da máquina) e consequentemente o modelo -COM-Poisson não se ajusta. +devido ao problema de \textit{overflow} (valores que ultrapassam o +limite de capacidade de cálculo da máquina) e, consequentemente, o +modelo COM-Poisson não se ajusta. <<loglik-whiteFly, include=FALSE>>= @@ -1538,24 +1615,19 @@ tab.ajuste <- data.frame( pred = rep(paste("Preditor", 1:2), 4), tab.ntot) -## ##---------------------------------------------------------------------- -## ## Copiar e colar o corpo do resultado na customização latex abaixo -## digits <- c(1, 1, 0, 2, 2, 2, 0, -2, 2) -## print(xtable(tab.ajuste, digits = digits)) - @ Nesse exemplo, os modelos COM-Poisson convergiram e seus resultados são -exibidos na tabela \ref{tab:ajuste-whiteFly} em conjunto com os +exibidos na \autoref{tab:ajuste-whiteFly} em conjunto com os resultados do ajuste dos modelos Poisson, Binomial Negativo e Quasi-Poisson. Todas as estimativas dos parâmetros extras nos modelos -concorrentes ao Poisson, $\hat{\phi}$, $\hat{\theta}$ e $\hat{\sigma^2}$ -indicam expressivamente a superdispersão os dados. Em benefÃcio dos +concorrentes ao Poisson $\hat{\phi}$, $\hat{\theta}$ e $\hat{\sigma^2}$ +indicam expressivamente superdispersão. Em benefÃcio dos modelos alternativos ao Poisson tem-se todas as medidas apresentadas indicando uma substancial melhora de ajuste quando flexibilizado o modelo. Destaque para a magnitude dessas evidências, em que, por exemplo, o AIC obtido dos modelos alternativos é em torno de 0,47 -vezes o obtido do Poisson. +vezes o AIC obtido do Poisson. \begin{table}[ht] \centering @@ -1567,37 +1639,68 @@ vezes o obtido do Poisson. \toprule Poisson & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & \\ \midrule - Preditor 1 & 12 & -922,98 & 1869,96 & & & & \\ - Preditor 2 & 27 & -879,23 & 1812,46 & 87,50 & 15 & 2,90E-12 & \\[0.3cm] +<<tab-whiteFly1, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 2, 0, -2, 3) +print.xtable(xtable(tab.ajuste[1:2, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento COM-Poisson & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & $\hat{\phi}$ \\ \midrule - Preditor 1 & 13 & -410,44 & 846,89 & & & & -3,08 \\ - Preditor 2 & 28 & -407,15 & 870,30 & 6,59 & 15 & 9,68E-01 & -2,95 \\[0.3cm] +<<tab-whiteFly2, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 2, 0, 4, 3) +print.xtable(xtable(tab.ajuste[3:4, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento Binomial Neg. & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & $\hat{\theta}$ \\ \midrule - Preditor 1 & 13 & -406,16 & 838,31 & & & & 3,44 \\ - Preditor 2 & 28 & -400,55 & 857,10 & 11,21 & 15 & 7,38E-01 & 3,99 \\[0.3cm] +<<tab-whiteFly3, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 2, 0, 4, 3) +print.xtable(xtable(tab.ajuste[5:6, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento Quase-Poisson & np & deviance & AIC & F & diff np & P(>F) & $\hat{\sigma}^2$ \\ \midrule - Preditor 1 & 12 & 1371,32 & & & & & 17,03 \\ - Preditor 2 & 27 & 1283,82 & & 0,31 & 15 & 9,93E-01 & 19,03 \\ +<<tab-whiteFly4, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 2, 0, 4, 3) +print.xtable(xtable(tab.ajuste[7:8, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ \bottomrule \end{tabular} \begin{tablenotes} \footnotesize -\item np, número de parâmetros, diff $\ell$, diferença entre - log-verossimilhanças, F, estatÃstica F baseada nas quasi-deviances, +\item np, número de parâmetros; diff $\ell$, diferença entre + log-verossimilhanças; F, estatÃstica F baseada nas quasi-deviances; diff np, diferença entre o np. \item Fonte: Elaborado pelo autor. \end{tablenotes} \end{table} -Para tomada de decisão, observa-se que o modelo Poisson é claramente -inadequado. Para avaliação dos preditores, na tabela -\ref{tab:ajuste-whiteFly}, o modelo Poisson indica (com uma -significância inferior a $1 \times 10^{-10}$) que há efeito de interação -entre os dias decorridos da primeira avaliação e as cultivares ao passo -que, nos modelos alternativos, esse efeito é marcadamente não +Para tomada de decisão quanto a significância dos efeitos, observa-se +que o modelo Poisson é claramente inadequado. Para avaliação dos +preditores, na \autoref{tab:ajuste-whiteFly}, o modelo Poisson indica +(com uma significância inferior a $1 \times 10^{-10}$) que há efeito de +interação entre os dias decorridos da primeira avaliação e as cultivares +ao passo que, nos modelos alternativos, esse efeito é marcadamente não significativo. Essa discordância se deve, conforme já discutido, ao fato de o modelo Poisson subestimar a variabilidade por sua restrição de equidispersão. Assim, com variâncias menores, qualquer efeito acrescido @@ -1605,9 +1708,9 @@ ao modelo passará por significativo. Enfatizando a superdispersão indicada pelo modelo COM-Poisson e considerando o preditor de efeitos aditivos, tem-se o perfil de -verossimilhança para o parâmetro $\phi$ apresentado na figura -\ref{fig:convergez-prof-whiteFly}. Pode-se observar que os limites -inferiores dos intervalos de confiança de 90, 95 e 99\% estão muito +verossimilhança para o parâmetro $\phi$ apresentado na +\autoref{fig:convergez-prof-whiteFly}. Pode-se observar que os limites +inferiores dos intervalos de confiança de 99, 95 e 90\% estão muito distantes do valor 0, sob o qual os modelos Poisson e COM-Poisson são equivalentes. Outra caracterÃstica desse gráfico é a leve assimetria à esquerda, indicando que haverá imperfeições para inferências baseadas na @@ -1630,8 +1733,8 @@ fonte("Fonte: Elaborado pelo autor.") @ As covariâncias entre os efeitos estimados pelo modelo COM-Poisson -também são apresentadas, conforme descrição do preditor 1 na figura -\ref{fig:corr-whiteFly}, sob a escala de correlação. Similarmente as +também são apresentadas, conforme descrição do preditor 1, na +\autoref{fig:corr-whiteFly}, sob a escala de correlação. Similarmente as análises anteriores observa-se a alta correlação entre $\hat{\phi}$ e os demais parâmetros de regressão. A soma dos valores absolutos das correlações observadas entre $\hat{\phi}$ e as demais estimativas é de @@ -1696,6 +1799,7 @@ pred.all <- pred.all[with(pred.all, order(cult, aval, modelo)), ] ## Gráfico key <- list(type = "o", divide = 1, + column = 2, lines = list(pch = 1:nlevels(pred.all$model) + 4, cex = 0.8), text = list(c("Poisson", "COM-Poisson", @@ -1732,25 +1836,25 @@ fonte.xy("Fonte: Elaborado pelo autor.") @ -As médias com intervalos de confiança calculadas para cada combinação -dos nÃveis de dias após a primeira avaliação e cultivar de soja +As médias, com intervalos de confiança, calculadas para cada combinação +dos nÃveis de dias após a primeira avaliação e cultivar de soja, considerando os modelos Poisson, COM-Poisson, Binomial-Negativo e -Quasi-Poisson, são apresentadas na figura \ref{fig:pred-whiteFly}. Para -o efeito de bloco foi considerado o efeito médio, para uma correta +Quasi-Poisson, são apresentadas na \autoref{fig:pred-whiteFly}. Para o +efeito de bloco foi considerado o efeito médio, para uma correta comparação. Pode-se observar que o intervalo de confiança descrito pelo modelo Poisson é quase imperceptÃvel quando comparados aos demais, mostrando novamente que seu uso é inadequado a esses dados. Já para as -outras alternativas não tivemos um comportamento padrão em todas as -cultivares. Os intervalos pelo modelos Quasi-Poisson e COM-Poisson foram -muito similares em todos os casos e os intervalos pelo modelo Binomial -Negativo mais amplos. Um fato interessante é que não necessariamente as -estimativas pontuais da média desses modelos alternativos serão iguais, -isso ocorre, por construção, somente para nos modelos Poisson e -Quasi-Poisson, esse exemplo ilustra na prática a constatação desse -fato. Para o modelo Binomial Negativo tivemos médias visivelmente -superiores que os demais para a cultivar BRS 239. Para o modelo -COM-Poisson as estimativas pontuais são visivelmente iguais as do modelo -Poisson. +outras alternativas não tivemos um comportamento razoavelmente similar +em todas as cultivares. Os intervalos pelos modelos Quasi-Poisson e +COM-Poisson foram muito similares em todos os casos e os intervalos pelo +modelo Binomial Negativo mais amplos. Um fato interessante é que não +necessariamente as estimativas pontuais da média desses modelos +alternativos serão iguais. Isso ocorre, por construção, somente para nos +modelos Poisson e Quasi-Poisson. Esse exemplo ilustra na prática a +constatação desse fato. Para o modelo Binomial Negativo tivemos médias +visivelmente superiores que os demais para a cultivar BRS 239. Para o +modelo COM-Poisson as estimativas pontuais são aproximadamente iguais as +do modelo Poisson. \section{Análise de captura de peixes em um parque estadual} \label{sec:analise-fish} @@ -1783,14 +1887,14 @@ m2HC <- hurdlecmp(f2, data = fish) Nesse exemplo ilustra-se a análise de um estudo observacional em que aparentemente há uma quantidade excessiva de contagens nulas (veja a -seção \ref{sec:fish}). O estudo tem por objetivo a modelagem do número -de peixes capturados por grupos de visitantes em um Parque Estadual. As -covariáveis mensuradas foram (\texttt{np}), o número de pessoas no -grupo, (\texttt{nc}), o número de crianças e (\texttt{ca}) variável -binária que indica a presença ou não de um campista no grupo. +\autoref{sec:fish}). O estudo tem por objetivo a modelagem do número de +peixes capturados por grupos de visitantes em um Parque Estadual. As +covariáveis mensuradas foram o número de pessoas no grupo (\texttt{np}), +o número de crianças (\texttt{nc}) e a indicação da presença ou não de +um campista no grupo (\texttt{ca}, 0: se não presente e 1: se presente). Como já antecipado pela visualização e apresentação dos dados, modelos -estruturados de forma convencional, que pressupõe apenas um processo +estruturados de forma convencional, que pressupõem apenas um processo estocástico na geração de dados, não se ajustaram adequadamente. A seguir a alternativa de inclusão de um efeito de barreira para acomodar a quantidade excessiva de valores zero é apresentada. Os modelos Poisson, @@ -1799,8 +1903,8 @@ comparados. O número de peixes capturados é modelado em duas partes, as contagens nulas e as não nulas, conforme descrito na seção -\ref{cap02:zeros}. Abaixo define-se os preditores considerados para as -duas partes +\autoref{cap02:zeros}. Abaixo define-se os preditores considerados para +as duas partes \noindent Preditor 1: \quad $ @@ -1825,7 +1929,7 @@ Preditor 2: \quad $ \noindent sendo $g(\mu)$ e $\textrm{logit}(\pi)$ as funções de ligação que relacionam os preditores lineares com as médias dos modelos para -contagens não nulas e contagens zero respectivamente. Os preditores +contagens não nulas e contagens zero, respectivamente. Os preditores lineares foram propostos de forma aninhada. No primeiro considera-se os efeitos aditivos de todas as covariáveis mensuradas para a parte das contagens nulas e efeitos aditivos do número de pessoas e de crianças @@ -1849,10 +1953,6 @@ dispersions <- c(NA, NA, dispHB, dispHC) tab <- data.frame(pred = rep(paste("Preditor", 1:2), 3)) tab <- cbind(tab, rbind(anHP, anHB, anHC), dispersions) -## ## Gerando o código latex -## digits <- c(1, 1, 0, 2, 2, 2, 0, -2, 2) -## xtable(tab, digits = digits) - @ \begin{table}[ht] @@ -1865,40 +1965,85 @@ tab <- cbind(tab, rbind(anHP, anHB, anHC), dispersions) \toprule Poisson & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & \\ \midrule - Preditor 1 & 7 & -857,48 & 1728,96 & & & & \\ - Preditor 2 & 10 & -744,58 & 1509,17 & 225,79 & 3 & 1,12E-48 & \\[0.3cm] +<<tab-fish1, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 2, 0, -2, 3) +print.xtable(xtable(tab[1:2, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento Binomial Negativo & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & $\hat{\theta}$ \\ \midrule - Preditor 1 & 8 & -399,79 & 815,58 & & & & 0,20 \\ - Preditor 2 & 11 & -393,72 & 809,44 & 12,14 & 3 & 6,91E-03 & 0,37 \\[0.3cm] +<<tab-fish2, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 2, 0, 4, 3) +print.xtable(xtable(tab[3:4, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento COM-Poisson & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & $\hat{\phi}$ \\ \midrule - Preditor 1 & 8 & -409,85 & 835,71 & & & & -8,77 \\ - Preditor 2 & 11 & -402,30 & 826,59 & 15,12 & 3 & 1,72E-03 & -3,77 \\ +<<tab-fish3, results="asis">>= + +digits <- c(1, 1, 0, 2, 2, 2, 0, 3, 3) +print.xtable(xtable(tab[5:6, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ +\specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento \bottomrule \end{tabular} \begin{tablenotes} \footnotesize -\item np, número de parâmetros, diff $\ell$, diferença entre - log-verossimilhanças, F, estatÃstica F baseada nas quasi-deviances, +\item np, número de parâmetros; diff $\ell$, diferença entre + log-verossimilhanças; F, estatÃstica F baseada nas quasi-deviances; diff np, diferença entre o np. \\[0.1cm] \item Fonte: Elaborado pelo autor. \end{tablenotes} \end{table} -Na tabela \ref{tab:ajuste-fish} as medidas de ajuste dos modelos -Poisson, Binomial Negativo e COM-Poisson são apresentadas para -comparação dos resultados. Observa-se pelas log-verossimilhanças -maximizadas que o modelo Poisson não se ajustou adequadamente quando -comparado aos demais. Isso se deve ao fato discutido na seção -\ref{cap02:zeros}, que mesmo modelando os zeros pode-se ter diferentes -nÃveis de dispersão para as contagens nulas. Nesse exemplo as contagens -não nulas são superdispersas, conforme visto pelas estimativas dos -parâmetros extras do modelo Binomial Negativo e COM-Poisson. Indicado -pelos nÃveis descritivos dos TRV's aplicados nos modelos encaixados há -evidências de que o modelo com efeitos de interação é distinto do modelo -com efeitos aditivos definidos no preditor 1. +Na \autoref{tab:ajuste-fish} as medidas de ajuste dos modelos Poisson, +Binomial Negativo e COM-Poisson são apresentadas para comparação dos +resultados. Observa-se pelas log-verossimilhanças maximizadas que o +modelo Poisson não se ajustou adequadamente quando comparado aos +demais. Isso se deve ao fato discutido na \autoref{cap02:zeros}, que +mesmo modelando os zeros pode-se ter diferentes nÃveis de dispersão para +as contagens não nulas. Nesse exemplo as contagens não nulas são +superdispersas, visto pelas estimativas dos parâmetros extras do modelo +Binomial Negativo e COM-Poisson. Indicado pelos nÃveis descritivos dos +TRV's aplicados nos modelos encaixados, há evidências de que o modelo com +efeitos de interação é distinto do modelo com efeitos aditivos definido +no preditor 1. +As estimativas dos parâmetros para cada especificação de modelos são +exibidas na \autoref{tab:coef-fish}. Observe, primeiramente, que as +estimativas dos parâmetros $\gamma_i$, $i = 0, 1, 2, 3, 4$ são +idênticas, independentemente do modelo adotado. Esse resultado é +esperado, pois na construção dos modelos com componente de barreira, a +modelagem da parte que contempla os valores zero é realizada via +distribuição Bernoulli com parâmetro $\pi = \textrm{logit}(Z\gamma)$. As +diferenças entre os modelos ocorre na distribuição considerada para a +parte das contagens não nulas. + +\begin{table}[ht] +\centering +\caption{Estimativas dos parâmetros e razões entre as estimativa e erro + padrão para os três modelos em estudo} +\label{tab:coef-fish} +\begin{tabular}{lcccccc} + \toprule + & \multicolumn{2}{c}{Poisson} & \multicolumn{2}{c}{Binomial Negativo} & \multicolumn{2}{c}{COM-Poisson} \\ + \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} + Parâmetro & Estimativa & Est/EP & Estimativa & Est/EP & Estimativa & Est/EP \\ + \midrule <<coef-fish, results="asis">>= ##====================================================================== @@ -1920,46 +2065,21 @@ coHB <- with(summary(m2HB)$coef, { coHC <- summary(m2HC)@coef[, c(1, 3)] ## Empilha -pnames <- c("phi", paste("beta", 0:4, sep = "_"), - paste("gamma", 0:4, sep = "_")) -tab <- data.frame(pnames, cbind(coHP, coHB, coHC)) +pnames <- c("$\\sigma^2,\\,\\phi$", + paste("$\\beta_{", 0:4, "}$"), + paste("$\\gamma_{", 0:4, "}$")) -## xtable(tab) - -@ +tab <- data.frame(cbind(coHP, coHB, coHC)) +rownames(tab) <- pnames -As estimativas dos parâmetros para cada especificação de modelos são -exibidas na tabela \ref{tab:coef-fish}. Observe, primeiramente, que as -estimativas dos parâmetros $\gamma_i$, $i = 0, 1, 2, 3, 4$ são -idênticas, independentemente do modelo adotado. Esse resultado é -esperado, pois na construção dos modelos com componente de barreira, a -modelagem da parte que contempla os valores zero é realizada via -distribuição Bernoulli com parâmetro $\pi = \textrm{logit}(Z\gamma)$. As -diferenças entre os modelos ocorre na distribuição considerada para a -parte das contagens não nulas. +print.xtable(xtable(tab), + include.rownames = TRUE, + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE, + sanitize.text.function = function(x) x) -\begin{table}[ht] -\centering -\caption{Estimativas dos parâmetros e razões entre as estimativa e erro - padrão para os três modelos em estudo} -\label{tab:coef-fish} -\begin{tabular}{lcccccc} - \toprule - & \multicolumn{2}{c}{Poisson} & \multicolumn{2}{c}{Binomial Negativo} & \multicolumn{2}{c}{COM-Poisson} \\ - \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} - Parâmetro & Estimativa & Est/EP & Estimativa & Est/EP & Estimativa & Est/EP \\ - \midrule - $\phi,\,\theta$ & & & 0,37 & -2,08 & -3,77 & -9,52 \\ - $\beta_0$ & -1,01 & -5,44 & -1,75 & -2,90 & -0,62 & -29,74 \\ - $\beta_1$ & 0,74 & 7,88 & 0,41 & 1,23 & 0,10 & 29,20 \\ - $\beta_2$ & 0,89 & 18,55 & 1,05 & 6,41 & 0,14 & 21,86 \\ - $\beta_3$ & 0,49 & 1,11 & -0,06 & -0,05 & -0,33 & -17,53 \\ - $\beta_4$ & -0,45 & -3,69 & -0,32 & -0,90 & 0,04 & 33,41 \\ - $\gamma_0$ & -2,58 & -5,08 & -2,58 & -5,08 & -2,59 & -5,09 \\ - $\gamma_1$ & 0,98 & 3,00 & 0,98 & 3,00 & 1,00 & 3,04 \\ - $\gamma_2$ & 1,25 & 5,60 & 1,25 & 5,60 & 1,26 & 5,61 \\ - $\gamma_3$ & -0,93 & -1,05 & -0,93 & -1,05 & -0,93 & -1,06 \\ - $\gamma_4$ & -0,41 & -1,41 & -0,41 & -1,41 & -0,41 & -1,41 \\ +@ \bottomrule \end{tabular} \begin{tablenotes} @@ -1969,8 +2089,8 @@ parte das contagens não nulas. \end{table} Nos efeitos estimados para a parte da modelagem dos valores não nulos -têm-se algumas diferenças consideráveis. Destaca-se que o valor das -estimativas dos modelos Poisson e Binomial Negativo são comparáveis +têm-se algumas diferenças consideráveis. Destaca-se que as estimativas +dos parâmetros dos modelos Poisson e Binomial Negativo são comparáveis entre si, pois modelam a média da distribuição, mas não comparáveis com as estimativas do modelo COM-Poisson, pois este modela um parâmetro que não representa, diretamente, a média. Contudo, independente da @@ -1978,20 +2098,25 @@ distribuição o sinal dos efeitos deve ser o mesmo. Isso não ocorre nas estimativas dos parâmetros $\beta_3$, positiva no modelo Poisson e negativa nos demais e $\beta_4$, positiva no modelo COM-Poisson e negativa nos demais. Porém, esses efeitos não tem impacto significativo -para definição dos parâmetros das distribuições, conforme visto na -figura \ref{fig:pred-fish} que exibe as médias calculadas com base nas -três formulações. A seguir uma discussão sobre os valores apresentados -para os erros padrão dessas estimativas é feita. - -Calculando a magnitude desses efeitos quando escalonados pelo seu erro -padrão, obtido pelo negativo do inverso da matriz hessiana, há -diferenças substanciais. O modelo COM-Poisson indica erros padrões das +para definição dos parâmetros das distribuições, conforme pode ser visto +na \autoref{fig:pred-fish}, que exibe as médias calculadas com base nas +três formulações. A seguir uma discussão sobre os erros padrão dessas +estimativas é feita. + +Considerando a magnitude dos efeitos estimados nos modelos Hurdle, +quando escalonados pelo seu erro padrão, obtido pelo negativo do inverso +da matriz hessiana, há diferenças substanciais entre o Poisson, Binomial +Negativo e COM-Poisson. O modelo COM-Poisson indica erros padrões das estimativas muito menores que os apresentados no modelo Binomial Negativo. Sob investigações do problema, encontrou-se que este resultado se deve por inconsistências no procedimento numérico para determinação da matriz hessiana por diferenças finitas no modelo COM-Poisson. Portanto, os erros padrão sob o modelo COM-Poisson -apresentados estão incorretos. +apresentados estão incorretos. Essa impossibilidade para realização de +testes do tipo Wald no modelo Hurdle COM-Poisson foi particular da +análise desse conjunto de dados, uma possÃvel causa seja a notável +superdispersão das contagens não nulas, $\hat{\theta}=$\Sexpr{coHB[1,1]} +e $\hat{\phi}=$\Sexpr{coHC[1,1]}. <<pred-fish, fig.height=4, fig.width=6.5, fig.cap="Valores preditos do número de peixes capturados considerando o número de crianças e pessoas no grupo e a presença de um campista.">>= @@ -2081,8 +2206,8 @@ fonte.xy("Fonte: Elaborado pelo autor.") Embora tenha-se constatado problemas nos algoritmos numéricos para determinar a curvatura da log-verossimilhança, as estimativas pontuais -são coerentes com os demais modelos, conforme visto na figura -\ref{fig:pred-fish} onde são apresentadas as médias calculadas com +são coerentes com os demais modelos, conforme visto na +\autoref{fig:pred-fish} onde são apresentadas as médias calculadas com base nos três modelos estudados. Observa-se em todos os modelos a mesma tendência. @@ -2112,6 +2237,7 @@ f2 <- nema ~ log(off) + (1|cult) m1PM <- glmer(f1, data = nematodes, family = poisson) m2PM <- glmer(f2, data = nematodes, family = poisson) +##------------------------------------------- ## ## Ajuste dos mixed COM-Poisson ## m1CM <- mixedcmp(f1, data = nematodes, sumto = 50) ## m2CM <- mixedcmp(f2, data = nematodes, sumto = 50) @@ -2121,21 +2247,21 @@ load("mixedcmp_models.rda") Nessa última aplicação apresentada no trabalho a extensão dos modelos de contagem para inclusão de efeitos aleatórios é ilustrada. Os modelos em -competição são o Poisson e o COM-Poisson com efeitos aleatórios. O -conjunto de dados se refere ao número de nematoides em cultivares -medidas em soluções \texttt{sol} compostas da massa fresca de raizes -diluÃdas em água, mensuradas em gramas$ \cdot$ ml$^{-1}$ conforme -apresentado na seção \ref{sec:nematodes}. Considera-se para os modelos +considerados para análise são o Poisson e o COM-Poisson com efeitos +aleatórios. O conjunto de dados se refere ao número de nematoides, +mensurados em soluções (\texttt{sol}) compostas da massa fresca de +raizes diluÃdas em água, para diferentes cultivares, conforme +apresentado na \autoref{sec:nematodes}. Considera-se para os modelos em competição, os seguintes preditores: \noindent -Preditor 2: $g(\mu) = \beta_0 + b_j$ \\ +Preditor 1: $g(\mu) = \beta_0 + b_j$ \\ Preditor 2: $g(\mu) = \beta_0 + \beta_1 \log(\textrm{sol})_i + b_j$ \noindent em que $i = 1, 2, \cdots, \Sexpr{nrow(nematodes)}$ (número de observações) e $j$ varia nos nÃveis da cultivar de feijão ($j =$ A, B, -C, $\cdots$, S representando o efeito aleatório, realização de uma +C, $\cdots$, S) representando o efeito aleatório, realização de uma variável aleatória Normal de média 0 e variância $\sigma^2$. Assim, nos modelos propostos têm-se a variabilidade entre as cultivares explicada por uma distribuição Normal e a variabilidade dentro das cultivares @@ -2159,29 +2285,27 @@ tab <- data.frame( rbind(cbind(anP, cbind(c(NA, NA)), c(NA, NA)), cbind(anC, phi, pvs))) -## digits <- c(1, 1, 0, 2, 2, 2, 0, -2, 2, -2) -## xtable(tab, digits = digits) +digits <- c(1, 1, 0, 2, 2, 2, 0, 4, 3, 4) @ O ajuste dos modelos com a inclusão de efeitos aleatórios requer a solução de uma integral que, em geral, é resolvida numericamente. Isso torna o procedimento de ajuste computacionalmente intensivo e bastante -suscetÃvel a problemas numéricos. Para o ajuste dos modelos COM-Poisson -de efeitos mistos algumas iterações do algoritmo de estimação propuseram -valores para os parâmetros que resultaram em somas $Z(\lambda_i, \phi)$ -que não puderam ser representados pela máquina, -\textit{overflow}. Porém, o algoritmo dispõe de procedimentos que evitam -sua interrupção, propondo novos valores mesmo quando a função objetivo -não puder ser calculada, alcançando o máximo da -log-verossimilhança. Para o modelo Poisson de efeito aleatório -utilizou-se das programações em R providas pelo pacote \texttt{lme4} -\cite{Bates2015}, que trabalham com matrizes esparsas para os efeitos -aleatórios e otimização em linguagem de baixo nÃvel, minimizando os -problemas numéricos. +suscetÃvel a problemas numéricos. Em algumas iterações durante o +algoritmo de estimação dos parâmetros dos modelos COM-Poisson de efeitos +mistos, os valores considerados para os parâmetros resultaram em somas +$Z(\lambda_i, \phi)$ que não puderam ser representados pela +máquina. Porém, o algoritmo dispõe de procedimentos que evitam sua +interrupção, propondo novos valores mesmo quando a função objetivo não +puder ser calculada, alcançando o máximo da log-verossimilhança. Para o +modelo Poisson de efeito aleatório utilizou-se das programações em R +providas pelo pacote \texttt{lme4} \cite{Bates2015}, que trabalham com +matrizes esparsas para os efeitos aleatórios e otimização em linguagem +de baixo nÃvel, minimizando os problemas numéricos. Os resultados do ajuste para avaliação e comparação dos modelos são -apresentados na tabela \ref{tab:ajuste-nematodes}. Os valores na tabela +apresentados na \autoref{tab:ajuste-nematodes}. Os valores indicam que os modelos Poisson e COM-Poisson se ajustaram de forma equivalente, os valores da log-verossimilhança foram muito próximos. Essa equivalência também é apontada pelos AIC's, que foram @@ -2201,19 +2325,31 @@ significativo para explicação do número de nematoides. \toprule Poisson & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & \\ \midrule - Preditor 1 & 2 & -237,20 & 478,40 & & & & & \\ - Preditor 2 & 3 & -234,66 & 475,32 & 5,07 & 1 & 2,43E-02 & & \\[0.3cm] +<<tab-nematodes1, results="asis">>= + +print(xtable(tab[1:2, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ + \specialrule{0em}{0.5em}{0em} %% Apenas para espaçamento COM-Poisson & np & $\ell$ & AIC & 2(diff $\ell$) & diff np & P($>\rchi^2$) & $\hat{\phi}$ & P($>\rchi^2$) \\ \midrule - Preditor 1 & 3 & -236,85 & 479,71 & & & & 0,15 & 4,06E-01 \\ - Preditor 2 & 4 & -233,86 & 475,72 & 5,99 & 1 & 1,44E-02 & 0,23 & 2,05E-01 \\ +<<tab-nematodes2, results="asis">>= + +print(xtable(tab[3:4, ], digits = digits), + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE) + +@ \bottomrule \end{tabular} \begin{tablenotes} \footnotesize -\item np, número de parâmetros, diff $\ell$, diferença entre - log-verossimilhanças, - diff np, diferença entre o np. +\item np, número de parâmetros; diff $\ell$, diferença entre + log-verossimilhanças; diff np, diferença entre o np. \item Fonte: Elaborado pelo autor. \end{tablenotes} \end{table} @@ -2231,39 +2367,24 @@ fonte.xy("Fonte: Elaborado pelo autor.") Permanecendo com o segundo preditor, com o efeito do logaritmo da solução, as estimativas dos parâmetros do modelo são apresentadas na -tabela \ref{tab:coef-nematodes} em conjunto com seu erro padrão, -calculado sob aproximação quadrática da verossimilhança, ou seja via -inversão da matriz hessiana. Novamente, os resultados entre os modelos -são similares. Lembre-se que, desta tabela o único resultado comparável -diretamente é a razão entre estimativa e erro padrão do parâmetro -$\beta_1$. O parâmetro $\sigma$ é a variância da distribuição dos -efeitos aleatórios, que no modelo Poisson são somados aos efeitos fixos -para composição de $\mu$ e na COM-Poisson para composição de -$\lambda$. Outro resultado interessante dessa tabela é a estimativa do -parâmetro $\phi$ da COM-Poisson, que positiva indica uma subdispersão -moderada nesse conjunto de dados. Uma vantagem do modelo misto -COM-Poisson é que pode-se distinguir a variabilidade da contagem com a -variabilidade do efeito do grupo no experimento. Nesse exemplo tem-se -uma variabilidade do efeito aleatório maior, $\sigma$ estimado no caso +\autoref{tab:coef-nematodes}, em conjunto os erros padrão das +estimativas, calculado sob aproximação quadrática da verossimilhança, ou +seja via inversão da matriz hessiana. Novamente, os resultados dos +modelos são similares. Lembre-se que, dessa tabela, o único resultado +comparável diretamente é a razão entre estimativa e erro padrão do +parâmetro $\beta_1$. O parâmetro $\sigma$ é a variância da distribuição +dos efeitos aleatórios, que no modelo Poisson são somados aos efeitos +fixos para composição de $\mu$ e na COM-Poisson para composição de +$\lambda$. Outro resultado interessante é a estimativa do parâmetro +$\phi$ da COM-Poisson que, positiva, indica uma subdispersão moderada +nesse conjunto de dados. Uma vantagem do modelo misto COM-Poisson é que +pode-se distinguir a variabilidade da contagem da variabilidade induzida +pelo efeito do grupo no experimento. Nesse exemplo tem-se uma +variabilidade do efeito aleatório maior, $\sigma$ estimado no caso COM-Poisson maior que no caso Poisson, porém essa variabilidade extra capturada pelo efeito aleatório é compensada pela subdispersão capturada pelo parâmetro $\phi$. -<<coef-nematodes, include=FALSE>>= - -tabMP <- rbind(lsigma = c(sqrt(VarCorr(m2PM)$cult), NA, NA), - summary(m2PM)$coef[, -4], - phi = c(NA, NA, NA)) - -tabMC <- rbind(lsigma = c(exp(m2CM@coef[2]), NA, NA), - summary(m2CM)@coef[-(1:2) ,-4], - summary(m2CM)@coef[1 ,-4]) - -## print(xtable(cbind(tabMP, tabMC), digits = 2), -## include.rownames = TRUE) - -@ - \begin{table}[ht] \centering \caption{Estimativas dos parâmetros e razões entre as estimativa e erro @@ -2275,17 +2396,38 @@ tabMC <- rbind(lsigma = c(exp(m2CM@coef[2]), NA, NA), \cmidrule(lr){2-4} \cmidrule(lr){5-7} Parâmetro & Estimativa & E. Padrão & Est/EP & Estimativa & E. Padrão & Est/EP \\ \midrule - $\sigma$ & 0,75 & & & 0,93 & & \\ - $\beta_0$ & 1,62 & 0,19 & 8,50 & 2,08 & 0,45 & 4,59 \\ - $\beta_1$ & 1,27 & 0,56 & 2,29 & 1,58 & 0,68 & 2,33 \\ - $\phi$ & & & & 0,23 & 0,18 & 1,33 \\ +<<coef-nematodes, results="asis">>= + +pnames <- paste0("$", c("\\phi", "\\sigma", "\\beta_0", "\\beta_1"), "$") + +tabMP <- rbind( + phi = c(NA, NA, NA), + lsigma = c(sqrt(VarCorr(m2PM)$cult), NA, NA), + summary(m2PM)$coef[, -4]) + +tabMC <- rbind( + summary(m2CM)@coef[1 ,-4], + lsigma = c(exp(m2CM@coef[2]), NA, NA), + summary(m2CM)@coef[-(1:2) ,-4]) + +tab <- cbind(tabMP, tabMC) +rownames(tab) <- pnames + +print.xtable(xtable(tab), + include.rownames = TRUE, + include.colnames = FALSE, + hline.after = NULL, + only.contents = TRUE, + sanitize.text.function = function(x) x) + +@ \bottomrule \end{tabular} \end{table} -Como resultados complementares a tabela \ref{tab:coef-nematodes}, tem-se -os perfis de verossimilhança com intervalos de confianças de nÃveis 90, -95 e 99\% apresentados na figura \ref{fig:prof-nematodes}. Observa-se um +Como resultados complementares à \autoref{tab:coef-nematodes}, tem-se +os perfis de verossimilhança com intervalos de confianças de nÃveis 99, +95 e 90\% apresentados na \autoref{fig:prof-nematodes}. Observa-se um comportamento razoavelmente simétrico para todos os parâmetros, apenas com uma assimetria levemente destacada para o parâmetro $\beta_0$. Isso traz mais segurança na interpretação dos resultados baseados na @@ -2350,7 +2492,7 @@ parâmetros $\phi$, da distribuição considerada para a variável de contagem condicional aos efeitos aleatórios e as covariáveis e $\sigma$, da distribuição considerada para os efeitos aleatórios são conjuntamente responsáveis pela explicação da variabilidade do processo em estudo. Na -figura \ref{fig:corr-nematodes} apresentados as covariâncias entre os +\autoref{fig:corr-nematodes} são apresentados as covariâncias entre os parâmetros do modelo, na escala de correlação, a fim de verificar, principalmente, a correlação entre $\sigma$ e $\phi$. Observa-se que, conforme esperado, estes parâmetros apresentam uma forte correlação e @@ -2359,17 +2501,17 @@ que não de forma acentuada. Nota-se também que a caracterÃstica de não ortogonalidade entre os parâmetros de locação e $\phi$ se mantém com a inclusão de efeitos aleatórios. -Na figura \ref{fig:pred-nematodes} são apresentados as predições do -efeito aleatório em cada modelo, à direita e as contagem preditas para -cada cultivar e para o comportamento médio, à esquerda. A distribuição +Na \autoref{fig:pred-nematodes} são apresentados as predições do efeito +aleatório em cada modelo (à direita) e as contagem preditas para cada +cultivar e para o comportamento médio (à esquerda). A distribuição empÃrica dos efeitos aleatórios, gráfico à direita, está de acordo com -os parâmetros estimados para $\sigma$, vistos na tabela -\ref{tab:coef-nematodes}. Têm-se a ordenação dos efeitos aleatórios +os parâmetros estimados para $\sigma$, vistos na +\autoref{tab:coef-nematodes}. Têm-se a ordenação dos efeitos aleatórios idêntica em ambos os modelos, porém valores mais dispersos no caso COM-Poisson. Devido ao parâmetro adicional $\phi$ do modelo COM-Poisson, -que indica subdispersão, tem-se os valores preditos por esse modelo muito -similares aos preditos pelo modelo Poisson, conforme observa-se no -gráfico à direita da figura \ref{fig:pred-nematodes}. A soma das +que indica subdispersão, tem-se os valores preditos por esse modelo +muito similares aos preditos pelo modelo Poisson, conforme observa-se no +gráfico à direita da \autoref{fig:pred-nematodes}. A soma das diferenças ao quadrado, entre valores preditos pelos dois modelos foi de \Sexpr{sum((predCM$y - predPM$y)^2)}, o que mostra que ambos os modelos levam ao mesmo resultado. @@ -2391,7 +2533,7 @@ xy1 <- densityplot( key <- list( column = 1, lines = list(lty = c(1, 2), lwd = c(3, 1)), - text = list(c("Perfil Médio", "Perfil por cultivar"))) + text = list(c("Perfil médio", "Perfil por cultivar"))) ## Faz o gráfico nematodes2 <- rbind(nematodes, nematodes) @@ -2440,7 +2582,7 @@ Nos quatro primeiros conjuntos de dados, em que modelou-se as contagens via modelos de regressão de efeitos fixos, observou-se resultados dos modelos COM-Poisson equivalentes a abordagem semi-paramétrica via quasi-verossimilhança, quanto a significância dos efeitos e predição com -bandas de confiança. Porém ressalta-se que na abordagem por +bandas de confiança. Porém, ressalta-se que na abordagem por quasi-verossimilhança, com a especificação de apenas dois momentos, i) não se pode representar a distribuição de probabilidades da variável em estudo, ii) a informação a respeito da média é igual ou inferior a uma @@ -2449,13 +2591,13 @@ excesso de zeros e modelagem do parâmetro de dispersão não são imediatas. Nos casos de superdispersão explorou-se também os resultados dos modelos baseados na distribuição Binomial Negativa e nessa abordagem tem-se o inconveniente de somente a caracterÃstica de superdispersão ser -contemplada. Nos estudos de caso os modelos Binomial Negativo +contemplada. Nos estudos de caso, os modelos Binomial Negativo proporcionaram resultados, com relação a significância dos efeitos, -equivalentes ao COM-Poisson e Quasi-Poisson. Porém, em um dos estudos de -caso com acentuada superdispersão, os valores preditos pontuais e -intervalares nessa abordagem diferiram dos modelos COM-Poisson e -Quasi-Poisson, isso devido a forma da relação média e variância dessa -distribuição, figura \ref{fig:mv-binomneg}. +equivalentes ao COM-Poisson e Quasi-Poisson. Porém, em um dos estudos +com acentuada superdispersão, os valores preditos pontuais e +intervalares obtidos do modelo Binomial Negativo, diferiram dos modelos +COM-Poisson e Quasi-Poisson, isso devido a forma da relação média e +variância dessa distribuição, \autoref{fig:mv-binomneg}. Nas extensões propostas para o modelo COM-Poisson obteve-se resultados satisfatórios. No caso da inclusão de um componente de barreira para @@ -2470,26 +2612,26 @@ em que acomoda-se efeitos aleatórios, os procedimentos computacionalmente intensivos que são empregados no algoritmo de estimação ganham destaque. A aplicação se deu a um experimento que apresentou contagens com um grau não significativo de -subdispersão. Nessa aplicação os modelos em competição foram o Poisson e -o COM-Poisson de efeitos mistos e todos os resultados em questões -inferenciais foram equivalentes em ambos os modelos, com poder de teste -maior para o modelo COM-Poisson. +subdispersão. Nessa aplicação os modelos empregados foram o Poisson e o +COM-Poisson de efeitos mistos e todos os resultados em questões, +inferenciais, foram equivalentes em, mas com poder de teste maior para o +modelo COM-Poisson. Nas aplicações, em geral, pode-se notar caracterÃsticas que permearam a todos os modelos baseados na distribuição COM-Poisson. A primeira delas, e talvez a mais difÃcil de se contornar, é a determinação da constante de normalização, pois essa depende do parâmetro que está associado a um -preditor linear assim deve-se calcular $n$ constantes a cada iteração +preditor linear, assim deve-se calcular $n$ constantes a cada iteração do algoritmo de estimação. Em casos de contagens altas e superdispersão o cálculo dessa constante é extremamente demorado. Outra caracterÃstica que se manisfestou em todas as aplicações foi a não ortogonalidade entre os parâmetros de regressão e o parâmetro adicional $\phi$, observada -pelas correlações calculadas a partir da matriz hessiana. O que torna as +pelas correlações calculadas a partir da matriz hessiana, o que torna as inferências dependentes. Em pesquisas não relatadas nesse trabalho verificou-se que a reparametrização do parâmetro $\lambda$, adotando a -aproximação para média contorna essa caracterÃstica com o preço de se +aproximação para média, contorna essa caracterÃstica com o preço de se ter uma distribuição aproximada. Nas aplicações explorou-se também os perfis de verossimilhança para o parâmetro $\phi$ da COM-Poisson e o -comportamento aproximadamente simétrico em todos casos induz que +comportamento aproximadamente simétrico, em todos casos, induz que aproximações quadráticas da verossimilhança podem ter desempenhos satisfatórios. diff --git a/docs/cap05_consideracoes-finais.Rnw b/docs/cap05_consideracoes-finais.Rnw index 81e58dc6fddee97b9576209fe96fdc1ef650c511..ada799a81396919fac07fc99b67c9e29715706eb 100644 --- a/docs/cap05_consideracoes-finais.Rnw +++ b/docs/cap05_consideracoes-finais.Rnw @@ -3,46 +3,48 @@ % ------------------------------------------------------------------------ Os objetivos nesse trabalho foram a exploração, extensão e aplicação da -distribuição COM-Poisson na análise de dados de contagem cujo foram +distribuição COM-Poisson na análise de dados de contagem, cujo foram atendidos com a apresentação de seis aplicações dos modelos COM-Poisson -à conjuntos de dados reais que exibem equidispersão, subdispersão, +a conjuntos de dados reais que exibem equidispersão, subdispersão, superdispersão, contagens altas, excesso de zeros e efeito aleatório, mostrando a flexibilidade do modelo COM-Poisson. Das análises realizadas destaca-se a caracterÃstica restritiva do modelo -Poisson, que na maioria dos casos não se ajustou adequadamente devido a +Poisson, que na maioria dos casos não se ajustou adequadamente devido à suposição de equidispersão. Para os modelos de regressão de efeitos fixos, os resultados obtidos com as abordagens via modelo COM-Poisson, Quasi-Poisson e Binomial Negativo (para os casos de superdispersão) -foram bastante similares quanto a significância dos efeitos e predição +foram bastante similares quanto à significância dos efeitos e predição com bandas de confiança. Resultados satisfatórios também foram obtidos -para nos modelos COM-Poisson para modelagem de excesso de zeros e +nos modelos COM-Poisson com modelagem de excesso de zeros e inclusão de efeitos aleatórios. Nessas extensões, há dificuldade -computacional de ajuste dos modelos, principalmente devido ao cálculo +computacional para ajuste dos modelos, principalmente devido ao cálculo das constantes de normalização, que mesmo nos modelos de efeitos fixos -ainda são problemáticas. +se mostram como dificuldades a serem superadas. -Em todas as aplicações observou-se a não ortonalidade empÃrica na matriz -hessiana, o que se mostra como caracterÃstica da distribuição. Outra -caracterÃstica observada na análise de dados é a simetria nos perfis de -verossimilhança para o parâmetro $\phi$, indicando que aproximações -quadráticas da verossimilhança podem ter bons desempenhos. +Em todas as aplicações observou-se a não ortonalidade empÃrica, via +matriz hessiana, o que se mostra como caracterÃstica da +distribuição. Outra caracterÃstica observada na análise de dados é a +simetria nos perfis de verossimilhança para o parâmetro $\phi$, +indicando que aproximações quadráticas da verossimilhança podem ter bons +desempenhos. De forma geral, sugere-se a aplicação dos modelos COM-Poisson na análise -de dados de contagem, pois devido a sua flexibilidade, seus resultados +de dados de contagem, pois devido à sua flexibilidade, seus resultados se equivalem a abordagem semi-paramétrica via quasi-verossimilhança, porém com todos os benefÃcios da inferência totalmente paramétrica. Dado o escopo do trabalho foram vários os tópicos levantados para pesquisas futuras. Estudo de reparametrizações que tornem os parâmetros $\lambda$ e $\nu$ ortogonais no modelo COM-Poisson podem ser de grande -valia, pois tornaram as inferências entre eles independentes, além de +valia, pois tornarão as inferências entre eles independentes, além de possivelmente permitir a fatoração da verossimilhança com estimação concentrada. Para acelerar o algoritmo de estimação aproximações da constante normalização podem resultar em ajustes satisfatórios. Estudos de simulação para verificar a robustez do modelo à má especificação da distribuição da variável resposta. Implementação da modelagem de excesso -de zeros via mistura de distribuições. Inclusão de efeitos aleatórios -dependentes no modelo misto COM-Poisson. São algumas das muitas -possibilidades para pesquisa envolvendo dados de contagem subdispersos -ou superdispersos modelados com a distribuição COM-Poisson. +de zeros via mistura de distribuições. Expansão do modelo misto +COM-Poisson com diferentes fontes de efeito aleatório e efeitos +aleatórios dependentes. São algumas das muitas possibilidades para +pesquisa envolvendo dados de contagem subdispersos ou superdispersos +modelados com a distribuição COM-Poisson. diff --git a/docs/capA_codigostcc.Rnw b/docs/capA_codigostcc.Rnw index 98fd5b98732eea39ddf92add0d75dab222ac7c92..3f6c173f66b5caa93b402d7f59c1c31f2ab73c65 100644 --- a/docs/capA_codigostcc.Rnw +++ b/docs/capA_codigostcc.Rnw @@ -5,7 +5,7 @@ Todos os resultados apresentados são realizados com o \textit{software} R, cujo códigos para ajuste dos modelos COM-Poisson de efeito fixo, aleatório e com componente de barreira foram disponibilizados em formato -de pacote no endereço \url{github.com/jreduardo/tccPackage}. Nesse +de pacote no endereço \url{github.com/jreduardo/cmpreg}. Nesse apêndice são apresentados os códigos, que utilizam as funções do pacote, para produzir os resultados da seção \ref{sec:analise-cottonBolls2} (modelos de regressão de efeitos fixos). Todavia, os códigos que @@ -15,15 +15,15 @@ visualizados no complemento online <<code-cottonBolls2, echo=TRUE, eval=FALSE>>= ##---------------------------------------------------------------------- -## Instalando o pacote tccPackage, elaborado no trabalho +## Instalando o pacote cmpreg, elaborado no trabalho library(devtools) -install_git("git@github.com:JrEduardo/tccPackage.git") +install_git("git@github.com:JrEduardo/cmpreg.git") ##---------------------------------------------------------------------- ## Análise de dados apresentados na seção ... (v.a. número de nós) ## Carrega o pacote no workspace -library(tccPackage) +library(cmpreg) ## Dados data(cottonBolls2)