library(tidyverse)
## graphics
theme_set(theme_bw()+theme(panel.spacing=grid::unit(0,"lines")))
## modeling/coef plots
library(lme4)
library(broom)
library(broom.mixed)
library(dotwhisker)
library(stargazer) ## coef tables
## predictions/effects plots
library(emmeans)
library(effects)
library(ggeffects)
## also: margins, marginaleffects, sjPlot, ...
library(remef) ## remotes::install_github("hohenstein/remef")
library(faux)

• make it easy to compare magnitudes
• use appropriate scales
• which comparisons do you want to make?
• typically drop intercept (why?)

• old and busted: arm::coefplot, coefplot2
• new hotness:
  • broom[.mixed] + ggplot2
  • dotwhisker

data(Contraception,package="mlmRev")
Contraception <- Contraception %>%
    mutate(ch = factor(livch != 0, labels = c("N", "Y")))
m3 <- glmer(use ~ age * ch + I(age^2) + urban + (1 | urban:district),
            data=Contraception, family=binomial)

gg0 <- dotwhisker::dwplot(m3, by_2sd=TRUE)
print(gg0)

gg1 <- gg0 + geom_vline(xintercept=0,lty=2)
print(gg1)

• need to scale parameters appropriately to compare magnitude
• Gelman (2008); Schielzeth (2010)
• arm::standardize, arm::rescale (handy but inefficient)

binary.inputs: options for standardizing binary variables, default is ‘center’; ‘0/1’ keeps original scale; ‘-0.5,0.5’ rescales 0 as -0.5 and 1 as 0.5; ‘center’ subtracts the mean; and ‘full’ subtracts the mean and divides by 2 sd.

• mixed models: standard deviations have same scale as corresponding parameter

Look at results if we turn off auto-scaling (default)

dotwhisker::dwplot(m3, by_2sd = FALSE)

Contraception <- Contraception %>%
    mutate(sc_age = drop(scale(age)))
m3_sc <- update(m3, 
      . ~ sc_age * ch + I(sc_age^2) + urban + (1 | urban:district))

dotwhisker::dwplot(m3_sc, effects="fixed", by_2sd = FALSE) +
    geom_vline(xintercept=0,lty=2)

• scaling makes units harder to interpret
• separate sets of predictors into facets according to their units
• ??

cc <- broom.mixed::tidy(m3_sc,effects="fixed")
print(cc,digits=3)

## # A tibble: 6 × 6
##   effect term        estimate std.error statistic  p.value
##   <chr>  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 fixed  (Intercept)   -1.34     0.222      -6.04 1.56e- 9
## 2 fixed  sc_age        -0.416    0.198      -2.10 3.57e- 2
## 3 fixed  chY            1.21     0.210       5.79 7.18e- 9
## 4 fixed  I(sc_age^2)   -0.457    0.0690     -6.62 3.51e-11
## 5 fixed  urbanY         0.787    0.172       4.57 4.88e- 6
## 6 fixed  sc_age:chY     0.599    0.231       2.59 9.57e- 3

• sometimes gets “confused” with complex models
• might want to post-process output of tidy() (edit variable/term names, etc.)
• use purrr::map_dfr(list(mod1=, mod2=, .id = "model") to process a list of models and collapse to a tibble

m3_res <- (purrr::map_dfr(list(with_re = m3_sc, no_re=m3_fixed),
                  tidy,
                  conf.int = TRUE,
                  .id = "model")
  %>% mutate(term=fct_inorder(term))
  %>% filter(term!="(Intercept)")
)

(purrr::map_dfr == “run function on elements of a list, bind results into a single data frame, add .id column based on names of the list elements”)

m3_res_order <- mutate(m3_res,term=reorder(term,estimate))
gg5 %+% m3_res_order

• Useful when plotting many parallel/similar coefficients
• e.g.:
  • random effects in a mixed model
  • bioinformatics: genes, microbial species, etc..

Reorders levels, allowing for multiple terms per level

caterpillar_levels <- function(x, order_term = "(Intercept)") {
    lev_order <- (x
        %>% filter(term==order_term)
        %>% group_by(level)
        %>% summarise(across(estimate, mean, na.rm = TRUE))
        %>% arrange(estimate)
        %>% pull(level)
    )
    x <- dplyr::mutate(x, across(level, factor, levels = lev_order))
    return(x)
}

• useful in comparing lots of models;
  models with many parallel parameters
• separate parameter types with facets
• sort by magnitude of estimate

• can often use a wide aspect ratio (don’t need much height)
• may want to replace/improve coefficient names (faux package?)
• Bayesian models: violin plots of posterior samples? Or 50% + 95% intervals?

e1 <- emmeans(m3, ~age*ch, at=list(age=-14:20), type="response")
plot(e1)

plot(ae <- allEffects(m3, xlevels = list(age=51)))

aed <- as.data.frame(ae)[["age:ch"]]
gg1 %+% aed + aes(y=fit, ymin=lower, ymax=upper)

ag <- ggeffect(m3, terms=c("age [all]", "ch"))
plot(ag)

‘Correct’ data for non-focal terms:

\[ \beta_0 + \sum_{i \in \substack{\textrm{focal} \\ \textrm{params}}} \beta_i x_i + \textrm{resid}_i \]

or

\[ Y_i - \sum_{i \in \substack{\textrm{non-focal} \\ \textrm{params}}} \beta_i x_i \]

fm1 <- lmer(Reaction ~ 1 + Days + (1 + Days | Subject), sleepstudy)
ag <- ggeffect(fm1, terms="Days [all]")
gg2 <- plot(ag) + geom_point(data=sleepstudy,
                             aes(x=Days, y = Reaction, colour=Subject)) +
  colorspace::scale_colour_discrete_qualitative(guide="none")

mm <- (mtcars
  %>% mutate(across(cyl, contr_code_treatment))
)
m1 <- lm(mpg ~ disp*cyl, mm)
dwplot(m1, by_2sd=TRUE) + geom_vline(xintercept = 0, lty=2)