\documentclass{article}

\usepackage[latin1]{inputenc}
\usepackage[brazil,brazilian]{babel}

\usepackage{graphicx,Rd,harvard}

\title{Modelo Autologístico: Alguns Detalhes Na \\
Estima\c{c}\~{a}o E Infer\^{e}ncia}
\author{Elias Teixeira Krainski e Paulo Justiniano Ribeiro Junior}

\begin{document}

\maketitle

\section*{Resumo}

O modelo autologístico tem sido muito utilizado para modelagem de 
dados de incidencia de doen\c{c}as em plantas. Neste trabalho, 
consideramos a corre\c{c}\~{a}o de vi\'{e}s de David Firth para a 
fun\c{c}\~{a}o de pseudo-verossimilhan\c{c}a. Tamb\'{e}m avaliamos o 
m\'{e}todo de Monte Carlo infer\^encia sobre os par\^ametros e 
para estima\c{c}\~{a}o da vari\^{a}ncia das estimativas. Fizemos um 
estudo de simula\c{c}\~{a}o para comparar esses m\'{e}todos. 

\section{Introdu\c{c}ão}

A modelagem de dados bin\'{a}rios de imagens, presen\c{c}a de 
doen\c{c}as em plantas, mapas epidemiol\'{o}gicos, etc., \'{e} vasto 
campo de pesquisa com freq\"{u}\^{e}ntes avan\c{c}os na busca de 
formas de modelagem e principalmente de metodologias de estima\c{c}\~{a}o. 
As duas abordagens mais comuns est\~{a}o associadas ao modelo de 
regress\~{a}o log\'{i}stica. Para dados de \'{a}reas geogr\'{a}ficas, 
em dados epidemiol\'{o}gicos por exemplo, o mais comum \'{e} considerar 
a depend\^{e}ncia espacial como um efeito aleat\'{o}rio. Isso \'{e} 
feito considerando utilizando esse efeito aleat\'{o}rio espacial como 
um campo aleat\'{o}rio markoviano. No caso de dados de \textit{lattices} 
regulares, dados de imagens e doen\c{c}as de plantas por exemplo, o 
mais comum \'{e} considerar o modelo autolog\'{i}stico. Neste caso, 
a probabilidade de uma planta estar contaminada é modelada em 
fun\c{c}ão da presença ou ausência de sucesso em plantas vizinhas. 

O modelo autolog\'{i}stico foi proposto por \cite{besag:72}. 
Posteriormente, \cite{besag:74} considerou um constexto mais amplo 
de auto-modelos. A estima\c{c}\~{a}o do modelo autolog\'{i}stico por 
m\'{a}xima verossimilhan\c{c}a \'e intrat\'avel na pr\'atica. 
\cite{besag:77} sugeriu maximizar a pseudo-verossimilhan\c{c}a - PL, 
de \textit{pseudo-likelihood}. A fun\c{c}\~{a}o de PL \'{e} simplesmente 
o produto da densidade das distribui\c{c}\~{o}es condicionais de cada 
observa\c{c}\~{a}o. Nessas condi\c{c}\~{o}es a estima\c{c}\~{a}o dos 
coeficientes de regress\~{a}o do modelo \'{e} obtida simplesmente 
utilizando o m\'{e}todo comum de estima\c{c}\~{a}o em modelos lineares 
generalizado, considerando a fun\c{c}\~ao de PL como se fosse a 
fun\c{c}\~ao de verossimilhan\c{c}a dos dados.

A estima\c{c}\~ao utilizando PL n\~{a}o apresenta problemas para a 
estima\c{c}\~{a}o dos coeficientes de regress\c{c}\~{a}o do modelo. 
Por\'{e}m as estimativas de vari\^{a}ncia dos coeficientes s\~ao 
subestimadas, induzindo erros ao se fazer infer\^{e}encia. 
Uma proposta para resolver esse problema, \'{e} utilizar um 
procedimento de reamostragem \cite{gumpertz:97}. Esse procedimento 
bootstrap s\'o \'e poss\'ivel via amostrador de gibbs, onde cada 
amostra \'e obtida de sua distribui\c{c}\~ao condicional.

