# 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"))
# 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)
# 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))
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)