Assessment of Change in CBT Knowledge

Repeated Measures Analysis in Diploma Program Students

1 Introduction

1.1 Study Description

This document presents the complete analysis of a study evaluating changes in Cognitive Behavioral Therapy (CBT) knowledge among students in a diploma program. Participants were assessed at three time points:

• Time 1: Program start
• Time 2: Mid-term evaluation
• Time 3: Program completion

1.2 Hypotheses

A significant improvement in CBT knowledge over time is expected, manifested in:

1. Increase in the total score on the instrument
2. Improvement in each of the dimensions evaluated
3. Possible differential effect based on prior clinical experience

1.3 Report Structure

This report includes:

• Complete descriptive statistics
• Repeated measures analysis (ANOVA)
• Post-hoc comparions
• Statistical power analysis
• Visualizations
• Interpretation of results

---

2 Data Preparation

2.1 Load Packages

# Data manipulation
library(readxl)
library(tidyverse)
library(here)

# Statistical analysis
library(lme4)
library(lmerTest)
library(emmeans)
library(car)
library(effectsize)
library(pwr)

# APA Tables
library(apaTables)
library(flextable)
library(officer)
library(gt)

# Visualization
library(ggplot2)
library(patchwork)
library(scales)

cat("✓ Packages loaded successfully\n")

## ✓ Packages loaded successfully

2.2 Configure Plot Theme

# Custom theme for plots
tema_custom <- theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14, hjust = 0),
    plot.subtitle = element_text(size = 11, color = "gray30"),
    axis.title = element_text(face = "bold"),
    legend.position = "bottom",
    panel.grid.minor = element_blank(),
    strip.text = element_text(face = "bold", size = 11)
  )

theme_set(tema_custom)

2.3 Define Instrument Dimensions

# Dimensions and their items
dimensiones <- list(
  "Mastering_CBT" = c(40, 36, 30, 32, 31, 17, 26, 39, 44, 37, 50, 45, 22, 33, 23, 41, 34),
  "Generic_Competencies" = c(21, 46, 51, 16, 6, 25, 20, 27, 47, 48, 15, 49, 1, 2, 13),
  "CBT_Techniques" = c(9, 14, 18, 12, 3, 38, 5),
  "Formulating_Problems" = c(19, 4, 29, 10, 42)  # Reverse-scored items
)

# Reverse-scored items
items_reversos <- c(19, 4, 29, 10, 42)

# Maximum scores
puntajes_max <- data.frame(
  Dimension = c("Total", names(dimensiones)),
  N_Items = c(44, sapply(dimensiones, length)),
  Min = c(44, sapply(dimensiones, length)),
  Max = c(220, sapply(dimensiones, length) * 5)
)

flextable(puntajes_max) %>%
  set_caption("Table 1. Instrument Structure and Scoring") %>%
  theme_booktabs() %>%
  autofit()

Table 1. Instrument Structure and Scoring

Dimension	N_Items	Min	Max
Total	44	44	220
Mastering_CBT	17	17	85
Generic_Competencies	15	15	75
CBT_Techniques	7	7	35
Formulating_Problems	5	5	25

2.4 Load Data

# ADJUST THIS PATH according to your file location
ruta_archivo <- "C:/Users/marco/Downloads/Datos_corregidos.xlsx"

# Read data
datos <- read_excel(ruta_archivo)

# Rename main columns
colnames(datos)[1:3] <- c("ID", "Time", "Clinical_Exp")

# Show structure
cat("Data dimensions:", nrow(datos), "rows ×", ncol(datos), "columns\n")

## Data dimensions: 93 rows × 47 columns

cat("Unique participants:", length(unique(datos$ID)), "\n")

## Unique participants: 31

cat("Measurements per participant:", table(table(datos$ID)), "\n")

## Measurements per participant: 31

2.5 Process Data

# 1. Reverse score items
for(item in items_reversos) {
  col_name <- as.character(item)
  if(col_name %in% colnames(datos)) {
    datos[[col_name]] <- 6 - datos[[col_name]]
  }
}

cat("✓ Reverse-scored items processed:", paste(items_reversos, collapse = ", "), "\n")

## ✓ Reverse-scored items processed: 19, 4, 29, 10, 42

# 2. Calculate dimension scores
for(dim in names(dimensiones)) {
  items <- as.character(dimensiones[[dim]])
  items_disponibles <- items[items %in% colnames(datos)]
  datos[[dim]] <- rowSums(datos[, items_disponibles, drop = FALSE], na.rm = TRUE)
}

cat("✓ Dimension scores calculated\n")

## ✓ Dimension scores calculated

# 3. Calculate total score
item_cols <- as.character(1:51)
item_cols_disponibles <- item_cols[item_cols %in% colnames(datos)]
datos$Total <- rowSums(datos[, item_cols_disponibles, drop = FALSE], na.rm = TRUE)

cat("✓ Total score calculated\n")

## ✓ Total score calculated

# 4. Factorize variables
datos <- datos %>%
  mutate(
    ID = as.factor(ID),
    Time = factor(Time, levels = c(1, 2, 3),
                      labels = c("Time 1", "Time 2", "Time 3")),
    Clinical_Exp = factor(Clinical_Exp, 
                         levels = c(0, 1),
                         labels = c("No clinical experience", "Clinical experience"))
  )

cat("✓ Variables factorized\n")

## ✓ Variables factorized

---

3 Descriptive Statistics

3.1 Total Score

# Calculate descriptive statistics
desc_total <- datos %>%
  group_by(Time, Clinical_Exp) %>%
  summarise(
    N = n(),
    M = mean(Total),
    SD = sd(Total),
    SE = SD / sqrt(N),
    CI_lower = M - 1.96 * SE,
    CI_upper = M + 1.96 * SE,
    Min = min(Total),
    Max = max(Total),
    .groups = "drop"
  )

# APA Table using flextable
desc_total %>%
  mutate(
    `M (SD)` = sprintf("%.2f (%.2f)", M, SD),
    `95% IC` = sprintf("[%.2f, %.2f]", CI_lower, CI_upper),
    Range = sprintf("[%.0f, %.0f]", Min, Max)
  ) %>%
  select(Time, Clinical_Exp, N, `M (SD)`, `95% IC`, Range) %>%
  flextable() %>%
  set_caption("Table 2. Descriptive Statistics for Total Score by Time and Group") %>%
  merge_v(j = "Time") %>%
  theme_booktabs() %>%
  autofit() %>%
  align(align = "center", part = "all") %>%
  align(j = 1:2, align = "left", part = "body")

Table 2. Descriptive Statistics for Total Score by Time and Group

Time	Clinical_Exp	N	M (SD)	95% IC	Range
Time 1	No clinical experience	9	123.67 (27.10)	[105.96, 141.37]	[73, 148]
	Clinical experience	22	158.32 (26.34)	[147.31, 169.33]	[89, 207]
Time 2	No clinical experience	9	147.11 (37.38)	[122.69, 171.54]	[78, 214]
	Clinical experience	22	172.59 (23.10)	[162.94, 182.24]	[122, 213]
Time 3	No clinical experience	9	163.56 (26.70)	[146.11, 181.00]	[135, 220]
	Clinical experience	22	179.59 (19.41)	[171.48, 187.70]	[144, 215]

Note. M = Mean; SD = Standard Deviation; CI = 95% Confidence Interval. Total score: theoretical range = 44-220.

3.2 Statistics by Dimension

