library(MASS) # Has to come first because of dplyr::select
library(tidyverse)
library(broom)
library(haven)
example_df <- read_stata("http://stats.idre.ucla.edu/stat/data/ologit.dta") %>%
mutate(apply = factor(apply, levels = 1:3,
labels = c("unlikely", "somewhat likely", "very likely"),
ordered = TRUE))
model <- polr(apply ~ pared + public + gpa,
data = example_df, Hess = TRUE)
model_results <- tidy(model, exponentiate = TRUE, conf.int = TRUE) %>%
mutate(p_value = pnorm(abs(statistic), lower.tail = FALSE) * 2)
model_results
#> # A tibble: 5 x 8
#> term estimate std.error statistic conf.low conf.high coefficient_type p_value
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
#> 1 gpa 1.62e0 0.504 0.955 0.604 4.39 coefficient 0.340
#> 2 pared 1.65e0 0.406 1.23 0.729 3.62 coefficient 0.219
#> 3 publ… 2.22e0 0.479 1.67 0.848 5.63 coefficient 0.0955
#> 4 some… 5.45e8 721. 0.0279 NA NA zeta 0.978
#> 5 unli… 2.03e1 1.54 1.95 NA NA zeta 0.0507
newdata <- expand_grid(pared = c(0, 1),
public = c(0, 1),
gpa = seq(1.9, 4, by = 0.1))
predicted_probs <- cbind(newdata,
predict(model, newdata = newdata, type = "probs")) %>%
pivot_longer(cols = c("unlikely", "somewhat likely", "very likely"),
names_to = "apply", values_to = "probs") %>%
mutate(apply = factor(apply, levels = levels(example_df$apply), ordered = TRUE))
ggplot(predicted_probs, aes(x = gpa, y = probs, color = apply)) +
geom_line() +
facet_grid(pared ~ public)
Created on 2020-06-02 by the reprex package (v0.3.0)