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) }
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
