Purpose
Purpose
1. Setup
2. Load and prepare data
3. Overall learning outcomes
- 3.1 Primary paired tests
- 3.2 Figure 1: Student gains
4. Mastery of laboratory skills
5. Affective factors associated with student success
6. Preparedness gaps
7. Gender differences in change in grit
8. Export analysis tables
9. Session information

Purpose

This file is the streamlined, reproducible analysis for the manuscript. It retains only:

Overall learning outcomes.
Mastery progression and attainment.
Perseverance and first-pass mastery.
Growth in consistency of interests and experimental-design competency.
Relative importance of baseline competency, consistency growth, and high school GPA.
Preparedness-gap analyses.

Where a figure was present in POBRebootAnalysis_edited.Rmd, its plotting code has been copied verbatim so that its appearance is unchanged.

Important: The edited source still does not contain the code used to create the poster figure titled Students Maintained Confidence. A clearly marked placeholder is included below rather than silently reconstructing a different figure.

Purpose

This streamlined analysis file contains only the data preparation, statistical tests, and figures supporting the current paper:

Mastery of laboratory skills.
Affective factors associated with mastery and experimental-design competency.
Reduction of preparedness gaps.

The plotting code for the retained figures is preserved from the original analysis so that their appearance does not change.

1. Setup

library(dplyr)
library(tidyr)
library(readr)
library(stringr)
library(forcats)
library(ggplot2)
library(ggalluvial)
library(patchwork)
library(scales)
library(broom)
library(relaimpo)

theme_paper  <- function(base_size = 16){

  theme_classic(base_size = base_size) +

    theme(

      panel.background = element_rect(
        fill = "transparent",
        color = NA
      ),

      plot.background = element_rect(
        fill = "transparent",
        color = NA
      ),

      legend.background = element_rect(
        fill = "transparent",
        color = NA
      ),

      legend.key = element_rect(
        fill = "transparent",
        color = NA
      ),

      strip.background = element_rect(
        fill = "white",
        color = "black"
      ),

      strip.text = element_text(
        face = "bold",
        color = "#353C47"
      ),

      axis.text = element_text(
        color = "#353C47",
        face = "bold"
      ),

      axis.title = element_text(
        color = "#353C47",
        face = "bold"
      ),

      plot.title = element_text(
        color = "#353C47",
        face = "bold",
        hjust = .5
      ),

      plot.subtitle = element_text(
        color = "#353C47"
      ),

      plot.caption = element_text(
        color = "gray40",
        face = "italic"
      ),

      legend.title = element_text(
        color = "#353C47",
        face = "bold"
      ),

      legend.text = element_text(
        color = "#353C47"
      ),

      axis.line = element_line(
        color = "#353C47"
      ),

      axis.ticks = element_line(
        color = "#353C47"
      )
    )
}

first_nonmissing <- function(x) {
  x <- x[!is.na(x)]
  if (length(x) == 0) NA else x[1]
}

attempt_cols <- c(
  "Pipette.Pass",
  "Excel.Pass",
  "Microscope.Pass",
  "PCR.Pass",
  "Gel.Pass",
  "Sequencing.Pass"
)

2. Load and prepare data

This file expects the following source files in the working directory:

POBRebootCombined.csv
POB1_StudentIDs.csv
StudentDemo.csv
SAlong.csv

2.1 Import and merge

dat <- read.csv("POBRebootCombined.csv")
names_file <- read.csv("POB1_StudentIDs.csv")
demo <- read.csv("StudentDemo.csv")
skills <- read.csv("SAlong.csv")

dat_full <- merge(dat, names_file, by = "StudentID", all.x = TRUE)

demo <- demo %>%
  rename(StudentID = Student.ID)

dat_full <- merge(dat_full, demo, by = "StudentID", all.x = TRUE)

dat_full <- dat_full %>%
  rename(Cohort.TR = Cohort) %>%
  mutate(
    Cohort = substr(Cohort.TR, 1, 4),
    Transfer = ifelse(grepl("TR", Cohort.TR), 1, 0),
    Gender = case_when(
      Gender == "Female" ~ "Female",
      Gender == "Male" ~ "Male",
      TRUE ~ "Other"
    ),
    Race = case_when(
      Race == "White" ~ "White",
      TRUE ~ "PEER"
    )
  )

2.2 Academic preparation variables

dat_full <- dat_full %>%
  mutate(
    across(
      c(HS.GPA_IR, SATmath, SATeng, ACT),
      ~ na_if(.x, "TR student")
    )
  )

convert_sat_to_act <- function(sat_total) {
  case_when(
    sat_total >= 1570 ~ 36,
    sat_total >= 1530 ~ 35,
    sat_total >= 1490 ~ 34,
    sat_total >= 1450 ~ 33,
    sat_total >= 1420 ~ 32,
    sat_total >= 1390 ~ 31,
    sat_total >= 1360 ~ 30,
    sat_total >= 1330 ~ 29,
    sat_total >= 1290 ~ 28,
    sat_total >= 1250 ~ 27,
    sat_total >= 1210 ~ 26,
    sat_total >= 1170 ~ 25,
    sat_total >= 1130 ~ 24,
    sat_total >= 1090 ~ 23,
    sat_total >= 1050 ~ 22,
    sat_total >= 1010 ~ 21,
    sat_total >= 970  ~ 20,
    sat_total >= 930  ~ 19,
    sat_total >= 890  ~ 18,
    sat_total >= 850  ~ 17,
    sat_total >= 810  ~ 16,
    sat_total >= 770  ~ 15,
    sat_total >= 730  ~ 14,
    sat_total >= 690  ~ 13,
    sat_total >= 650  ~ 12,
    sat_total >= 620  ~ 11,
    sat_total >= 590  ~ 10,
    sat_total >= 560  ~ 9,
    TRUE ~ NA_real_
  )
}

dat_full <- dat_full %>%
  mutate(
    SATmath = as.numeric(SATmath),
    SATeng = as.numeric(SATeng),
    SAT_total = SATmath + SATeng,
    ACT_equiv = convert_sat_to_act(SAT_total),
    ACT_numeric = suppressWarnings(as.numeric(ACT)),
    ACT.complete = coalesce(ACT_equiv, ACT_numeric),
    HS.GPA_IR = suppressWarnings(as.numeric(HS.GPA_IR))
  ) %>%
  dplyr::select(-ACT_numeric)

## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `SATmath = as.numeric(SATmath)`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.

2.3 Score grit

grit_scale_map <- c(
  "Not at all like me" = 1,
  "Not much like me" = 2,
  "Somewhat like me" = 3,
  "Mostly like me" = 4,
  "Very much like me" = 5
)

grit_items <- paste0("Grit", 1:12)

dat_full[grit_items] <- lapply(
  dat_full[grit_items],
  function(col) as.numeric(grit_scale_map[as.character(col)])
)

reverse_items <- c("Grit2", "Grit3", "Grit5", "Grit7", "Grit8", "Grit11")
dat_full[reverse_items] <- lapply(
  dat_full[reverse_items],
  function(col) 6 - col
)

perseverance_items <- c("Grit1", "Grit4", "Grit6", "Grit9", "Grit10", "Grit12")
consistency_items <- c("Grit2", "Grit3", "Grit5", "Grit7", "Grit8", "Grit11")

dat_full <- dat_full %>%
  mutate(
    Grit_Total = rowMeans(across(all_of(grit_items)), na.rm = TRUE),
    Grit_Perseverance = rowMeans(across(all_of(perseverance_items)), na.rm = TRUE),
    Grit_Consistency = rowMeans(across(all_of(consistency_items)), na.rm = TRUE)
  )

2.4 Score STEP-U, RSQ, and EDCI

step_map <- c(
  "Not Important" = 1,
  "Slightly important" = 2,
  "Fairly important" = 3,
  "Important" = 4,
  "Very Important" = 5
)

step_items <- paste0("STEP", 1:14)

dat_full[step_items] <- lapply(
  dat_full[step_items],
  function(col) as.numeric(step_map[as.character(col)])
)

rsq_map <- c(
  "Somewhat Confident" = 1,
  "Moderately Confident" = 2,
  "Confident" = 3,
  "Very Confident" = 4
)

rsq_items <- paste0("RSQ", 1:12)

