Algoritmo utilizando quadraturas gaussianas para a obtenção das probabilidades do teste bilateral de Dunnett para dados balanceados

Algoritmo utilizando quadraturas gaussianas para a obtenção das probabilidades do teste bilateral de Dunnett para dados balanceados^*

S.C. Broch^I,**; D.F. Ferreira^II

^IProfessora do Instituto Federal Farroupilha, IFFarroupilha, 98130-000 Santa Maria, RS, Brasil

^IIProfessor associado 4, Bolsista CNPq, Departamento de Ciências Exatas, DEX, UFLA - Universidade Federal de Lavras, 37200-000 Lavras, MG, Brasil. E-mail: [e-mail]

RESUMO

O teste de Dunnett é um teste de comparações múltiplas em que se confrontam as médias de r tratamentos em teste com a média de um tratamento controle ou testemunha. Este teste controla simultaneamente a taxa de erro tipo I por experimento num nível de significância nominal α. A limitação para o uso do teste de Dunnett é a dificuldade de obter as probabilidades da distribuição e os valores dos quantis das estatísticas do teste, pois as correlações possíveis entre os tratamentos têm larga amplitude. Neste artigo é apresentado um algoritmo que utiliza o método numérico das quadratura gaussiana para resolver as integrais e obter as probabilidades do teste bilateral de Dunnett balanceado. O algoritmo mostrou resultados precisos quando comparados com valores obtidos na literatura e em relação aos valores obtidos nos três programas analisados.

Palavras-chave: Procedimentos de comparações múltiplas, integração numérica, estatística experimental.

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 α. 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.

Keywords: multiple comparison procedures, numerical integration, experimental statistic.

1 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 em teste 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 num valor específico α [6]. A estatística deste teste também e utilizada para determinar os intervalos de confiança dos verdadeiros valores das diferenças entre a média de cada um dos tratamentos em teste e o tratamento controle, com um valor 1 - α de coeficiente de confiança conjunta.

O teste de Dunnett pode ser realizado em situações de experimentos balanceados, em que todos os tratamentostêmo mesmo número derepetições ou em que todos os tratamento sem teste 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 e 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 tabela dos valores para os casos de correlações diferentes de 0, 5.

Neste artigo é apresentado um algoritmo que utiliza quadratura gaussianas para resolver as integrais e obter as probabilidades do teste bilateral de Dunnett balanceado. A seção 2 apresenta aspectos teóricos básicos da estatística do teste bilateral de Dunnett para dados balanceados. Também são citados os principais algoritmos encontrados na literatura que buscaram resolver o problema em estudo. A seção 3 apresenta a teoria das quadraturas gaussianas, que é a ferramenta matemática numérica utilizada na resolução das integrais envolvidas na estatística do teste. A seção 4 apresenta o pseudo-código do algoritmo implementado em ambiente R. Posteriormente são mostrados resultados obtidos com o algoritmo desenvolvido em comparação com outros algoritmos existentes na literatura (implementados por nós no R), e outros já disponíveis no próprio R.

2 TESTE DE DUNNETT

Considerando Y_j1, ..., Y_jnj 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, ..., n_j.Da mesma forma, os Y_ja's são independentemente distribuídos para os diferentes valores de j = 1, 2, ..., r + 1. Assumindo normalidade e homogeneidade de variâncias da amostra aleatória, tem-se o modelo

Y_ja = µ_j + є_ja, j = 1, ..., r + 1, a = 1, ..., n_j,

em que µ_j média do j-ésimo tratamento , e є₁₁,...,є são independentes e identicamentes distribuídos como uma normal com média 0 e variância σ² desconhecida.

Usando a notação

para a média amostral do j-ésimo tratamento, e

para a variância residual amostral combinada com ν graus de liberdade.

A distribuição da estatística para comparar r médias com a testemunha, considera ν 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 _ρj , j = 1, 2, ..., 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 ρ_j = , em que é a variância da média do tratamento controle e é 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 σ². Assim, o coeficiente de correlação é ρ_j =n_j /(n_j + n_r+1).

