# Clear work space
rm(list = ls())

# Load packages
suppressPackageStartupMessages({
    library(tidyverse)
    library(ggplot2)
    library(ggpubr)
    library(Hmisc)
    library(merTools)
})

# Load preprocessed data with invMR
dat <- readRDS(here::here("Data/data_imr.rds"))

# Load fitted meta-analytic models
ls3 <- readRDS(here::here("Data/ls3.rds"))
hp3 <- readRDS(here::here("Data/hp3.rds"))

1 Plot well-being by day of week

# M and SD by day of the week
mns <- dat %>%
        group_by(dow) %>%
        summarise(n = n(),
                  mnls = wtd.mean(ls, wgt),
                  mnhp = wtd.mean(happy, wgt),
                  sdls = wtd.var(ls, wgt)^.5,
                  sdhp = wtd.var(happy, wgt)^.5,
                  .groups = "keep") %>%
        pivot_longer(c(mnls, mnhp, sdls, sdhp),
                     names_to = c(".value", "name"),
                     names_pattern = "^(.{2})(.{2})$") %>%
        mutate(se = sd / n^.5,
               name = factor(name, c("ls", "hp")))

# Histogram for mean well-being ratings
ggplot(mns, aes(x = dow, y = mn, fill = name)) +
    geom_bar(position = position_dodge(), stat = "identity") +
    #scale_fill_grey(start = 0.8, end = 0.6) +
    # geom_errorbar(aes(ymin = mn - se, ymax = mn + se),
    #               size = .3,    # Thinner lines
    #               width = .2,
    #               position = position_dodge(.9)) +
    geom_text(position = position_dodge(width = 1),
              size = 3,
              aes(y = 8.0, label = round(mn, 2)),
              vjust = 0) +
    geom_text(position = position_dodge(width = 1),
              size = 3,
              aes(y = 7.6, label = paste0("(", round(sd, 2), ")")),
              vjust = 0) +
    geom_text(position = position_dodge(width = 1),
              size = 3,
              aes(y = .5, label = format(n, big.mark = ","), fill = NULL),
              vjust = 0) +
    ylim(0, 10) +
    xlab("Day of the week") +
    ylab("Mean subjective well-being rating") +
    scale_fill_manual(values = c("gray73", "gray53"),
                      name = "",
                      labels = c("Life satisfaction", "Happiness")) +
    ggsave(here::here("Plots/Figure 1.png"),
           device = "png",
           width = 200,
           height = 100,
           units = "mm",
           dpi = 300)

rm(mns)

2 Plot random effects

2.1 Empirical Bayes estimates

# EB for life satisfaction
if (file.exists(here::here("Data/lsre.rds"))) {
    lsre <- readRDS(here::here("Data/lsre.rds"))
} else {
    lsre <- REsim(ls3, 100, seed = 243)
    saveRDS(lsre, here::here("Data/lsre.rds"))
}

# EB for happinness
if (file.exists(here::here("Data/hpre.rds"))) {
    hpre <- readRDS(here::here("Data/hpre.rds"))
} else {
    hpre <- REsim(hp3, 100, seed = 243)
    saveRDS(hpre, here::here("Data/hpre.rds"))
}

# Descriptives for life satisfaction
lsre %>%
    filter(term == "isSu2") %>%
    summarise(lb = quantile(median, .05),
              ub = quantile(median, .95),
              min = min(median),
              max = max(median))
# Descriptives for happiness
hpre %>%
    filter(term == "isSu2") %>%
    summarise(lb = quantile(median, .05),
              ub = quantile(median, .95),
              min = min(median),
              max = max(median))

2.2 Plot for Sunday

There are over 30,000 random effects for interviewers. Because these cannot be easily displayed in a single plot, a random sample of 100 interviewers is drawn and plotted.

# Select interviewers for life satisfaction
set.seed(342)
lsre_su <- lsre %>%
             filter(term == "isSu2") %>%
             slice_sample(n = 100, replace = FALSE) %>%
             arrange(median) %>%
             mutate(groupID = factor(groupID, groupID))

# Select interviewers for happiness
set.seed(342)
hpre_su <- hpre %>%
            filter(term == "isSu2") %>%
            slice_sample(n = 100, replace = FALSE) %>%
            arrange(median) %>%
            mutate(groupID = factor(groupID, groupID))

# Create plot for life satisfaction
p1 <- ggplot(lsre_su, aes(x = groupID, y = median)) +
        geom_hline(yintercept = 0, size = 1, color = "gray62") +
        geom_point() +
        geom_errorbar(aes(ymin = median - sd, ymax = median + sd),
                      size = .3,    # Thinner lines
                      width = .2) +
        ylab("Random effect") +
        xlab("Interviewer") +
        ggtitle("Life satisfaction") +
        scale_x_discrete(expand = expansion(mult = c(.05, .05))) +
        scale_y_continuous(breaks = seq(-.2, .2, .1), limits = c(-.25, .25)) +
        theme(axis.line.x = element_blank(),
              axis.line.y = element_line(),
              axis.text.x = element_blank(),
              axis.ticks.x = element_blank(),
              #axis.title.x = element_blank(),
              panel.background = element_blank(),
              panel.border = element_blank(),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              plot.title = element_text(lineheight = .8, face = "bold"))

# Create plot for happiness
p2 <- ggplot(hpre_su, aes(x = groupID, y = median)) +
        geom_hline(yintercept = 0, size = 1, color = "gray62") +
        geom_point() +
        geom_errorbar(aes(ymin = median - sd, ymax = median + sd),
                      size = .3,    # Thinner lines
                      width = .2) +
        ylab("Random effect") +
        xlab("Interviewer") +
        scale_x_discrete(expand = expansion(mult = c(.05, .05))) +
        scale_y_continuous(breaks = seq(-.2, .2, .1), limits = c(-.25, .25)) +
        ggtitle("Happiness") +
        theme(axis.line.x = element_blank(),
              axis.line.y = element_line(),
              axis.text.x = element_blank(),
              axis.ticks.x = element_blank(),
              #axis.title.x = element_blank(),
              panel.background = element_blank(),
              panel.border = element_blank(),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              plot.title = element_text(lineheight = .8, face = "bold"))

# Save combined plot
ggarrange(p1, p2, ncol = 2, nrow = 1) +
    ggsave(here::here("Plots/Figure 2.png"),
           device = "png",
           width = 200,
           height = 100,
           units = "mm",
           dpi = 300)

rm(p1, p2)