You estimated a model, but you probably want to report an `emmeans`.

Whether a paper is observational or experimental, they almost always end with a set of statistical models. But gosh, they can be a bear to map a set of coefficients to your intuition.

An indicative example:

library(tidyverse)
library(magrittr)
library(broom)

titanic <- "https://hbiostat.org/data/repo/titanic3.csv" %>% 
  read_csv

g1 <- glm(
  factor(survived) ~ age * sex * factor(pclass), 
  data = titanic, 
  family = binomial()
)

g1 %>% 
  broom::tidy()

## # A tibble: 12 × 5
##    term                        estimate std.error statistic p.value
##    <chr>                          <dbl>     <dbl>     <dbl>   <dbl>
##  1 (Intercept)                  2.90       1.24       2.34   0.0193
##  2 age                          0.00960    0.0326     0.294  0.769 
##  3 sexmale                     -2.02       1.35      -1.50   0.134 
##  4 factor(pclass)2              0.595      1.53       0.388  0.698 
##  5 factor(pclass)3             -2.61       1.28      -2.03   0.0423
##  6 age:sexmale                 -0.0470     0.0350    -1.34   0.179 
##  7 age:factor(pclass)2         -0.0546     0.0414    -1.32   0.187 
##  8 age:factor(pclass)3         -0.0274     0.0353    -0.776  0.438 
##  9 sexmale:factor(pclass)2     -0.433      1.73      -0.251  0.802 
## 10 sexmale:factor(pclass)3      0.965      1.43       0.675  0.499 
## 11 age:sexmale:factor(pclass)2 -0.0201     0.0500    -0.403  0.687 
## 12 age:sexmale:factor(pclass)3  0.0309     0.0399     0.776  0.438

Interpreting these results is a Gordian Knot!

Other than the constant, only one significant covariate (p < .05),
having to interpret the effects in terms of log odds
Two variables are categorical, one is continuous

Let’s do the legwork of generating some out of sample fitted values, scaling to probabilities.

nd <- expand_grid(
  age = seq(1, 70, by = .25),
  sex = titanic$sex %>% unique,
  pclass = titanic$pclass %>% factor %>% levels
)

nd %>% 
  mutate(
    fit = g1 %>% 
      predict(
        newdata = nd,
        se.fit = F,
        type = "response"
      ),
    se = g1 %>% 
      predict(
        newdata = nd,
        se.fit = T,
        type = "response"
      ) %>% 
      extract2("se.fit"),
    lo = fit %>% 
      subtract(
        se %>% multiply_by(1.96)
      ),
    hi = fit %>% 
      add(
        se %>% multiply_by(1.96)
      )
  ) %>% 
  ggplot(
    aes(
      ymin = lo, ymax = hi, y = fit, x = age, fill = sex, color = sex
    )
  ) +
  geom_ribbon(
    alpha = .5,
    linewidth = .3
  ) +
  geom_line(
    linewidth = .3
  ) +
  facet_wrap(~pclass, nrow = 1) +
  scale_fill_brewer(
    palette = "Set1"
  ) +
  scale_color_brewer(
    palette = "Set1"
  )

Nuts:

.2 difference in survival probability for young first class passengers,
a .7 gender difference for old first class passengers, and
gender difference shrinking by about 3/4rs when the elderly passengers are in the 3rd class.

Is this what you expected when you only had a set of coefficients to eyeball?

Even this doesn’t necessarily address the questions we typically have for a set of results:

Where exactly are the gender contrasts (the survival differences between men and women) significant, and at what ⍺?
What are the average gender or age or class effects, averaging over the other covariates?

These are the questions that we have after we’re done estimating a model. To address them, we use emmeans (for Estimated Marginals Means).

The `emmeans` package

Forestry science and agronomy play a huge role in the development of the modern statistical toolkit, especially at anglosphere universities in the early 20th century. When you do experiments in a lab which is measured in hundreds of square kilometers, only some of which get rain and sun, orthogonal to whatever you’ve randomly assigned, you need an ability to remove these sources of statistical noise to get a better measure of your causal estimand.

Least Squares Means were introduced into SAS in the 1970s, and ported by Russell Lenth (an emeritus statistician at Iowa) to R in 2010. emmeans superseded lsmeans in 2018.

What’s an estimated marginal mean? To the author:

