SEED <- 123; set.seed(SEED)
library(readxl)
library(dplyr)
##
## 载入程辑包：'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(stringr)
library(lme4)
## 载入需要的程辑包：Matrix
##
## 载入程辑包：'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lmerTest)
##
## 载入程辑包：'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(broom.mixed)
library(purrr)
path <- "C:/Users/97380/Documents/LFF_data process/Scale/BADS/BADS_fanxiang.xlsx"
bads_raw <- read_excel(path, sheet = "Sheet1")
names(bads_raw) <- trimws(names(bads_raw))
# 把数据从 wide 变 long：提取 time 与 量表维度
bads_long <- bads_raw %>%
pivot_longer(
cols = matches("^(pre|post|fu)_"),
names_to = c("time","scale"),
names_pattern = "^(pre|post|fu)_(.+)$",
values_to = "score"
) %>%
mutate(
time = factor(time, levels = c("pre","post","fu")),
group = factor(group, levels = c(1,2), labels = c("控制组","实验组"))
)
mk_fit <- function(df, times, label){
d <- df %>%
dplyr::filter(time %in% times) %>%
droplevels() %>%
dplyr::mutate(time = factor(time, levels = times),
ID = factor(ID))
fit <- lmer(score ~ group * time + (1|ID), data = d, REML = TRUE)
list(fit = fit, label = label, times = times)
}
## 三个模型的“固定效应摘要”·按分量表分别输出
## 三个模型的“固定效应摘要”·按分量表分别输出（无CI，含 SE 与 SE_d）
fixed_tbl <- bads_long %>%
split(.$scale) %>% # 按分量表拆分
purrr::imap_dfr(function(df_sc, sc_name){ # sc_name 为分量表名
fits <- list(
mk_fit(df_sc, c("pre","post"), "pre-post"),
mk_fit(df_sc, c("pre","fu"), "pre-fu"),
mk_fit(df_sc, c("post","fu"), "post-fu")
)
purrr::map_dfr(fits, function(x){
m <- x$fit; lab <- x$label
sig <- sigma(m)
co <- as.data.frame(coef(summary(m)))
co$term <- rownames(co)
co <- co %>%
dplyr::mutate(
d = Estimate / sig, # 标准化效应
SE_d = `Std. Error` / sig # d 的标准误
) %>%
dplyr::select(
term, Estimate, `Std. Error`, df, `t value`, `Pr(>|t|)`, d, SE_d
)
tibble::tibble(scale = sc_name, model = lab) %>%
dplyr::bind_cols(co)
})
})
fixed_tbl
## # A tibble: 48 × 10
## scale model term Estimate `Std. Error` df `t value` `Pr(>|t|)` d
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BADS_工… pre-… (Int… 7.67 0.697 84.8 11.0 5.44e-18 3.20
## 2 BADS_工… pre-… grou… 0.167 0.986 84.8 0.169 8.66e- 1 0.0695
## 3 BADS_工… pre-… time… -0.500 0.619 58.0 -0.808 4.23e- 1 -0.209
## 4 BADS_工… pre-… grou… -0.0667 0.876 58.0 -0.0761 9.40e- 1 -0.0278
## 5 BADS_工… pre-… (Int… 7.67 0.665 91.2 11.5 1.71e-19 3.04
## 6 BADS_工… pre-… grou… 0.167 0.941 91.2 0.177 8.60e- 1 0.0661
## 7 BADS_工… pre-… time… -0.908 0.651 58.0 -1.40 1.68e- 1 -0.360
## 8 BADS_工… pre-… grou… 1.10 0.920 58.0 1.19 2.38e- 1 0.436
## 9 BADS_工… post… (Int… 7.17 0.672 88.2 10.7 1.52e-17 2.94
## 10 BADS_工… post… grou… 0.100 0.950 88.2 0.105 9.16e- 1 0.0410
## # ℹ 38 more rows
## # ℹ 1 more variable: SE_d <dbl>
## === 2) 仅做 pre–post 与 pre–fu 的“组内简单效应”
within_two <- bads_long %>%
split(.$scale) %>%
purrr::imap_dfr(function(df_sc, sc_name){
fits <- list(
mk_fit(df_sc, c("pre","post"), "pre-post"),
mk_fit(df_sc, c("pre","fu"), "pre-fu")
)
purrr::map_dfr(fits, function(x){
m <- x$fit; lab <- x$label; sig <- sigma(m)
emm <- emmeans(m, ~ time | group)
contr_name <- if (lab == "pre-post") "post - pre" else "fu - pre"
w <- summary(
contrast(emm, method = setNames(list(c(-1,1)), contr_name),
by = "group", adjust = "none"),
infer = c(TRUE, TRUE)
) %>% as.data.frame()
tibble::tibble(
scale = sc_name,
model = lab,
group = w$group,
contrast = w$contrast,
estimate = w$estimate,
ci_low = w$lower.CL,
ci_high = w$upper.CL,
p = w$p.value,
d = w$estimate / sig,
d_low = d + qt(0.025, w$df) * (w$SE / sig),
d_high = d + qt(0.975, w$df) * (w$SE / sig)
)
})
})
within_two
## # A tibble: 16 × 11
## scale model group contrast estimate ci_low ci_high p d d_low
## <chr> <chr> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BADS_工… pre-… 控制… post - … -0.500 -1.74 0.739 4.23e-1 -0.209 -0.725
## 2 BADS_工… pre-… 实验… post - … -0.567 -1.81 0.673 3.64e-1 -0.236 -0.753
## 3 BADS_工… pre-… 控制… fu - pre -0.908 -2.21 0.395 1.68e-1 -0.360 -0.877
## 4 BADS_工… pre-… 实验… fu - pre 0.190 -1.11 1.49 7.71e-1 0.0754 -0.441
## 5 BADS_回… pre-… 控制… post - … 1.03 -0.843 2.91 2.75e-1 0.285 -0.232
## 6 BADS_回… pre-… 实验… post - … 1 -0.876 2.88 2.90e-1 0.275 -0.241
## 7 BADS_回… pre-… 控制… fu - pre 0.821 -1.58 3.22 4.96e-1 0.177 -0.340
## 8 BADS_回… pre-… 实验… fu - pre 1.07 -1.33 3.47 3.76e-1 0.230 -0.286
## 9 BADS_激活 pre-… 控制… post - … 0.533 -1.13 2.20 5.23e-1 0.166 -0.351
## 10 BADS_激活 pre-… 实验… post - … 2.43 0.770 4.10 4.85e-3 0.756 0.239
## 11 BADS_激活 pre-… 控制… fu - pre 0.703 -1.26 2.67 4.76e-1 0.185 -0.332
## 12 BADS_激活 pre-… 实验… fu - pre 3.51 1.54 5.47 7.19e-4 0.922 0.406
## 13 BADS_社… pre-… 控制… post - … -1.23 -2.73 0.267 1.05e-1 -0.425 -0.942
## 14 BADS_社… pre-… 实验… post - … 2.71 1.21 4.21 6.26e-4 0.934 0.417
## 15 BADS_社… pre-… 控制… fu - pre 0.467 -1.17 2.10 5.70e-1 0.148 -0.369
## 16 BADS_社… pre-… 实验… fu - pre 2.62 0.991 4.26 2.13e-3 0.830 0.313
## # ℹ 1 more variable: d_high <dbl>
# 从 bads_long 汇总出“总分”
bads_total <- bads_long %>%
dplyr::group_by(ID, group, time) %>%
dplyr::summarise(score = sum(score, na.rm = TRUE), .groups = "drop") %>%
dplyr::mutate(scale = "BADS_总分")
fits_total <- list(
mk_fit(bads_total, c("pre","post"), "pre-post"),
mk_fit(bads_total, c("pre","fu"), "pre-fu"),
mk_fit(bads_total, c("post","fu"), "post-fu")
)
fixed_total <- purrr::map_dfr(fits_total, function(x){
m <- x$fit; lab <- x$label
sig <- sigma(m)
co <- as.data.frame(coef(summary(m)))
co$term <- rownames(co)
co <- co %>%
dplyr::mutate(
d = Estimate / sig,
SE_d = `Std. Error` / sig
) %>%
dplyr::select(term, Estimate, `Std. Error`, df, `t value`, `Pr(>|t|)`, d, SE_d)
tibble::tibble(scale = "BADS_总分", model = lab) %>%
dplyr::bind_cols(co)
})
fixed_total
## # A tibble: 12 × 10
## scale model term Estimate `Std. Error` df `t value` `Pr(>|t|)` d
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BADS_总… pre-… (Int… 61.4 2.61 75.6 23.5 1.86e-36 8.28
## 2 BADS_总… pre-… grou… 4.20 3.70 75.6 1.14 2.60e- 1 0.566
## 3 BADS_总… pre-… time… -0.167 1.92 58.0 -0.0870 9.31e- 1 -0.0225
## 4 BADS_总… pre-… grou… 5.74 2.71 58.0 2.12 3.83e- 2 0.774
## 5 BADS_总… pre-… (Int… 61.4 2.49 82.6 24.6 7.73e-40 7.45
## 6 BADS_总… pre-… grou… 4.20 3.52 82.6 1.19 2.37e- 1 0.510
## 7 BADS_总… pre-… time… 1.08 2.13 58.0 0.509 6.13e- 1 0.131
## 8 BADS_总… pre-… grou… 6.31 3.01 58.0 2.10 4.05e- 2 0.765
## 9 BADS_总… post… (Int… 61.2 2.53 81.1 24.2 8.06e-39 7.53
## 10 BADS_总… post… grou… 9.94 3.58 81.1 2.78 6.82e- 3 1.22
## 11 BADS_总… post… time… 1.25 2.10 58.0 0.595 5.54e- 1 0.154
## 12 BADS_总… post… grou… 0.562 2.97 58.0 0.189 8.51e- 1 0.0691
## # ℹ 1 more variable: SE_d <dbl>
# === 2) “总分”的组内简单效应（只做 pre-post / pre-fu；含 95%CI、p、d）===
within_total <- purrr::map_dfr(fits_total[1:2], function(x){
m <- x$fit; lab <- x$label; sig <- sigma(m)
emm <- emmeans(m, ~ time | group)
contr_name <- if (lab == "pre-post") "post - pre" else "fu - pre"
w <- summary(
contrast(emm, method = setNames(list(c(-1,1)), contr_name),
by = "group", adjust = "none"),
infer = c(TRUE, TRUE)
) %>% as.data.frame()
tibble::tibble(
scale = "BADS_总分",
model = lab,
group = w$group,
contrast = w$contrast, # post - pre / fu - pre
estimate = w$estimate,
ci_low = w$lower.CL,
ci_high = w$upper.CL,
p = w$p.value,
d = w$estimate / sig,
d_low = d + qt(0.025, w$df) * (w$SE / sig),
d_high = d + qt(0.975, w$df) * (w$SE / sig)
)
})
within_total
## # A tibble: 4 × 11
## scale model group contrast estimate ci_low ci_high p d d_low
## <chr> <chr> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BADS_总分 pre-p… 控制… post - … -0.167 -4.00 3.67 9.31e-1 -0.0225 -0.539
## 2 BADS_总分 pre-p… 实验… post - … 5.58 1.74 9.41 5.10e-3 0.752 0.235
## 3 BADS_总分 pre-fu 控制… fu - pre 1.08 -3.18 5.34 6.13e-1 0.131 -0.385
## 4 BADS_总分 pre-fu 实验… fu - pre 7.39 3.13 11.6 9.82e-4 0.897 0.380
## # ℹ 1 more variable: d_high <dbl>