dat_full[rsq_items] <- lapply(
  dat_full[rsq_items],
  function(col) as.numeric(rsq_map[as.character(col)])
)

correct_answers <- c(
  "No - mice should face a random direction and the enclosure should rotate randomly",
  "It should lead to mice with variable characteristics being distributed fairly evenly",
  "The other two methods do not directly assess weight changes in mice",
  "It provides more chances to compare activity rate across different temperatures",
  "Other variables as well as temperature will differ between treatment groups",
  "Mice in each treatment group should not vary in size, age or sex (small, young, females only)",
  "The 14C water temperature treatment group",
  "Take measurements from 14 fish in one temperature group but take measurements from all 15 fish in the other temperature groups (total n=44)",
  "The weights of fish in each treatment group were too variable",
  "Sex, age, health and weight at start",
  "Use a set of weighing scales that is able to measure smaller weights",
  "You can test all four of these hypotheses",
  "Spray plain water on another group of 150 tomato plants to compare results with P1, P2 and P3",
  "Potting soil type and temperature",
  "If you are able to support a hypothesis (either the alternate or null hypothesis) in all three experiments",
  "Any of the above sources could affect data more than the others in any given experiment",
  "Splitting the treatment groups into thirds is unfair as field conditions might vary",
  "Count aphids on 1 randomly selected leaf from 100 plants"
)

question_cols <- paste0("ED", 1:18)

for (i in seq_along(question_cols)) {
  col <- question_cols[i]
  student_answer <- trimws(iconv(dat_full[[col]], from = "", to = "UTF-8"))
  correct_answer <- trimws(correct_answers[i])

  dat_full[[paste0(col, "_correct")]] <- case_when(
    is.na(student_answer) | student_answer == "" ~ NA_integer_,
    student_answer == correct_answer ~ 1L,
    TRUE ~ 0L
  )
}

ed_correct_cols <- paste0(question_cols, "_correct")

dat_full <- dat_full %>%
  mutate(
    STEP_Total = rowMeans(across(all_of(step_items)), na.rm = TRUE),
    RSQ_Total = rowMeans(across(all_of(rsq_items)), na.rm = TRUE),
    EDCI_Total_Score = rowSums(across(all_of(ed_correct_cols)), na.rm = TRUE),
    EDCI_Proportion_Correct = rowMeans(across(all_of(ed_correct_cols)), na.rm = TRUE)
  )

2.5 Merge mastery-assessment outcomes

dat_full <- dat_full %>%
  left_join(skills, by = c("StudentID" = "Student_ID"))

## Warning in left_join(., skills, by = c(StudentID = "Student_ID")): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 49 of `x` matches multiple rows in `y`.
## ℹ Row 108 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.

mastery_levels <- c("First", "Second", "Third", "No Mastery")

dat_full <- dat_full %>%
  mutate(
    across(
      all_of(attempt_cols),
      ~ factor(.x, levels = mastery_levels, ordered = TRUE)
    )
  )

# Create totals only if they are not already supplied in SAlong.csv.
if (!"Total.Pass" %in% names(dat_full)) {
  dat_full <- dat_full %>%
    mutate(
      Total.Pass = rowSums(
        across(
          all_of(attempt_cols),
          ~ !is.na(.x) & as.character(.x) != "No Mastery"
        )
      )
    )
}

if (!"Total.NoMastery" %in% names(dat_full)) {
  dat_full <- dat_full %>%
    mutate(
      Total.NoMastery = rowSums(
        across(
          all_of(attempt_cols),
          ~ as.character(.x) == "No Mastery"
        ),
        na.rm = TRUE
      )
    )
}

2.6 Create paired student-level dataset

dat_paired <- dat_full %>%
  group_by(StudentID) %>%
  filter(all(c("Pre", "Post") %in% Timepoint)) %>%
  group_by(StudentID, Timepoint) %>%
  slice_max(Progress, n = 1, with_ties = FALSE) %>%
  ungroup()

scores_wide <- dat_paired %>%
  dplyr::select(
    StudentID,
    Timepoint,
    EDCI_Total_Score,
    STEP_Total,
    RSQ_Total,
    Grit_Total,
    Grit_Perseverance,
    Grit_Consistency
  ) %>%
  pivot_wider(
    id_cols = StudentID,
    names_from = Timepoint,
    values_from = c(
      EDCI_Total_Score,
      STEP_Total,
      RSQ_Total,
      Grit_Total,
      Grit_Perseverance,
      Grit_Consistency
    )
  )

predictors <- dat_paired %>%
  group_by(StudentID) %>%
  summarise(
    Gender = first_nonmissing(Gender),
    Race = first_nonmissing(Race),
    Transfer = first_nonmissing(Transfer),
    HS.GPA_IR = first_nonmissing(HS.GPA_IR),
    ACT.complete = first_nonmissing(ACT.complete),
    Total.Pass = first_nonmissing(Total.Pass),
    Total.NoMastery = first_nonmissing(Total.NoMastery),
    .groups = "drop"
  )

paired_wide <- scores_wide %>%
  left_join(predictors, by = "StudentID") %>%
  mutate(
    HS.GPA_IR = as.numeric(HS.GPA_IR),

    EDCI_Gain = EDCI_Total_Score_Post - EDCI_Total_Score_Pre,
    STEP_Gain = STEP_Total_Post - STEP_Total_Pre,
    RSQ_Gain = RSQ_Total_Post - RSQ_Total_Pre,

    Grit_Total_Gain = Grit_Total_Post - Grit_Total_Pre,
    Perseverance_Gain = Grit_Perseverance_Post - Grit_Perseverance_Pre,
    Consistency_Gain = Grit_Consistency_Post - Grit_Consistency_Pre,

    # Preserve names used by the original figure code.
    EDCI_Gain.x = EDCI_Gain,
    STEP_Gain.x = STEP_Gain,
    RSQ_Gain.x = RSQ_Gain
  ) %>%
  filter(
  !(STEP_Total_Pre == 1 | STEP_Total_Post == 1),
  !(EDCI_Total_Score_Post == 0)
)

3. Overall learning outcomes

3.1 Primary paired tests

paired_tests <- list(
  EDCI = t.test(
    paired_wide$EDCI_Total_Score_Post,
    paired_wide$EDCI_Total_Score_Pre,
    paired = TRUE
  ),
  RSQ = t.test(
    paired_wide$RSQ_Total_Post,
    paired_wide$RSQ_Total_Pre,
    paired = TRUE
  ),
  STEP = t.test(
    paired_wide$STEP_Total_Post,
    paired_wide$STEP_Total_Pre,
    paired = TRUE
  )
)

paired_tests

## $EDCI
## 
##  Paired t-test
## 
## data:  paired_wide$EDCI_Total_Score_Post and paired_wide$EDCI_Total_Score_Pre
## t = 2.0002, df = 55, p-value = 0.05042
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.002036708  2.144893850
## sample estimates:
## mean difference 
##        1.071429 
## 
## 
## $RSQ
## 
##  Paired t-test
## 
## data:  paired_wide$RSQ_Total_Post and paired_wide$RSQ_Total_Pre
## t = 4.0965, df = 53, p-value = 0.0001443
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.1596604 0.4659942
## sample estimates:
## mean difference 
##       0.3128273 
## 
## 
## $STEP
## 
##  Paired t-test
## 
## data:  paired_wide$STEP_Total_Post and paired_wide$STEP_Total_Pre
## t = 0.95252, df = 55, p-value = 0.345
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.06984441  0.19638156
## sample estimates:
## mean difference 
##      0.06326858

3.2 Figure 1: Student gains

The following block is copied directly from the edited source file.

Figure 1: (A) EDCI pre/post, (B) RSQ pre/post, (C) STEP pre/post

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ lubridate 1.9.5     ✔ tibble    3.3.1
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ scales::col_factor() masks readr::col_factor()
## ✖ purrr::discard()     masks scales::discard()
## ✖ Matrix::expand()     masks tidyr::expand()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ Matrix::pack()       masks tidyr::pack()
## ✖ MASS::select()       masks dplyr::select()
## ✖ Matrix::unpack()     masks tidyr::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggplot2)
library(patchwork)

