将 p 值添加到 polr 模型用于模型摘要

问题描述

我知道 polr 不会给出 p 值,因为它们不是很可靠。尽管如此,我还是想将它们添加到我的 modelsummary (Vignette) 输出中。我知道要获得如下值:

library(MASS)
polr_res <- polr(as.ordered(rep77) ~ foreign + length + mpg,Hess=TRUE,data=fullauto);summary(polr_res)

Call:
polr_res(formula = as.ordered(rep77) ~ foreign + length + mpg,data = fullauto,Hess = TRUE)

## coefficient test
library("AER")
coeftest(polr_res)

模型摘要

因为 polr 没有 p 值,所以我不能在我的 modelsummary(models,stars=TRUE)调用 models(其中包括其他有 p 值并且我想显示星星的模型) .

library(modelsummary)
models <- list(
"Ordinal Probit" = polr_res,)
# model_names <- c("OLS","")
modelsummary(models,stars=TRUE)

我首先尝试简单地将 p 值添加tidy 对象,但我无法将该对象添加models 的列表中。

polr_pval <- coeftest(polr)[,4]
polr_pval <- as.data.frame(polr_pval)
tidy_polr <- tidy(polr)
tidy_polr[,5] <- polr_pval

小插图描述我可以制作一个适应 polr自定义类,但我不明白如何: https://vincentarelbundock.github.io/modelsummary/articles/modelsummary.html#customizing-existing-models-part-i-

https://vincentarelbundock.github.io/modelsummary/articles/modelsummary.html#customizing-existing-models-part-ii-

谁能帮我解决这个问题?

编辑:

我发布了一个编辑,显示了我在使用 vincent 的答案时遇到的问题,带有 R version 3.6.1 (2019-07-05)。如果您遇到此问题,(最好)更新到 R 版本 4.0.0 或从 Github 下载 modelsummary 的更新(另请参阅下面 vincent评论)。:

library(remotes)
remotes::install_github('vincentarelbundock/modelsummary')

输出

enter image description here

R 的数据

