forked from tidymodels/broom
-
Notifications
You must be signed in to change notification settings - Fork 0
/
mass-polr-tidiers.R
133 lines (124 loc) · 3.84 KB
/
mass-polr-tidiers.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#' @templateVar class polr
#' @template title_desc_tidy
#'
#' @param x A `polr` object returned from [MASS::polr()].
#' @template param_confint
#' @template param_exponentiate
#' @param p.values Logical. Should p-values be returned,
#' based on chi-squared tests from [MASS::dropterm()]. Defaults to FAlSE
#' @template param_unused_dots
#'
#' @examples
#'
#' library(MASS)
#'
#' fit <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
#'
#' tidy(fit, exponentiate = TRUE, conf.int = TRUE)
#'
#' glance(fit)
#' augment(fit, type.predict = "class")
#'
#' fit2 <- polr(factor(gear) ~ am + mpg + qsec, data = mtcars)
#' tidy(fit, p.values = TRUE)
#'
#' @evalRd return_tidy(regression = TRUE)
#'
#' @details In `broom 0.7.0` the `coefficient_type` column was renamed to
#' `coef.type`, and the contents were changed as well. Now the contents
#' are `coefficient` and `scale`, rather than `coefficient` and `zeta`.
#'
#' Calculating p.values with the `dropterm()` function is the approach
#' suggested by the MASS package author
#' (https://r.789695.n4.nabble.com/p-values-of-plor-td4668100.html). This
#' approach is computationally intensive, so that p.values are only
#' returned if requested explicitly. Additionally, it only works for
#' models containing no variables with more than two categories. If this
#' condition is not met, a message is shown and NA is returned instead of
#' p-values.
#'
#' @aliases polr_tidiers
#' @export
#' @seealso [tidy], [MASS::polr()]
#' @family ordinal tidiers
tidy.polr <- function(x, conf.int = FALSE, conf.level = 0.95,
exponentiate = FALSE, p.values = FALSE, ...) {
ret <- as_tibble(coef(summary(x)), rownames = "term")
colnames(ret) <- c("term", "estimate", "std.error", "statistic")
if (conf.int) {
ci <- broom_confint_terms(x, level = conf.level)
ret <- dplyr::left_join(ret, ci, by = "term")
}
if (exponentiate) {
ret <- exponentiate(ret)
if (p.values) {
sig <- MASS::dropterm(x, test = "Chisq")
p <- sig %>% dplyr::select(`Pr(Chi)`) %>% dplyr::pull() %>% .[-1]
terms <- purrr::map(rownames(sig)[-1], function(x)
ret$term[stringr::str_detect(ret$term, stringr::fixed(x))]) %>% unlist()
if (length(p) == length(terms)) {
ret <- dplyr::left_join(ret, tibble::tibble(term = terms, p.value = p), by = "term")
} else {
message("p-values can presently only be returned for models that contain
no categorical variables with more than two levels")
ret$p.value <- NA
}
}
}
mutate(
ret,
coef.type = if_else(term %in% names(x$zeta), "scale", "coefficient")
)
}
#' @templateVar class polr
#' @template title_desc_glance
#'
#' @inherit tidy.polr params examples
#'
#' @evalRd return_glance(
#' "edf",
#' "logLik",
#' "AIC",
#' "BIC",
#' "deviance",
#' "df.residual",
#' "nobs"
#' )
#'
#' @export
#' @seealso [tidy], [MASS::polr()]
#' @family ordinal tidiers
glance.polr <- function(x, ...) {
tibble(
edf = x$edf,
logLik = as.numeric(stats::logLik(x)),
AIC = stats::AIC(x),
BIC = stats::BIC(x),
deviance = stats::deviance(x),
df.residual = stats::df.residual(x),
nobs = stats::nobs(x)
)
}
#' @templateVar class polr
#' @template title_desc_augment
#'
#' @inherit tidy.polr params examples
#' @template param_data
#' @template param_newdata
#'
#' @param type.predict Which type of prediction to compute,
#' passed to `MASS:::predict.polr()`. Only supports `"class"` at
#' the moment.
#'
#' @export
#' @seealso [tidy()], [MASS::polr()]
#' @family ordinal tidiers
#'
augment.polr <- function(x, data = model.frame(x), newdata = NULL,
type.predict = c("class"), ...) {
type <- rlang::arg_match(type.predict)
df <- if (is.null(newdata)) data else newdata
df <- as_broom_tibble(df)
df$.fitted <- predict(object = x, newdata = df, type = type)
df
}