Estimated marginal means… are derived by using a model to make predictions over a regular grid of predictor combinations (called a reference grid). These predictions may possibly be averaged (typically with equal weights) over one or more of the predictors. Such marginally-averaged predictions are useful for describing the results of fitting a model, particularly in presenting the effects of factors.

This package provides an encapsulated and general way to automate the bespoke analysis done above–we estimate a model, and then to better explain how the factors condition the DV. Here’s how we’d use it for the titanic example:

library(emmeans)

g1 %>% 
  emmeans(
    ~ sex | age + pclass, 
    type = "response",
    at = list(
      age = c(1, 15, 35, 55, 75)
    )
  )

## age =  1, pclass = 1:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.948103 0.059413 Inf  6.31e-01   0.99489
##  male   0.699067 0.108440 Inf  4.58e-01   0.86450
## 
## age = 15, pclass = 1:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.954331 0.035011 Inf  8.12e-01   0.99019
##  male   0.579027 0.086011 Inf  4.08e-01   0.73309
## 
## age = 35, pclass = 1:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.962005 0.016671 Inf  9.12e-01   0.98410
##  male   0.394149 0.043725 Inf  3.12e-01   0.48225
## 
## age = 55, pclass = 1:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.968432 0.024116 Inf  8.67e-01   0.99310
##  male   0.235306 0.048557 Inf  1.53e-01   0.34306
## 
## age = 75, pclass = 1:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.973802 0.035055 Inf  7.16e-01   0.99818
##  male   0.127052 0.054604 Inf  5.25e-02   0.27641
## 
## age =  1, pclass = 2:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.969110 0.026341 Inf  8.48e-01   0.99435
##  male   0.717133 0.116967 Inf  4.50e-01   0.88700
## 
## age = 15, pclass = 2:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.943518 0.030051 Inf  8.47e-01   0.98056
##  male   0.345133 0.069192 Inf  2.24e-01   0.48988
## 
## age = 35, pclass = 2:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.871616 0.036734 Inf  7.81e-01   0.92816
##  male   0.052926 0.020965 Inf  2.40e-02   0.11258
## 
## age = 55, pclass = 2:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.733987 0.126098 Inf  4.38e-01   0.90727
##  male   0.005891 0.005082 Inf  1.08e-03   0.03144
## 
## age = 75, pclass = 2:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.528612 0.278144 Inf  1.12e-01   0.90907
##  male   0.000628 0.000848 Inf  4.44e-05   0.00881
## 
## age =  1, pclass = 3:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.567301 0.080839 Inf  4.07e-01   0.71429
##  male   0.311056 0.070648 Inf  1.91e-01   0.46281
## 
## age = 15, pclass = 3:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.505324 0.047301 Inf  4.13e-01   0.59681
##  male   0.219254 0.030950 Inf  1.65e-01   0.28584
## 
## age = 35, pclass = 3:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.416973 0.058496 Inf  3.09e-01   0.53404
##  male   0.124733 0.023414 Inf  8.56e-02   0.17829
## 
## age = 55, pclass = 3:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.333651 0.106176 Inf  1.64e-01   0.56077
##  male   0.067441 0.028314 Inf  2.91e-02   0.14876
## 
## age = 75, pclass = 3:
##  sex    response       SE  df asymp.LCL asymp.UCL
##  female 0.259566 0.142056 Inf  7.61e-02   0.59880
##  male   0.035400 0.024211 Inf  9.06e-03   0.12838
## 
## Unknown transformation "factor": no transformation done 
## Confidence level used: 0.95

Or contrasts (ie, differences between men and women)

g1 %>% 
  emmeans(
    ~ sex | age + pclass, 
    # type = "response",
    at = list(
      age = c(1, 15, 35, 55, 75)
    )
  ) %>% 
  regrid %>% 
  pairs

