logo标签列表

从rms拟合中提取所有模型统计信息？

rstatistics
7

rms包含大量有用的统计函数。然而，我无法找到从拟合对象中提取某些适配统计数据的正确方法。考虑以下示例：

library(pacman)
p_load(rms, stringr, readr)

#fit
> (fit = rms::ols(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species, data = iris))
Linear Regression Model

 rms::ols(formula = Sepal.Length ~ Sepal.Width + Petal.Length + 
     Petal.Width + Species, data = iris)

                Model Likelihood     Discrimination    
                   Ratio Test           Indexes        
 Obs     150    LR chi2    302.96    R2       0.867    
 sigma0.3068    d.f.            5    R2 adj   0.863    
 d.f.    144    Pr(> chi2) 0.0000    g        0.882    

 Residuals

       Min        1Q    Median        3Q       Max 
 -0.794236 -0.218743  0.008987  0.202546  0.731034 


                    Coef    S.E.   t     Pr(>|t|)
 Intercept           2.1713 0.2798  7.76 <0.0001 
 Sepal.Width         0.4959 0.0861  5.76 <0.0001 
 Petal.Length        0.8292 0.0685 12.10 <0.0001 
 Petal.Width        -0.3152 0.1512 -2.08 0.0389  
 Species=versicolor -0.7236 0.2402 -3.01 0.0031  
 Species=virginica  -1.0235 0.3337 -3.07 0.0026

因此，适合的print函数打印了许多有用的内容，包括标准误差和调整后的R2。不幸的是，如果我们检查模型拟合对象，这些值似乎没有出现在任何地方。
> str(fit)
List of 19
 $ coefficients     : Named num [1:6] 2.171 0.496 0.829 -0.315 -0.724 ...
  ..- attr(*, "names")= chr [1:6] "Intercept" "Sepal.Width" "Petal.Length" "Petal.Width" ...
 $ residuals        : Named num [1:150] 0.0952 0.1432 -0.0731 -0.2894 -0.0544 ...
  ..- attr(*, "names")= chr [1:150] "1" "2" "3" "4" ...
 $ effects          : Named num [1:150] -71.5659 -1.1884 9.1884 -1.3724 -0.0587 ...
  ..- attr(*, "names")= chr [1:150] "Intercept" "Sepal.Width" "Petal.Length" "Petal.Width" ...
 $ rank             : int 6
 $ fitted.values    : Named num [1:150] 5 4.76 4.77 4.89 5.05 ...
  ..- attr(*, "names")= chr [1:150] "1" "2" "3" "4" ...
 $ assign           :List of 4
  ..$ Sepal.Width : int 2
  ..$ Petal.Length: int 3
  ..$ Petal.Width : int 4
  ..$ Species     : int [1:2] 5 6
 $ qr               :List of 5
  ..$ qr   : num [1:150, 1:6] -12.2474 0.0816 0.0816 0.0816 0.0816 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:150] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:6] "Intercept" "Sepal.Width" "Petal.Length" "Petal.Width" ...
  ..$ qraux: num [1:6] 1.08 1.02 1.11 1.02 1.02 ...
  ..$ pivot: int [1:6] 1 2 3 4 5 6
  ..$ tol  : num 1e-07
  ..$ rank : int 6
  ..- attr(*, "class")= chr "qr"
 $ df.residual      : int 144
 $ var              : num [1:6, 1:6] 0.07828 -0.02258 -0.00198 0.01589 -0.02837 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:6] "Intercept" "Sepal.Width" "Petal.Length" "Petal.Width" ...
  .. ..$ : chr [1:6] "Intercept" "Sepal.Width" "Petal.Length" "Petal.Width" ...
 $ stats            : Named num [1:6] 150 302.964 5 0.867 0.882 ...
  ..- attr(*, "names")= chr [1:6] "n" "Model L.R." "d.f." "R2" ...
 $ linear.predictors: Named num [1:150] 5 4.76 4.77 4.89 5.05 ...
  ..- attr(*, "names")= chr [1:150] "1" "2" "3" "4" ...
 $ call             : language rms::ols(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species, data = iris)
 $ terms            :Classes 'terms', 'formula'  language Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species
  .. ..- attr(*, "variables")= language list(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, Species)
  .. ..- attr(*, "factors")= int [1:5, 1:4] 0 1 0 0 0 0 0 1 0 0 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:5] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" ...
  .. .. .. ..$ : chr [1:4] "Sepal.Width" "Petal.Length" "Petal.Width" "Species"
  .. ..- attr(*, "term.labels")= chr [1:4] "Sepal.Width" "Petal.Length" "Petal.Width" "Species"
  .. ..- attr(*, "order")= int [1:4] 1 1 1 1
  .. ..- attr(*, "intercept")= num 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, Species)
  .. ..- attr(*, "dataClasses")= Named chr [1:5] "numeric" "numeric" "numeric" "numeric" ...
  .. .. ..- attr(*, "names")= chr [1:5] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" ...
  .. ..- attr(*, "formula")=Class 'formula'  language Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
 $ Design           :List of 12
  ..$ name        : chr [1:4] "Sepal.Width" "Petal.Length" "Petal.Width" "Species"
  ..$ label       : chr [1:4] "Sepal.Width" "Petal.Length" "Petal.Width" "Species"
  ..$ units       : Named chr [1:4] "" "" "" ""
  .. ..- attr(*, "names")= chr [1:4] "Sepal.Width" "Petal.Length" "Petal.Width" "Species"
  ..$ colnames    : chr [1:5] "Sepal.Width" "Petal.Length" "Petal.Width" "Species=versicolor" ...
  ..$ mmcolnames  : chr [1:5] "Sepal.Width" "Petal.Length" "Petal.Width" "Speciesversicolor" ...
  ..$ assume      : chr [1:4] "asis" "asis" "asis" "category"
  ..$ assume.code : int [1:4] 1 1 1 5
  ..$ parms       :List of 1
  .. ..$ Species: chr [1:3] "setosa" "versicolor" "virginica"
  ..$ limits      : list()
  ..$ values      : list()
  ..$ nonlinear   :List of 4
  .. ..$ Sepal.Width : logi FALSE
  .. ..$ Petal.Length: logi FALSE
  .. ..$ Petal.Width : logi FALSE
  .. ..$ Species     : logi [1:2] FALSE FALSE
  ..$ interactions: NULL
 $ non.slopes       : num 1
 $ na.action        : NULL
 $ scale.pred       : chr "Sepal.Length"
 $ fail             : logi FALSE
 $ sformula         :Class 'formula'  language Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
 - attr(*, "class")= chr [1:3] "ols" "rms" "lm"

