Senhores,<div><br></div><div>Dia ou outro me vejo com a necessidade de aplicar o teste de Bartlett (função bartllet.test(), do pacote stats).</div><div><br></div><div>Em um desses usos regulares notei certo "comportamento" anômalo na mesma. O uso que dou a ela era para testar a homogeneidade de variâncias entre modelos ajustados, ou seja, aplicava-a sobre uma lista de saídas da função lm(), da forma:</div>
<div><br></div><div>bartlett.test(list(x = lm(...), z = lm(...), w = lm(...), y = lm(...))), De acordo com a sintaxe apresentada no help() da mesma (<a href="http://stat.ethz.ch/R-manual/R-patched/library/stats/html/bartlett.test.html">http://stat.ethz.ch/R-manual/R-patched/library/stats/html/bartlett.test.html</a>).</div>
<div><br></div><div>Para investigar o porque desse comportamento anômalo, além do help() indicado no link, busquei a citação do help, de autoria próprio Bartlett (<a href="http://www.jstor.org/stable/96803">http://www.jstor.org/stable/96803</a>). E a fonte da função disponível no tar.gz do R 2.15.2, que baixei do CRAN (<a href="http://cran.r-project.org/src/base/R-2/R-2.15.2.tar.gz">http://cran.r-project.org/src/base/R-2/R-2.15.2.tar.gz</a>, acesso hoje 10 de dezembro de 2012).</div>
<div><br></div><div>Transcrevo no final da mensagem um trecho da fonte, que se encontra em: ~/R-2.15.2/src/library/stats/R/bartlett.test.R.</div><div><br></div><div>Pois bem, lendo a código fonte, o artigo do Bartlett e consultando o Livro Steel, Torrie and Dickey "Principles and Procedures of Statistics: A biometrial Approach", terceira edição, capítulo 19, página 480. Pude perceber que todo o fundamento do teste aplica-se a variâncias. Entretanto ao observar a implementação da função (linha destacada) observei que ao "capturar" sum(obj$residuals^2), a função estaria atendo-se a soma de quadrados do resíduo, e não a variância. Obviamente, uma soma de quadrados, dividido pelo seu respectivo grau de liberdade é uma quadrado médio e por sua vez uma variância, entretanto, isso não é realizado!</div>
<div><br></div><div>Pelos testes que fiz, ao utilizar a função com diferentes modelos ajustados com o mesmo grau de liberdade dos resíduos, pela proporcionalidade das somas de quadrados o resultado do teste é satisfatório, entretanto, para diferentes modelos, com diferentes graus de liberdades o resultado do teste não confere com o obtido calculado a mão!</div>
<div><br></div><div>Diante do exposto, gostaria de ouvir a opinião dos senhores quanto a esse fato, estou certo? Em caso positivo, qual seria um meio elegante de entrar em contato com o R CORE TEAM no sentido de apresentar esse fato?</div>
<div><br></div><div>Agradeço qualquer explicação, sugestão, etc.</div><div><br></div><div>att,</div><div>FH</div><div><br></div><div><div>bartlett.test.default <-</div><div>function(x, g, ...)</div><div>{</div><div> LM <- FALSE</div>
<div> if (is.list(x)) {</div><div> if (length(x) < 2L)</div><div> stop("'x' must be a list with at least 2 elements")</div><div> DNAME <- deparse(substitute(x))</div><div>
if (all(sapply(x, function(obj) inherits(obj, "lm"))))</div><div> LM <- TRUE</div><div> else</div><div> x <- lapply(x, function(x) x <- x[is.finite(x)])</div><div> k <- length(x)</div>
<div> }</div><div> else {</div><div> if (length(x) != length(g))</div><div> stop("'x' and 'g' must have the same length")</div><div> DNAME <- paste(deparse(substitute(x)), "and",</div>
<div> deparse(substitute(g)))</div><div> OK <- complete.cases(x, g)</div><div> x <- x[OK]</div><div> g <- factor(g[OK])</div><div> k <- nlevels(g)</div><div> if (k < 2)</div>
<div> stop("all observations are in the same group")</div><div> x <- split(x, g)</div><div> }</div><div><br></div><div> if (LM) {</div><div> n <- sapply(x, function(obj) obj$df.resid)</div>
<div> v <- sapply(x, function(obj) sum(obj$residuals^2)) ***** DESTAQUE ***** DESTAQUE ***** DESTAQUE *****</div><div> } else {</div><div> n <- sapply(x, "length") - 1</div><div> if (any(n <= 0))</div>
<div> stop("there must be at least 2 observations in each group")</div><div> v <- sapply(x, "var")</div><div> }</div><div><br></div><div> n.total <- sum(n)</div><div> v.total <- sum(n * v) / n.total</div>
<div> STATISTIC <- ((n.total * log(v.total) - sum(n * log(v))) /</div><div> (1 + (sum(1 / n) - 1 / n.total) / (3 * (k - 1))))</div><div> PARAMETER <- k - 1</div><div> PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)</div>
<div> names(STATISTIC) <- "Bartlett's K-squared"</div><div> names(PARAMETER) <- "df"</div><div><br></div><div> RVAL <- list(statistic = STATISTIC,</div><div> parameter = PARAMETER,</div>
<div> p.value = PVAL,</div><div> <a href="http://data.name">data.name</a> = DNAME,</div><div> method = "Bartlett test of homogeneity of variances")</div><div> class(RVAL) <- "htest"</div>
<div> return(RVAL)</div><div>}</div></div>