## age =  1, pclass = 1:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.249 0.1236 Inf   2.014  0.0440
## 
## age = 15, pclass = 1:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.375 0.0929 Inf   4.041  0.0001
## 
## age = 35, pclass = 1:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.568 0.0468 Inf  12.135  <.0001
## 
## age = 55, pclass = 1:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.733 0.0542 Inf  13.522  <.0001
## 
## age = 75, pclass = 1:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.847 0.0649 Inf  13.049  <.0001
## 
## age =  1, pclass = 2:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.252 0.1199 Inf   2.102  0.0356
## 
## age = 15, pclass = 2:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.598 0.0754 Inf   7.932  <.0001
## 
## age = 35, pclass = 2:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.819 0.0423 Inf  19.356  <.0001
## 
## age = 55, pclass = 2:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.728 0.1262 Inf   5.769  <.0001
## 
## age = 75, pclass = 2:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.528 0.2781 Inf   1.898  0.0577
## 
## age =  1, pclass = 3:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.256 0.1074 Inf   2.387  0.0170
## 
## age = 15, pclass = 3:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.286 0.0565 Inf   5.061  <.0001
## 
## age = 35, pclass = 3:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.292 0.0630 Inf   4.638  <.0001
## 
## age = 55, pclass = 3:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.266 0.1099 Inf   2.423  0.0154
## 
## age = 75, pclass = 3:
##  contrast      estimate     SE  df z.ratio p.value
##  female - male    0.224 0.1441 Inf   1.556  0.1198

Whoa, p values, the whole lot! We could even tidy it and plot these!

g1 %>% 
  emmeans(
    consec ~ sex | age + pclass, 
    # type = "response",
    at = list(
      age = c(1, 15, 35, 55, 75)
    ), 
  ) %>% 
  regrid %>% 
  pairs %>% 
  tidy %>% 
  mutate(
    lo = estimate %>% 
      subtract(
        std.error %>% 
          multiply_by(1.96)
      ),
    hi = estimate %>% 
      add(
        std.error %>% 
          multiply_by(1.96)
      )
  ) %>% 
  ggplot(
    aes(
      factor(age), ymin = lo, ymax = hi, y = estimate, color = factor(pclass)
    )
  ) +
  geom_line(
    aes(group = factor(pclass)),
    position = position_dodge(width = .4)
  ) +
  geom_pointrange(
    position = position_dodge(width = .4)
  )

Don’t worry about the details–just luxuriate in the power and compactness of the formula object:

\[ contrast \sim x1 | x2 \]

The left hand side specifies what contrast you’re estimating–all the factor levels against a common baseline (trt.vs.ctrl), all the factor levels against the one that preceded it (consec) or each pairwise set of levels (pairwise), or all the factor levels against the group mean (eff). If you omit the stuff on the left hand side of the tilde, you won’t estimate a contrast.

The right hand side indicates the factors for which you’d like to estimate means and contrasts.

let’s do a fun

Applied Example

TESS grant recipients have their data posted on Roper within 2 years of their collection. Here we’re going to test a set of survey experiments seeing if whites’ racial attitudes condition their Covid-19 public health attitudes. The codebook is here.

library(haven)

t2 <- "https://github.com/thomasjwood/code_lab/raw/main/data/31120363.por" %>% 
  read_por

t2 %<>% 
  modify_if(
    t2 %>% 
      map_lgl(
        ~attr(., "class") %>% 
          str_detect("label") %>% 
          any
      ),
    as_factor
  )

It’s a two by one vignette experiment:

We’re going to check the outcomes on seven outcome variables: Q14:20

t2 %>% 
  select(
    Q14:Q20
  ) %>% 
  map(
    ~attr(., "label")
  )

## $Q14
## [1] "The US should ease up on measures aimed at slowing the spread of the coronavirus"
## 
## $Q15
## [1] "The United States should take measures aimed at slowing the spread of the corona"
## 
## $Q16
## [1] "State and local directives that ask people to 'shelter in place' or to be 'safer"
## 
## $Q17
## [1] "People who participate in large-scale protests and rallies to call for the reope"
## 
## $Q18
## [1] "How important is it that people stay home rather than participating in protests"
## 
## $Q19
## [1] "There should be targeted aid packages, including monetary grants, to Black commu"
## 
## $Q20
## [1] "How important is it for people to wear a mask when coming close to people outsid"

t2[, 
   str_c(
     "Q", 14:20, "_num"
   )
] <- t2 %>% 
  select(
    Q14:Q20
  ) %>% 
  map(
    \(i)
    i %>% 
      as.character %>% 
      case_match(
        c("Strongly Agree",
          "Extremely Important") ~ 5,
        c("Agree",
          "Very Important") ~ 4,
        c("Neither Agree Nor Disagree",
          "Somewhat Important")  ~ 3,
        c("Disagree",
          "Not Very Important") ~ 2,
        c("Strongly Disagree",
          "Not At All Important") ~ 1,
        c("DON'T KNOW", "SKIPPED ON WEB", "REFUSED") ~ NA_integer_,
      ) %>% 
      as.character %>% 
      as.numeric
  )

