Generate heatmaps of mixed effect model estimates using a model of poverty regressed on age, gender and random effects for education and sexual orientation by gender

NHANESre.R

library(NHANES)
data(NHANES)
library(tidyverse)
library(extrafont)
library(lme4)
library(RColorBrewer)

NHANES <- mutate(NHANES, 
                 SexOrientation = fct_relevel(SexOrientation, "Heterosexual"),
                 Education      = fct_relevel(Education, "High School"))

model  <- lmer(data = NHANES, Poverty ~ Age + Gender + (Education + SexOrientation + 1 | Gender) )

pal    <- brewer.pal(n=3, name = "RdYlBu")

labels <- c(SexOrientation = "Sexual\nOrientation",
            Education      = "Education",
            Intercept      = "")

if (any(grepl(pattern = "^Lato$", x = fonts()))){
    base_font <- "Lato"
} else {
    base_font <- "Arial Narrow"
}

p <- coef(model)$Gender %>%
    rownames_to_column(var = "Gender") %>%
    select(-Age, -Gendermale) %>% 
    gather(key, value, -`Gender`) %>% 
    mutate(key = sub(pattern = "Education", replacement = "Education:", x = key),
           key = sub(pattern = "SexOrientation", replacement = "SexOrientation:", key)) %>%
    separate(key, into = c("key", "label"), ":") %>%
    mutate(label = ifelse(key == "(Intercept)", "Intercept", label),
           key = ifelse(key == "(Intercept)", "Intercept", key)) %>%
    mutate(label = factor(label, 
                          levels = c("Homosexual",                                    
                                     "Bisexual",
                                     "Intercept",
                                     "8th Grade",
                                     "9 - 11th Grade",
                                     "Some College",
                                     "College Grad")),
           key = factor(key, levels = c("Education", "Intercept", "SexOrientation"))) %>%
    ggplot(data = ., aes(x = Gender, y = label)) +
    geom_tile(aes(fill = value)) +
    scale_fill_gradient2(low = pal[1], mid = pal[2], high = pal[3], midpoint = 0,
                         name = "Median value of\nrandom effect") +
    theme_minimal() +
    facet_grid(key ~ ., scales = "free_y", space = "free", 
               labeller = labeller(key = labels)) + 
    theme(text = element_text(family = base_font),
          axis.title.y = element_blank(),
          panel.grid = element_blank()) +
    labs(title = "Ratio of family income to poverty guidelines\nRandom effects model fit to NHANES data",
         subtitle = "Adjusted for age and gender\nIntercept = heterosexual male with high school education")

ggsave(filename = "NHANESre.png", plot = p,  width = 5, height = 5, units = "in", dpi = 300)