Rückgabevarianz aus Bartlett-Test der Varianzhomogenität in R

Question

Bartlett-Test ermöglicht es Ihnen zu testen, ob die Varianz in verschiedenen Gruppen gleich ist oder nicht.

Das stats Paket in R hat die Funktion bartlett.test. Dies ist ein Beispiel eines Datensatzes verwenden, die in R.

bt <- bartlett.test(count ~ spray, data = InsectSprays) 

verfügbar ist, wie Sie die tatsächliche Abweichung von bartlett.test bekommen Sie?

Ich kann nicht das scheint in dem Objekt zu finden bt

names(bt) 
[1] "data.name" "method" "p.value" "parameter" "statistic" 

Sie können die Varianz berechnen sich var() verwenden. Eine Möglichkeit dazu besteht in der Verwendung von summaryBy.

library(doBy) 
summaryBy(count~spray, data=InsectSprays, FUN=var) 

jedoch würden Sie bartlett.test erwarten, dass die Varianz pro Gruppe zur Verfügung zu stellen. Ähnlich berechnet die Berechnung eines T.tests in R auch den Mittelwert pro Gruppe. Also, können wir die Varianz pro Gruppe von bartlett.test in R extrahieren, und wie?

Answer 1

Nein, leider können Sie nicht.

Varianz erscheint nicht in der Struktur des zurückgegebenen Objekts.

Ich lese die Quelle der Funktion und Sie könnten es extrahieren, indem Sie das ein bisschen neu schreiben, aber das wäre viel mehr benutzerdefinierte Arbeit als die Lösung, die Sie bereits an Ort und Stelle haben.

können Sie sehen, was ich meine hier:

# File src/library/stats/R/bartlett.test.R 
# Part of the R package, https://www.R-project.org 
# 
# Copyright (C) 1995-2015 The R Core Team 
# 
# This program is free software; you can redistribute it and/or modify 
# it under the terms of the GNU General Public License as published by 
# the Free Software Foundation; either version 2 of the License, or 
# (at your option) any later version. 
# 
# This program is distributed in the hope that it will be useful, 
# but WITHOUT ANY WARRANTY; without even the implied warranty of 
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU General Public License for more details. 
# 
# A copy of the GNU General Public License is available at 
# https://www.R-project.org/Licenses/ 

bartlett.test <- function(x, ...) UseMethod("bartlett.test") 

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)) 
    } 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) 
} 

bartlett.test.formula <- 
function(formula, data, subset, na.action, ...) 
{ 
    if(missing(formula) || (length(formula) != 3L)) 
     stop("'formula' missing or incorrect") 
    m <- match.call(expand.dots = FALSE) 
    if(is.matrix(eval(m$data, parent.frame()))) 
     m$data <- as.data.frame(data) 
    ## need stats:: for non-standard evaluation 
    m[[1L]] <- quote(stats::model.frame) 
    mf <- eval(m, parent.frame()) 
    if(length(mf) != 2L) 
     stop("'formula' should be of the form response ~ group") 
    DNAME <- paste(names(mf), collapse = " by ") 
    names(mf) <- NULL 
    y <- do.call("bartlett.test", as.list(mf)) 
    y$data.name <- DNAME 
    y 
} 

Answer 2

Dies ist leicht mehr hackable als Hack-R schlägt vor, aber sie sind richtig, dass so etwas wie sapply(split(InsectSprays,spray),function(x) var(x$count)) (tun sie alle in der Basis R) könnte einfacher sein.

Die hier gezeigte Technik ist zerbrechlich, weil es auf der genauen Form der eingebauten Funktion beruht; es würde aufhören zu arbeiten, wenn es sogar kleine Änderungen an der Funktion in zukünftigen Versionen von R. Safer zu dump() die gesamte Funktion und ändern Sie es nach Ihren Wünschen, dann source() die Ergebnisse.

bb <- stats:::bartlett.test.default 
bb2 <- body(bb) 
## add a line to save the variances 
bb2[[12]] <- quote(ESTIMATE <- v) 
## add the variances to the return list 
bb2[[13]] <- quote(RVAL <- list(statistic = STATISTIC, parameter = PARAMETER, estimate = ESTIMATE, p.value = PVAL, data.name = DNAME, method = "Bartlett test of homogeneity of variances")) 
## restore the rest of the function 
bb2[14:15] <- body(bb)[13:14] 
body(bb) <- bb2 

legte es nun im stats Namespace zurück:

assignInNamespace("bartlett.test.default",bb,pos="package:stats") 

Test:

(bt <- bartlett.test(count~spray,data=InsectSprays)) 
## Bartlett test of homogeneity of variances 
## 
## data: count by spray 
## Bartlett's K-squared = 25.96, df = 5, p-value = 9.085e-05 
## sample estimates: 
##   A   B   C   D   E   F 
## 22.272727 18.242424 3.901515 6.265152 3.000000 38.606061 

Sie können die Werte über bt$estimate abrufen.

Es könnte sich lohnen, dies auf r-devel als eine Erweiterung zu bartlett.test vorzuschlagen: das einzige Gegenargument, das ich mir vorstellen kann, ist, dass die Ausgabe unhandlich wäre, wenn viele Gruppen getestet würden.