Блог

Я использую функцию пакета Anova ‘car’ для некоторого статистического тестирования.

Это дает следующий результат:

     Y = cbind(curdata$V1, curdata$V2, curdata$V3)
    mymdl = lm(Y ~ curdata$V4   curdata$V5)
    myanova = Anova(mymdl)
    
Type II MANOVA Tests: Pillai test statistic
           Df test stat approx F num Df den Df  Pr(>F)  
curdata$V4  1   0.27941   2.9728      3     23 0.05280 .
curdata$V5  1   0.33570   3.8743      3     23 0.02228 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 

Я хотел бы извлечь значения из столбца ‘Pr (> F)’, чтобы я мог поместить эти p-значения в другую матрицу для последующей коррекции множественных сравнений.

Я пробовал использовать unlist, но он по-прежнему не предоставляет p-значения, найденные в столбце.

Любая помощь в этом была бы весьма признательна.

Если у нас есть несколько переменных ответа, это a Manova . Мы могли бы захватить выходные данные и использовать регулярное выражение

 as.numeric(sub(".*\s*(\d \.[0-9e-] )\s*[*.]*", "\1", capture.output(out)[4:5]))
#[1] 8.836e-06 2.200e-16
 

данные

  mymdl <- lm(cbind(Sepal.Length, Sepal.Width) ~ Petal.Width   
     Petal.Length, data = iris)

 out <- Anova(mymdl)
 

Возможно, это не самый практичный способ, но вы можете поиграть со столбцами, используя separate() from tidyr :

 library(car)
library(dplyr)
library(tidyr)
#Code
v1 <- data.frame(capture.output(myanova))
v1 <- v1[3:5,,drop=F]
names(v1)<-'v1'
v2 <- separate(v1,v1,c(paste0('v',1:21)),sep = '\s')
v2 <- v2[-1,]
 

Вывод:

 as.numeric(v2$v21)
[1] 8.836e-06 2.200e-16
 

Предупреждение: при необходимости вам потребуется изменить 1:21 , если в действии захвата присутствует больше столбцов.

TLDR:

 # define helper:
get_summary_for_print <- car:::print.Anova.mlm
body(get_summary_for_print) <- local({tmp <- body(get_summary_for_print);tmp[-(length(tmp)-(0:1))]})
#use it:
get_summary_for_print(Anova(mymdl))$`Pr(>F)`
 

К сожалению, нет определенного способа. Но вы можете посмотреть источник car:::print.Anova.mlm (введя это в R консоли), чтобы узнать, как он получает нужные вам значения:

 function (x, ...) 
{
    if ((!is.null(x$singular)) amp;amp; x$singular) 
        stop("singular error SSP matrix; multivariate tests unavailablentry summary(object, multivariate=FALSE)")
    test <- x$test
    repeated <- x$repeated
    ntests <- length(x$terms)
    tests <- matrix(NA, ntests, 4)
    if (!repeated) 
        SSPE.qr <- qr(x$SSPE)
    for (term in 1:ntests) {
        eigs <- Re(eigen(qr.coef(if (repeated) qr(x$SSPE[[term]]) else SSPE.qr, 
            x$SSP[[term]]), symmetric = FALSE)$values)
        tests[term, 1:4] <- switch(test, Pillai = Pillai(eigs, 
            x$df[term], x$error.df), Wilks = Wilks(eigs, x$df[term], 
            x$error.df), `Hotelling-Lawley` = HL(eigs, x$df[term], 
            x$error.df), Roy = Roy(eigs, x$df[term], x$error.df))
    }
    ok <- tests[, 2] >= 0 amp; tests[, 3] > 0 amp; tests[, 4] > 0
    ok <- !is.na(ok) amp; ok
    tests <- cbind(x$df, tests, pf(tests[ok, 2], tests[ok, 3], 
        tests[ok, 4], lower.tail = FALSE))
    rownames(tests) <- x$terms
    colnames(tests) <- c("Df", "test stat", "approx F", "num Df", 
        "den Df", "Pr(>F)")
    tests <- structure(as.data.frame(tests), heading = paste("nType ", 
        x$type, if (repeated) 
            " Repeated Measures", " MANOVA Tests: ", test, " test statistic", 
        sep = ""), class = c("anova", "data.frame"))
    print(tests, ...)
    invisible(x)
}
<bytecode: 0x56032ea80990>
<environment: namespace:car>
 

В этом случае для вычисления p-значений требуется довольно много строк кода. Однако мы можем легко создать модифицированную версию print функции для возврата таблицы ( tests ) вместо того, чтобы только печатать ее ( print(tests, ...) ) и возвращать исходный объект ( invisible(x) ):

 get_summary_for_print <- car:::print.Anova.mlm # copy the original print function (inclusive environment)
body(get_summary_for_print) <- # replace the code of our copy
    local({ # to avoid pollution of environment by tmp
        tmp <- body(get_summary_for_print) # to avoid code duplication
        tmp[-(length(tmp)-(0:1))] # remove the last two code lines of the function
    })
 

И используйте его, например, так:

 library(car)
#> Loading required package: carData
res <- Anova(lm(cbind(Sepal.Width, Sepal.Length, Petal.Width) ~ Species   Petal.Length, iris))
res
#> 
#> Type II MANOVA Tests: Pillai test statistic
#>              Df test stat approx F num Df den Df    Pr(>F)    
#> Species       2   0.70215   26.149      6    290 < 2.2e-16 ***
#> Petal.Length  1   0.63487   83.461      3    144 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
str(get_summary_for_print(res))
#> Classes 'anova' and 'data.frame':    2 obs. of  6 variables:
#>  $ Df       : num  2 1
#>  $ test stat: num  0.702 0.635
#>  $ approx F : num  26.1 83.5
#>  $ num Df   : num  6 3
#>  $ den Df   : num  290 144
#>  $ Pr(>F)   : num  7.96e-25 2.41e-31
#>  - attr(*, "heading")= chr "nType II MANOVA Tests: Pillai test statistic"