Chapter 8 Latent Variable Models
# Parallel processing settings
plan(multisession(workers = parallel::detectCores(logical = FALSE)))8.1 Data
mtf <- read_rds("data/mtf.rds")
us <- read_rds("data/us.rds")
yrbs <- read_rds("data/yrbs.rds")8.2 Models
fit <- function(data, items, x, y, name, missing = "ml") {
# Center year
data <- mutate(data, Year = Year - 2017)
# Contrast code sex
data <- mutate(data, Sex = ifelse(Sex=="Male", -0.5, 0.5))
# Drop rows with missing predictor
data <- drop_na(data, all_of(x))
# Drop rows where all outcome items are missing
data <- drop_na(data, all_of(y))
# Mean-center predictors
data <- data %>%
mutate(
across(
all_of(x),
~as.numeric(scale(., center = TRUE, scale = FALSE))
)
)
# Ordered?
if (name=="YRBS") {
data <- mutate(data, across(sad_lonely:suicide_3, ordered))
missing = "listwise"
}
# Create interaction terms because lavaan doesn't know how to
newdata <- model.matrix(
as.formula(str_glue("{y} ~ Sex * Year * {x}")),
data = data
)[,-1] %>% # Take out the intercept column because it causes lavaan to break
as.data.frame()
# Interaction term breaks lavaan so change to dot
names(newdata) <- str_replace_all(names(newdata), ":", ".")
# return(newdata)
newdata <- cbind(data[,items], newdata)
# Combine names of items to a string for lavaan model
items_all <- paste0(items, collapse = " + ")
# Model strings
sem0 <- str_glue("{y} =~ {items_all}\n{y} ~ Sex + Year + Sex.Year")
sem1 <- str_glue("{y} =~ {items_all}\n{y} ~ Sex + Year + {x} + Sex.Year + Sex.{x} + Year.{x} + Sex.Year.{x}")
ml0 <- sem(sem0, data = newdata, missing = missing)
ml1 <- sem(sem1, data = newdata, missing = missing)
tibble(
data = name,
Technology = x,
Outcome = y,
ml0 = list(ml0),
ml1 = list(ml1)
)
}x1 %<-% fit(yrbs, c("sad_lonely", paste0("suicide_", 1:3)), "TV", "Suicide", "YRBS")
x2 %<-% fit(yrbs, c("sad_lonely", paste0("suicide_", 1:3)), "DV", "Suicide", "YRBS")
x3 %<-% fit(mtf, paste0("D_B_", 1:6), "TV", "Depression", "MTF")
x4 %<-% fit(mtf, paste0("D_B_", 1:6), "SM", "Depression", "MTF")
sdq_con <- c("sdqe", "sdqg", "sdql", "sdqr", "sdqv")
sdq_emo <- c("sdqc", "sdqh", "sdqm", "sdqp", "sdqx")
x5 %<-% fit(us, sdq_con, "TV", "Conduct", "US")
x6 %<-% fit(us, sdq_con, "SM", "Conduct", "US")
x7 %<-% fit(us, sdq_emo, "TV", "Emotion", "US")
x8 %<-% fit(us, sdq_emo, "SM", "Emotion", "US")
# Rename variables
fits <- bind_rows(x1,x2,x3,x4,x5,x6,x7,x8)
fits <- fits %>%
mutate(
Technology = ifelse(
Technology %in% c("SM", "DV"),
"Social media / device",
"Television"
)
) %>%
arrange(Outcome, Technology)8.3 Results
coefs <- fits %>%
mutate(p = map(ml1, ~tidy(., conf.int = TRUE))) %>%
unnest(p) %>%
filter(op == "~") %>%
separate(term, c("lhs", "rhs"), sep = " ~ ") %>%
mutate(N = map_dbl(ml1, nobs)) %>%
mutate(N = scales::comma(N, accuracy = 1))
coefs %>%
mutate(Parameter = str_replace(rhs, "SM|DV|TV", "Technology")) %>%
mutate(Parameter = str_replace_all(Parameter, "\\.", ":")) %>%
mutate(Parameter = fct_inorder(Parameter)) %>%
mutate(Outcome = fct_rev(Outcome)) %>%
ggplot(aes(estimate, Outcome, shape = Technology, fill = p.value < 0.05)) +
scale_shape_manual(values = c(21, 22)) +
scale_fill_manual(values = c("white", "black"), guide = FALSE) +
scale_x_continuous(
"Parameter estimate",
breaks = scales::pretty_breaks(),
expand = expansion(.25)
) +
geom_vline(xintercept = 0, lty = 2, size = .25) +
geom_linerangeh(
aes(xmin = conf.low, xmax = conf.high), size = .25,
position = position_dodge2v(.4), show.legend = FALSE
) +
geom_point(
size = 2, position = position_dodge2v(.4),
) +
# geom_text(
# aes(label = N), size = 2, vjust = 2,
# position = position_dodge2v(.4)
# ) +
facet_wrap("Parameter", scales = "free_x", nrow = 2) +
theme(
legend.position = c(.875, .25),
axis.title.y = element_blank(),
panel.spacing.x = unit(12, "pt")
)
# Numbers
coefs %>%
mutate(coef = str_glue("{lhs} {op} {rhs}")) %>%
select(data:Outcome, coef, estimate:conf.high)coefs %>%
distinct(data, Technology, Outcome, N)Just time
fits %>%
mutate(p = map(ml0, ~tidy(., conf.int = TRUE))) %>%
unnest(p) %>%
filter(op == "~") %>%
separate(term, c("lhs", "rhs"), sep = " ~ ") %>%
# mutate(
# Parameter = factor(rhs, levels = c("Year", "Technology", "Year x Technology"))
# ) %>%
mutate(Outcome = fct_rev(Outcome)) %>%
ggplot(aes(estimate, Outcome, shape = Technology, fill = p.value < 0.05)) +
scale_shape_manual(values = c(21, 22)) +
scale_fill_manual(values = c("white", "black"), guide = FALSE) +
scale_x_continuous(
"Parameter estimate",
breaks = scales::pretty_breaks(),
expand = expansion(.25)
) +
geom_vline(xintercept = 0, lty = 2, size = .25) +
geom_linerangeh(
aes(xmin = conf.low, xmax = conf.high), size = .25,
position = position_dodge2v(.4), show.legend = FALSE
) +
geom_point(
size = 2, position = position_dodge2v(.4),
) +
facet_wrap("rhs", scales = "free_x", nrow = 1) +
theme(
legend.position = "bottom",
axis.title.y = element_blank(),
panel.spacing.x = unit(12, "pt")
)
options(width = 120)
library(sessioninfo)
session_info()## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.3 (2020-10-10)
## os macOS Big Sur 10.16
## system x86_64, darwin17.0
## ui X11
## language (EN)
## collate en_GB.UTF-8
## ctype en_GB.UTF-8
## tz Europe/London
## date 2021-03-01
##
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
## package * version date lib source
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.0)
## backports 1.2.1 2020-12-09 [1] CRAN (R 4.0.2)
## bookdown 0.21.6 2021-03-01 [1] Github (rstudio/bookdown@ca0145f)
## broom * 0.7.5.9000 2021-03-01 [1] Github (tidymodels/broom@0b3528b)
## bslib 0.2.4 2021-01-25 [1] CRAN (R 4.0.3)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.0.0)
## cli 2.3.1 2021-02-23 [1] CRAN (R 4.0.3)
## codetools 0.2-18 2020-11-04 [1] CRAN (R 4.0.2)
## colorspace 2.0-0 2020-11-11 [1] CRAN (R 4.0.2)
## crayon 1.4.1 2021-02-08 [1] CRAN (R 4.0.3)
## DBI 1.1.1 2021-01-15 [1] CRAN (R 4.0.2)
## dbplyr 2.1.0 2021-02-03 [1] CRAN (R 4.0.2)
## digest 0.6.27 2020-10-24 [1] CRAN (R 4.0.2)
## dplyr * 1.0.4 2021-02-02 [1] CRAN (R 4.0.2)
## ellipsis 0.3.1 2020-05-15 [1] CRAN (R 4.0.0)
## evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.0)
## fansi 0.4.2 2021-01-15 [1] CRAN (R 4.0.2)
## farver 2.0.3 2020-01-16 [1] CRAN (R 4.0.0)
## forcats * 0.5.1 2021-01-27 [1] CRAN (R 4.0.2)
## fs 1.5.0 2020-07-31 [1] CRAN (R 4.0.2)
## future * 1.21.0 2020-12-10 [1] CRAN (R 4.0.2)
## generics 0.1.0 2020-10-31 [1] CRAN (R 4.0.2)
## ggplot2 * 3.3.3 2020-12-30 [1] CRAN (R 4.0.2)
## ggstance * 0.3.5 2020-12-17 [1] CRAN (R 4.0.2)
## globals 0.14.0 2020-11-22 [1] CRAN (R 4.0.2)
## glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 4.0.0)
## haven 2.3.1 2020-06-01 [1] CRAN (R 4.0.0)
## highr 0.8 2019-03-20 [1] CRAN (R 4.0.0)
## hms 1.0.0 2021-01-13 [1] CRAN (R 4.0.2)
## htmltools 0.5.1.1 2021-01-22 [1] CRAN (R 4.0.2)
## httr 1.4.2 2020-07-20 [1] CRAN (R 4.0.2)
## jquerylib 0.1.3 2020-12-17 [1] CRAN (R 4.0.2)
## jsonlite 1.7.2 2020-12-09 [1] CRAN (R 4.0.2)
## knitr * 1.31 2021-01-27 [1] CRAN (R 4.0.2)
## lavaan * 0.6-7 2020-07-31 [1] CRAN (R 4.0.2)
## lifecycle 1.0.0 2021-02-15 [1] CRAN (R 4.0.2)
## listenv 0.8.0 2019-12-05 [1] CRAN (R 4.0.0)
## lubridate 1.7.9.2 2020-11-13 [1] CRAN (R 4.0.2)
## magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.0.2)
## mnormt 2.0.2 2020-09-01 [1] CRAN (R 4.0.2)
## modelr 0.1.8 2020-05-19 [1] CRAN (R 4.0.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.0.0)
## pacman 0.5.1 2019-03-11 [1] CRAN (R 4.0.0)
## parallelly 1.23.0 2021-01-04 [1] CRAN (R 4.0.2)
## pbivnorm 0.6.0 2015-01-23 [1] CRAN (R 4.0.0)
## pillar 1.5.0 2021-02-22 [1] CRAN (R 4.0.3)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.0)
## purrr * 0.3.4 2020-04-17 [1] CRAN (R 4.0.0)
## R6 2.5.0 2020-10-28 [1] CRAN (R 4.0.2)
## Rcpp 1.0.6 2021-01-15 [1] CRAN (R 4.0.2)
## readr * 1.4.0 2020-10-05 [1] CRAN (R 4.0.2)
## readxl 1.3.1 2019-03-13 [1] CRAN (R 4.0.0)
## reprex 1.0.0 2021-01-27 [1] CRAN (R 4.0.2)
## rlang 0.4.10 2020-12-30 [1] CRAN (R 4.0.2)
## rmarkdown 2.7.2 2021-03-01 [1] Github (rstudio/rmarkdown@9bfaf4a)
## rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.0.2)
## rvest 0.3.6 2020-07-25 [1] CRAN (R 4.0.2)
## sass 0.3.1 2021-01-24 [1] CRAN (R 4.0.2)
## scales 1.1.1 2020-05-11 [1] CRAN (R 4.0.0)
## sessioninfo * 1.1.1 2018-11-05 [1] CRAN (R 4.0.0)
## stringi 1.5.3 2020-09-09 [1] CRAN (R 4.0.2)
## stringr * 1.4.0 2019-02-10 [1] CRAN (R 4.0.0)
## tibble * 3.1.0 2021-02-25 [1] CRAN (R 4.0.2)
## tidyr * 1.1.2 2020-08-27 [1] CRAN (R 4.0.2)
## tidyselect 1.1.0 2020-05-11 [1] CRAN (R 4.0.0)
## tidyverse * 1.3.0 2019-11-21 [1] CRAN (R 4.0.0)
## tmvnsim 1.0-2 2016-12-15 [1] CRAN (R 4.0.0)
## utf8 1.1.4 2018-05-24 [1] CRAN (R 4.0.0)
## vctrs 0.3.6 2020-12-17 [1] CRAN (R 4.0.2)
## withr 2.4.1 2021-01-26 [1] CRAN (R 4.0.2)
## xfun 0.21 2021-02-10 [1] CRAN (R 4.0.2)
## xml2 1.3.2 2020-04-23 [1] CRAN (R 4.0.0)
## yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.0)
##
## [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library