# The following packages might need to be installed onto your version
# of R prior to the running of the code below.
# Package names
packages <- c("MASS", "lavaan", "processR", "corrplot", "tidytext", "tidyverse", "haven", "jmv", "Hmisc", "cluster", "kableExtra", "factoextra")
# Install packages not yet installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages])
}
# Packages loading
invisible(lapply(packages, library, character.only = TRUE))
# We load the original Czech dataset (in SPSS format) from a local directory.
data <- zap_labels(haven::read_sav(file = "INTERNATIONAL DATASET_12 Oct_final.sav"))
# For use in correlation analysis, we duplicate the dataset under name data_corr
data_corr <- data
# We also try to limit the decimals to three significant figures
options(digits = 3, scipen = 999)
# Calculating PH8
data <- data %>%
mutate(Q49_1_rec = recode(Q49_1,
`1` = 0,
`2` = 1,
`3` = 2,
`4` = 3),
Q49_2_rec = recode(Q49_2,
`1` = 0,
`2` = 1,
`3` = 2,
`4` = 3),
Q49_3_rec = recode(Q49_3,
`1` = 0,
`2` = 1,
`3` = 2,
`4` = 3),
Q49_4_rec = recode(Q49_4,
`1` = 0,
`2` = 1,
`3` = 2,
`4` = 3),
Q49_5_rec = recode(Q49_5,
`1` = 0,
`2` = 1,
`3` = 2,
`4` = 3),
Q49_6_rec = recode(Q49_6,
`1` = 0,
`2` = 1,
`3` = 2,
`4` = 3),
Q49_7_rec = recode(Q49_7,
`1` = 0,
`2` = 1,
`3` = 2,
`4` = 3),
Q49_8_rec = recode(Q49_8,
`1` = 0,
`2` = 1,
`3` = 2,
`4` = 3),
PHQ8 = Q49_1_rec + Q49_2_rec + Q49_3_rec + Q49_4_rec + Q49_5_rec + Q49_6_rec + Q49_7_rec + Q49_8_rec)
# Renaming variables to identical ones as before
data_full <- data %>%
transmute(country = recode_factor(as_factor(Country),
`1` = "IT",
`2` = "FR",
`3` = "PL",
`4` = "UK",
`5` = "CZ",
`6` = "SW"),
q01_gender = recode_factor(as_factor(Q3),
`1` = "male",
`2` = "female",
`3` = "other"),
q02_age = Q4,
q02_age_group = recode_factor(as_factor(Q4_AGE_r),
`1` = "16-29 years",
`2` = "30-49 years",
`3` = "50-64 years",
`4` = "65+"),
q03_relationship_type = recode_factor(as_factor(Q5),
`1` = "single",
`2` = "relationship",
`3` = "married",
`4` = "divorced",
`5` = "widowed"),
q04_children = recode_factor(as_factor(Q6),
`1` = "yes",
`2` = "no"),
q11_education = recode_factor(as_factor(Q13R),
`1` = "low",
`2` = "medium",
`3` = "high"),
q18_02_soc_media = recode_factor(as_factor(replace_na(Q21_2, 0)),
`0` = "no",
`1` = "yes"),
q20_public_info = recode_factor(as_factor(Q23),
`1` = "yes",
`2` = "no",
`99` = "do_not_know"),
q34_02_face_mask = recode_factor(as_factor(Q37_2),
`1` = "yes",
`2` = "no"),
q34_07_hand_washing = recode_factor(as_factor(Q37_7),
`1` = "yes",
`2` = "no"),
q35_01_contact_close_family = recode_factor(as_factor(Q38_1),
`1` = "less_often",
`2` = "as_before",
`3` = "more_often"),
q35_03_contact_friends = recode_factor(as_factor(Q38_3),
`1` = "less_often",
`2` = "as_before",
`3` = "more_often"),
q36_econ_worry = recode_factor(as_factor(Q39),
`1` = "very_serious",
`2` = "serious",
`3` = "limited"),
q38_alcohol = recode_factor(as_factor(Q42),
`1` = "yes",
`2` = "no"),
q40_smoking = recode_factor(as_factor(Q44),
`1` = "yes",
`2` = "no"),
q42_sport = recode_factor(as_factor(Q46),
`1` = "yes",
`2` = "no"),
q47_self_reporting_health = recode_factor(as_factor(Q50),
`1` = "excellent",
`2` = "good",
`3` = "neutral",
`4` = "bad",
`5` = "very_bad"),
q48_chronic_illness = recode_factor(as_factor(Q51),
`1` = "yes",
`2` = "no"),
q49_health_limitations = recode_factor(as_factor(Q52),
`1` = "limits",
`2` = "partially_limits",
`3` = "no_limits"),
q30_concern_infection_covid = Q33,
q31_concern_infection_friends = Q34,
q33_01_concern_situation = Q36_1,
q33_02_concern_low_control = Q36_2,
q33_03_concern_survival_covid = Q36_3,
q33_04_concern_change_employment = Q36_4,
q33_05_concern_infecting_others = Q36_5,
PHQ8 = PHQ8)
# Filtering out Czech sample for analyses
data <- data_full %>%
filter(country != "CZ")
data_full[data_full$country %in% c("IT", "FR", "PL") & !is.na(data_full$q02_age),] %>%
ggplot(aes(x = q02_age,
fill = q01_gender)) +
geom_bar(data = data_full[data_full$country %in% c("IT", "FR", "PL") & data_full$q01_gender == "female" & !is.na(data_full$q02_age),]) +
geom_bar(data = data_full[data_full$country %in% c("IT", "FR", "PL") & data_full$q01_gender == "male" & !is.na(data_full$q02_age),],
aes(y = ..count..*(-1))) +
scale_y_continuous(breaks = seq(-250, 250, 50),
labels = abs(seq(-250, 250, 50)),
limits = c(-100, 100)) +
scale_x_continuous(breaks = seq(15, 100, 5),
labels = seq(15, 100, 5),
limits = c(15, 100)) +
coord_flip() +
theme_bw() +
theme(legend.title = element_blank()) +
facet_grid(cols = vars(country)) +
scale_fill_manual(values = c("#E64B35CC", "#00A087CC")) +
xlab("Age") +
ylab("Count")
data_full[data_full$country %in% c("UK", "CZ", "SW") & !is.na(data_full$q02_age),] %>%
ggplot(aes(x = q02_age,
fill = q01_gender)) +
geom_bar(data = data_full[data_full$country %in% c("UK", "CZ", "SW") & data_full$q01_gender == "female" & !is.na(data_full$q02_age),]) +
geom_bar(data = data_full[data_full$country %in% c("UK", "CZ", "SW") & data_full$q01_gender == "male" & !is.na(data_full$q02_age),],
aes(y = ..count..*(-1))) +
scale_y_continuous(breaks = seq(-250, 250, 50),
labels = abs(seq(-250, 250, 50)),
limits = c(-100, 100)) +
scale_x_continuous(breaks = seq(15, 100, 5),
labels = seq(15, 100, 5),
limits = c(15, 100)) +
coord_flip() +
theme_bw() +
theme(legend.title = element_blank()) +
facet_grid(cols = vars(country)) +
scale_fill_manual(values = c("#E64B35CC", "#00A087CC")) +
xlab("Age") +
ylab("Count")
# To summarize the dependent continuous variable, we use the descriptives()
# function from the jmv package.
descriptives <- jmv::descriptives(
data = data,
vars = "PHQ8",
freq = TRUE,
box = TRUE,
median = FALSE,
range = TRUE,
sd = TRUE,
pc = TRUE)
descriptives$descriptives
Descriptives
──────────────────────────────
PHQ8
──────────────────────────────
N 7804
Missing 187
Mean 6.91
Standard deviation 5.45
Range 24.0
Minimum 0.00
Maximum 24.0
25th percentile 3.00
50th percentile 6.00
75th percentile 10.0
──────────────────────────────
data %>%
filter(PHQ8 != "NA") %>%
ggplot(aes(y = PHQ8)) +
geom_boxplot(fill = "#4DBBD5CC") +
theme_bw() +
scale_y_continuous(breaks = seq(0, 24, 1),
labels = seq(0, 24, 1),
limits = c(0, 24)) +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank())
data %>%
filter(PHQ8 != "NA") %>%
ggplot(aes(x = as_factor(PHQ8))) +
geom_bar(fill = "#4DBBD5CC") +
scale_y_continuous(breaks = seq(0, 700, 50),
labels = seq(0, 700, 50),
limits = c(0, 700)) +
theme_bw() +
theme(legend.position = "none") +
xlab("PHQ8")
# To summarize the key demographic variables, we use the descriptives()
# function from the jmv package.
demo_descriptives <- jmv::descriptives(
data = data,
vars = vars("q01_gender",
"q02_age_group",
"q03_relationship_type",
"q04_children",
"q11_education"),
bar = TRUE,
freq = TRUE,
missing = FALSE,
mean = FALSE,
median = FALSE,
sd = FALSE,
min = FALSE,
max = FALSE)
demo_descriptives$frequencies
##
## FREQUENCIES
##
## Frequencies of q01_gender
## ──────────────────────────────────────────────────
## Levels Counts % of Total Cumulative %
## ──────────────────────────────────────────────────
## male 1738 21.8 21.8
## female 6242 78.2 99.9
## other 7 8.76e-4 100.0
## ──────────────────────────────────────────────────
##
##
## Frequencies of q02_age_group
## ───────────────────────────────────────────────────────
## Levels Counts % of Total Cumulative %
## ───────────────────────────────────────────────────────
## 16-29 years 1655 20.7 20.7
## 30-49 years 4177 52.3 73.0
## 50-64 years 1683 21.1 94.1
## 65+ 471 5.9 100.0
## ───────────────────────────────────────────────────────
##
##
## Frequencies of q03_relationship_type
## ────────────────────────────────────────────────────────
## Levels Counts % of Total Cumulative %
## ────────────────────────────────────────────────────────
## single 1816 22.8 22.8
## relationship 1930 24.2 46.9
## married 3490 43.7 90.7
## divorced 621 7.8 98.4
## widowed 125 1.6 100.0
## ────────────────────────────────────────────────────────
##
##
## Frequencies of q04_children
## ──────────────────────────────────────────────────
## Levels Counts % of Total Cumulative %
## ──────────────────────────────────────────────────
## yes 4745 59.4 59.4
## no 3242 40.6 100.0
## ──────────────────────────────────────────────────
##
##
## Frequencies of q11_education
## ──────────────────────────────────────────────────
## Levels Counts % of Total Cumulative %
## ──────────────────────────────────────────────────
## low 1011 17.3 17.3
## medium 1712 29.3 46.6
## high 3116 53.4 100.0
## ──────────────────────────────────────────────────
data %>%
filter(q01_gender %in% c("male", "female"),
q02_age_group != "NA") %>%
ggplot(aes(x = q01_gender)) +
geom_bar(aes(fill = q01_gender)) +
scale_y_continuous(breaks = seq(0, 3500, 500),
labels = seq(0, 3500, 500),
limits = c(0, 3500)) +
theme_bw() +
theme(legend.position = "none", axis.title.x = element_blank()) +
facet_grid(cols = vars(q02_age_group)) +
scale_fill_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC"))
data %>%
filter(q01_gender %in% c("male", "female"),
q03_relationship_type != "NA") %>%
ggplot(aes(x = q01_gender)) +
geom_bar(aes(fill = q01_gender)) +
scale_y_continuous(breaks = seq(0, 2500, 250),
labels = seq(0, 2500, 250),
limits = c(0, 2500)) +
theme_bw() +
theme(legend.position = "none", axis.title.x = element_blank()) +
facet_grid(cols = vars(q03_relationship_type),
rows = vars(q04_children)) +
scale_fill_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC"))
data %>%
filter(q01_gender %in% c("male", "female"),
q11_education != "NA",
q02_age_group != "NA") %>%
ggplot(aes(x = q01_gender)) +
geom_bar(aes(fill = q01_gender)) +
scale_y_continuous(breaks = seq(0, 1600, 250),
labels = seq(0, 1600, 250),
limits = c(0, 1600)) +
theme_bw() +
theme(legend.position = "none", axis.title.x = element_blank()) +
facet_grid(cols = vars(q11_education),
rows = vars(q02_age_group)) +
scale_fill_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC"))
data_full %>%
filter(q01_gender != "other",
PHQ8 != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.x = element_blank()) +
facet_grid(cols = vars(country)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Countries")
data %>%
filter(q01_gender != "other",
q02_age_group != "NA",
PHQ8 != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.x = element_blank()) +
facet_grid(cols = vars(q02_age_group)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Age groups")
# For reproducibility of results
set.seed(2021)
data %>%
filter(q01_gender != "other",
q02_age_group != "NA",
PHQ8 != "NA") %>%
ggplot(aes(x = q02_age,
y = PHQ8)) +
geom_jitter(na.rm = TRUE,
size = 1,
alpha = 0.1) +
stat_smooth(method = "gam",
formula = y ~ x,
na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.title = element_blank()) +
scale_x_continuous(limits = c(16, 92),
breaks = seq(from = 16,
to = 92,
by = 8)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
xlab("Age") +
ggtitle("Age groups")
data %>%
filter(q01_gender != "other",
q03_relationship_type != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.x = element_blank()) +
facet_grid(cols = vars(q03_relationship_type)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Relationship type")
data %>%
filter(q01_gender != "other",
q11_education != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.x = element_blank()) +
facet_grid(cols = vars(q11_education)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Educational attainment")
data %>%
filter(q01_gender != "other") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.x = element_blank()) +
facet_grid(cols = vars(q04_children)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Children")
data %>%
filter(q01_gender != "other",
q35_01_contact_close_family != "NA",
q02_age_group != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.y = element_blank()) +
facet_grid(cols = vars(q35_01_contact_close_family), rows = vars(q02_age_group)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Frequency of physical contact - close family") +
coord_flip()
data %>%
filter(q01_gender != "other",
q35_03_contact_friends != "NA",
q02_age_group != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.y = element_blank()) +
facet_grid(cols = vars(q35_03_contact_friends), rows = vars(q02_age_group)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Frequency of physical contact - friends") +
coord_flip()
data %>%
filter(q01_gender != "other",
q38_alcohol != "NA",
q02_age_group != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.y = element_blank()) +
facet_grid(cols = vars(q38_alcohol), rows = vars(q02_age_group)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Alcohol") +
coord_flip()
data %>%
filter(q01_gender != "other",
q40_smoking != "NA",
q02_age_group != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.y = element_blank()) +
facet_grid(cols = vars(q40_smoking), rows = vars(q02_age_group)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Smoking") +
coord_flip()
data %>%
filter(q01_gender != "other",
q42_sport != "NA",
q02_age_group != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.y = element_blank()) +
facet_grid(cols = vars(q42_sport ), rows = vars(q02_age_group)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Sport") +
coord_flip()
data %>%
filter(q01_gender != "other",
q47_self_reporting_health != "NA",
q02_age_group != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.y = element_blank()) +
facet_grid(cols = vars(q47_self_reporting_health), rows = vars(q02_age_group)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Self-reported Health") +
coord_flip()
data %>%
filter(q01_gender != "other",
q48_chronic_illness != "NA",
q02_age_group != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.y = element_blank()) +
facet_grid(cols = vars(q48_chronic_illness), rows = vars(q02_age_group)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Chronic illness") +
coord_flip()
data %>%
filter(q01_gender != "other",
q49_health_limitations != "NA",
q02_age_group != "NA") %>%
ggplot(aes(y = PHQ8,
x = q01_gender)) +
geom_boxplot(na.rm = TRUE,
aes(color = q01_gender)) +
theme_bw() +
theme(legend.position = "none", axis.title.y = element_blank()) +
facet_grid(cols = vars(q49_health_limitations ), rows = vars(q02_age_group)) +
scale_color_manual(values = c("female" = "#E64B35CC", "male" = "#00A087CC")) +
ggtitle("Health limitations") +
coord_flip()
# While the dataset has been already imported, the values of factor variables
# were changed from numerics to text strings, therefore that dataset is unsuitable
# for correlation analysis. To solve this, we create a parallel dataset,
# again renaming the key variables to a more understandable form.
data_corr <- data_corr %>%
transmute(country = Country,
q01_gender = Q3,
q02_age = Q4,
q03_relationship_type = Q5,
q04_children = Q6,
q11_education = Q13R,
q18_02_soc_media = replace_na(Q21_2, 0),
q20_public_info = Q23,
q34_02_face_mask = Q37_2,
q34_07_hand_washing = Q37_7,
q36_econ_worry = Q39,
q35_01_contact_close_family = Q38_1,
q35_03_contact_friends = Q38_3,
q38_alcohol = Q42,
q40_smoking = Q44,
q42_sport = Q46,
q47_self_reporting_health = Q50,
q48_chronic_illness = Q51,
q49_health_limitations = Q52) %>%
filter(country != 5) %>%
select(-country)
data_corr <- cbind(data_corr, PHQ8_t)
res1 <- cor.mtest(data_corr, conf.level = .95)
#Correlation matrix using Spearman coefficient (values with p > 0.05 are crossed)
corrplot(cor(data_corr,
method = "spearman",
use = "complete.obs"),
method = "circle",
title = "Correlation Matrix - Spearman Coefficient",
type = "lower",
p.mat = res1$p,
sig.level = .05,
mar = c(0,0,1,0))
linreg_theory <- jmv::linReg(
data = data,
dep = "PHQ8_t",
covs = "q02_age",
factors = vars("q01_gender",
"q03_relationship_type",
"q04_children",
# "q11_education",
"q18_02_soc_media",
"q20_public_info",
"q34_02_face_mask",
"q34_07_hand_washing",
"q35_01_contact_close_family",
"q35_03_contact_friends",
"q36_econ_worry",
"q38_alcohol",
"q40_smoking",
"q42_sport",
"q47_self_reporting_health",
"q48_chronic_illness",
"q49_health_limitations"),
blocks = list(
list(
"q01_gender",
"q02_age",
"q03_relationship_type",
"q04_children"),
# "q11_education"),
list(
"q18_02_soc_media",
"q20_public_info",
"q34_02_face_mask",
"q34_07_hand_washing",
"q36_econ_worry"),
list(
"q40_smoking",
"q42_sport",
"q38_alcohol",
"q35_01_contact_close_family",
"q35_03_contact_friends"),
list(
"q47_self_reporting_health",
"q48_chronic_illness",
"q49_health_limitations")),
refLevels = list(
list(
var = "q01_gender",
ref = "female"),
list(
var = "q04_children",
ref = "no"),
list(
var = "q20_public_info",
ref = "no"),
list(
var = "q34_02_face_mask",
ref = "no"),
list(
var = "q34_07_hand_washing",
ref = "no"),
list(
var = "q36_econ_worry",
ref = "very_serious"),
list(
var = "q42_sport",
ref = "no"),
list(
var = "q40_smoking",
ref = "yes"),
list(
var = "q38_alcohol",
ref = "yes"),
list(
var = "q35_01_contact_close_family",
ref = "less_often"),
list(
var = "q35_03_contact_friends",
ref = "less_often"),
list(
var = "q18_02_soc_media",
ref = "yes"),
list(
var = "q03_relationship_type",
ref = "single"),
list(
var = "q47_self_reporting_health",
ref = "very_bad"),
list(
var = "q49_health_limitations",
ref = "limits"),
# list(
# var = "q11_education",
# ref = "low"),
list(
var = "q48_chronic_illness",
ref = "yes")),
r2Adj = TRUE,
aic = TRUE,
bic = TRUE,
rmse = TRUE,
modelTest = TRUE,
anova = TRUE,
ci = TRUE,
stdEst = TRUE,
ciStdEst = TRUE,
durbin = TRUE,
collin = TRUE)
linreg_theory$modelFit
Model Fit Measures
──────────────────────────────────────────────────────────────────────────────────────────────────────
Model R R² Adjusted R² AIC BIC RMSE F df1 df2 p
──────────────────────────────────────────────────────────────────────────────────────────────────────
1 0.274 0.0749 0.0740 27219 27289 1.43 77.3 8 7633 < .001
2 0.356 0.1267 0.1249 26794 26912 1.39 73.7 15 7626 < .001
3 0.369 0.1362 0.1337 26724 26891 1.39 54.6 22 7619 < .001
4 0.488 0.2384 0.2355 25775 25990 1.30 82.2 29 7612 < .001
──────────────────────────────────────────────────────────────────────────────────────────────────────
linreg_theory$modelComp
Model Comparisons
────────────────────────────────────────────────────────────────────
Model Model ΔR² F df1 df2 p
────────────────────────────────────────────────────────────────────
1 - 2 0.05173 64.5 7 7626 < .001
2 - 3 0.00951 12.0 7 7619 < .001
3 - 4 0.10228 146.0 7 7612 < .001
────────────────────────────────────────────────────────────────────
linreg_theory$models
MODEL SPECIFIC RESULTS
MODEL 1
Omnibus ANOVA Test
──────────────────────────────────────────────────────────────────────────────────────
Sum of Squares df Mean Square F p
──────────────────────────────────────────────────────────────────────────────────────
q01_gender 393.14 2 196.57 95.44 < .001
q02_age 579.47 1 579.47 281.35 < .001
q03_relationship_type 59.23 4 14.81 7.19 < .001
q04_children 4.12 1 4.12 2.00 0.157
Residuals 15720.82 7633 2.06
──────────────────────────────────────────────────────────────────────────────────────
Note. Type 3 sum of squares
Model Coefficients - PHQ8_t
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Predictor Estimate SE Lower Upper t p Stand. Estimate Lower Upper
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Intercept ᵃ 3.9502 0.05909 3.8344 4.0661 66.851 < .001
q01_gender:
male – female -0.5557 0.04049 -0.6351 -0.4763 -13.725 < .001 -0.3726 -0.4258 -0.3194
other – female 0.9239 0.64223 -0.3351 2.1829 1.439 0.150 0.6195 -0.2247 1.4637
q02_age -0.0238 0.00142 -0.0266 -0.0210 -16.774 < .001 -0.2129 -0.2378 -0.1880
q03_relationship_type:
relationship – single -0.1705 0.05057 -0.2696 -0.0714 -3.371 < .001 -0.1143 -0.1808 -0.0479
married – single -0.1534 0.05337 -0.2580 -0.0488 -2.874 0.004 -0.1028 -0.1730 -0.0327
divorced – single 0.0992 0.07715 -0.0521 0.2504 1.285 0.199 0.0665 -0.0349 0.1679
widowed – single 0.0817 0.14991 -0.2122 0.3756 0.545 0.586 0.0548 -0.1423 0.2518
q04_children:
yes – no -0.0585 0.04135 -0.1396 0.0225 -1.415 0.157 -0.0392 -0.0936 0.0151
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
ᵃ Represents reference level
ASSUMPTION CHECKS
Durbin–Watson Test for Autocorrelation
────────────────────────────────────────────
Autocorrelation DW Statistic p
────────────────────────────────────────────
0.0325 1.93 0.008
────────────────────────────────────────────
Collinearity Statistics
──────────────────────────────────────────────
VIF Tolerance
──────────────────────────────────────────────
q01_gender 1.01 0.992
q02_age 1.15 0.867
q03_relationship_type 1.08 0.930
q04_children 1.24 0.808
──────────────────────────────────────────────
MODEL 2
Omnibus ANOVA Test
─────────────────────────────────────────────────────────────────────────────────────────
Sum of Squares df Mean Square F p
─────────────────────────────────────────────────────────────────────────────────────────
q01_gender 333.56504 2 166.78252 85.69636 < .001
q02_age 435.39966 1 435.39966 223.71749 < .001
q03_relationship_type 45.57174 4 11.39293 5.85393 < .001
q04_children 9.08554 1 9.08554 4.66834 0.031
q18_02_soc_media 54.79480 1 54.79480 28.15472 < .001
q20_public_info 132.17343 2 66.08671 33.95674 < .001
q34_02_face_mask 0.00213 1 0.00213 0.00109 0.974
q34_07_hand_washing 20.29233 1 20.29233 10.42663 0.001
q36_econ_worry 578.14533 2 289.07266 148.53161 < .001
Residuals 14841.74441 7626 1.94620
─────────────────────────────────────────────────────────────────────────────────────────
Note. Type 3 sum of squares
Model Coefficients - PHQ8_t
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Predictor Estimate SE Lower Upper t p Stand. Estimate Lower Upper
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Intercept ᵃ 4.21742 0.11257 3.9967 4.43809 37.4644 < .001
q01_gender:
male – female -0.51516 0.03966 -0.5929 -0.43741 -12.9887 < .001 -0.34544 -0.3976 -0.29330
other – female 0.95339 0.62457 -0.2709 2.17773 1.5265 0.127 0.63929 -0.1817 1.46026
q02_age -0.02095 0.00140 -0.0237 -0.01820 -14.9572 < .001 -0.18728 -0.2118 -0.16274
q03_relationship_type:
relationship – single -0.17703 0.04920 -0.2735 -0.08058 -3.5981 < .001 -0.11870 -0.1834 -0.05403
married – single -0.12680 0.05199 -0.2287 -0.02489 -2.4391 0.015 -0.08503 -0.1534 -0.01669
divorced – single 0.06902 0.07504 -0.0781 0.21611 0.9198 0.358 0.04628 -0.0524 0.14491
widowed – single -0.00889 0.14585 -0.2948 0.27702 -0.0610 0.951 -0.00596 -0.1977 0.18575
q04_children:
yes – no -0.08694 0.04024 -0.1658 -0.00806 -2.1606 0.031 -0.05830 -0.1112 -0.00541
q18_02_soc_media:
no – yes -0.18010 0.03394 -0.2466 -0.11357 -5.3061 < .001 -0.12077 -0.1654 -0.07615
q20_public_info:
yes – no -0.27747 0.03528 -0.3466 -0.20832 -7.8653 < .001 -0.18606 -0.2324 -0.13969
do_not_know – no -0.03903 0.05864 -0.1540 0.07592 -0.6656 0.506 -0.02617 -0.1032 0.05091
q34_02_face_mask:
yes – no -0.00115 0.03486 -0.0695 0.06718 -0.0331 0.974 -7.73e-4 -0.0466 0.04504
q34_07_hand_washing:
yes – no 0.30711 0.09511 0.1207 0.49355 3.2290 0.001 0.20593 0.0809 0.33095
q36_econ_worry:
serious – very_serious -0.31430 0.04251 -0.3976 -0.23095 -7.3926 < .001 -0.21075 -0.2666 -0.15487
limited – very_serious -0.72547 0.04344 -0.8106 -0.64031 -16.6990 < .001 -0.48646 -0.5436 -0.42935
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
ᵃ Represents reference level
ASSUMPTION CHECKS
Durbin–Watson Test for Autocorrelation
─────────────────────────────────────────────
Autocorrelation DW Statistic p
─────────────────────────────────────────────
0.0359 1.93 < .001
─────────────────────────────────────────────
Collinearity Statistics
──────────────────────────────────────────────
VIF Tolerance
──────────────────────────────────────────────
q01_gender 1.01 0.988
q02_age 1.17 0.855
q03_relationship_type 1.08 0.928
q04_children 1.24 0.808
q18_02_soc_media 1.02 0.982
q20_public_info 1.01 0.992
q34_02_face_mask 1.02 0.982
q34_07_hand_washing 1.01 0.986
q36_econ_worry 1.01 0.991
──────────────────────────────────────────────
MODEL 3
Omnibus ANOVA Test
─────────────────────────────────────────────────────────────────────────────────────────────
Sum of Squares df Mean Square F p
─────────────────────────────────────────────────────────────────────────────────────────────
q01_gender 309.977 2 154.988 80.440 < .001
q02_age 411.513 1 411.513 213.577 < .001
q03_relationship_type 44.404 4 11.101 5.761 < .001
q04_children 7.994 1 7.994 4.149 0.042
q18_02_soc_media 55.925 1 55.925 29.025 < .001
q20_public_info 128.831 2 64.415 33.432 < .001
q34_02_face_mask 0.267 1 0.267 0.139 0.710
q34_07_hand_washing 12.984 1 12.984 6.739 0.009
q36_econ_worry 519.640 2 259.820 134.848 < .001
q40_smoking 8.914 1 8.914 4.627 0.032
q42_sport 23.841 1 23.841 12.374 < .001
q38_alcohol 0.348 1 0.348 0.181 0.671
q35_01_contact_close_family 98.680 2 49.340 25.608 < .001
q35_03_contact_friends 21.667 2 10.833 5.623 0.004
Residuals 14680.048 7619 1.927
─────────────────────────────────────────────────────────────────────────────────────────────
Note. Type 3 sum of squares
Model Coefficients - PHQ8_t
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Predictor Estimate SE Lower Upper t p Stand. Estimate Lower Upper
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Intercept ᵃ 4.4676 0.11994 4.2325 4.70274 37.2482 < .001
q01_gender:
male – female -0.5013 0.03977 -0.5793 -0.42335 -12.6045 < .001 -0.33616 -0.3884 -0.28388
other – female 0.8139 0.62170 -0.4048 2.03255 1.3091 0.191 0.54572 -0.2715 1.36291
q02_age -0.0207 0.00142 -0.0235 -0.01794 -14.6143 < .001 -0.18524 -0.2101 -0.16039
q03_relationship_type:
relationship – single -0.1828 0.04914 -0.2791 -0.08645 -3.7196 < .001 -0.12256 -0.1872 -0.05797
married – single -0.1302 0.05202 -0.2322 -0.02820 -2.5024 0.012 -0.08729 -0.1557 -0.01891
divorced – single 0.0556 0.07472 -0.0909 0.20205 0.7436 0.457 0.03726 -0.0610 0.13548
widowed – single -0.0283 0.14524 -0.3130 0.25646 -0.1946 0.846 -0.01895 -0.2099 0.17196
q04_children:
yes – no -0.0822 0.04033 -0.1612 -0.00309 -2.0369 0.042 -0.05509 -0.1081 -0.00207
q18_02_soc_media:
no – yes -0.1823 0.03383 -0.2486 -0.11595 -5.3875 < .001 -0.12222 -0.1667 -0.07775
q20_public_info:
yes – no -0.2746 0.03514 -0.3435 -0.20575 -7.8158 < .001 -0.18415 -0.2303 -0.13797
do_not_know – no -0.0411 0.05838 -0.1555 0.07338 -0.7032 0.482 -0.02753 -0.1043 0.04921
q34_02_face_mask:
yes – no -0.0131 0.03505 -0.0818 0.05565 -0.3724 0.710 -0.00875 -0.0548 0.03732
q34_07_hand_washing:
yes – no 0.2479 0.09550 0.0607 0.43512 2.5959 0.009 0.16624 0.0407 0.29177
q36_econ_worry:
serious – very_serious -0.3011 0.04243 -0.3843 -0.21790 -7.0954 < .001 -0.20189 -0.2577 -0.14611
limited – very_serious -0.6937 0.04358 -0.7791 -0.60824 -15.9159 < .001 -0.46514 -0.5224 -0.40785
q40_smoking:
no – yes -0.0835 0.03884 -0.1597 -0.00740 -2.1509 0.032 -0.05601 -0.1071 -0.00497
q42_sport:
yes – no -0.1140 0.03241 -0.1776 -0.05048 -3.5176 < .001 -0.07645 -0.1191 -0.03385
q38_alcohol:
no – yes -0.0145 0.03406 -0.0812 0.05230 -0.4249 0.671 -0.00970 -0.0545 0.03507
q35_01_contact_close_family:
as_before – less_often -0.1982 0.03490 -0.2667 -0.12982 -5.6795 < .001 -0.13293 -0.1788 -0.08705
more_often – less_often 0.1040 0.05265 7.94e-4 0.20721 1.9754 0.048 0.06974 5.33e-4 0.13895
q35_03_contact_friends:
as_before – less_often -0.2596 0.07753 -0.4116 -0.10760 -3.3481 < .001 -0.17407 -0.2760 -0.07215
more_often – less_often 0.0138 0.17948 -0.3381 0.36560 0.0767 0.939 0.00923 -0.2267 0.24515
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
ᵃ Represents reference level
ASSUMPTION CHECKS
Durbin–Watson Test for Autocorrelation
────────────────────────────────────────────
Autocorrelation DW Statistic p
────────────────────────────────────────────
0.0331 1.93 0.006
────────────────────────────────────────────
Collinearity Statistics
────────────────────────────────────────────────────
VIF Tolerance
────────────────────────────────────────────────────
q01_gender 1.02 0.984
q02_age 1.19 0.840
q03_relationship_type 1.08 0.924
q04_children 1.25 0.802
q18_02_soc_media 1.02 0.980
q20_public_info 1.01 0.991
q34_02_face_mask 1.03 0.971
q34_07_hand_washing 1.02 0.977
q36_econ_worry 1.01 0.986
q40_smoking 1.03 0.972
q42_sport 1.02 0.981
q38_alcohol 1.02 0.979
q35_01_contact_close_family 1.02 0.978
q35_03_contact_friends 1.02 0.983
────────────────────────────────────────────────────
MODEL 4
Omnibus ANOVA Test
──────────────────────────────────────────────────────────────────────────────────────────────
Sum of Squares df Mean Square F p
──────────────────────────────────────────────────────────────────────────────────────────────
q01_gender 250.5010 2 125.2505 73.6681 < .001
q02_age 683.5154 1 683.5154 402.0204 < .001
q03_relationship_type 36.5492 4 9.1373 5.3742 < .001
q04_children 2.5845 1 2.5845 1.5201 0.218
q18_02_soc_media 36.6038 1 36.6038 21.5291 < .001
q20_public_info 76.5542 2 38.2771 22.5133 < .001
q34_02_face_mask 0.0257 1 0.0257 0.0151 0.902
q34_07_hand_washing 10.7997 1 10.7997 6.3520 0.012
q36_econ_worry 366.1595 2 183.0797 107.6812 < .001
q40_smoking 2.6140 1 2.6140 1.5375 0.215
q42_sport 1.4776 1 1.4776 0.8690 0.351
q38_alcohol 13.3845 1 13.3845 7.8723 0.005
q35_01_contact_close_family 80.7180 2 40.3590 23.7378 < .001
q35_03_contact_friends 20.4692 2 10.2346 6.0196 0.002
q47_self_reporting_health 1144.1397 4 286.0349 168.2360 < .001
q48_chronic_illness 1.1393 1 1.1393 0.6701 0.413
q49_health_limitations 66.7325 2 33.3663 19.6249 < .001
Residuals 12941.9291 7612 1.7002
──────────────────────────────────────────────────────────────────────────────────────────────
Note. Type 3 sum of squares
Model Coefficients - PHQ8_t
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Predictor Estimate SE Lower Upper t p Stand. Estimate Lower Upper
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Intercept ᵃ 6.36415 0.24142 5.8909 6.83739 26.362 < .001
q01_gender:
male – female -0.45310 0.03741 -0.5264 -0.37976 -12.111 < .001 -0.30382 -0.3530 -0.2546
other – female 0.42854 0.58586 -0.7199 1.57699 0.731 0.465 0.28735 -0.4827 1.0574
q02_age -0.02735 0.00136 -0.0300 -0.02468 -20.050 < .001 -0.24453 -0.2684 -0.2206
q03_relationship_type:
relationship – single -0.14746 0.04619 -0.2380 -0.05690 -3.192 0.001 -0.09888 -0.1596 -0.0382
married – single -0.06930 0.04894 -0.1652 0.02663 -1.416 0.157 -0.04647 -0.1108 0.0179
divorced – single 0.11024 0.07023 -0.0274 0.24792 1.570 0.117 0.07392 -0.0184 0.1662
widowed – single -0.02409 0.13651 -0.2917 0.24352 -0.176 0.860 -0.01615 -0.1956 0.1633
q04_children:
yes – no -0.04687 0.03801 -0.1214 0.02765 -1.233 0.218 -0.03143 -0.0814 0.0185
q18_02_soc_media:
no – yes -0.14764 0.03182 -0.2100 -0.08527 -4.640 < .001 -0.09900 -0.1408 -0.0572
q20_public_info:
yes – no -0.21014 0.03310 -0.2750 -0.14526 -6.349 < .001 -0.14091 -0.1844 -0.0974
do_not_know – no -0.02047 0.05491 -0.1281 0.08718 -0.373 0.709 -0.01372 -0.0859 0.0585
q34_02_face_mask:
yes – no 0.00407 0.03310 -0.0608 0.06895 0.123 0.902 0.00273 -0.0408 0.0462
q34_07_hand_washing:
yes – no 0.22636 0.08981 0.0503 0.40242 2.520 0.012 0.15178 0.0337 0.2698
q36_econ_worry:
serious – very_serious -0.28713 0.03989 -0.3653 -0.20893 -7.198 < .001 -0.19253 -0.2450 -0.1401
limited – very_serious -0.59292 0.04108 -0.6734 -0.51239 -14.434 < .001 -0.39758 -0.4516 -0.3436
q40_smoking:
no – yes -0.04532 0.03655 -0.1170 0.02633 -1.240 0.215 -0.03039 -0.0784 0.0177
q42_sport:
yes – no 0.02873 0.03081 -0.0317 0.08913 0.932 0.351 0.01926 -0.0212 0.0598
q38_alcohol:
no – yes -0.09008 0.03211 -0.1530 -0.02714 -2.806 0.005 -0.06040 -0.1026 -0.0182
q35_01_contact_close_family:
as_before – less_often -0.18729 0.03281 -0.2516 -0.12298 -5.709 < .001 -0.12559 -0.1687 -0.0825
more_often – less_often 0.07432 0.04948 -0.0227 0.17132 1.502 0.133 0.04984 -0.0152 0.1149
q35_03_contact_friends:
as_before – less_often -0.25137 0.07289 -0.3943 -0.10850 -3.449 < .001 -0.16856 -0.2644 -0.0728
more_often – less_often 0.04463 0.16879 -0.2862 0.37550 0.264 0.791 0.02993 -0.1919 0.2518
q47_self_reporting_health:
excellent – very_bad -2.14982 0.22440 -2.5897 -1.70994 -9.580 < .001 -1.44154 -1.7365 -1.1466
good – very_bad -1.60660 0.22207 -2.0419 -1.17129 -7.235 < .001 -1.07730 -1.3692 -0.7854
neutral – very_bad -1.00010 0.22224 -1.4358 -0.56444 -4.500 < .001 -0.67061 -0.9627 -0.3785
bad – very_bad -0.44890 0.22843 -0.8967 -0.00112 -1.965 0.049 -0.30101 -0.6013 -7.53e-4
q48_chronic_illness:
no – yes 0.03105 0.03793 -0.0433 0.10539 0.819 0.413 0.02082 -0.0290 0.0707
q49_health_limitations:
partially_limits – limits -0.27710 0.07033 -0.4150 -0.13923 -3.940 < .001 -0.18581 -0.2783 -0.0934
no_limits – limits -0.41416 0.06857 -0.5486 -0.27975 -6.040 < .001 -0.27771 -0.3678 -0.1876
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
ᵃ Represents reference level
ASSUMPTION CHECKS
Durbin–Watson Test for Autocorrelation
────────────────────────────────────────────
Autocorrelation DW Statistic p
────────────────────────────────────────────
0.0261 1.95 0.020
────────────────────────────────────────────
Collinearity Statistics
────────────────────────────────────────────────────
VIF Tolerance
────────────────────────────────────────────────────
q01_gender 1.02 0.982
q02_age 1.22 0.820
q03_relationship_type 1.08 0.924
q04_children 1.25 0.799
q18_02_soc_media 1.02 0.979
q20_public_info 1.01 0.990
q34_02_face_mask 1.03 0.966
q34_07_hand_washing 1.02 0.976
q36_econ_worry 1.02 0.983
q40_smoking 1.03 0.970
q42_sport 1.03 0.969
q38_alcohol 1.02 0.976
q35_01_contact_close_family 1.02 0.978
q35_03_contact_friends 1.02 0.982
q47_self_reporting_health 1.05 0.951
q48_chronic_illness 1.19 0.838
q49_health_limitations 1.12 0.897
────────────────────────────────────────────────────
linreg_stepwise2 <- jmv::linReg(
data = data,
dep = "PHQ8_t",
covs = "q02_age",
factors = vars("q01_gender",
"q03_relationship_type",
"q04_children",
"q18_02_soc_media",
"q20_public_info",
"q34_02_face_mask",
"q36_econ_worry",
"q38_alcohol",
"q40_smoking",
"q47_self_reporting_health",
"q48_chronic_illness",
"q49_health_limitations"),
blocks = list(
list(
"q01_gender",
"q02_age",
"q04_children",
"q36_econ_worry",
"q18_02_soc_media",
"q47_self_reporting_health",
"q49_health_limitations"),
list(
"q03_relationship_type",
"q20_public_info",
"q34_02_face_mask",
"q38_alcohol",
"q40_smoking",
"q48_chronic_illness")),
refLevels = list(
list(
var = "q01_gender",
ref = "female"),
list(
var = "q04_children",
ref = "no"),
list(
var = "q20_public_info",
ref = "no"),
list(
var = "q34_02_face_mask",
ref = "no"),
list(
var = "q36_econ_worry",
ref = "very_serious"),
list(
var = "q40_smoking",
ref = "yes"),
list(
var = "q38_alcohol",
ref = "yes"),
list(
var = "q18_02_soc_media",
ref = "yes"),
list(
var = "q03_relationship_type",
ref = "single"),
list(
var = "q47_self_reporting_health",
ref = "very_bad"),
list(
var = "q49_health_limitations",
ref = "limits"),
list(
var = "q48_chronic_illness",
ref = "yes")),
r2Adj = TRUE,
aic = TRUE,
bic = TRUE,
rmse = TRUE,
modelTest = TRUE,
anova = TRUE,
ci = TRUE,
stdEst = TRUE,
ciStdEst = TRUE,
durbin = TRUE,
collin = TRUE)
linreg_stepwise2$modelFit
Model Fit Measures
────────────────────────────────────────────────────────────────────────────────────────────────────
Model R R² Adjusted R² AIC BIC RMSE F df1 df2 p
────────────────────────────────────────────────────────────────────────────────────────────────────
1 0.472 0.223 0.222 26085 26190 1.31 170 13 7684 < .001
2 0.481 0.231 0.229 26027 26200 1.31 100 23 7674 < .001
────────────────────────────────────────────────────────────────────────────────────────────────────
linreg_stepwise2$modelComp
Model Comparisons
───────────────────────────────────────────────────────────────────
Model Model ΔR² F df1 df2 p
───────────────────────────────────────────────────────────────────
1 - 2 0.00791 7.90 10 7674 < .001
───────────────────────────────────────────────────────────────────
linreg_stepwise2$models
MODEL SPECIFIC RESULTS
MODEL 1
Omnibus ANOVA Test
──────────────────────────────────────────────────────────────────────────────────────────
Sum of Squares df Mean Square F p
──────────────────────────────────────────────────────────────────────────────────────────
q01_gender 282.57 2 141.29 81.63 < .001
q02_age 838.16 1 838.16 484.24 < .001
q04_children 8.79 1 8.79 5.08 0.024
q36_econ_worry 451.29 2 225.65 130.36 < .001
q18_02_soc_media 38.98 1 38.98 22.52 < .001
q47_self_reporting_health 1273.34 4 318.34 183.91 < .001
q49_health_limitations 69.10 2 34.55 19.96 < .001
Residuals 13300.19 7684 1.73
──────────────────────────────────────────────────────────────────────────────────────────
Note. Type 3 sum of squares
Model Coefficients - PHQ8_t
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Predictor Estimate SE Lower Upper t p Stand. Estimate Lower Upper
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Intercept ᵃ 6.3406 0.22475 5.9000 6.78116 28.212 < .001
q01_gender:
male – female -0.4706 0.03696 -0.5430 -0.39813 -12.733 < .001 -0.3155 -0.3641 -0.26696
other – female 0.5513 0.59053 -0.6063 1.70894 0.934 0.351 0.3697 -0.4065 1.14591
q02_age -0.0268 0.00122 -0.0292 -0.02445 -22.005 < .001 -0.2404 -0.2618 -0.21896
q04_children:
yes – no -0.0738 0.03276 -0.1381 -0.00962 -2.254 0.024 -0.0495 -0.0926 -0.00645
q36_econ_worry:
serious – very_serious -0.3098 0.03985 -0.3879 -0.23169 -7.774 < .001 -0.2077 -0.2601 -0.15536
limited – very_serious -0.6423 0.04050 -0.7217 -0.56289 -15.860 < .001 -0.4307 -0.4839 -0.37744
q18_02_soc_media:
no – yes -0.1510 0.03182 -0.2134 -0.08863 -4.745 < .001 -0.1013 -0.1431 -0.05943
q47_self_reporting_health:
excellent – very_bad -2.1832 0.22508 -2.6244 -1.74199 -9.700 < .001 -1.4639 -1.7598 -1.16807
good – very_bad -1.6334 0.22307 -2.0707 -1.19617 -7.323 < .001 -1.0953 -1.3885 -0.80208
neutral – very_bad -1.0295 0.22358 -1.4678 -0.59121 -4.605 < .001 -0.6903 -0.9842 -0.39643
bad – very_bad -0.4859 0.22983 -0.9364 -0.03538 -2.114 0.035 -0.3258 -0.6279 -0.02372
q49_health_limitations:
partially_limits – limits -0.2810 0.07025 -0.4187 -0.14330 -4.000 < .001 -0.1884 -0.2808 -0.09609
no_limits – limits -0.4063 0.06722 -0.5380 -0.27449 -6.044 < .001 -0.2724 -0.3608 -0.18406
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
ᵃ Represents reference level
ASSUMPTION CHECKS
Durbin–Watson Test for Autocorrelation
────────────────────────────────────────────
Autocorrelation DW Statistic p
────────────────────────────────────────────
0.0289 1.94 0.004
────────────────────────────────────────────
Collinearity Statistics
──────────────────────────────────────────────────
VIF Tolerance
──────────────────────────────────────────────────
q01_gender 1.01 0.991
q02_age 1.09 0.921
q04_children 1.07 0.932
q36_econ_worry 1.01 0.994
q18_02_soc_media 1.02 0.984
q47_self_reporting_health 1.04 0.965
q49_health_limitations 1.07 0.932
──────────────────────────────────────────────────
MODEL 2
Omnibus ANOVA Test
───────────────────────────────────────────────────────────────────────────────────────────
Sum of Squares df Mean Square F p
───────────────────────────────────────────────────────────────────────────────────────────
q01_gender 284.252 2 142.126 82.848 < .001
q02_age 686.861 1 686.861 400.386 < .001
q04_children 4.143 1 4.143 2.415 0.120
q36_econ_worry 399.464 2 199.732 116.428 < .001
q18_02_soc_media 31.564 1 31.564 18.400 < .001
q47_self_reporting_health 1181.145 4 295.286 172.129 < .001
q49_health_limitations 66.702 2 33.351 19.441 < .001
q03_relationship_type 42.350 4 10.587 6.172 < .001
q20_public_info 77.084 2 38.542 22.467 < .001
q34_02_face_mask 1.246 1 1.246 0.726 0.394
q38_alcohol 17.978 1 17.978 10.480 0.001
q40_smoking 0.704 1 0.704 0.410 0.522
q48_chronic_illness 0.923 1 0.923 0.538 0.463
Residuals 13164.721 7674 1.715
───────────────────────────────────────────────────────────────────────────────────────────
Note. Type 3 sum of squares
Model Coefficients - PHQ8_t
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Predictor Estimate SE Lower Upper t p Stand. Estimate Lower Upper
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Intercept ᵃ 6.4542 0.22777 6.00775 6.9007 28.337 < .001
q01_gender:
male – female -0.4747 0.03700 -0.54729 -0.4022 -12.830 < .001 -0.3183 -0.36698 -0.2697
other – female 0.5557 0.58825 -0.59746 1.7088 0.945 0.345 0.3726 -0.40062 1.1458
q02_age -0.0268 0.00134 -0.02941 -0.0242 -20.010 < .001 -0.2400 -0.26346 -0.2164
q04_children:
yes – no -0.0588 0.03782 -0.13292 0.0154 -1.554 0.120 -0.0394 -0.08913 0.0103
q36_econ_worry:
serious – very_serious -0.2970 0.03974 -0.37491 -0.2191 -7.475 < .001 -0.1992 -0.25139 -0.1469
limited – very_serious -0.6136 0.04085 -0.69363 -0.5335 -15.019 < .001 -0.4114 -0.46511 -0.3577
q18_02_soc_media:
no – yes -0.1365 0.03182 -0.19887 -0.0741 -4.289 < .001 -0.0915 -0.13335 -0.0497
q47_self_reporting_health:
excellent – very_bad -2.1392 0.22501 -2.58024 -1.6981 -9.507 < .001 -1.4344 -1.73015 -1.1386
good – very_bad -1.5847 0.22272 -2.02129 -1.1481 -7.115 < .001 -1.0626 -1.35536 -0.7699
neutral – very_bad -0.9846 0.22292 -1.42160 -0.5476 -4.417 < .001 -0.6602 -0.95324 -0.3672
bad – very_bad -0.4319 0.22906 -0.88091 0.0171 -1.885 0.059 -0.2896 -0.59069 0.0115
q49_health_limitations:
partially_limits – limits -0.2769 0.07011 -0.41431 -0.1394 -3.949 < .001 -0.1857 -0.27781 -0.0935
no_limits – limits -0.4119 0.06836 -0.54592 -0.2779 -6.025 < .001 -0.2762 -0.36606 -0.1863
q03_relationship_type:
relationship – single -0.1507 0.04606 -0.24096 -0.0604 -3.271 0.001 -0.1010 -0.16157 -0.0405
married – single -0.0633 0.04878 -0.15888 0.0324 -1.297 0.195 -0.0424 -0.10654 0.0217
divorced – single 0.1300 0.07002 -0.00727 0.2673 1.856 0.063 0.0872 -0.00488 0.1792
widowed – single 0.0150 0.13537 -0.25034 0.2804 0.111 0.912 0.0101 -0.16786 0.1880
q20_public_info:
yes – no -0.2100 0.03310 -0.27487 -0.1451 -6.344 < .001 -0.1408 -0.18431 -0.0973
do_not_know – no -0.0204 0.05500 -0.12818 0.0875 -0.370 0.711 -0.0137 -0.08595 0.0586
q34_02_face_mask:
yes – no 0.0280 0.03287 -0.03642 0.0924 0.852 0.394 0.0188 -0.02442 0.0620
q38_alcohol:
no – yes -0.1035 0.03196 -0.16613 -0.0408 -3.237 0.001 -0.0694 -0.11140 -0.0274
q40_smoking:
no – yes -0.0232 0.03627 -0.09434 0.0479 -0.641 0.522 -0.0156 -0.06326 0.0321
q48_chronic_illness:
no – yes 0.0278 0.03793 -0.04653 0.1022 0.734 0.463 0.0187 -0.03120 0.0685
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
ᵃ Represents reference level
ASSUMPTION CHECKS
Durbin–Watson Test for Autocorrelation
────────────────────────────────────────────
Autocorrelation DW Statistic p
────────────────────────────────────────────
0.0288 1.94 0.014
────────────────────────────────────────────
Collinearity Statistics
──────────────────────────────────────────────────
VIF Tolerance
──────────────────────────────────────────────────
q01_gender 1.01 0.988
q02_age 1.20 0.835
q04_children 1.24 0.804
q36_econ_worry 1.01 0.986
q18_02_soc_media 1.02 0.980
q47_self_reporting_health 1.05 0.955
q49_health_limitations 1.11 0.898
q03_relationship_type 1.08 0.925
q20_public_info 1.01 0.990
q34_02_face_mask 1.03 0.975
q38_alcohol 1.02 0.980
q40_smoking 1.02 0.979
q48_chronic_illness 1.19 0.839
──────────────────────────────────────────────────
Secondly, we conduct a Reliability Analysis of the Covid-19 concern factor. We use a cutoff value of 0.7 for both McDonald’s Omega and Cronbach’s Alpha. The scale passes this cutoff and the statistics would not be improved if any of the items were dropped.
# To conduct the reliability analysis, we use the reliability() function from the
# jmv package on the "test" data set (as opposed to the "train" dataset used for EFA).
jmv::reliability(
data = data,
vars = vars("q30_concern_infection_covid",
"q31_concern_infection_friends",
"q33_01_concern_situation",
"q33_02_concern_low_control",
"q33_03_concern_survival_covid",
"q33_05_concern_infecting_others"),
omegaScale = TRUE,
alphaItems = TRUE,
omegaItems = TRUE)
RELIABILITY ANALYSIS
Scale Reliability Statistics
─────────────────────────────────────────
Cronbach's α McDonald's ω
─────────────────────────────────────────
scale 0.801 0.815
─────────────────────────────────────────
Item Reliability Statistics
───────────────────────────────────────────────────────────────────
Cronbach's α McDonald's ω
───────────────────────────────────────────────────────────────────
q30_concern_infection_covid 0.738 0.750
q31_concern_infection_friends 0.751 0.768
q33_01_concern_situation 0.745 0.766
q33_02_concern_low_control 0.820 0.833
q33_03_concern_survival_covid 0.763 0.784
q33_05_concern_infecting_others 0.794 0.811
───────────────────────────────────────────────────────────────────
# To conduct the CFA, we use the cfa() function from the jmv package on the "test"
# data set (as opposed to the "train" dataset used for EFA).
jmv::cfa(
data = data,
factors = list(
list(
label = "Concern",
vars = c("q30_concern_infection_covid",
"q31_concern_infection_friends",
"q33_01_concern_situation",
"q33_02_concern_low_control",
"q33_03_concern_survival_covid",
"q33_05_concern_infecting_others"))),
resCov = list(),
ci = TRUE,
stdEst = TRUE,
factCovEst = FALSE,
fitMeasures = c("cfi", "tli", "rmsea", "srmr"),
corRes = TRUE)
CONFIRMATORY FACTOR ANALYSIS
Factor Loadings
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Factor Indicator Estimate SE Lower Upper Z p Stand. Estimate
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Concern q30_concern_infection_covid 2.270 0.0255 2.220 2.321 88.9 < .001 0.862
q31_concern_infection_friends 1.862 0.0256 1.812 1.912 72.6 < .001 0.755
q33_01_concern_situation 2.059 0.0282 2.004 2.114 72.9 < .001 0.746
q33_02_concern_low_control 0.855 0.0317 0.793 0.917 26.9 < .001 0.319
q33_03_concern_survival_covid 1.847 0.0291 1.790 1.904 63.4 < .001 0.672
q33_05_concern_infecting_others 1.419 0.0352 1.350 1.488 40.3 < .001 0.467
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
MODEL FIT
Test for Exact Fit
────────────────────────
χ² df p
────────────────────────
1194 9 < .001
────────────────────────
Fit Measures
───────────────────────────────────────────────────────
CFI TLI SRMR RMSEA Lower Upper
───────────────────────────────────────────────────────
0.926 0.877 0.0432 0.129 0.122 0.135
───────────────────────────────────────────────────────
POST-HOC MODEL PERFORMANCE
Residuals for Observed Correlation Matrix
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
q30_concern_infection_covid q31_concern_infection_friends q33_01_concern_situation q33_02_concern_low_control q33_03_concern_survival_covid q33_05_concern_infecting_others
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
q30_concern_infection_covid 0.017 -0.007 -0.055 0.036 -0.073
q31_concern_infection_friends -0.007 -0.041 -0.078 0.103
q33_01_concern_situation 0.058 -0.003 0.027
q33_02_concern_low_control 0.064 0.116
q33_03_concern_survival_covid -0.012
q33_05_concern_infecting_others
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
# Creating the Covid-19-related concern/anxiety index, consisting of the average of
# the values of the multiple variables selected through factor analysis to
# represent the underlying construct.
concern_index <- apply(cbind(data$q30_concern_infection_covid,
data$q31_concern_infection_friends,
data$q33_01_concern_situation,
data$q33_02_concern_low_control,
data$q33_03_concern_survival_covid,
data$q33_05_concern_infecting_others), 1, mean)
#Adding the vector as an column to the existing dataset.
data <- cbind(data, concern_index)
data_corr <- cbind(data_corr, concern_index)
#To summarize the concern_index variable, we use the descriptives()
# function from the jmv package.
anx_index_descriptives <- jmv::descriptives(data = data,
missing = TRUE,
vars = "concern_index",
sd = TRUE,
median = FALSE,
pc = TRUE,
range = TRUE,
box = TRUE)
# Function to get the result from the correlation matrix into a data frame
flattenCorrMatrix <- function(cormat, pmat) {
ut <- upper.tri(cormat)
data.frame(
row = rownames(cormat)[row(cormat)[ut]],
column = rownames(cormat)[col(cormat)[ut]],
cor = (cormat)[ut],
p = pmat[ut]
)
}
#Correlation matrix using Spearman coefficient
corr_mtx <- rcorr(as.matrix(data_corr), type = "spearman")
# Selecting only significant correlates for PHQ8 (values with p>0.05 are excluded)
flattenCorrMatrix(corr_mtx$r, corr_mtx$P) %>% filter(p <= 0.05,
column %in% c("PHQ8_t")) %>%
arrange(desc(abs(cor)))
row column cor p
1 q47_self_reporting_health PHQ8_t 0.3090 0.0000000000000
2 q02_age PHQ8_t -0.2162 0.0000000000000
3 q36_econ_worry PHQ8_t -0.2121 0.0000000000000
4 q49_health_limitations PHQ8_t -0.1651 0.0000000000000
5 q01_gender PHQ8_t 0.1460 0.0000000000000
6 q20_public_info PHQ8_t 0.1340 0.0000000000000
7 q18_02_soc_media PHQ8_t 0.1123 0.0000000000000
8 q03_relationship_type PHQ8_t -0.0998 0.0000000000000
9 q48_chronic_illness PHQ8_t -0.0949 0.0000000000000
10 q04_children PHQ8_t 0.0783 0.0000000000044
11 q42_sport PHQ8_t 0.0630 0.0000000256406
12 q11_education PHQ8_t -0.0546 0.0000383164221
13 q40_smoking PHQ8_t -0.0536 0.0000021447198
14 q35_03_contact_friends PHQ8_t -0.0392 0.0005412876720
15 q34_02_face_mask PHQ8_t -0.0337 0.0029777622136
# Selecting only significant correlates for concern index (values with p>0.05 are excluded)
flattenCorrMatrix(corr_mtx$r, corr_mtx$P) %>% filter(p <= 0.05,
column %in% c("concern_index")) %>%
arrange(desc(abs(cor)))
row column cor p
1 PHQ8_t concern_index 0.3296 0.00000000000000000
2 q47_self_reporting_health concern_index 0.2600 0.00000000000000000
3 q01_gender concern_index 0.1922 0.00000000000000000
4 q49_health_limitations concern_index -0.1483 0.00000000000000000
5 q34_07_hand_washing concern_index -0.1483 0.00000000000000000
6 q48_chronic_illness concern_index -0.1410 0.00000000000000000
7 q20_public_info concern_index 0.1237 0.00000000000000000
8 q35_03_contact_friends concern_index -0.1215 0.00000000000000000
9 q36_econ_worry concern_index -0.1131 0.00000000000000000
10 q34_02_face_mask concern_index -0.0911 0.00000000000000488
11 q42_sport concern_index 0.0770 0.00000000003795764
12 q11_education concern_index 0.0684 0.00000057803401488
13 q38_alcohol concern_index 0.0548 0.00000255020021722
14 q04_children concern_index -0.0520 0.00000783470995414
15 q18_02_soc_media concern_index 0.0454 0.00009800170357166
16 q35_01_contact_close_family concern_index -0.0453 0.00010456661621516
17 q03_relationship_type concern_index 0.0296 0.01110354779772171
data %>%
filter(concern_index != "NA") %>%
ggplot(aes(y = concern_index)) +
geom_boxplot(fill = "#4DBBD5CC") +
theme_bw() +
scale_y_continuous(breaks = seq(0, 10, 1),
labels = seq(0, 10, 1),
limits = c(0, 10)) +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank()) +
ylab("Concern Index")
data %>%
filter(concern_index != "NA") %>%
ggplot(aes(x = concern_index)) +
geom_density(fill = "#4DBBD5CC") +
scale_x_continuous(breaks = seq(0, 10, 1),
labels = seq(0, 10, 1),
limits = c(0, 10)) +
theme_bw() +
theme(legend.position = "none") +
xlab("Concern Index")
anx_index_descriptives$descriptives
Descriptives
───────────────────────────────────────
concern_index
───────────────────────────────────────
N 7371
Missing 620
Mean 5.47
Standard deviation 1.92
Range 9.00
Minimum 1.00
Maximum 10.0
25th percentile 4.17
50th percentile 5.50
75th percentile 6.83
───────────────────────────────────────
# Before running the model, we need to transform the social media string
# dummy (yes/no) back to its numeric form, with similar operation for gender.
levels(data$q18_02_soc_media) <- list("1" = "yes", "0" = "no")
levels(data$q01_gender) <- list("0" = "male", "1" = "female", "2" = "other")
data$q01_gender <- as.numeric(as.character(data$q01_gender))
data$q18_02_soc_media <- as.numeric(as.character(data$q18_02_soc_media))
# Centering continuous variables with scaling
data_sem <- data %>%
filter(q01_gender != 2, !is.na(concern_index)) %>%
mutate(concern_index.c = scale(concern_index, scale = TRUE),
PHQ8.c = scale(PHQ8_t, scale = TRUE),
q02_age.c = scale(q02_age, scale = TRUE))
# Labels for diagrams
labels_H76 <- list(X = "Social Media",
M = "Concern",
Y = "Depression",
W = "Age",
Z = "Gender")
pmacroModel(no = 76,
labels = labels_H76,
xmargin = 0,
rady = 0.047,
radx = 0.09,
ylim = c(0.15, 0.8))
statisticalDiagram(no = 76,
labels = labels_H76,
whatLabel = "name",
xmargin = 0.01,
rady = 0.03,
radx = 0.11,
ylim = c(0.06, 0.95),
xlim = c(0.01, 1))
In the second step, we specify the key pathways and run the analysis, while bootstrapping the confidence intervals.
# Mediation-moderation analysis (path analysis framework, SEM) using lavaan package.
# First, we specify the model pathways
spec_mod <- "
# Regressions
concern_index.c ~ a1*q18_02_soc_media + a2*q02_age.c + a3*q01_gender + a4*q18_02_soc_media:q02_age.c + a5*q18_02_soc_media:q01_gender
PHQ8.c ~ c1*q18_02_soc_media + c2*q02_age.c + c3*q01_gender + c4*q18_02_soc_media:q02_age.c + c5*q18_02_soc_media:q01_gender + b1*concern_index.c + b2*concern_index.c:q02_age.c + b3*concern_index.c:q01_gender
#Mean and variance of age and gender moderators
q02_age.c ~ q02_age.c.mean*1
q02_age.c ~~ q02_age.c.var*q02_age.c
q01_gender ~ q01_gender.mean*1
q01_gender ~~ q01_gender.var*q01_gender
# Effect specifications
XonM := a1 + a4*q02_age.c.mean + a5*q01_gender.mean
MonY := b1 + b2*q02_age.c.mean + b3*q01_gender.mean
indirect := (a1 + a4*q02_age.c.mean + a5*q01_gender.mean)*(b1 + b2*q02_age.c.mean + b3*q01_gender.mean)
direct := c1 + c4*q02_age.c.mean + c5*q01_gender.mean
total := direct + indirect
prop.mediated := indirect / total
# Component effects conditional on moderators (X = Social Media, M = Concern, Y = Depression, W = Age, Z = Gender)
XonM.mean.male := a1 + a4*q02_age.c.mean + a5*0
XonM.mean.female := a1 + a4*q02_age.c.mean + a5*1
XonM.blw.male := a1 + a4*(q02_age.c.mean - sqrt(q02_age.c.var)) + a5*0
XonM.blw.female := a1 + a4*(q02_age.c.mean - sqrt(q02_age.c.var)) + a5*1
XonM.blw.avg := a1 + a4*(q02_age.c.mean - sqrt(q02_age.c.var)) + a5*q01_gender.mean
XonM.abv.male := a1 + a4*(q02_age.c.mean + sqrt(q02_age.c.var)) + a5*0
XonM.abv.female := a1 + a4*(q02_age.c.mean + sqrt(q02_age.c.var)) + a5*1
XonM.abv.avg := a1 + a4*(q02_age.c.mean + sqrt(q02_age.c.var)) + a5*q01_gender.mean
MonY.mean.male := b1 + b2*q02_age.c.mean + b3*0
MonY.mean.female := b1 + b2*q02_age.c.mean + b3*1
MonY.blw.male := b1 + b2*(q02_age.c.mean - sqrt(q02_age.c.var)) + b3*0
MonY.blw.female := b1 + b2*(q02_age.c.mean - sqrt(q02_age.c.var)) + b3*1
MonY.blw.avg := b1 + b2*(q02_age.c.mean - sqrt(q02_age.c.var)) + b3*q01_gender.mean
MonY.abv.male := b1 + b2*(q02_age.c.mean + sqrt(q02_age.c.var)) + b3*0
MonY.abv.female := b1 + b2*(q02_age.c.mean + sqrt(q02_age.c.var)) + b3*1
MonY.abv.avg := b1 + b2*(q02_age.c.mean + sqrt(q02_age.c.var)) + b3*q01_gender.mean
# Indirect effects conditional on moderators
indirect.mean.male := (a1 + a4*q02_age.c.mean + a5*0)*(b1 + b2*q02_age.c.mean + b3*0)
indirect.mean.female := (a1 + a4*q02_age.c.mean + a5*1)*(b1 + b2*q02_age.c.mean + b3*1)
indirect.blw.male := (a1 + a4*(q02_age.c.mean - sqrt(q02_age.c.var)) + a5*0)*(b1 + b2*(q02_age.c.mean - sqrt(q02_age.c.var)) + b3*0)
indirect.blw.female := (a1 + a4*(q02_age.c.mean - sqrt(q02_age.c.var)) + a5*1)*(b1 + b2*(q02_age.c.mean - sqrt(q02_age.c.var)) + b3*1)
indirect.blw.avg := (a1 + a4*(q02_age.c.mean - sqrt(q02_age.c.var)) + a5*q01_gender.mean)*(b1 + b2*(q02_age.c.mean - sqrt(q02_age.c.var)) + b3*q01_gender.mean)
indirect.abv.male := (a1 + a4*(q02_age.c.mean + sqrt(q02_age.c.var)) + a5*0)*(b1 + b2*(q02_age.c.mean + sqrt(q02_age.c.var)) + b3*0)
indirect.abv.female := (a1 + a4*(q02_age.c.mean + sqrt(q02_age.c.var)) + a5*1)*(b1 + b2*(q02_age.c.mean + sqrt(q02_age.c.var)) + b3*1)
indirect.abv.avg := (a1 + a4*(q02_age.c.mean + sqrt(q02_age.c.var)) + a5*q01_gender.mean)*(b1 + b2*(q02_age.c.mean + sqrt(q02_age.c.var)) + b3*q01_gender.mean)
# Direct effects conditional on moderators
direct.mean.male := c1 + c4*q02_age.c.mean + c5*0
direct.mean.female := c1 + c4*q02_age.c.mean + c5*1
direct.blw.male := c1 + c4*(q02_age.c.mean - sqrt(q02_age.c.var)) + c5*0
direct.blw.female := c1 + c4*(q02_age.c.mean - sqrt(q02_age.c.var)) + c5*1
direct.blw.avg := c1 + c4*(q02_age.c.mean - sqrt(q02_age.c.var)) + c5*q01_gender.mean
direct.abv.male := c1 + c4*(q02_age.c.mean + sqrt(q02_age.c.var)) + c5*0
direct.abv.female := c1 + c4*(q02_age.c.mean + sqrt(q02_age.c.var)) + c5*1
direct.abv.avg := c1 + c4*(q02_age.c.mean + sqrt(q02_age.c.var)) + c5*q01_gender.mean
# Total effects conditional on moderators
total.mean.male := direct.mean.male + indirect.mean.male
total.mean.female := direct.mean.female + indirect.mean.female
total.blw.male := direct.blw.male + indirect.blw.male
total.blw.female := direct.blw.female + indirect.blw.female
total.blw.avg := direct.blw.avg + indirect.blw.avg
total.abv.male := direct.abv.male + indirect.abv.male
total.abv.female := direct.abv.female + indirect.abv.female
total.abv.avg := direct.abv.avg + indirect.abv.avg
# Proportion mediated conditional on moderators
prop.med.mean.male := indirect.mean.male / total.mean.male
prop.med.mean.female := indirect.mean.female / total.mean.female
prop.med.blw.male := indirect.blw.male / total.blw.male
prop.med.blw.female := indirect.blw.female / total.blw.female
prop.med.blw.avg := indirect.blw.avg / total.blw.avg
prop.med.abv.male := indirect.abv.male / total.abv.male
prop.med.abv.female := indirect.abv.female / total.abv.male
prop.med.abv.avg := indirect.abv.avg / total.abv.avg"
# For reproducibility of results (using bootstrap)
set.seed(2021)
# Secondly, we fit/estimate the model and we use bootstrap for robustness.
fit_mod <- lavaan::sem(model = spec_mod,
data = data_sem,
se = "bootstrap",
bootstrap = 1000)
# Labels for statistical diagrams
labels_stats_H76 <- list(X = "q18_02_soc_media",
M = "concern_index.c",
Y = "PHQ8.c",
W = "q02_age.c",
Z = "q01_gender")
statisticalDiagram(no = 76,
labels = labels_stats_H76,
fit = fit_mod,
whatLabel = "est",
xmargin = 0.01,
rady = 0.03,
radx = 0.158,
ylim = c(0.06, 0.95),
xlim = c(0.01, 1))
statisticalDiagram(no = 76,
labels = labels_stats_H76,
fit = fit_mod,
whatLabel = "std",
xmargin = 0.01,
rady = 0.03,
radx = 0.158,
ylim = c(0.06, 0.95),
xlim = c(0.01, 1))
lavaan::summary(fit_mod,
rsquare = TRUE,
ci = TRUE,
fit.measures = TRUE,
standardize = TRUE)
lavaan 0.6-9 ended normally after 44 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 21
Used Total
Number of observations 7239 7365
Model Test User Model:
Test statistic 18410.525
Degrees of freedom 13
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 20100.115
Degrees of freedom 26
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.084
Tucker-Lewis Index (TLI) -0.833
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -33838.618
Loglikelihood unrestricted model (H1) -24633.356
Akaike (AIC) 67719.236
Bayesian (BIC) 67863.868
Sample-size adjusted Bayesian (BIC) 67797.135
Root Mean Square Error of Approximation:
RMSEA 0.442
90 Percent confidence interval - lower 0.437
90 Percent confidence interval - upper 0.448
P-value RMSEA <= 0.05 0.000
Standardized Root Mean Square Residual:
SRMR 0.162
Parameter Estimates:
Standard errors Bootstrap
Number of requested bootstrap draws 1000
Number of successful bootstrap draws 1000
Regressions:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
concern_index.c ~
q18_02_s_ (a1) -0.027 0.050 -0.552 0.581 -0.128 0.070
q02_age.c (a2) 0.010 0.014 0.708 0.479 -0.017 0.038
q01_gendr (a3) 0.397 0.032 12.448 0.000 0.332 0.457
q18_02__: (a4) 0.023 0.026 0.879 0.379 -0.030 0.072
q18_02__: (a5) 0.160 0.056 2.865 0.004 0.051 0.269
PHQ8.c ~
q18_02_s_ (c1) 0.145 0.050 2.883 0.004 0.046 0.240
q02_age.c (c2) -0.204 0.014 -14.813 0.000 -0.230 -0.176
q01_gendr (c3) 0.201 0.033 6.152 0.000 0.136 0.269
q18_02__: (c4) -0.032 0.024 -1.341 0.180 -0.079 0.018
q18_02__: (c5) -0.024 0.056 -0.424 0.671 -0.135 0.090
cncrn_nd. (b1) 0.328 0.026 12.525 0.000 0.275 0.377
cn_.:02_. (b2) 0.026 0.012 2.226 0.026 0.003 0.047
cnc_.:01_ (b3) -0.009 0.029 -0.326 0.745 -0.065 0.045
Std.lv Std.all
-0.027 -0.013
0.010 0.010
0.397 0.166
0.023 0.014
0.160 0.072
0.145 0.070
-0.204 -0.204
0.201 0.083
-0.032 -0.019
-0.024 -0.011
0.328 0.327
0.026 0.026
-0.009 -0.008
Intercepts:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
q02_g.c (q02_) -0.001 0.012 -0.049 0.961 -0.023 0.023
q01_gnd (q01_) 0.780 0.005 156.255 0.000 0.769 0.789
.cncrn_. -0.341 0.028 -12.069 0.000 -0.395 -0.283
.PHQ8.c -0.204 0.029 -7.162 0.000 -0.264 -0.147
Std.lv Std.all
-0.001 -0.001
0.780 1.880
-0.341 -0.343
-0.204 -0.205
Variances:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
q02_g.c (q02_) 0.997 0.015 67.692 0.000 0.969 1.026
q01_gnd (q01_) 0.172 0.003 61.727 0.000 0.167 0.177
.cncrn_. 0.956 0.014 69.225 0.000 0.928 0.981
.PHQ8.c 0.824 0.012 66.280 0.000 0.799 0.848
Std.lv Std.all
0.997 1.000
0.172 1.000
0.956 0.969
0.824 0.830
R-Square:
Estimate
concern_indx.c 0.031
PHQ8.c 0.170
Defined Parameters:
Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
XonM 0.097 0.025 3.825 0.000 0.047 0.148
MonY 0.321 0.012 27.188 0.000 0.298 0.344
indirect 0.031 0.008 3.782 0.000 0.015 0.048
direct 0.126 0.022 5.620 0.000 0.083 0.173
total 0.157 0.024 6.608 0.000 0.110 0.205
prop.mediated 0.198 0.052 3.792 0.000 0.101 0.310
XonM.mean.male -0.027 0.050 -0.552 0.581 -0.128 0.069
XonM.mean.feml 0.132 0.029 4.634 0.000 0.076 0.187
XonM.blw.male -0.050 0.051 -0.977 0.329 -0.152 0.050
XonM.blw.femal 0.109 0.035 3.159 0.002 0.039 0.177
XonM.blw.avg 0.074 0.031 2.366 0.018 0.011 0.137
XonM.abv.male -0.005 0.060 -0.077 0.939 -0.130 0.110
XonM.abv.femal 0.155 0.042 3.678 0.000 0.071 0.234
XonM.abv.avg 0.120 0.041 2.949 0.003 0.040 0.196
MonY.mean.male 0.328 0.026 12.518 0.000 0.275 0.377
MonY.mean.feml 0.319 0.013 24.703 0.000 0.294 0.344
MonY.blw.male 0.302 0.029 10.451 0.000 0.247 0.356
MonY.blw.femal 0.293 0.018 16.526 0.000 0.258 0.326
MonY.blw.avg 0.295 0.017 17.410 0.000 0.263 0.327
MonY.abv.male 0.354 0.029 12.423 0.000 0.298 0.410
MonY.abv.femal 0.345 0.017 20.128 0.000 0.311 0.380
MonY.abv.avg 0.347 0.016 21.269 0.000 0.314 0.378
indirect.mn.ml -0.009 0.016 -0.552 0.581 -0.041 0.023
indirct.mn.fml 0.042 0.009 4.542 0.000 0.024 0.060
indirct.blw.ml -0.015 0.016 -0.969 0.333 -0.046 0.016
indrct.blw.fml 0.032 0.010 3.071 0.002 0.011 0.053
indirct.blw.vg 0.022 0.009 2.330 0.020 0.003 0.041
indirect.bv.ml -0.002 0.021 -0.077 0.939 -0.043 0.040
indirct.bv.fml 0.053 0.015 3.642 0.000 0.024 0.081
indirect.bv.vg 0.042 0.014 2.935 0.003 0.013 0.067
direct.mean.ml 0.145 0.050 2.882 0.004 0.047 0.240
direct.men.fml 0.121 0.025 4.820 0.000 0.073 0.173
direct.blw.mal 0.177 0.053 3.343 0.001 0.072 0.281
direct.blw.fml 0.153 0.033 4.643 0.000 0.087 0.213
direct.blw.avg 0.158 0.030 5.213 0.000 0.095 0.217
direct.abv.mal 0.113 0.058 1.931 0.053 0.002 0.225
direct.abv.fml 0.089 0.036 2.438 0.015 0.020 0.163
direct.abv.avg 0.094 0.035 2.670 0.008 0.027 0.165
total.mean.mal 0.136 0.054 2.535 0.011 0.035 0.241
total.mean.fml 0.163 0.027 6.132 0.000 0.111 0.216
total.blw.male 0.162 0.056 2.887 0.004 0.057 0.270
total.blw.feml 0.185 0.035 5.360 0.000 0.114 0.248
total.blw.avg 0.180 0.032 5.656 0.000 0.114 0.241
total.abv.male 0.111 0.063 1.770 0.077 -0.005 0.231
total.abv.feml 0.142 0.039 3.652 0.000 0.064 0.222
total.abv.avg 0.136 0.038 3.600 0.000 0.065 0.213
prop.med.mn.ml -0.066 1.150 -0.058 0.954 -0.640 0.156
prop.md.mn.fml 0.258 0.061 4.247 0.000 0.146 0.389
prop.md.blw.ml -0.094 0.484 -0.194 0.846 -0.505 0.091
prp.md.blw.fml 0.173 0.060 2.878 0.004 0.070 0.312
prop.md.blw.vg 0.121 0.053 2.303 0.021 0.020 0.238
prop.med.bv.ml -0.015 1.938 -0.008 0.994 -1.877 0.558
prop.md.bv.fml 0.482 3.956 0.122 0.903 -1.201 4.940
prop.med.bv.vg 0.307 0.134 2.281 0.023 0.105 0.630
Std.lv Std.all
0.097 0.123
0.321 0.312
0.031 0.038
0.126 0.049
0.157 0.088
0.198 0.437
-0.027 -0.013
0.132 0.059
-0.050 -0.027
0.109 0.046
0.074 0.109
-0.005 0.000
0.155 0.073
0.120 0.136
0.328 0.327
0.319 0.319
0.302 0.301
0.293 0.293
0.295 0.285
0.354 0.353
0.345 0.345
0.347 0.338
-0.009 -0.004
0.042 0.019
-0.015 -0.008
0.032 0.013
0.022 0.031
-0.002 0.000
0.053 0.025
0.042 0.046
0.145 0.070
0.121 0.059
0.177 0.089
0.153 0.078
0.158 0.068
0.113 0.051
0.089 0.040
0.094 0.030
0.136 0.065
0.163 0.078
0.162 0.081
0.185 0.091
0.180 0.100
0.111 0.051
0.142 0.065
0.136 0.076
-0.066 -0.066
0.258 0.243
-0.094 -0.100
0.173 0.146
0.121 0.313
-0.015 0.002
0.482 0.495
0.307 0.603
estimates <- parameterEstimates(fit_mod, standardized = TRUE) %>%
filter(op == "~") %>%
select(-c(std.nox))
p_adj <- p.adjust(estimates$pvalue, method = "holm")
estimates <- cbind(estimates, p_adj)
kbl(estimates) %>%
kable_classic(full_width = FALSE, lightable_options = c("striped")) %>%
row_spec(row = which(estimates$p_adj < 0.05), bold = TRUE)
lhs | op | rhs | label | est | se | z | pvalue | ci.lower | ci.upper | std.lv | std.all | p_adj |
---|---|---|---|---|---|---|---|---|---|---|---|---|
concern_index.c | ~ | q18_02_soc_media | a1 | -0.027 | 0.050 | -0.552 | 0.581 | -0.128 | 0.070 | -0.027 | -0.013 | 1.000 |
concern_index.c | ~ | q02_age.c | a2 | 0.010 | 0.014 | 0.708 | 0.479 | -0.017 | 0.038 | 0.010 | 0.010 | 1.000 |
concern_index.c | ~ | q01_gender | a3 | 0.397 | 0.032 | 12.448 | 0.000 | 0.332 | 0.457 | 0.397 | 0.166 | 0.000 |
concern_index.c | ~ | q18_02_soc_media:q02_age.c | a4 | 0.023 | 0.026 | 0.879 | 0.379 | -0.030 | 0.072 | 0.023 | 0.014 | 1.000 |
concern_index.c | ~ | q18_02_soc_media:q01_gender | a5 | 0.160 | 0.056 | 2.865 | 0.004 | 0.051 | 0.269 | 0.160 | 0.072 | 0.035 |
PHQ8.c | ~ | q18_02_soc_media | c1 | 0.145 | 0.050 | 2.883 | 0.004 | 0.046 | 0.240 | 0.145 | 0.070 | 0.035 |
PHQ8.c | ~ | q02_age.c | c2 | -0.204 | 0.014 | -14.813 | 0.000 | -0.230 | -0.176 | -0.204 | -0.204 | 0.000 |
PHQ8.c | ~ | q01_gender | c3 | 0.201 | 0.033 | 6.152 | 0.000 | 0.136 | 0.269 | 0.201 | 0.083 | 0.000 |
PHQ8.c | ~ | q18_02_soc_media:q02_age.c | c4 | -0.032 | 0.024 | -1.341 | 0.180 | -0.079 | 0.018 | -0.032 | -0.019 | 1.000 |
PHQ8.c | ~ | q18_02_soc_media:q01_gender | c5 | -0.024 | 0.056 | -0.424 | 0.671 | -0.135 | 0.090 | -0.024 | -0.011 | 1.000 |
PHQ8.c | ~ | concern_index.c | b1 | 0.328 | 0.026 | 12.525 | 0.000 | 0.275 | 0.377 | 0.328 | 0.327 | 0.000 |
PHQ8.c | ~ | concern_index.c:q02_age.c | b2 | 0.026 | 0.012 | 2.226 | 0.026 | 0.003 | 0.047 | 0.026 | 0.026 | 0.182 |
PHQ8.c | ~ | concern_index.c:q01_gender | b3 | -0.009 | 0.029 | -0.326 | 0.745 | -0.065 | 0.045 | -0.009 | -0.008 | 1.000 |
parameters <- parameterEstimates(fit_mod, standardized = TRUE) %>%
filter(op == ":=") %>%
select(-c(op, lhs, rhs, std.nox))
p_adj <- p.adjust(parameters$pvalue, method = "holm")
parameters <- cbind(parameters, p_adj)
kbl(parameters) %>%
kable_classic(full_width = FALSE, lightable_options = c("striped")) %>%
row_spec(row = which(parameters$p_adj < 0.05), bold = TRUE)
label | est | se | z | pvalue | ci.lower | ci.upper | std.lv | std.all | p_adj |
---|---|---|---|---|---|---|---|---|---|
XonM | 0.097 | 0.025 | 3.825 | 0.000 | 0.047 | 0.148 | 0.097 | 0.123 | 0.004 |
MonY | 0.321 | 0.012 | 27.188 | 0.000 | 0.298 | 0.344 | 0.321 | 0.312 | 0.000 |
indirect | 0.031 | 0.008 | 3.782 | 0.000 | 0.015 | 0.048 | 0.031 | 0.038 | 0.005 |
direct | 0.126 | 0.022 | 5.620 | 0.000 | 0.083 | 0.173 | 0.126 | 0.049 | 0.000 |
total | 0.157 | 0.024 | 6.608 | 0.000 | 0.110 | 0.205 | 0.157 | 0.088 | 0.000 |
prop.mediated | 0.198 | 0.052 | 3.792 | 0.000 | 0.101 | 0.310 | 0.198 | 0.437 | 0.005 |
XonM.mean.male | -0.027 | 0.050 | -0.552 | 0.581 | -0.128 | 0.069 | -0.027 | -0.013 | 1.000 |
XonM.mean.female | 0.132 | 0.029 | 4.634 | 0.000 | 0.076 | 0.187 | 0.132 | 0.059 | 0.000 |
XonM.blw.male | -0.050 | 0.051 | -0.977 | 0.329 | -0.152 | 0.050 | -0.050 | -0.027 | 1.000 |
XonM.blw.female | 0.109 | 0.035 | 3.159 | 0.002 | 0.039 | 0.177 | 0.109 | 0.046 | 0.041 |
XonM.blw.avg | 0.074 | 0.031 | 2.366 | 0.018 | 0.011 | 0.137 | 0.074 | 0.109 | 0.288 |
XonM.abv.male | -0.005 | 0.060 | -0.077 | 0.939 | -0.130 | 0.110 | -0.005 | 0.000 | 1.000 |
XonM.abv.female | 0.155 | 0.042 | 3.678 | 0.000 | 0.071 | 0.234 | 0.155 | 0.073 | 0.007 |
XonM.abv.avg | 0.120 | 0.041 | 2.949 | 0.003 | 0.040 | 0.196 | 0.120 | 0.136 | 0.076 |
MonY.mean.male | 0.328 | 0.026 | 12.518 | 0.000 | 0.275 | 0.377 | 0.328 | 0.327 | 0.000 |
MonY.mean.female | 0.319 | 0.013 | 24.703 | 0.000 | 0.294 | 0.344 | 0.319 | 0.319 | 0.000 |
MonY.blw.male | 0.302 | 0.029 | 10.451 | 0.000 | 0.247 | 0.356 | 0.302 | 0.301 | 0.000 |
MonY.blw.female | 0.293 | 0.018 | 16.526 | 0.000 | 0.258 | 0.326 | 0.293 | 0.293 | 0.000 |
MonY.blw.avg | 0.295 | 0.017 | 17.410 | 0.000 | 0.263 | 0.327 | 0.295 | 0.285 | 0.000 |
MonY.abv.male | 0.354 | 0.029 | 12.423 | 0.000 | 0.298 | 0.410 | 0.354 | 0.353 | 0.000 |
MonY.abv.female | 0.345 | 0.017 | 20.128 | 0.000 | 0.311 | 0.380 | 0.345 | 0.345 | 0.000 |
MonY.abv.avg | 0.347 | 0.016 | 21.269 | 0.000 | 0.314 | 0.378 | 0.347 | 0.338 | 0.000 |
indirect.mean.male | -0.009 | 0.016 | -0.552 | 0.581 | -0.041 | 0.023 | -0.009 | -0.004 | 1.000 |
indirect.mean.female | 0.042 | 0.009 | 4.542 | 0.000 | 0.024 | 0.060 | 0.042 | 0.019 | 0.000 |
indirect.blw.male | -0.015 | 0.016 | -0.969 | 0.333 | -0.046 | 0.016 | -0.015 | -0.008 | 1.000 |
indirect.blw.female | 0.032 | 0.010 | 3.071 | 0.002 | 0.011 | 0.053 | 0.032 | 0.013 | 0.053 |
indirect.blw.avg | 0.022 | 0.009 | 2.330 | 0.020 | 0.003 | 0.041 | 0.022 | 0.031 | 0.297 |
indirect.abv.male | -0.002 | 0.021 | -0.077 | 0.939 | -0.043 | 0.040 | -0.002 | 0.000 | 1.000 |
indirect.abv.female | 0.053 | 0.015 | 3.642 | 0.000 | 0.024 | 0.081 | 0.053 | 0.025 | 0.008 |
indirect.abv.avg | 0.042 | 0.014 | 2.935 | 0.003 | 0.013 | 0.067 | 0.042 | 0.046 | 0.077 |
direct.mean.male | 0.145 | 0.050 | 2.882 | 0.004 | 0.047 | 0.240 | 0.145 | 0.070 | 0.085 |
direct.mean.female | 0.121 | 0.025 | 4.820 | 0.000 | 0.073 | 0.173 | 0.121 | 0.059 | 0.000 |
direct.blw.male | 0.177 | 0.053 | 3.343 | 0.001 | 0.072 | 0.281 | 0.177 | 0.089 | 0.022 |
direct.blw.female | 0.153 | 0.033 | 4.643 | 0.000 | 0.087 | 0.213 | 0.153 | 0.078 | 0.000 |
direct.blw.avg | 0.158 | 0.030 | 5.213 | 0.000 | 0.095 | 0.217 | 0.158 | 0.068 | 0.000 |
direct.abv.male | 0.113 | 0.058 | 1.931 | 0.053 | 0.002 | 0.225 | 0.113 | 0.051 | 0.641 |
direct.abv.female | 0.089 | 0.036 | 2.438 | 0.015 | 0.020 | 0.163 | 0.089 | 0.040 | 0.251 |
direct.abv.avg | 0.094 | 0.035 | 2.670 | 0.008 | 0.027 | 0.165 | 0.094 | 0.030 | 0.144 |
total.mean.male | 0.136 | 0.054 | 2.535 | 0.011 | 0.035 | 0.241 | 0.136 | 0.065 | 0.203 |
total.mean.female | 0.163 | 0.027 | 6.132 | 0.000 | 0.111 | 0.216 | 0.163 | 0.078 | 0.000 |
total.blw.male | 0.162 | 0.056 | 2.887 | 0.004 | 0.057 | 0.270 | 0.162 | 0.081 | 0.085 |
total.blw.female | 0.185 | 0.035 | 5.360 | 0.000 | 0.114 | 0.248 | 0.185 | 0.091 | 0.000 |
total.blw.avg | 0.180 | 0.032 | 5.656 | 0.000 | 0.114 | 0.241 | 0.180 | 0.100 | 0.000 |
total.abv.male | 0.111 | 0.063 | 1.770 | 0.077 | -0.005 | 0.231 | 0.111 | 0.051 | 0.844 |
total.abv.female | 0.142 | 0.039 | 3.652 | 0.000 | 0.064 | 0.222 | 0.142 | 0.065 | 0.008 |
total.abv.avg | 0.136 | 0.038 | 3.600 | 0.000 | 0.065 | 0.213 | 0.136 | 0.076 | 0.009 |
prop.med.mean.male | -0.066 | 1.150 | -0.058 | 0.954 | -0.640 | 0.156 | -0.066 | -0.066 | 1.000 |
prop.med.mean.female | 0.258 | 0.061 | 4.247 | 0.000 | 0.146 | 0.389 | 0.258 | 0.243 | 0.001 |
prop.med.blw.male | -0.094 | 0.484 | -0.194 | 0.846 | -0.505 | 0.091 | -0.094 | -0.100 | 1.000 |
prop.med.blw.female | 0.173 | 0.060 | 2.878 | 0.004 | 0.070 | 0.312 | 0.173 | 0.146 | 0.085 |
prop.med.blw.avg | 0.121 | 0.053 | 2.303 | 0.021 | 0.020 | 0.238 | 0.121 | 0.313 | 0.298 |
prop.med.abv.male | -0.015 | 1.938 | -0.008 | 0.994 | -1.877 | 0.558 | -0.015 | 0.002 | 1.000 |
prop.med.abv.female | 0.482 | 3.956 | 0.122 | 0.903 | -1.201 | 4.940 | 0.482 | 0.495 | 1.000 |
prop.med.abv.avg | 0.307 | 0.134 | 2.281 | 0.023 | 0.105 | 0.630 | 0.307 | 0.603 | 0.298 |
We are conducting clustering with the K-Medoids algorithm to find distinct groups in our dataset. We are only using a subset of the dataset, which contains values that we used in the previous step within the path model - age, gender, social media, concern index and PHQ8. The reasoning is to supplement the insights from the path model, while also limiting ourselves only to the core variables, most of them showing relatively strong effects in the regression models.
data_scaled <- data %>%
select(q01_gender, q02_age, q18_02_soc_media, PHQ8, concern_index) %>%
na.omit() %>%
scale()
set.seed(2021)
pam3 <- pam(data_scaled, 3)
# res3 <- hcut(data_scaled, k = 3, stand = TRUE, method = "median")
# fviz_dend(res3, rect = TRUE, cex = 0.5,
# k_colors = "simpsons")
# cluster_sum <- as_tibble(data_scaled) %>%
# mutate(cluster = as_factor(pam3$clustering)) %>%
# group_by(cluster) %>%
# summarise_all("mean")
fviz_cluster(pam3, data = data_scaled, ellipse.type = "norm", palette = c("#E64B3599", "#4DBBD599", "#00A08799"), ggtheme = theme_bw(), main = "Cluster plot - K-Medoids algorithm")
data_scaled_df <- as_tibble(data_scaled)
data_scaled_df$cluster <- as_factor(pam3$clustering)
data_scaled_df_long <- data_scaled_df %>% pivot_longer(cols = q01_gender:concern_index, names_to = "variable", values_to = "value")
ggplot(data_scaled_df_long, aes(x = variable, y = value, group = cluster, colour = cluster)) +
stat_summary(fun = mean, geom = "pointrange", size = 1, aes(shape = cluster)) +
stat_summary(geom = "line") +
scale_color_manual(values = c("#E64B35CC", "#4DBBD5CC", "#00A087CC")) +
ggtitle("Average value of selected variables per cluster") +
theme_bw() +
ylab("relative value")