knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE,
                      cache = TRUE)
options(digits = 4)

# load libraries
libraries <- c("rio",          # import export
               "kableExtra",   # kable styling
               "reshape2",     # reshape data wide long
               "cowplot",      # plot formatting
               "mice",         # missingness pattern
               "summarytools", # descriptives
               "flextable",    # reporting tables
               "officer",      # export tables
               "apaTables",    # apa cor table
               "lavaan",       # RICLPM
               "semPlot",      # RICLPM plot
               "lmerTest",     # change analyses
               "tidyverse")    # general wrangling

lapply(libraries, require, character.only = TRUE)

# load the process macro
# https://www.processmacro.org/download.html
source(file = "process.R") 
# list of loaded packages and versions
si <- devtools::session_info()[[2]]
rownames(si) <- NULL
si %>% 
  select(package, loadedversion, date, source) %>% 
  
  #red bold the called packages
  mutate(package = 
           cell_spec(package, 
                     color = ifelse(package %in% libraries, "red", "black"),
                     bold = ifelse(package %in% libraries, TRUE, FALSE))) %>% 
  knitr::kable(escape = F, caption = "All loaded packages. Bolded in red are those loaded explicitly with <code>library()</code>") %>% 
  kable_styling() %>% 
  scroll_box(height = "300px")
