## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", 
  message = FALSE,
  warning = FALSE
)
library(ggplot2) #
#### 0. Initialization ####
# remotes::install_github('davidaarmstrong/VizTest')
library(VizTest)

# Datasets from
library(wooldridge) #
library(carData) #

# For wrangling and analysis
library(dplyr) #
library(tidyr) #
library(ggplot2)
library(patchwork) #
library(lme4) #
library(marginaleffects) #
library(multcomp) #

## -----------------------------------------------------------------------------
data(chickwts)
chickwts$feed <- reorder(chickwts$feed, 
                         chickwts$weight, 
                         FUN=mean)
chick_preds <- chickwts %>% 
  group_by(feed) %>% 
  summarise(ave_weight = mean(weight),
            n = n(),
            sd = sd(weight)) %>%
  mutate(se = sd/sqrt(n))
chick_est <- chick_preds$ave_weight
names(chick_est) <- chick_preds$feed

## -----------------------------------------------------------------------------
chick_vt_data <- make_vt_data(estimates = chick_est, variances = chick_preds$se^2, type="est_var")

## -----------------------------------------------------------------------------
chick_vt <- viztest(chick_vt_data, include_zero = FALSE)
chick_vt

## ----chick-plot1, fig.height=4, fig.width=6, out.width="75%", fig.cap="Chick Weights with Inferential COnfidence Intervals (wide gray bars) and standard 95% confidence interval (thin black line).  If the gray bars overlap, it indicates the pair of estimates whose intervals overlap are not statistically different from each other.  If the gray bars do not overlap, the pair of estimates are statistically different from each other."----
plot(chick_vt)

## -----------------------------------------------------------------------------
data(gpa1, package="wooldridge")
model1 <- lm(colGPA~skipped+alcohol+PC+male+car+job20,data=gpa1)
summary(model1)


## -----------------------------------------------------------------------------
gpa_vt <- viztest(model1, test_level = 0.05, 
                  range_levels = c(0.25,0.99), 
                  level_increment = 0.01, 
                  include_zero=FALSE)

gpa_vt


## ----gpa1, fig.height=6.5, fig.width=6, out.width="75%", fig.cap="GPA Model Coefficients with Inferential Confidence Intervals (80%, wide gray bars) and standard 95% confidence interval (thin black line).  If the gray bars overlap, it indicates the pair of estimates whose intervals overlap are not statistically different from each other.  If the gray bars do not overlap, the pair of estimates are statistically different from each other."----
plot(gpa_vt, make_plot = FALSE) %>%
  mutate(vbl = factor(vbl), 
         vbl = reorder(vbl, est)) %>%
  ggplot(aes(x = est, y = vbl)) +
  geom_vline(xintercept = 0, linetype="dashed", color="black") +
  geom_linerange(aes(xmin = lwr, xmax = upr, color="Optimal Visual Testing Intervals (80%)", linewidth="Optimal Visual Testing Intervals (80%)")) +
  geom_linerange(aes(xmin = lwr_add, xmax = upr_add, color="Standard Confidence Interval (95%)", linewidth="Standard Confidence Interval (95%)")) +
  geom_point(size=3) +
  scale_color_manual(values=c("gray75", "black")) + 
  scale_linewidth_manual(values=c(3.5, .5)) + 
  labs(x = "Coefficient Estimate", y = "",
       color = "", linewidth="") +
  theme_classic() + 
  theme(legend.position="top", 
        legend.box="vertical")

## -----------------------------------------------------------------------------
data(TitanicSurvival, package="carData")
NewTitanicSurvival <- TitanicSurvival %>% 
  mutate(ageCat = case_when(age <= 10 ~ "0-10",
                            age > 10 & age <=18 ~ "11-18",
                            age > 18 & age <=30 ~ "19-30",
                            age > 30 & age <=40 ~ "31-40",
                            age > 40 & age <=50 ~ "41-50",
                            age >50 ~ "51+"))
# The model
model3 <- glm(survived~sex*ageCat+passengerClass,
              data=NewTitanicSurvival,family = binomial(link="logit")) 

## -----------------------------------------------------------------------------
## all pairwise tests
mes <- avg_predictions(model3, 
                variables = list(ageCat = levels(NewTitanicSurvival$ageCat),
                                 sex=levels(NewTitanicSurvival$sex)))

