Dia ou outro me vejo com a necessidade de aplicar o teste de Bartlett (função bartllet.test(), do pacote stats).

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:

bartlett.test(list(x = lm(...), z = lm(...), w = lm(...), y = lm(...))), De acordo com a sintaxe apresentada no help() da mesma (http://stat.ethz.ch/R-manual/R-patched/library/stats/html/bartlett.test.html).

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 (http://www.jstor.org/stable/96803). E a fonte da função disponível no tar.gz do R 2.15.2, que baixei do CRAN (http://cran.r-project.org/src/base/R-2/R-2.15.2.tar.gz, acesso hoje 10 de dezembro de 2012).

Transcrevo no final da mensagem um trecho da fonte, que se encontra em: ~/R-2.15.2/src/library/stats/R/bartlett.test.R.

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!

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!

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?

Agradeço qualquer explicação, sugestão, etc.

bartlett.test.default <-

function(x, g, ...)

{

LM <- FALSE

if (is.list(x)) {

if (length(x) < 2L)

stop("'x' must be a list with at least 2 elements")

DNAME <- deparse(substitute(x))

if (all(sapply(x, function(obj) inherits(obj, "lm"))))

LM <- TRUE

else

x <- lapply(x, function(x) x <- x[is.finite(x)])

k <- length(x)

}

else {

if (length(x) != length(g))

stop("'x' and 'g' must have the same length")

DNAME <- paste(deparse(substitute(x)), "and",

deparse(substitute(g)))

OK <- complete.cases(x, g)

x <- x[OK]

g <- factor(g[OK])

k <- nlevels(g)

if (k < 2)

stop("all observations are in the same group")

x <- split(x, g)

}

if (LM) {

n <- sapply(x, function(obj) obj$df.resid)

v <- sapply(x, function(obj) sum(obj$residuals^2)) ***** DESTAQUE ***** DESTAQUE ***** DESTAQUE *****

} else {

n <- sapply(x, "length") - 1

if (any(n <= 0))

stop("there must be at least 2 observations in each group")

v <- sapply(x, "var")

}

n.total <- sum(n)

v.total <- sum(n * v) / n.total

STATISTIC <- ((n.total * log(v.total) - sum(n * log(v))) /

(1 + (sum(1 / n) - 1 / n.total) / (3 * (k - 1))))

PARAMETER <- k - 1

PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)

names(STATISTIC) <- "Bartlett's K-squared"

names(PARAMETER) <- "df"

RVAL <- list(statistic = STATISTIC,

parameter = PARAMETER,

p.value = PVAL,

data.name = DNAME,

method = "Bartlett test of homogeneity of variances")

class(RVAL) <- "htest"

return(RVAL)

}