All loaded packages. Bolded in red are those loaded explicitly with library()
package loadedversion date source
abind 1.4-5 2016-07-21 CRAN (R 4.1.0)
apaTables 2.0.8 2021-01-04 CRAN (R 4.1.2)
arm 1.12-2 2021-10-15 CRAN (R 4.1.2)
assertthat 0.2.1 2019-03-21 CRAN (R 4.1.0)
backports 1.4.1 2021-12-13 CRAN (R 4.1.2)
base64enc 0.1-3 2015-07-28 CRAN (R 4.1.0)
boot 1.3-28 2021-05-03 CRAN (R 4.1.2)
broom 0.7.12 2022-01-28 CRAN (R 4.1.2)
bslib 0.3.0 2021-09-02 CRAN (R 4.1.1)
cachem 1.0.6 2021-08-19 CRAN (R 4.1.1)
callr 3.7.0 2021-04-20 CRAN (R 4.1.0)
carData 3.0-4 2020-05-22 CRAN (R 4.1.0)
cellranger 1.1.0 2016-07-27 CRAN (R 4.1.0)
checkmate 2.0.0 2020-02-06 CRAN (R 4.1.0)
cli 3.0.1 2021-07-17 CRAN (R 4.1.1)
cluster 2.1.2 2021-04-17 CRAN (R 4.1.2)
coda 0.19-4 2020-09-30 CRAN (R 4.1.0)
codetools 0.2-18 2020-11-04 CRAN (R 4.1.2)
colorspace 2.0-2 2021-06-24 CRAN (R 4.1.0)
corpcor 1.6.10 2021-09-16 CRAN (R 4.1.1)
cowplot 1.1.1 2020-12-30 CRAN (R 4.1.0)
crayon 1.5.0 2022-02-14 CRAN (R 4.1.2)
curl 4.3.2 2021-06-23 CRAN (R 4.1.0)
data.table 1.14.2 2021-09-27 CRAN (R 4.1.1)
DBI 1.1.1 2021-01-15 CRAN (R 4.1.0)
dbplyr 2.1.1 2021-04-06 CRAN (R 4.1.0)
desc 1.4.0 2021-09-28 CRAN (R 4.1.1)
devtools 2.4.2 2021-06-07 CRAN (R 4.1.0)
digest 0.6.27 2020-10-24 CRAN (R 4.1.0)
dplyr 1.0.7 2021-06-18 CRAN (R 4.1.0)
ellipsis 0.3.2 2021-04-29 CRAN (R 4.1.0)
evaluate 0.14 2019-05-28 CRAN (R 4.1.0)
fansi 0.5.0 2021-05-25 CRAN (R 4.1.0)
fastmap 1.1.0 2021-01-25 CRAN (R 4.1.0)
fdrtool 1.2.16 2021-01-06 CRAN (R 4.1.0)
flextable 0.7.0 2022-03-06 CRAN (R 4.1.2)
forcats 0.5.1 2021-01-27 CRAN (R 4.1.0)
foreign 0.8-81 2020-12-22 CRAN (R 4.1.2)
Formula 1.2-4 2020-10-16 CRAN (R 4.1.0)
fs 1.5.0 2020-07-31 CRAN (R 4.1.0)
gdtools 0.2.3 2021-01-06 CRAN (R 4.1.1)
generics 0.1.2 2022-01-31 CRAN (R 4.1.2)
ggplot2 3.3.5 2021-06-25 CRAN (R 4.1.2)
glasso 1.11 2019-10-01 CRAN (R 4.1.0)
glue 1.4.2 2020-08-27 CRAN (R 4.1.0)
gridExtra 2.3 2017-09-09 CRAN (R 4.1.0)
gtable 0.3.0 2019-03-25 CRAN (R 4.1.0)
gtools 3.9.2 2021-06-06 CRAN (R 4.1.0)
haven 2.4.3 2021-08-04 CRAN (R 4.1.1)
Hmisc 4.5-0 2021-02-28 CRAN (R 4.1.1)
hms 1.1.1 2021-09-26 CRAN (R 4.1.1)
htmlTable 2.2.1 2021-05-18 CRAN (R 4.1.0)
htmltools 0.5.2 2021-08-25 CRAN (R 4.1.1)
htmlwidgets 1.5.4 2021-09-08 CRAN (R 4.1.1)
httr 1.4.2 2020-07-20 CRAN (R 4.1.0)
igraph 1.2.6 2020-10-06 CRAN (R 4.1.0)
jpeg 0.1-9 2021-07-24 CRAN (R 4.1.0)
jquerylib 0.1.4 2021-04-26 CRAN (R 4.1.0)
jsonlite 1.7.2 2020-12-09 CRAN (R 4.1.0)
kableExtra 1.3.4 2021-02-20 CRAN (R 4.1.0)
knitr 1.36 2021-09-29 CRAN (R 4.1.1)
kutils 1.70 2020-04-29 CRAN (R 4.1.0)
lattice 0.20-45 2021-09-22 CRAN (R 4.1.2)
latticeExtra 0.6-29 2019-12-19 CRAN (R 4.1.0)
lavaan 0.6-10 2022-01-25 CRAN (R 4.1.2)
lifecycle 1.0.1 2021-09-24 CRAN (R 4.1.1)
lisrelToR 0.1.4 2013-05-08 CRAN (R 4.1.0)
lme4 1.1-27.1 2021-06-22 CRAN (R 4.1.0)
lmerTest 3.1-3 2020-10-23 CRAN (R 4.1.2)
lubridate 1.7.10 2021-02-26 CRAN (R 4.1.0)
magick 2.7.3 2021-08-18 CRAN (R 4.1.1)
magrittr 2.0.1 2020-11-17 CRAN (R 4.1.0)
MASS 7.3-54 2021-05-03 CRAN (R 4.1.2)
Matrix 1.3-4 2021-06-01 CRAN (R 4.1.2)
matrixcalc 1.0-5 2021-07-28 CRAN (R 4.1.0)
matrixStats 0.61.0 2021-09-17 CRAN (R 4.1.1)
memoise 2.0.0 2021-01-26 CRAN (R 4.1.0)
mi 1.0 2015-04-16 CRAN (R 4.1.0)
mice 3.13.0 2021-01-27 CRAN (R 4.1.0)
minqa 1.2.4 2014-10-09 CRAN (R 4.1.0)
mnormt 2.0.2 2020-09-01 CRAN (R 4.1.0)
modelr 0.1.8 2020-05-19 CRAN (R 4.1.0)
munsell 0.5.0 2018-06-12 CRAN (R 4.1.0)
nlme 3.1-153 2021-09-07 CRAN (R 4.1.2)
nloptr 1.2.2.2 2020-07-02 CRAN (R 4.1.0)
nnet 7.3-16 2021-05-03 CRAN (R 4.1.2)
numDeriv 2016.8-1.1 2019-06-06 CRAN (R 4.1.0)
officer 0.4.1 2021-11-14 CRAN (R 4.1.2)
OpenMx 2.19.8 2021-09-06 CRAN (R 4.1.1)
openxlsx 4.2.4 2021-06-16 CRAN (R 4.1.0)
pander 0.6.4 2021-06-13 CRAN (R 4.1.0)
pbapply 1.5-0 2021-09-16 CRAN (R 4.1.1)
pbivnorm 0.6.0 2015-01-23 CRAN (R 4.1.0)
pillar 1.7.0 2022-02-01 CRAN (R 4.1.2)
pkgbuild 1.2.0 2020-12-15 CRAN (R 4.1.0)
pkgconfig 2.0.3 2019-09-22 CRAN (R 4.1.0)
pkgload 1.2.2 2021-09-11 CRAN (R 4.1.1)
plyr 1.8.6 2020-03-03 CRAN (R 4.1.0)
png 0.1-7 2013-12-03 CRAN (R 4.1.0)
prettyunits 1.1.1 2020-01-24 CRAN (R 4.1.0)
processx 3.5.2 2021-04-30 CRAN (R 4.1.0)
pryr 0.1.5 2021-07-26 CRAN (R 4.1.1)
ps 1.6.0 2021-02-28 CRAN (R 4.1.0)
psych 2.1.9 2021-09-22 CRAN (R 4.1.1)
purrr 0.3.4 2020-04-17 CRAN (R 4.1.0)
qgraph 1.6.9 2021-01-28 CRAN (R 4.1.0)
R6 2.5.1 2021-08-19 CRAN (R 4.1.1)
rapportools 1.0 2014-01-07 CRAN (R 4.1.0)
RColorBrewer 1.1-2 2014-12-07 CRAN (R 4.1.0)
Rcpp 1.0.7 2021-07-07 CRAN (R 4.1.1)
RcppParallel 5.1.4 2021-05-04 CRAN (R 4.1.0)
readr 2.0.2 2021-09-27 CRAN (R 4.1.1)
readxl 1.3.1 2019-03-13 CRAN (R 4.1.0)
regsem 1.8.0 2021-06-03 CRAN (R 4.1.0)
remotes 2.4.1 2021-09-29 CRAN (R 4.1.1)
reprex 2.0.1 2021-08-05 CRAN (R 4.1.1)
reshape2 1.4.4 2020-04-09 CRAN (R 4.1.0)
rio 0.5.27 2021-06-21 CRAN (R 4.1.0)
rlang 1.0.2 2022-03-04 CRAN (R 4.1.2)
rmarkdown 2.11 2021-09-14 CRAN (R 4.1.1)
rockchalk 1.8.144 2019-03-08 CRAN (R 4.1.0)
rpart 4.1-15 2019-04-12 CRAN (R 4.1.2)
rprojroot 2.0.2 2020-11-15 CRAN (R 4.1.0)
Rsolnp 1.16 2015-12-28 CRAN (R 4.1.0)
rstudioapi 0.13 2020-11-12 CRAN (R 4.1.0)
rvest 1.0.1 2021-07-26 CRAN (R 4.1.0)
sass 0.4.0 2021-05-12 CRAN (R 4.1.0)
scales 1.1.1 2020-05-11 CRAN (R 4.1.0)
sem 3.1-11 2020-05-19 CRAN (R 4.1.0)
semPlot 1.1.2 2019-08-20 CRAN (R 4.1.0)
sessioninfo 1.1.1 2018-11-05 CRAN (R 4.1.0)
stringi 1.7.4 2021-08-25 CRAN (R 4.1.1)
stringr 1.4.0 2019-02-10 CRAN (R 4.1.0)
summarytools 1.0.0 2021-07-28 CRAN (R 4.1.1)
survival 3.2-13 2021-08-24 CRAN (R 4.1.2)
svglite 2.0.0 2021-02-20 CRAN (R 4.1.0)
systemfonts 1.0.2 2021-05-11 CRAN (R 4.1.0)
testthat 3.0.4 2021-07-01 CRAN (R 4.1.0)
tibble 3.1.4 2021-08-25 CRAN (R 4.1.1)
tidyr 1.1.4 2021-09-27 CRAN (R 4.1.1)
tidyselect 1.1.2 2022-02-21 CRAN (R 4.1.2)
tidyverse 1.3.1 2021-04-15 CRAN (R 4.1.0)
tmvnsim 1.0-2 2016-12-15 CRAN (R 4.1.0)
truncnorm 1.0-8 2018-02-27 CRAN (R 4.1.0)
tzdb 0.1.2 2021-07-20 CRAN (R 4.1.1)
usethis 2.0.1 2021-02-10 CRAN (R 4.1.0)
utf8 1.2.2 2021-07-24 CRAN (R 4.1.1)
uuid 0.1-4 2020-02-26 CRAN (R 4.1.0)
vctrs 0.3.8 2021-04-29 CRAN (R 4.1.0)
viridisLite 0.4.0 2021-04-13 CRAN (R 4.1.0)
webshot 0.5.2 2019-11-22 CRAN (R 4.1.0)
withr 2.5.0 2022-03-03 CRAN (R 4.1.2)
xfun 0.26 2021-09-14 CRAN (R 4.1.1)
XML 3.99-0.8 2021-09-17 CRAN (R 4.1.1)
xml2 1.3.2 2020-04-23 CRAN (R 4.1.0)
xtable 1.8-4 2019-04-21 CRAN (R 4.1.0)
yaml 2.2.1 2020-02-01 CRAN (R 4.1.0)
zip 2.2.0 2021-05-31 CRAN (R 4.1.0) 
# load in data
demo  <- rio::import(file = "Demographics.sav")
full  <- rio::import(file = "Full data set 196 participants.sav")
rrs10 <- rio::import(file = "RRS-10 196 participants.sav")

# center mindfulness
full <- full %>%
  mutate(MINDc = scale(MIND, scale = FALSE),
         MIND_T1c = scale(MIND_T1, scale = FALSE),
         MIND_T2c = scale(MIND_T2, scale = FALSE),
         MIND_T3c = scale(MIND_T3, scale = FALSE),
         MIND_T4c = scale(MIND_T4, scale = FALSE))