Quando os experimentos são não balanceados, n_r+1 ≠ n_j e tem-se ρ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, ρ = 0, 5; se n_j =n para todo j mas n≠ n_r+1, tem-se o caso com apenas um parâmetro de correlação dado por ρ = 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|, segundo Hsu [6], como a solução da equação:

sendo

•

a função densidade da normal padrão;

•

a função de distribuição da normal padrão

• a função densidade da variável aleatória S definida como S =

2.1 Algoritmos propostos na literatura

Dunlap et al. (1981)[1] 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 e Duffy (1975 apud [1]) ou, para valores extremos, pela aproximação de Zelen e Severo (1965 apud [1]). 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 e Hendrickson (1971 apud [1]) era da ordem da quarta casa decimal.

Schervish [9] 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 e Schervish (1981) (apud [9]). 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 correlação deve ser não singular e passível de ser escrita como uma matriz triangular inferior.

Dunnett publica em 1989 [2] 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 [2])(Hsu, 1984 apud [2]). O pacote "mvtnormpcs" do software 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 e Bretz [4] apresentam um novo método para calcular probabilidades da distribuição t multivariada a partir de uma integral q-variada transformando-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 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. [1] 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.

3 QUADRATURAS GAUSSIANAS

A quadratura gaussiana consiste em tomar x₁< x₂ < x₃ < ... < x_n , um conjunto de n pontos distintos no intervalo [a, b], de modo que a aproximação

cujos coeficientes w₁, w₂, ..., w_n são determinados sob os pontos x₁, x₂, ..., x_n, minimize o erro esperado no cálculo da integral (adaptado de Khuri [7]). São escolhidos os w_i's e os x_i's tais que a aproximação seja exata quando a função a ser integrada g(x) é um polinômio de grau menor ou igual a 2n -1. Nos casos de polinômios de grau maior ou igual a 2n ou quando g(x) é outra função não polinomial, tem-se um erro de aproximação dado por

Algumas extensões dessa aproximação são obtidas tomando-se ƒ(x) = g(x)/λ(x) como uma nova função. Assim

em que λ(x) é uma particular função positiva denominada de função peso. A escolha dos x_i's depende da forma de λ(x). Os valores dos x_i's são as raízes do polinômio de grau n, pertencente a sequência de polinômios em [a, b] em relação a λ(x) [7]. A finalidade de uma função peso é atribuir graus variados de importância a aproximações em certas partes do intervalo.

Na Tabela 1 estão sintetizados alguns dos polinômios ortogonais clássicos, suas respectivas função peso, intervalo e simbologia usualmente utilizada para representá-los. Cada função peso e seu polinômio ortogonal dá origem a um tipo de quadratura gaussiana, cujo nome está relacionado com o polinômio interpolador.

Para resolver a integral (2.1) utilizando quadraturas, optou-se por não utilizar a quadratura Gauss-Hermite e transformar a integral original numa integral com intervalo [-1, 1] eavaliá-la por uma quadratura de Gauss-Legendre. Esta decisão se deveu principalmente por dois fatores:

• na implementação das quadraturas Gauss-Hermite e Gauss-Laguerre surgem grandes valores (ou valores muito pequenos) dos nós das respectivas quadraturas. Nos computadores, a precisão é finita. Assim, determinadas funções apresentam imagens que computacionalmente dão resultados infinitos ou são arredondados erroneamente para zero, quando computadas para argumentos muito grandes ou muito pequenos. Isso decorre da precisão finita dos números reais representados por máquinas. Esses erros são denotados por overflow ou underflow. Determinados casos das distribuições da estatística do teste de Dunnett, quando submetidos a esses tipos de quadraturas, apresentaram problemas dessa natureza, o que conduziu a uma menor acurácia da integração, quando comparada com a quadratura Gauss-Legendre.

• as funções que serão integradas não são definidas para x < 0, o que pode acarretar dificuldade para avaliar estas integrais em intervalos infinitos.

Para aplicar a quadratura Gauss-Legendre, em que λ(x) = 1, quando os limites de integração em não são a =-1 e b = 1, pode-se converter a integral para estes limites fazendo uma transformação da forma

em que os z_i's são as raízes dos polinômios de Legendre de grau n.

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, ∞) pode ser transformado inicialmente para um intervalo finito [0, 1], dado por

