\documentclass[a4paper,titlepage]{article} \usepackage[brazilian]{babel} \usepackage{url} \usepackage[ansinew]{inputenc} \usepackage{amsthm,amsfonts} \usepackage{graphicx} \usepackage{placeins} \newcommand{\stackunder}[2]{ \renewcommand{\arraystretch}{0.3} \displaystyle \begin{array}[t]{l} {#1}\\_{#2}\end{array} \renewcommand{\arraystretch}{1}} \newcommand{\e}{\textrm{ e }} \title{MAE0418 - Estatística Documentária\\Lista 3\\Equilíbrio de Hardy-Weinberg e o sistema ABO\\Manual de Uso} \author{Fernando Henrique Ferraz Pereira da Rosa\\ Guilherme Miguel Mitne\\ Vagner Aparecido Pedro Junior} \date{24 de outubro de 2004} \begin{document} \maketitle \begin{abstract} Descrevemos no presente trabalho a modelagem através de técnicas de análise categorizada de dados estudadas no curso de Estatística Documentária, do equilíbrio de Hardy-Weinberg para tabelas do sistema de classificação ABO. Estimamos os parâmetros do modelo através de máxima verossimilhança, utilizando procedimentos iterativos. Essa estimação foi implementada no pacote estatístico R\cite{R}. Descrevemos o uso dessa implementação e discutimos alguns exemplos. \end{abstract} \section{Metodologia} \subsection{O equilíbrio de Hardy-Weinberg} Seja uma população com muitos indivíduos em que cada um possue um par particular de genes. Cada gene desse par pode ser de um dado tipo (por exemplo, a ou A). Suponhamos ainda que a os acasalamentos entre os indivíduos dessa população sejam aleatórios ou não seletivos, que a população não esteja sujeita a migrações nem a mutações constantes. Segundo o equilíbrio de Hardy-Weinberg\cite{uzunian3}, as freqüências e razões genotípicas nessa população serão constantes de geração para geração. Notemos que as suposições para a validade da lei são muito restritivas, e que na maioria dos casos portanto ela só pode se aplicada a populações teóricas. Entretanto, ela é uma importante ferramenta para os geneticistas no estudo de populações naturais. Através das fugas do equilíbrio de Hardy-Weinberg nessas populações, é possível estudar seu comportamento e modificações. É portanto de grande interesse a possibilidade de se delinear um método quantitativo para verificar se uma certa população está de acordo com a lei de Hardy-Weinberg para um dado par de genes. \subsection{O modelo probabilístico considerado} Consideremos o sistema ABO de classificação de tipo sangüíneo. Existem quatro fenótipos diferentes, ``A'', ``B'', ``O'' e ``AB'', e três genes diferentes envolvidos: $i$, $I_a$ e $I_b$. Os pares de genes e os fenótipos se relacionam de acordo com a Tabela \ref{gene-fenotipo}. \begin{table}[ht!] \centering \begin{tabular}{ll} \hline Fenótipo & Genótipos \\ \hline O & $ii$ \\ A & $I_aI_a$, $I_ai$ \\ B & $I_bI_b$, $I_bi$ \\ AB& $I_aI_b$ \\ \hline \end{tabular} \caption{Relação entre fenótipos e genótipos} \label{gene-fenotipo} \end{table} Seja $N$ o tamanho da população considerada. Vamos modelar a probabilidade de ocorrência de cada um dos fenótipos através de uma distribuição multinomial com parâmetros $N, \theta_O,\theta_A,\theta_B,\theta_{AB}$. \subsection{O modelo estrutural considerado} Sejam agora $\beta_A, \beta_B$ e $\beta_O$ as freqüências relativas na população dos genes $I_a$, $I_b$ e $i$ respectivamente. Sob o equilíbrio de Hardy-Weinberg, temos que os parâmetros $\theta_x, x = A, B, AB, O$ do modelo probabilístico considerado podem ser expressos em função dos parâmetros $\beta_x, x = A, B, O$. Consideremos então o modelo estrutural: \begin{displaymath} H : \stackunder{\theta}{\sim} = \stackunder{\theta}{\sim}\stackunder{(\beta)}{\sim} \end{displaymath} com: \begin{displaymath} H : \left( \begin{array}{l} \theta_O \\ \theta_A \\ \theta_B \\ \theta_{AB} \end{array} \right) = \left( \begin{array}{l} \beta_O^2 \\ \beta_A^2 + 2\beta_A\beta_O \\ \beta_B^2 + 2\beta_B\beta_O \\ 2\beta_A\beta_B \end{array} \right) \end{displaymath} \noindent onde, dada a restrição natural em $\stackunder{\beta}{\sim}$, podemos ainda escrever $\beta_O = 1 - \beta_A - \beta_B$. A verossimilhança sob o equilíbrio de Hardy-Weinberg é dada por: \begin{eqnarray} L(\stackunder{\beta}{\sim} | \stackunder{n}{\sim}) &=& \frac{N!}{n_O! n_A! n_B! n_{AB}!} \theta_O^{n_O} \theta_A^{n_A} \theta_B^{n_B} \theta_{AB}^{n_{AB}} \nonumber \\ &\propto& \beta_O^{2n_O} (\beta_A^2 + 2\beta_A\beta_O)^{n_A} (\beta_B^2 + 2\beta_B\beta_O)^{n_B} (2\beta_A\beta_B)^{n_{AB}} \label{verossim} \end{eqnarray} \subsection{Descrição dos métodos iterativos} A partir da verossimilhança dada em \ref{verossim}, podemos então obter os estimadores de máxima verossimilhança de $\stackunder{\beta}{\sim}$. Pelo princípio da invariância, podemos obter a partir desse estimador o estimador de máxima verossimilhança de $\stackunder{\theta}{\sim}$, de acordo com o modelo estrutural. Não é possível obter os estimadores de máxima verossimilhança de forma explícita. Consideraremos então três métodos iterativos para obtenção de $\hat{\stackunder{\beta}{\sim}}$. \subsubsection{Multiplicadores de Lagrange} Usando o método dos multiplicadores de Lagrange, obtemos a seguinte relação de recorrência para encontrar $\hat{\stackunder{\beta}{\sim}}$: \begin{displaymath} \left\{ \begin{array}{lcl} \hat{\beta_j} &=& p_j \frac{\hat{\beta_j} + \hat{\beta_O}}{\hat{\beta_j} +2 \hat{\beta_O}} + \frac{p_{AB}}{2}, \quad j= A,B \\ \hat{\beta_O} &=& p_0 + \sum_{j=A,B} \frac{p_j \hat{\beta_O}}{\hat{\beta_j} + 2\hat{\beta_O}} \end{array}\right. \end{displaymath} \noindent onde $\stackunder{p}{\sim} = \stackunder{n}{\sim}/N$. \subsubsection{Newton-Raphson} Seja $H(\stackunder{\beta}{\sim},\stackunder{n}{\sim})$ a matriz Hessiana calculada em $\stackunder{\beta}{\sim}$ e $\stackunder{n}{\sim}$, e $U(\stackunder{\beta}{\sim},\stackunder{n}{\sim})$ a função score, calculada também nos mesmos parâmetros. O método de Newton-Raphson nos provê a seguinte relação de recorrência para encontrar $\hat{\stackunder{\beta}{\sim}}$: \begin{displaymath} \stackunder{\beta^{(k)}}{\sim} = \stackunder{\beta^{(k-1)}}{\sim} - \left[ H(\stackunder{\beta^{(k-1)}}{\sim},\stackunder{n}{\sim})\right]^{-1} U(\stackunder{\beta^{(k-1)}}{\sim},\stackunder{n}{\sim}) . \end{displaymath} \subsubsection{Scoring de Fisher} Seja $\mathcal{I}(\stackunder{\beta}{\sim}) = E(-H(\stackunder{\beta}{\sim},\stackunder{n}{\sim}))$ a matriz de informação de Fisher. O método de Scoring de Fisher, define a seguinte relação iterativa para encontrarmos o estimador de máxima verossimilhança: \begin{displaymath} \stackunder{\beta^{(k)}}{\sim} = \stackunder{\beta^{(k-1)}}{\sim} + \left[ \mathcal{I}(\stackunder{\beta^{(k-1)}}{\sim})\right]^{-1} U(\stackunder{\beta^{(k-1)}}{\sim},\stackunder{n}{\sim}) . \end{displaymath} \subsubsection{Chutes Iniciais} Os três métodos propostos precisam de um valor inicial para que seja começado o processo iterativo. Nos três casos, consideraremos o estimador dos momentos para obte esse chute. Temos: \begin{eqnarray*} \beta_A^{(0)} &=& 1 - \sqrt{p_O + p_B} \\ \beta_B^{(0)} &=& 1 - \sqrt{p_O + p_A} \\ \beta_O^{(0)} &=& 1 - \beta_A^{(0)} - \beta_B^{(0)} \end{eqnarray*} \section{Uso do Programa} <>= source("hw.R") @ Implementamos o programa através de um conjunto de funções no R, na qual a interface principal para análise dos dados é a função \texttt{hardy.weinberg()}. Ela tem os parâmetros: <<>>= args(hardy.weinberg) @ O parâmetro \texttt{dados} é o vetor com os dados (na ordem O,A,B,AB) a serem analisados. Os parâmetros \texttt{eps.min} e \texttt{n.max} controlam os critérios de convergência. Eles são por padrão 0.0001 e 30. Por fim o parâmetro \texttt{method} controla qual dos três métodos iterativos será usado para obter as estimativas de máxima verossimilhança. Ao executar a função com os parâmetros desejados é criado um objeto da classe hardy.weinberg, que pode ser mostrado na tela ou guardado para futuras análises. O objeto criado contém todas as informações a respeito da análise, como o valor das estatisticas, os valores escolhidos de \texttt{n.max} e \texttt{eps.min} e o número de iterações realizadas. Para entrar com os dados, basta criar um vetor com os valores de $\stackunder{n}{\sim}$, na ordem O, A, B, AB: <<>>= sangue1 <- c(4578,4219,890,313) @ Por padrão, no caso de a função ser chamada somente com o argumento \texttt{dados}, é utilizado o método de Multiplicadores de Lagrange. <>= m1 <- hardy.weinberg(sangue1) @ O comando acima vai fazer a análise dos dados contidos em \texttt{sangue1}, definido conforme a linha anterior, utilizando o método de Multiplicadores de Lagrange com $n$ e $\epsilon$ padrões. O resultado fica guardado no objeto \texttt{m1}. Podemos então ver o resultado da análise imprimindo esse objeto na tela: <<>>= m1 @ Criamos também um método \texttt{plot} para a classe \texttt{hardy.weinberg}, de forma que além do resultado poder ser mostrado na tela, pode ser feito um gráfico a partir dele (vide Figura \ref{fig:m1}). O gráfico contém os gráficos de barra dos valores esperados e observados para cada grupo, os P-valores para as estatísticas $Q_V, Q_P$ e $Q_N$, com o gráfico da distribuição Qui-Quadrado com 1 grau de liberdade e a área preenchida indicando o P-Valor. \begin{figure}[ht!] \centering <>= plot(m1) @ \caption{Gráficos para a análise dos dados \texttt{sangue1}} \label{fig:m1} \end{figure} Por fim, criamos um método \texttt{latex}, que permite que seja gerado um relatório em \LaTeX{} a partir da análise, automaticamente. O output por padrão é impresso na tela, e daí pode ser copiado e colado em um documento em \LaTeX{} existente. Para salvar o resultado em um arquivo, basta usar-se o parâmetro \texttt{file}. O comando abaixo gera o relatório e salva-o no rquivo \texttt{m1.tex}, no diretório corrente. Pode-se então inserir esse documento dentro de um documento \LaTeX{}: <<>>= latex(m1,file="m1.tex") @ Note-se que deve ser definido no cabeçalho do documento o seguinte comando, para que a notação vetorial seja mostrada corretamente: \begin{verbatim} \newcommand{\stackunder}[2]{ \renewcommand{\arraystretch}{0.3} \displaystyle \begin{array}[t]{l} {#1}\\_{#2}\end{array} \renewcommand{\arraystretch}{1}} \end{verbatim} Inserindo os comandos diretamente no presente documento temos: <>= latex(m1) @ \section{Exemplos de Uso} Nos exemplos abaixo, assume-se que o script \texttt{hw.R} esteja no diretório em que o R esteja sendo rodado. Caso não, deve ser especificado o seu caminho completo como parâmetro da função \texttt{source()}. Por exemplo: \texttt{source("a://hw.R")} ou \texttt{source("c://home/hw.R")}. <<>>= source("hw.R") duodenal <- c(298,214,39,13) @ Com os comandos acima lemos o script na sessão corrente do R e guardamos os dados da Tabela de distribuição no vetor duodenal. Para fazer a análise através do método de Multiplicadores de Lagrange, fazemos: <>= ex1 <- hardy.weinberg(duodenal,method="lagrange") latex(ex1) @ Para usarmos Newton-Raphson, com \texttt{eps=0.00001}, fazemos\footnote{Note-se que estamos trabalhando diretamente com a saída em \LaTeX{} do programa, pois o Manual está escrito em nessa linguagem. Caso deseje-se, é possível ver o output que seria gerado na tela do R, com o comando \texttt{print(obj)} ou ainda simplesmente \texttt{obj}, numa sessão interativa.}: <>= ex2 <- hardy.weinberg(duodenal,eps.min=0.00001,method="newton") latex(ex2) @ Por fim, para usarmos o método de Scoring de Fisher: <>= ex3 <- hardy.weinberg(duodenal,method="fisher") latex(ex3) @ É possível também obter componentes individuais ao invés da saída completa. Por exemplo, para comparar os valores estimados de beta entre os três métodos: <<>>= ex1$beta ex2$beta ex3$beta @ Para ter uma lista dos componentes que podem ser obtidos: <<>>= str(ex1) @ Podemos então obter somente as estatisticas $Q_V$ de cada um dos métodos, fazendo: <<>>= ex1$estatisticas["Qv",] ex2$estatisticas["Qv",] ex3$estatisticas["Qv",] @ Podemos ainda obter os gráficos para qualquer uma das análises. Na Figura \ref{fig:ex2}, temos os gráficos para o caso da análise \texttt{ex2}. \begin{figure}[ht!] \centering <>= plot(ex2) @ \caption{Gráficos para o \texttt{ex2}} \label{fig:ex2} \end{figure} Tomemos agora o conjunto de dados: <<>>= outros <- c(25,424,844,12) @ Fazendo a análise por Fisher-Scoring, temos: <>= ex4 <- hardy.weinberg(outros,method="fisher") latex(ex4) @ \noindent e o respectivo gráfico na Figura \ref{fig:ex4}. Observamos que rejeitamos a Hipótese de Equilíbrio de Hardy-Weinberg para essa população. Notemos ainda como o número de iterações realizadas foi maior até se obter a convergência. \begin{figure}[ht!] \centering <>= plot(ex4) @ \caption{Gráficos para o \texttt{ex4}.} \label{fig:ex4} \end{figure} \FloatBarrier \bibliography{ref} \bibliographystyle{abnt-num} \input{/home/mentus/utils/latex/copyright_listas_br} \end{document}