coefficient plots, etc.

setup

load packages

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)

coefficient plots

principles

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

nuts and bolts

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

example

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)

coefficient table

stargazer(m3,type="html")


	Dependent variable:

	use

age	-0.046^**
	(0.022)

chY	1.213^***
	(0.210)

I(age2)	-0.006^***
	(0.001)

urbanY	0.787^***
	(0.172)

age:chY	0.066^***
	(0.026)

Constant	-1.341^***
	(0.222)


Observations	1,934
Log Likelihood	-1,177.237
Akaike Inf. Crit.	2,368.475
Bayesian Inf. Crit.	2,407.446

Note:	p<0.1; p<0.05; p<0.01

Also see broom/broom.mixed + huxtable packages

coefficient plot

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

add reference line

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

scaling

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

dwplot without auto-scaling

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

dotwhisker::dwplot(m3, by_2sd = FALSE)

alternative: explicitly scale age parameter

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

plot

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

other alternatives

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

manually tidying

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

GLM for comparison

m3_fixed <-  glm(
    use ~ sc_age * ch + I(sc_age^2) + urban,
    data=Contraception, family=binomial)
dotwhisker::dwplot(list(sc_GLMM=m3_sc, GLM=m3_fixed))+
    geom_vline(xintercept=0,lty=2)+ scale_colour_brewer(palette="Dark2")

dotwhisker limitations

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

alternative: data prep

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

alternative: plot

pd <- position_dodge(width=0.5)
(gg5 <- ggplot(m3_res,aes(x=estimate,y=term,colour=model))
  + geom_pointrange(aes(xmin=conf.low,xmax=conf.high),
                    position=pd)
  + labs(y="")
  + scale_colour_brewer(palette="Dark2") + geom_vline(xintercept=0,lty=2))

reorder terms

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

caterpillar plots

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

caterpillar plots (utility)

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

caterpillar plot

cat_data  <- (m3
  %>% broom.mixed::tidy(m3, effects="ran_vals")
  %>% select(level, term, estimate, std.error)
  %>% caterpillar_levels()
)
ggplot(cat_data, aes(level, estimate,
                     ymin = estimate-2*std.error,
                     ymax = estimate+2*std.error)) +
  geom_pointrange() +
  guides(x = guide_axis(n.dodge = 3))

summary

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

other thoughts

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?

marginal/effects/prediction plots

emmeans

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

as.data.frame → ggplot

e1d <- as.data.frame(e1)
print(gg1 <- ggplot(e1d,
              aes(age, prob, colour=ch, fill=ch,
                  ymin=asymp.LCL, ymax=asymp.UCL)) +
  geom_line() +
  geom_ribbon(colour=NA, alpha=0.2))

effects

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

as.data.frame

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

ggeffects

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

partial effects removal

‘Correct’ data for non-focal terms:

\[ \\beta\_0 + \\sum\_{i \\in \\substack{\\textrm{focal} \\\\ \\textrm{params}}} \\beta\_i x\_i + \\textrm{resid}\_i \]

\[ Y\_i - \\sum\_{i \\in \\substack{\\textrm{non-focal} \\\\ \\textrm{params}}} \\beta\_i x\_i \]

example

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

example

example (cont.)

sleepstudy$re <- remef(fm1, ran = list(Subject = "Days"))
gg3 <- plot(ag) + geom_point(data=sleepstudy,
                      aes(x=Days, y = re, colour=Subject)) +
  colorspace::scale_colour_discrete_qualitative(guide="none")
cowplot::plot_grid(gg2, gg3)

faux: nice contrast names

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)

references

Gelman, Andrew. 2008. “Scaling Regression Inputs by Dividing by Two Standard Deviations.” Statistics in Medicine 27 (15): 2865–73. https://doi.org/10.1002/sim.3107.

Schielzeth, Holger. 2010. “Simple Means to Improve the Interpretability of Regression Coefficients.” Methods in Ecology and Evolution 1: 103–13. https://doi.org/10.1111/j.2041-210X.2010.00012.x.