gain_long <- paired_wide %>%
  dplyr::select(
    StudentID,
    EDCI_Total_Score_Pre, EDCI_Total_Score_Post,
    RSQ_Total_Pre, RSQ_Total_Post,
    STEP_Total_Pre, STEP_Total_Post
  ) %>%
  pivot_longer(
    -StudentID,
    names_to = c("Measure", "Time"),
    names_pattern = "(.*)_(Pre|Post)",
    values_to = "Score"
  ) %>%
  mutate(
    Measure = recode(
      Measure,
      "EDCI_Total_Score" = "Experimental Design\nCompetency",
      "RSQ_Total" = "Research\nSelf-Efficacy",
      "STEP_Total" = "Science\nAttitudes"
    ),
    Time = factor(Time, levels = c("Pre", "Post"))
  )

fig1 <- ggplot(gain_long, aes(Time, Score, group = StudentID)) +
  geom_line(alpha = 0.25) +
  geom_point(alpha = 0.45, size = 1.8) +
  stat_summary(aes(group = 1), fun = mean, geom = "line", linewidth = 1.2) +
  stat_summary(aes(group = 1), fun = mean, geom = "point", size = 3) +
  facet_wrap(~Measure, scales = "free_y") +
  theme_classic(base_size = 13) +
  labs(x = NULL, y = "Score")

fig1

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 2 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("Figure1_Student_Gains.png", fig1, width = 9, height = 4.5, dpi = 600)

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

Figure 2: (A) EDCI gain by baseline EDCI quartile, (B) RSQ gain by baseline RSQ quartile

4. Mastery of laboratory skills

4.1 Descriptive mastery statistics

mastery_student <- dat_full %>%
  dplyr::select(StudentID, Total.Pass, Total.NoMastery, all_of(attempt_cols)) %>%
  distinct()

mastery_student %>%
  summarise(
    N = sum(!is.na(Total.Pass)),
    Mean_Skills_Mastered = mean(Total.Pass, na.rm = TRUE),
    SD_Skills_Mastered = sd(Total.Pass, na.rm = TRUE),
    Percent_Mastered_5_or_6 = mean(Total.Pass >= 5, na.rm = TRUE) * 100,
    Mean_No_Mastery = mean(Total.NoMastery, na.rm = TRUE)
  )

##     N Mean_Skills_Mastered SD_Skills_Mastered Percent_Mastered_5_or_6
## 1 148             4.777027            1.33412                64.18919
##   Mean_No_Mastery
## 1        1.222973

4.2 Figure: Distribution of mastery and mean attempts

The following blocks are copied directly from the edited source file.

attempt_cols <- c(
  "Pipette.Pass",
  "Excel.Pass",
  "Microscope.Pass",
  "PCR.Pass",
  "Gel.Pass",
  "Sequencing.Pass"
)

attempt_long <- dat_full %>%
  dplyr::select(all_of(attempt_cols)) %>%
  pivot_longer(
    everything(),
    names_to = "Skill",
    values_to = "Attempt"
  ) %>%
  filter(!is.na(Attempt)) %>%
  mutate(
    Skill = str_remove(Skill, "\\.Pass"),
    Skill = factor(
      Skill,
      levels = c("Pipette", "Excel", "Microscope", "PCR", "Gel", "Sequencing")
    ),
    Attempt = factor(
      Attempt,
      levels = c("First", "Second", "Third", "No Mastery")
    )
  )

fig3 <- ggplot(
  attempt_long,
  aes(x = Skill, fill = Attempt)
) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  theme_classic(base_size = 13) +
  labs(
    x = "Skills Assessment",
    y = "Percent of Students",
    fill = "Mastery Attempt",
    title = "Most students achieved mastery, but differed in number of attempts"
  ) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

fig3

ggsave("Figure3_Mastery_Attempts.png", fig3, width = 7, height = 4.5, dpi = 600)

Options presented here

# ============================================================
# FIGURE 3B. Mean attempts required versus perseverance
# ============================================================
attempt_map <- c(
  "First" = 1,
  "Second" = 2,
  "Third" = 3,
  "No Mastery" = 4
)

attempt_num <- dat_full %>%
  dplyr::select(
    StudentID,
    Grit_Perseverance,
    all_of(attempt_cols)
  ) %>%
  distinct()

attempt_num[attempt_cols] <- lapply(
  attempt_num[attempt_cols],
  function(x) as.numeric(attempt_map[as.character(x)])
)

attempt_num <- attempt_num %>%
  mutate(
    Number_First_Pass = rowSums(
      across(all_of(attempt_cols), ~ .x == 1),
      na.rm = TRUE
    ),
    Mean_Attempt = rowMeans(
      across(all_of(attempt_cols)),
      na.rm = TRUE
    ),
    Perseverance_Quartile = ntile(Grit_Perseverance, 4),
    Perseverance_Quartile_Label = factor(
      Perseverance_Quartile,
      levels = 1:4,
      labels = c("Lowest", "Q2", "Q3", "Highest")
    )
  ) %>%
  filter(
    !is.na(Grit_Perseverance),
    !is.na(Number_First_Pass)
  )
attempt_num <- attempt_num %>%
  mutate(
    Perseverance_Quartile = ntile(Grit_Perseverance, 4)
  )

ggplot(
  attempt_num,
  aes(
    factor(Perseverance_Quartile),
    Mean_Attempt
  )
) +
  geom_boxplot() +
  geom_jitter(
    width = .15,
    alpha = .4
  ) +
  theme_classic() +
  labs(
    x = "Perseverance Quartile",
    y = "Mean Attempts Required"
  )

## Warning: Removed 7 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).

attempt_cols <- c(
  "Pipette.Pass",
  "Excel.Pass",
  "Microscope.Pass",
  "PCR.Pass",
  "Gel.Pass",
  "Sequencing.Pass"
)

attempt_map <- c(
  "First" = 1,
  "Second" = 2,
  "Third" = 3,
  "No Mastery" = 4
)

attempt_num <- dat_full %>%
  dplyr::select(
    StudentID,
    Grit_Perseverance,
    all_of(attempt_cols)
  ) %>%
  distinct()

attempt_num[attempt_cols] <- lapply(
  attempt_num[attempt_cols],
  function(x) {
    as.numeric(attempt_map[as.character(x)])
  }
)

attempt_num <- attempt_num %>%
  mutate(
    Mean_Attempt = rowMeans(
      across(all_of(attempt_cols)),
      na.rm = TRUE
    ),

    Number_First_Pass = rowSums(
      across(all_of(attempt_cols), ~ .x == 1),
      na.rm = TRUE
    ),

    Perseverance_Quartile = ntile(
      Grit_Perseverance,
      4
    )
  ) %>%
  filter(
    !is.na(Grit_Perseverance),
    !is.na(Number_First_Pass)
  )

ggplot(
  attempt_num,
  aes(
    Grit_Perseverance,
    Number_First_Pass
  )
) +
  geom_jitter(
    width = .1,
    height = .1,
    alpha = .5
  ) +
  geom_smooth(
    method = "lm"
  ) +
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

ggplot(
  attempt_num,
  aes(
    factor(ntile(Grit_Perseverance,4)),
    Number_First_Pass
  )
) +
  geom_boxplot()

ggplot(
  attempt_num,
  aes(
    factor(ntile(Grit_Perseverance,4)),
    Number_First_Pass
  )
) +
  geom_boxplot() +
  geom_jitter(width=.15)

4.3 Figure: Mastery progression

The following block is copied directly from the edited source file.

# ============================================================
# Sankey / alluvial diagram: mastery flow stops after passing
# ============================================================

library(dplyr)
library(tidyr)
library(ggplot2)
library(ggalluvial)
library(stringr)
library(scales)
library(ggalluvial)

attempt_cols <- c(
  "Pipette.Pass",
  "Excel.Pass",
  "Microscope.Pass",
  "PCR.Pass",
  "Gel.Pass",
  "Sequencing.Pass"
)

# Use one row per student to avoid duplicated pre/post rows
mastery_wide <- dat_full %>%
  dplyr::select(StudentID, all_of(attempt_cols)) %>%
  distinct()

# One row per student-skill assessment
mastery_long <- mastery_wide %>%
  pivot_longer(
    cols = all_of(attempt_cols),
    names_to = "Skill",
    values_to = "Final_Attempt"
  ) %>%
  filter(!is.na(Final_Attempt)) %>%
  mutate(
    Skill = str_remove(Skill, "\\.Pass"),
    Final_Attempt = as.character(Final_Attempt)
  )

