Differences in slopes or residual variance

library(tidyverse)
library(brms)

Data-generating model

set.seed(201910052)
n <- 1000
people <- tibble(
  id = 1:n,
  days_WFH = rep(0:5, length.out = n),
  parent = id > n/2,
  sleep = rnorm(n, if_else(parent, 6, 8), if_else(parent, 5, 1)),
  RandomVariation = rnorm(n)
)

table(round(people$sleep))

## 
##  -7  -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6   7   8   9  10  11  12 
##   1   3   2   4   6  11  17  17  18  36  28  38  39  89 159 231 153  53  29  17 
##  13  14  15  16  17  18  21 
##  19  11   6   3   2   6   2

xtabs(~ parent + days_WFH, people)

##        days_WFH
## parent   0  1  2  3  4  5
##   FALSE 84 84 83 83 83 83
##   TRUE  83 83 84 84 83 83

table(people$days_WFH)

## 
##   0   1   2   3   4   5 
## 167 167 167 167 166 166

Variance in outcome increases with moderator, predictor variance is lower

DGM

people <- people %>% mutate(
  Productivity = 20 + sleep + if_else(parent, 0.31, 0.1) * days_WFH + 1.75 * RandomVariation
)
ggplot(data = people) + geom_histogram(aes(sleep, fill = parent), position = position_identity(), alpha = 0.4)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(data = people) + geom_histogram(aes(Productivity, fill = parent), position = position_identity(), alpha = 0.4)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Fit LM

model <- brm(Productivity ~ parent * days_WFH, data = people, cores = 4, backend = "cmdstanr", file = "models/productivity_parents_mean")
model2 <- brm(bf(Productivity ~ parent * days_WFH, 
                 sigma ~ parent), data = people, cores = 4, backend = "cmdstanr", file = "models/productivity_parents_mean_sigma")
people[,c("estimate__", "std.error__", "lower__", "upper__")] <- fitted(model, re_formula = NA)

Inspect correlations

(cors <- people %>% group_by(parent) %>%
  summarise(cor = cor(days_WFH, Productivity),
            slope = broom::tidy(lm(Productivity ~ days_WFH))[2, "estimate"][[1]],
            sd_resid = sd(resid(lm(Productivity ~ days_WFH))),
            sd_prod = sd(Productivity),
            mean_prod = mean(Productivity),
            max_js = max(days_WFH),
            mean_js = mean(days_WFH),
            mean_est = mean(estimate__),
            sd_js = sd(days_WFH))
)

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 2 x 10
##   parent   cor slope sd_resid sd_prod mean_prod max_js mean_js mean_est sd_js
##   <lgl>  <dbl> <dbl>    <dbl>   <dbl>     <dbl>  <int>   <dbl>    <dbl> <dbl>
## 1 FALSE  0.131 0.157     2.03    2.05      28.2      5    2.49     28.2  1.71
## 2 TRUE   0.132 0.398     5.10    5.15      27.0      5    2.5      27.0  1.71

Scatter plot

m6 <- ggplot(people, aes(days_WFH, Productivity)) +
  facet_wrap(~ parent, labeller = labeller(parent = c("TRUE" = "parent", "FALSE" = "non-parent"))) +
  geom_point(alpha = 0.1, color = "#6E9181") +
  # geom_line(aes(color = factor(if_else(id > 25, id - 24.5, as.numeric(id)))), stat = "smooth", method = 'lm', alpha = 0.3) +
  geom_smooth(aes(days_WFH, Productivity), method = 'lm', color = 'black') +
  geom_text(aes(max_js - 3, mean_est + 1.1, label = sprintf('b==plain("%.2f")', slope)), data = cors, hjust = 0, parse = TRUE, size = 3) +

  scale_color_viridis_d(guide = F, option = "B")+
  geom_text(aes(label = sprintf('italic("r")==plain("%.2f")', cor)), x= 0.5, y=39, data = cors, hjust = 0, parse = TRUE, size = 3) +

  geom_text(aes(max_js + 0.3, y = mean_prod - 2, label = sprintf('sigma[y]==plain("%.2f")', sd_prod)),, data = cors, hjust = 0, parse = TRUE, size = 3, angle = 90) +
  geom_segment(aes(x = max_js + 0.55, y = mean_prod - sd_prod/2, xend = max_js + 0.55, yend = mean_prod + sd_prod/2), data = cors, size = 1, color = "#219FD8") +

  geom_text(aes(max_js + 0.9, y = mean_prod - 2, label = sprintf('sigma[res(y)]==plain("%.2f")', sd_resid)),, data = cors, hjust = 0, parse = TRUE, size = 3, angle = 90) +
  geom_segment(aes(x = max_js + 1.15, y = mean_prod - sd_resid/2, xend = max_js + 1.15, yend = mean_prod + sd_resid/2), data = cors, size = 1, color = "#219FD8") +

  geom_text(aes(label = sprintf('sigma[x]==plain("%.1f")', sd_js),  x= mean_js - .7, y = 16), data = cors, hjust = 0, parse = TRUE, size = 3) +
  geom_segment(aes(x = mean_js - sd_js/2, y = 15, xend = mean_js + sd_js/2, yend = 15), data = cors, size = 1, color = "#219FD8") +
  

  # coord_cartesian(ylim = c(-3.5, 3.5)) +
  scale_x_continuous("Days worked from home", breaks = 0:5) +
  ylab("Productivity (h)") +
  ggthemes::theme_clean()
m6

## `geom_smooth()` using formula 'y ~ x'

ggsave("plots/parents_slope_vs_corr.pdf", width = 5.5, height = 4)

## `geom_smooth()` using formula 'y ~ x'

ggsave("plots/parents_slope_vs_corr.png", width = 4.5, height = 3)

## `geom_smooth()` using formula 'y ~ x'