环境与依赖
# 自动安装并加载常用包
pkgs <- c(
"tidyverse","rio","readxl","janitor","stringr","psych",
"lavaan","semPlot","lme4","lmerTest","BayesFactor",
"smacof","GGally","corrplot","jmv"
)
to_install <- pkgs[!sapply(pkgs, requireNamespace, quietly = TRUE)]
if (length(to_install)) install.packages(to_install, repos = "https://cloud.r-project.org")
invisible(lapply(pkgs, library, character.only = TRUE))
knitr::opts_chunk$set(message = FALSE, warning = FALSE, fig.align = "center")
set.seed(123)
R 与 RStudio 获取
cat("CRAN:", "https://cran.r-project.org/", "\n")
## CRAN: https://cran.r-project.org/
cat("RStudio (Posit):", "https://posit.co/download/rstudio-desktop/", "\n")
## RStudio (Posit): https://posit.co/download/rstudio-desktop/
R
的优点与生态（心理学常用）
- 可复现：Rmd/Quarto 一键生成报告
- tidyverse：数据清洗与可视化（dplyr/ggplot2
等）
- jmv：jamovi 的 R 后端，菜单式常用分析
- psych（心理统计/测量）、lavaan（SEM）、lme4/lmerTest（混合效应）、BayesFactor（贝叶斯）、smacof（MDS）
RStudio 常用快捷键
cat("运行当前行/选中：Ctrl/Cmd + Enter\n",
"插入管道：Ctrl/Cmd + Shift + M\n",
"注释/取消注释：Ctrl/Cmd + Shift + C\n",
"清屏：Ctrl/Cmd + L；撤销：Ctrl/Cmd + Z\n")
## 运行当前行/选中：Ctrl/Cmd + Enter
## 插入管道：Ctrl/Cmd + Shift + M
## 注释/取消注释：Ctrl/Cmd + Shift + C
## 清屏：Ctrl/Cmd + L；撤销：Ctrl/Cmd + Z
基础语法与数据结构
赋值与向量化
x <- 1:5
mean(x)
## [1] 3
x^2 + 1
## [1] 2 5 10 17 26
向量（vector）
a <- c(1,2,3,4,5)
b <- c("one","two","three")
gender <- c("male","female","male","male")
a[3]
## [1] 3
a[c(1,3,5)]
## [1] 1 3 5
a[2:4]
## [1] 2 3 4
矩阵（matrix）
y <- matrix(1:20, nrow = 5, ncol = 4)
y
## [,1] [,2] [,3] [,4]
## [1,] 1 6 11 16
## [2,] 2 7 12 17
## [3,] 3 8 13 18
## [4,] 4 9 14 19
## [5,] 5 10 15 20
cells <- c(1,26,24,68)
rnames <- c("R1","R2")
cnames <- c("C1","C2")
mymat <- matrix(cells, nrow = 2, ncol = 2, byrow = TRUE,
dimnames = list(rnames, cnames))
mymat
## C1 C2
## R1 1 26
## R2 24 68
x <- matrix(1:10, nrow = 2)
x[2, ]
## [1] 2 4 6 8 10
x[, 2]
## [1] 3 4
x[1, 4]
## [1] 7
数组（array，三维及以上）
dim1 <- c("A1","A2")
dim2 <- c("B1","B2","B3")
dim3 <- c("C1","C2","C3","C4")
Z <- array(1:24, dim = c(2,3,4), dimnames = list(dim1, dim2, dim3))
Z[,"B2",]
## C1 C2 C3 C4
## A1 3 9 15 21
## A2 4 10 16 22
数据框（data.frame）与因子（factor）
stuID <- c(1,2,3,4)
age <- c(25,34,28,52)
personality <- c("TYPE1","TYPE2","TYPE3","TYPE1")
study <- c("poor","improved","excellent","poor")
studata <- data.frame(stuID, age, personality, study, stringsAsFactors = FALSE)
# 选取
studata[1:2]
studata[c("personality","study")]
studata$age
## [1] 25 34 28 52
# 有序因子与重编码
studata$study <- factor(studata$study,
levels = c("poor","improved","excellent"),
ordered = TRUE)
studata <- studata |>
dplyr::mutate(age_group = dplyr::case_when(
age <= 20 ~ "低",
age <= 30 ~ "中",
TRUE ~ "高"
))
table(studata$age_group)
##
## 中 高
## 2 2
列表（list）
g <- "my first list"; h <- c(25,26,18,39); j <- matrix(1:10, nrow = 5)
k <- c("one","two","three")
mylist <- list(title = g, ages = h, mat = j, labels = k)
str(mylist)
## List of 4
## $ title : chr "my first list"
## $ ages : num [1:4] 25 26 18 39
## $ mat : int [1:5, 1:2] 1 2 3 4 5 6 7 8 9 10
## $ labels: chr [1:3] "one" "two" "three"
数据导入与导出（rio）
getwd()
## [1] "/Users/mac/Documents/小红书/心理学学习笔记"
# setwd("/your/path")
# 读写（自动识别后缀）：
# dat <- rio::import("data.xlsx") # 读
# rio::export(dat, "data_out.sav") # 写
数据清洗与管理（dplyr +
janitor）
df <- airquality |>
janitor::clean_names() |>
tibble::as_tibble()
# 缺失统计与处理
sum(is.na(df$ozone))
## [1] 37
df2 <- df |> tidyr::drop_na(ozone)
# 新变量与重命名
df3 <- df2 |>
dplyr::mutate(ozone_cat = cut(ozone, breaks = c(-Inf, 50, 100, Inf),
labels = c("低","中","高"))) |>
dplyr::rename(temp_f = temp)
# 排序与选择
df4 <- df3 |>
dplyr::arrange(dplyr::desc(temp_f)) |>
dplyr::select(ozone, ozone_cat, wind, temp_f)
head(df4)
描述统计与频率表
desc_tab <- psych::describe(df4 |> dplyr::select(ozone, wind, temp_f))
round(desc_tab, 3)
table(df3$ozone_cat)
##
## 低 中 高
## 82 27 7
prop.table(table(df3$ozone_cat))
##
## 低 中 高
## 0.70689655 0.23275862 0.06034483
xt <- table(df3$ozone_cat, cut(df3$wind, breaks = 3))
addmargins(xt)
##
## (2.28,8.43] (8.43,14.6] (14.6,20.7] Sum
## 低 18 51 13 82
## 中 19 8 0 27
## 高 7 0 0 7
## Sum 44 59 13 116
可视化（ggplot2 /
corrplot）
library(showtext)
font_add("Heiti TC Light", regular = "STHeiti Light.ttc")
showtext_auto()
ggplot(df3, aes(ozone)) +
geom_histogram(bins = 20) +
geom_density(linewidth = 0.7) +
labs(title = "Ozone 分布", x = "Ozone", y = "频数/密度")
ggplot(df3, aes(temp_f, ozone)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE) +
labs(title = "温度与臭氧", x = "温度 (F)", y = "Ozone")
num_df <- df3 |> dplyr::select(ozone, wind, temp_f) |> tidyr::drop_na()
Rmat <- cor(num_df)
corrplot::corrplot(Rmat, method = "color", addCoef.col = "black", number.cex = .7)
信度分析（Cronbach’s
α）
data(bfi, package = "psych")
bfi_sub <- bfi |> dplyr::select(A1:A5) |> tidyr::drop_na()
alpha_res <- psych::alpha(bfi_sub)
## Some items ( A1 ) were negatively correlated with the first principal component and
## probably should be reversed.
## To do this, run the function again with the 'check.keys=TRUE' option
alpha_res$total$raw_alpha
## [1] 0.4306169
alpha_res
##
## Reliability analysis
## Call: psych::alpha(x = bfi_sub)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.43 0.46 0.53 0.14 0.84 0.017 4.2 0.74 0.32
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.4 0.43 0.46
## Duhachek 0.4 0.43 0.46
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## A1 0.72 0.73 0.67 0.397 2.64 0.0089 0.0066 0.375
## A2 0.28 0.30 0.39 0.096 0.42 0.0223 0.1106 0.079
## A3 0.17 0.20 0.31 0.060 0.25 0.0253 0.1015 0.079
## A4 0.25 0.30 0.44 0.098 0.43 0.0233 0.1614 0.103
## A5 0.21 0.23 0.36 0.071 0.31 0.0242 0.1322 0.093
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## A1 2709 0.055 0.021 -0.40 -0.31 2.4 1.4
## A2 2709 0.631 0.665 0.58 0.37 4.8 1.2
## A3 2709 0.728 0.743 0.72 0.48 4.6 1.3
## A4 2709 0.689 0.661 0.50 0.37 4.7 1.5
## A5 2709 0.701 0.718 0.64 0.45 4.6 1.3
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## A1 0.33 0.30 0.14 0.12 0.08 0.03 0
## A2 0.02 0.05 0.05 0.20 0.37 0.31 0
## A3 0.03 0.06 0.08 0.20 0.36 0.27 0
## A4 0.05 0.08 0.07 0.16 0.24 0.41 0
## A5 0.02 0.07 0.09 0.22 0.35 0.25 0
探索性因素分析（EFA）
R <- cor(bfi_sub)
psych::KMO(R)
## Kaiser-Meyer-Olkin factor adequacy
## Call: psych::KMO(r = R)
## Overall MSA = 0.76
## MSA for each item =
## A1 A2 A3 A4 A5
## 0.78 0.76 0.73 0.83 0.76
psych::cortest.bartlett(R, n = nrow(bfi_sub))
## $chisq
## [1] 2530.363
##
## $p.value
## [1] 0
##
## $df
## [1] 10
fa.parallel(bfi_sub, fm = "pa", fa = "fa")
## Parallel analysis suggests that the number of factors = 2 and the number of components = NA
fit_efa <- fa(bfi_sub, nfactors = 2, fm = "pa", rotate = "varimax")
print(fit_efa, cut = .35, sort = TRUE)
## Factor Analysis using method = pa
## Call: fa(r = bfi_sub, nfactors = 2, rotate = "varimax", fm = "pa")
## Standardized loadings (pattern matrix) based upon correlation matrix
## item PA1 PA2 h2 u2 com
## A3 3 0.68 0.58 0.42 1.4
## A5 5 0.64 0.44 0.56 1.2
## A4 4 0.43 0.23 0.77 1.5
## A2 2 0.44 0.58 0.53 0.47 1.9
## A1 1 -0.48 0.25 0.75 1.2
##
## PA1 PA2
## SS loadings 1.28 0.76
## Proportion Var 0.26 0.15
## Cumulative Var 0.26 0.41
## Proportion Explained 0.63 0.37
## Cumulative Proportion 0.63 1.00
##
## Mean item complexity = 1.4
## Test of the hypothesis that 2 factors are sufficient.
##
## df null model = 10 with the objective function = 0.94 with Chi Square = 2530.36
## df of the model are 1 and the objective function was 0
##
## The root mean square of the residuals (RMSR) is 0.01
## The df corrected root mean square of the residuals is 0.03
##
## The harmonic n.obs is 2709 with the empirical chi square 5.31 with prob < 0.021
## The total n.obs was 2709 with Likelihood Chi Square = 6.7 with prob < 0.0096
##
## Tucker Lewis Index of factoring reliability = 0.977
## RMSEA index = 0.046 and the 90 % confidence intervals are 0.018 0.081
## BIC = -1.2
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## PA1 PA2
## Correlation of (regression) scores with factors 0.78 0.66
## Multiple R square of scores with factors 0.61 0.43
## Minimum correlation of possible factor scores 0.23 -0.14
fa.diagram(fit_efa, digits = 2)
验证性因素分析（CFA,
lavaan）
data("HolzingerSwineford1939", package = "lavaan")
cfa_model <- '
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
'
fit_cfa <- lavaan::cfa(cfa_model, data = HolzingerSwineford1939)
summary(fit_cfa, fit.measures = TRUE, standardized = TRUE, rsquare = TRUE)
## lavaan 0.6-19 ended normally after 35 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 21
##
## Number of observations 301
##
## Model Test User Model:
##
## Test statistic 85.306
## Degrees of freedom 24
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 918.852
## Degrees of freedom 36
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.931
## Tucker-Lewis Index (TLI) 0.896
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3737.745
## Loglikelihood unrestricted model (H1) -3695.092
##
## Akaike (AIC) 7517.490
## Bayesian (BIC) 7595.339
## Sample-size adjusted Bayesian (SABIC) 7528.739
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.092
## 90 Percent confidence interval - lower 0.071
## 90 Percent confidence interval - upper 0.114
## P-value H_0: RMSEA <= 0.050 0.001
## P-value H_0: RMSEA >= 0.080 0.840
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.065
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## visual =~
## x1 1.000 0.900 0.772
## x2 0.554 0.100 5.554 0.000 0.498 0.424
## x3 0.729 0.109 6.685 0.000 0.656 0.581
## textual =~
## x4 1.000 0.990 0.852
## x5 1.113 0.065 17.014 0.000 1.102 0.855
## x6 0.926 0.055 16.703 0.000 0.917 0.838
## speed =~
## x7 1.000 0.619 0.570
## x8 1.180 0.165 7.152 0.000 0.731 0.723
## x9 1.082 0.151 7.155 0.000 0.670 0.665
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## visual ~~
## textual 0.408 0.074 5.552 0.000 0.459 0.459
## speed 0.262 0.056 4.660 0.000 0.471 0.471
## textual ~~
## speed 0.173 0.049 3.518 0.000 0.283 0.283
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.549 0.114 4.833 0.000 0.549 0.404
## .x2 1.134 0.102 11.146 0.000 1.134 0.821
## .x3 0.844 0.091 9.317 0.000 0.844 0.662
## .x4 0.371 0.048 7.779 0.000 0.371 0.275
## .x5 0.446 0.058 7.642 0.000 0.446 0.269
## .x6 0.356 0.043 8.277 0.000 0.356 0.298
## .x7 0.799 0.081 9.823 0.000 0.799 0.676
## .x8 0.488 0.074 6.573 0.000 0.488 0.477
## .x9 0.566 0.071 8.003 0.000 0.566 0.558
## visual 0.809 0.145 5.564 0.000 1.000 1.000
## textual 0.979 0.112 8.737 0.000 1.000 1.000
## speed 0.384 0.086 4.451 0.000 1.000 1.000
##
## R-Square:
## Estimate
## x1 0.596
## x2 0.179
## x3 0.338
## x4 0.725
## x5 0.731
## x6 0.702
## x7 0.324
## x8 0.523
## x9 0.442
semPlot::semPaths(fit_cfa, what = "path", layout = "tree2",
whatLabels = "std", style = "ram", nDigits = 2)
多维尺度分析（MDS）
dist_mat <- dist(scale(attitude)) # 欧氏距离
mds_fit <- cmdscale(dist_mat, k = 2) # 经典 MDS
mds_df <- tibble::as_tibble(mds_fit) |> stats::setNames(c("Dim1","Dim2"))
ggplot(mds_df, aes(Dim1, Dim2)) + geom_point() + labs(title = "经典 MDS 二维解")
# # 非度量 MDS（需要 smacof 包）
# smacof_fit <- smacof::mds(dist_mat, type = "ordinal", ndim = 2)
# plot(smacof_fit)
混合/多层线性模型（lme4 + lmerTest）
data("sleepstudy", package = "lme4")
m1 <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
summary(m1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Reaction ~ Days + (Days | Subject)
## Data: sleepstudy
##
## REML criterion at convergence: 1743.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9536 -0.4634 0.0231 0.4634 5.1793
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 612.10 24.741
## Days 35.07 5.922 0.07
## Residual 654.94 25.592
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 251.405 6.825 17.000 36.838 < 2e-16 ***
## Days 10.467 1.546 17.000 6.771 3.26e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.138
贝叶斯示例（BayesFactor）
pg <- PlantGrowth
bf <- BayesFactor::ttestBF(
x = pg$weight[pg$group == "ctrl"],
y = pg$weight[pg$group == "trt1"]
)
bf
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 0.650041 ±0%
##
## Against denominator:
## Null, mu1-mu2 = 0
## ---
## Bayes factor type: BFindepSample, JZS
jmv：描述/检验分析快速统计
# 示例（注：jmv 多以“菜单函数”形式组织）：
jmv::descriptives(data = iris, vars = vars(Sepal.Length, Sepal.Width))
jmv::ttestIS(data = ToothGrowth, vars = "len", group = "supp")
jmv::anova(data = iris, deps = "Sepal.Length", factors = "Species")
常见问题与排查
cat("• 包函数冲突：用 pkg::fun()，如 psych::alpha()\n",
"• 缺失与类型：is.na()/drop_na()/as.numeric()/as.factor()\n",
"• 复现：set.seed() + sessionInfo()\n",
"• Knit 失败：逐块运行定位问题包/数据\n")
## • 包函数冲突：用 pkg::fun()，如 psych::alpha()
## • 缺失与类型：is.na()/drop_na()/as.numeric()/as.factor()
## • 复现：set.seed() + sessionInfo()
## • Knit 失败：逐块运行定位问题包/数据
参考与延伸
cat("R for Data Science:", "https://r4ds.hadley.nz/", "\n",
"psych 包:", "https://cran.r-project.org/package=psych", "\n",
"lavaan:", "https://lavaan.ugent.be/", "\n",
"MDS 原理与实现（示例）:", "https://www.jianshu.com/p/39332c72828e", "\n")
## R for Data Science: https://r4ds.hadley.nz/
## psych 包: https://cran.r-project.org/package=psych
## lavaan: https://lavaan.ugent.be/
## MDS 原理与实现（示例）: https://www.jianshu.com/p/39332c72828e
sessionInfo()
## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Asia/Shanghai
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] showtext_0.9-7 showtextdb_3.0 sysfonts_0.8.9
## [4] jmv_2.7.7 corrplot_0.95 GGally_2.4.0
## [7] smacof_2.1-7 e1071_1.7-16 colorspace_2.1-1
## [10] plotrix_3.8-4 BayesFactor_0.9.12-4.7 coda_0.19-4.1
## [13] lmerTest_3.1-3 lme4_1.1-37 Matrix_1.7-4
## [16] semPlot_1.1.7 lavaan_0.6-19 psych_2.5.6
## [19] janitor_2.2.1 readxl_1.4.5 rio_1.2.3
## [22] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.2
## [25] dplyr_1.1.4 purrr_1.1.0 readr_2.1.5
## [28] tidyr_1.3.1 tibble_3.3.0 ggplot2_4.0.0
## [31] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 rstudioapi_0.17.1 jsonlite_2.0.0
## [4] shape_1.4.6.1 magrittr_2.0.4 jomo_2.7-6
## [7] farver_2.1.2 nloptr_2.2.1 rmarkdown_2.29
## [10] vctrs_0.6.5 minqa_1.2.8 base64enc_0.1-3
## [13] htmltools_0.5.8.1 polynom_1.4-1 weights_1.1.2
## [16] broom_1.0.10 cellranger_1.1.0 Formula_1.2-5
## [19] mitml_0.4-5 sass_0.4.10 bslib_0.9.0
## [22] htmlwidgets_1.6.4 plyr_1.8.9 cachem_1.1.0
## [25] igraph_2.1.4 lifecycle_1.0.4 iterators_1.0.14
## [28] pkgconfig_2.0.3 R6_2.6.1 fastmap_1.2.0
## [31] rbibutils_2.3 snakecase_0.11.1 digest_0.6.37
## [34] numDeriv_2016.8-1.1 OpenMx_2.21.13 fdrtool_1.2.18
## [37] ellipse_0.5.0 Hmisc_5.2-3 labeling_0.4.3
## [40] timechange_0.3.0 nnls_1.6 gdata_3.0.1
## [43] mgcv_1.9-3 abind_1.4-8 compiler_4.4.2
## [46] proxy_0.4-27 withr_3.0.2 doParallel_1.0.17
## [49] glasso_1.11 htmlTable_2.4.3 S7_0.2.0
## [52] backports_1.5.0 carData_3.0-5 ggstats_0.11.0
## [55] pan_1.9 MASS_7.3-65 corpcor_1.6.10
## [58] gtools_3.9.5 tools_4.4.2 pbivnorm_0.6.0
## [61] foreign_0.8-90 zip_2.3.3 nnet_7.3-20
## [64] glue_1.8.0 quadprog_1.5-8 nlme_3.1-168
## [67] lisrelToR_0.3 grid_4.4.2 checkmate_2.3.3
## [70] cluster_2.1.8.1 reshape2_1.4.4 generics_0.1.4
## [73] gtable_0.3.6 tzdb_0.5.0 class_7.3-23
## [76] hms_1.1.3 data.table_1.17.8 sem_3.1-16
## [79] foreach_1.5.2 pillar_1.11.1 rockchalk_1.8.157
## [82] splines_4.4.2 lattice_0.22-7 survival_3.8-3
## [85] kutils_1.73 tidyselect_1.2.1 pbapply_1.7-4
## [88] knitr_1.50 reformulas_0.4.1 gridExtra_2.3
## [91] stats4_4.4.2 xfun_0.53 qgraph_1.9.8
## [94] jmvcore_2.7.7 arm_1.14-4 stringi_1.8.7
## [97] yaml_2.3.10 boot_1.3-32 evaluate_1.0.5
## [100] codetools_0.2-20 wordcloud_2.6 mi_1.2
## [103] cli_3.6.5 RcppParallel_5.1.11-1 rpart_4.1.24
## [106] xtable_1.8-4 Rdpack_2.6.4 jquerylib_0.1.4
## [109] Rcpp_1.1.0 png_0.1-8 XML_3.99-0.19
## [112] parallel_4.4.2 MatrixModels_0.5-4 jpeg_0.1-11
## [115] glmnet_4.1-10 mvtnorm_1.3-3 scales_1.4.0
## [118] openxlsx_4.2.8 rlang_1.1.6 mnormt_2.1.1
## [121] mice_3.18.0