# Create pathways where flow stops once students pass
mastery_flow <- mastery_long %>%
  mutate(
    Attempt_1 = case_when(
      Final_Attempt == "First" ~ "Pass",
      Final_Attempt %in% c("Second", "Third", "No Mastery") ~ "Fail",
      TRUE ~ NA_character_
    ),

    Attempt_2 = case_when(
      Final_Attempt == "First" ~ NA_character_,
      Final_Attempt == "Second" ~ "Pass",
      Final_Attempt %in% c("Third", "No Mastery") ~ "Fail",
      TRUE ~ NA_character_
    ),

    Attempt_3 = case_when(
      Final_Attempt %in% c("First", "Second") ~ NA_character_,
      Final_Attempt == "Third" ~ "Pass",
      Final_Attempt == "No Mastery" ~ "No Mastery",
      TRUE ~ NA_character_
    )
  )

# Convert to alluvial-ready long format
flow_long <- mastery_flow %>%
  mutate(
    Pathway = row_number()
  ) %>%
  dplyr::select(Pathway, Attempt_1, Attempt_2, Attempt_3) %>%
  pivot_longer(
    cols = starts_with("Attempt"),
    names_to = "Attempt",
    values_to = "Outcome"
  ) %>%
  filter(!is.na(Outcome)) %>%
  mutate(
    Attempt = recode(
      Attempt,
      "Attempt_1" = "Attempt 1",
      "Attempt_2" = "Attempt 2",
      "Attempt_3" = "Attempt 3"
    ),
    Attempt = factor(
      Attempt,
      levels = c("Attempt 1", "Attempt 2", "Attempt 3")
    ),
    Outcome = factor(
      Outcome,
      levels = c("Pass", "Fail", "No Mastery")
    )
  )

# Summarize counts for checking
flow_counts <- flow_long %>%
  count(Attempt, Outcome) %>%
  group_by(Attempt) %>%
  mutate(
    Percent = n / sum(n)
  )

flow_counts

## # A tibble: 6 × 4
## # Groups:   Attempt [3]
##   Attempt   Outcome        n Percent
##   <fct>     <fct>      <int>   <dbl>
## 1 Attempt 1 Pass         408   0.570
## 2 Attempt 1 Fail         308   0.430
## 3 Attempt 2 Pass         236   0.766
## 4 Attempt 2 Fail          72   0.234
## 5 Attempt 3 Pass          63   0.875
## 6 Attempt 3 No Mastery     9   0.125

# Plot
fig_sankey <- ggplot(
  flow_long,
  aes(
    x = Attempt,
    stratum = Outcome,
    alluvium = Pathway,
    y = 1,
    fill = Outcome,
    label = Outcome
  )
) +
 geom_alluvium(
  alpha = 0.75
) +
  geom_stratum(
    width = 0.22,
    color = "black"
  ) +
  geom_text(
    stat = "stratum",
    aes(
      label = after_stat(
        paste0(
          stratum,
          "\n",
          n,
          " (",
          percent(n / tapply(n, x, sum)[as.character(x)], accuracy = 1),
          ")"
        )
      )
    ),
    size = 3.3
  ) +
  scale_fill_manual(
    values = c(
      "Pass" = "#33A02C",
      "Fail" = "#E31A1C",
      "No Mastery" = "#6A3D9A"
    )
  ) +
  theme_classic(base_size = 13) +
  labs(
    x = NULL,
    y = "Number of student-skill attempts",
    fill = "Outcome",
    title = "Mastery progression across skills assessments"
  ) +
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(face = "bold", size = 12),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

fig_sankey

## Warning in setup_data(...): Some differentiation aesthetics vary within alluvia, and will be diffused by their first value.
## Consider using `geom_flow()` instead.

ggsave(
  "Figure_Mastery_Sankey_StopAfterPass.png",
  fig_sankey,
  width = 8.5,
  height = 5,
  dpi = 600
)

## Warning in setup_data(...): Some differentiation aesthetics vary within alluvia, and will be diffused by their first value.
## Consider using `geom_flow()` instead.

Figure 4: Forest Plot. EDCI, RSQ, STEP (Predictors Baseline score, Grit, GPA, Mastery, and Collaboration)

4.4 Confidence across mastery attempts

The code for the poster panel Students Maintained Confidence was not present in either uploaded R Markdown file. Insert the original code here once located.

# ------------------------------------------------------------
# Confidence across final mastery-attempt categories
# One observation = one student × one laboratory skill
# ------------------------------------------------------------

# One row per student with post-course confidence
post_confidence <- paired_wide %>%
  dplyr::select(
    StudentID,
    RSQ_Total_Post
  ) %>%
  filter(!is.na(RSQ_Total_Post))

# One row per student × laboratory skill
confidence_attempt_data <- dat_full %>%
  dplyr::select(
    StudentID,
    all_of(attempt_cols)
  ) %>%
  distinct() %>%
  pivot_longer(
    cols = all_of(attempt_cols),
    names_to = "Skill",
    values_to = "Final_Mastery_Attempt"
  ) %>%
  left_join(
    post_confidence,
    by = "StudentID"
  ) %>%
  filter(
    !is.na(RSQ_Total_Post),
    !is.na(Final_Mastery_Attempt),
    Final_Mastery_Attempt != "No Mastery"   # <-- remove these observations
  ) %>%
  mutate(
    Final_Mastery_Attempt = factor(
      Final_Mastery_Attempt,
      levels = c("First", "Second", "Third")
    )
  )


confidence_attempt_summary <- confidence_attempt_data %>%
  group_by(Final_Mastery_Attempt) %>%
  summarise(
    n = n(),
    Mean = mean(RSQ_Total_Post, na.rm = TRUE),
    SD = sd(RSQ_Total_Post, na.rm = TRUE),
    SE = SD / sqrt(n),
    .groups = "drop"
  )

confidence_attempt_summary

confidence_attempt_plot <- ggplot(
  confidence_attempt_data,
  aes(
    x = Final_Mastery_Attempt,
    y = RSQ_Total_Post
  )
) +

  # Distribution of confidence scores
  geom_boxplot(
    outlier.shape = NA,
    width = 0.62,
    fill = "#CFF4EF",
    color = "#353C47",
    linewidth = 0.9
  ) +

  # Individual student-skill observations
  geom_jitter(
    width = 0.10,
    height = 0.025,
    size = 2.5,
    alpha = 0.80,
    color = "#F4A261"
  ) +

  # Mean confidence
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 21,
    size = 4.5,
    stroke = 0.8,
    fill = "#29B897",
    color = "#353C47"
  ) +

  # Mean ± standard error
  stat_summary(
    fun.data = mean_se,
    geom = "errorbar",
    width = 0.16,
    linewidth = 0.9,
    color = "#353C47"
  ) +

  scale_y_continuous(
    breaks = 1:4,
    limits = c(0.8, 4.25),
    expand = expansion(mult = c(0, 0.02))
  ) +

  labs(
    title = "Students Maintained Confidence",
    x = "Final Mastery Attempt",
    y = "Post-Course Confidence"
  ) +

  theme_paper(base_size = 16) +

  theme(
    plot.title = element_text(
      face = "bold",
      size = 18,
      hjust = 0.5,
      color = "#353C47"
    ),

    axis.title = element_text(
      face = "bold",
      color = "#353C47"
    ),

    axis.text = element_text(
      face = "bold",
      color = "#353C47"
    ),

    axis.text.x = element_text(
      size = 11
    ),

    panel.grid = element_blank(),

    plot.background = element_rect(
      fill = "transparent",
      color = NA
    ),

    panel.background = element_rect(
      fill = "transparent",
      color = NA
    )
  )

confidence_attempt_plot

ggsave(
  "Confidence_by_Final_Mastery_Attempt.png",
  plot = confidence_attempt_plot,
  width = 6,
  height = 4.5,
  dpi = 600,
  bg = "transparent"
)

library(ordinal)

confidence_attempt_data <- confidence_attempt_data %>%
  mutate(
    Attempt_Ordinal = ordered(
      Final_Mastery_Attempt,
      levels = c("First", "Second", "Third")
    )
  )