all_vt <- viztest(mes, include_zero=FALSE)
plot_data <- plot(all_vt, make_plot=FALSE) %>% 
  dplyr::select(est, lwr, upr, lwr_add, upr_add) 
mes <- mes %>% 
  as.data.frame() %>% 
  left_join(plot_data, by = join_by(estimate == est)) 

## ----titanic1, fig.height=8, fig.width=6, out.width="75%", fig.cap="Predicted Probabilities of Survival on the Titanic with Inferential Confidence Intervals (91%, wide gray bars) and standard 95% confidence interval (thin black line).  If the gray bars overlap, it indicates the pair of estimates whose intervals overlap are not statistically different from each other.  If the gray bars do not overlap, the pair of estimates are statistically different from each other. Tests were done using normal-theory on the response scale with standard errors computed by the delta method.  Comparisons both within and across panels are valid. "----
ggplot(mes, aes(x=estimate, y = ageCat)) + 
  geom_linerange(aes(xmin = lwr, xmax = upr, color="Optimal Visual Testing Intervals (91%)", linewidth="Optimal Visual Testing Intervals (91%)")) +
  geom_linerange(aes(xmin = lwr_add, xmax = upr_add, color="Standard Confidence Intervals (95%)", linewidth="Standard Confidence Intervals (95%)")) +
  geom_point(size=3) + 
  scale_color_manual(values=c("gray75", "black")) +
  scale_linewidth_manual(values=c(3.5, .5)) + 
  facet_grid(sex~.) +
  theme_classic() + 
  theme(legend.position="top") + 
  labs(x="Predicted Pr(Survival)", y="Age Category",
       color="", linewidth="")

## -----------------------------------------------------------------------------
## all pairwise tests
mes_link <- avg_predictions(model3, 
                variables = list(ageCat = levels(NewTitanicSurvival$ageCat),
                                 sex=levels(NewTitanicSurvival$sex)), 
                type="link")

all_vt <- viztest(mes_link, include_zero=FALSE)

## -----------------------------------------------------------------------------
plot_data <- plot(all_vt, make_plot=FALSE) %>% 
  dplyr::select(est, lwr, upr, lwr_add, upr_add) 
mes_resp <- mes_link %>% 
  as.data.frame() %>% 
  left_join(plot_data, by = join_by(estimate == est)) %>%
  mutate(across(c(estimate, lwr, upr, lwr_add, upr_add), plogis))

## ----titanic2, fig.height=8, fig.width=6, out.width="75%", fig.cap="Predicted Probabilities of Survival on the Titanic with Inferential Confidence Intervals (88%, wide gray bars) and standard 95% confidence interval (thin black line).  If the gray bars overlap, it indicates the pair of estimates whose intervals overlap are not statistically different from each other.  If the gray bars do not overlap, the pair of estimates are statistically different from each other. Tests were done using normal-theory on the link scale.  Comparisons both within and across panels are valid. "----
ggplot(mes_resp, aes(x=estimate, y = ageCat)) + 
  geom_linerange(aes(xmin = lwr, xmax = upr, color="Optimal Visual Testing Intervals (88%)", linewidth="Optimal Visual Testing Intervals (88%)")) +
  geom_linerange(aes(xmin = lwr_add, xmax = upr_add, color="Standard Confidence Intervals (95%)", linewidth="Standard Confidence Intervals (95%)")) +
  geom_point(size=3) + 
  scale_color_manual(values=c("gray75", "black")) +
  scale_linewidth_manual(values=c(3.5, .5)) + 
  facet_grid(sex~.) +
  theme_classic() + 
  theme(legend.position="top") + 
  labs(x="Predicted Pr(Survival)", y="Age Category",
       color="", linewidth="")


## -----------------------------------------------------------------------------
mes_sim <- avg_predictions(model3, 
                variables = list(ageCat = levels(NewTitanicSurvival$ageCat),
                                 sex=levels(NewTitanicSurvival$sex))) %>% 
  inferences(method="simulation", R=1000)

## -----------------------------------------------------------------------------
post <- posterior_draws(mes_sim) %>% dplyr::select(drawid, draw, ageCat, sex) %>% 
  pivot_wider(names_from = "drawid", values_from = "draw", names_prefix="draw")
sim <- post %>% dplyr::select(starts_with("draw")) %>% as.matrix() %>% t()
colnames(sim) <- paste0("b", 1:ncol(sim))
sim_vt_data <- make_vt_data(sim, type="sim")
all_vt_sim <- viztest(sim_vt_data, include_zero=FALSE)
all_vt_sim