# Prepare long data
datos_largo <- datos %>%
  select(ID, Time, Clinical_Exp, all_of(names(dimensiones)), Total) %>%
  pivot_longer(
    cols = c(Total, all_of(names(dimensiones))),
    names_to = "Dimension",
    values_to = "Score"
  ) %>%
  mutate(
    Dimension = factor(Dimension,
                      levels = c("Total", names(dimensiones)),
                      labels = c("Total Score", "Mastering CBT theory and skills",
                                "Generic psychotherapy competencies", "Using CBT techniques properly",
                                "Formulating clients' problems and plans"))
  )

# Calculate descriptives
desc_dim <- datos_largo %>%
  group_by(Dimension, Time) %>%
  summarise(
    N = n(),
    M = mean(Score),
    SD = sd(Score),
    .groups = "drop"
  ) %>%
  pivot_wider(
    names_from = Time,
    values_from = c(M, SD),
    names_glue = "{Time}_{.value}"
  )

# Format as APA table
desc_dim %>%
  mutate(
    `Time 1` = sprintf("%.2f (%.2f)", `Time 1_M`, `Time 1_SD`),
    `Time 2` = sprintf("%.2f (%.2f)", `Time 2_M`, `Time 2_SD`),
    `Time 3` = sprintf("%.2f (%.2f)", `Time 3_M`, `Time 3_SD`)
  ) %>%
  select(Dimension, `Time 1`, `Time 2`, `Time 3`) %>%
  flextable() %>%
  set_caption("Table 3. Means and Standard Deviations by Dimension and Time") %>%
  theme_booktabs() %>%
  autofit() %>%
  align(align = "center", part = "all") %>%
  align(j = 1, align = "left", part = "body") %>%
  add_footer_lines("Note. Values are presented as M (SD).")

Table 3. Means and Standard Deviations by Dimension and Time

Dimension	Time 1	Time 2	Time 3
Total Score	148.26 (30.61)	165.19 (29.74)	174.94 (22.55)
Mastering CBT theory and skills	47.39 (15.51)	58.39 (14.12)	64.26 (10.22)
Generic psychotherapy competencies	56.81 (11.11)	59.84 (10.46)	62.23 (8.44)
Using CBT techniques properly	23.68 (5.36)	25.90 (4.92)	26.90 (3.87)
Formulating clients' problems and plans	20.39 (2.60)	21.06 (2.93)	21.55 (2.78)
Note. Values are presented as M (SD).

---

4 Descriptive Visualizations

4.1 Evolution of Total Score

ggplot(datos, aes(x = Time, y = Total, group = ID, color = Clinical_Exp)) +
  geom_line(alpha = 0.3, linewidth = 0.5) +
  geom_point(alpha = 0.3, size = 2) +
  stat_summary(aes(group = Clinical_Exp), fun = mean, geom = "line", linewidth = 2) +
  stat_summary(aes(group = Clinical_Exp), fun = mean, geom = "point", 
               size = 4, shape = 21, fill = "white", stroke = 1.5) +
  stat_summary(aes(group = Clinical_Exp), fun.data = mean_se, 
               geom = "errorbar", width = 0.2, linewidth = 1) +
  geom_hline(yintercept = 220, linetype = "dashed", color = "gray40", alpha = 0.6) +
  scale_color_manual(values = c("#E69F00", "#56B4E9")) +
  scale_y_continuous(limits = c(44, 220), breaks = seq(50, 220, 25)) +
  facet_wrap(~ Clinical_Exp) +
  labs(
    title = "Individual Trajectories of Total Score",
    subtitle = "Thin lines: individual participants | Thick line: group mean | Dashed line: theoretical maximum",
    x = "Assessment Time",
    y = "Total Score (range: 44-220)",
    color = "Clinical Experience"
  ) +
  theme(legend.position = "none")

Figure 1. Evolution of total score over time by clinical experience group.

4.2 Evolution by Dimension

# Prepare data with limits
limites_ejes <- data.frame(
  Dimension = c("Total Score", "Mastering CBT theory and skills",
                "Generic psychotherapy competencies", "Using CBT techniques properly",
                "Formulating clients' problems and plans"),
  Min = c(44, 17, 15, 7, 5),
  Max = c(220, 85, 75, 35, 25)
)

stats_dim <- datos_largo %>%
  group_by(Dimension, Time, Clinical_Exp) %>%
  summarise(
    Mean = mean(Score),
    SE = sd(Score) / sqrt(n()),
    .groups = "drop"
  ) %>%
  left_join(limites_ejes, by = "Dimension")

ggplot(stats_dim, aes(x = Time, y = Mean, 
                      color = Clinical_Exp, group = Clinical_Exp)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), 
                width = 0.2, linewidth = 0.8) +
  geom_hline(aes(yintercept = Max), linetype = "dashed", 
             color = "gray40", alpha = 0.5) +
  geom_blank(aes(y = Min)) +
  geom_blank(aes(y = Max)) +
  facet_wrap(~ Dimension, scales = "free_y", ncol = 2) +
  scale_color_manual(values = c("#E69F00", "#56B4E9")) +
  labs(
    title = "Evolution of CBT Knowledge by Dimension",
    subtitle = "Dashed line indicates maximum possible score",
    x = "Assessment Time",
    y = "Average Score",
    color = "Clinical Experience",
    caption = "Error bars represent standard errors of the mean"
  )

Figure 2. Evolution of knowledge by dimension.

---

5 Inferential Analysis

5.1 Repeated Measures ANOVA - Total Score

# Mixed model
modelo_total <- lmer(Total ~ Time * Clinical_Exp + (1|ID),
                     data = datos)

# ANOVA Table
anova_resultado <- Anova(modelo_total, type = 3)

# Format as table
anova_resultado %>%
  as.data.frame() %>%
  rownames_to_column("Effect") %>%
  mutate(
    Sig = case_when(
      `Pr(>Chisq)` < 0.001 ~ "***",
      `Pr(>Chisq)` < 0.01 ~ "**",
      `Pr(>Chisq)` < 0.05 ~ "*",
      TRUE ~ "ns"
    )
  ) %>%
  flextable() %>%
  set_caption("Table 4. Repeated Measures ANOVA for Total Score") %>%
  colformat_double(j = c("Chisq", "Pr(>Chisq)"), digits = 4) %>%
  theme_booktabs() %>%
  autofit() %>%
  add_footer_lines("Note. ***p < .001, **p < .01, *p < .05")

Table 4. Repeated Measures ANOVA for Total Score

Effect	Chisq	Df	Pr(>Chisq)	Sig
(Intercept)	212.1210	1	0.0000	***
Time	39.1513	2	0.0000	***
Clinical_Exp	11.8191	1	0.0006	***
Time:Clinical_Exp	5.9907	2	0.0500	ns
Note. ***p < .001, **p < .01, *p < .05

5.1.1 Effect Sizes

# Calculate partial eta squared
eta <- eta_squared(modelo_total, partial = TRUE)

eta %>%
  as.data.frame() %>%
  flextable() %>%
  set_caption("Table 5. Effect Sizes (Partial Eta Squared)") %>%
  colformat_double(j = c("Eta2_partial", "CI_low", "CI_high"), digits = 3) %>%
  theme_booktabs() %>%
  autofit() %>%
  add_footer_lines("Note. Small: η² < .06, Medium: η² < .14, Large: η² ≥ .14 (Cohen, 1988)")

Table 5. Effect Sizes (Partial Eta Squared)

Parameter	Eta2_partial	CI	CI_low	CI_high
Time	0.532	0.95	0.379	1.000
Clinical_Exp	0.213	0.95	0.035	1.000
Time:Clinical_Exp	0.094	0.95	0.000	1.000
Note. Small: η² < .06, Medium: η² < .14, Large: η² ≥ .14 (Cohen, 1988)

5.2 Post-Hoc Comparisons