fullauto <- structure(list(make = structure(c(1,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23),label = "Make",format.stata = "%8.0g",class = c("haven_labelled","vctrs_vctr","double"),labels = c(AMC = 1,Audi = 2,BMW = 3,Buick = 4,Cad. = 5,Chev. = 6,Datsun = 7,Dodge = 8,Fiat = 9,Ford = 10,Honda = 11,Linc. = 12,Mazda = 13,Merc. = 14,Olds = 15,Peugeot = 16,Plym. = 17,Pont. = 18,Renault = 19,Subaru = 20,Toyota = 21,VW = 22,Volvo = 23)),model = structure(c(1,5000,320,200,210,510,810,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,98,604,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,260),label = "Model",labels = c(Concord = 1,Pacer = 2,Spirit = 3,Fox = 4,Century = 5,Electra = 6,LeSabre = 7,Opel = 8,Regal = 9,Riviera = 10,Skylark = 11,Deville = 12,Eldrado = 13,Seville = 14,Chevette = 15,Impala = 16,Malibu = 17,MCarlo = 18,Monza = 19,Nova = 20,Colt = 21,Diplomat = 22,Magnum = 23,StRegis = 24,STrada = 25,Fiesta = 26,Mustang = 27,Accord = 28,Civic = 29,Cntntl = 30,`Mark V` = 31,Vrsills = 32,GLC = 33,Bobcat = 34,Cougar = 35,`XR-7` = 36,Marquis = 37,Monarch = 38,Zephyr = 39,Cutlass = 40,CutlSupr = 41,`Delta 88` = 42,Omega = 43,Starfire = 44,Toronado = 45,Arrow = 46,Champ = 47,Horizon = 48,Sapporo = 49,Volare = 50,Catalina = 51,Firebird = 52,GranPrix = 53,`Le Mans` = 54,Phoenix = 55,Sunbird = 56,`Le Car` = 57,Subaru = 58,Celica = 59,Corolla = 60,Corona = 61,Rabbit = 62,Diesel = 63,Scirocco = 64,Dasher = 65)),price = structure(c(4099,4749,3799,6295,9690,9735,4816,7827,5788,4453,5189,10372,4082,11385,14500,15906,3299,5705,4504,5104,3667,3955,6229,4589,5079,8129,3984,4010,5886,6342,4296,4389,4187,5799,4499,11497,13594,13466,3995,3829,5379,6303,6165,4516,3291,4733,5172,4890,4181,4195,10371,8814,12990,4647,4425,4482,6486,4060,5798,4934,5222,4723,4424,4172,3895,3798,5899,3748,5719,4697,5397,6850,7140,11995),label = "Price",format.stata = "%8.0g"),mpg = structure(c(22,17),label = "Mileage (mpg)",rep78 = structure(c(3,NA,5),label = "Repair Record 1978",format.stata = "%9.0g",labels = c(Poor = 1,Fair = 2,Average = 3,Good = 4,Excellent = 5)),rep77 = structure(c(2,3),label = "Repair Record 1977",hdroom = structure(c(2.5,2.5,4.5,3.5,1.5,2.5),label = "Headroom (in.)",format.stata = "%6.1f"),rseat = structure(c(27.5,25.5,18.5,31.5,30.5,28.5,29.5,23.5,26.5,21.5,24.5,37.5,29.5),label = "Rear Seat (in.)",trunk = structure(c(11,14),label = "Trunk space (cu. ft.)",weight = structure(c(2930,3350,2640,2070,2830,2650,3250,4080,3670,2230,3280,3880,3400,4330,3900,4290,2110,3690,3180,3220,2750,3430,2370,2020,2280,2120,3600,3740,2130,1800,2240,1760,4840,4720,3830,1980,2580,4130,3720,3370,3300,3310,2730,4030,3420,3260,2200,2520,3330,3700,3470,3210,3200,2690,1830,2050,2410,2670,1930,2040,1990,2160,3170),label = "Weight (lbs.)",length = structure(c(186,173,168,174,189,177,196,222,218,170,207,221,204,163,212,193,179,197,165,184,206,220,161,147,172,149,233,230,201,154,169,217,198,195,180,192,157,182,214,199,203,142,164,175,155,156,193),label = "Length (in.)",turn = structure(c(40,37),label = "Turn Circle (ft.) ",displ = structure(c(121,258,121,97,131,350,231,304,425,250,151,119,85,146,318,225,105,140,107,91,400,302,86,79,134,89,90,163),label = "displacement (cu. in.)",gratio = structure(c(3.57999992370605,2.52999997138977,3.07999992370605,3.70000004768372,3.20000004768372,3.64000010490417,2.9300000667572,2.41000008583069,2.73000001907349,2.86999988555908,2.27999997138977,2.19000005722046,2.24000000953674,2.55999994277954,3.89000010490417,3.53999996185303,3.54999995231628,2.47000002861023,2.94000005722046,3.36999988555908,3.15000009536743,3.04999995231628,3.29999995231628,3.73000001907349,2.75,2.25999999046326,2.4300000667572,3.57999992370605,2.97000002861023,3.23000001907349,3.72000002861023,3.80999994277954,3.05999994277954,3.21000003814697,3.77999997138977,3.74000000953674,2.98000001907349),label = "Gear Ratio",format.stata = "%6.2f"),order = structure(c(1,66,67,68,69,70,71,72,73,74),label = "Original order",foreign = structure(c(0,1),label = "Foreign",labels = c(Domestic = 0,Foreign = 1
    )),wgtd = structure(c(2930,NA),format.stata = "%9.0g"),wgtf = structure(c(NA,format.stata = "%9.0g")),label = "Automobile Models",row.names = c(NA,-74L),class = c("tbl_df","tbl","data.frame"))

解决方法

我认为实现这一点的最简单方法是定义一个 tidy_custom.polr method as described here in the documentation.。例如,你可以这样做:

library(MASS)
library(AER)
library(modelsummary)

tidy_custom.polr <- function(x,...) {
  s <- coeftest(x)
  out <- data.frame(
    term = row.names(s),p.value = s[,"Pr(>|z|)"])
  out
}

mod = list(
  "LM" = lm(gear ~ hp + mpg,data = mtcars),"POLR" = polr(as.ordered(gear) ~ hp + mpg,data = mtcars))

modelsummary(mod,stars = TRUE)

enter image description here