quarto-article/code.R at main · perlatex/quarto-article · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
library(tidyverse)
library(lavaan)
library(modelsummary)
library(marginaleffects)
library(flextable)


###############################################################################
# Dataset kindly shared by Dr. Ingrid Galán for educational use
# SRS   = Autistic traits  AT
# SPAI  = Social anxiety   SA
# TOMI  = Theory of mind   TOM
# ABAS  = Social adaptive behavior  SAB

d <- readxl::read_excel("./data/data.xlsx")
sjPlot::view_df(d)
###############################################################################


###############################################################################
d %>%
  select(AT, SA, TOM,  SAB) %>%
  psych::describe() %>%
  as.data.frame() %>%
  rownames_to_column("name") %>%
  select(name, Mean = mean, SD = sd)


d %>%
  select(AT, SA, TOM, SAB) %>%
  corrr::correlate() %>%
  corrr::shave()
###############################################################################


# 结构模型检验
#################################################################################
model <- "

  TOM ~ a*AT
  SAB ~ d*AT

  SA  ~ b*TOM + cprime*AT + e*SAB

"

fit <- sem(model, data = d, se = "bootstrap", bootstrap = 2000)

fit %>%
  parameterEstimates(standardized = TRUE) %>%
  filter(op == "~") %>%
  mutate(stars = gtools::stars.pval(pvalue)) %>%
  select(label, est, se, z, pvalue, stars, ci.lower, ci.upper) %>%
  flextable::flextable() %>%
  flextable::colformat_double(digits = 3) %>%
  flextable::autofit()
###################################################################################


###################################################################################
library(modsem)

model <- "

  TOM ~ a*AT
  SAB ~ d*AT

  SA  ~ b*TOM + cprime*AT + e*SAB

"

fit_modsem <- modsem(model, data = d, se = "bootstrap", bootstrap = 2000)
summary(fit_modsem)

orders <- c("a", "d", "b", "cprime", "e")

fit_modsem %>%
  parameter_estimates() %>%
  filter(op %in% c("~", ":=")) %>%
  filter(label %in% orders) %>%
  arrange(factor(label, levels = orders)) %>%
  select(label, est, se, pvalue, ci.lower, ci.upper)
###################################################################################