# descriptives
full %>% select(DAS_D:RRS) %>%descr(stats = "common", order = "p")

Descriptive Statistics
full
N: 196

DAS_D MIND COVID RRS
Mean 1.80 2.56 3.11 2.08
Std.Dev 0.60 0.47 0.72 0.59
Min 1.00 1.21 1.23 1.06
Median 1.71 2.55 3.10 2.06
Max 3.86 3.75 5.00 3.91
N.Valid 196.00 196.00 196.00 196.00
Pct.Valid 100.00 100.00 100.00 100.00 
cormat <- round(x = cor(select(full, c("COVID", "RRS", "MIND", "DAS_D"))),
                digits = 2)
cormat[lower.tri(cormat)] <- NA
cormat <- melt(cormat, na.rm = TRUE) %>%
  # delete the perfect correlations
  filter(value != 1.00)
ggplot(data = cormat, aes(Var2, Var1, fill = value)) +
 geom_tile(color = "white" ) +
 scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                      midpoint = 0, limit = c(-1,1), space = "Lab", 
                      name = "Pearson\nCorrelation") +
  theme_minimal()+ 
  labs(x = NULL, y = NULL) +
  coord_fixed() +
  geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.major = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    axis.ticks = element_blank(),
    legend.justification = c(1, 0),
    legend.position = c(0.4, 0.7),
    legend.direction = "horizontal")+
  guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
                               title.position = "top", title.hjust = 0.5))

I. Missingness analyses

# create data frame that has missingness indicator and ID
missing <- data.frame(Respondent_ID = full$Respondent_ID)
missing$missing <- apply(X = select(full, starts_with(c("COVID_", "DAS_D_",
                                                        "RRS_", "MIND_"))),
                         MARGIN = 1,
                         FUN = function(x) as.numeric(any(is.na(x))))
 
# merge in demographic info
missing <- merge(missing, demo)

# merge in average scores
scale <- c("Respondent_ID", "DAS_D", "COVID", "MIND", "RRS")
missing <- merge(missing, 
                 select(full, all_of(scale)))     

# pull out missing ID
missing_ID <- missing %>% 
  filter(missing == 1) %>% 
  pull(Respondent_ID)

# recode demo variables in missing data frame
missing <- missing %>%
  mutate(gender_f = ifelse(Gender == 1, 0, # 0 = female
                    ifelse(Gender == 2, 1, # 1 = male
                    NA)),                  
         race_f   = as.numeric(Race != 5), # 0 = white, 1 = non-white
         employ_f = ifelse(Employment_T1 %in% c(3, 4, 5, 6, 7), 0, # 0 = no work
                    ifelse(Employment_T1 %in% c(1, 2), 1,          # 1 = work
                           NA)))

# general missing patterns for the main 4 variables
full %>% select(COVID_T1:COVID_T4) %>% md.pattern()

##     COVID_T1 COVID_T2 COVID_T3 COVID_T4    
## 150        1        1        1        1   0
## 12         1        1        1        0   1
## 9          1        1        0        0   2
## 25         1        0        0        0   3
##            0       25       34       46 105
full %>% select(DAS_D_T1:DAS_D_T4) %>% md.pattern()

##     DAS_D_T1 DAS_D_T2 DAS_D_T3 DAS_D_T4    
## 149        1        1        1        1   0
## 12         1        1        1        0   1
## 1          1        1        0        1   1
## 10         1        1        0        0   2
## 24         1        0        0        0   3
##            0       24       35       46 105
full %>% select(RRS_T1:RRS_T4) %>% md.pattern()

##     RRS_T1 RRS_T2 RRS_T3 RRS_T4    
## 149      1      1      1      1   0
## 12       1      1      1      0   1
## 1        1      1      0      1   1
## 10       1      1      0      0   2
## 24       1      0      0      0   3
##          0     24     35     46 105
full %>% select(MIND_T1:MIND_T4) %>% md.pattern()

##     MIND_T1 MIND_T2 MIND_T3 MIND_T4    
## 150       1       1       1       1   0
## 12        1       1       1       0   1
## 10        1       1       0       0   2
## 24        1       0       0       0   3
##           0      24      34      46 104
# compare age
summary(
  glm(missing ~ Age, family = "binomial", data = missing))$coefficients %>%
  as.data.frame() %>%
  knitr::kable(
    caption = "Age (continuous) predicting missingness") %>%
  kable_styling()
Age (continuous) predicting missingness
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.2277 1.059 -0.2150 0.8297
Age -0.0440 0.050 -0.8793 0.3792 
# compare gender
summary(
  glm(missing ~ gender_f, family = "binomial", data = missing))$coefficients %>%
  as.data.frame() %>%
  knitr::kable(
    caption = "Gender (0 = female, 1 = male) predicting missingness") %>%
  kable_styling()
Gender (0 = female, 1 = male) predicting missingness
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.2040 0.1900 -6.3357 0.000
gender_f 0.2877 0.4197 0.6855 0.493 
# compare race
summary(
  glm(missing ~ race_f, family = "binomial", data = missing))$coefficients %>%
  as.data.frame() %>%
  knitr::kable(
    caption = "Race (0 = white, 1 = non-white) predicting missingness") %>%
  kable_styling()