confidence_attempt_model <- clmm(
  Attempt_Ordinal ~
    RSQ_Total_Post +
    (1 | StudentID) +
    (1 | Skill),
  data = confidence_attempt_data,
  link = "logit",
  Hess = TRUE
)

summary(confidence_attempt_model)

confidence_beta <- coef(summary(confidence_attempt_model))[
  "RSQ_Total_Post",
  "Estimate"
]

confidence_or <- exp(confidence_beta)

confidence_or

confidence_ci <- confint(
  confidence_attempt_model,
  parm = "RSQ_Total_Post"
)

exp(confidence_ci)

confidence_p <- coef(summary(confidence_attempt_model))[
  "RSQ_Total_Post",
  "Pr(>|z|)"
]

confidence_p_label <- paste0(
  "p = ",
  format.pval(confidence_p, digits = 2, eps = 0.001)
)


confidence_attempt_plot <- confidence_attempt_plot +
  annotate(
    "text",
    x = Inf,
    y = Inf,
    label = confidence_p_label,
    hjust = 1.15,
    vjust = 1.5,
    fontface = "bold",
    size = 5,
    color = "#29B897"
  )

confidence_attempt_plot

5. Affective factors associated with student success

5.1 Perseverance and first-pass mastery

attempt_map <- c(
  "First" = 1,
  "Second" = 2,
  "Third" = 3,
  "No Mastery" = 4
)

attempt_num <- dat_full %>%
  dplyr::select(
    StudentID,
    Grit_Perseverance,
    all_of(attempt_cols)
  ) %>%
  distinct()

attempt_num[attempt_cols] <- lapply(
  attempt_num[attempt_cols],
  function(x) as.numeric(attempt_map[as.character(x)])
)

attempt_num <- attempt_num %>%
  mutate(
    Number_First_Pass = rowSums(
      across(all_of(attempt_cols), ~ .x == 1),
      na.rm = TRUE
    ),
    Mean_Attempt = rowMeans(
      across(all_of(attempt_cols)),
      na.rm = TRUE
    ),
    Perseverance_Quartile = ntile(Grit_Perseverance, 4),
    Perseverance_Quartile_Label = factor(
      Perseverance_Quartile,
      levels = 1:4,
      labels = c("Lowest", "Q2", "Q3", "Highest")
    )
  ) %>%
  filter(
    !is.na(Grit_Perseverance),
    !is.na(Number_First_Pass)
  )

first_pass_model <- lm(
  Number_First_Pass ~ Grit_Perseverance,
  data = attempt_num
)

summary(first_pass_model)

## 
## Call:
## lm(formula = Number_First_Pass ~ Grit_Perseverance, data = attempt_num)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7200 -1.2721 -0.3588  1.2077  3.7712 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         1.7086     0.4389   3.893 0.000121 ***
## Grit_Perseverance   0.2167     0.1142   1.897 0.058700 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.453 on 315 degrees of freedom
## Multiple R-squared:  0.0113, Adjusted R-squared:  0.00816 
## F-statistic:   3.6 on 1 and 315 DF,  p-value: 0.0587

5.2 Figure: First-pass mastery by perseverance quartile

The following block is copied directly from the edited source file.

# ============================================================
# FIGURE 2. First-pass mastery by perseverance quartile
# ============================================================

library(dplyr)
library(ggplot2)

attempt_cols <- c(
  "Pipette.Pass",
  "Excel.Pass",
  "Microscope.Pass",
  "PCR.Pass",
  "Gel.Pass",
  "Sequencing.Pass"
)

attempt_map <- c(
  "First" = 1,
  "Second" = 2,
  "Third" = 3,
  "No Mastery" = 4
)

attempt_num <- dat_full %>%
  dplyr::select(
    StudentID,
    Grit_Perseverance,
    all_of(attempt_cols)
  ) %>%
  distinct()

attempt_num[attempt_cols] <- lapply(
  attempt_num[attempt_cols],
  function(x) {
    as.numeric(attempt_map[as.character(x)])
  }
)

attempt_num <- attempt_num %>%
  mutate(
    Number_First_Pass = rowSums(
      across(all_of(attempt_cols), ~ .x == 1),
      na.rm = TRUE
    ),
    Mean_Attempt = rowMeans(
      across(all_of(attempt_cols)),
      na.rm = TRUE
    ),
    Perseverance_Quartile = ntile(Grit_Perseverance, 4)
  ) %>%
  filter(
    !is.na(Grit_Perseverance),
    !is.na(Number_First_Pass)
  )

fig2_explor <- ggplot(
  attempt_num,
  aes(
    x = factor(Perseverance_Quartile),
    y = Number_First_Pass
  )
) +
  geom_boxplot(
    outlier.shape = NA,
    width = 0.6
  ) +
  geom_jitter(
    width = 0.12,
    alpha = 0.55,
    size = 2
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3
  ) +
  theme_classic(base_size = 13) +
  labs(
    x = "Perseverance of Effort Quartile",
    y = "Skills Mastered on First Attempt"
  )

fig2_explor

ggsave(
  "BriefReport_Figure2_FirstPass_by_Perseverance.png",
  fig2_explor,
  width = 6,
  height = 4.5,
  dpi = 600
)

Figure 3:

5.3 Consistency of interests and EDCI outcome

consistency_model_data <- paired_wide %>%
  dplyr::select(
    EDCI_Total_Score_Post,
    EDCI_Total_Score_Pre,
    HS.GPA_IR,
    Consistency_Gain
  ) %>%
  filter(complete.cases(.))

consistency_model <- lm(
  EDCI_Total_Score_Post ~
    EDCI_Total_Score_Pre +
    HS.GPA_IR +
    Consistency_Gain,
  data = consistency_model_data
)

summary(consistency_model)

## 
## Call:
## lm(formula = EDCI_Total_Score_Post ~ EDCI_Total_Score_Pre + HS.GPA_IR + 
##     Consistency_Gain, data = consistency_model_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.816 -1.634 -0.157  1.796  6.421 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)           -1.1046     4.6327  -0.238  0.81265   
## EDCI_Total_Score_Pre   0.3837     0.1419   2.704  0.00970 **
## HS.GPA_IR              1.9018     1.2281   1.549  0.12864   
## Consistency_Gain       1.7144     0.6276   2.731  0.00904 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.897 on 44 degrees of freedom
## Multiple R-squared:  0.2939, Adjusted R-squared:  0.2458 
## F-statistic: 6.105 on 3 and 44 DF,  p-value: 0.001449

confint(consistency_model)

##                             2.5 %    97.5 %
## (Intercept)          -10.44117577 8.2319386
## EDCI_Total_Score_Pre   0.09772889 0.6696065
## HS.GPA_IR             -0.57322334 4.3769082
## Consistency_Gain       0.44940867 2.9792990

5.4 Figure: adjusted relationship between consistency and EDCI

The following plotting code is retained from the original file, apart from assigning the plot to an object so that it can be saved explicitly.

plot_data <- consistency_model_data

y_model <- lm(
  EDCI_Total_Score_Post ~
    EDCI_Total_Score_Pre +
    HS.GPA_IR,
  data = plot_data
)

x_model <- lm(
  Consistency_Gain ~
    EDCI_Total_Score_Pre +
    HS.GPA_IR,
  data = plot_data
)

plot_data <- plot_data %>%
  mutate(
    x_resid = resid(x_model),
    y_resid = resid(y_model)
  )

full_model <- lm(
  EDCI_Total_Score_Post ~
    EDCI_Total_Score_Pre +
    HS.GPA_IR +
    Consistency_Gain,
  data = plot_data
)

consistency_p <- summary(full_model)$coefficients[
  "Consistency_Gain",
  "Pr(>|t|)"
]

p_label <- paste0(
  "p = ",
  format.pval(consistency_p, digits = 2, eps = 0.001)
)