e convertendo esta nova integral para os limites [-1, 1], obtendo-se

Para resolver a expressão (2.1) utilizando quadratura gaussiana a integral interna original foi transformada numa integral com intervalo [-1, 1] e avaliada por uma quadratura de Gauss-Legendre. Para resolver a integral externa optou-se por não utilizar a quadratura Gauss-Laguerre mas dividir a integral original na soma de duas integrais, uma com intervalo [0, 1] e a outra com intervalo [1, ∞]. Ambas as integrais sofreram mudanças de variáveis para serem resolvida no intervalo [-1, 1] por uma quadratura de Gauss-Legendre.

ALGORITMO - PSEUDOCÓDIGO

O algoritmo descrito calcula a probabilidade do teste bilateral de Dunnett por meio da expressão (2.1) para o quantil desejado. Ele foi implementadas no software R por ser um programa de livre distribuição e código fonte aberto e utiliza algumas rotinas já existentes no programa.

No pacote "statmod " encontra-se a função "gauss.quad" que possibilita a obtenção dos nós (x_i's) e pesos (w_i's) da quadratura desejada fornecendo como argumentos:

• o número de pontos que se deseja realizar a quadratura (n);

• o tipo de quadratura (kind =), que pode ser: "legendre"; "chebyshev1"; "chebyshev2"; "hermite" que neste programa ómios de Hermite físicos, ou seja, com λ (x) = e^-x; "jacobi" que neste caso e necessário α e de β:e "laguerre" cujo padrão é α = 0 mas pode ser especificado outro valor de α, contanto que seja maior do que -1;

• valor de α que por padrão é zero, nesse caso não sendo necessário especificar;

• valor de β que por padrão é zero, nesse caso não sendo necessário especificar.

Outras rotinas utilizadas para calcular as principais distribuições de probabilidades foram:

• "pnorm" possibilita a obtenção do valor da função de distribuição da normal padrão fornecendo o quantil desejado;

• "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.

A.1 Função auxiliar para computar os valores da função da integral interna da expressão (2.1) (g x DunI n fty (x, q, c, r))

1. Recebe x: real (peso da quadratura Gauss-Hermite), q =|d|s[i], c e r;

2. calcular

3. retornar ƒ

x

.

A. 2 Quadratura Gauss-Hermite para obtenção do valor da integral interna da expressão (2.1) (GH Dun(q, c, r, n))

1. Recebe q =|d|s[i], c, r e n;

2. calcular os nós e os pesos da quadratura Gauss-Hermite por xalpha =gauss.quad(n, kind = "hermite") que corresponde aos vetores xalpha$nodes(n ×1) dos nós e xalpha$weights(n ×1) dos pesos;

3. chamar a função gxDunInfty e criar um vetor denotado por dun(n ×1) e calcular o seu i-ésimo componente por dun[i] =gxDunI n f t y(xalpha$nodes[i]), q, c, r) com i = 1, 2, ..., n;

4. calcular aux = Σdun[i] × xalpha$weights[i];

5. retornar aux.

A.3 Função original da integral interna da expressão (2.1) ( f 0dunnet t(s, |d|, c, r, ν))

1. Recebe s: vetor de pesos da quadratura Gauss-Legendre de dimensão (n ×1), |d|, c, r e ν;

2. calcular o vetor gsdf(n ×1), sendo seu i-ésimo elemento dado por

com i =1, 2, ..., n;

3. chamar a função GH Dun e calcular o vetor gs(n ×1) por gs[i] =GH Dun(s[i]× |d|, c, r, 32) com i =1, 2, ..., n;

4. obter o vetor auxiliar aux(n ×1)por aux[i] =gs[i] × g s d f[i]com i =1, 2, ..., n;

5. retornar aux.

A.4 Função para transformar a escala de 0 a 1 para 0 a ∞ numericamente (ƒ1 dunnett (s, |d|, c, r, ν))

1. Recebe s(n ×1), |d|, c, r e ν;

2. calcular o vetor trans(n ×1), sendo seu i-ésimo elemento dado por trans[i]=-ln(s[i]), com i =1, 2, ..., n;

3. chamar a função ƒ0dunnett e calcular os elementos do vetor f1(n ×1) por f 1[i]= × f 0dunnett(trans[i], |d|, c, r, ν) com i =1, 2, ..., n;