Race (0 = white, 1 = non-white) predicting missingness
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.1575 0.2444 -4.7355 0.0000
race_f 0.0069 0.3353 0.0205 0.9836 
# compare employment
summary(
  glm(missing ~ employ_f, family = "binomial", data = missing))$coefficients %>%
  as.data.frame() %>%
  knitr::kable(
    caption = "Employment at wave 1 (0 = not currently working, 1 = working) 
    predicting missingness") %>%
  kable_styling()
Employment at wave 1 (0 = not currently working, 1 = working) predicting missingness
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3705 0.2800 -4.896 0.0000
employ_f 0.3737 0.3501 1.068 0.2857 
# compare mindfulness
summary(
  glm(missing ~ MIND, family = "binomial", data = missing))$coefficients %>%
  as.data.frame() %>%
  knitr::kable(
    caption = "Mindfulness average scores (FMI; Walach et al., 2006) 
    predicting missingness") %>%
  kable_styling()
Mindfulness average scores (FMI; Walach et al., 2006) predicting missingness
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.8629 0.9172 -0.9407 0.3468
MIND -0.1141 0.3545 -0.3217 0.7476 
# compare covid stress
summary(
  glm(missing ~ COVID, family = "binomial", data = missing))$coefficients %>%
  as.data.frame() %>%
  knitr::kable(
    caption = "COVID stress average scores predicting missingness") %>%
  kable_styling()
COVID stress average scores predicting missingness
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.1992 0.7827 -2.810 0.0050
COVID 0.3317 0.2393 1.386 0.1657 
# compare depression
summary(
  glm(missing ~ DAS_D, family = "binomial", data = missing))$coefficients %>%
  as.data.frame() %>%
  knitr::kable(
    caption = "Depression average scores (DASS; Lovibond & Lovibond, 1995) 
    predicting missingness") %>%
  kable_styling()
Depression average scores (DASS; Lovibond & Lovibond, 1995) predicting missingness
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.7110 0.5318 -3.217 0.0013
DAS_D 0.3045 0.2717 1.121 0.2625 
# compare rumination
summary(
  glm(missing ~ RRS, family = "binomial", data = missing))$coefficients %>%
  as.data.frame() %>%
  knitr::kable(
    caption = "Rumination average scores (RRS; Nolen-Hoeksema & Morrow, 1991) 
    predicting missingness") %>%
  kable_styling()
Rumination average scores (RRS; Nolen-Hoeksema & Morrow, 1991) predicting missingness
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.304 0.6411 -3.594 0.0003
RRS 0.542 0.2854 1.899 0.0576

There are 47 participants with at least 1 missing data point and 149 with complete data. There was no significant differences between participants with full vs. missing data, either in demographic variables or measured scale responses.

II. Robustness analyses

# create functions to create and export tables for mediation/mod-mediation
med_table <- function(a, b, c, indirect,
                      labels = c("a", "b", "c (direct)", "indirect"),
                      footer = NULL) {
  df <- rbind(a,b,c,indirect) %>% as.data.frame()
  names(df) <- c("B", "t", "p", "LLCI", "ULCI")
  df <- cbind(effects = labels,
              df)
  ft <- flextable(df)
  ft <- theme_vanilla(ft)
  if (!is.null(footer)) {
    ft <- add_footer_lines(ft, footer)
  }
  ft
}

modmed_table <- function(a1, a2, a3, cp, b1,
                         R2a, R2b, Fpa, Fpb,
                         labels = c("X", "W", "XW", "M",
                                    "R-squared", "F p-value"),
                         footer = NULL) {
  header <- c("Predictors", "B", "t", "p", " ", "B ", "t ", "p ", "  ")
  X <- c(a1, "", cp, "")
  W <- c(a2, rep(NA, 5))
  XW <- c(a3, rep(NA, 5))
  M <- c(rep(NA, 4), b1, "")
  R2 <- c(rep(NA, 3), R2a, rep(NA, 3), R2b)
  pvalue <- c(rep(NA, 3), Fpa, rep(NA, 3), Fpb)
  
  df <- rbind(X, W, XW, M, R2, pvalue) %>% as.data.frame()
  df <- cbind(Predictors = labels,
              df)
  names(df) <- header
  ft <- flextable(df)
  if (!is.null(footer)) {
    ft <- add_footer_lines(ft, footer)
  }
  ft <- theme_vanilla(ft)
  ft <- vline(ft, j = c("Predictors", " "))
  ft <- hline(ft, i = 4)
  ft <- add_header_row(ft,
                       colwidths = c(1, 4, 4),
                       values = c("", "Model 1 (Rumination)", 
                                  "Model 2 (Depression)"))
  ft
}

modmed_condl <- function(low, mean, high) {
  header <- c("MIND", "Effect", "BootSE", "LLCI", "ULCI")
  df <- rbind(low, mean, high) %>% as.data.frame()
  names(df) <- header
  ft <- flextable(df)
  ft <- theme_vanilla(ft)
  ft
}

II.1. Using data from different waves

Replicating SPSS analyses in R

Mediation - Average

process(data = full,
        y = "DAS_D",
        x = "COVID",
        m = "RRS",
        model = 4,
        boot = 5000,
        seed = 20220308)
a <- c(0.3182, 5.8810, "<0.001", 0.2115, 0.4249)
b <- c(0.7213, 13.3337, "<0.001", 0.6146, 0.8280)
c <- c(0.0615, 1.3901, 0.1661, -0.0258, 0.1488)
indirect <- c(0.2295, NA, NA, 0.1506, 0.3132)

(med_avg <- med_table(a = a, b = b, c = c, indirect = indirect))

Moderated mediation - Average

process(data = full,
        y = "DAS_D",
        x = "COVID",
        m = "RRS",
        w = "MIND",
        model = 7,
        moments = 1,
        center = 2, # continuous variables in interaction terms
        boot = 5000,
        seed = 20220309)
a1 <- c(0.3021, 6.3427, "<0.001")
a2 <- c(-0.5682,  -7.8077, "<0.001")
a3 <- c(-0.2197, -2.5726, 0.0108)
cp <- c(0.0615, 1.3901, 0.1661)
b1 <- c(0.7213, 13.3337, "<0.001")
R2a <- 0.3592
R2b <- 0.5424
Fpa <- "<0.001"
Fpb <- "<0.001"
low <- c(-0.4742,    0.2931,    0.0497,    0.1991,    0.3936)
mean <- c(0.0000,    0.2179,    0.0394,    0.1416,    0.2962)
high <- c(0.4742,    0.1428,    0.0438,    0.0578,    0.2322)

(modmed_avg <- modmed_table(a1 = a1, a2 = a2, a3 = a3, cp = cp, b1 = b1,
                           R2a = R2a, R2b = R2b, Fpa = Fpa, Fpb = Fpb)) 
(modmed_avg_condl <- modmed_condl(low = low, mean = mean, high = high))

Wave 1 data only

Mediation - Wave 1

process(data = full,
        y = "DAS_D_T1",
        x = "COVID_T1",
        m = "RRS_T1",
        model = 4,
        boot = 5000,
        seed = 20220308)
a <- c(0.3388, 5.8085, "<0.001",0.2238, 0.4539)
b <- c(0.6159, 9.9135, "<0.001", 0.4934, 0.7384)
c <- c(0.1062, 1.9416, 0.0536, -0.0017, 0.2141)
indirect <- c(0.2087, NA, NA, 0.1377, 0.2837)

(med_w1 <- med_table(a = a, b = b, c = c, indirect = indirect))

Moderated mediation - Wave 1

process(data = full,
        y = "DAS_D_T1",
        x = "COVID_T1",
        m = "RRS_T1",
        w = "MIND_T1",
        model = 7,
        moments = 1,
        center = 2, # continuous variables in interaction terms
        boot = 5000,
        seed = 20220309)
a1 <- c(0.3373, 6.1311, "<0.001")
a2 <- c(-0.3887, -5.0186, "<0.001")
a3 <- c(-0.1475, -1.5675, 0.1186)
cp <- c(0.1062, 1.9416, 0.0536)
b1 <- c(0.6159, 9.9135, "<0.001")
R2a <- 0.2507
R2b <- 0.4155
Fpa <- "<0.001"
Fpb <- "<0.001"
low <- c(-0.5137,    0.2544,    0.0466,    0.1670,    0.3511)
mean <- c(0.0000,    0.2077,    0.0372,    0.1390,    0.2852)
high <- c(0.5137,    0.1611,    0.0404,    0.0842,    0.2447)

(modmed_w1 <- modmed_table(a1 = a1, a2 = a2, a3 = a3, cp = cp, b1 = b1,
                           R2a = R2a, R2b = R2b, Fpa = Fpa, Fpb = Fpb))
(modmed_w1_condl <- modmed_condl(low = low, mean = mean, high = high))

Wave 2 data only

Mediation - Wave 2

process(data = full,
        y = "DAS_D_T2",
        x = "COVID_T2",
        m = "RRS_T2",
        model = 4,
        boot = 5000,
        seed = 20220308)
a <- c(0.3358, 5.8736, "<0.001",0.2229, 0.4486)
b <- c(0.6376, 10.2228, "<0.001", 0.5145, 0.7608)
c <- c(0.1143, 2.2472, 0.0259, 0.0139, 0.2147)
indirect <- c(0.2141, NA, NA, 0.1243, 0.3212)

(med_w2 <- med_table(a = a, b = b, c = c, indirect = indirect))

Moderated mediation - Wave 2

process(data = full,
        y = "DAS_D_T2",
        x = "COVID_T2",
        m = "RRS_T2",
        w = "MIND_T2",
        model = 7,
        moments = 1,
        center = 2, # continuous variables in interaction terms
        boot = 5000,
        seed = 20220309)
a1 <- c(0.3223, 6.0988, "<0.001")
a2 <- c(-0.4312, -5.5346, "<0.001")
a3 <- c(-0.1896, -2.1841, 0.0303)
cp <- c(0.1143, 2.2472, 0.0259)
b1 <- c( 0.6376, 10.2228, "<0.001")
R2a <- 0.3035
R2b <- 0.4794
Fpa <- "<0.001"
Fpb <- "<0.001"
low <- c(-0.5389,    0.2707,    0.0605,    0.1565,    0.3918)
mean <- c(0.0000,    0.2055,    0.0478,    0.1175,    0.3062)
high <- c(0.5389,    0.1404,    0.0634,    0.0184,    0.2692)
  
(modmed_w2 <- modmed_table(a1 = a1, a2 = a2, a3 = a3, cp = cp, b1 = b1,
                           R2a = R2a, R2b = R2b, Fpa = Fpa, Fpb = Fpb))
(modmed_w2_condl <- modmed_condl(low = low, mean = mean, high = high))

Wave 3 data only

Mediation - Wave 3

process(data = full,
        y = "DAS_D_T3",
        x = "COVID_T3",
        m = "RRS_T3",
        model = 4,
        boot = 5000,
        seed = 20220308)
a <- c(0.2257, 3.6596, "<0.001",0.1039 , 0.3475)
b <- c(0.7379, 10.8904, "<0.001", 0.6040, 0.8717)
c <- c(0.0159, 0.2890, 0.7730, -0.0925, 0.1242)
indirect <- c(0.1665, NA, NA, 0.0803, 0.2625)

(med_w3 <- med_table(a = a, b = b, c = c, indirect = indirect))

Moderated mediation - Wave 3

process(data = full,
        y = "DAS_D_T3",
        x = "COVID_T3",
        m = "RRS_T3",
        w = "MIND_T3",
        model = 7,
        moments = 1,
        center = 2, # continuous variables in interaction terms
        boot = 5000,
        seed = 20220309)
a1 <- c(0.2043, 3.4614, "<0.001")
a2 <- c(-0.3514, -4.2507, "<0.001")
a3 <- c(-0.0709, -0.7906, 0.4304)
cp <- c(0.0159, 0.2890, 0.7730)
b1 <- c(0.7379, 10.8904, "<0.001")
R2a <- 0.1729
R2b <- 0.4525
Fpa <- "<0.001"
Fpb <- "<0.001"
low <- c(-0.5629,    0.1802,    0.0664,    0.0442,    0.3072)
mean <- c(0.0000,    0.1508,    0.0472,    0.0592,    0.2465)
high <- c(0.5629,    0.1213,    0.0529,    0.0181,    0.2287)

(modmed_w3 <- modmed_table(a1 = a1, a2 = a2, a3 = a3, cp = cp, b1 = b1,
                           R2a = R2a, R2b = R2b, Fpa = Fpa, Fpb = Fpb))
(modmed_w3_condl <- modmed_condl(low = low, mean = mean, high = high))

Wave 4 data only

Mediation - Wave 4

process(data = full,
        y = "DAS_D_T4",
        x = "COVID_T4",
        m = "RRS_T4",
        model = 4,
        boot = 5000,
        seed = 20220308)
a <- c(0.1777, 2.8369, 0.0052, 0.0539 , 0.3015)
b <- c(0.7617, 12.2976, "<0.001", 0.6393, 0.8841)
c <- c(0.0677, 1.3974, 0.1644, -0.0281, 0.1635)
indirect <- c(0.1354, NA, NA, 0.0302, 0.2482)

(med_w4 <- med_table(a = a, b = b, c = c, indirect = indirect))

Moderated mediation - Wave 4

process(data = full,
        y = "DAS_D_T4",
        x = "COVID_T4",
        m = "RRS_T4",
        w = "MIND_T4",
        model = 7,
        moments = 1,
        center = 2, # continuous variables in interaction terms
        boot = 5000,
        seed = 20220309)
a1 <- c(0.1845, 3.2802, 0.0013)
a2 <- c(-0.4931, -5.8144, "<0.001")
a3 <- c(-0.2766, -3.0717, 0.0025)
cp <- c(0.0677, 1.3974, 0.1644)
b1 <- c(0.7617, 12.2976, "<0.001")
R2a <- 0.2560
R2b <- 0.5359
Fpa <- "<0.001"
Fpb <- "<0.001"
low <- c(-0.5520,    0.2568,    0.0592,    0.1421,    0.3788)
mean <- c(0.0000,    0.1405,    0.0505,    0.0455,    0.2426)
high <- c(0.5520,    0.0242,    0.0665,   -0.1168,    0.1471)

(modmed_w4 <- modmed_table(a1 = a1, a2 = a2, a3 = a3, cp = cp, b1 = b1,
                           R2a = R2a, R2b = R2b, Fpa = Fpa, Fpb = Fpb))
(modmed_w4_condl <- modmed_condl(low = low, mean = mean, high = high))

II.2. Excluding depression subscale in Rumination (RRS-10)

Mediation - RRS-10

process(data = rrs10,
        y = "DAS_D",
        x = "COVID",
        m = "RRS",
        model = 4,
        boot = 5000,
        seed = 20220308)
a <- c(0.2998, 5.5227, "<0.001", 0.1927 , 0.4069)
b <- c(0.5915, 9.6296, "<0.001", 0.4704, 0.7127)
c <- c(0.1137, 2.2761, 0.0239, 0.0152, 0.2123)
indirect <- c(0.1773, NA, NA, 0.1060, 0.2581)

(med_rrs10 <- med_table(a = a, b = b, c = c, indirect = indirect))

Moderated mediation - RRS-10

process(data = rrs10,
        y = "DAS_D",
        x = "COVID",
        m = "RRS",
        w = "MIND",
        model = 7,
        moments = 1,
        center = 2, # continuous variables in interaction terms
        boot = 5000,
        seed = 20220309)
a1 <- c(0.2876, 5.6432, "<0.001")
a2 <- c( -0.4329, -5.5601, "<0.001")
a3 <- c(-0.1682, -1.8407,  0.0672)
cp <- c(0.1137, 2.2761, 0.0239)
b1 <- c(0.5915, 9.6296, "<0.001")
R2a <- 0.2580
R2b <- 0.4062
Fpa <- "<0.001"
Fpb <- "<0.001"
low <- c(-0.4742,    0.2173,    0.0494,    0.1245,    0.3188)
mean <- c(0.0000,    0.1701,    0.0380,    0.1013,    0.2491)
high <- c(0.4742,    0.1229,    0.0406,    0.0470,    0.2053)

(modmed_rrs10 <- modmed_table(a1 = a1, a2 = a2, a3 = a3, cp = cp, b1 = b1,
                              R2a = R2a, R2b = R2b, Fpa = Fpa, Fpb = Fpb))
(modmed_rrs10_condl <- modmed_condl(low = low, mean = mean, high = high))

III. Change analyses

# select only needed data
DAS_D <- full %>%
  select(Respondent_ID, DAS_D_T1:DAS_D_T4)
MIND <- full %>%
  select(Respondent_ID, MIND_T1:MIND_T4)
RRS <- full %>%
  select(Respondent_ID, RRS_T1:RRS_T4)
COVID <- full %>%
  select(Respondent_ID, COVID_T1:COVID_T4)

# melt wide to long
DAS_D <- melt(DAS_D, id.vars = c("Respondent_ID")) %>%
  mutate(time = ifelse(grepl(pattern = "T1", x = variable), 0,
                ifelse(grepl(pattern = "T2", x = variable), 1,
                ifelse(grepl(pattern = "T3", x = variable), 2, 3)))) %>%
  select(-variable) %>%
  rename(DAS_D = value)
MIND <- melt(MIND, id.vars = c("Respondent_ID")) %>%
  mutate(time = ifelse(grepl(pattern = "T1", x = variable), 0,
                ifelse(grepl(pattern = "T2", x = variable), 1,
                ifelse(grepl(pattern = "T3", x = variable), 2, 3)))) %>%
  select(-variable) %>%
  rename(MIND = value)
RRS <- melt(RRS, id.vars = c("Respondent_ID")) %>%
  mutate(time = ifelse(grepl(pattern = "T1", x = variable), 0,
                ifelse(grepl(pattern = "T2", x = variable), 1,
                ifelse(grepl(pattern = "T3", x = variable), 2, 3)))) %>%
  select(-variable) %>%
  rename(RRS = value)
COVID <- melt(COVID, id.vars = c("Respondent_ID")) %>%
  mutate(time = ifelse(grepl(pattern = "T1", x = variable), 0,
                ifelse(grepl(pattern = "T2", x = variable), 1,
                ifelse(grepl(pattern = "T3", x = variable), 2, 3)))) %>%
  select(-variable) %>%
  rename(COVID = value)

# merge all variables to a long data frame
long <- merge(DAS_D, MIND) %>%
  merge(RRS) %>%
  merge(COVID)

Spaghetti plots

# random sample for plots
set.seed(20200308)
sample <- sample(full$Respondent_ID, size = 75)
sample <- long %>% filter(Respondent_ID %in% c(sample))

# COVID stress
p1 <- ggplot(data = sample,
             aes(x = time, y = COVID, group = Respondent_ID)) + 
  geom_line(linetype = "dashed") +
  theme_bw() +
  stat_smooth(aes(data = sample$COVID, group = 1),
              method = "lm",
              formula = y ~ poly(x, 2),
              lwd = 1.5, color = "red") +
  labs(x = "Time in weeks",
       y = "COVID stress")

# Depression DAS_D
p2 <- ggplot(data = sample,
             aes(x = time, y = DAS_D, group = Respondent_ID)) + 
  geom_line(linetype = "dashed") +
  theme_bw() +
  stat_smooth(aes(data = sample$DAS_D, group = 1),
              method = "lm",
              formula = y ~ poly(x, 2),
              lwd = 1.5, color = "red") +
  labs(x = "Time in weeks",
       y = "Depression (DAS_D)")

# Rumination
p3 <- ggplot(data = sample,
             aes(x = time, y = RRS, group = Respondent_ID)) + 
  geom_line(linetype = "dashed") +
  theme_bw() +
  stat_smooth(aes(data = sample$RRS, group = 1),
              method = "lm",
              formula = y ~ poly(x, 2),
              lwd = 1.5, color = "red") +
  labs(x = "Time in weeks",
       y = "Rumination (RRS)")

# Mindfulness
p4 <- ggplot(data = sample,
             aes(x = time, y = MIND, group = Respondent_ID)) + 
  geom_line(linetype = "dashed") +
  theme_bw() +
  stat_smooth(aes(data = sample$MIND, group = 1),
              method = "lm",
              formula = y ~ poly(x, 2),
              lwd = 1.5, color = "red") +
  labs(x = "Time in weeks",
       y = "Mindfulness (FMI)")

# all 4 plots in 4 panels
plot_grid(p1, p2, p3, p4, 
          nrow = 2)

Mean-level change

# covid stress
linear.covid <- lmer(COVID ~ time + (1 + time | Respondent_ID),
                     data = long)
summary(linear.covid)$coefficients %>%
  as.data.frame() %>%
  knitr::kable(caption = "COVID stress: 
               linear change with random intercept and random slope") %>%
  kable_styling()
COVID stress: linear change with random intercept and random slope
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.285 0.0506 194.9 64.928 0
time -0.146 0.0159 149.5 -9.205 0 
# depression
linear.dasd <- lmer(DAS_D ~ time + (1 + time | Respondent_ID),
                    data = long)
summary(linear.dasd)$coefficients %>%
  as.data.frame() %>%
  knitr::kable(caption = "Depression: 
               linear change with random intercept and random slope") %>%
  kable_styling()
Depression: linear change with random intercept and random slope
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.8342 0.0463 193.6 39.63 0.0000
time -0.0283 0.0166 164.2 -1.71 0.0892 
# rumination
linear.rrs <- lmer(RRS ~ time + (1 + time | Respondent_ID),
                   data = long)
summary(linear.rrs)$coefficients %>%
  as.data.frame() %>%
  knitr::kable(caption = "Rumination: 
               linear change with random intercept and random slope") %>%
  kable_styling()
Rumination: linear change with random intercept and random slope
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.1275 0.0446 194.4 47.697 0e+00
time -0.0454 0.0132 161.9 -3.429 8e-04 
# mindfulness
linear.mind <- lmer(MIND ~ time + (1 + time | Respondent_ID),
                    data = long)
summary(linear.mind)$coefficients %>%
  as.data.frame() %>%
  knitr::kable(caption = "Mindfulness: 
               linear change with random intercept and random slope") %>%
  kable_styling()
Mindfulness: linear change with random intercept and random slope
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.5747 0.0354 195.5 72.716 0.0000
time -0.0123 0.0112 160.5 -1.103 0.2716

We fit linear mixed-effects models with random intercepts and random slopes to test for linear change in the four main variables over time. There was a significant change in COVID stress and Rumination: COVID stress decreased by 0.146 points per week (\(t_{(149.5)} = -9.205, p < .001\)) and Rumination decreased by 0.045 points per week (\(t_{(161.9)} = -3.429, p = .0008\)). There was no significant change in depression symptoms or mindfulness at an alpha level of 0.05.

Individual differences in change

# covid stress
linear.covid.int <- lmer(COVID ~ time + (1 | Respondent_ID),
                         data = long)
anova(linear.covid.int, linear.covid) %>%
  as.data.frame() %>%
  knitr::kable(caption = "COVID stress: 
               Model comparison for individual differences in change") %>%
  kable_styling()
COVID stress: Model comparison for individual differences in change
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
linear.covid.int 4 1085 1103 -538.5 1077 NA NA NA
linear.covid 6 1067 1094 -527.4 1055 22.15 2 0 
# depression
linear.dasd.int <- lmer(DAS_D ~ time + (1 | Respondent_ID),
                        data = long)
anova(linear.dasd.int, linear.dasd) %>%
  as.data.frame() %>%
  knitr::kable(caption = "Depression: 
               Model comparison for individual differences in change") %>%
  kable_styling()
Depression: Model comparison for individual differences in change
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
linear.dasd.int 4 1110 1128 -551.1 1102 NA NA NA
linear.dasd 6 1107 1134 -547.4 1095 7.479 2 0.0238 
# rumination
linear.rrs.int <- lmer(RRS ~ time + (1 | Respondent_ID),
                       data = long)
anova(linear.rrs.int, linear.rrs) %>%
  as.data.frame() %>%
  knitr::kable(caption = "Rumination: 
               Model comparison for individual differences in change") %>%
  kable_styling()
Rumination: Model comparison for individual differences in change
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
linear.rrs.int 4 885.8 903.9 -438.9 877.8 NA NA NA
linear.rrs 6 883.0 910.1 -435.5 871.0 6.841 2 0.0327 
# mindfulness
linear.mind.int <- lmer(MIND ~ time + (1 | Respondent_ID),
                        data = long)
anova(linear.mind.int, linear.mind) %>%
  as.data.frame() %>%
  knitr::kable(caption = "Mindfulness: 
               Model comparison for individual differences in change") %>%
  kable_styling()
Mindfulness: Model comparison for individual differences in change
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
linear.mind.int 4 726.1 744.2 -359.1 718.1 NA NA NA
linear.mind 6 727.4 754.5 -357.7 715.4 2.75 2 0.2528

To test for individual differences in change, we fit a separate set of models without random slopes. These models were compared to previous models that included both random intercepts and slopes. A significant likelihood ratio test signified significant individual differences in change. This was found for COVID stress, depression, and rumination.

Random-intercept cross-lagged panel model

RICLPM <- '
  
  #random intercepts
  RIx =~ 1*COVID_T1 + 1*COVID_T2 + 1*COVID_T3 + 1*COVID_T4
  RIm =~ 1*RRS_T1 + 1*RRS_T2 + 1*RRS_T3 + 1*RRS_T4
  RIy =~ 1*DAS_D_T1 + 1*DAS_D_T2 + 1*DAS_D_T3 + 1*DAS_D_T4
  
  #within-subject
  wx1 =~ 1*COVID_T1
  wx2 =~ 1*COVID_T2
  wx3 =~ 1*COVID_T3
  wx4 =~ 1*COVID_T4
  wy1 =~ 1*DAS_D_T1
  wy2 =~ 1*DAS_D_T2
  wy3 =~ 1*DAS_D_T3
  wy4 =~ 1*DAS_D_T4 
  wm1 =~ 1*RRS_T1
  wm2 =~ 1*RRS_T2
  wm3 =~ 1*RRS_T3
  wm4 =~ 1*RRS_T4
  
  #cross-lagged effects
  wx2 + wm2 + wy2 ~ wx1 + wm1 + wy1
  wx3 + wm3 + wy3 ~ wx2 + wm2 + wy2
  wx4 + wm4 + wy4 ~ wx3 + wm3 + wy3
  
  #covariances
  wx1 ~~ wy1
  wy1 ~~ wm1
  wx1 ~~ wm1
  wx2 ~~ wy2
  wy2 ~~ wm2
  wx2 ~~ wm2
  wx3 ~~ wy3
  wy3 ~~ wm3
  wx3 ~~ wm3
  wx4 ~~ wy4
  wy4 ~~ wm4
  wx4 ~~ wm4
  RIx ~~ RIy
  RIy ~~ RIm
  RIx ~~ RIm
  
  #variances and residual variances
  RIx ~~ RIx
  RIy ~~ RIy
  RIm ~~ RIm
  wx1 ~~ wx1
  wy1 ~~ wy1 
  wm1 ~~ wm1
  wx2 ~~ wx2
  wy2 ~~ wy2 
  wm2 ~~ wm2
  wx3 ~~ wx3
  wy3 ~~ wy3 
  wm3 ~~ wm3
  wx4 ~~ wx4
  wy4 ~~ wy4 
  wm4 ~~ wm4
  '
  
RICLPM.fit <- lavaan(RICLPM, data = full, missing = 'ML', 
                     meanstructure = T, int.ov.free = T) 

fitMeasures(RICLPM.fit, c("chisq", "df", "pvalue", "cfi", "rmsea", "srmr"), 
            output = "text")
## 
## Model Test User Model:
## 
##   Test statistic                                25.560
##   Degrees of freedom                                21
##   P-value                                        0.224
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.997
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.033
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.028
as.data.frame(unclass(standardizedSolution(RICLPM.fit))) %>% 
  filter(op == "~") %>% 
  knitr::kable(caption = "Cross-lagged paths", digits = 3) %>% 
  kable_styling() %>% 
  kableExtra::footnote("x = COVID stress, y = Depression, m = Rumination")
Cross-lagged paths
lhs op rhs est.std se z pvalue ci.lower ci.upper
wx2 ~ wx1 -0.041 0.391 -0.104 0.917 -0.806 0.725
wx2 ~ wm1 0.353 0.244 1.448 0.148 -0.125 0.832
wx2 ~ wy1 0.090 0.193 0.466 0.642 -0.289 0.469
wm2 ~ wx1 0.407 0.247 1.647 0.100 -0.077 0.892
wm2 ~ wm1 -0.012 0.201 -0.061 0.951 -0.406 0.381
wm2 ~ wy1 0.090 0.147 0.613 0.540 -0.198 0.379
wy2 ~ wx1 0.058 0.218 0.268 0.788 -0.368 0.485
wy2 ~ wm1 -0.172 0.156 -1.097 0.273 -0.478 0.135
wy2 ~ wy1 0.405 0.122 3.331 0.001 0.167 0.643
wx3 ~ wx2 0.347 0.152 2.279 0.023 0.049 0.645
wx3 ~ wm2 0.084 0.157 0.534 0.594 -0.224 0.392
wx3 ~ wy2 -0.130 0.138 -0.938 0.348 -0.401 0.141
wm3 ~ wx2 0.168 0.163 1.031 0.303 -0.151 0.487
wm3 ~ wm2 -0.074 0.188 -0.392 0.695 -0.442 0.295
wm3 ~ wy2 0.038 0.174 0.221 0.825 -0.302 0.379
wy3 ~ wx2 -0.026 0.126 -0.205 0.838 -0.273 0.222
wy3 ~ wm2 -0.170 0.140 -1.212 0.226 -0.444 0.105
wy3 ~ wy2 0.369 0.145 2.551 0.011 0.086 0.653
wx4 ~ wx3 0.689 0.073 9.419 0.000 0.546 0.833
wx4 ~ wm3 -0.014 0.091 -0.155 0.876 -0.193 0.164
wx4 ~ wy3 -0.057 0.087 -0.657 0.511 -0.227 0.113
wm4 ~ wx3 -0.161 0.160 -1.004 0.315 -0.474 0.153
wm4 ~ wm3 0.009 0.196 0.045 0.964 -0.375 0.393
wm4 ~ wy3 0.171 0.163 1.048 0.294 -0.148 0.489
wy4 ~ wx3 -0.050 0.132 -0.381 0.703 -0.309 0.209
wy4 ~ wm3 -0.096 0.145 -0.663 0.507 -0.381 0.188
wy4 ~ wy3 0.267 0.141 1.889 0.059 -0.010 0.543
Note:
x = COVID stress, y = Depression, m = Rumination 
semPaths(RICLPM.fit, what = "col", 
         structural = TRUE,
         residuals = FALSE,
         whatLabels = "est", intercepts = T)

IV. Adding change to mediation analyses

# extract individual slopes for each variable as a data frame
slope_COVID <- data.frame("id" = rownames(coef(linear.covid)$Respondent_ID),
                          "slope_COVID" = coef(linear.covid)$Respondent_ID[,2])
slope_DASD <- data.frame("id" = rownames(coef(linear.dasd)$Respondent_ID),
                         "slope_DASD" = coef(linear.dasd)$Respondent_ID[,2])
slope_RRS <- data.frame("id" = rownames(coef(linear.rrs)$Respondent_ID),
                        "slope_RRS" = coef(linear.rrs)$Respondent_ID[,2])
slope_MIND <- data.frame("id" = rownames(coef(linear.mind)$Respondent_ID),
                         "slope_MIND" = coef(linear.mind)$Respondent_ID[,2])

# merge all slopes into one data frame to use in mediation analysis
slope <- merge(slope_COVID, slope_DASD) %>%
  merge(slope_RRS) %>%
  merge(slope_MIND)

# bivariate slope correlation
slope %>% select(-id) %>% apa.cor.table(file = "correlated slopes.docx")
## 
## 
## Means, standard deviations, and correlations with confidence intervals
##  
## 
##   Variable       M     SD   1           2            3          
##   1. slope_COVID -0.15 0.08                                     
##                                                                 
##   2. slope_DASD  -0.03 0.06 .14*                                
##                             [.00, .28]                          
##                                                                 
##   3. slope_RRS   -0.05 0.05 .06         .46**                   
##                             [-.08, .20] [.34, .56]              
##                                                                 
##   4. slope_MIND  -0.01 0.02 -.06        -.16*        -.06       
##                             [-.20, .08] [-.29, -.02] [-.20, .08]
##                                                                 
## 
## Note. M and SD are used to represent mean and standard deviation, respectively.
## Values in square brackets indicate the 95% confidence interval.
## The confidence interval is a plausible range of population correlations 
## that could have caused the sample correlation (Cumming, 2014).
##  * indicates p < .05. ** indicates p < .01.
## 
cormat <- round(x = cor(select(slope, c("slope_COVID", "slope_RRS", 
                                        "slope_MIND", "slope_DASD"))),
                digits = 2)
cormat[lower.tri(cormat)] <- NA
cormat <- melt(cormat, na.rm = TRUE) %>%
  # delete the perfect correlations
  filter(value != 1.00)
ggplot(data = cormat, aes(Var2, Var1, fill = value)) +
 geom_tile(color = "white" ) +
 scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                      midpoint = 0, limit = c(-1,1), space = "Lab", 
                      name = "Pearson\nCorrelation") +
  theme_minimal()+ 
  labs(x = NULL, y = NULL) +
  coord_fixed() +
  geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.major = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    axis.ticks = element_blank(),
    legend.justification = c(1, 0),
    legend.position = c(0.4, 0.7),
    legend.direction = "horizontal")+
  guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
                               title.position = "top", title.hjust = 0.5))

Mediation - Slope

process(data = slope,
        y = "slope_DASD",
        x = "slope_COVID",
        m = "slope_RRS",
        model = 4,
        boot = 5000,
        seed = 20220308)
a <- c(0.0358, 0.8373, 0.4034, -0.0485 , 0.1200)
b <- c(0.5758, 7.0970, "<0.001", 0.4158, 0.7358)
c <- c(0.0890, 1.8409, 0.0672,  -0.0064, 0.1843)
indirect <- c(0.0206, NA, NA, -0.0267, 0.0719)

(med_slope <- med_table(a = a, b = b, c = c, indirect = indirect))

Moderated mediation - Slope

process(data = slope,
        y = "slope_DASD",
        x = "slope_COVID",
        m = "slope_RRS",
        w = "slope_MIND",
        model = 7,
        moments = 1,
        center = 2, # continuous variables in interaction terms
        boot = 5000,
        seed = 20220309)
a1 <- c( 0.0485,  1.1183, 0.2649)
a2 <- c(-0.1652, -0.8977, 0.3704)
a3 <- c(-3.4274, -1.8297, 0.0689)
cp <- c(0.0890, 1.8409, 0.0672)
b1 <- c(0.5758, 7.0970, "<0.001")
R2a <- 0.0236 
R2b <- 0.2234
Fpa <- 0.2045
Fpb <- "<0.001"
low <- c(-0.0184,    0.0643,    0.0356,   -0.0044,    0.1383)
mean <- c(0.0000,    0.0279,    0.0245,   -0.0192,    0.0791)
high <- c(0.0184,   -0.0085,    0.0278,   -0.0556,    0.0556)

(modmed_slope <- modmed_table(a1 = a1, a2 = a2, a3 = a3, cp = cp, b1 = b1,
                              R2a = R2a, R2b = R2b, Fpa = Fpa, Fpb = Fpb))
(modmed_slope_condl <- modmed_condl(low = low, mean = mean, high = high))
# export all mediation and moderated/mediation tables to docx
save_as_docx(`Average - Mediation` = med_avg, 
             `Average - Moderated/Mediation` = modmed_avg,
             `Average - Moderated/Mediation Conditional` = modmed_avg_condl,
             `Wave 1 - Mediation` = med_w1, 
             `Wave 1 - Moderated/Mediation` = modmed_w1,
             `Wave 1 - Moderated/Mediation Conditional` = modmed_w1_condl,
             `Wave 2 - Mediation` = med_w2, 
             `Wave 2 - Moderated/Mediation` = modmed_w2,
             `Wave 2 - Moderated/Mediation Conditional` = modmed_w2_condl,
             `Wave 3 - Mediation` = med_w3, 
             `Wave 3 - Moderated/Mediation` = modmed_w3,
             `Wave 3 - Moderated/Mediation Conditional` = modmed_w3_condl,
             `Wave 4 - Mediation` = med_w4, 
             `Wave 4 - Moderated/Mediation` = modmed_w4,
             `Wave 4 - Moderated/Mediation Conditional` = modmed_w4_condl,
             `Shortened RRS-10 - Mediation` = med_rrs10, 
             `Shortened RRS-10 - Moderated/Mediation` = modmed_rrs10,
             `Shortened RRS-10 - Moderated/Mediation Conditional` = modmed_rrs10_condl,
             `Slope - Mediation` = med_slope, 
             `Slope - Moderated/Mediation` = modmed_slope,
             `Slope - Moderated/Mediation Conditional` = modmed_slope_condl,
             path = "table.docx")