Consistency_EDCI <- ggplot(
  plot_data,
  aes(x = x_resid, y = y_resid)
) +
  geom_point(
    size = 3,
    alpha = 0.85,
    color = "#353C47"
  ) +
  geom_smooth(
    method = "lm",
    se = TRUE,
    linewidth = 1.3,
    color = "#29B897",
    fill = "#353C47",
    alpha = 0.45
  ) +
  annotate(
    "text",
    x = Inf,
    y = Inf,
    label = p_label,
    hjust = 1.15,
    vjust = 1.5,
    size = 5.5,
    fontface = "bold",
    color = "#29B897"
  ) +
  labs(
    x = "Consistency of Interests Gain",
    y = "Experimental Design Competency",
    caption = "Note. Relationship shown after adjusting for baseline experimental design competency and high school GPA."
  ) +
  theme_classic(base_size = 16) +
  theme(
    panel.background = element_rect(fill = "transparent", color = NA),
    plot.background = element_rect(fill = "transparent", color = NA),
    legend.background = element_rect(fill = "transparent", color = NA),
    legend.key = element_rect(fill = "transparent", color = NA),
    axis.line = element_line(color = "#353C47", linewidth = 0.8),
    axis.ticks = element_line(color = "#353C47"),
    axis.text = element_text(
      color = "#353C47",
      face = "bold"
    ),
    axis.title = element_text(
      color = "#353C47",
      face = "bold",
      size = 16
    ),
    plot.caption = element_text(
      size = 10,
      color = "gray40",
      hjust = 0,
      face = "italic",
      margin = margin(t = 8)
    )
  )

Consistency_EDCI

## `geom_smooth()` using formula = 'y ~ x'

ggsave(
  "Consistency_EDCI.png",
  plot = Consistency_EDCI,
  width = 4,
  height = 5,
  dpi = 600,
  bg = "transparent"
)

## `geom_smooth()` using formula = 'y ~ x'

5.5 Relative importance analysis

relative_importance <- calc.relimp(
  consistency_model,
  type = "lmg",
  rela = TRUE
)

relative_importance

## Response variable: EDCI_Total_Score_Post 
## Total response variance: 11.12722 
## Analysis based on 48 observations 
## 
## 3 Regressors: 
## EDCI_Total_Score_Pre HS.GPA_IR Consistency_Gain 
## Proportion of variance explained by model: 29.39%
## Metrics are normalized to sum to 100% (rela=TRUE). 
## 
## Relative importance metrics: 
## 
##                            lmg
## EDCI_Total_Score_Pre 0.4361402
## HS.GPA_IR            0.1285899
## Consistency_Gain     0.4352699
## 
## Average coefficients for different model sizes: 
## 
##                             1X      2Xs       3Xs
## EDCI_Total_Score_Pre 0.4162404 0.400223 0.3836677
## HS.GPA_IR            1.8651038 1.884843 1.9018424
## Consistency_Gain     1.8218110 1.767828 1.7143538

relative_importance$lmg

## EDCI_Total_Score_Pre            HS.GPA_IR     Consistency_Gain 
##            0.4361402            0.1285899            0.4352699

set.seed(123)

relative_boot <- boot.relimp(
  consistency_model,
  b = 1000,
  type = "lmg",
  rank = TRUE,
  diff = TRUE,
  rela = TRUE
)

booteval.relimp(
  relative_boot,
  bty = "perc"
)

## Response variable: EDCI_Total_Score_Post 
## Total response variance: 11.12722 
## Analysis based on 48 observations 
## 
## 3 Regressors: 
## EDCI_Total_Score_Pre HS.GPA_IR Consistency_Gain 
## Proportion of variance explained by model: 29.39%
## Metrics are normalized to sum to 100% (rela=TRUE). 
## 
## Relative importance metrics: 
## 
##                            lmg
## EDCI_Total_Score_Pre 0.4361402
## HS.GPA_IR            0.1285899
## Consistency_Gain     0.4352699
## 
## Average coefficients for different model sizes: 
## 
##                             1X      2Xs       3Xs
## EDCI_Total_Score_Pre 0.4162404 0.400223 0.3836677
## HS.GPA_IR            1.8651038 1.884843 1.9018424
## Consistency_Gain     1.8218110 1.767828 1.7143538
## 
##  
##  Confidence interval information ( 1000 bootstrap replicates, bty= perc ): 
## Relative Contributions with confidence intervals: 
##  
##                                         Lower  Upper
##                          percentage 0.95 0.95   0.95  
## EDCI_Total_Score_Pre.lmg 0.4361     ABC  0.0378 0.8673
## HS.GPA_IR.lmg            0.1286     ABC  0.0031 0.5433
## Consistency_Gain.lmg     0.4353     ABC  0.0343 0.8528
## 
## Letters indicate the ranks covered by bootstrap CIs. 
## (Rank bootstrap confidence intervals always obtained by percentile method) 
## CAUTION: Bootstrap confidence intervals can be somewhat liberal. 
## 
##  
##  Differences between Relative Contributions: 
##  
##                                                           Lower   Upper
##                                           difference 0.95 0.95    0.95   
## EDCI_Total_Score_Pre-HS.GPA_IR.lmg         0.3076         -0.3485  0.8208
## EDCI_Total_Score_Pre-Consistency_Gain.lmg  0.0009         -0.7581  0.7994
## HS.GPA_IR-Consistency_Gain.lmg            -0.3067         -0.8024  0.3497
## 
## * indicates that CI for difference does not include 0. 
## CAUTION: Bootstrap confidence intervals can be somewhat liberal.

5.6 Figure: relative contribution to explained variance

The following plotting code is retained from the original file.

importance_df <- tibble(
  Predictor = factor(
    c(
      "Baseline\nCompetency",
      "Growth in\nConsistency",
      "High School\nGPA"
    ),
    levels = c(
      "Baseline\nCompetency",
      "Growth in\nConsistency",
      "High School\nGPA"
    )
  ),
  Importance = c(43.6, 43.5, 12.8)
)

Relative_Importance_EDCI <- ggplot(
  importance_df,
  aes(
    x = Predictor,
    y = Importance
  )
) +
  geom_col(
    width = 0.6,
    fill = "#353C47"
  ) +
  geom_text(
    aes(label = paste0(round(Importance), "%")),
    vjust = -0.6,
    color = "#29B897",
    fontface = "bold",
    size = 6
  ) +
  scale_y_continuous(
    limits = c(0, 60),
    breaks = seq(0, 60, 20),
    expand = expansion(mult = c(0, 0.05))
  ) +
  labs(
    y = "Contribution to\nExplained Variance (%)",
    x = NULL,
    caption = "Relative importance estimates from regression model."
  ) +
  theme_classic(base_size = 18) +
  theme(
    panel.background = element_rect(fill = "transparent", color = NA),
    plot.background = element_rect(fill = "transparent", color = NA),
    axis.text.x = element_text(
      face = "bold",
      color = "#353C47",
      size = 12
    ),
    axis.text.y = element_text(
      face = "bold",
      color = "#353C47"
    ),
    axis.title.y = element_text(
      face = "bold",
      color = "#353C47"
    ),
    axis.line.x = element_blank(),
    axis.ticks.x = element_blank(),
    plot.caption = element_text(
      hjust = 0,
      size = 9,
      face = "italic",
      color = "gray40"
    )
  )

Relative_Importance_EDCI

ggsave(
  "Relative_Importance_EDCI.png",
  plot = Relative_Importance_EDCI,
  width = 4.5,
  height = 6.2,
  dpi = 600,
  bg = "transparent"
)

6. Preparedness gaps

6.1 Regression models

edci_preparedness_model <- lm(
  EDCI_Gain ~ EDCI_Total_Score_Pre,
  data = paired_wide
)

rsq_preparedness_model <- lm(
  RSQ_Gain ~ RSQ_Total_Pre,
  data = paired_wide
)

step_preparedness_model <- lm(
  STEP_Gain ~ STEP_Total_Pre,
  data = paired_wide
)

summary(edci_preparedness_model)

## 
## Call:
## lm(formula = EDCI_Gain ~ EDCI_Total_Score_Pre, data = paired_wide)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7140 -1.8256 -0.1217  2.2473  6.2473 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            6.7527     1.0572   6.387 4.06e-08 ***
## EDCI_Total_Score_Pre  -0.7039     0.1201  -5.862 2.83e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.162 on 54 degrees of freedom
## Multiple R-squared:  0.3889, Adjusted R-squared:  0.3776 
## F-statistic: 34.37 on 1 and 54 DF,  p-value: 2.832e-07

summary(rsq_preparedness_model)