# Estimated marginal means
emm_medicion <- emmeans(modelo_total, ~ Time)

# Pairwise comparisons
comparaciones <- pairs(emm_medicion, adjust = "bonferroni")

# Format results
comp_df <- summary(comparaciones) %>%
  as.data.frame() %>%
  mutate(
    Sig = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ "ns"
    ),
    `d de Cohen` = estimate / sqrt(2 * SE^2)  # Approximation
  )

comp_df %>%
  select(contrast, estimate, SE, df, t.ratio, p.value, `d de Cohen`, Sig) %>%
  flextable() %>%
  set_caption("Table 6. Pairwise Comparisons Between Assessment Times") %>%
  colformat_double(j = c("estimate", "SE", "t.ratio", "p.value", "d de Cohen"), 
                  digits = 3) %>%
  theme_booktabs() %>%
  autofit() %>%
  add_footer_lines("Note. Bonferroni correction applied (α = .0167 for three comparisons).")

Table 6. Pairwise Comparisons Between Assessment Times

contrast	estimate	SE	df	t.ratio	p.value	d de Cohen	Sig
Time 1 - Time 2	-18.859	3.803	58	-4.959	0.000	-3.506	***
Time 1 - Time 3	-30.581	3.803	58	-8.041	0.000	-5.686	***
Time 2 - Time 3	-11.722	3.803	58	-3.082	0.009	-2.180	**
Note. Bonferroni correction applied (α = .0167 for three comparisons).

# Calculate OBSERVED MEANS
medias_reales <- datos %>%
  group_by(Time) %>%
  summarise(
    N = n(),
    `Observed M` = mean(Total),
    `SD` = sd(Total),
    `SE Real` = SD / sqrt(N),
    .groups = "drop"
  )

# Get estimated marginal means from model
medias_estimadas <- summary(emm_medicion) %>%
  as.data.frame() %>%
  select(Time, `Estimated M` = emmean, `Estimated SE` = SE, 
         `CI_Lower` = lower.CL, `CI_Upper` = upper.CL)

# Combine
medias_combinadas <- medias_reales %>%
  left_join(medias_estimadas, by = "Time") %>%
  mutate(
    Difference = `Estimated M` - `Observed M`
  )

# Table with both means
medias_combinadas %>%
  mutate(
    `Observed Mean` = sprintf("%.2f (%.2f)", `Observed M`, SD),
    `Estimated Mean` = sprintf("%.2f (%.2f)", `Estimated M`, `Estimated SE`),
    Dif = sprintf("%+.2f", Difference)
  ) %>%
  select(Time, N, `Observed Mean`, `Estimated Mean`, Dif) %>%
  flextable() %>%
  set_caption("Table 7. Observed Means vs Estimated Marginal Means") %>%
  theme_booktabs() %>%
  autofit() %>%
  add_footer_lines(c(
    "Note. Observed Mean = M (SD) directly observed in the data.",
    "Estimated Mean = Estimated marginal mean (SE) from mixed model.",
    "Dif = Difference between estimated and observed mean. Minimal differences indicate good model fit."
  ))

Table 7. Observed Means vs Estimated Marginal Means

Time	N	Observed Mean	Estimated Mean	Dif
Time 1	31	148.26 (30.61)	140.99 (5.04)	-7.27
Time 2	31	165.19 (29.74)	159.85 (5.04)	-5.34
Time 3	31	174.94 (22.55)	171.57 (5.04)	-3.36
Note. Observed Mean = M (SD) directly observed in the data.
Estimated Mean = Estimated marginal mean (SE) from mixed model.
Dif = Difference between estimated and observed mean. Minimal differences indicate good model fit.

Explanation of Observed vs Estimated Marginal Means:

• Observed Mean: Simple arithmetic average of observed data (sum/n)
• Estimated Mean Marginal: Mean adjusted by the statistical model, controlling for random effects and other variables (clinical experience)
• Why do they differ?: EMMs “average” across experience groups and consider the repeated measures structure
• For reporting: Use actual means in descriptive tables, and EMMs in statistical comparisons
  (consistent with ANOVA p-values)

5.3 ANOVA by Dimension

# Function to analyze one dimension
analizar_dimension <- function(datos, variable) {
  modelo <- lmer(as.formula(paste(variable, "~ Time * Clinical_Exp + (1|ID)")),
                data = datos)
  anova_res <- Anova(modelo, type = 3)
  eta_res <- eta_squared(modelo, partial = TRUE)
  
  list(
    anova = anova_res,
    eta = eta_res,
    modelo = modelo
  )
}

# Analyze all dimensions
resultados_dim <- lapply(names(dimensiones), function(dim) {
  analizar_dimension(datos, dim)
})
names(resultados_dim) <- names(dimensiones)

# Create summary table
resumen_dim <- data.frame()

for(dim in names(dimensiones)) {
  anova_res <- resultados_dim[[dim]]$anova
  eta_res <- resultados_dim[[dim]]$eta
  
  # Extract values for Time effect
  f_val <- anova_res["Time", "Chisq"]
  p_val <- anova_res["Time", "Pr(>Chisq)"]
  eta_val <- eta_res$Eta2_partial[eta_res$Parameter == "Time"]
  
  resumen_dim <- rbind(resumen_dim, data.frame(
    Dimension = gsub("_", " ", dim),
    F = f_val,
    p = p_val,
    `η²` = eta_val,
    check.names = FALSE
  ))
}

resumen_dim %>%
  mutate(
    Sig = case_when(
      p < 0.001 ~ "***",
      p < 0.01 ~ "**",
      p < 0.05 ~ "*",
      TRUE ~ "ns"
    )
  ) %>%
  flextable() %>%
  set_caption("Table 8. ANOVA by Dimension (Effect of Time)") %>%
  colformat_double(j = c("F", "p", "η²"), digits = 3) %>%
  theme_booktabs() %>%
  autofit() %>%
  add_footer_lines("Note. ***p < .001, **p < .01, *p < .05")

Table 8. ANOVA by Dimension (Effect of Time)

Dimension	F	p	η²	Sig
Mastering CBT	46.648	0.000	0.603	***
Generic Competencies	25.049	0.000	0.325	***
CBT Techniques	14.250	0.001	0.275	***
Formulating Problems	0.082	0.960	0.060	ns
Note. ***p < .001, **p < .01, *p < .05

5.4 Post-Hoc Comparisons by Dimension

Pairwise comparisons between assessment times for each dimension are presented below.

5.4.1 Mastering CBT theory and skills

# Marginal means and comparisons
emm_mastering <- emmeans(resultados_dim$Mastering_CBT$modelo, ~ Time)
comp_mastering <- pairs(emm_mastering, adjust = "bonferroni")

# Comparison table
summary(comp_mastering) %>%
  as.data.frame() %>%
  mutate(
    Sig = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ "ns"
    ),
    `d de Cohen` = estimate / sqrt(2 * SE^2)
  ) %>%
  select(contrast, estimate, SE, df, t.ratio, p.value, `d de Cohen`, Sig) %>%
  flextable() %>%
  set_caption("Table 9. Post-Hoc Comparisons: Mastering CBT theory and skills") %>%
  colformat_double(j = c("estimate", "SE", "t.ratio", "p.value", "d de Cohen"), 
                  digits = 3) %>%
  theme_booktabs() %>%
  autofit() %>%
  add_footer_lines("Note. Bonferroni correction (α = .0167). d = Cohen's effect size.")

Table 9. Post-Hoc Comparisons: Mastering CBT theory and skills

