library(tidyverse)
theme_set(theme_bw(base_size = 14))
library(lubridate)
library(tsibble)
library(cmdstanr)
register_knitr_engine(override = FALSE)
library(posterior)
library(bayesplot)
color_scheme_set("blue")
Exam Study Analysis by Three-Level Mixed Effect Model with Two Autoregressive Processes
1 exam_3l-lmm_ZARdARmHm_Seed20250616 results
2 Setup
First, we need to load the necessary packages.
3 Load data
<- read_rds("data/exam_data.rds")
data
<- params$model_name
model_name # model_name <- "exam_3l-lmm_ZHm_Seed20250616"
<- str_detect(model_name, "ARd")
is_ARd <- str_detect(model_name, "Hd")
is_Hd <- str_detect(model_name, "ARm")
is_ARm <- str_detect(model_name, "Hm") is_Hm
4 MCMC check
<- as_cmdstan_fit(list.files(str_glue("stan/draws/{model_name}"),
lmm_fit full.names = TRUE))
# write_rds(lmm_fit, str_glue("stan/{model_name}.rds"))
# lmm_fit <- read_rds("stan/{model_name}.rds")
# lmm_summary <- lmm_fit$summary()
# write_csv(lmm_summary, str_glue("stan/summary/{model_name}_summary.csv"))
<- lmm_fit$draws(format = "df") |> thin_draws(thin = 10)
lmm_draws
<- c(5, 32, 38, 47, 58, 66, 71, 99, 101) selected_subj
Fixed and random effects
mcmc_trace(lmm_draws, pars = "beta")
|> subset_draws(variable = paste0("s[", selected_subj, "]")) |> mcmc_trace() lmm_draws
|> subset_draws(variable = paste0("b[", selected_subj, "]")) |> mcmc_trace() lmm_draws
mcmc_trace(lmm_draws, regex_pars = "^d\\[8,\\d+\\]")
mcmc_trace(lmm_draws, pars = c("psi_s", "psi_b"))
mcmc_trace(lmm_draws, regex_pars = "psi_d")
Autoregressive parameters
<- vector("character")
pars_phi <- vector("character")
pars_tau if (is_ARd) {
<- c(pars_phi, "phi_d")
pars_phi <- c(pars_tau, "tau2_d")
pars_tau
}if (is_ARm) {
<- c(pars_phi, "phi_m")
pars_phi <- c(pars_tau, "tau2_m")
pars_tau
}if (is_ARd || is_ARm) {
mcmc_trace(lmm_draws, pars = pars_phi)
}
if (is_ARd || is_ARm) {
mcmc_trace(lmm_draws, pars = pars_tau)
}
Measurement errors
mcmc_trace(lmm_draws, regex_pars = "sigma_epsilon")
Reliability
mcmc_trace(lmm_draws, pars = c("rel_T"))
5 Fitting results
mcmc_intervals(lmm_draws, regex_pars = "psi_[bsd]")
if (is_ARd || is_ARm) {
mcmc_intervals(lmm_draws, regex_pars = pars_phi)
}
if (is_ARd || is_ARm) {
mcmc_intervals(lmm_draws, regex_pars = pars_tau)
}
mcmc_intervals(lmm_draws, regex_pars = "psi_d")
mcmc_intervals(lmm_draws, regex_pars = "sigma_epsilon")
mcmc_intervals(lmm_draws, regex_pars = "rel_T")
<- lmm_fit$summary("y_hat", mean, median, quantile2)
.y_hat_summary_lmm
<- .y_hat_summary_lmm |>
y_hat_summary_lmm mutate(Indices = str_extract_all(variable, "\\d+"),
Participant = map_dbl(Indices, \(x) as.integer(x[1])),
Day = map_dbl(Indices, \(x) as.integer(x[2])),
Moment = map_dbl(Indices, \(x) as.integer(x[3])))
<- data |>
data_predict_lmm left_join(y_hat_summary_lmm)
|>
data_predict_lmm filter(Participant %in% selected_subj) |>
mutate(Date_time = as_datetime(days(Day) + Time)) |>
ggplot(aes(x = Date_time, y = Neg_aff)) +
geom_line() + geom_point() +
geom_line(aes(y = mean), linetype = "dashed") +
geom_ribbon(aes(ymin = q5, ymax = q95), alpha = 0.25) +
# coord_cartesian(ylim = c(-20, 100)) +
scale_x_datetime(breaks = as_datetime(1:9 * 86400),
labels = paste("Day", 1:9)) +
facet_grid(Participant ~ ., space = "free_y")