## 
## Call:
## lm(formula = RSQ_Gain ~ RSQ_Total_Pre, data = paired_wide)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.05025 -0.32036  0.01007  0.28045  1.20265 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.0351     0.2583   4.007 0.000197 ***
## RSQ_Total_Pre  -0.3391     0.1166  -2.909 0.005317 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5254 on 52 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:   0.14,  Adjusted R-squared:  0.1235 
## F-statistic: 8.465 on 1 and 52 DF,  p-value: 0.005317

summary(step_preparedness_model)

## 
## Call:
## lm(formula = STEP_Gain ~ STEP_Total_Pre, data = paired_wide)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.02810 -0.19274  0.03976  0.25297  0.97333 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.9453     0.5236   3.715 0.000483 ***
## STEP_Total_Pre  -0.4600     0.1271  -3.618 0.000654 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.45 on 54 degrees of freedom
## Multiple R-squared:  0.1951, Adjusted R-squared:  0.1802 
## F-statistic: 13.09 on 1 and 54 DF,  p-value: 0.0006537

6.2 EDCI preparedness-gap figure

The following block is copied directly from the edited source file.

paired_wide <- paired_wide %>%
  mutate(
    EDCI_Quartile = ntile(EDCI_Total_Score_Pre, 4),
    EDCI_Quartile_Label = factor(
      EDCI_Quartile,
      levels = 1:4,
      labels = c("Lowest", "Q2", "Q3", "Highest")
    )
  )

fig2 <- ggplot(
  paired_wide,
  aes(
    x = factor(EDCI_Quartile),
    y = EDCI_Gain.x
  )
) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_boxplot(outlier.shape = NA, width = 0.6) +
  geom_jitter(width = 0.12, alpha = 0.55, size = 2) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  theme_classic(base_size = 13) +
  labs(
    x = "Baseline EDCI Quartile",
    y = "EDCI Gain",
    title = "Students with the lowest baseline competency showed the largest gains"
  )

fig2

ggsave("Figure2_EDCI_Baseline_Quartile.png", fig2, width = 6.5, height = 4.5, dpi = 600)

6.3 RSQ preparedness-gap figure

The following block is copied directly from the edited source file.

# ============================================================
# FIGURE 2B. RSQ gain by baseline RSQ quartile
# ============================================================

library(dplyr)
library(ggplot2)

paired_wide <- paired_wide %>%
  mutate(
    RSQ_Quartile = ntile(
      RSQ_Total_Pre,
      4
    )
  )


rsq_quartile_summary <- paired_wide %>%
  group_by(RSQ_Quartile) %>%
  summarise(
    N = n(),
    Mean_Pre = mean(RSQ_Total_Pre, na.rm = TRUE),
    Mean_Gain = mean(RSQ_Gain.x, na.rm = TRUE),
    SD_Gain = sd(RSQ_Gain.x, na.rm = TRUE),
    Mean_Post = mean(RSQ_Total_Post, na.rm = TRUE),
    .groups = "drop"
  )

rsq_quartile_summary

## # A tibble: 5 × 6
##   RSQ_Quartile     N Mean_Pre Mean_Gain SD_Gain Mean_Post
##          <int> <int>    <dbl>     <dbl>   <dbl>     <dbl>
## 1            1    14     1.45    0.669    0.526      2.12
## 2            2    14     1.95    0.402    0.547      2.35
## 3            3    13     2.16    0.0736   0.459      2.23
## 4            4    13     3.03    0.0728   0.525      3.10
## 5           NA     2   NaN     NaN       NA          2.29

fig2b <- ggplot(
  paired_wide,
  aes(
    x = factor(RSQ_Quartile),
    y = RSQ_Gain.x
  )
) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed",
    linewidth = 0.6,
    color = "gray40"
  ) +
  geom_boxplot(
    outlier.shape = NA,
    width = 0.6
  ) +
  geom_jitter(
    width = 0.12,
    alpha = 0.55,
    size = 2
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3
  ) +
  theme_classic(base_size = 13) +
  labs(
    x = "Baseline Research Self-Efficacy Quartile",
    y = "Research Self-Efficacy Gain"
  )

fig2b

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(
  "Figure2B_RSQ_Gain_by_Baseline_Quartile.png",
  fig2b,
  width = 6,
  height = 4.5,
  dpi = 600
)

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(
  "Figure2B_RSQ_Gain_by_Baseline_Quartile.pdf",
  fig2b,
  width = 6,
  height = 4.5
)

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

paired_wide <- paired_wide %>%
  mutate(
    STEP_Quartile = ntile(STEP_Total_Pre, 4),
    STEP_Quartile_Label = factor(
      STEP_Quartile,
      levels = 1:4,
      labels = c("Lowest", "Q2", "Q3", "Highest")
    )
  )

fig_step <- ggplot(
  paired_wide,
  aes(
    x = STEP_Quartile_Label,
    y = STEP_Gain.x
  )
) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed"
  ) +
  geom_boxplot(
    outlier.shape = NA,
    width = 0.6
  ) +
  geom_jitter(
    width = 0.12,
    alpha = 0.55,
    size = 2
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    size = 3
  ) +
  theme_classic(base_size = 13) +
  labs(
    x = "Baseline STEP-U Quartile",
    y = "STEP-U Gain",
    title = "Students with the lowest baseline science perceptions showed the largest gains"
  )

fig_step

step_gap_model <- lm(
  STEP_Gain ~ STEP_Total_Pre,
  data = paired_wide
)

summary(step_gap_model)

## 
## Call:
## lm(formula = STEP_Gain ~ STEP_Total_Pre, data = paired_wide)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.02810 -0.19274  0.03976  0.25297  0.97333 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.9453     0.5236   3.715 0.000483 ***
## STEP_Total_Pre  -0.4600     0.1271  -3.618 0.000654 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.45 on 54 degrees of freedom
## Multiple R-squared:  0.1951, Adjusted R-squared:  0.1802 
## F-statistic: 13.09 on 1 and 54 DF,  p-value: 0.0006537

Figure 3: (A) Distribution of total skills mastered, (B) Mean attempts required versus perseverance

7. Gender differences in change in grit