contrast	estimate	SE	df	t.ratio	p.value	d de Cohen	Sig
Time 1 - Time 2	-12.051	2.013	58	-5.985	0.000	-4.232	***
Time 1 - Time 3	-18.649	2.013	58	-9.262	0.000	-6.549	***
Time 2 - Time 3	-6.598	2.013	58	-3.277	0.005	-2.317	**
Note. Bonferroni correction (α = .0167). d = Cohen's effect size.

# Marginal means
summary(emm_mastering) %>%
  as.data.frame() %>%
  mutate(Type = "Estimated (Model)") %>%
  select(Time, Type, M = emmean, SE, CI_Lower = lower.CL, CI_Upper = upper.CL) %>%
  bind_rows(
    datos %>%
      group_by(Time) %>%
      summarise(
        Type = "Observed",
        M = mean(Mastering_CBT),
        SE = sd(Mastering_CBT) / sqrt(n()),
        CI_Lower = M - 1.96 * SE,
        CI_Upper = M + 1.96 * SE,
        .groups = "drop"
      )
  ) %>%
  flextable() %>%
  set_caption("Means by Time: Mastering CBT theory and skills (Observed vs Estimated)") %>%
  colformat_double(j = c("M", "SE", "CI_Lower", "CI_Upper"), digits = 2) %>%
  theme_booktabs() %>%
  merge_v(j = "Time") %>%
  autofit()

Means by Time: Mastering CBT theory and skills (Observed vs Estimated)

Time	Type	M	SE	CI_Lower	CI_Upper
Time 1	Estimated (Model)	44.38	2.53	39.29	49.48
Time 2	Estimated (Model)	56.43	2.53	51.34	61.53
Time 3	Estimated (Model)	63.03	2.53	57.93	68.13
Time 1	Observed	47.39	2.79	41.93	52.85
Time 2	Observed	58.39	2.54	53.42	63.36
Time 3	Observed	64.26	1.84	60.66	67.86

5.4.2 Generic Psychotherapy Competencies

# Marginal means and comparisons
emm_generic <- emmeans(resultados_dim$Generic_Competencies$modelo, ~ Time)
comp_generic <- pairs(emm_generic, adjust = "bonferroni")

# Comparison table
summary(comp_generic) %>%
  as.data.frame() %>%
  mutate(
    Sig = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ "ns"
    ),
    `d de Cohen` = estimate / sqrt(2 * SE^2)
  ) %>%
  select(contrast, estimate, SE, df, t.ratio, p.value, `d de Cohen`, Sig) %>%
  flextable() %>%
  set_caption("Table 10. Post-Hoc Comparisons: Generic Psychotherapy Competencies") %>%
  colformat_double(j = c("estimate", "SE", "t.ratio", "p.value", "d de Cohen"), 
                  digits = 3) %>%
  theme_booktabs() %>%
  autofit() %>%
  add_footer_lines("Note. Bonferroni correction (α = .0167). d = Cohen's effect size.")

Table 10. Post-Hoc Comparisons: Generic Psychotherapy Competencies

contrast	estimate	SE	df	t.ratio	p.value	d de Cohen	Sig
Time 1 - Time 2	-3.876	1.370	58	-2.828	0.019	-2.000	*
Time 1 - Time 3	-7.232	1.370	58	-5.277	0.000	-3.732	***
Time 2 - Time 3	-3.356	1.370	58	-2.449	0.052	-1.732	ns
Note. Bonferroni correction (α = .0167). d = Cohen's effect size.

# Marginal means
summary(emm_generic) %>%
  as.data.frame() %>%
  mutate(Type = "Estimated (Model)") %>%
  select(Time, Type, M = emmean, SE, CI_Lower = lower.CL, CI_Upper = upper.CL) %>%
  bind_rows(
    datos %>%
      group_by(Time) %>%
      summarise(
        Type = "Observed",
        M = mean(Generic_Competencies),
        SE = sd(Generic_Competencies) / sqrt(n()),
        CI_Lower = M - 1.96 * SE,
        CI_Upper = M + 1.96 * SE,
        .groups = "drop"
      )
  ) %>%
  flextable() %>%
  set_caption("Means by Time: Generic psychotherapy competencies (Observed vs Estimated)") %>%
  colformat_double(j = c("M", "SE", "CI_Lower", "CI_Upper"), digits = 2) %>%
  theme_booktabs() %>%
  merge_v(j = "Time") %>%
  autofit()

Means by Time: Generic psychotherapy competencies (Observed vs Estimated)

Time	Type	M	SE	CI_Lower	CI_Upper
Time 1	Estimated (Model)	53.97	1.80	50.34	57.61
Time 2	Estimated (Model)	57.85	1.80	54.22	61.49
Time 3	Estimated (Model)	61.21	1.80	57.57	64.84
Time 1	Observed	56.81	1.99	52.90	60.72
Time 2	Observed	59.84	1.88	56.16	63.52
Time 3	Observed	62.23	1.52	59.26	65.20

5.4.3 Using CBT techniques properly

# Marginal means and comparisons
emm_techniques <- emmeans(resultados_dim$CBT_Techniques$modelo, ~ Time)
comp_techniques <- pairs(emm_techniques, adjust = "bonferroni")

# Comparison table
summary(comp_techniques) %>%
  as.data.frame() %>%
  mutate(
    Sig = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ "ns"
    ),
    `d de Cohen` = estimate / sqrt(2 * SE^2)
  ) %>%
  select(contrast, estimate, SE, df, t.ratio, p.value, `d de Cohen`, Sig) %>%
  flextable() %>%
  set_caption("Table 11. Post-Hoc Comparisons: Using CBT techniques properly") %>%
  colformat_double(j = c("estimate", "SE", "t.ratio", "p.value", "d de Cohen"), 
                  digits = 3) %>%
  theme_booktabs() %>%
  autofit() %>%
  add_footer_lines("Note. Bonferroni correction (α = .0167). d = Cohen's effect size.")

Table 11. Post-Hoc Comparisons: Using CBT techniques properly

contrast	estimate	SE	df	t.ratio	p.value	d de Cohen	Sig
Time 1 - Time 2	-2.422	0.823	58	-2.944	0.014	-2.082	*
Time 1 - Time 3	-3.816	0.823	58	-4.638	0.000	-3.280	***
Time 2 - Time 3	-1.394	0.823	58	-1.694	0.287	-1.198	ns
Note. Bonferroni correction (α = .0167). d = Cohen's effect size.

# Marginal means
summary(emm_techniques) %>%
  as.data.frame() %>%
  mutate(Type = "Estimated (Model)") %>%
  select(Time, Type, M = emmean, SE, CI_Lower = lower.CL, CI_Upper = upper.CL) %>%
  bind_rows(
    datos %>%
      group_by(Time) %>%
      summarise(
        Type = "Observed",
        M = mean(CBT_Techniques),
        SE = sd(CBT_Techniques) / sqrt(n()),
        CI_Lower = M - 1.96 * SE,
        CI_Upper = M + 1.96 * SE,
        .groups = "drop"
      )
  ) %>%
  flextable() %>%
  set_caption("Means by Time: Using CBT techniques properly (Observed vs Estimated)") %>%
  colformat_double(j = c("M", "SE", "CI_Lower", "CI_Upper"), digits = 2) %>%
  theme_booktabs() %>%
  merge_v(j = "Time") %>%
  autofit()

Means by Time: Using CBT techniques properly (Observed vs Estimated)

Time	Type	M	SE	CI_Lower	CI_Upper
Time 1	Estimated (Model)	22.59	0.88	20.83	24.35
Time 2	Estimated (Model)	25.01	0.88	23.25	26.77
Time 3	Estimated (Model)	26.41	0.88	24.65	28.17
Time 1	Observed	23.68	0.96	21.79	25.56
Time 2	Observed	25.90	0.88	24.17	27.64
Time 3	Observed	26.90	0.69	25.54	28.26