Alguns outros m\'etodos de estima\c{c}\~ao foram propostos para a 
estima\c{c}\~ao do modelo autolog\'istico, tais como: m\'{e}todo de PL 
generalizado, aproxima\c{c}\~oes para a verossimilhan\c{c}a exata via 
monte carlo e aproxima\c{c}\~ao estoc\'{a}stica por monte carlo via 
cadeia de markov, \cite{he+zhou+zhu,03}, \cite{friel+pettitt,04} e 
\cite{sherman+apanasovich+carroll,06}. Todos esses m\'etodos foram 
omparados com o m\'etodo PL. 
O m\'etodo de PL \'e o mais simples, por\'em inadequado para se fazer 
infer\^encia sobre os par\^ametros, por subestimar suas vari\^ancias. 
O m\'etodo de PL generalizado n\~ao \'e muito complicado, mas produz 
resultados similares ao m\'etodo PL simples. Os demais m\'etodos tem 
custo computacional elevado, devido ao m\'etodo de simula\c{a}o 
necess\'ario para a obten\c{c}ao de realiza\c{c}\~oes do processo. 

Todos os m\'etodos de estima\c{c}\~ao citados envolvem a maximiza\c{c}\~ao 
da fun\c{c}\~ao de PL. Esta fun\c{c}\~ao \'e a mesma para um modelo de 
regress\~ao log\'istica, considerando dados independentes. Por\'em, pode 
n\~ao existir estimativas\cite{silvapulle,1981}, fato \'e bastante comum 
em regress\~ao log\'istica. Isso \'e devido \`a separa\c{c}\~ao completa 
ou quase-completa, fato que ocorre quando h\'a uma covari\'avel que 
prediz perfeitamente a resposta. A separa\c{c}\~o quase-completa pode 
ocorrer se, por exemplo, 
todas as plantas n\~ao tenham vizinhas doentes ou tenham apenas uma planta 
vizinha doente. Al\'em disso, todas as plantas doentes tem uma planta 
vizinha doente. Neste caso, uma tabela de freq\"u\^encias $2\times 2$ 
apresentar\'a uma casela nula. Este fato \'e um m\'etodo simples para 
detectar o problema de separa\c{c}\~ao. Duas caselas nulas indica 
separa\c{c}\~ao completa. Os m\'etodos de simula\c{c}\~ao propostos 
para a estima\c{c}\~ao do modelo autolog\'istico envolvem 
Uma solu\c{c}\~ao para esse problema \'e o 
uso da corre\c{c}\~ao de vi\'es de David Firth, \cite{firth:92} e 
\cite{heinze+schemper:02}. 

Se o problema de separa\c{c}\~ao ocorre nos dados, as estimativas de PS 
n\~ao podem ser obtidas e torna-se necess\'ario uma modifica\c{c}\~ao 
nos algor\'itmos para os demais m\'etodos. Por\'em, utilizando os outros 
m\'etodos de estima\c{c}\~ao, o problema de separa\c{c}\~ao pode ocorrer 
em algum momento. A chance disso ocorrer aumenta se a incid\^encia \'e 
baixa e ou se a \textit{lattice} \'e pequena.

Neste artigo avaliamos duas propostas para a estima\c{c}\~ao e infer\^encia 
do modelo autolog\'istico: 1) O uso da corre\c{c}\~ao de vi\'es; e 
2) O uso do m\'etodo de monte carlo. Na se\c{c}\~ao~\ref{sec:model} 
apresentamos o modelo autolog\'istico e os m\'etodos de estima\c{a}o. 
Na sec\c{c}\~ao~\ref{sec:aplic} aplicamos os m\'etodos a um conjunto de 
dados da literatura e na se\c{c}\~ao~\ref{sec:simul} descrevemos o estudo 
de simula\c{c}ao feito e os resultados obtidos. As conclus\~{o}es s\~ao 
feitas na se\c{c}\~ao~\ref{sec:conc}

