[R-br] Implementação Função

[FHRB Toledo]

Senhores,

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.

att,
FH

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)
}
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20121210/b2f06111/attachment.html>