### School start times project ### ## Script: Generative model (preregistration) ## # Author: Theresa Lemke # Last edit: 06-07-2026 # Set-up ------------------------------------------------------------ rm(list = ls(all.names = TRUE)) # delete all variables in the environment window try(invisible(dev.off()), silent = TRUE) # remove all previous plots, and some other things cat("\014") # hides all text in the console window library(dplyr) library(tidyr) library(ggplot2) library(brms) library(lme4) library(bayesplot) datadir <- file.path(dirname(dirname(box::file())), "5) Data") outputdir <- file.path(dirname(dirname(box::file())), "6) Analysis/Generative model") # Bayesian Workflow ------------------------------------------------------- # 1. Define the generative model (data-generating process) # 2. Simulate dataset(s) using the generative model # 3. Specify priors # 4. Prior predictive checks # 5. Fit model on simulated data # 6. Posterior predictive checks # 1. Generative Model ----------------------------------------------------- # Model: 3-level multilevel model # Level 1: repeated measurements (t: T1, T2, T3, T4) nested within students - residual error (epsilon) at the observation level (deviation from a student's expected value) # Level 2: students nested within classes, with a student random intercept (student-level deviation from the class-level expected value/between-student differences) # Level 3: classes with SST as class-level predictor and a class random intercept (class-level deviation from grand mean beyond what's explained by SST) # # Rationale for excluding school fixed effects: # - SST varies within each school (half delayed, half not) # - School membership adds little information beyond class-level SST # - Only 3 schools: insufficient to estimate school effects precisely # - Can be added as sensitivity check in secondary analysis # # Model formula: # # Model (1) # weekday sleep duration = alpha [intercept] Grand mean across all observations. # + beta_sst * (MeanSST_c - 8) [fixed effect] MeanSST_c is the weekly average SST for class c. Centered around SST of 08:00. # + beta_occasion * Occasion [fixed effect] Effect of the survey time on sleep. Included primarily to improve imputation of missing # values during the model fitting. # + v_class [class random intercept] Residual class-level differences after accounting for fixed effects. # + u_student [student random intercept] Residual student-level differences after accounting for fixed effects. # + e [measurement residual, i.e., student's sleep fluctuation between measurements] # # # Model (2) - add variance in school start times # weekday sleep duration = alpha # + beta_sst * (MeanSST_c - 8) # + beta_occasion * Occasion # + beta_sst_variation * (SD_SST_c) [fixed effect] SD_SST_c is the standard deviation of the weekly average SST for class c (reflecting the variation in start times during the week). # + v_class # + u_student # + e # # # # Model (3) - add gender effect # weekday sleep duration = alpha # + beta_sst * (MeanSST_c - 8) # + beta_occasion * Occasion # + beta_sst_variation * (SD_SST_c) # + beta_gender * Gender_i [fixed effect] Add fixed-effect for gender # + v_class # + u_student # + e # # Where: # v_class ~ Normal(0, sigma_v) class random intercept # (sigma_v being the standard deviation, distribution of the deviation, i.e., of how much class intercepts vary from overall average) # u_student ~ Normal(0, sigma_u) student random intercept # (sigma_u being the standard deviation, distribution of the deviation, i.e., of how much class intercepts vary from class average) # e ~ Normal(0, sigma_e) measurement-level variance # (sigma_e being the standard deviation of within-student fluctuations across occasions) # # Parameter interpretations: # alpha = expected weekday sleep duration (hours) for a boy (Gender = 0), # in a class with mean SST = 08:00 (centered sst) and SST variability = 0, # for a student and class with random effects equal to zero # (i.e., at the average class and average student) # beta_sst = change in weekday sleep duration (hours) per minute of later SST beyond 08:00, holding all other variables constant # beta_occasion = change in weekday sleep duration (hours) associated with the measurement occasion # beta_sst_variation = expected change in weekday sleep duration (hours) associated with a one-unit increase in the # variability (standard deviation) of start times throughout the week within a class, holding other variables constant. # beta_gender = expected difference in weekday sleep duration between girls (Gender = 1) and boys (Gender = 0) # sigma_v = SD of class-level random intercept (how much average sleep duration varies between classes) # sigma_u = SD of student-level random intercept (how much average sleep duration varies between students in the same class) # sigma_e = SD of measurement-level residual (within-student variability in sleep duration across repeated measurements # (measurement noise and day-to-day fluctuations)) ### 2. SIMULATE DATA ### # ============================================================================== # 2.1 Load observed schedule data - used ONLY to calibrate the simulated SST # distribution (mean, SD, range), not to supply exposure values directly. # This keeps the generative model generalisable beyond this one dataset. # ============================================================================== schedule <- read.csv(file.path(datadir, "School start times.csv"), sep = ";") names(schedule) names(schedule)[names(schedule) == "Grupp..1.Intervention..2..Controle."] <- "Group" days <- c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday") # Transform SST (Monday - Friday) from minutes to hours schedule [days] <- schedule[days]/60 # Dataset with averaged schedules across both schedule periods (i.e., a classes average school start time on Mondays across both schedule periods) schedule_avg <- schedule %>% group_by(Class, School, Group) %>% summarise( across(all_of(days), ~ mean(.x, na.rm = TRUE)), .groups = "drop" ) # Calculate weekly average time per class schedule_avg$sst_week <- rowMeans( schedule_avg [, days], na.rm = TRUE) # standard deviation during the week per class (reflects variation during the week) schedule_avg$sst_sd <- apply( schedule_avg [, days], 1, sd, na.rm = TRUE) # ============================================================================== # 2.2 Extract parameters from observed school start times in the study. # # These are used to simulate SST in the generative model below, # rather than using the fixed 20 real class values. # This allows to increase the generalisability of the generative model and # test the model on a different range of school start times. # # ============================================================================== # Weekly average start time: distribution across the 20 classes total_sst_mean <- mean(schedule_avg$sst_week) # total average weekly start time (across the 20 classes) total_sst_sd <- sd(schedule_avg$sst_week) # total standard deviation of the classes' weekly average start time total_sst_min <- min(schedule_avg$sst_week) # total min weekly average start time total_sst_max <- max(schedule_avg$sst_week) # total max weekly average start time # Variability of school start times during the week: distribution across the 20 classes total_sst_sd_mean <- mean(schedule_avg$sst_sd) # average standard deviation in school start times (the typical within-class SD of start times) total_sst_sd_sd <- sd(schedule_avg$sst_sd) # how much the within-class SD of start times varies between classes total_sst_sd_min <- min(schedule_avg$sst_sd) total_sst_sd_max <- max(schedule_avg$sst_sd) # ============================================================================== # 2.3 Specify true parameter values for the simulated data # ============================================================================== sst_centering <- 8 # centering constant: 08:00 in hours alpha <- 8 # intercept: 8 h - weekday sleep duration (in hours) for boys with 08:00 sst and no within-week variation (sst_sd=0) # (weekday sleep duration for boys was 485 minutes/8.08 hours in Lemke et al. 2023) beta_sst <- 0.5 # 0.5 h - increase in sleep (in hours) per later SST (in hours) (literature: 30min/60min) beta_sst_sd <- -0.25 # 15 minutes decrease in sleep duration per 1-hour increase in SD of sst during the week beta_gender <- -0.4 # -0.4 h - girls sleep 24 min less on weekdays than boys (Lemke et al. 2023) beta_occasion <- 0.1 # Small improvement in sleep over time. ## Variance components ## # Derived from a total SD (sigma_total) + class-level ICC. # # sigma_total is empirically grounded: Lemke et al. (2023) report an overall # weekday sleep duration SD of ~70 min (~68 min boys, ~73 min girls). But this was # a single cross-sectional SD, not a multilevel decomposition. # # icc_class is empirically grounded: Parker et al. (2023) report an average # school-level ICC of 0.04 (95% CI 0.01-0.08) for socio-emotional and health # outcomes in European cluster trials. This is used here as an estimate for the # class-level ICC. # # student_residual is the variance remaining after accounting for the the class-level share (1 - icc_class), # this specifies what proportion is stable between-student differences vs # occasion-to-occasion fluctuation. This remains an assumption and is informed by the general expectation that # repeated self-report measures show more stable between-person variance than pure measurement noise. # # icc_class = proportion of TOTAL variance that is between- # class (Parker et al. 2023, school-level ICC) # student_residual = of the REMAINING variance (1 - icc_class), # proportion that is stable between-student # differences (sigma_u^2); the rest is # occasion-level residual (sigma_e^2) sigma_total <- 70/60 # total SD in hours: ~70 minute weekday sleep duration SD in Lemke et al. (2023) icc_class <- 0.04 # Parker et al. (2023), school-level ICC, here used to approximate class-level ICC student_residual <- 0.65 # assumption sigma_v <- sqrt(icc_class * sigma_total^2) sigma_u <- sqrt(student_residual * (1 - icc_class) * sigma_total^2) sigma_e <- sqrt((1 - student_residual) * (1 - icc_class) * sigma_total^2) # Sanity check: components should add to ~sigma_total (1.166) implied_sigma_total <- sqrt(sigma_v^2 + sigma_u^2 + sigma_e^2) # correct, sigma_v, sigma_u, and sigma_e split the variance correctly according to icc_class and student_share_of_residual # Parameters for data structure (number of schools, classes, time points) n_classes <- 20 n_occasions <- 4 school_names <- c("School 1", "School 2", "School 3") classes_per_school <- c(4, 6, 10) set.seed(2026); students_per_class <- pmax(20, round(rnorm(n_classes, mean = 30, sd = 5))) # simulate a class size for each of the 20 classes, drawn from N(30, 5), but never smaller than 20 dropout_prob <- c(0.10, 0.20, 0.30, 0.40) # ============================================================================== # 2.4 Simulate data (single simulation run) # ============================================================================== sim_dat <- function(seed) { set.seed(seed) # ---- Step A: pre-allocate empty CLASS-level data frame ---- n_class_rows <- n_classes classes <- data.frame( class_id = rep(NA_integer_, n_class_rows), school = rep(NA_character_, n_class_rows), sst_week = rep(NA_real_, n_class_rows), sst_week_centered = rep(NA_real_, n_class_rows), sst_sd = rep(NA_real_, n_class_rows), v_class = rep(NA_real_, n_class_rows) # class random intercept - will be propagated to every student later ) # ---- Fill in class-level columns ---- classes$class_id <- 1:n_classes classes$school <- rep(school_names, times = classes_per_school) # Simulates average weekly start time for each 20 classes # Each simulation run gets a fresh plausible set of 20 class mean start times, anchored in the observed distribution shape but not identical to the observed 20 classes classes$sst_week <- rnorm(n_classes, mean = total_sst_mean, sd = total_sst_sd) # Centers the mean start time around 08:00 classes$sst_week_centered <- classes$sst_week - sst_centering # Simulates day-to-day variability in start times for each class, drawn from a distribution anchored in the observed distribution of start times across the week classes$sst_sd <- abs(rnorm(n_classes, mean = total_sst_sd_mean, sd = total_sst_sd_sd)) classes$v_class <- rnorm(n_classes, 0, sigma_v) # ---- Step B: pre-allocate empty STUDENT-level data frame ---- n_students <- sum(students_per_class) students <- data.frame( student_id = rep(NA_integer_, n_students), class_id = rep(NA_integer_, n_students), gender = rep(NA_integer_, n_students), u_student = rep(NA_real_, n_students) ) # ---- Fill in student-level columns ---- students$student_id <- seq_len(n_students) students$class_id <- rep(1:n_classes, times = students_per_class) students$gender <- rbinom(n_students, 1, prob = 0.50) # 0 = boy, 1 = girl students$u_student <- rnorm(n_students, 0, sigma_u) # student random intercept # Merge students dataframe with classes dataframe # Propagates sst, variation in sst, and class random intercept down to the student level students <- merge(students, classes, by = "class_id", all.x = TRUE) # ---- Step C: pre-allocate empty OCCASION-level (long-format) data frame ---- n_long_rows <- n_students * n_occasions # one row per student per measurement occasion long_data <- data.frame( student_id = rep(NA_integer_, n_long_rows), class_id = rep(NA_integer_, n_long_rows), school = rep(NA_character_, n_long_rows), sst_week = rep(NA_real_, n_long_rows), sst_week_centered = rep(NA_real_, n_long_rows), sst_sd = rep(NA_real_, n_long_rows), v_class = rep(NA_real_, n_long_rows), gender = rep(NA_integer_, n_long_rows), u_student = rep(NA_real_, n_long_rows), occasion = rep(NA_integer_, n_long_rows), mu = rep(NA_real_, n_long_rows), sleep = rep(NA_real_, n_long_rows) ) # ---- Fill in occasion-level columns ---- # Expand each student to n_occasions rows # class-level information (SST, school, v_class) and student-level information (gender, u_student) are propagated to every observation expand_idx <- rep(1:n_students, each = n_occasions) long_data$student_id <- students$student_id[expand_idx] long_data$class_id <- students$class_id[expand_idx] long_data$school <- students$school[expand_idx] long_data$sst_week <- students$sst_week[expand_idx] long_data$sst_week_centered <- students$sst_week_centered[expand_idx] long_data$sst_sd <- students$sst_sd[expand_idx] long_data$v_class <- students$v_class[expand_idx] long_data$gender <- students$gender[expand_idx] long_data$u_student <- students$u_student[expand_idx] long_data$occasion <- rep(1:n_occasions, times = n_students) # occasion number # ---- Linear predictor (for Model 3 - full model) ---- # Model 1 and 2 are simplified fits to this same data long_data$mu <- alpha + beta_sst * long_data$sst_week_centered + # mean SST effect beta_sst_sd * long_data$sst_sd + # SST variability effect beta_gender * long_data$gender + # gender effect beta_occasion * long_data$occasion + # occasion effect long_data$v_class + # class random intercept long_data$u_student # student random intercept # ---- Observed sleep duration (hours) with occasion-level noise ---- # Draws the actual observed sleep duration for each row from a normal distribution centred at that row's mu with standard deviation sigma_e. long_data$sleep_raw <- rnorm(n_long_rows, mean = long_data$mu, sd = sigma_e) # Limit to plausible range of sleep duration: 3 to 12 hours long_data$sleep_raw <- pmax(3, pmin(12, long_data$sleep_raw)) # ---- Missing data (increasing dropout over occasions) ---- long_data$missing <- rbinom(n_long_rows, 1, prob = dropout_prob[long_data$occasion]) long_data$sleep <- ifelse(long_data$missing == 1, NA, long_data$sleep_raw) # Return both levels for inspection list(classes = classes, long_data = long_data) } # Simulation sim <- sim_dat(seed = 2026) classes_sim <- sim$classes # simulated dataset with class-level data only long_data_sim <- sim$long_data # simulated dataset that will be used for the following prior model # Save simulated data # saveRDS(long_data_sim, file.path(outputdir, "simulated_data.rds")) # write.csv(long_data_sim, file.path(outputdir, "simulated_data.csv"), row.names = FALSE) # ============================================================================== # 3. Priors # ============================================================================== # alpha ~ Normal(8, 1) Expected ~8h sleep at 08:00 SST. # Weekday sleep duration for boys was 485 minutes/8.08 hours in Lemke et al. 2023. # SD = 1h - allows ~2 SD range (95% of the data) of 6-10h. # beta_sst_week_centered ~ Normal(0.5, 0.25) Centered at estimate from the literature (30 min per # 60 min delay = 0.5). Central 95% include 0.01 to 0.99 hours of extra sleep # per hour of later school start time. The bottom 2.5% include values close to zero and slightly # negative (later school start decreasing sleep). # beta_sst_sd ~ Normal(-0.25, 0.15) Expected effect of 1-hour SD in start times on sleep duration: 15 minutes decrease # Central 95% interval include -33 to +3. 4.8% prior probability of a positive effect. # Assumption, no literature about this. But circadian disruption seems likely. # beta_gender ~ Normal(-0.4, 0.25) Centered at literature estimate (girls sleep ~24 min / 0.4h less than boys). # 95% of values will fall within 2 SD, i.e., 50 minutes (girls sleep 54 min less to girls sleep 6 min more) # beta_occasion ~ Normal(0 , 0.5) Expected effect of measurement occasion on sleep. Mean=0 (no assumed effect of change over time). # SD = 0.5h allows up to 1h (2SD) change in sleep duration per occasion in either direction. # Included primarily to improve imputation of missing values. # sigma_v ~ Normal(0.5, 0.17) class-level variance. mean=30 min, SD=10 minutes # sigma_u ~ Normal(1.0, 0.33) student-level variance. mean=1 hour, SD=20 minutes # sigma_e ~ Normal(1.0, 0.33) measurement-level residual. mean=1 hour, SD=20 minutes # Check what priors brms requires and what brms uses internally as grouping factor names for random effects # get_prior(sleep ~ sst_week_centered + occasion + sst_sd + gender + (1 | class_id) + (1 | class_id:student_id), data = long_data1) # --- Shared priors used by ALL three models --- prior_base <- c( prior(normal(8, 1), class = "Intercept"), prior(normal(0.5, 0.25), class = "b", coef = "sst_week_centered"), prior(normal(0, 0.5), class = "b", coef = "occasion"), prior(normal(0.5, 0.17), class = "sd", group = "class_id"), prior(normal(1.0, 0.33), class = "sd", group = "class_id:student_id"), prior(normal(1.0, 0.33), class = "sigma") ) # --- Model 1: SST mean + occasion only --- prior_m1 <- prior_base # --- Model 2: + SST variability --- prior_m2 <- c( prior_base, prior(normal(-0.25, 0.15), class = "b", coef = "sst_sd") ) # --- Model 3: + gender --- prior_m3 <- c( prior_m2, prior(normal(-0.4, 0.25), class = "b", coef = "gender") ) # ============================================================================== # 4. Prior predictive checks # Are the priors reasonable? Do the priors include reasonable data? # ============================================================================== # Sample the prior predictive distribution (sample_prior = "only") # The prior-only model didn't converge when doing the imputation, so the prior predictive checks # or done using complete-cases only. Then later when fitting the model on the data, the model # includes the imputation. # MODEL 1 # Each chain contributes 2000 (iter - warmup) to the total samples (8000) fit_m1_prior <- brm( sleep ~ sst_week_centered + occasion + (1 | class_id) + (1 | class_id:student_id) , data = long_data_sim, prior = prior_m1, sample_prior = "only", chains = 4, iter = 4000, warmup = 1000, # MCMC starts a chain with some initial parameters, runs it for a while to get into our region of interest (warm up) seed = 2026, silent = 2 ) # MCMC diagnostic plots plot(fit_m1_prior) # chains are converged and mixing well # Plot the predicted y density (sleep duration predicted from the prior predictive distribution) across the data's observed range pp_check(fit_m1_prior, ndraws = 100) # MODEL 2 # add SST variability to the model fit_m2_prior <- brm( sleep ~ sst_week_centered + occasion + sst_sd + (1 | class_id) + (1 | class_id:student_id) , data = long_data_sim, prior = prior_m2, sample_prior = "only", chains = 4, iter = 4000, warmup = 1000, seed = 2026, silent = 2 ) # MCMC diagnostic plots plot(fit_m2_prior) # Plot the predicted y density (sleep duration predicted from the prior predictive distribution) across the data's observed range pp_check(fit_m2_prior, ndraws = 100) # MODEL 3 fit_m3_prior <- brm( sleep ~ sst_week_centered + occasion + sst_sd + gender + (1 | class_id) + (1 | class_id:student_id) , data = long_data_sim, prior = prior_m3, sample_prior = "only", chains = 4, iter = 4000, warmup = 1000, seed = 2026, silent = 2 ) # MCMC diagnostic plots plot(fit_m3_prior) # Plot the predicted y density (sleep duration predicted from the prior predictive distribution) across the data's observed range pp_check(fit_m3_prior, ndraws = 100) # Save final prior model (model 3) # Extract and save all prior samples variables(fit_m3_prior) samples_prior <- as_draws_df(fit_m3_prior) # saveRDS(samples_prior, file.path(outputdir, "prior_samples.rds")) # Save prior model # saveRDS(fit_m3_prior, file.path(outputdir, "prior_model.rds")) # Save the standardized model (correlation coefficients) # Standardised model # Compute SDs from the simulated dataset (complete cases) sd_sst_week <- sd(long_data_sim$sst_week_centered, na.rm = TRUE) sd_sst_sd <- sd(long_data_sim$sst_sd, na.rm = TRUE) sd_gender <- sd(long_data_sim$gender, na.rm = TRUE) sd_occasion <- sd(long_data_sim$occasion, na.rm = TRUE) sd_sleep <- sd(long_data_sim$sleep, na.rm = TRUE) # Standardize by both the predictor SD and the outcome SD # Expresses each coefficient as "standard deviations of sleep per standard deviation of the predictor" — making all predictors directly comparable regardless of their original units std_samples_prior <- samples_prior |> mutate( std_b_sst_week_centered = b_sst_week_centered * (sd_sst_week / sd_sleep), std_b_sst_sd = b_sst_sd * (sd_sst_sd / sd_sleep), std_b_gender = b_gender * (sd_gender / sd_sleep), std_b_occasion = b_occasion * (sd_occasion / sd_sleep) ) head(std_samples_prior) # Save # saveRDS(std_samples_prior, file.path(outputdir, "prior_samples_standardized.rds")) # ============================================================================== # 5. Fit models (incl imputation) # ============================================================================== # Imputation of missing data: see https://paulbuerkner.com/brms/articles/brms_missings.html # Fit M1 fit_m1 <- brm(formula = bf( sleep | mi() ~ sst_week_centered + occasion + (1 | class_id) + (1 | class_id:student_id)) , data = long_data_sim, prior = prior_m1, family = gaussian(), sample_prior = "yes", # to allow for plotting prior vs posterior from the fitted model later on chains = 4, iter = 10000, warmup = 3000, cores = 4, control = list(adapt_delta = .99), seed = 2026, silent = 2 ) # Convergence diagnostics plot(fit_m1) # Fit M2 fit_m2 <- brm(formula = bf( sleep | mi() ~ sst_week_centered + occasion + sst_sd + (1 | class_id) + (1 | class_id:student_id)) , data = long_data_sim, prior = prior_m2, family = gaussian(), sample_prior = "yes", chains = 4, iter = 10000, warmup = 3000, cores = 4, control = list(adapt_delta = .99), seed = 2026, silent = 2 ) # Convergence diagnostics plot(fit_m2) # Fit M3 fit_m3 <- brm(formula = bf( sleep | mi() ~ sst_week_centered + occasion + sst_sd + gender + (1 | class_id) + (1 | class_id:student_id)) , data = long_data_sim, prior = prior_m3, family = gaussian(), sample_prior = "yes", chains = 4, iter = 10000, warmup = 3000, cores = 4, control = list(adapt_delta = .99), seed = 2026, silent = 2 ) # Convergence diagnostics plot(fit_m3) # Save final posterior model (model 3) # Extract and save posterior samples samples_posterior <- as_draws_df(fit_m3) # saveRDS(samples_posterior, file.path(outputdir, "posterior_samples.rds")) # Save posterior model # saveRDS(fit_m3, file.path(outputdir, "posterior_model.rds")) # Save the standardized model (correlation coefficients) # Standardised model # Compute SDs from the simulated dataset (complete cases) # sd_sst_week <- sd(long_data_sim$sst_week_centered, na.rm = TRUE) # sd_sst_sd <- sd(long_data_sim$sst_sd, na.rm = TRUE) # sd_gender <- sd(long_data_sim$gender, na.rm = TRUE) # sd_occasion <- sd(long_data_sim$occasion, na.rm = TRUE) # sd_sleep <- sd(long_data_sim$sleep, na.rm = TRUE) # Standardize by both the predictor SD and the outcome SD # Expresses each coefficient as "standard deviations of sleep per standard deviation of the predictor" — making all predictors directly comparable regardless of their original units std_samples_posterior <- samples_posterior |> mutate( std_b_sst_week_centered = b_sst_week_centered * (sd_sst_week / sd_sleep), std_b_sst_sd = b_sst_sd * (sd_sst_sd / sd_sleep), std_b_gender = b_gender * (sd_gender / sd_sleep), std_b_occasion = b_occasion * (sd_occasion / sd_sleep) ) head(std_samples_posterior) # Save # saveRDS(std_samples_posterior, file.path(outputdir, "posterior_samples_standardized.rds")) # ============================================================================== # 6. Posterior predictive checks and posterior distribution # ============================================================================== # MODEL 3 summary(fit_m3) # 50 sleep duration density curves drawn from the posterior predictive, overlaid on the observed # sleep duration density (dark blue) pp_check(fit_m3, type = "dens_overlay", ndraws = 100) # Distribution of the posterior intercept mcmc_areas(fit_m3, pars = c("b_Intercept"), prob = 0.95) # Distribution of posterior parameter estimate for school start times mcmc_areas(fit_m3, pars = "b_sst_week_centered", prob = 0.95) # Posterior distribution for variation in school start times throughout the week mcmc_areas(fit_m3, pars = "b_sst_sd", prob = 0.95) # Posterior distribution for gender mcmc_areas(fit_m3, pars = "b_gender", prob = 0.95) # Create a prior posterior plot # samples_prior # samples from the prior distribution # plot the prior distribution ggplot() + geom_density(data = samples_prior, aes(x = b_sst_week_centered)) + labs(title = "Prior Distribution") + theme_classic() # plot posterior distribution # samples_posterior # extracted posterior samples ggplot() + geom_density(data = samples_posterior, aes(x = b_sst_week_centered), color = "blue") + labs(title = "Posterior Distribution") + theme_classic() # show prior and posterior distributions in same plot ggplot() + geom_density(data = samples_prior, aes(x = b_sst_week_centered)) + geom_density(data = samples_posterior, aes(x = b_sst_week_centered), color = "blue") + geom_vline(xintercept = beta_sst, color = "red", linetype = "dashed", linewidth = 1) + labs(title = "Prior Distribution vs. Posterior Distribution") + theme_classic() # Plot the fitted SST effect (marginal effects) # Shows the posterior mean sleep duration, 95% credible interval, effect of SST ce <- conditional_effects( fit_m3, effects = "sst_week_centered", re_formula = NA # include all group-level effects/random effects ) plot(ce)