library(lme4)
#> Loading required package: Matrix
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(ggplot2)
m <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

# Make a dataframe with the fitted effects
ab_lines <- coef(m)[["Subject"]] %>% 
  tibble::rownames_to_column("Subject") %>% 
  rename(intercept = `(Intercept)`)
ab_lines
#>    Subject intercept       Days
#> 1      308  253.6637 19.6662579
#> 2      309  211.0065  1.8475828
#> 3      310  212.4449  5.0184060
#> 4      330  275.0956  5.6529547
#> 5      331  273.6653  7.3973914
#> 6      332  260.4446 10.1951153
#> 7      333  268.2455 10.2436615
#> 8      334  244.1725 11.5418620
#> 9      335  251.0714 -0.2848732
#> 10     337  286.2955 19.0955699
#> 11     349  226.1950 11.6407002
#> 12     350  238.3351 17.0814910
#> 13     351  255.9829  7.4520288
#> 14     352  272.2687 14.0032993
#> 15     369  254.6806 11.3395026
#> 16     370  225.7922 15.2897506
#> 17     371  252.2121  9.4791309
#> 18     372  263.7196 11.7513157

ggplot(sleepstudy) + 
  aes(x = Days, y = Reaction) + 
  # Add the lines
  geom_abline(aes(intercept = intercept, slope = Days), data = ab_lines) +
  # Map the points
  geom_point() +
  # Tile
  facet_wrap("Subject")

Update

Now, I’m curious about the shrinkage effect of the model.

ab_lines <- ab_lines %>% 
  tibble::add_column(Model = "partial pooling")

ab_lines2 <- lmList(Reaction ~ Days | Subject, sleepstudy) %>% 
  coef() %>% 
  tibble::rownames_to_column("Subject") %>% 
  rename(intercept = `(Intercept)`) %>% 
  tibble::add_column(Model = "no pooling")

ab_lines3 <- data_frame(
  Subject = unique(sleepstudy$Subject),
  intercept = coef(lm(Reaction ~ Days, sleepstudy))[1],
  Days = coef(lm(Reaction ~ Days, sleepstudy))[2],
  Model = "total pooling"
)

all_lines <- bind_rows(ab_lines, ab_lines2, ab_lines3)

ggplot(sleepstudy) + 
  aes(x = Days, y = Reaction) + 
  # Add the lines
  geom_abline(aes(intercept = intercept, slope = Days, color = Model), data = all_lines) +
  # Map the points
  geom_point() +
  # Tile
  facet_wrap("Subject") + 
  theme(legend.position = c(.9, .1))

ggplot(bind_rows(ab_lines2, ab_lines)) + 
  aes(x = intercept, y = Days, color = Model) + 
  geom_point() + 
  geom_point(data = ab_lines3) + 
  geom_path(aes(group = Subject), 
            arrow = arrow(length = unit(.02, "npc"))) + 
  theme(legend.position = "bottom") + 
  ggtitle("Shrinkage of regression parameters")