Значения P из функции регрессии fastbw пакета rms

Я пытаюсь использовать функцию fastbw пакета rms для обратной регрессии следующим образом (используя набор данных mtcars):

> mod = ols(mpg~am+vs+cyl+drat+wt+gear, mtcars) 
> mod

Linear Regression Model

ols(formula = mpg ~ am + vs + cyl + drat + wt + gear, data = mtcars)

                Model Likelihood     Discrimination    
                   Ratio Test           Indexes        
Obs       32    LR chi2     58.26    R2       0.838    
sigma 2.7008    d.f.            6    R2 adj   0.799    
d.f.      25    Pr(> chi2) 0.0000    g        6.383    

Residuals

    Min      1Q  Median      3Q     Max 
-4.3807 -1.4314 -0.5405  1.5828  5.4703 

          Coef    S.E.   t     Pr(>|t|)
Intercept 39.9804 8.8745  4.51 0.0001  
am         1.5981 1.9927  0.80 0.4301  
vs         0.8011 1.9201  0.42 0.6801  
cyl       -1.3163 0.7033 -1.87 0.0730  
drat       0.3488 1.6201  0.22 0.8313  
wt        -3.0390 0.9510 -3.20 0.0038  
gear      -1.1450 1.1420 -1.00 0.3256  


> modbw = fastbw(mod)
> modbw

 Deleted Chi-Sq d.f. P      Residual d.f. P      AIC   R2   
 drat    0.05   1    0.8296 0.05     1    0.8296 -1.95 0.838
 vs      0.17   1    0.6800 0.22     2    0.8974 -3.78 0.837
 am      0.58   1    0.4473 0.79     3    0.8509 -5.21 0.833
 gear    0.42   1    0.5194 1.21     4    0.8766 -6.79 0.830

Approximate Estimates after Deleting Factors

            Coef   S.E. Wald Z          P
Intercept 39.686 1.8040 21.999 0.00000000
cyl       -1.508 0.4362 -3.457 0.00054706
wt        -3.191 0.7962 -4.008 0.00006128

Factors in Final Model

[1] cyl wt

Ниже приведена структура этой модели:

> str(modbw)
List of 10
 $ result         : num [1:4, 1:8] 0.0463 0.1701 0.5775 0.4152 1 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:4] "drat" "vs" "am" "gear"
  .. ..$ : chr [1:8] "Chi-Sq" "d.f." "P" "Residual" ...
 $ names.kept     : chr [1:2] "cyl" "wt"
 $ factors.kept   : int [1:2] 3 5
 $ factors.deleted: int [1:4] 4 2 1 6
 $ parms.kept     : int [1:3] 1 4 6
 $ parms.deleted  : int [1:4] 5 3 2 7
 $ coefficients   : Named num [1:3] 39.69 -1.51 -3.19
  ..- attr(*, "names")= chr [1:3] "Intercept" "cyl" "wt"
 $ var            : num [1:3, 1:3] 3.254 -0.303 -0.358 -0.303 0.19 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "Intercept" "cyl" "wt"
  .. ..$ : chr [1:3] "Intercept" "cyl" "wt"
 $ Coefficients   : num [1:4, 1:7] 41.26 43.17 42.39 39.69 1.68 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:7] "Intercept" "am" "vs" "cyl" ...
 $ force          : NULL
 - attr(*, "class")= chr "fastbw"

Ниже приведены структуры выходных данных сводки (функция summary.lm не работает на этой модели):

> summary(modbw)
                Length Class  Mode     