t3 <- t2 %>% 
  select(
    CASEID, P_PARTYI, EDUC4, P_EXP, ends_with("_num")
  ) %>% 
  mutate(
    P_PARTYI = P_PARTYI %>% 
      str_extract(
        "Democrat|Independent|Republican"
      ) %>% 
      factor(
        c("Democrat", "Independent", "Republican")
      ),
    EDUC4 = EDUC4 %>% 
      case_match(
        c("No HS diploma", "HS graduate or equivalent") ~ "HS or less",
        .default = EDUC4
        ) %>% 
      factor(
        c("HS or less",
          "Some college",
          "BA or above")
      )
    ) %>%  
  pivot_longer(
    names_to = "dv", 
    values_to = "dv_val", 
    ends_with("_num"), 
    values_drop_na = T 
  ) %>% 
  mutate(
    dv = dv %>% 
      str_remove("_num")
  ) %>% 
  left_join(
    tibble(
      dv = str_c("Q", 14:20),
      dv_lab = c(
        "The US should ease up on measures aimed at slowing the spread of the coronavirus soon, in order to open businesses and get the US economy going again, even if that means more people would get coronavirus and could die",
        "The United States should take measures aimed at slowing the spread of the coronavirus while more widespread testing becomes available, even if that means many businesses will have to stay closed",
        "State and local directives that ask people to 'shelter in place' or to be 'safer at home' are a threat to individual rights and freedom.",
        "People who participate in large-scale protests and rallies to call for the reopening of the states should be fined or ticketed, because these protests and rallies may further enable the spread of the coronavirus.",
        "How important is it that people stay home rather than participating in protests and rallies to pressure their governors to reopen their states?",
        "There should be targeted aid packages, including monetary grants, to Black communities that have been disproportionately affected by the coronavirus.",
        "How important is it for people to wear a mask when coming close to people outside of their home (because of the coronavirus outbreak)?"
      )
    )
  ) %>% 
  group_by(dv) %>% 
  nest %>% 
  full_join(
    expand_grid(
      dv = str_c(
        "Q", 14:20
      ),
      frm = c(
        "dv_val ~ P_EXP",
        "dv_val ~ P_EXP * P_PARTYI",
        "dv_val ~ P_EXP * EDUC4",
        "dv_val ~ P_EXP * EDUC4 * P_PARTYI"
      )
    )
  )

With our nice long nested tibble, estimating the 21 models is a snap

t3$mods <- t3$data %>% 
  map2(
    t3$frm,
    \(i, j)
    j %>% 
      lm(i)
  )

Did any of these DVs move with the vignettes?

library(modelsummary)

t3 %>% 
  filter(
    frm %>% equals("dv_val ~ P_EXP")
  ) %>% 
  use_series(mods) %>% 
  set_attr(
    "names", 
    t3 %>% 
      filter(
        frm %>% equals("dv_val ~ P_EXP")
      ) %>% 
      use_series(dv)
  ) %>% 
  modelsummary(
    gof_map = "nobs", 
    stars = T
  )

	Q14	Q15	Q16	Q17	Q18	Q19	Q20
(Intercept)	2.917***	3.199***	2.615***	2.948***	3.540***	2.920***	3.721***
	(0.079)	(0.075)	(0.079)	(0.078)	(0.076)	(0.071)	(0.073)
P_EXPExperimental group B	−0.065	0.003	−0.057	−0.014	0.024	0.050	−0.038
	(0.110)	(0.104)	(0.110)	(0.108)	(0.106)	(0.099)	(0.102)
Num.Obs.	591	588	587	588	590	590	590
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Hmmm. Well gosh. so much for q racial covid backlash. Maybe… party id differences?

t3 %>%
  filter(
    frm %>% 
      equals("dv_val ~ P_EXP * P_PARTYI")
  ) %>%
  use_series(mods) %>%
  set_attr(
    "names",
    t3 %>%
      filter(
        frm %>% equals("dv_val ~ P_EXP")
      ) %>%
      use_series(dv)
  ) %>%
  modelsummary(
    gof_map = "nobs",
    stars = T
  )

	Q14	Q15	Q16	Q17	Q18	Q19	Q20