4. retornar ƒ1.

A.5 Função para transformar a escala de -1 a 1 para 0 a 1 numericamente (dunnet t(s, |d|, c, r, ν))

1. Recebe s(n ×1), |d|, c, r e ν;

2. calcular o vetor trans(n×1), sendo seu i-ésimo elemento dado por trans[i]= ×s[i]+ , com i =1, 2, ..., n;

3. chamar a função f 1 dunnett e calcular os elementos do vetor f2 (n ×1)por f 2[i]= × f 1dunnett(trans[i], |d|, c, r, ν)com i =1, 2, ..., n;

4. retornar ƒ2.

A.6 Função principal para calcular a expressão (2.1) (GLdunnet t (|d|, c, r, ν , n))

1. Recebe |d|, c, r, ν e n;

2. se (df > 1000)&(n < 200) faça n = 1000 senão se (df > 10000)&(n < 500) faça n =1000;

3. se df =∞, então retorne GH Dun(|d|, r, c, 32);senão váparao passo 4;

4. 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");

5. chamar a função dunnett para cada nóe calcular o vetor dunn, sendo seu i-ésimo elemento dado por

dunn[i] =dunnett(xalpha$nodes[i], |d|, r, c, ν),

com i =1, 2, ..., n;

6. calcular aux =dunn[i] × xalpha$weights[i] que é a probabilidade desejada.

5 RESULTADOS

Para aplicar o teste de Dunnett deve-se seguir os seguintes passos:

1. Para cada contraste determinar a estimativa

2. entrar com as informações como argumentos do algoritmo implementado: o valor 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 ν e o número de pontos n para a quadratura;

3. a rotina fornece a probabilidade conjunta estimada 1 - p, da qual se obtém o valor-p do teste;

4. compara-se o valor-p obtido no contraste com o nível de significância nominal α do teste;

5. conclui-se sobre a significância da diferença de médias observada da seguinte forma: se o valor-p for inferior ao α, 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.

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 2 estão alguns valores de probabilidades do teste bilateral de Dunnett balanceado, calculados pelo presente algoritmo, em função de |d|, ρ, ν e k, e também daqueles obtidos em tabelas fornecidas em Hsu [6], calculados pelo algoritmo de Dunlap et al. [1], pelos pacotes "mvtnormpcs" [2] e "mvtnorm" [4] do software R.

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 [6] 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 ao 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.

6 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.

Recebido em 28 setembro, 2012

Aceito em 16 agosto, 2013

• [1] W.P. Dunlap, M.S. Marx & G.J. Agamy. FORTRAN IV functions for calculating probabilities associated with Dunnett's test. Behavior Research Methods & Instrumentation, 13(3)(1981),363-366.
• [2] C. Dunnett. Algorithm AS 251: Multivariate normal probabilities integrals with product correlation structure. Appl. Statist., 38(1989),564-579.
• [3] C.W. Dunnett. New Tables for Multiple Comparisons with a Control. Biometrika, 20(3)(1964),482-491, Sept.
• [4] A. Genz & F. Bretz. Numerical computation of the multivariate t-probabilities with applications to power calculation of multiple contrasts. Journal of Statistical Computation and Simulation, 63(4)(1999),361-378.
• [5] Y. Hochberg & A.C. Tamhane. "Multiple Comparisons Procedures", John Wiley & Sons, Canadá, (1987).
• [6] J.C. Hsu. "Multiple Comparisons - Theory and methods", Chapman & Hall, USA, (1999).
• [7] A.I. Khuri. "Advanced Calculus with Applications in Statistics", Second Edition, John Wiley & Sons, USA, (2003).
• [8] 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.
• [9] M.J. Schervish. Algorithm AS 195: Multivariate normal probabilities with error bound. Appl. Statist., 33(1984),81-94.

*

Trabalho apresentado no XXXIV CNMAC (2012) com apoio da Capes e Fapemig.

Datas de Publicação

• Publicação nesta coleção
  11 Out 2013
• Data do Fascículo
  Ago 2013

Histórico

• Recebido
  28 Set 2012
• Aceito
  16 Ago 2013