这里有一个关于R语言帮助的7年老问题，包创建者在其中解释了如何解决以下问题：

2010年8月11日，david dav写道：

嗨， 我想提取逻辑回归的系数（估计值和标准误差），与glm中的summary(fit.glm)$coef相同。

谢谢 David

coef(fit) sqrt(diag(vcov(fit)))

但是，除非一切都是线性的，没有任何交互作用，并且因子只有两个水平，否则这些内容并不会很有帮助。

Frank

根据作者的说法，这个解决方案并不是最优的。这让人们想知道显示的值是如何计算出来的。查找代码后发现需要在未记录的软件包代码（该软件包代码在Github上）中进行搜索。也就是说，我们从print.ols()开始：

> rms:::print.ols
function (x, digits = 4, long = FALSE, coefs = TRUE, title = "Linear Regression Model", 
    ...) 
{
    latex <- prType() == "latex"
    k <- 0
    z <- list()
    if (length(zz <- x$na.action)) {
        k <- k + 1
        z[[k]] <- list(type = paste("naprint", class(zz)[1], 
            sep = "."), list(zz))
    }
    stats <- x$stats
...

阅读更多后，我们发现例如 R2 adj. 是在打印函数中计算的。
rsqa <- 1 - (1 - r2) * (n - 1) / rdf

我们还发现了一些标准误差的计算，但没有 p 值。
  se <- sqrt(diag(x$var))
  z[[k]] <- list(type='coefmatrix',
                 list(coef    = x$coefficients,
                      se      = se,
                      errordf = rdf))

所有结果都会进一步传递给prModFit()。我们可以查找并找到p值的计算等。不幸的是，print命令返回NULL，因此这些值无法在程序中重复使用：
> x = print((fit = rms::ols(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species, data = iris)))
#printed output...
> x
NULL

如何获得所有统计数据？
我不确定“所有的统计数据”是什么意思。它是否指下面受访者提供的基本屏幕截取？（“打印命令返回NULL”的意思是什么……用完整的代码请说明是哪个命令。） - IRTFM
提取模型信息的最有效方法是使用“broom”包。请查看此处：https://cran.r-project.org/web/packages/broom/index.html。我将使用这种方法发布答案... - AntoniosK
@42，计算出的所有内容，包括没有在屏幕上显示但必要的统计数据。至于打印，请执行x = print((fit = rms::ols(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species, data = iris))) ， > x NULL。求解该问题（即我）的回答者提供了黑客解决方案。 - CoderGuy123
2个回答
5
这里是一个hack的解决方案，我们捕获print命令的输出结果。
#parser
get_model_stats = function(x, precision=60) {

  # remember old number formatting function
  # (which would round and transforms p-values to formats like "<0.01")
  old_format_np = rms::formatNP
  # substitute it with a function which will print out as many digits as we want
  assignInNamespace("formatNP", function(x, ...) formatC(x, format="f", digits=precision), "rms")

  # remember old width setting
  old_width = options('width')$width
  # substitute it with a setting making sure the table will not wrap
  options(width=old_width + 4 * precision)

  # actually print the data and capture it
  cap = capture.output(print(x))

  # restore original settings
  options(width=old_width)
  assignInNamespace("formatNP", old_format_np, "rms")
  
  #model stats
  stats = c()
  stats$R2.adj = str_match(cap, "R2 adj\\s+ (\\d\\.\\d+)") %>% na.omit() %>% .[, 2] %>% as.numeric()
  
  #coef stats lines
  coef_lines = cap[which(str_detect(cap, "Coef\\s+S\\.E\\.")):(length(cap) - 1)]
  
  #parse
  coef_lines_table = suppressWarnings(readr::read_table(coef_lines %>% stringr::str_c(collapse = "\n")))
  colnames(coef_lines_table)[1] = "Predictor"
  
  list(
    stats = stats,
    coefs = coef_lines_table
  )
}

例子：

> get_model_stats(fit)
$stats
$stats$R2.adj
[1] 0.86


$coefs
# A tibble: 6 x 5
           Predictor  Coef  S.E.     t `Pr(>|t|)`
               <chr> <dbl> <dbl> <dbl>      <chr>
1          Intercept  2.17 0.280   7.8    <0.0001
2        Sepal.Width  0.50 0.086   5.8    <0.0001
3       Petal.Length  0.83 0.069  12.1    <0.0001
4        Petal.Width -0.32 0.151  -2.1     0.0389
5 Species=versicolor -0.72 0.240  -3.0     0.0031
6  Species=virginica -1.02 0.334  -3.1     0.0026

这仍存在问题，例如p值未返回为数字且仅有4位数字，在某些情况下可能会导致问题。更新的代码应该提取任意精度的数字。

在使用长变量名时要特别小心，因为这些变量名可能会将表格分成多行，并在输出中引入缺失值（NA），尽管统计数据已经在其中！

要将p值作为数字获取并且不进行四舍五入，请在调用get_model_stats之前添加assignInNamespace('formatNP', function(x, ...) x, 'rms')。只要prModFit中的代码使用formatNP来格式化p值，就可以正常工作 :) - krassowski
1

