## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", 
  message = FALSE, 
  warning = FALSE
)
library(dplyr)

## -----------------------------------------------------------------------------
data(chickwts)
chickwts$feed <- reorder(chickwts$feed, 
                         chickwts$weight, 
                         FUN=mean)
chick_mod <- lm(weight~ feed, data=chickwts)

## -----------------------------------------------------------------------------
library(emmeans)
chick_emm <- emmeans(chick_mod, ~ feed)
pwpm(chick_emm, 
     by = NULL, 
     adjust = "none", 
     type = "response")

## ----chick-fp, fig.height=6.5, fig.width=6, out.width="75%", fig.align="center", fig.cap="Heatmap of Pairwise Differences in Chick Weights by Feed Type using `factorplot()`."----
library(factorplot)
chick_fp <- factorplot(chick_mod, factor.var="feed", order="size")
plot(chick_fp)

## -----------------------------------------------------------------------------
chick_df1 <- fp_to_df(chick_fp) 
head(chick_df1)

## ----chick-fp2, fig.height=6.5, fig.width=6, out.width="75%", fig.align="center", fig.cap="Custom Heatmap of Pairwise Differences in Chick Weights by Feed Type Encoding Difference with Color. White asterisks indicate statistical significance at the 0.05 level (two-tailed)."----
ggplot(chick_df1, aes(x=row, y=column)) + 
  geom_tile(aes(fill = difference), color="black", linewidth=.25) + 
  scale_fill_viridis_c(option = "D", na.value = "transparent") + 
  geom_text(aes(label = ifelse(p_value < .05, "*", "")), 
            color = "white", size = 6) + 
  geom_text(data = filter(chick_df1, !is.na(estimate)), aes(label = sprintf("%.2f", estimate))) + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle=45, hjust=1), 
        legend.position = "top") + 
  labs(x="", y="", fill="Difference")

## -----------------------------------------------------------------------------
library(marginaleffects)
chick_preds <- predictions(chick_mod, variables="feed", by="feed")
chick_b <- coef(chick_preds)
names(chick_b) <- chick_preds$feed
chick_fp2 <- factorplot(chick_b, var=vcov(chick_preds), order="size")
chick_df1b <- fp_to_df(chick_fp2) 
head(chick_df1b)

## ----chick-fp3, fig.height=6.5, fig.width=6, out.width="75%", fig.align="center", fig.cap="Custom Heatmap of Pairwise Differences in Chick Weights by Feed Type Encoding Difference with Color, Absolute Estimates.White asterisks indicate statistical significance at the 0.05 level (two-tailed)."----
ggplot(chick_df1b, aes(x=row, y=column)) + 
  geom_tile(aes(fill = difference), color="black", linewidth=.25) + 
  scale_fill_viridis_c(option = "D", na.value = "transparent") + 
  geom_text(aes(label = ifelse(p_value < .05, "*", "")), 
            color = "white", size = 6) + 
  geom_text(data = filter(chick_df1b, !is.na(estimate)), aes(label = sprintf("%.2f", estimate))) + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle=45, hjust=1), 
        legend.position = "top") + 
  labs(x="", y="", fill="Difference")

## -----------------------------------------------------------------------------
chick_df2 <- fp_to_df(chick_fp2, type="both_tri")
head(chick_df2)

## ----chick-fp4, fig.height=7, fig.width=6, out.width="75%", fig.align="center", fig.cap="Custom Heatmap of Pairwise Differences in Chick Weights by Feed Type Encoding Difference and p-values with Color in Alternating Triangles."----
library(ggnewscale)
chick_df2 <- chick_df2 %>% 
  mutate(row = factor(as.character(row), levels=rev(levels(chickwts$feed))), 
         column = factor(as.character(column), levels=levels(chickwts$feed)), 
         significant = as.factor(ifelse(p_value < .05, "Significant", "Not Significant")))
ggplot(chick_df2, aes(x=row, y=column)) + 
  geom_tile(fill="transparent", color="black", linewidth=.25) + 
  geom_tile(data= filter(chick_df2, type=="difference"), 
            aes(fill = difference), color="black", linewidth=.25) +
  scale_fill_viridis_c(option = "D", na.value = "transparent") + 
  labs(fill = "Difference") + 
  new_scale_fill() +
  geom_tile(data= filter(chick_df2, type=="pval"), 
            aes(fill = significant), color="black", linewidth=.25) +
  scale_fill_manual(values=c("white", "gray50")) + 
  labs(fill="p-value", x="", y="") + 
  geom_text(data = filter(chick_df2, !is.na(estimate)), aes(label = sprintf("%.2f", estimate))) +
  theme_classic() + 
  theme(legend.position = "top", legend.box="vertical") 
  

## -----------------------------------------------------------------------------
data(UCBAdmissions)

ucb_dat <- as.data.frame(UCBAdmissions) %>% 
  tidyr::pivot_wider(names_from = "Admit", values_from = "Freq") %>% 
  mutate(Dept = as.factor(Dept))

ucb_mod <- glm(cbind(Admitted, Rejected) ~ Gender * Dept, 
               data=ucb_dat, family=binomial)

## -----------------------------------------------------------------------------
gen_list <- c("Male", "Female")
pred_list <- lapply(gen_list, \(g)predictions(ucb_mod, 
            variables = list("Dept" = levels(ucb_dat$Dept), 
                             "Gender" = g), 
            by="Dept"))
names(pred_list) <- gen_list
df_list <- lapply(pred_list, \(x){
  b <- coef(x)
  names(b) <- x$Dept
  v <- vcov(x)
  fp <- factorplot(b, var=v)
  fp_to_df(fp, type="upper_tri")
})
ucb_df <- bind_rows(df_list, .id="Gender")
ucb_df <- ucb_df %>% 
  mutate(across(c(row, column), ~as.factor(as.character(.x))))

## ----ucb-mf, fig.height=6, fig.width=10, out.width="100%", fig.align="center", fig.cap="Heatmaps of Predicted Admission Rates to UCB by Gender and Department with Pairwise Differences Encoded with Color. White asterisks indicate statistical significance at the 0.05 level (two-tailed)."----
ggplot(ucb_df, aes(x=row, y=forcats::fct_rev(column), fill=difference)) + 
  geom_tile(color="black", linewidth=.25) + 
  geom_text(aes(label = ifelse(p_value < .05, "*", "")), 
            color = "white", size = 6) +
  geom_text(data=filter(ucb_df, !is.na(estimate)), 
            aes(label = sprintf("%.2f", estimate))) +
  scale_fill_viridis_c(na.value = "transparent") + 
  theme_classic() + 
  facet_wrap(~Gender, nrow=1) + 
  theme(legend.position="top") + 
  labs(x="", y="", fill="Difference")