【问题标题】:Convert part of a statistical function's output into a data.frame将统计函数的部分输出转换为 data.frame
【发布时间】:2021-11-26 05:47:04
【问题描述】:

我想知道是否有办法将resres2 对象的输出的以下部分转换为data.frame?

注意:下面的答案适用于res,但不适用于res2

感谢您提供实用的答案,因为以下数据只是玩具。

library(metafor)

dat <- dat.konstantopoulos2011
res <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat)

#== OUTPUT (CAN WE TURN ONLY BELOW PART INTO A data.frame?):

#Variance Components:

#            estim    sqrt  nlvls  fixed           factor 
#sigma^2.1  0.0651  0.2551     11     no         district 
#sigma^2.2  0.0327  0.1809     56     no  district/school 

#Test for Heterogeneity:
#Q(df = 55) = 578.8640, p-val < .0001

# AND

res2 <- rma.mv(yi, vi, random = ~ factor(school) | district, data=dat)

#== OUTPUT (CAN WE TURN ONLY BELOW PART INTO A data.frame?):
#Variance Components:

#outer factor: district       (nlvls = 11)
#inner factor: factor(school) (nlvls = 11)

#            estim    sqrt  fixed 
#tau^2      0.0978  0.3127     no 
#rho        0.6653             no 

#Test for Heterogeneity:
#Q(df = 55) = 578.8640, p-val < .0001

【问题讨论】:

  • 你可以检查它是否兼容扫帚包
  • @ViníciusFélix,如果您可以发布解决方案,我很乐意尝试一下。

标签: r dataframe function dplyr tidyverse


【解决方案1】:

如果没有提取数据的默认/标准方法,那么您可以使用capture.output 操作输出。

return_data <- function(res) {
  tmp <- capture.output(res)
  #data start from second line after "Variance Components:"
  start <- which(tmp == "Variance Components:") + 2
  index <- which(tmp == "") 
  #Data ends before the empty line after "Variance Components:"
  end <- index[which.max(index > start)] - 1
  data <- read.table(text = paste0(tmp[start:end], collapse = '\n'), header = T)
  heterogeneity_index <- which(tmp == "Test for Heterogeneity:") + 1
  list(data = data, heterogeneity = tmp[heterogeneity_index])
}

res <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat)
return_data(res)

#$data
#           estim   sqrt nlvls fixed          factor
#sigma^2.1 0.0651 0.2551    11    no        district
#sigma^2.2 0.0327 0.1809    56    no district/school

#$heterogeneity
#[1] "Q(df = 55) = 578.8640, p-val < .0001"

【讨论】:

  • 谢谢!我们能不能把它变成一个函数来与其他模型一起工作,比如:res2 &lt;- rma.mv(yi, vi, random = ~ factor(school) | district, data=dat)
  • 我看到resres2 完全不同。不能真正概括它们并将它们放在一个函数中。
  • 好的,+2-1 是什么?如果我们假设所有对象都像res,但输出的行数更多,我们可以把它变成一个函数吗?
  • 我在答案中添加了一些 cmets 并使其成为一个函数。
  • 在问题中添加了res2。我还为res2 使用了+2,但可以正确捕获方差分量部分。有什么办法解决这个问题吗?
【解决方案2】:

这符合您的目的吗? 'Test for Heterogeneity' 并不真正适合数据框,所以我将它添加为一个单独的列,结果它被重复了。我不知道你还能怎么做。

library(tidyverse)
#install.packages("metafor")
library(metafor)
#> Loading required package: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack
#> 
#> Loading the 'metafor' package (version 3.0-2). For an
#> introduction to the package please type: help(metafor)

dat <- dat.konstantopoulos2011
res <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat)
res
#> 
#> Multivariate Meta-Analysis Model (k = 56; method: REML)
#> 
#> Variance Components:
#> 
#>             estim    sqrt  nlvls  fixed           factor 
#> sigma^2.1  0.0651  0.2551     11     no         district 
#> sigma^2.2  0.0327  0.1809     56     no  district/school 
#> 
#> Test for Heterogeneity:
#> Q(df = 55) = 578.8640, p-val < .0001
#> 
#> Model Results:
#> 
#> estimate      se    zval    pval   ci.lb   ci.ub 
#>   0.1847  0.0846  2.1845  0.0289  0.0190  0.3504  * 
#> 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

vc <- cbind(estim = res$sigma2,
            sqrt = res$sigma,
            nlvls = res$s.nlevels, 
            fixed = ifelse(res$vc.fix$sigma2, "yes", "no"), 
            factor = res$s.names,
            R = ifelse(res$Rfix, "yes", "no"),
            Test_for_heterogeneity = paste0("Q(df = ", res$k - res$p, ") = ", metafor:::.fcf(res$QE, res$digits[["test"]]), ", p-val ", metafor:::.pval(res$QEp, 
                                                                                                                               res$digits[["pval"]], showeq = TRUE, sep = " "))
            
)
rownames(vc) <- c("sigma^2.1", "sigma^2.2")
result <- as.data.frame(vc)
result
#>      estim                nlvls fixed factor            R    Test_for_heterogeneity                
#> sigma^2.1 "0.0650619442753117" "11"  "no"  "district"        "no" "Q(df = 55) = 578.8640, p-val < .0001"
#> sigma^2.2 "0.0327365170279351" "56"  "no"  "district/school" "no" "Q(df = 55) = 578.8640, p-val < .0001"

reprex package (v2.0.1) 于 2021-10-06 创建

【讨论】:

  • 非常感谢您的时间和精力。但正如我在帖子中指出的那样,我正在寻找一种功能性(并且可能有效)的解决方案。我的一部分说capture.output() 是解决方案。
  • 我认为@ronak-shah 的解决方案正是您想要的 :)
猜你喜欢
  • 1970-01-01
  • 2011-02-06
  • 1970-01-01
  • 2013-04-23
  • 2016-01-13
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多