5.4.4 Formulating Clients’ Problems and Plans

# Marginal means and comparisons
emm_formulating <- emmeans(resultados_dim$Formulating_Problems$modelo, ~ Time)
comp_formulating <- pairs(emm_formulating, adjust = "bonferroni")

# Comparison table
summary(comp_formulating) %>%
  as.data.frame() %>%
  mutate(
    Sig = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ "ns"
    ),
    `d de Cohen` = estimate / sqrt(2 * SE^2)
  ) %>%
  select(contrast, estimate, SE, df, t.ratio, p.value, `d de Cohen`, Sig) %>%
  flextable() %>%
  set_caption("Table 12. Post-Hoc Comparisons: Formulating clients' problems and plans and Plans") %>%
  colformat_double(j = c("estimate", "SE", "t.ratio", "p.value", "d de Cohen"), 
                  digits = 3) %>%
  theme_booktabs() %>%
  autofit() %>%
  add_footer_lines("Note. Bonferroni correction (α = .0167). d = Cohen's effect size.")

Table 12. Post-Hoc Comparisons: Formulating clients' problems and plans and Plans

contrast	estimate	SE	df	t.ratio	p.value	d de Cohen	Sig
Time 1 - Time 2	-0.510	0.462	58	-1.105	0.821	-0.781	ns
Time 1 - Time 3	-0.884	0.462	58	-1.915	0.181	-1.354	ns
Time 2 - Time 3	-0.374	0.462	58	-0.810	1.000	-0.572	ns
Note. Bonferroni correction (α = .0167). d = Cohen's effect size.

# Marginal means
summary(emm_formulating) %>%
  as.data.frame() %>%
  mutate(Type = "Estimated (Model)") %>%
  select(Time, Type, M = emmean, SE, CI_Lower = lower.CL, CI_Upper = upper.CL) %>%
  bind_rows(
    datos %>%
      group_by(Time) %>%
      summarise(
        Type = "Observed",
        M = mean(Formulating_Problems),
        SE = sd(Formulating_Problems) / sqrt(n()),
        CI_Lower = M - 1.96 * SE,
        CI_Upper = M + 1.96 * SE,
        .groups = "drop"
      )
  ) %>%
  flextable() %>%
  set_caption("Means by Time: Formulating clients' problems and plans (Observed vs Estimated)") %>%
  colformat_double(j = c("M", "SE", "CI_Lower", "CI_Upper"), digits = 2) %>%
  theme_booktabs() %>%
  merge_v(j = "Time") %>%
  autofit()

Means by Time: Formulating clients' problems and plans (Observed vs Estimated)

Time	Type	M	SE	CI_Lower	CI_Upper
Time 1	Estimated (Model)	20.04	0.51	19.02	21.07
Time 2	Estimated (Model)	20.55	0.51	19.53	21.58
Time 3	Estimated (Model)	20.93	0.51	19.90	21.95
Time 1	Observed	20.39	0.47	19.47	21.30
Time 2	Observed	21.06	0.53	20.03	22.10
Time 3	Observed	21.55	0.50	20.57	22.53

5.5 Visual Summary of Changes by Dimension

# Prepare change data
cambios_data <- data.frame()

for(dim in names(dimensiones)) {
  emm <- emmeans(resultados_dim[[dim]]$modelo, ~ Time)
  comp <- summary(pairs(emm, adjust = "bonferroni"))
  
  for(i in 1:nrow(comp)) {
    cambios_data <- rbind(cambios_data, data.frame(
      Dimension = gsub("_", " ", dim),
      Comparacion = as.character(comp$contrast[i]),
      Estimado = comp$estimate[i],
      SE = comp$SE[i],
      p_valor = comp$p.value[i],
      Sig = ifelse(comp$p.value[i] < 0.0167, "Sig", "ns")
    ))
  }
}

# Plot
ggplot(cambios_data, aes(x = Comparacion, y = Estimado, fill = Sig)) +
  geom_col(color = "black", linewidth = 0.5) +
  geom_errorbar(aes(ymin = Estimado - SE, ymax = Estimado + SE),
                width = 0.3, linewidth = 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
  facet_wrap(~ Dimension, scales = "free_y", ncol = 2) +
  scale_fill_manual(values = c("Sig" = "#56B4E9", "ns" = "#E0E0E0"),
                    labels = c("Sig" = "Significant (p < .0167)", 
                              "ns" = "Not significant")) +
  coord_flip() +
  labs(
    title = "Magnitude of Changes by Dimension",
    subtitle = "Differences between assessment times with Bonferroni correction",
    x = "",
    y = "Estimated Difference (points)",
    fill = ""
  ) +
  theme(legend.position = "bottom")

Figure 3. Changes between time by dimension.

5.6 Simple Effects Analysis (Interaction)

When the Time × Clinical Experience interaction is significant, we analyze simple effects to understand how the pattern of change differs between groups.

# Function to analyze simple effects
analizar_efectos_simples <- function(modelo, nombre_dim) {
  
  # Interaction test
  anova_res <- Anova(modelo, type = 3)
  p_interaccion <- anova_res["Time:Clinical_Exp", "Pr(>Chisq)"]
  
  # Only proceed if interaction is significant or marginally significant
  if(p_interaccion < 0.10) {
    
    cat("\n", rep("-", 70), "\n", sep = "")
    cat(nombre_dim, "\n")
    cat("Interaction p =", sprintf("%.4f", p_interaccion), "\n")
    cat(rep("-", 70), "\n", sep = "")
    
    # Simple effects: Change in each group
    emm_interaccion <- emmeans(modelo, ~ Time | Clinical_Exp)
    comp_interaccion <- pairs(emm_interaccion, adjust = "bonferroni")
    
    resultado <- summary(comp_interaccion) %>%
      as.data.frame() %>%
      mutate(
        Sig = case_when(
          p.value < 0.001 ~ "***",
          p.value < 0.01 ~ "**",
          p.value < 0.05 ~ "*",
          p.value < 0.10 ~ "†",
          TRUE ~ "ns"
        ),
        `d de Cohen` = estimate / sqrt(2 * SE^2),
        Dimension = nombre_dim
      ) %>%
      select(Dimension, Clinical_Exp, contrast, estimate, SE, t.ratio, p.value, `d de Cohen`, Sig)
    
    return(resultado)
  } else {
    return(NULL)
  }
}

# Analyze all dimensions
efectos_simples_todas <- list()

for(dim in names(dimensiones)) {
  nombre_display <- gsub("_", " ", dim)
  resultado <- analizar_efectos_simples(resultados_dim[[dim]]$modelo, nombre_display)
  
  if(!is.null(resultado)) {
    efectos_simples_todas[[dim]] <- resultado
  }
}

## 
## ----------------------------------------------------------------------
## Generic Competencies 
## Interaction p = 0.0068 
## ----------------------------------------------------------------------

# Combine todos los resultados
if(length(efectos_simples_todas) > 0) {
  
  efectos_simples_df <- bind_rows(efectos_simples_todas)
  
  efectos_simples_df %>%
    flextable() %>%
    set_caption("Table 13. Effects Simples: Cambio en Cada Grupo de Clinical Experience") %>%
    colformat_double(j = c("estimate", "SE", "t.ratio", "p.value", "d de Cohen"), 
                    digits = 3) %>%
    theme_booktabs() %>%
    autofit() %>%
    merge_v(j = c("Dimension", "Clinical_Exp")) %>%
    add_footer_lines("Note. Bonferroni correction. ***p < .001, **p < .01, *p < .05, †p < .10")
  
  cat("\n\nSignificant or marginal interactions were found in:\n")
  cat(paste("  -", unique(efectos_simples_df$Dimension), collapse = "\n"))
  
} else {
  cat("\nNo significant interactions were found (p < .10) in any dimension.\n")
  cat("This indicates that the pattern of change is similar between groups with and without clinical experience.\n")
}

## 
## 
## Significant or marginal interactions were found in:
##   - Generic Competencies

# Only generate if there are significant interactions
if(length(efectos_simples_todas) > 0) {
  
  # Prepare data for dimensions with interaction
  dims_con_interaccion <- names(efectos_simples_todas)
  
  datos_interaccion <- datos %>%
    select(ID, Time, Clinical_Exp, all_of(dims_con_interaccion)) %>%
    pivot_longer(
      cols = all_of(dims_con_interaccion),
      names_to = "Dimension",
      values_to = "Score"
    ) %>%
    mutate(Dimension = gsub("_", " ", Dimension))
  
  # Calculate means and SE
  stats_interaccion <- datos_interaccion %>%
    group_by(Dimension, Time, Clinical_Exp) %>%
    summarise(
      Mean = mean(Score),
      SE = sd(Score) / sqrt(n()),
      .groups = "drop"
    )
  
  ggplot(stats_interaccion, aes(x = Time, y = Mean, 
                                 color = Clinical_Exp, group = Clinical_Exp)) +
    geom_line(linewidth = 1.5) +
    geom_point(size = 4, shape = 21, fill = "white", stroke = 1.5) +
    geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE),
                  width = 0.2, linewidth = 1) +
    facet_wrap(~ Dimension, scales = "free_y", ncol = 2) +
    scale_color_manual(values = c("#E69F00", "#56B4E9")) +
    labs(
      title = "Time × Clinical Experience Interaction",
      subtitle = "Dimensions with significant interaction (p < .10)",
      x = "Assessment Time",
      y = "Average Score",
      color = "Clinical Experience"
    ) +
    theme(legend.position = "bottom")
}

