%% Antes de processar este arquivo LaTeX (LaTeX2e) deve %% verificar que o arquivo TEMA.cls estah no mesmo %% diretorio. O arquivo TEMA.cls pode ser obtido do %% endereco www.sbmac.org.br/tema. \documentclass{TEMA} \usepackage[brazil]{babel} % para texto em Português %\usepackage[english]{babel} % para texto em Inglês \usepackage[latin1]{inputenc} % para acentuação em Português %\input{P-margin.inf} \usepackage[dvips]{graphics} \usepackage{subfigure} \usepackage{graphicx} \usepackage{epsfig} \usepackage{icomma} \newcommand{\B}{{\tt\symbol{92}}} \newcommand{\til}{{\tt\symbol{126}}} \newcommand{\chap}{{\tt\symbol{94}}} \newcommand{\agud}{{\tt\symbol{13}}} \newcommand{\crav}{{\tt\symbol{18}}} \begin{document} %******************************************************** \title {Algoritmo utilizando quadraturas gaussianas para a obtenção das probabilidades do teste bilateral de Dunnett para dados balanceados\thanks{Agradecimentos à Capes, ao CNPq e à FAPEMIG; Trabalho apresentado no 34º Congresso Nacional de Matemática Aplicada e Computacional.}} \author {S. C. BROCH \thanks{siomarabroch@jc.iffarroupilha.edu.br; doutoranda do Programa de Estatística e Experimentação Agropecuária.}\,, Departamento de Ciências Exatas, DEX, UFLA - Universidade Federal de Lavras, 37200-000 Lavras, MG, Brasil \\ \\ D. F. FERREIRA% \thanks{danielff@dex.ufla.br; Professor associado 4, Bolsista CNPq.}\,, Departamento de Ciências Exatas, DEX, UFLA - Universidade Federal de Lavras, 37200-000 Lavras, MG, Brasil.} \criartitulo \runningheads {Broch e Ferreira}{Teste bilateral de Dunnett balanceado} \begin{abstract} {\bf Resumo}. O teste de Dunnett é um teste de comparações múltiplas em que se confrontam as médias de $r$ novos tratamentos com a média de um tratamento testemunha controlando simultaneamente a taxa de erro tipo I por experimento num valor específico $\alpha$. A limitação para o seu uso é a dificuldade de obter as probabilidades da distribuição e os valores dos quantis da estatística do teste, pois as correlações possíveis entre os tratamentos têm larga amplitude. Neste trabalho é apresentado um algoritmo para obter probabilidades relacionadas ao teste de Dunnett bilateral para dados balanceados utilizando para resolver as integrais métodos numéricos de quadratura gaussiana. O algoritmo apresentou resultados precisos quando comparados com valores das tabelas divulgadas na literatura e em relação aos valores obtidos nos três programas analisados. {\bf Palavras-chave}. Procedimentos de comparações múltiplas, integração numérica, estatística experimental. \end{abstract} %******************************************************** \newsec{Introdução} O teste de Dunnett, apresentado em 1955, é um teste de comparações múltiplas em que se confrontam as médias de $r$ tratamentos novos ou de interesse com a média de um tratamento controle ou testemunha. Ele tem a propriedade de controlar simultaneamente a taxa de erro tipo I por experimento ($FWE$) num valor específico $\alpha$ \cite{HSU}. A estatística deste teste também é utilizada para determinar os intervalos de confiança dos verdadeiros valores das diferenças entre a média de cada um dos novos tratamentos e o tratamento controle, com um valor $1 - \alpha$ de coeficiente de confiança conjunta. O teste de Dunnett pode ser realizado em situações de experimentos balanceados, em que todos os tratamentos têm o mesmo número de repetições ou em que todos os novos tratamentos têm o mesmo número de repetições mesmo que diferente das repetições do tratamento controle, e não balanceados, em que os tratamentos têm diferentes números de repetições. Além disso, o teste pode ser bilateral ou unilateral. Cada situação utiliza diferentes distribuições de probabilidade para a estatística do teste. A limitação para o uso do teste de Dunnett é a dificuldade de obter as probabilidades da distribuição e os valores dos quantis da estatística do teste, pois as correlações possíveis entre os tratamentos têm larga amplitude e geralmente não se encontram tabelados valores para os casos de correlações diferentes de $0,5$. Neste trabalho é apresentado um algoritmo que utiliza quadratura gaussianas para resolver as integrais e obter as probabilidades do teste bilateral de Dunnett balanceado. \newsec{Teste de Dunnett} Considerando $Y_{j1}$, $\ldots$, $Y_{jn_j}$ uma amostra aleatória do $j$-ésimo tratamento (população) de tamanho $n_j$, em que os $Y_{ja}$'s são independentemente distribuídos para $a$ = $1$, $2$, $\ldots$, $n_j$. Da mesma forma, os $Y_{ja}$'s são independentemente distribuídos para os diferentes valores de $j$ = $1$, $2$, $\ldots$, $r+1$. Assumindo normalidade e homogeneidade de variâncias da amostra aleatória, tem-se o modelo $$Y_{ja}=\mu_j + \epsilon_{ja}, \ \ \ j = 1,\ldots,r+1, \ \ \ a=1,\ldots, n_j,$$ em que $\mu_j$ é a média do $j$-ésimo tratamento, e $\epsilon_{11},\ldots, \epsilon_{(r+1)(n_{r+1})}$ são independentes e identicamente distribuídos como uma normal com média $0$ e variância $\sigma^2$ desconhecida. Usando a notação $$\widehat{\mu}_j = \overline{Y}_j = \sum_{a=1}^{n_j} \frac{Y_{ja}}{n_j}$$ para a média amostral do j-ésimo tratamento, e $$\widehat{\sigma}^2= QME =\sum_{j=1}^{r+1}\sum_{a=1}^{n_j}\frac{(Y_{ja}-\overline{Y}_j)^2}{(n_j-1)}$$ para a variância residual amostral combinada. A distribuição da estatística para comparar $r$ médias com a testemunha, considera $\nu$ graus de liberdade associados à variância residual amostral conforme o delineamento amostral utilizado, com distribuição independente à das médias, e com diferentes correlações entre os novos tratamentos e o tratamento controle. O coeficiente de correlação entre o j-ésimo novo tratamento e o tratamento controle é representado por $\rho_j$, $j = 1,2,\ldots, r$. Sendo $n_j$ o número de repetições do $j$-ésimo novo tratamento e $n_{r+1}$ o número de repetições do tratamento controle, tem-se $\rho_j=(\sigma_{r+1}^2)/(\sigma_{r+1}^2 + \sigma_j^2)$, em que $\sigma_{r+1}^2 = \sigma_{\mu_{r+1}}^2 = \sigma^2/n_{r+1}$ é a variância da média do tratamento controle e $\sigma_j^2 = \sigma_{\mu_j}^2 = \sigma^2/n_j$ é a variância da média de cada novo tratamento, pois é considerado o caso homocedástico para o resíduo sendo a variância comum dada por $\sigma^2$. Assim, o coeficiente de correlação é $\rho_j = n_j/(n_j + n_{r+1})$. Quando os experimentos são não balanceados, $n_{r+1}\neq n_j$ e tem-se $\rho_j$ variável para cada $j$. No caso de experimentos balanceados, se $n_j = n_{r+1}=n$ para todo $j$, tem-se o caso equicorrelacionado com apenas um parâmetro, $\rho = 0,5$; se $n_j=n$ para todo $j$ mas $n \neq n_{r+1}$, tem-se o caso com apenas um parâmetro de correlação dado por $\rho = n/(n + n_{r+1})$. Para determinar quais tratamentos regulares diferem do tratamento controle, o teste bilateral de Dunnett para dados balanceados define a sua estatística do teste $|d|$ como a solução da equação: \begin{equation}\label{1} \int_0^\infty \left\{\int_{-\infty}^{+\infty}\left[\Phi\left(\frac{\sqrt{\rho} z + |d|s}{\sqrt{1 - \rho}}\right) - \Phi\left(\frac{\sqrt{\rho} z - |d|s}{\sqrt{1 - \rho}}\right)\right]^r \phi(z)dz\right\}\gamma(s) ds = 1 - \alpha, \end{equation} sendo \begin{itemize} \item $\phi(z)=\frac{1}{\sqrt{2\pi}}e^{-\frac{z^2}{2}}$ a função densidade da normal padrão; \item $\Phi(x)=\int_{-\infty}^x \frac{1}{\sqrt{2\pi}}e^{-\frac{w^2}{2}}dw$ a função de distribuição da normal padrão; \item $f(s;\nu) = \frac{\nu^{\nu/2}}{\Gamma(\nu/2)2^{\nu/2-1}}s^{\nu-1}e^{-\nu s^2/2}$, com $s\geq0$, a função densidade de $S = \widehat{\sigma}/\sigma$. \end{itemize} \subsection{Algoritmos propostos na literatura} Dunlap et al. (1981)\cite{DUN} descreveram uma função em FORTRAN IV para obtenção da função de distribuição da estatística do teste bilateral de Dunnett para dados balanceados. Nesse trabalho, as probabilidades acumuladas da normal padrão são resolvidas por aproximações fornecidas por Dunlap \& Duffy (1975 apud \cite{DUN}) ou, para valores extremos, pela aproximação de Zelen \& Severo (1965 apud \cite{DUN}). As demais integrais envolvidas no cálculo da distribuição de probabilidade de Dunnett foram resolvidas utilizando o método numérico de Newton-Côtes denominado de regra de Simpson com 20 a 30 pontos. Os autores comentam que a precisão de seus valores quando comparados aos valores das tabelas do teste de Dunnett publicadas em Hanhn and Hendrickson (1971 apud \cite{DUN}) era da ordem da quarta casa decimal. Schervish \cite{SCH} apresenta em 1983 o algoritmo AS 195 para calcular probabilidades da normal multivariada dentro de uma faixa de erro estabelecida utilizando para resolver todas as integrais envolvidas também os métodos numéricos de Newton-Côtes com sequências de subintervalos de 2 a 5 pontos. O autor baseou-se em estudos e algoritmos anteriores dados por Milton (1972), Bohrer et al. (1982) e Bohrer and Schervish (1981) (apud \cite{SCH}). A vantagem desse algoritmo sobre os anteriores foi a melhora do tempo de execução com o aumento da precisão requerida para os valores. Porém a restrição para o uso deste algoritmo é que calcula para no máximo $7$ dimensões, em função do tempo de processamento. A matriz de corre\-la\-ção deve ser não singular e passível de ser escrita como uma matriz triangular inferior. Dunnett publica em 1989 \cite{DUNNETT} o algoritmo AS 251 que possibilita calcular integrais de probabilidade normal multivariada com estrutura de correlação produto entre as variáveis, com dimensão máxima de $50$. Este algoritmo baseia-se no de Schervish, exceto pela particularidade de assumir uma estrutura de correlação especial, o que reduz consideravelmente o tempo de cálculo e viabiliza a aplicação do método com um número maior de variáveis. A restrição da estrutura de correlação produto é aplicada em muitos problemas de comparações múltiplas como o teste de Dunnett (Dunnett, 1955 apud \cite{DUNNETT})(Hsu, 1984 apud \cite{DUNNETT}). O pacote ``mvtnormpcs'' do $software$ \textsf{R} desenvolvido por Duane Currie (2006) baseou-se neste algoritmo para o cálculo das funções de distribuição $t$ e normal multivariadas. Em 1999 Genz and Bretz \cite{GENZ} apresentam um novo método para calcular probabilidades da distribuição $t$ multivariada a partir de uma integral $q$-variada transfor\-man\-do-a em um hipercubo $q-1$ dimensional que é resolvido por métodos padrões de integração numérica. Genz et al. desenvolvem e aperfeiçoam o pacote ``mvtnorm'' do $software$ \textsf{R} para o cálculo de probabilidades das distribuições $t$ e normal multivariadas, cuja última atualização data de 2012. O presente artigo foi baseado na proposta de Dunlap et al \cite{DUN} para calcular probabilidades da t-multivariada relacionadas ao teste de Dunnett bilateral para dados balanceados utilizando para resolver as integrais métodos numéricos de quadratura gaussiana ao invés da regra de Simpson. %******************************************************** \newsec{Quadraturas gaussianas} A quadratura gaussiana consiste em tomar $x_1 < x_2 < x_3 < ... $-1$ & \\ & & \\ Chebyshev do $1^\circ$ tipo & $T_n(x)$ & $\frac{1}{\sqrt{1-x^2}}$ & [-1,1] \\ & & \\ Chebyshev do $2^\circ$ tipo & $U_n(x)$ & $\sqrt{1-x^2}$ & [-1,1] \\ & & \\ Hermite físicos & $H_{nf}(x)$ & $e^{-x^2}$ & [$-\infty,\infty$] \\ & & \\ Hermite probabilísticos & $H_{np}(x)$ & $e^{-\frac{x^2}{2}}$ & [$-\infty,\infty$] \\ & & \\ Laguerre & $L_n^{(\alpha)}(x)$ & $e^{-x}x^\alpha$, $\alpha>-1$ & [$0,\infty$] \\ \hline \end{tabular} \end{small} \end{center} \end{table} Para aplicar a quadratura Gauss-Legendre, em que $\lambda(x) = 1$, quando os limites de integração em $\int_a^b \lambda(x)f(x)dx$ não são $a=-1$ e $b=1$, pode-se converter a integral para estes limites fazendo uma transformação da forma \begin{equation}\label{mvgeral} \int_a^b \lambda(x)f(x)dx \approx \frac{(b-a)}{2} \sum_{i=1}^n w_if\left(\frac{(b-a)z_i + (a+b)}{2}\right), \end{equation} em que os $z_i's$ são as raízes dos polinômios de Legendre de grau $n$. Nas quadraturas Gauss-Laguerre e Gauss-Hermite caso o integrando não apresente a função $\lambda(x) = e^{-x}$ ou $\lambda(x) = e^{-x^2}$, respectivamente, pode-se proceder ao artifício de multiplicar o integrando pelos termos $(e^{-x}e^x)$ ou $(e^{-x^2}e^{x^2})$ obtendo-se: $$\int_0^\infty g(x)dx = \int_0^\infty e^{-x} e^x g(x)dx = \int_0^\infty e^{-x} f(x) dx \approx \sum_{i=1}^n w_if(x_i),$$ ou \begin{equation}\label{2} \int_{-\infty}^\infty g(x)dx = \int_{-\infty}^\infty e^{-x^2} e^{x^2} g(x)dx = \int_{-\infty}^\infty e^{-x^2} f(x) dx \approx \sum_{i=1}^n w_if(x_i). \end{equation} Uma alternativa à aplicação das quadraturas de Gauss-Laguerre e Gauss-Hermite é transformar o intervalo de integração infinito ou semi-infinito para o intervalo $[-1,1]$ e aplicar uma quadratura Gauss-Legendre. No caso do intervalo $[0,\infty)$ pode ser transformado inicialmente para um intervalo finito $[0,1]$, dado por \begin{equation}\label{3} \int_0^\infty g(x)dx = \int_0^1 g(-ln(t))\frac{dt}{t} \end{equation} e convertendo esta nova integral para os limites $[-1,1]$, obtendo-se \begin{equation}\label{4} \int_0^1 g(-ln(t))\frac{dt}{t} = \frac{1}{2}\int_{-1}^1 \frac{1}{t} g\left(\frac{1-ln(t)}{2}\right)dt. \end{equation} %******************************************************** \newsec{Algoritmo - pseudocódigo} O algoritmo descrito calcula a probabilidade do teste bilateral de Dunnett por meio da expressão ({\bf \ref{1}}) para o quantil desejado. Ele foi implementadas no $software$ \textsf{R} por ser um programa de livre distribuição e código fonte aberto e utiliza algumas rotinas já existentes no programa. No pacote ``\textit{statmod}'' encontra-se a função ``\textit{gauss.quad}'' que possibilita a obtenção dos nós ($x_i$'$s$) e pesos ($w_i$'$s$) da quadratura desejada fornecendo como argumentos: \begin{itemize} \item o número de pontos que se deseja realizar a quadratura ($n$); \item o tipo de quadratura ($kind=$), que pode ser: ``legendre''; ``chebyshev1''; ``che\-by\-shev2''; ``hermite'' que neste programa é com os polinômios de Hermite físicos, ou seja, com $\lambda(x) = e^{-x^2}$; ``jacobi'' que neste caso é necessário especificar o valor de $\alpha$ e de $\beta$: e ``laguerre'' cujo padrão é $\alpha=0$ mas pode ser especificado outro valor de $\alpha$, contanto que seja maior do que $-1$; \item valor de $\alpha$ que por padrão é zero, nesse caso não sendo necessário especificar; \item valor de $\beta$ que por padrão é zero, nesse caso não sendo necessário especificar. \end{itemize} Outras rotinas utilizadas para calcular as principais distribuições de probabilidades foram: \begin{itemize} \item ``\textit{pnorm}'' possibilita a obtenção do valor da função de distribuição da normal padrão fornecendo o quantil desejado; \item ``\textit{lgamma}'' possibilita a obtenção do logaritmo do valor da função de densidade da probabilidade Gamma fornecendo o valor dos graus de liberdade desejado. \end{itemize} \noindent\large{\textbf{A.1. Função auxiliar para computar os valores da função da integral interna da expressão ({\bf \ref{1}}) ($gxDunInfty(x,q,c,r)$).} \begin{enumerate} \item Recebe $x$: real (peso da quadratura Gauss-Hermite), $q=|d|s[i]$, $c$ e $r$; \item calcular $$f(x) = \left[\Phi\left(\frac{\sqrt{c} x + q}{\sqrt{1 - c}}\right) - \Phi\left(\frac{\sqrt{c} x - q}{\sqrt{1 - c}}\right)\right]^r \frac{1}{\sqrt{2\pi}} \times e^{\frac{x^2}{2}}$$ \item retornar $f_x$. \end{enumerate} \noindent\large{\textbf{A.2. Quadratura Gauss-Hermite para obtenção do valor da integral interna da expressão ({\bf \ref{1}}) ($GHDun(q,c,r,n)$).} \begin{enumerate} \item Recebe $q=|d|s[i]$, $c$, $r$ e $n$; \item calcular os nós e os pesos da quadratura Gauss-Hermite por $xalpha$ $=$ $gauss.quad$ $(n, kind$ $=$ ``hermite'') que corresponde aos vetores $\mathbf{xalpha\$nodes}$ $(n \times 1)$ dos nós e $\mathbf{xalpha\$weights}$ $(n \times 1)$ dos pesos; \item chamar a função $gxDunInfty$ e criar um vetor denotado por $\textbf{dun}$ $(n \times 1)$ e calcular o seu i-ésimo componente por $dun[i] = gxDunInfty(xalpha\$nodes[i]), q, c, r)$ com $i$ = $1,$ $2,$ $\cdots,$ $n$; \item calcular $aux = \sum_{i=1}^n dun[i] \times xalpha\$weights[i]$; \item retornar $aux$. \end{enumerate} \noindent\large{\textbf{A.3. Função original da integral interna da expressão ({\bf \ref{1}}) ($f0dunnett(\textbf{s},|d|,c,r,\nu)$).} \begin{enumerate} \item Recebe $\textbf{s}$: vetor de pesos da quadratura Gauss-Legendre de dimensão $(n \times 1)$, $|d|$, $c$, $r$ e $\nu$; \item calcular o vetor $\textbf{gsdf}$ $(n \times 1)$, sendo seu i-ésimo elemento dado por $$gsdf[i] = \frac{\nu^{\nu/2}}{\Gamma(\nu/2)2^{\nu/2-1}}s[i]^{\nu-1}e^{-\nu s[i]^2/2},$$ com $i$ = $1,$ $2,$ $\cdots,$ $n$; \item chamar a função $GHDun$ e calcular o vetor $\textbf{gs}$ $(n \times 1)$ por $gs[i]=GHDun$ $(s[i] \times |d|,c,r,32)$ com $i$ = $1,$ $2,$ $\cdots,$ $n$; \item obter o vetor auxiliar $\textbf{aux}$ $(n \times 1)$ por $aux[i] = gs[i] \times gsdf[i]$ com $i$ = $1,$ $2,$ $\cdots,$ $n$; \item retornar $aux$. \end{enumerate} \noindent\large{\textbf{A.4. Função para transformar a escala de $0$ a $1$ para $0$ a $\infty$ numericamente ($f1dunnett(\textbf{s},|d|,c,r,\nu)$).} \begin{enumerate} \item Recebe $\textbf{s}$ $(n \times 1)$, $|d|$, $c$, $r$ e $\nu$; \item calcular o vetor $\textbf{trans}$ $(n \times 1)$, sendo seu i-ésimo elemento dado por $trans[i]=-ln(s[i])$, com $i$ = $1,$ $2,$ $\cdots,$ $n$; \item chamar a função $f0dunnett$ e calcular os elementos do vetor $\textbf{f1}$ $(n \times 1)$ por $f1[i]=\frac{1}{s} \times f0dunnett(trans[i],|d|,c,r,\nu)$ com $i$ = $1,$ $2,$ $\cdots,$ $n$; \item retornar $f1$. \end{enumerate} \noindent\large{\textbf{A.5. Função para transformar a escala de $-1$ a $1$ para $0$ a $1$ numericamente ($dunnett(\textbf{s},|d|,c,r,\nu)$).} \begin{enumerate} \item Recebe $\textbf{s}$ $(n \times 1)$, $|d|$, $c$, $r$ e $\nu$; \item calcular o vetor $\textbf{trans}$ $(n \times 1)$, sendo seu i-ésimo elemento dado por $trans[i]=\frac{1}{2} \times s[i] + \frac{1}{2}$, com $i$ = $1,$ $2,$ $\cdots,$ $n$; \item chamar a função $f1dunnett$ e calcular os elementos do vetor $\textbf{f2}$ $(n \times 1)$ por $f2[i]=\frac{1}{2} \times f1dunnett(trans[i],|d|,c,r,\nu)$ com $i$ = $1,$ $2,$ $\cdots,$ $n$; \item retornar $f2$. \end{enumerate} \noindent\large{\textbf{A.6. Função principal para calcular a expressão ({\bf \ref{1}}) ($GLdunnett$ $(|d|$, $c$, $r$, $\nu$, $n)$).} \begin{enumerate} \item Recebe $|d|$, $c$, $r$, $\nu$ e $n$; \item se $(df>1000)\&(n<200)$ faça $n=1000$ senão se $(df>10000)\&(n<500)$ faça $n=1000$; \item se $df=\infty$, então retorne $GHDun(|d|,r,c,32)$; senão vá para o passo 4; \item calcular os nós e os pesos da quadratura Gauss-Legendre usando a função ``gauss.quad'' do R por $xalpha=gauss.quad(n,kind=$ ``legendre''$)$; \item chamar a função $dunnett$ para cada nó e calcular o vetor $\textbf{dunn}$, sendo seu i-ésimo elemento dado por $$dunn[i]=dunnett(xalpha\$nodes[i],|d|,r,c,\nu),$$ com $i$ = $1,$ $2,$ $\cdots,$ $n$; \item calcular $aux = dunn[i] \times xalpha\$weights[i]$ que é a probabilidade desejada. \end{enumerate} \newsec{Resultados} Para aplicar o teste de Dunnett deve-se seguir os seguintes passos: \begin{enumerate} \item Para cada contraste determinar a estimativa $\widehat{d}_i =\frac{\widehat{\mu}_j- \widehat{\mu}_{r+1}}{\widehat{\sigma}^2\sqrt{\frac{1}{n}+\frac{1}{n_{r+1}}}}$; \item entrar com as informações como argumentos do algoritmo implementado: o valor $\widehat{d}_i$, o valor da correlação $c = n/(n + n_{r+1})$, o número de tratamentos $r$, excluindo-se o controle, os graus de liberdade $\nu$ e o número de pontos $n$ para a quadratura; \item a rotina fornece a probabilidade conjunta estimada $1-p$, da qual se obtém o valor-$p$ do teste; \item compara-se o valor-$p$ obtido no contraste com o nível de significância nominal $\alpha$ do teste; \item conclui-se sobre a significância da diferença de médias observada da seguinte forma: se o valor-$p$ for inferior ao $\alpha$, considera-se que o novo tratamento e o controle diferem; caso contrário, considera-se que não há diferença significativa entre os dois tratamentos. \end{enumerate} Repete-se este procedimento para cada um dos $r$ contrastes. Assim, determina-se quais dos $r$ novos tratamentos apresentam média significativamente diferente da média do tratamento controle, em um nível de significância nominal pré-especificado para as $r$ comparações simultaneamente. Na Tabela \ref{tab:compar} estão alguns valores de probabilidades do teste bilateral de Dunnett balanceado, calculados pelo presente algoritmo, em função de $|d|$, $\rho$, $\nu$ e $k$, e também daqueles obtidos em tabelas fornecidas em HSU \cite{HSU}, calculados pelo algoritmo de Dunlap et al. \cite{DUN}, pelos pacotes ``mvtnormpcs'' \cite{DUNNETT} e ``mvtnorm'' \cite{GENZ} do $software$ \textsf{R}. \begin{table}[h] \caption{Probabilidade do teste bilateral de Dunnett balanceado em função de $|d|$, $\rho$, $\nu$ e $k$}\label{tab:compar} \begin{footnotesize} \begin{center} \begin{tabular}{r|r|r|r|r|r} \hline $\rho$ & $\nu$ & \multicolumn{4}{c}{Probabilidade da t-multivariada no intervalo $(-d, d)$}\\ \cline{3-6} & & alg. proposto & alg. Dunlap et al. & ``mvtnormpcs'' & ``mvtnorm''\\ \hline \hline \multicolumn{6}{c}{$k=3$; $|d|=3,391$}\\ \hline 0,5 & 16 & $^{(**)}$0,9900151883803 & 0.9900126692757 & 0,9900016784668 & 0,9900441397330\\ 0,3 & 24 & 0,9930189750833 & 0,9930163962220 & 0,9930208325386 & 0,9930297267354\\ \hline \multicolumn{6}{c}{$k=3$; $|d|=2,000$}\\ \hline 0,5 & 30& 0,8661810780890 & 0,8661786202006 & 0,8661742210388 & 0,8661298672577\\ \hline \multicolumn{6}{c}{$k=5$; $|d|=3,052$}\\ \hline 0,5 & 9 & $^{(*)}$0,9500095172420 & 0,9500061659667 & 0,9500006437302 & 0,9500525922825\\ 0,6 & 29 & 0,9812681036455 & 0,9812654899493 & 0,9812633395195 & 0,9812400176518\\ \hline \multicolumn{6}{c}{$k=12$; $|d|=3,504$}\\ \hline 0,5 & 45 & $^{(**)}$0,9900110737664 & 0,9900084317938 & 0,9900038242340 & 0,9899819900391\\ \hline \multicolumn{6}{c}{$k=16$; $|d|=2,881$}\\ \hline 0,5 & 160 & $^{(*)}$0,9500110840683 & 0,9500084900263 & 0,9500097036362 & 0,9499497774682\\ 0,8 & 200 & 0,9735943018411 & 0,9735914613436 & 0,9735981822014 & 0,9733991646750\\ \hline \multicolumn{6}{c}{$k=48$; $|d|=2,000$}\\ \hline 0,5 & 900 & 0,4203154083059 & 0,4203149222736 & 0,4203275144100 & $ ^{(***)}$0,4195540530385\\ \hline \multicolumn{6}{c}{$k=60$; $|d|=3,000$}\\ \hline 0,5 & 900 & 0,9161963127005 & 0,9161937957534 & NC & 0,9161224918813\\ \hline \end{tabular} \end{center} Nota: segundo tabelas do HSU (*) $p=0,95$ e (**) $p=0,99$; (***) valor com indicação de erro maior do que o valor absoluto do eps mínimo aceito na rotina. \end{footnotesize} \end{table} Quando comparados os valores da probabilidade obtidos no algoritmo proposto com valores tabelados e em relação aos valores obtidos nos três programas analisados o algoritmo proposto utilizando quadraturas gaussianas na resolução das integrais forneceu resultados semelhantes na quarta decimal; É difícil falar em precisão de valores das probabilidades pois não se tem o valor exato dessas probabilidades. Os valores obtidos nas tabelas HSU \cite{HSU} podem ser imprecisos, pois dependem da qualidade da implementação: algoritmos antigos e sem uso de computador. Uma quadratura de $500$ nós fornece os resultados semelhantes na precisão quanto o uso de quadratura de $50$ nós. Considerando que as quadraturas são mais precisas do que métodos de simulação Monte-Carlo ou quasi-Monte Carlo e métodos numéricos baseados em Newton-Côtes, espera-se que os resultados obtidos pelo algoritmo proposto seja mais preciso do que os dos concorrentes. O fato de que a quadratura aumenta a precisão conforme o aumento dos pontos fortalece esta conjectura. \newsec{Conclusões} As probabilidades da $t$-multivariada relacionadas ao teste de Dunnett bilateral para dados balanceados foram obtidas com sucesso por meio de métodos numéricos de quadratura gaussiana. As comparações com outros métodos, mostram que a acurácia das quadraturas Gaussianas neste contexto são elevadas e comparáveis com todos os métodos discutidos no presente trabalho. Essa acurácia pode ser aumentada com o aumento do número de pontos, se isso for demandado. \vspace*{.3cm} \begin{abstract} {\bf Abstract}. Dunnett's test is appropriate for multiple comparisons where the means of $r$ new treatments are compared with a control treatment mean. this test controls the experimentwise type I error rate for a nominal significance level $\alpha$. The limitation in their use is the difficulty of obtaining computationally cumulative distribution function and quantiles of the test statistic, since the possible correlations between treatments have broad amplitude. This paper presents an algorithm for probabilities related to bilateral Dunnett test for balanced data using numerical methods to solve the integrals by Gaussian quadrature. The algorithm showed accurate results when compared with values reported in the literature and with values obtained in three computational programs.\\ {\bf Keywords}. Multiple comparison procedures, numerical integration, experimental statistic. \end{abstract} \begin{thebibliography}{9} \bibitem{DUN} W. P. Dunlap; M. S. Marx; G. J. Agamy, FORTRAN IV functions for calculating probabilities associated with Dunnett's test, {\em Behavior Research Methods \& Instrumentation }, {\bf 13(3)}, (1981) 363-366. \bibitem{DUNNETT} C. Dunnett, Algorithm AS 251: Multivariate normal probabilities integrals with product correlation structure, {\em Appl. Statist.}, {\bf 38}, (1989), 564-579. \bibitem{DUNNET} C. W. Dunnett; New Tables for Multiple Comparisons with a Control {\em Biometrika}, {\bf 20(3)}, (Sep., 1964), p. 482-491, 1964. \bibitem{GENZ} A. Genz; F. Bretz, Numerical computation of the multivariate t-probabilities with applications to power calculation of multiple contrasts, {\em Journal of Statistical Computation and Simulation }, {\bf 63(4)}, (1999), 361-378. \bibitem{HOC} Y. Hochberg; A. C. Tamhane, ``Multiple Comparisons Procedures'', John Wiley \& Sons, Canadá, 1987. \bibitem{HSU} J. C. Hsu, ``Multiple Comparisons - Theory and methods'', Chapman \& Hall, USA, 1999. \bibitem{KHU} A. I. Khuri, ``Advanced Calculus with Applications in Statistics'', Second Edition, John Wiley \& Sons, USA, 2003. \bibitem{R.2.13} R Development Core Team. R: a languague and environment for statistical computing. Vienna: R Foundation for Statistical Computing, 2011. Disponível em http://www.R-project.org > Acesso em: 20 mar. 2011. \bibitem{SCH} M. J. Schervish, Algorithm AS 195: Multivariate normal probabilities with error bound. {\em Appl. Statist.}, 33, (1984), 81-94. \end{thebibliography} \end{document} \newpage $ \ \ $ \thispagestyle{myheadings} \markboth{ }{ }