r - 从 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 fit 函数打印了很多有用的东西，包括标准误差和调整后的 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包创建者解释了获取这些的解决方案:

On Wed, 11 Aug 2010, david dav wrote:

Hi, I would like to extract the coefficients of a logistic regression (estimates and standard error as well) in lrm as in glm with

summary(fit.glm)$coef

Thanks David

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 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() .我们可以 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

这是一个黑客解决方案，我们捕获 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)，即使统计数据在那里!