Figure 4. Change profiles by clinical experience group.

---

6 Individual Items Analysis

6.1 Assessment of Ceiling Effects by Item

This analysis examines each item to identify possible ceiling effects, which is important for interpreting results and planning future program improvements.

# Prepare item data
items_largo <- datos %>%
  select(ID, Time, Clinical_Exp, any_of(as.character(1:51))) %>%
  pivot_longer(
    cols = any_of(as.character(1:51)),
    names_to = "Item",
    values_to = "Respuesta"
  ) %>%
  mutate(
    Item_Num = as.numeric(Item),
    # Classify by dimension
    Dimension = case_when(
      Item_Num %in% dimensiones$Mastering_CBT ~ "Mastering CBT theory and skills",
      Item_Num %in% dimensiones$Generic_Competencies ~ "Generic psychotherapy competencies",
      Item_Num %in% dimensiones$CBT_Techniques ~ "Using CBT techniques properly",
      Item_Num %in% dimensiones$Formulating_Problems ~ "Formulating clients' problems and plans",
      TRUE ~ "Otro"
    )
  ) %>%
  filter(Dimension != "Otro")

# Calculate statistics by item and time
stats_items <- items_largo %>%
  group_by(Dimension, Item_Num, Time) %>%
  summarise(
    N = n(),
    Mean = mean(Respuesta, na.rm = TRUE),
    SD = sd(Respuesta, na.rm = TRUE),
    Prop_Max = mean(Respuesta == 5, na.rm = TRUE) * 100,  # % in maximum
    .groups = "drop"
  )

# Calculate change from T1 to T3
cambio_items <- stats_items %>%
  filter(Time %in% c("Time 1", "Time 3")) %>%
  select(Dimension, Item_Num, Time, Mean, Prop_Max) %>%
  pivot_wider(
    names_from = Time,
    values_from = c(Mean, Prop_Max)
  ) %>%
  mutate(
    Cambio_Mean = `Mean_Time 3` - `Mean_Time 1`,
    Cambio_Prop_Max = `Prop_Max_Time 3` - `Prop_Max_Time 1`
  )

# Identify items with ceiling effect
items_techo <- cambio_items %>%
  filter(`Prop_Max_Time 3` > 50 | `Mean_Time 3` > 4.5) %>%
  arrange(desc(`Prop_Max_Time 3`))

cat("\nItems with possible ceiling effect (>50% at maximum or M>4.5 at T3):\n")

## 
## Items with possible ceiling effect (>50% at maximum or M>4.5 at T3):

cat("Total:", nrow(items_techo), "items\n\n")

## Total: 7 items

6.1.1 Table of Items with Ceiling Effect

if(nrow(items_techo) > 0) {
  items_techo %>%
    mutate(
      `% Maximum T1` = sprintf("%.1f%%", `Prop_Max_Time 1`),
      `% Maximum T3` = sprintf("%.1f%%", `Prop_Max_Time 3`),
      `M T1` = sprintf("%.2f", `Mean_Time 1`),
      `M T3` = sprintf("%.2f", `Mean_Time 3`),
      `Mean Change` = sprintf("%+.2f", Cambio_Mean)
    ) %>%
    select(Dimension, Item = Item_Num, `M T1`, `M T3`, `Mean Change`, 
           `% Maximum T1`, `% Maximum T3`) %>%
    flextable() %>%
    set_caption("Table 14. Items with Possible Ceiling Effect at Time 3") %>%
    theme_booktabs() %>%
    autofit() %>%
    merge_v(j = "Dimension") %>%
    add_footer_lines("Note. M = Mean on 1-5 scale. % Maximum = percentage of responses with value 5.")
} else {
  cat("No items with marked ceiling effect were identified.\n")
}

Table 14. Items with Possible Ceiling Effect at Time 3

Dimension	Item	M T1	M T3	Mean Change	% Maximum T1	% Maximum T3
Formulating clients' problems and plans	4	4.39	4.65	+0.26	51.6%	74.2%
Generic psychotherapy competencies	20	4.06	4.61	+0.55	48.4%	71.0%
Formulating clients' problems and plans	29	4.39	4.55	+0.16	51.6%	67.7%
Generic psychotherapy competencies	46	4.06	4.45	+0.39	45.2%	58.1%
	6	4.29	4.42	+0.13	58.1%	54.8%
Formulating clients' problems and plans	42	3.77	4.32	+0.55	19.4%	51.6%
Generic psychotherapy competencies	2	4.23	4.42	+0.19	29.0%	51.6%
Note. M = Mean on 1-5 scale. % Maximum = percentage of responses with value 5.

6.1.2 Items with Largest and Smallest Change

# Top 10 items with largest improvement
top_mejora <- cambio_items %>%
  arrange(desc(Cambio_Mean)) %>%
  head(10)

# Top 10 items with smallest improvement (or deterioration)
top_menor <- cambio_items %>%
  arrange(Cambio_Mean) %>%
  head(10)

# Table of largest improvement
top_mejora %>%
  mutate(
    `M T1` = sprintf("%.2f", `Mean_Time 1`),
    `M T3` = sprintf("%.2f", `Mean_Time 3`),
    `Cambio` = sprintf("%+.2f", Cambio_Mean)
  ) %>%
  select(Dimension, Item = Item_Num, `M T1`, `M T3`, `Cambio`) %>%
  flextable() %>%
  set_caption("Table 15. Items with Largest Improvement (Top 10)") %>%
  theme_booktabs() %>%
  autofit() %>%
  merge_v(j = "Dimension") %>%
  add_footer_lines("Note. Change = M T3 - M T1. Positive values indicate improvement.")

Table 15. Items with Largest Improvement (Top 10)

