What does the following code do?

source("preproc.R")
set.seed(123)
################################################################################
#1----------- Get empty models to calculate "unadjusted" ICCs with bootstrap CI

iccs <- data.table("Score" = names(score_pairs), "iccMedian" = NaN, 
                   "iccLower" = NaN, "iccUpper" = NaN)

# Loop over scores
for (name in names(score_pairs)) {
  
  empty_model <- lme4::lmer(get(name) ~ 1 + (1|pseudonym) + (1|praxis),
                            data = longdata_full)
  
  boot_lme <- bootMer(empty_model,
                      FUN  = icc_fun_lme,
                      nsim = 5000,
                      type = "parametric")
  
  iccs[Score == name, iccMedian := median(boot_lme$t)]
  iccs[Score == name, iccLower := quantile(boot_lme$t, probs = 0.025)]
  iccs[Score == name, iccUpper := quantile(boot_lme$t, probs = 0.975)]
}

################################################################################
#2----------- Estimate random intercept LMMs for each outcome with bootstrap CI

lme_models <- sapply(names(score_pairs), function(name){
  
  model <- lme4::lmer(get(name) ~ 1 + group + time + group:time + 
                            (1|pseudonym) + (1|praxis),
                          data = longdata_full)
  
}, simplify = FALSE, USE.NAMES = TRUE)

################################################################################
#--------------- Result table with REML-treatment effect estimation and
# Bootstrap CIs plus Model "adjusted" ICC

results <- rbindlist(
  lapply(lme_models, function(model){
    
    fe <- fixef(model)["group:time"]
    ci <- confint(model, method = "boot", parm = "group:time")
    
    boot_lme <- bootMer(model,
                        FUN  = icc_fun_lme,
                        nsim = 5000,
                        type = "parametric")
    
    data.table("TreatmentEffect" = fe,
               "teLower" = ci[1], 
               "teUpper" = ci[2], 
               "iccMedian" = median(boot_lme$t), 
               "iccLower" = quantile(boot_lme$t, probs = 0.025), 
               "iccUpper" = quantile(boot_lme$t, probs = 0.975))
  }
  ), idcol = TRUE)

################################################################################
#-------------- get Data for plot. Population prediction without random effects

plotData <- rbindlist(
  lapply(lme_models, function(model){
    
    # Bootstrap
    merBoot <- bootMer(model, predict_lmer, nsim = 1000, re.form = NA)
    preds <- data.table(expand.grid(pseudonym = NA,
                                    praxis = NA,
                                    group = 0:1,
                                    time = 0:1)
    )
    preds[, median := apply(merBoot$t, 2, median)]
    preds[, conf.low := apply(merBoot$t, 2, quantile, probs = 0.025)]
    preds[, conf.high := apply(merBoot$t, 2, quantile, probs = 0.975)]
    preds[, group := factor(group, levels = c(1, 0), labels = c("IG", "CG"))]
    
    return(preds)
    
  }), idcol = TRUE)


prediction_plot <- ggplot(data = plotData, 
       aes(x = time, y = median, ymin = conf.low, ymax = conf.high, 
           color = group, fill = group)) +
  geom_ribbon(alpha = 1/3, linewidth = 0) +
  geom_line(linewidth = 1) +
  scale_fill_manual(values = c("IG" = "red", "CG" = "blue")) +
  scale_color_manual(values = c("IG" = "red", "CG" = "blue")) +
  scale_x_continuous(breaks = 0:1, labels = c("Pre", "Post"), 
                     expand = c(0.1, 0.1)) +
  scale_y_continuous("Score", sec.axis = dup_axis(name = NULL)) +
  labs(color = "Group", fill = "Group") +  # Proper group label
  theme_minimal() +
  theme(
    axis.title.x = element_text(margin = margin(t = 10)),  
    axis.title.y = element_text(margin = margin(r = 10)),  
    legend.position = "top",  # Move legend for clarity
    panel.spacing = unit(2, "lines")  # Increase space between facets
  ) + 
  ggtitle("Population-mean trajectories",
subtitle = "Bootstrap Predictions with 95% Prediction Interval") +
  facet_wrap(~.id, scales = "free_y", ncol = 2, labeller = as_labeller(c(
    "PHQ9" = "PHQ9",
    "GAD7" = "GAD7",
    "PHQ15" = "PHQ15",
    "PID5BFM" = "PID5BF+M",
    "EBQ" = "EBQ",
    "ERQ_CR" = "ERQ_CR",
    "ERQ_ES" = "ERQ_ES",
    "BEAQ" = "BEAQ"
)))


################################################################################
#-------------------------------- Save Results
write.csv2(iccs,
           file = "Results/Bootstrap_ICCs.csv")

write.csv2(results,
           file = "Results/Bootstrap_CIs.csv")

ggsave("Results/Prediction_Plot.pdf", plot = prediction_plot,
       widt = 10, height = 8)

#-------- check ICC plot

ggplot(iccs) + aes(x = Score, y = iccMedian, ymin = iccLower, ymax = iccUpper) +
  geom_point() + geom_errorbar(aes(colour = "red"))

##########

Also my statistician said that the above script calls the following script. Was does he mean by that:

library(data.table)
library(lme4)
library(ggplot2)
library(marginaleffects)
library(readxl)


# Inserting data
STARC_merged_data <- read_excel("STARC_merged_data.xlsx")

# Ensure data is available
daten <- as.data.table(copy(STARC_merged_data))

# Add numeric group indicator (IG = 1, CG = 0)
daten[, group := 2 - as.numeric(as.factor(gruppe))]


# List of outcome score pairs
score_pairs <- list(
  PHQ9 = c("PHQ9_score_pre", "PHQ9_score_post"),
  GAD7 = c("GAD7_score_pre", "GAD7_score_post"),
  PHQ15 = c("PHQ15_score_pre", "PHQ15_score_post"),
  EBQ = c("EBQ_score_pre", "EBQ_score_post"),
  ERQ_CR = c("ERQ_CR_score_pre", "ERQ_CR_score_post"),
  ERQ_ES = c("ERQ_ES_score_pre", "ERQ_ES_score_post"),
  BEAQ = c("BEAQ_score_pre", "BEAQ_score_post"),
  PID5BFM = c("PID5BFM_score_pre", "PID5BFM_score_post")
)

# Long data set
longdata_full <- copy(daten[, c("group", "praxis", "pseudonym", 
                           unlist(score_pairs)), with = FALSE])

longdata_full <- melt(longdata_full, id.vars = c("pseudonym", "praxis", "group"), 
                 measure.vars = patterns("^PHQ9", "^GAD7", "^PHQ15",
                                         "^EBQ", "^ERQ_CR", "^ERQ_ES", "^BEAQ", 
                                         "^PID5BFM"),
                 variable.name = "time", value.name = names(score_pairs))

longdata_full[, time := as.numeric(time) - 1] # Convert time to numeric

# Define Bootstrap function for icc calculation
icc_fun_lme <- function(fit) {
  vc <- as.data.table(VarCorr(fit))
  vp <- vc[grp == "praxis" & var1 == "(Intercept)", vcov]
  vr <- vc[grp == "Residual", vcov]
  vs <- vc[grp == "pseudonym" & var1 == "(Intercept)", vcov]
  vp / (vp + vr + vs)
}

# Define Bootstrap function for prediction intervals
predict_lmer <- function(fit){
  predict(fit, newdata = expand.grid(pseudonym = NA,
                                     praxis = NA,
                                     group = 0:1,
                                     time = 0:1), 
          re.form = NA)
}