result          32     -none- numeric  
names.kept       2     -none- character
factors.kept     2     -none- numeric  
factors.deleted  4     -none- numeric  
parms.kept       3     -none- numeric  
parms.deleted    4     -none- numeric  
coefficients     3     -none- numeric  
var              9     -none- numeric  
Coefficients    28     -none- numeric  
force            0     -none- NULL     
> 
> summary.lm(modbw)
Error in if (p == 0) { : argument is of length zero

Но я не могу найти значения P ни в одном из них. Как я могу получить список значений P для окончательной модели функции fastbw?


r regression linear-regression rms
person rnso    schedule 06.06.2015    source источник
comment
Я почти уверен, что fastbw не предназначен для статистического вывода... у вас есть доступ к книге Харрелла Стратегии регрессионного моделирования?   -  person Ben Bolker    schedule 06.06.2015
comment
Значение P появляется на выходе. Так что это должно быть разумно важно и здесь. Почему его не видно нигде в структурах?   -  person rnso    schedule 06.06.2015

Ответы (2)


arrow_upward
2
arrow_downward

Они вычисляются на лету. Покопавшись внутри rms:::print.fastbw (метод print для объектов класса fastbw) можно найти:

   cof <- coef(x)
   vv <- if (length(cof) > 1) diag(x$var) else x$var
   z <- cof/sqrt(vv)
   stats <- cbind(cof, sqrt(vv), z, 1 - pchisq(z^2, 1))

(если вам нужны более точные маленькие p-значения, было бы лучше заменить pchisq(z^2,1,lower.tail=FALSE) на 1-pchisq(z^2,1))

Если x является вашим объектом fastbw, то последний столбец stats содержит ваши p-значения.

person Ben Bolker    schedule 06.06.2015

arrow_upward
2
arrow_downward

Вычисление p-значений происходит в функции print.fastbw и по какой-то причине они не возвращаются из функции. Я собирался использовать исходный код print.fastbw, чтобы пересчитать их самостоятельно, но обнаружил, что переписывание моей собственной функции print.fastbw для возврата p-значений происходит намного быстрее.

Вот переработанная функция (обратите внимание, что print2 не является универсальной):

print2.fastbw <- function (x, digits = 4, estimates = TRUE, ...) 
{
  res <- x$result
  fd <- x$factors.deleted
  if (length(fd)) {
    cres <- cbind(dimnames(res)[[1]], format(round(res[, 
                                                       1], 2)), format(res[, 2]), format(round(res[, 3], 
                                                                                               4)), format(round(res[, 4], 2)), format(res[, 5]), 
                  format(round(res[, 6], 4)), format(round(res[, 7], 
                                                           2)), if (ncol(res) > 7) 
                                                             format(round(res[, 8], 3)))
    dimnames(cres) <- list(rep("", nrow(cres)), c("Deleted", 
                                                  dimnames(res)[[2]]))
    cat("\n")
    if (length(x$force)) 
      cat("Parameters forced into all models:\n", paste(x$force, 
                                                        collapse = ", "), "\n\n")
    print(cres, quote = FALSE)
    if (estimates && length(x$coef)) {
      cat("\nApproximate Estimates after Deleting Factors\n\n")
      cof <- coef(x)
      vv <- if (length(cof) > 1) 
        diag(x$var)
      else x$var
      z <- cof/sqrt(vv)
      stats <- cbind(cof, sqrt(vv), z, 1 - pchisq(z^2, 
                                                  1))
      dimnames(stats) <- list(names(cof), c("Coef", "S.E.", 
                                            "Wald Z", "P"))
      return(stats)
    }
  }
  else cat("\nNo Factors Deleted\n")
  cat("\nFactors in Final Model\n\n")
  nk <- x$names.kept
  if (length(nk)) 
    print(nk, quote = FALSE)
  else cat("None\n")
}

Выход:

> results <- print2.fastbw(modbw)

 Deleted Chi-Sq d.f. P      Residual d.f. P      AIC   R2   
 drat    0.05   1    0.8296 0.05     1    0.8296 -1.95 0.838
 vs      0.17   1    0.6800 0.22     2    0.8974 -3.78 0.837
 am      0.58   1    0.4473 0.79     3    0.8509 -5.21 0.833
 gear    0.42   1    0.5194 1.21     4    0.8766 -6.79 0.830

Approximate Estimates after Deleting Factors

> results
               Coef      S.E.    Wald Z            P
Intercept 39.686261 1.8039853 21.999216 0.000000e+00
cyl       -1.507795 0.4362091 -3.456588 5.470608e-04
wt        -3.190972 0.7961871 -4.007817 6.128261e-05

И p-значения:

> results[,4]
   Intercept          cyl           wt 
0.000000e+00 5.470608e-04 6.128261e-05
person LyzandeR    schedule 06.06.2015