inst/doc/compact_letter_displays.R

In VizTest: Optimal Confidence Intervals for Visual Testing

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

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

## -----------------------------------------------------------------------------
library(VizTest)
library(multcomp)
chick_gl <- glht(chick_mod, linfct=mcp(feed = "Tukey"))

## -----------------------------------------------------------------------------
lmat_gl <- get_letters(chick_gl, test=adjusted(type="none"))
lmat_gl

## -----------------------------------------------------------------------------
library(emmeans)
chick_em <- emmeans(chick_mod, "feed")
lmat_em <- get_letters(chick_em, adjust="none")
lmat_em

## -----------------------------------------------------------------------------
library(marginaleffects)
chick_preds <- predictions(chick_mod, variables="feed", by="feed")
chick_b <- coef(chick_preds)
names(chick_b) <- chick_preds$feed
lmat_df <- get_letters(list(est=chick_b, var=vcov(chick_preds)), adjust="none")
lmat_df

## -----------------------------------------------------------------------------

library(dplyr)
chick_preds_dat <- chick_preds %>% 
  as.data.frame() %>% 
  rename(x = feed, predicted=estimate)

## ----chick-letter, fig.width=6, fig.height=4, fig.align="center", fig.cap="Predicted chick weights by feed type with 95% confidence intervals and a compact letter display"----
library(psre)
library(ggplot2)
letter_plot(chick_preds_dat, lmat_gl)

## -----------------------------------------------------------------------------
## Load Data
data(Ornstein, package="carData")

## Estimate Model
imod <- glm(interlocks ~ log2(assets) + sector + nation, data=Ornstein, 
    family=poisson)

## Generate Predictions
ipreds <- predictions(imod, variables = "sector", by = "sector")

## -----------------------------------------------------------------------------
## merge predictions and data
Ornstein <- Ornstein %>% 
  left_join(ipreds %>% 
              as.data.frame() %>% 
              dplyr::select(sector, estimate)) %>% 
  ## re-order factor 
  mutate(sector = reorder(sector, estimate, mean)) %>% 
  dplyr::select(-estimate)

## update model
imod <- update(imod, data=Ornstein)

## make predictions again
ipreds <- predictions(imod, variables = "sector", by = "sector")

## -----------------------------------------------------------------------------
est <- ipreds$estimate
names(est) <- ipreds$sector
i_lets <- get_letters(list(est=est, var = vcov(ipreds)), adjust="none")

## ----orn-lets, fig.width=6, fig.height=6, fig.align="center", fig.cap="Predicted chick weights by feed type with 95% confidence intervals and a compact letter display"----
ipreds <- ipreds %>% 
  as.data.frame() %>% 
  rename(x = sector, predicted=estimate)
letter_plot(ipreds, i_lets)

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

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

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

## ----ucb-lets, fig.width=10, fig.height=6, fig.align="center", fig.cap="Probability of being admitted to UC Berkeley by Gender and Department", out.width="100%"----
## make predictions for both variables
ucb_pred <- predictions(ucb_mod, variables=c("Gender", "Dept"), by=c("Gender", "Dept"))

## save coefficients and variance-covariance matrix
ucb_v <- vcov(ucb_pred)
ucb_b <- coef(ucb_pred)

## make joint-name of the group and y-axis variables (in this case, Gender and Dept)
## We also rename the identifying variable "x" and the estimate variable to "predicted" 
## as is expected by letter_plot().
ucb_pred <- ucb_pred %>% 
  as.data.frame() %>% 
  mutate(gd = sprintf("%s-%s", Gender, Dept)) %>% 
  rename(x=gd, predicted=estimate)

## set the names of the estimates as the "x" variable. 
names(ucb_b) <- ucb_pred$x

## Get the letters
ucb_lets <- get_letters(list(est=ucb_b, var=ucb_v), adjust="none")

## subset the letters by looking for "Male" in the rownames
m_lets <- ucb_lets[grepl("Male", rownames(ucb_lets)), ]

## subset the letters by looking for "Female" in the rownames
f_lets <- ucb_lets[grepl("Female", rownames(ucb_lets)), ]

## Subset the predictions by gender
m_preds <- ucb_pred %>% filter(Gender == "Male")
f_preds <- ucb_pred %>% filter(Gender == "Female")

## Make the letter plots - when done, remove "Male" and "Female" from the y-axis labels. 
lp_m <- letter_plot(m_preds, m_lets) +
  scale_y_discrete(labels = ~gsub("Male-", "", .x)) + 
  facet_wrap(~"Male")

lp_f <- letter_plot(f_preds, f_lets) +
  scale_y_discrete(labels = ~gsub("Female-", "", .x)) + 
  facet_wrap(~"Female")

## Put the plots together. 
library(patchwork)
(lp_m + lp_f) + plot_layout(ncol=2)



## ----ucb-lets2, fig.width=10, fig.height=6, fig.align="center", fig.cap="Probability of being admitted to UC Berkeley by Gender and Department", out.width="100%"----
## Estimate subset models. 
ucb_mod_m <- glm(cbind(Admitted, Rejected) ~ Dept, 
               data=subset(ucb_dat, Gender == "Male"), family=binomial)
ucb_mod_f <- glm(cbind(Admitted, Rejected) ~ Dept, 
               data=subset(ucb_dat, Gender == "Female"), family=binomial)

## generate subset model predictions
ucb_pred_m <- predictions(ucb_mod_m, variables="Dept", by="Dept")
ucb_pred_f <- predictions(ucb_mod_f, variables="Dept", by="Dept")

## save and name coefficients; save variance-covariance matrices
ucbb_m <- coef(ucb_pred_m)
ucbb_f <- coef(ucb_pred_f)
names(ucbb_m) <- names(ucbb_f) <- ucb_pred_m$Dept
ucbv_m <- vcov(ucb_pred_m)
ucbv_f <- vcov(ucb_pred_f)

## get letters for each set of predictions and variances
ucb_lets_m <- get_letters(list(est=ucbb_m, var=ucbv_m), adjust="none")
ucb_lets_f <- get_letters(list(est=ucbb_f, var=ucbv_f), adjust="none")

## rename relevant variabels. 
ucb_pred_m <- ucb_pred_m %>% 
  as.data.frame() %>% 
  rename("x" = "Dept", "predicted" = "estimate")

ucb_pred_f <- ucb_pred_f %>% 
  as.data.frame() %>% 
  rename("x" = "Dept", "predicted" = "estimate")

## make plots
lp_m <- letter_plot(ucb_pred_m, ucb_lets_m) + facet_wrap(~"Male")
lp_f <- letter_plot(ucb_pred_f, ucb_lets_f) + facet_wrap(~"Female")

## print plots together 
(lp_m + lp_f) + plot_layout(ncol=2)