如何从 R 中的 rms 包的 ols() 模型中获取 p.values 系数
How to get p.values of coefficients from ols() model of rms package in R
我运行回归使用mtcars数据集和rms包的ols()函数如下:
> library(rms)
> mod_ols = ols(mpg~wt+qsec+am+hp+vs+drat+cyl, mtcars)
>
> mod_ols
Linear Regression Model
ols(formula = mpg ~ wt + qsec + am + hp + vs + drat + cyl, data = mtcars)
Model Likelihood Discrimination
Ratio Test Indexes
Obs 32 LR chi2 62.73 R2 0.859
sigma 2.5703 d.f. 7 R2 adj 0.818
d.f. 24 Pr(> chi2) 0.0000 g 6.519
Residuals
Min 1Q Median 3Q Max
-3.3615 -1.7407 -0.2537 1.0893 4.5929
Coef S.E. t Pr(>|t|)
Intercept 15.7714 16.5491 0.95 0.3501
wt -3.0826 1.0888 -2.83 0.0092
qsec 0.7426 0.6632 1.12 0.2739
am 2.5945 1.8802 1.38 0.1803
hp -0.0182 0.0162 -1.12 0.2728
vs 0.2321 2.0291 0.11 0.9099
drat 0.6387 1.5186 0.42 0.6778
cyl 0.0309 0.9091 0.03 0.9732
生成的模型结构如下:
> str(mod_ols)
List of 18
$ coefficients : Named num [1:8] 15.7714 -3.0826 0.7426 2.5945 -0.0182 ...
..- attr(*, "names")= chr [1:8] "Intercept" "wt" "qsec" "am" ...
$ residuals : Named num [1:32] -2.193 -1.823 -3.362 0.714 1.81 ...
..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
$ effects : Named num [1:32] -113.65 -29.12 -9.1 5.12 3.04 ...
..- attr(*, "names")= chr [1:32] "Intercept" "wt" "qsec" "am" ...
$ rank : int 8
$ fitted.values : Named num [1:32] 23.2 22.8 26.2 20.7 16.9 ...
..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
$ assign :List of 7
..$ wt : int 2
..$ qsec: int 3
..$ am : int 4
..$ hp : int 5
..$ vs : int 6
..$ drat: int 7
..$ cyl : int 8
$ qr :List of 5
..$ qr : num [1:32, 1:8] -5.657 0.177 0.177 0.177 0.177 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
.. .. ..$ : chr [1:8] "Intercept" "wt" "qsec" "am" ...
.. ..- attr(*, "assign")= int [1:8] 0 1 2 3 4 5 6 7
..$ qraux: num [1:8] 1.18 1.05 1.08 1.18 1.03 ...
..$ pivot: int [1:8] 1 2 3 4 5 6 7 8
..$ tol : num 0.0000001
..$ rank : int 8
..- attr(*, "class")= chr "qr"
$ df.residual : int 24
$ var : num [1:8, 1:8] 273.872 3.919 -9.225 -13.12 -0.054 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:8] "Intercept" "wt" "qsec" "am" ...
.. ..$ : chr [1:8] "Intercept" "wt" "qsec" "am" ...
$ stats : Named num [1:6] 32 62.733 7 0.859 6.519 ...
..- attr(*, "names")= chr [1:6] "n" "Model L.R." "d.f." "R2" ...
$ linear.predictors: Named num [1:32] 23.2 22.8 26.2 20.7 16.9 ...
..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
$ call : language ols(formula = mpg ~ wt + qsec + am + hp + vs + drat + cyl, data = mtcars)
$ terms :Classes 'terms', 'formula' length 3 mpg ~ wt + qsec + am + hp + vs + drat + cyl
.. ..- attr(*, "variables")= language list(mpg, wt, qsec, am, hp, vs, drat, cyl)
.. ..- attr(*, "factors")= int [1:8, 1:7] 0 1 0 0 0 0 0 0 0 0 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:8] "mpg" "wt" "qsec" "am" ...
.. .. .. ..$ : chr [1:7] "wt" "qsec" "am" "hp" ...
.. ..- attr(*, "term.labels")= chr [1:7] "wt" "qsec" "am" "hp" ...
.. ..- attr(*, "order")= int [1:7] 1 1 1 1 1 1 1
.. ..- attr(*, "intercept")= num 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(mpg, wt, qsec, am, hp, vs, drat, cyl)
.. ..- attr(*, "dataClasses")= Named chr [1:8] "numeric" "numeric" "numeric" "numeric" ...
.. .. ..- attr(*, "names")= chr [1:8] "mpg" "wt" "qsec" "am" ...
.. ..- attr(*, "formula")=Class 'formula' length 3 mpg ~ wt + qsec + am + hp + vs + drat + cyl
.. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
$ Design :List of 11
..$ name : chr [1:7] "wt" "qsec" "am" "hp" ...
..$ label : chr [1:7] "wt" "qsec" "am" "hp" ...
..$ units : Named chr [1:7] "" "" "" "" ...
.. ..- attr(*, "names")= chr [1:7] "wt" "qsec" "am" "hp" ...
..$ colnames : chr [1:7] "wt" "qsec" "am" "hp" ...
..$ assume : chr [1:7] "asis" "asis" "asis" "asis" ...
..$ assume.code : int [1:7] 1 1 1 1 1 1 1
..$ parms : list()
..$ limits : list()
..$ values : list()
..$ nonlinear :List of 7
.. ..$ wt : logi FALSE
.. ..$ qsec: logi FALSE
.. ..$ am : logi FALSE
.. ..$ hp : logi FALSE
.. ..$ vs : logi FALSE
.. ..$ drat: logi FALSE
.. ..$ cyl : logi FALSE
..$ interactions: list()
$ non.slopes : num 1
$ na.action : NULL
$ scale.pred : chr "mpg"
$ fail : logi FALSE
- attr(*, "class")= chr [1:3] "ols" "rms" "lm"
>
但是,我无法在 mod_ols 结构中找到 p.values。如何从 rms 包的 ols() 函数生成的模型中获取系数的 P 值?感谢您的帮助。
使用summary.lm
:
summary.lm(mod_ols)$coefficients[ , "Pr(>|t|)"]
# Intercept wt qsec am hp vs
#0.350087341 0.009234825 0.273857810 0.180345204 0.272787068 0.909889126
# drat cyl
#0.677806300 0.973163043
以下每一行都可以解决问题:
summary.lm(mod_ols)$coefficients[ , 4]
coef(summary.lm(mod_ols))[ , 4]
coef(summary.lm(mod_ols))[ , "Pr(>|t|)"]
我运行回归使用mtcars数据集和rms包的ols()函数如下:
> library(rms)
> mod_ols = ols(mpg~wt+qsec+am+hp+vs+drat+cyl, mtcars)
>
> mod_ols
Linear Regression Model
ols(formula = mpg ~ wt + qsec + am + hp + vs + drat + cyl, data = mtcars)
Model Likelihood Discrimination
Ratio Test Indexes
Obs 32 LR chi2 62.73 R2 0.859
sigma 2.5703 d.f. 7 R2 adj 0.818
d.f. 24 Pr(> chi2) 0.0000 g 6.519
Residuals
Min 1Q Median 3Q Max
-3.3615 -1.7407 -0.2537 1.0893 4.5929
Coef S.E. t Pr(>|t|)
Intercept 15.7714 16.5491 0.95 0.3501
wt -3.0826 1.0888 -2.83 0.0092
qsec 0.7426 0.6632 1.12 0.2739
am 2.5945 1.8802 1.38 0.1803
hp -0.0182 0.0162 -1.12 0.2728
vs 0.2321 2.0291 0.11 0.9099
drat 0.6387 1.5186 0.42 0.6778
cyl 0.0309 0.9091 0.03 0.9732
生成的模型结构如下:
> str(mod_ols)
List of 18
$ coefficients : Named num [1:8] 15.7714 -3.0826 0.7426 2.5945 -0.0182 ...
..- attr(*, "names")= chr [1:8] "Intercept" "wt" "qsec" "am" ...
$ residuals : Named num [1:32] -2.193 -1.823 -3.362 0.714 1.81 ...
..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
$ effects : Named num [1:32] -113.65 -29.12 -9.1 5.12 3.04 ...
..- attr(*, "names")= chr [1:32] "Intercept" "wt" "qsec" "am" ...
$ rank : int 8
$ fitted.values : Named num [1:32] 23.2 22.8 26.2 20.7 16.9 ...
..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
$ assign :List of 7
..$ wt : int 2
..$ qsec: int 3
..$ am : int 4
..$ hp : int 5
..$ vs : int 6
..$ drat: int 7
..$ cyl : int 8
$ qr :List of 5
..$ qr : num [1:32, 1:8] -5.657 0.177 0.177 0.177 0.177 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
.. .. ..$ : chr [1:8] "Intercept" "wt" "qsec" "am" ...
.. ..- attr(*, "assign")= int [1:8] 0 1 2 3 4 5 6 7
..$ qraux: num [1:8] 1.18 1.05 1.08 1.18 1.03 ...
..$ pivot: int [1:8] 1 2 3 4 5 6 7 8
..$ tol : num 0.0000001
..$ rank : int 8
..- attr(*, "class")= chr "qr"
$ df.residual : int 24
$ var : num [1:8, 1:8] 273.872 3.919 -9.225 -13.12 -0.054 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:8] "Intercept" "wt" "qsec" "am" ...
.. ..$ : chr [1:8] "Intercept" "wt" "qsec" "am" ...
$ stats : Named num [1:6] 32 62.733 7 0.859 6.519 ...
..- attr(*, "names")= chr [1:6] "n" "Model L.R." "d.f." "R2" ...
$ linear.predictors: Named num [1:32] 23.2 22.8 26.2 20.7 16.9 ...
..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
$ call : language ols(formula = mpg ~ wt + qsec + am + hp + vs + drat + cyl, data = mtcars)
$ terms :Classes 'terms', 'formula' length 3 mpg ~ wt + qsec + am + hp + vs + drat + cyl
.. ..- attr(*, "variables")= language list(mpg, wt, qsec, am, hp, vs, drat, cyl)
.. ..- attr(*, "factors")= int [1:8, 1:7] 0 1 0 0 0 0 0 0 0 0 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:8] "mpg" "wt" "qsec" "am" ...
.. .. .. ..$ : chr [1:7] "wt" "qsec" "am" "hp" ...
.. ..- attr(*, "term.labels")= chr [1:7] "wt" "qsec" "am" "hp" ...
.. ..- attr(*, "order")= int [1:7] 1 1 1 1 1 1 1
.. ..- attr(*, "intercept")= num 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(mpg, wt, qsec, am, hp, vs, drat, cyl)
.. ..- attr(*, "dataClasses")= Named chr [1:8] "numeric" "numeric" "numeric" "numeric" ...
.. .. ..- attr(*, "names")= chr [1:8] "mpg" "wt" "qsec" "am" ...
.. ..- attr(*, "formula")=Class 'formula' length 3 mpg ~ wt + qsec + am + hp + vs + drat + cyl
.. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
$ Design :List of 11
..$ name : chr [1:7] "wt" "qsec" "am" "hp" ...
..$ label : chr [1:7] "wt" "qsec" "am" "hp" ...
..$ units : Named chr [1:7] "" "" "" "" ...
.. ..- attr(*, "names")= chr [1:7] "wt" "qsec" "am" "hp" ...
..$ colnames : chr [1:7] "wt" "qsec" "am" "hp" ...
..$ assume : chr [1:7] "asis" "asis" "asis" "asis" ...
..$ assume.code : int [1:7] 1 1 1 1 1 1 1
..$ parms : list()
..$ limits : list()
..$ values : list()
..$ nonlinear :List of 7
.. ..$ wt : logi FALSE
.. ..$ qsec: logi FALSE
.. ..$ am : logi FALSE
.. ..$ hp : logi FALSE
.. ..$ vs : logi FALSE
.. ..$ drat: logi FALSE
.. ..$ cyl : logi FALSE
..$ interactions: list()
$ non.slopes : num 1
$ na.action : NULL
$ scale.pred : chr "mpg"
$ fail : logi FALSE
- attr(*, "class")= chr [1:3] "ols" "rms" "lm"
>
但是,我无法在 mod_ols 结构中找到 p.values。如何从 rms 包的 ols() 函数生成的模型中获取系数的 P 值?感谢您的帮助。
使用summary.lm
:
summary.lm(mod_ols)$coefficients[ , "Pr(>|t|)"]
# Intercept wt qsec am hp vs
#0.350087341 0.009234825 0.273857810 0.180345204 0.272787068 0.909889126
# drat cyl
#0.677806300 0.973163043
以下每一行都可以解决问题:
summary.lm(mod_ols)$coefficients[ , 4]
coef(summary.lm(mod_ols))[ , 4]
coef(summary.lm(mod_ols))[ , "Pr(>|t|)"]