Dimension	Item	M T1	M T3	Cambio
Mastering CBT theory and skills	34	2.16	4.26	+2.10
	26	2.29	3.84	+1.55
	31	2.55	3.71	+1.16
	30	3.00	4.13	+1.13
	39	2.58	3.71	+1.13
	37	2.74	3.84	+1.10
	17	2.10	3.19	+1.10
	40	2.61	3.71	+1.10
	50	2.97	4.06	+1.10
Using CBT techniques properly	18	2.71	3.77	+1.06
Note. Change = M T3 - M T1. Positive values indicate improvement.

# Table of smallest improvement
top_menor %>%
  mutate(
    `M T1` = sprintf("%.2f", `Mean_Time 1`),
    `M T3` = sprintf("%.2f", `Mean_Time 3`),
    `Cambio` = sprintf("%+.2f", Cambio_Mean)
  ) %>%
  select(Dimension, Item = Item_Num, `M T1`, `M T3`, `Cambio`) %>%
  flextable() %>%
  set_caption("Table 16. Items with Smallest Improvement (Bottom 10)") %>%
  theme_booktabs() %>%
  autofit() %>%
  merge_v(j = "Dimension") %>%
  add_footer_lines("Note. Items showing less change or even deterioration.")

Table 16. Items with Smallest Improvement (Bottom 10)

Dimension	Item	M T1	M T3	Cambio
Formulating clients' problems and plans	10	3.94	3.94	+0.00
Generic psychotherapy competencies	1	4.00	4.13	+0.13
	6	4.29	4.42	+0.13
Using CBT techniques properly	5	3.97	4.10	+0.13
Generic psychotherapy competencies	15	3.81	3.97	+0.16
Formulating clients' problems and plans	29	4.39	4.55	+0.16
	19	3.90	4.10	+0.19
Generic psychotherapy competencies	2	4.23	4.42	+0.19
	51	3.81	4.00	+0.19
Using CBT techniques properly	3	3.45	3.65	+0.19
Note. Items showing less change or even deterioration.

6.1.3 Visualization of Items by Dimension

ggplot(stats_items, aes(x = Time, y = Mean, group = Item_Num)) +
  geom_line(alpha = 0.3, color = "gray60") +
  geom_point(alpha = 0.3, size = 1, color = "gray60") +
  stat_summary(aes(group = 1), fun = mean, geom = "line", 
               color = "#D55E00", linewidth = 2) +
  stat_summary(aes(group = 1), fun = mean, geom = "point",
               color = "#D55E00", size = 4, shape = 21, fill = "white", stroke = 1.5) +
  geom_hline(yintercept = 4.5, linetype = "dashed", color = "red", alpha = 0.5) +
  geom_hline(yintercept = 5, linetype = "solid", color = "red", alpha = 0.3) +
  facet_wrap(~ Dimension, ncol = 2) +
  scale_y_continuous(limits = c(1, 5), breaks = 1:5) +
  labs(
    title = "Evolución de Items Individuales by Dimension",
    subtitle = "Gray lines: individual items | Orange line: average | Red lines: thresholds for ceiling effect",
    x = "Assessment Time",
    y = "Average Score (scale 1-5)",
    caption = "Dashed line: M = 4.5 (possible ceiling) | Solid line: M = 5.0 (absolute ceiling)"
  )

Figure 5. Distribution of item means by Dimension and Time.

# Calculate % at maximum for T3
prop_max_t3 <- stats_items %>%
  filter(Time == "Time 3") %>%
  mutate(Item = paste0("Item ", Item_Num))

ggplot(prop_max_t3, aes(x = reorder(Item, Prop_Max), y = Prop_Max, fill = Dimension)) +
  geom_col() +
  geom_hline(yintercept = 50, linetype = "dashed", color = "red", linewidth = 1) +
  coord_flip() +
  scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#CC79A7")) +
  labs(
    title = "Percentage of Participants with Maximum Score by Item (Time 3)",
    subtitle = "Red line indicates 50% threshold (possible ceiling effect)",
    x = "",
    y = "% of Responses with Value 5 (maximum)",
    fill = "Dimension"
  ) +
  theme(
    axis.text.y = element_text(size = 6),
    legend.position = "bottom"
  )

Figure 6. Percentage of responses at the maximum per item.

6.1.4 Interpretation of Item Analysis

# Statistical summary
n_items_techo <- nrow(items_techo)
n_items_total <- length(unique(items_largo$Item_Num))
prop_techo <- 100 * n_items_techo / n_items_total

cat("\n**Item Analysis Summary:**\n\n")

Item Analysis Summary:

cat(sprintf("- **Total items analyzed:** %d\n", n_items_total))

• Total items analyzed: 44

cat(sprintf("- **Items with possible ceiling effect:** %d (%.1f%%)\n", n_items_techo, prop_techo))

• Items with possible ceiling effect: 7 (15.9%)

if(n_items_techo > 0) {
  # By dimension
  techo_por_dim <- items_techo %>%
    group_by(Dimension) %>%
    summarise(N = n(), .groups = "drop")
  
  cat("\n- **Distribution of ceiling effects by dimension:**\n")
  for(i in 1:nrow(techo_por_dim)) {
    cat(sprintf("  - %s: %d items\n", techo_por_dim$Dimension[i], techo_por_dim$N[i]))
  }
  
  cat("\n**Implicaciones:**\n\n")
  cat("Items with ceiling effects suggest areas where students achieved\n")
  cat("a high level of competence by the end of the program. These items could:\n\n")
  cat("1. Reflect fundamental aspects well mastered by all students\n")
  cat("2. Be too easy for this level of training\n")
  cat("3. Require revision to include more advanced content\n\n")
} else {
  cat("\n**Implicaciones:**\n\n")
  cat("The absence of widespread ceiling effects indicates that the instrument has\n")
  cat("an adequate range of difficulty and can detect improvements even at levels\n")
  cat("of advanced competence.\n\n")
}

• Distribution of ceiling effects by dimension:
  • Formulating clients’ problems and plans: 3 items
  • Generic psychotherapy competencies: 4 items

Implicaciones:

Items with ceiling effects suggest areas where students achieved a high level of competence by the end of the program. These items could:

1. Reflect fundamental aspects well mastered by all students
2. Be too easy for this level of training
3. Require revision to include more advanced content

# Items with little or no change
items_sin_cambio <- cambio_items %>%
  filter(abs(Cambio_Mean) < 0.2)

if(nrow(items_sin_cambio) > 0) {
  cat(sprintf("\n- **Items with minimal change** (|Δ| < 0.2): %d\n", nrow(items_sin_cambio)))
  cat("\nThese items might represent:\n")
  cat("- Content not emphasized in the program\n")
  cat("- Prior knowledge already consolidated\n")
  cat("- Aspects requiring more time or practice\n")
}

• Items with minimal change (|Δ| < 0.2): 11

These items might represent: - Content not emphasized in the program - Prior knowledge already consolidated - Aspects requiring more time or practice

---

7 Statistical Power Analysis

7.1 A Priori Power (MDES)

# Study parameters
n_participants <- length(unique(datos$ID))
n_mediciones <- 3
alpha <- 0.05

# Calculate correlation between Times
datos_wide <- datos %>%
  select(ID, Time, Total) %>%
  pivot_wider(names_from = Time, values_from = Total)

# Column names will be "Time 1", "Time 2", "Time 3"
# Create version with simple names
datos_wide_simple <- datos_wide %>%
  rename(T1 = `Time 1`, T2 = `Time 2`, T3 = `Time 3`)

cor_mediciones <- cor(datos_wide_simple[, c("T1", "T2", "T3")], use = "complete.obs")
cor_promedio <- mean(cor_mediciones[lower.tri(cor_mediciones)])

