## ----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)