## -----------------------------------------------------------------------------
plot_data <- plot(all_vt_sim, make_plot=FALSE) %>% 
  dplyr::select(vbl, est, lwr, upr, lwr_add, upr_add) %>% 
  arrange(as.numeric(gsub("b", "", vbl)))
mes_sim <- mes_sim %>% 
  as.data.frame() %>% 
  bind_cols(plot_data) 

## ----titanic3, fig.height=8, fig.width=6, out.width="75%", fig.cap="Predicted Probabilities of Survival on the Titanic with Inferential Credible Intervals (89%, wide gray bars) and standard 95% credible interval (thin black line).  If the gray bars overlap, it indicates the pair of estimates whose intervals overlap are not credibly different from each other.  If the gray bars do not overlap, the pair of estimates are credibly different from each other. Tests were done using quasi-Bayesian simulation with quantile credible intervals on the response scale. Comparisons both within and across panels are valid. "----
ggplot(mes_sim, aes(x=estimate, y = ageCat)) + 
  geom_linerange(aes(xmin = lwr, xmax = upr, color="Optimal Visual Testing Intervals (89%)", linewidth="Optimal Visual Testing Intervals (89%)")) +
  geom_linerange(aes(xmin = lwr_add, xmax = upr_add, color="Standard Credible Intervals (95%)", linewidth="Standard Credible Intervals (95%)")) +
  geom_point(size=3) + 
  scale_color_manual(values=c("gray75", "black")) +
  scale_linewidth_manual(values=c(3.5, .5)) + 
  facet_grid(sex~.) +
  theme_classic() + 
  theme(legend.position="top") + 
  labs(x="Predicted Pr(Survival)", y="Age Category",
       color="", linewidth="")

## -----------------------------------------------------------------------------
sex_vals <- c("female", "male")

## -----------------------------------------------------------------------------
preds <- lapply(sex_vals, \(s){
  avg_predictions(model3, newdata=datagrid(sex = s, grid_type="counterfactual"), 
                  variables="ageCat")
})

## -----------------------------------------------------------------------------
names(preds) <- sex_vals
group_vt <- lapply(preds, \(p){
  est <- p$estimate
  names(est) <- p$ageCat
  d <- make_vt_data(est, vcov(p), type="est_var")
  viztest(d, include_zero=FALSE)
})

## -----------------------------------------------------------------------------
group_vt

## -----------------------------------------------------------------------------
plot_data <- lapply(names(group_vt), \(g){
  d <- plot(group_vt[[g]], make_plot=FALSE, level = ifelse(g == "female", "median", "ce")) %>% 
    dplyr::select(vbl, est, lwr, upr, lwr_add, upr_add) 
  d
})
names(plot_data) <- sex_vals
plot_data <- bind_rows(plot_data, .id="sex") %>% 
  rename(ageCat = vbl)

preds <- preds %>% 
  bind_rows(.id="sex") %>% 
  as.data.frame() %>% 
  left_join(plot_data)

## ----titanic4, fig.height=8, fig.width=6, out.width="75%", fig.cap="Predicted Probabilities of Survival on the Titanic with Inferential Confidence Intervals (wide gray bars, 91% for males and 73% for females) and standard 95% confidence interval (thin black line).  If the gray bars overlap, it indicates the pair of estimates whose intervals overlap are not statistically different from each other.  If the gray bars do not overlap, the pair of estimates are statistically different from each other. Tests were done using normal-theory on the response scale with standard errors computed by the delta method.  Intervals were built so comparisons within, but not across panels are valid. "----
ggplot(preds, aes(x=estimate, y = ageCat)) + 
  geom_linerange(aes(xmin = lwr, xmax = upr, color="Optimal Visual Testing Intervals (91%)", linewidth="Optimal Visual Testing Intervals (91%)")) +
  geom_linerange(aes(xmin = lwr_add, xmax = upr_add, color="Standard Confidence Intervals (95%)", linewidth="Standard Confidence Intervals (95%)")) +
  geom_point(size=3) + 
  scale_color_manual(values=c("gray75", "black")) +
  scale_linewidth_manual(values=c(3.5, .5)) + 
  facet_grid(sex~.) +
  theme_classic() + 
  theme(legend.position="top") + 
  labs(x="Predicted Pr(Survival)", y="Age Category",
       color="", linewidth="")