(Intercept)	2.038***	3.962***	1.792***	3.764***	4.349***	3.528***	4.425***
	(0.108)	(0.105)	(0.109)	(0.112)	(0.108)	(0.102)	(0.108)
P_EXPExperimental group B	0.032	−0.032	−0.021	−0.220	−0.182	0.171	−0.100
	(0.150)	(0.146)	(0.152)	(0.156)	(0.150)	(0.142)	(0.150)
P_PARTYIIndependent	0.637**	−0.450*	0.636**	−0.578**	−0.563**	−0.435*	−0.820***
	(0.201)	(0.195)	(0.205)	(0.208)	(0.203)	(0.190)	(0.201)
P_PARTYIRepublican	1.722***	−1.518***	1.643***	−1.579***	−1.597***	−1.168***	−1.304***
	(0.146)	(0.143)	(0.149)	(0.152)	(0.147)	(0.139)	(0.147)
P_EXPExperimental group B × P_PARTYIIndependent	0.168	−0.280	0.242	0.333	0.322	−0.164	0.395
	(0.286)	(0.278)	(0.291)	(0.297)	(0.288)	(0.272)	(0.286)
P_EXPExperimental group B × P_PARTYIRepublican	−0.260	0.162	−0.200	0.370+	0.387+	−0.188	0.081
	(0.203)	(0.198)	(0.206)	(0.211)	(0.204)	(0.193)	(0.204)
Num.Obs.	565	562	561	562	564	564	564
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

Mmmm I struggle to do mental arithmetic, summing three decimal places. or a dreaded three way interaction!

t3 %>%
  filter(
    frm %>% 
      equals("dv_val ~ P_EXP * EDUC4 * P_PARTYI")
  ) %>%
  use_series(mods) %>%
  set_attr(
    "names",
    t3 %>%
      filter(
        frm %>% equals("dv_val ~ P_EXP * EDUC4 * P_PARTYI")
      ) %>%
      use_series(dv)
  ) %>%
  modelsummary(
    coef_omit = "^(?!.*P_EXPE).*$",
    gof_map = "nobs",
    stars = T
  )

	Q14	Q15	Q16	Q17	Q18	Q19	Q20
P_EXPExperimental group B	0.092	0.275	−0.987*	0.118	0.065	0.810+	0.451
	(0.455)	(0.443)	(0.454)	(0.481)	(0.458)	(0.429)	(0.454)
P_EXPExperimental group B × EDUC4Some college	−0.248	−0.112	0.969+	−0.353	−0.091	−0.405	−0.589
	(0.510)	(0.498)	(0.510)	(0.538)	(0.514)	(0.481)	(0.509)
P_EXPExperimental group B × EDUC4BA or above	0.003	−0.526	1.056*	−0.387	−0.428	−0.937*	−0.576
	(0.504)	(0.491)	(0.504)	(0.532)	(0.508)	(0.476)	(0.503)
P_EXPExperimental group B × P_PARTYIIndependent	0.565	−0.446	1.643*	0.458	−0.176	−1.063+	−0.310
	(0.673)	(0.655)	(0.672)	(0.707)	(0.678)	(0.635)	(0.672)
P_EXPExperimental group B × P_PARTYIRepublican	−0.268	−0.404	1.071+	−0.182	0.230	−0.553	−0.366
	(0.549)	(0.539)	(0.549)	(0.581)	(0.554)	(0.518)	(0.548)
P_EXPExperimental group B × EDUC4Some college × P_PARTYIIndependent	−0.122	−0.027	−1.372+	−0.088	0.700	0.843	0.935
	(0.790)	(0.769)	(0.790)	(0.829)	(0.797)	(0.745)	(0.788)
P_EXPExperimental group B × EDUC4BA or above × P_PARTYIIndependent	−1.104	0.392	−1.990*	−0.300	0.428	1.273	0.518
	(0.859)	(0.836)	(0.857)	(0.901)	(0.865)	(0.810)	(0.857)
P_EXPExperimental group B × EDUC4Some college × P_PARTYIRepublican	0.213	0.320	−1.180+	0.714	0.010	0.058	0.404
	(0.629)	(0.617)	(0.629)	(0.664)	(0.634)	(0.594)	(0.628)
