This script reproduces the student composite-score analyses reported in Tables 2 and 3 of the manuscript. The workflow is intentionally limited to the analyses required for those tables and the associated post hoc comparisons.
required_packages <- c("dplyr", "readxl", "psych", "lavaan", "knitr")
missing_packages <- required_packages[
!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages) > 0) {
stop(
"Install the following packages before running this analysis: ",
paste(missing_packages, collapse = ", ")
)
}
options(digits = 6)
# Bonferroni-adjusted significance threshold used in the manuscript for
# the three composite-score comparisons within each subgroup analysis.
alpha_bonf <- 0.05 / 3find_student_file <- function() {
file_name <- "AI4DS_Student_Survey_October_5_ 2025.xlsx"
candidates <- c(
file_name,
file.path("data", file_name),
file.path("..", "data", file_name)
)
found <- candidates[file.exists(candidates)]
if (length(found) == 0) {
stop(
"Student survey Excel file not found. Expected the file in the repo root, ",
"data/, or ../data/ relative to the working directory."
)
}
found[[1]]
}
load_student_data <- function(path, min_duration = 120, min_progress = 25) {
# The Qualtrics export contains short variable names in row 1,
# question labels in row 2, and respondent data beginning in row 3.
header <- readxl::read_excel(path, n_max = 1, col_names = FALSE)
short_names <- as.character(unlist(header[1, ], use.names = FALSE))
short_names <- short_names[!is.na(short_names) & short_names != ""]
student_df <- readxl::read_excel(
path,
skip = 2,
col_names = short_names,
.name_repair = "minimal"
)
student_df |>
dplyr::filter(
`Duration (in seconds)` >= min_duration,
Progress >= min_progress
)
}
student_data_path <- find_student_file()
student_data <- load_student_data(student_data_path)
cat("Data file:", normalizePath(student_data_path), "\n")## Data file: C:\Users\amhasan1\OneDrive - North Carolina A&T State University\Documents\GitHub\AI4DS\Data & RCode\data\AI4DS_Student_Survey_October_5_ 2025.xlsxcat("Data MD5:", unname(tools::md5sum(student_data_path)), "\n")## Data MD5: 81957b81331eabab6a7e93ca8c39c0f7cat("Analytic sample after quality filters:", nrow(student_data), "students\n")## Analytic sample after quality filters: 119 students# The distributed repository dataset is the finalized analytic file.
# The manuscript reports n = 119 valid student responses.
if (nrow(student_data) != 119L) {
stop(
"Unexpected analytic sample size. Expected 119 students but found ",
nrow(student_data),
". Verify the data-file version and filtering criteria."
)
}
required_variables <- c(
"Q3", "Q5",
paste0("Q10_", 1:5),
paste0("Q18_", 1:5),
"Q12_1", "Q12_2",
"Q21_1", "Q21_2", "Q21_4"
)
missing_variables <- setdiff(required_variables, names(student_data))
if (length(missing_variables) > 0) {
stop(
"Required variables are missing from the input file: ",
paste(missing_variables, collapse = ", ")
)
}# In the initial survey version, "Graduate Student" was inadvertently
# omitted from Q5. Review of respondents' reported majors established that
# respondents with missing Q5 values in this dataset were graduate students.
# These missing Q5 values are therefore recoded as "Graduate Student".
student_data <- student_data |>
dplyr::mutate(
Q5 = dplyr::coalesce(as.character(Q5), "Graduate Student"),
class_group = dplyr::case_when(
Q5 %in% c("Freshman", "Sophomore") ~ "Lower-Division",
Q5 %in% c("Junior", "Senior") ~ "Upper-Division",
Q5 == "Graduate Student" ~ "Graduate",
TRUE ~ NA_character_
),
class_group = factor(
class_group,
levels = c("Lower-Division", "Upper-Division", "Graduate")
)
)
score_familiarity <- function(x) {
dplyr::recode(
as.character(x),
"Not at All Familiar" = 1,
"Not Very Familiar" = 2,
"Neutral" = 3,
"Somewhat Familiar" = 4,
"Very Familiar" = 5,
.default = NA_real_
)
}
score_awareness <- function(x) {
dplyr::recode(
as.character(x),
"Not at All Aware" = 1,
"Not Very Aware" = 2,
"Neutral" = 3,
"Somewhat Aware" = 4,
"Very Aware" = 5,
.default = NA_real_
)
}
score_agreement <- function(x) {
dplyr::recode(
as.character(x),
"Strongly Disagree" = 1,
"Disagree" = 2,
"Neutral" = 3,
"Agree" = 4,
"Strongly Agree" = 5,
.default = NA_real_
)
}
student_data <- student_data |>
dplyr::mutate(
dplyr::across(Q10_1:Q10_5, score_familiarity),
dplyr::across(Q18_1:Q18_5, score_awareness),
dplyr::across(
c(Q12_1, Q12_2, Q21_1, Q21_2, Q21_4),
score_agreement
)
)The three composite scores are arithmetic means of the available
selected items for each respondent. A composite is set to
NA only when all five constituent items are missing. The
Perceptions composite uses the five conceptually selected items Q12_1,
Q12_2, Q21_1, Q21_2, and Q21_4.
student_data <- student_data |>
dplyr::rowwise() |>
dplyr::mutate(
Q10_comp = ifelse(
all(is.na(dplyr::c_across(Q10_1:Q10_5))),
NA_real_,
mean(dplyr::c_across(Q10_1:Q10_5), na.rm = TRUE)
),
Q18_comp = ifelse(
all(is.na(dplyr::c_across(Q18_1:Q18_5))),
NA_real_,
mean(dplyr::c_across(Q18_1:Q18_5), na.rm = TRUE)
),
Qper_comp = ifelse(
all(is.na(dplyr::c_across(c(Q12_1, Q12_2, Q21_1, Q21_2, Q21_4)))),
NA_real_,
mean(
dplyr::c_across(c(Q12_1, Q12_2, Q21_1, Q21_2, Q21_4)),
na.rm = TRUE
)
)
) |>
dplyr::ungroup()
composite_summary <- dplyr::bind_rows(
data.frame(
Composite = "Familiarity",
n = sum(!is.na(student_data$Q10_comp)),
Mean = mean(student_data$Q10_comp, na.rm = TRUE),
SD = stats::sd(student_data$Q10_comp, na.rm = TRUE)
),
data.frame(
Composite = "Perceptions",
n = sum(!is.na(student_data$Qper_comp)),
Mean = mean(student_data$Qper_comp, na.rm = TRUE),
SD = stats::sd(student_data$Qper_comp, na.rm = TRUE)
),
data.frame(
Composite = "Awareness of Limitations",
n = sum(!is.na(student_data$Q18_comp)),
Mean = mean(student_data$Q18_comp, na.rm = TRUE),
SD = stats::sd(student_data$Q18_comp, na.rm = TRUE)
)
) |>
dplyr::mutate(
Mean = round(Mean, 3),
SD = round(SD, 3)
)
knitr::kable(composite_summary, caption = "Overall composite-score summary.")| Composite | n | Mean | SD |
|---|---|---|---|
| Familiarity | 109 | 2.800 | 0.973 |
| Perceptions | 108 | 3.814 | 0.578 |
| Awareness of Limitations | 102 | 4.112 | 0.866 |
CFA is conducted after the composite-score variables have been constructed. The CFA models use the five constituent items for each construct and explicitly use listwise deletion for missing CFA indicators. Cronbach’s alpha is calculated from the same item sets.
model_familiarity <- '
Familiarity =~ Q10_1 + Q10_2 + Q10_3 + Q10_4 + Q10_5
'
model_perceptions <- '
Perceptions =~ Q12_1 + Q12_2 + Q21_1 + Q21_2 + Q21_4
'
model_awareness <- '
Awareness =~ Q18_1 + Q18_2 + Q18_3 + Q18_4 + Q18_5
'
fit_familiarity <- lavaan::cfa(
model_familiarity,
data = student_data,
missing = "listwise"
)
fit_perceptions <- lavaan::cfa(
model_perceptions,
data = student_data,
missing = "listwise"
)
fit_awareness <- lavaan::cfa(
model_awareness,
data = student_data,
missing = "listwise"
)
alpha_familiarity <- psych::alpha(
student_data |> dplyr::select(Q10_1:Q10_5)
)
alpha_perceptions <- psych::alpha(
student_data |>
dplyr::select(Q12_1, Q12_2, Q21_1, Q21_2, Q21_4)
)
alpha_awareness <- psych::alpha(
student_data |> dplyr::select(Q18_1:Q18_5)
)
make_cfa_row <- function(label, fit, alpha_object) {
data.frame(
Composite = label,
Chi_Square = unname(lavaan::fitMeasures(fit, "chisq")),
df = as.integer(unname(lavaan::fitMeasures(fit, "df"))),
P_value = unname(lavaan::fitMeasures(fit, "pvalue")),
CFI = unname(lavaan::fitMeasures(fit, "cfi")),
TLI = unname(lavaan::fitMeasures(fit, "tli")),
RMSEA = unname(lavaan::fitMeasures(fit, "rmsea")),
Raw_Alpha = alpha_object$total$raw_alpha,
Standard_Alpha = alpha_object$total$std.alpha,
stringsAsFactors = FALSE
)
}
cfa_results <- dplyr::bind_rows(
make_cfa_row("Familiarity", fit_familiarity, alpha_familiarity),
make_cfa_row("Perceptions", fit_perceptions, alpha_perceptions),
make_cfa_row("Awareness of Limitations", fit_awareness, alpha_awareness)
)Table_2 <- cfa_results |>
dplyr::transmute(
Composite,
`Chi-square Stat (p-value)` = sprintf("%.3f (%.3f)", Chi_Square, P_value),
CFI = sprintf("%.3f", CFI),
TLI = sprintf("%.3f", TLI),
RMSEA = sprintf("%.3f", RMSEA),
`Raw Alpha` = sprintf("%.3f", Raw_Alpha),
`Standard Alpha` = sprintf("%.3f", Standard_Alpha)
)
knitr::kable(
Table_2,
align = c("l", rep("c", 6)),
caption = "Table 2. Model fit statistics from confirmatory factor analysis."
)| Composite | Chi-square Stat (p-value) | CFI | TLI | RMSEA | Raw Alpha | Standard Alpha |
|---|---|---|---|---|---|---|
| Familiarity | 8.875 (0.114) | 0.969 | 0.939 | 0.086 | 0.756 | 0.753 |
| Perceptions | 8.309 (0.140) | 0.971 | 0.942 | 0.081 | 0.766 | 0.767 |
| Awareness of Limitations | 3.320 (0.651) | 1.000 | 1.017 | 0.000 | 0.859 | 0.860 |
For validation, the CFA degrees of freedom and number of observations used by each listwise CFA are retained below.
cfa_audit <- data.frame(
Composite = c("Familiarity", "Perceptions", "Awareness of Limitations"),
df = c(
lavaan::fitMeasures(fit_familiarity, "df"),
lavaan::fitMeasures(fit_perceptions, "df"),
lavaan::fitMeasures(fit_awareness, "df")
),
N_used = c(
lavaan::lavInspect(fit_familiarity, "nobs"),
lavaan::lavInspect(fit_perceptions, "nobs"),
lavaan::lavInspect(fit_awareness, "nobs")
)
)
knitr::kable(cfa_audit, caption = "CFA audit information (not part of manuscript Table 2).")| Composite | df | N_used |
|---|---|---|
| Familiarity | 5 | 105 |
| Perceptions | 5 | 100 |
| Awareness of Limitations | 5 | 100 |
The subgroup analyses are intentionally performed only after all composite-score and CFA calculations are complete. Academic-level comparisons use the full eligible analytic sample. Gender comparisons use only respondents who selected Female or Male.
score_variables <- c(
Familiarity = "Q10_comp",
Perceptions = "Qper_comp",
`Awareness of Limitations` = "Q18_comp"
)
class_summary_list <- lapply(names(score_variables), function(composite_name) {
score_name <- score_variables[[composite_name]]
student_data |>
dplyr::filter(!is.na(class_group)) |>
dplyr::group_by(class_group) |>
dplyr::summarise(
n = sum(!is.na(.data[[score_name]])),
Mean = mean(.data[[score_name]], na.rm = TRUE),
SD = stats::sd(.data[[score_name]], na.rm = TRUE),
.groups = "drop"
) |>
dplyr::mutate(Composite = composite_name, .before = 1)
})
class_summary <- dplyr::bind_rows(class_summary_list)
class_anova <- lapply(names(score_variables), function(composite_name) {
score_name <- score_variables[[composite_name]]
fit <- stats::aov(
stats::reformulate("class_group", response = score_name),
data = student_data
)
data.frame(
Composite = composite_name,
F = unname(summary(fit)[[1]]["class_group", "F value"]),
P_value = unname(summary(fit)[[1]]["class_group", "Pr(>F)"]),
stringsAsFactors = FALSE
)
}) |>
dplyr::bind_rows()
class_anova## Composite F P_value
## 1 Familiarity 8.6524132 0.000330573
## 2 Perceptions 0.0862061 0.917469972
## 3 Awareness of Limitations 5.6971339 0.004550815get_tukey <- function(score_name, composite_name) {
fit <- stats::aov(
stats::reformulate("class_group", response = score_name),
data = student_data
)
tukey <- as.data.frame(stats::TukeyHSD(fit)$class_group)
tukey$Comparison <- rownames(tukey)
rownames(tukey) <- NULL
tukey |>
dplyr::transmute(
Composite = composite_name,
Comparison,
Difference = round(diff, 3),
Lower_CI = round(lwr, 3),
Upper_CI = round(upr, 3),
Adjusted_P = round(`p adj`, 3)
)
}
# The manuscript reports post hoc results only for constructs with a
# significant omnibus class-level ANOVA.
tukey_results <- dplyr::bind_rows(
get_tukey("Q10_comp", "Familiarity"),
get_tukey("Q18_comp", "Awareness of Limitations")
)
knitr::kable(
tukey_results,
caption = "Tukey HSD post hoc comparisons for significant class-level ANOVAs."
)| Composite | Comparison | Difference | Lower_CI | Upper_CI | Adjusted_P |
|---|---|---|---|---|---|
| Familiarity | Upper-Division-Lower-Division | 0.425 | -0.085 | 0.935 | 0.122 |
| Familiarity | Graduate-Lower-Division | 1.013 | 0.430 | 1.597 | 0.000 |
| Familiarity | Graduate-Upper-Division | 0.588 | 0.084 | 1.092 | 0.018 |
| Awareness of Limitations | Upper-Division-Lower-Division | 0.548 | 0.068 | 1.028 | 0.021 |
| Awareness of Limitations | Graduate-Lower-Division | 0.752 | 0.194 | 1.310 | 0.005 |
| Awareness of Limitations | Graduate-Upper-Division | 0.204 | -0.276 | 0.684 | 0.572 |
# Restrict only the gender-comparison analysis to respondents who selected
# Female or Male. The master student_data object remains unchanged.
gender_data <- student_data |>
dplyr::filter(Q3 %in% c("Female", "Male")) |>
dplyr::mutate(Q3 = factor(Q3, levels = c("Female", "Male")))
gender_summary_list <- lapply(names(score_variables), function(composite_name) {
score_name <- score_variables[[composite_name]]
gender_data |>
dplyr::group_by(Q3) |>
dplyr::summarise(
n = sum(!is.na(.data[[score_name]])),
Mean = mean(.data[[score_name]], na.rm = TRUE),
SD = stats::sd(.data[[score_name]], na.rm = TRUE),
.groups = "drop"
) |>
dplyr::mutate(Composite = composite_name, .before = 1)
})
gender_summary <- dplyr::bind_rows(gender_summary_list)
gender_tests <- lapply(names(score_variables), function(composite_name) {
score_name <- score_variables[[composite_name]]
test <- stats::t.test(
stats::reformulate("Q3", response = score_name),
data = gender_data,
var.equal = FALSE
)
data.frame(
Composite = composite_name,
t = unname(test$statistic),
df = unname(test$parameter),
P_value = test$p.value,
stringsAsFactors = FALSE
)
}) |>
dplyr::bind_rows()
gender_tests## Composite t df P_value
## 1 Familiarity -0.724643 68.5186 0.471138143
## 2 Perceptions -1.453337 76.7824 0.150205253
## 3 Awareness of Limitations -3.575149 96.2851 0.000549378Table 3 is generated directly from the subgroup summaries and hypothesis tests above. All displayed numeric values are rounded to three decimal places. Asterisks identify p-values below the Bonferroni-adjusted threshold of 0.0167 used in the manuscript.
format_mean_sd <- function(mean_value, sd_value) {
sprintf("%.3f (%.3f)", mean_value, sd_value)
}
format_p <- function(p_value) {
if (is.na(p_value)) return("")
paste0(
sprintf("%.3f", p_value),
ifelse(p_value < alpha_bonf, "*", "")
)
}
extract_summary <- function(summary_df, composite_name, group_column, group_name) {
row <- summary_df[
summary_df$Composite == composite_name &
as.character(summary_df[[group_column]]) == group_name,
,
drop = FALSE
]
if (nrow(row) != 1) {
stop("Could not uniquely identify summary row for ", composite_name, " / ", group_name)
}
list(
n = row$n,
mean_sd = format_mean_sd(row$Mean, row$SD)
)
}
get_test_p <- function(test_df, composite_name) {
test_df$P_value[test_df$Composite == composite_name][1]
}
# The manuscript displays p-values once per comparison: on the Male row for
# Welch gender tests and on the Graduate row for three-level ANOVA tests.
make_table3_row <- function(group_name, comparison_type) {
if (comparison_type == "gender") {
source_df <- gender_summary
group_col <- "Q3"
p_source <- gender_tests
show_p <- group_name == "Male"
} else {
source_df <- class_summary
group_col <- "class_group"
p_source <- class_anova
show_p <- group_name == "Graduate"
}
fam <- extract_summary(source_df, "Familiarity", group_col, group_name)
per <- extract_summary(source_df, "Perceptions", group_col, group_name)
awa <- extract_summary(source_df, "Awareness of Limitations", group_col, group_name)
data.frame(
Group = group_name,
Familiarity_n = fam$n,
`Familiarity Mean (SD)` = fam$mean_sd,
`Familiarity P-value` = if (show_p) format_p(get_test_p(p_source, "Familiarity")) else "",
Perceptions_n = per$n,
`Perceptions Mean (SD)` = per$mean_sd,
`Perceptions P-value` = if (show_p) format_p(get_test_p(p_source, "Perceptions")) else "",
Awareness_n = awa$n,
`Awareness Mean (SD)` = awa$mean_sd,
`Awareness P-value` = if (show_p) format_p(get_test_p(p_source, "Awareness of Limitations")) else "",
check.names = FALSE,
stringsAsFactors = FALSE
)
}
# Retain the manuscript row order while keeping the analysis code itself in
# the requested order: class-level analysis first, then gender analysis.
Table_3 <- dplyr::bind_rows(
make_table3_row("Female", "gender"),
make_table3_row("Male", "gender"),
make_table3_row("Lower-Division", "class"),
make_table3_row("Upper-Division", "class"),
make_table3_row("Graduate", "class")
)
knitr::kable(
Table_3,
align = c("l", rep("c", 9)),
caption = "Table 3. Composite score comparisons by gender and class levels."
)| Group | Familiarity_n | Familiarity Mean (SD) | Familiarity P-value | Perceptions_n | Perceptions Mean (SD) | Perceptions P-value | Awareness_n | Awareness Mean (SD) | Awareness P-value |
|---|---|---|---|---|---|---|---|---|---|
| Female | 69 | 2.695 (0.921) | 70 | 3.767 (0.591) | 65 | 3.920 (0.940) | |||
| Male | 36 | 2.836 (0.961) | 0.471 | 36 | 3.933 (0.539) | 0.150 | 35 | 4.463 (0.576) | 0.001* |
| Lower-Division | 27 | 2.330 (0.748) | 28 | 3.811 (0.579) | 25 | 3.648 (0.947) | |||
| Upper-Division | 54 | 2.755 (0.914) | 54 | 3.796 (0.556) | 52 | 4.196 (0.818) | |||
| Graduate | 28 | 3.343 (1.037) | 0.000* | 26 | 3.854 (0.639) | 0.917 | 25 | 4.400 (0.716) | 0.005* |
Table 3 note. Gender comparisons use Welch’s two-sample t-test. Class-level comparisons use one-way ANOVA. The manuscript applies a Bonferroni-adjusted significance threshold of 0.05/3 = 0.0167 across the three composite-score tests within each subgroup analysis. Tukey HSD is used for post hoc class-level comparisons following significant omnibus ANOVAs.
sessionInfo()## R version 4.5.2 (2025-10-31 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=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## loaded via a namespace (and not attached):
## [1] jsonlite_2.0.0 dplyr_1.1.4 compiler_4.5.2 tidyselect_1.2.1
## [5] psych_2.5.6 parallel_4.5.2 jquerylib_0.1.4 yaml_2.3.10
## [9] fastmap_1.2.0 readxl_1.4.5 lattice_0.22-7 pbivnorm_0.6.0
## [13] R6_2.6.1 generics_0.1.4 knitr_1.50 tibble_3.3.0
## [17] bslib_0.9.0 pillar_1.11.1 rlang_1.1.6 cachem_1.1.0
## [21] xfun_0.54 quadprog_1.5-8 sass_0.4.10 cli_3.6.5
## [25] withr_3.0.2 magrittr_2.0.4 digest_0.6.38 grid_4.5.2
## [29] rstudioapi_0.17.1 lifecycle_1.0.4 nlme_3.1-168 lavaan_0.6-21
## [33] vctrs_0.6.5 mnormt_2.1.1 evaluate_1.0.5 glue_1.8.0
## [37] cellranger_1.1.0 stats4_4.5.2 rmarkdown_2.30 tools_4.5.2
## [41] pkgconfig_2.0.3 htmltools_0.5.8.1