# Find minimum detectable f
f_min <- uniroot(
  function(f) {
    pwr.anova.test(k = n_mediciones, n = n_participants, 
                   f = f, sig.level = alpha)$power - 0.80
  },
  interval = c(0.01, 2),
  tol = 0.001
)$root

eta_min <- f_min^2 / (1 + f_min^2)

# For paired comparisons
d_min <- uniroot(
  function(d) {
    pwr.t.test(n = n_participants, d = d, sig.level = alpha/3, 
               type = "paired")$power - 0.80
  },
  interval = c(0.01, 3),
  tol = 0.001
)$root

A Priori Power Analysis:

With N = 31 participants, k = 3 measures, and α = 0.05, the study has 80% of power to detect:

• Repeated measures ANOVA: Cohen’s f ≥ 0.327 (equivalent to η² ≥ 0.097)
• Paired comparions (with Bonferroni correction): Cohen’s d ≥ 0.611

These values represent the Minimun Detectable Effect Size (MDES). Effects smaller than these thresholds may not be detectable with this sample.

# Create power table for different effect sizes
tamanos_efecto <- c(0.10, 0.25, 0.40, 0.50, f_min, 0.80)

poder_tabla <- data.frame()
for(f in tamanos_efecto) {
  poder <- pwr.anova.test(k = n_mediciones, n = n_participants, 
                         f = f, sig.level = alpha)$power
  
  poder_tabla <- rbind(poder_tabla, data.frame(
    `Cohen's f` = f,
    `η²` = f^2 / (1 + f^2),
    Power = poder,
    check.names = FALSE
  ))
}

poder_tabla %>%
  mutate(
    Category = case_when(
      `Cohen's f` < 0.25 ~ "Small",
      `Cohen's f` < 0.40 ~ "Medium",
      TRUE ~ "Large"
    ),
    `Power (%)` = sprintf("%.1f%%", Power * 100)
  ) %>%
  arrange(`Cohen's f`) %>%
  flextable() %>%
  set_caption("Table 9. Statistical Power for Different Effect Sizes") %>%
  colformat_double(j = c("Cohen's f", "η²", "Power"), digits = 3) %>%
  theme_booktabs() %>%
  autofit() %>%
  bold(i = ~ round(`Cohen's f`, 3) == round(f_min, 3)) %>%
  add_footer_lines(sprintf("Note. MDES (80%% power) highlighted in bold: f = %.3f", f_min))

Table 9. Statistical Power for Different Effect Sizes

Cohen's f	η²	Power	Category	Power (%)
0.100	0.010	0.124	Small	12.4%
0.250	0.059	0.555	Medium	55.5%
0.327	0.097	0.800	Medium	80.0%
0.400	0.138	0.935	Large	93.5%
0.500	0.200	0.993	Large	99.3%
0.800	0.390	1.000	Large	100.0%
Note. MDES (80% power) highlighted in bold: f = 0.327

---

8 Interpretation of Results

8.1 Main Findings

1. Significant Improvement Over Time: A highly significant main significant effect of time on CBT knowledge (p < .001), with an effect size large effect size (η² = 0.532).

2. Continuous Progress: All pairwise comparisons were statistically significant even after Bonferroni correction, indicating continuous improvement in each phase of the diploma.

3. Magnitude of Change: Students improved on average by 26.68 points (IC 95%: [18.25, 35.11]) from program start to completion.

4. Absence of Ceiling Effect: The average final score (174.9) represents approximately 79.5% of the theoretical maximum, indicating there was no ceiling effect that limited the detection of improvements.

8.2 Implications

Results demonstrate that the diploma program was effective in improving CBT knowledge. The observed effect size far exceeds the MDES, providing confidence that these findings represent real and educationally significant effects.

---

9 Session Information

sessionInfo()

## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Chile.utf8  LC_CTYPE=Spanish_Chile.utf8   
## [3] LC_MONETARY=Spanish_Chile.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Chile.utf8    
## 
## time zone: America/Santiago
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] scales_1.4.0     patchwork_1.3.2  gt_1.3.0         officer_0.7.3   
##  [5] flextable_0.9.10 apaTables_2.0.8  pwr_1.3-0        effectsize_1.0.1
##  [9] car_3.1-3        carData_3.0-5    emmeans_1.11.2-8 lmerTest_3.1-3  
## [13] lme4_1.1-37      Matrix_1.7-3     here_1.0.2       lubridate_1.9.4 
## [17] forcats_1.0.1    stringr_1.6.0    dplyr_1.1.4      purrr_1.2.0     
## [21] readr_2.1.6      tidyr_1.3.1      tibble_3.3.0     ggplot2_4.0.1   
## [25] tidyverse_2.0.0  readxl_1.4.5    
## 
## loaded via a namespace (and not attached):
##  [1] Rdpack_2.6.4            sandwich_3.1-1          rlang_1.1.6            
##  [4] magrittr_2.0.3          multcomp_1.4-28         compiler_4.5.1         
##  [7] systemfonts_1.3.1       vctrs_0.6.5             pkgconfig_2.0.3        
## [10] fastmap_1.2.0           backports_1.5.0         labeling_0.4.3         
## [13] rmarkdown_2.29          tzdb_0.5.0              nloptr_2.2.1           
## [16] ragg_1.5.0              xfun_0.53               cachem_1.1.0           
## [19] jsonlite_2.0.0          uuid_1.2-1              parallel_4.5.1         
## [22] broom_1.0.10            R6_2.6.1                bslib_0.9.0            
## [25] stringi_1.8.7           RColorBrewer_1.1-3      boot_1.3-31            
## [28] jquerylib_0.1.4         cellranger_1.1.0        numDeriv_2016.8-1.1    
## [31] estimability_1.5.1      Rcpp_1.1.0              knitr_1.50             
## [34] zoo_1.8-14              parameters_0.28.3       splines_4.5.1          
## [37] timechange_0.3.0        tidyselect_1.2.1        rstudioapi_0.17.1      
## [40] abind_1.4-8             yaml_2.3.10             codetools_0.2-20       
## [43] lattice_0.22-7          withr_3.0.2             bayestestR_0.17.0      
## [46] S7_0.2.0                askpass_1.2.1           coda_0.19-4.1          
## [49] evaluate_1.0.5          survival_3.8-3          zip_2.3.3              
## [52] xml2_1.4.0              pillar_1.11.0           reformulas_0.4.1       
## [55] insight_1.4.4           generics_0.1.4          rprojroot_2.1.1        
## [58] hms_1.1.3               minqa_1.2.8             glue_1.8.0             
## [61] gdtools_0.4.4           tools_4.5.1             data.table_1.17.8      
## [64] fs_1.6.6                mvtnorm_1.3-3           grid_4.5.1             
## [67] rbibutils_2.3           datawizard_1.3.0        nlme_3.1-168           
## [70] Formula_1.2-5           cli_3.6.5               textshaping_1.0.3      
## [73] fontBitstreamVera_0.1.1 gtable_0.3.6            sass_0.4.10            
## [76] digest_0.6.37           fontquiver_0.2.1        pbkrtest_0.5.5         
## [79] TH.data_1.1-4           farver_2.1.2            htmltools_0.5.8.1      
## [82] lifecycle_1.0.4         fontLiberation_0.1.0    openssl_2.3.4          
## [85] MASS_7.3-65

10 References

Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). Lawrence Erlbaum Associates.

Hoenig, J. M., & Heisey, D. M. (2001). The abuse of power: The pervasive fallacy of power calculations for data analysis. The American Statistician, 55(1), 19-24.

---

End of Report