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

Question

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"

包创建者在a 7 year old question on R help 处解释了获取这些的解决方案：

2010 年 8 月 11 日，星期三，大卫 dav 写道：

嗨， 我想提取逻辑回归的系数 （估计和标准误差）在 lrm 和 glm 中一样

总结(fit.glm)$coef

谢谢 大卫

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

但这些不会很有帮助，除非在微不足道的情况下 一切都是线性的，没有相互作用，因素有两个层次。

弗兰克

根据作者的说法，解决方案不是最优的。这让人想知道如何计算显示的值。跟踪代码会导致搜索未记录的包代码（包代码is on Github）。 IE。我们从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 调整。在打印函数中计算：

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()。我们可以look it up 并找到 p 值计算等。不幸的是，print 命令返回NULL，因此这些值无法在任何地方进行编程重用：

> x = print((fit = rms::ols(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species, data = iris)))
#printed output...
> x
NULL

如何获得所有统计数据？

Answer 1

这是一个 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)，即使统计信息在其中！

Answer 2

包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