P_EXPExperimental group B × EDUC4BA or above × P_PARTYIRepublican	−0.029	0.906	−1.552*	0.502	0.197	0.507	0.382
	(0.636)	(0.623)	(0.635)	(0.671)	(0.641)	(0.600)	(0.635)
Num.Obs.	565	562	561	562	564	564	564
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

now, the magic of an emmeans call:

t4 <- t3 %>%
  filter(
    frm %>% 
      equals("dv_val ~ P_EXP * EDUC4 * P_PARTYI")
  ) %>%
  use_series(mods) %>%
  map(
    \(i)
    emmeans(
      i, 
      consec ~ P_EXP | EDUC4 + P_PARTYI, 
      options = list(
        adjust = "none"
      ) 
    )
  ) %>% 
  map2(
    str_c("Q", 14:20),
    \(i, j)
    i %>% 
      extract2("contrasts") %>% 
      tidy %>% 
      mutate(
        dv = j
      )
    ) %>% 
  list_rbind %>% 
  na.omit %>% 
  mutate(
    across(
      c(EDUC4:term, dv),
      fct_inorder
    ),
    lo = estimate %>% 
      subtract(
        std.error %>% multiply_by(1.96)
        ),
    hi = estimate %>% 
      add(
        std.error %>% multiply_by(1.96)
        )
    )

which we can plot…

t4 %>% 
  ggplot(
    aes(
      EDUC4, ymin = lo, ymax = hi, y  = estimate, 
    )
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed", 
    linewidth  =.4
  ) +
  geom_pointrange() +
  facet_grid(
    P_PARTYI ~ dv
  ) +
  theme(
    axis.text.x = element_text(
      angle = 45/1.5, hjust = 1, vjust = 1
      )
  ) +
  labs(
    x = "Education",
    y = "Effect of nominating racial difference in Covid"
  )

Just to understand how you can use the formula function to estimate different contrasts

c(
  consec ~  P_PARTYI,
  pairwise ~ P_PARTYI,
  eff ~ P_PARTYI
  ) %>% 
  map(
    \(i)
    
    t3$mods[[2]] %>% 
      emmeans(i) %>% 
      pluck("contrasts")
    )

## [[1]]
##  contrast                 estimate    SE  df t.ratio p.value
##  Independent - Democrat      0.721 0.143 559   5.040  <.0001
##  Republican - Independent    0.872 0.140 559   6.234  <.0001
## 
## Results are averaged over the levels of: P_EXP 
## P value adjustment: mvt method for 2 tests 
## 
## [[2]]
##  contrast                 estimate    SE  df t.ratio p.value
##  Democrat - Independent     -0.721 0.143 559  -5.040  <.0001
##  Democrat - Republican      -1.592 0.102 559 -15.684  <.0001
##  Independent - Republican   -0.872 0.140 559  -6.234  <.0001
## 
## Results are averaged over the levels of: P_EXP 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## 
## [[3]]
##  contrast           estimate     SE  df t.ratio p.value
##  Democrat effect     -0.7711 0.0683 559 -11.292  <.0001
##  Independent effect  -0.0503 0.0880 559  -0.572  0.5677
##  Republican effect    0.8214 0.0661 559  12.433  <.0001
## 
## Results are averaged over the levels of: P_EXP 
## P value adjustment: fdr method for 3 tests

An exercise

Take the variables Q22A:Q22D which measure perceptions of different racial groups’ compliance with covid social distancing mandates. For each DV, regress it on the experimental indicator, and then the experimental indicator interacted with party identification. Report all experimental differences in a figure or a table.

Code Lab-Lab 5 (Estimated Marginal Means `emmeans`)

Thomas Wood (wood.1080@osu.edu)

December 1, 2023

You estimated a model, but you probably want to report an `emmeans`.

The `emmeans` package

Applied Example

An exercise

Code Lab-Lab 5 (Estimated Marginal Means emmeans)

Thomas Wood (wood.1080@osu.edu)

December 1, 2023

You estimated a model, but you probably want to report an emmeans.

The emmeans package

Applied Example

An exercise

Code Lab-Lab 5 (Estimated Marginal Means `emmeans`)

You estimated a model, but you probably want to report an `emmeans`.

The `emmeans` package