\section{Modelo autolog\'{i}stico}

O modelo de regress\~ao autolog\'istico \'e um modelo de regress\~ao 
log\'istica utilizando como covari\'aveis a informa\c{c}\~ao do 
\textit{status} das observa\c{c}\~oes vizinhas. A express\~{a}o do modelo 
pode ser dada por
\begin{equation}\label{eq:prob}\begin{array}{lcl}
logit(p_{ij}) & = & eta_{ij} = \beta_0 + \gamma_1(y_{i,j-1} + y_{i,j+1}) + 
   \gamma_2(y_{i-1,j} + y_{i+1,j})\\
p_{ij} & = & \fraq{exp\{eta_{ij}\}}{1+exp{eta_{ij}}}
\end{array}{lcl}
\end{equation}
em que $y_{i,j-1}$ e $y_{i,j+1}$ s\~ao os vizinhos na mesma linha; 
$y_{i,j-1}$ e $y_{i,j+1}$ s\~ao os vizinhos na linha adjacente; 
$\beta_0$ \'e uma constante desconhecida; $\gamma_1$ e $\gamma_2$ 
s\~ao coeficientes que indicam a import\^ancia do \textit{status} 
das plantas vizinhas. 

Esse modelo possui uma interpreta\c{c}\~ao bastante conhecida. 
Por exemplo, se $exp{\gamma_1}=2$, significa que uma planta tem o 
dobro de chance de ficar doente se houver uma planta vizinha doente 
na mesma linha em rela\c{c}\~ao a uma planta que n\~ao tem planta 
vizinha doente na mesma linha. 
O fato de considerar-mos duas covari\'aveis de vizinhan\c{c}a, 
diferenciando o efeito na linha e entre linha, torna o modelo mais 
flex\'ivel para avalizar se h\'a efeito direcional. Em muitas 
culturas \'e razo\'avel considerar essa abordagem, pois o espa\c{c}amento 
entre linhas \'e diferente do espa\c{c}amento dentro das linhas. 

A fun\c{c}\~ao de verossimilhan\c{c}a pode ser escrita na forma 
\begin{equation}
\frac{exp(Q(y))}{\display \sum_{D} exp(Q(z))}
\end{equation}
em que $Q(y)$ \'e uma soma nas densidades em cada $y_{ij}$ e $D$ \'e o 
conjunto de todas as configura\c{c}\~oes poss\'iveis. Este denominador 
\'e praticamente intrat\'avel para \textit{lattices} de grandes 
dimens\~oes. Uma possibilidade para estima\c{c}\~ao dos par\^ametros 
utilizando o m\'etodo de m\'axima verossimilhan\c{c}a \'e utilizar 
\'aproxima\c{c}\~ao utilizando um m\'etodo baseado em monte carlo via 
algoritmo de Newton-Raphson ou aproxima\c{c}\~ao estoc\'astica utilizando 
monte carlo via cadeias de markov\cite{gu+zhu,01}. 

Um m\'etodo simples para estimar os par\^ametros do modelo \'e 
maximizar a fun\c{c}\~ao de pseudo-verossimilhan\c{c}a, que \'e o 
produto das distribui\c{c}\~oes condicionais, ou, mais eficientemente 
o seu logaritmo
\begin{equation}
\display \sum_i^I \sum_j^J p_{ij}^y_{ij}(1-p_{ij})^{(1-y)}
\end{equation}
em que $p_{ij}$ \'e dado pela equa\c{c}\~ao~\ref{eq:prob}, $I$ e $J$ 
s\~ao o total de linhas e colunas da \textit{lattice}. 

O m\'etodo de PS n\~ao apresenta grandes problemas na estima\c{c}\~ao 
dos par\^ametros do modelo, por\'em n\~ao \'e adequado para estimar
a vari\^ancia das estimativas.
Uma alterativa \'e utilizar \cite{gumpertz:97} propõe o uso de 
reamostragem para estimar a vari\^ancia das estimativas. 

\section{Aplica\c{c}\~ao a dados de Pimentas de Sino}



\section{Estudo de simula\c{c}\~ao}


\section{Conclus\~oes}

A metodologia bootstrap, consiste em gerar amostras dos dados originais 
para estimar a quantidade de interesse com estas amostras. 
Porém, devido a configuração dos dados, a reamostragem só é 
possível utilizando o algoritmo amostrador de Gibbs.
A idéia é de retirar amostras da distribuição de cada 
observa\c{c}ão $y_{ij}$ condicionando a todas os outras $y_{(ij)}$.

O algoritmo de Gibbs para obter $N$ amostras $\hat{\gamma}$ é 
da seguinte forma: 
\begin{enumerate}
\item ajustar o modelo aos dados observados utilizando o método 
de pseudo-verossimilhan\c{c}a; 
\item fazer um ciclo em ordem aleatória pelas plantas atualizando 
cada planta de acordo com o modelo ajustado dos dados observados e o 
estado atual dos vizinhos; 
\item ajustar o modelo com os re-amostrados e guardar as estimativas 
  dos parâmetros; 
\item repetir os passos 2 e 3 $N$ vezes, 
\item obter as estimativas de interesse usando os $N$ vetores de 
parâmetros estimados.
\end{enumerate}

Esse método é bastante utilizado, porém pode apresentar 
problemas quando a incidência é baixa. Na próxima se\c{c}ão 
apresentamos o problema que pode ocorrer. 

\section{Problema}

Incialmente consideremos a aplica\c{c}ão do método 
apresentado na se\c{c}ão anterior ao conjunto de dados 
de Bell Pepper~\cite{gumpertz:97}. 
Esse conjunto de dados é um talhão com 20 linhas de plantas 
e em cada linha os dados foram agregados em 20 quadrat s. 
O modelo proposto é:
\begin{equation}
logit(P(Y_{ij}=1)) = \beta_{0} + \gamma_{1}R + \gamma_{2}C + 
\gamma_{3}dA + \gamma_{4}dB \;,
\end{equation}
em que $R=y_{i,j-1}+y_{i,j+1}$, $C=y_{i-1,j}+y_{i+1,j}$, 
$dA=y_{i-1,j-1}+y_{i+1,j+1}$ e $dA=y_{i+1,j+1}+y_{i-1,j-1}$. 

A análise foi feita utilizando-se o pacote \pkg{Rcitrus}, 
desenvolvido em \R{}~\cite{R:06}. 
Foram feitas 1100 re-amostras bootstrap. 
Na Figura~\ref{fig:bootps} vizualizamos o tra\c{c}o 
dos valores obtidos em cada simula\c{c}ão. 
Chamamos a aten\c{c}ão para o fato de alguns valores 
estimados serem muito discrepantes. 

\setkeys{Gin}{width=0.99\textwidth}
\begin{figure}[ht]
\begin{center}
\includegraphics{gl}
\label{fig:bootps}
\caption{Mapa dos dados de Bell Pepper (superior esquerdo) e 
tra\c{o} dos valores estimados com 1100 reamostragens}
\end{center}
\end{figure}

Nas 1100 simula\c{c}ões observamos 22 valores menores que -5 para 
o coeficiente associado à covariável de vizinhan\c{c}a na 
coluna $C$ e cinco valores para o coeficiente associado à 
covariável de vizihan\c{c}a na primeira diagonal $dA$.
\cite{gumpertz:97} fez esta análise utilizando 550 
simula\c{c}ões e descartou as primeiras 50 reamostras 
para fazer a inferência. Se nós descartarmos as 
primeiras 100 reamostras, teremos o sumário da 
Tabela~\ref{tab:sumar1} para as amostras dos coeficientes.  

\begin{table}[ht]
\begin{center}
\begin{tabular}{rrrrrr}
\hline
 & (Intercept) & R & C & dA & dB \\
\hline
Min. & -4.1540 & -2.3120 & -18.3500 & -16.8900 & -0.9619 \\
1st Qu. & -3.2050 & 0.9239 & -0.7711 & 0.1457 & 0.6832 \\
Median & -2.9770 & 1.2560 & -0.2393 & 0.5366 & 0.9930 \\
Mean & -2.9930 & 1.2590 & -0.6203 & 0.4070 & 1.0050 \\
3rd Qu. & -2.7620 & 1.5800 & 0.1912 & 0.8764 & 1.3330 \\
Max. & -1.8040 & 2.8950 & 1.6440 & 1.9570 & 2.9470 \\
Std. Error & 0.3321 & 0.5186 & 2.5017 & 1.2272 & 0.5194 \\
2.5\% & -3.6829 & 0.3000 & -1.9442 & -0.9353 & -0.0455 \\
97.5\% & -2.3785 & 2.2841 & 0.9410 & 1.4931 & 2.0169 \\
\hline
\end{tabular}
\label{tab:sumar1}
\caption{Sumário das dos coeficientes estimados 
nas 1000 simula\c{c}ões consideradas.}
\end{center}
\end{table}

Nota-se na Tabela~\ref{tab:sumar1}, que ocorreram valores de 
$-18.35$ para $\hat{\gamma_{2}}$ e $-16.89$ para $\hat{\gamma_{3}}$. 
A estimativa do erro-padrão de $\gamma_{2}$ é estimada como 
sendo 2.5017. Mas, como pode ser visto na Tabela~\ref{tab:sumar1b}, 
ao desconsiderar os 22 valores discrepantes 
obtidos, essa estimativa muda para 0.6550! 

\begin{table}[ht]
\begin{center}
\begin{tabular}{rrrrrr}
\hline
 & (Intercept) & R & C & dA & dB \\
\hline
Min. & -4.1540 & -2.3120 & -2.3790 & -2.7120 & -0.9619 \\
1st Qu. & -3.2020 & 0.9212 & -0.7395 & 0.1517 & 0.6839 \\
Median & -2.9710 & 1.2510 & -0.2211 & 0.5381 & 0.9909 \\
Mean & -2.9910 & 1.2510 & -0.2668 & 0.4753 & 1.0060 \\
3rd Qu. & -2.7570 & 1.5700 & 0.1981 & 0.8765 & 1.3330 \\
Max. & -1.8040 & 2.8950 & 1.6440 & 1.9570 & 2.9470 \\
Std. Error & 0.3329 & 0.5178 & 0.6550 & 0.6029 & 0.5161 \\
2.5\% & -3.6832 & 0.2987 & -1.6387 & -0.8900 & -0.0136 \\
97.5\% & -2.3644 & 2.2829 & 0.9471 & 1.4934 & 2.0189 \\
\hline
\end{tabular}
\label{tab:sumar1b}
\caption{Sumário das dos coeficientes estimados 
nas 978 simula\c{c}ões sem valores discrepantes consideradas.}
\end{center}
\end{table}

O procedimento ad hoc de utilizar apenas observa\c{c}ões não 
discrepantes pode ser útil para corrigir o 
problema de valores estimados muito discrepantes, porém, 
convém estudar a causa do problema. 

A regressão logística é comumente utilizada em diversas areas do 
conhecimento e é comum o problema de separa\c{c}ão 
completa e quase-completa~\cite{heinze+schemper:02}~\cite{zorn:05}. 
Uma das formas de identificar a ocorrência dese problema, 
é a identifica\c{c}ão de estimativa discrepante do coeficiente 
da variável causadora do problema de separa\c{c}ão, 
associada a valor muito grande do erro-padrão do coeficiente. 
Se a covariável é discreta, uma maneira muito simples de identificar, 
é utilizando uma tabela de contingência. Neste caso o problema é 
identificado pela ocorrência de zero em pelo menos uma das caselas. 
Na próxima se\c{c}ão apresentamos o método de estima\c{c}ão 
via verossimilhan\c{c}a penalizada que resolve o problema de 
separa\c{c}ão. 

\section{Metodologia}

\cite{firth:93} propos um método para reduzir o viés das estimativas 
de máxima verossimilhan\c{c}a. Esse método é baseado em uma 
penaliza\c{c}ão feita através da priori invariante de 
Jeffreys~\cite{firth:92}. Esta penaliza\c{c}ão na fun\c{c}ão 
de verossimilhan\c{c}a remove o termo $O(1/n$ do vício 
assintótico dos coeficientes estimados e encontra estimativas 
finitas para os coeficientes e erros-padrões. 

\section{Aplica\c{c}ão}

O procedimento de Firth está implementado para a regressão logística 
no pacote \pkg{brlr}, em \R. Esse procedimento foi incluido 
como op\c{c}ão no pacote \pkg{Rcitrus}. O modelo ajustado 
na se\c{c}ão anterior, foi também ajustado usando o método de 
verossimilhan\c{c}a penalizada, em logar do método de máxima 
verossimilhan\c{c}a. Para comparar os resultados, utilizamos 
o mesmo número de simula\c{c}ões e a mesma semente. 

\setkeys{Gin}{width=0.99\textwidth}
\begin{figure}[ht]
\begin{center}
\includegraphics{br}
\label{fig:bootbr}
\caption{Tra\c{o} dos valores estimados em 1100 reamostragens, 
utilizando o método de verossimilhan\c{c}a penalizada.}
\end{center}
\end{figure}

Na Figura~\ref{fig:bootbr}, notamos que não há valores discrepantes. 
nas estimativas dos coeficientes e na Tabela~\ref{tab:sumar2}. 

\begin{table}[ht]
\begin{center}
\begin{tabular}{rrrrrr}
\hline
 & (Intercept) & R & C & dA & dB \\
\hline
Min. & -3.8610 & -2.0970 & -3.2170 & -2.4780 & -0.7457 \\
1st Qu. & -3.0660 & 0.9196 & -0.6288 & 0.1848 & 0.6520 \\
Median & -2.8620 & 1.1840 & -0.1683 & 0.5437 & 0.9505 \\
Mean & -2.8730 & 1.2060 & -0.2344 & 0.4839 & 0.9600 \\
3rd Qu. & -2.6640 & 1.4880 & 0.1901 & 0.8439 & 1.2550 \\
Max. & -1.7830 & 2.6160 & 1.6380 & 1.8860 & 2.5950 \\
Std. Error & 0.3015 & 0.4589 & 0.6338 & 0.5500 & 0.4694 \\
2.5\% & -3.4869 & 0.3561 & -1.5340 & -0.8289 & 0.0097 \\
97.5\% & -2.3122 & 2.0922 & 0.8677 & 1.4014 & 1.9051 \\
\hline
\end{tabular}
\label{tab:sumar2}
\caption{Sumário das dos coeficientes estimados pelo 
método de verossimilhan\c{c}a penalizada 
nas 1000 simula\c{c}ões consideradas.}
\end{center}
\end{table}

\section{Conclusão}

A solu\c{c}ão para o problema bastante comum que ocorre 
em regressão logística é aplicada na metodologia de estima\c{c}ão 
do modelo autologístico. 

\bibliographystyle{./dcu}
\bibliography{./biblio}

\end{document}