# # ============================================================
# # Gender differences in change in grit
# # ============================================================
# 
# library(dplyr)
# library(tidyr)
# library(lme4)
# library(lmerTest)
# library(emmeans)
# library(broom.mixed)
# library(ggplot2)
# library(patchwork)
# 
# grit_gender_long <- dat_paired %>%
#   dplyr::select(
#     StudentID,
#     Timepoint,
#     Gender,
#     Grit_Total,
#     Grit_Perseverance,
#     Grit_Consistency
#   ) %>%
#   filter(
#     Timepoint %in% c("Pre", "Post"),
#     Gender %in% c("Female", "Male")
#   ) %>%
#   mutate(
#     Timepoint = factor(
#       Timepoint,
#       levels = c("Pre", "Post")
#     ),
#     Gender = factor(
#       Gender,
#       levels = c("Female", "Male")
#     )
#   )
# 
# grit_gender_long %>%
#   count(Gender, Timepoint)
# 
# model_grit_gender <- lmer(
#   Grit_Total ~
#     Timepoint * Gender +
#     (1 | StudentID),
#   data = grit_gender_long,
#   REML = FALSE
# )
# 
# anova(model_grit_gender)
# summary(model_grit_gender)
# 
# model_perseverance_gender <- lmer(
#   Grit_Perseverance ~
#     Timepoint * Gender +
#     (1 | StudentID),
#   data = grit_gender_long,
#   REML = FALSE
# )
# 
# anova(model_perseverance_gender)
# summary(model_perseverance_gender)
# 
# model_consistency_gender <- lmer(
#   Grit_Consistency ~
#     Timepoint * Gender +
#     (1 | StudentID),
#   data = grit_gender_long,
#   REML = FALSE
# )
# 
# anova(model_consistency_gender)
# summary(model_consistency_gender)
# 
# emm_grit <- emmeans(
#   model_grit_gender,
#   ~ Timepoint | Gender
# )
# 
# emm_grit
# pairs(emm_grit, adjust = "holm")
# 
# emm_perseverance <- emmeans(
#   model_perseverance_gender,
#   ~ Timepoint | Gender
# )
# 
# emm_perseverance
# pairs(emm_perseverance, adjust = "holm")
# 
# emm_consistency <- emmeans(
#   model_consistency_gender,
#   ~ Timepoint | Gender
# )
# 
# emm_consistency
# pairs(emm_consistency, adjust = "holm")
# 
# contrast(
#   emmeans(
#     model_grit_gender,
#     ~ Timepoint * Gender
#   ),
#   interaction = "revpairwise",
#   adjust = "holm"
# )
# 
# contrast(
#   emmeans(
#     model_perseverance_gender,
#     ~ Timepoint * Gender
#   ),
#   interaction = "revpairwise",
#   adjust = "holm"
# )
# 
# contrast(
#   emmeans(
#     model_consistency_gender,
#     ~ Timepoint * Gender
#   ),
#   interaction = "revpairwise",
#   adjust = "holm"
# )
# 
# grit_gender_summary <- grit_gender_long %>%
#   pivot_longer(
#     cols = c(
#       Grit_Total,
#       Grit_Perseverance,
#       Grit_Consistency
#     ),
#     names_to = "Grit_Construct",
#     values_to = "Score"
#   ) %>%
#   group_by(
#     Grit_Construct,
#     Gender,
#     Timepoint
#   ) %>%
#   summarise(
#     n = sum(!is.na(Score)),
#     Mean = mean(Score, na.rm = TRUE),
#     SD = sd(Score, na.rm = TRUE),
#     SE = SD / sqrt(n),
#     .groups = "drop"
#   )
# 
# grit_gender_summary
# 
# grit_colors <- c(
#   "Female" = "#F8766D",
#   "Male" = "#00BFC4"
# )
# 
# make_grit_gender_plot <- function(data, outcome, y_label, title) {
# 
#   ggplot(
#     data,
#     aes(
#       x = Timepoint,
#       y = .data[[outcome]],
#       group = Gender,
#       color = Gender
#     )
#   ) +
#     stat_summary(
#       fun = mean,
#       geom = "line",
#       linewidth = 1.2
#     ) +
#     stat_summary(
#       fun = mean,
#       geom = "point",
#       size = 3.5
#     ) +
#     stat_summary(
#       fun.data = mean_se,
#       geom = "errorbar",
#       width = 0.12,
#       linewidth = 0.8
#     ) +
#     scale_color_manual(values = grit_colors) +
#     labs(
#       title = title,
#       x = "Timepoint",
#       y = y_label,
#       color = "Gender"
#     ) +
#     theme_classic(base_size = 13) +
#     theme(
#       plot.title = element_text(
#         face = "bold",
#         hjust = 0
#       ),
#       axis.title = element_text(
#         face = "bold"
#       )
#     )
# }
# 
# plot_grit_total <- make_grit_gender_plot(
#   grit_gender_long,
#   outcome = "Grit_Total",
#   y_label = "Mean Grit Total",
#   title = "Grit: Overall Score"
# )
# 
# plot_grit_perseverance <- make_grit_gender_plot(
#   grit_gender_long,
#   outcome = "Grit_Perseverance",
#   y_label = "Mean Grit Perseverance",
#   title = "Grit: Perseverance"
# ) +
#   theme(legend.position = "none")
# 
# plot_grit_consistency <- make_grit_gender_plot(
#   grit_gender_long,
#   outcome = "Grit_Consistency",
#   y_label = "Mean Grit Consistency",
#   title = "Grit: Consistency"
# ) +
#   theme(legend.position = "none")
# 
# grit_gender_figure <-
#   plot_grit_total /
#   (plot_grit_perseverance + plot_grit_consistency) +
#   plot_layout(heights = c(1.15, 1))
# 
# grit_gender_figure
# 
# anova(lm(
#   EDCI_Total_Score_Post ~
#     EDCI_Total_Score_Pre +
#     HS.GPA_IR +
#     Consistency_Gain * Gender,
#   data = paired_wide
# ))

Figure 1:

8. Export analysis tables

# dir.create("Paper_Output", showWarnings = FALSE)
# 
# write.csv(
#   broom::tidy(consistency_model, conf.int = TRUE),
#   "Paper_Output/Consistency_EDCI_Model.csv",
#   row.names = FALSE
# )
# 
# write.csv(
#   data.frame(
#     Predictor = names(relative_importance$lmg),
#     Relative_Importance = as.numeric(relative_importance$lmg)
#   ),
#   "Paper_Output/Relative_Importance_EDCI.csv",
#   row.names = FALSE
# )
# 
# write.csv(
#   bind_rows(
#     broom::tidy(edci_preparedness_model, conf.int = TRUE) %>%
#       mutate(Outcome = "EDCI"),
#     broom::tidy(rsq_preparedness_model, conf.int = TRUE) %>%
#       mutate(Outcome = "RSQ"),
#     broom::tidy(step_preparedness_model, conf.int = TRUE) %>%
#       mutate(Outcome = "STEP-U")
#   ),
#   "Paper_Output/Preparedness_Gap_Models.csv",
#   row.names = FALSE
# )

9. Session information

sessionInfo()

## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.7.4
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] lubridate_1.9.5   purrr_1.2.2       tibble_3.3.1      tidyverse_2.0.0  
##  [5] relaimpo_2.2-7    mitools_2.4       survey_4.5        survival_3.8-3   
##  [9] Matrix_1.7-4      boot_1.3-32       MASS_7.3-65       broom_1.0.12     
## [13] scales_1.4.0      patchwork_1.3.2   ggalluvial_0.12.6 ggplot2_4.0.2    
## [17] forcats_1.0.1     stringr_1.6.0     readr_2.2.0       tidyr_1.3.2      
## [21] dplyr_1.2.1      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       xfun_0.56          bslib_0.10.0       lattice_0.22-7    
##  [5] tzdb_0.5.0         vctrs_0.7.1        tools_4.5.2        generics_0.1.4    
##  [9] pkgconfig_2.0.3    RColorBrewer_1.1-3 S7_0.2.1           lifecycle_1.0.5   
## [13] compiler_4.5.2     farver_2.1.2       textshaping_1.0.5  htmltools_0.5.9   
## [17] sass_0.4.10        yaml_2.3.12        pillar_1.11.1      jquerylib_0.1.4   
## [21] cachem_1.1.0       nlme_3.1-168       tidyselect_1.2.1   digest_0.6.39     
## [25] stringi_1.8.7      labeling_0.4.3     splines_4.5.2      fastmap_1.2.0     
## [29] cli_3.6.5          magrittr_2.0.4     utf8_1.2.6         corpcor_1.6.10    
## [33] withr_3.0.2        backports_1.5.0    timechange_0.4.0   rmarkdown_2.30    
## [37] otel_0.2.0         ragg_1.5.2         hms_1.1.4          evaluate_1.0.5    
## [41] knitr_1.51         viridisLite_0.4.3  mgcv_1.9-3         rlang_1.3.0       
## [45] Rcpp_1.1.2         glue_1.8.1         DBI_1.3.0          rstudioapi_0.18.0 
## [49] jsonlite_2.0.0     R6_2.6.1           systemfonts_1.3.2

Iterative Mastery-Based Laboratory Assessment: Paper Analysis

Abby Beatty

2026-08-05

Purpose

Purpose

1. Setup

2. Load and prepare data

2.1 Import and merge

2.2 Academic preparation variables

2.3 Score grit

2.4 Score STEP-U, RSQ, and EDCI

2.5 Merge mastery-assessment outcomes

2.6 Create paired student-level dataset

3. Overall learning outcomes

3.1 Primary paired tests

3.2 Figure 1: Student gains

4. Mastery of laboratory skills

4.1 Descriptive mastery statistics

4.2 Figure: Distribution of mastery and mean attempts

4.3 Figure: Mastery progression

4.4 Confidence across mastery attempts

5. Affective factors associated with student success

5.1 Perseverance and first-pass mastery

5.2 Figure: First-pass mastery by perseverance quartile

5.3 Consistency of interests and EDCI outcome

5.4 Figure: adjusted relationship between consistency and EDCI

5.5 Relative importance analysis

5.6 Figure: relative contribution to explained variance

6. Preparedness gaps

6.1 Regression models

6.2 EDCI preparedness-gap figure

6.3 RSQ preparedness-gap figure

7. Gender differences in change in grit

8. Export analysis tables

9. Session information