包 broom 是提取模型信息的好方法。

library(pacman)
library(rms)
library(broom)

fit = ols(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species, 
               data = iris)

tidy(summary.lm(fit))

#                 term   estimate  std.error statistic      p.value
# 1          Intercept  2.1712663 0.27979415  7.760227 1.429502e-12
# 2        Sepal.Width  0.4958889 0.08606992  5.761466 4.867516e-08
# 3       Petal.Length  0.8292439 0.06852765 12.100867 1.073592e-23
# 4        Petal.Width -0.3151552 0.15119575 -2.084418 3.888826e-02
# 5 Species=versicolor -0.7235620 0.24016894 -3.012721 3.059634e-03
# 6  Species=virginica -1.0234978 0.33372630 -3.066878 2.584344e-03

glance(fit)

#   r.squared adj.r.squared     sigma statistic      p.value df    logLik      AIC      BIC df.residual
# 1 0.8673123      0.862705 0.3068261   188.251 2.666942e-61  6 -32.55801 79.11602 100.1905         144

对象fit还包含一些易于访问的信息，您可以获取并存储在数据帧中：

fit$coefficients

# Intercept        Sepal.Width       Petal.Length        Petal.Width Species=versicolor  Species=virginica 
# 2.1712663          0.4958889          0.8292439         -0.3151552         -0.7235620         -1.0234978

fit$stats

#           n  Model L.R.        d.f.          R2           g       Sigma 
# 150.0000000 302.9635115   5.0000000   0.8673123   0.8820479   0.3068261
1
我注意到broom方法会跳过某些预测变量，因为它们被赋予了NA值。我有一个私人数据集，其中出现了这种情况。然而，使用coef获取系数时会包括这些NA值。 - CoderGuy123
@Deleet coef 返回一个数据框吗？如果不是，您可以创建一个包含所有预测变量名称的数据框，然后将其与 broom 输出连接起来。这样，您将获得缺失值的 NA，同时也会得到一个数据框。 - AntoniosK
2
broom 目前还不支持 rms 包中的对象，因此这种方法是一种 hack，我无法保证其正确性。rms 的支持正在（缓慢地）进行中，但目前我建议使用 capture.output() 方法。 - alexpghayes
2021年，broom是否支持rms？ - krassowski
1
他们在0.7.3中重新添加了对此的支持，但这不是正确的方式。https://github.com/harrelfe/rms/issues/93#issuecomment-762397451 - CoderGuy123
显示剩余5条评论
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接