Working memory (WM) is a limited-capacity system that temporarily stores, manipulates, and retrieves information necessary for ongoing cognitive processes. It is a core construct in cognitive psychology and neuroscience. Yet, its assessment remains methodologically inconsistent. Researchers employ a wide range of task combinations to assess this construct, including single-task measurements, employing multiple versions of the same paradigm (class of WM tasks), and administering diverse tasks drawn from different paradigms. This leads to divergent representations of WM across studies. This variability hinders comparability and limits the accumulation of a coherent evidence base.

A growing consensus suggests that WM should be assessed using a heterogeneous battery encompassing tasks from multiple paradigms. Such an approach captures broader WM-relevant variance and engages distinct functional components (e.g., simultaneous storage and processing, updating of goal-relevant information, and short-term binding), while minimizing paradigm-specific noise. Despite these recommendations, many recent studies continue to rely on single-task or single-paradigm assessments. This may be partly explain by the absence of an open, automated, and psychometrically validated WM battery.

To address this gap, we developed the OpenWMB, an automated working memory battery that includes seven tasks from three established paradigms. The full set of tasks is summarized below.

Paradigm Tasks
Complex span Reading Span, Operation Span, Symmetry Span
Updating Working Memory Updating Task, N-back Task
Binding Multimodal Span, Binding-and-Maintenance Task

These tasks were designed to capture the aforementioned facets of WM (simultaneous storage and processing, updating, and binding) and are freely available in multiple languages under an open license. The full battery can be accessed at https://github.com/MonteiroEtAl/OpenWMB.

In order to evaluate the psychometic properties of the OpenWMB, we conducted a large validation study in which 169 participants completed the entire battery in a single session. We examined three main aspects of psychometric quality. Internal consistency was assessed through reliability estimates for each task. Convergent validity was tested by analyzing the interrelations among the WM tasks at the latent level using exploratory and confirmatory factor analyses. Finally, concurrent validity was evaluated by modeling the association between a latent factor representing WM and another latent factor derived from three fluid intelligence (Gf) measures — Letter Series, Number Series, and Raven’s Advanced Progressive Matrices — using structural equation modeling (SEM).

This report presents the statistical procedures and R code used to examine the psychometric properties of the OpenWMB.

0. Libraries used to compute the analyses

library(purrr)
library(DT)
library(semPlot)
library(dplyr)
library(semTools)
library(MVN)
library(mice)
library(readxl)
library(tibble)
library(openxlsx)
library(simsem)
library(psych)
library(sjmisc)
library(tidyr)
library(knitr)

1. Importing and Structuring the Raw Data

In the first subsection, we import the raw data for the the cognitive tasks from the Excel workbook BD_Original.xlsx. Each worksheet in the file corresponds to a specific task. The first seven sheets contain data regarding the WM tasks (Reading Span, N-Back Task, Operation Span, Binding and Maintenance Task, Multimodal Span, Symmetry Span, and WM Updating Task). The last three sheets contain the fluid intelligence (Gf) tasks (Number Series, Raven’s Advanced Progressive Matrices, and Letter Series). Each worksheet is read into R as an individual dataframe for subsequent data processing and analysis.

In the second subsection, we compute a single performance score for each participant on every WM and Gf task. Performance scores were calculated as proportions (nr of correct responses / total number of trials). We then create dataframes that contain both sociodemographic information and task scores to generate participant-level datasets for both WM and Gf. Finally, we display these datasets as interactive tables.

1.1. Importing Raw Data from Excel Workbooks

## Working Memory (WM) tasks
# Each dataframe corresponds to a different WM task
df_RS <- read_excel("BD_Original.xlsx", sheet = 'Reading_Span')
df_NB <- read_excel("BD_Original.xlsx", sheet = 'N_Back')
df_OS <- read_excel("BD_Original.xlsx", sheet = 'Operation_Span')
df_BT <- read_excel("BD_Original.xlsx", sheet = 'Binding_Task')
df_MS <- read_excel("BD_Original.xlsx", sheet = 'Multimodal_Span')
df_SS <- read_excel("BD_Original.xlsx", sheet = 'Symmetry Span')
df_UT <- read_excel("BD_Original.xlsx", sheet = 'WMU Task')

## Gf (Fluid Intelligence) tasks
df_NS <- read_excel("BD_Original.xlsx", sheet = 'Number_Series')
df_RAPM <- read_excel("BD_Original.xlsx", sheet = 'RAPM')
df_LS <- read_excel("BD_Original.xlsx", sheet = 'Letter_Series')

1.2. Aggregating Sociodemographic Information and Performance Scores per Participant

#### WORKING MEMORY (WM) TASKS
# Select max Reading Span score per participant
df_RS_scores <- df_RS %>% 
  group_by(subject_nr) %>% 
  slice_max(score_reading_span, n = 1, with_ties = FALSE)
# Keep score column only
df_RS_scores <- df_RS_scores["score_reading_span"]
# Round scores to 2 decimals
df_RS_scores$score_reading_span <- round(df_RS_scores$score_reading_span, 2)

# Select max N-back score per participant
df_NB_scores <- df_NB %>% 
  group_by(subject_nr) %>% 
  slice_max(TotalCorrectTwoBack, n = 1, with_ties = FALSE)
# Convert raw correct count to proportion correct
df_NB_scores <- df_NB_scores["TotalCorrectTwoBack"] / 12
# Round scores to 2 decimals
df_NB_scores$TotalCorrectTwoBack <- round(df_NB_scores$TotalCorrectTwoBack, 2)

# Select max Operation Span score per participant
df_OS_scores <- df_OS %>% 
  group_by(subject_nr) %>% 
  slice_max(score_operation_span, n = 1, with_ties = FALSE)
# Keep score column only
df_OS_scores <- df_OS_scores["score_operation_span"]
# Round scores to 2 decimals
df_OS_scores$score_operation_span <- round(df_OS_scores$score_operation_span, 2)

# Select max Binding Task score per participant
df_BT_scores <- df_BT %>% 
  group_by(subject_nr) %>% 
  slice_max(BindingRawScore, n = 1, with_ties = FALSE)
# Convert to proportion correct
df_BT_scores <- df_BT_scores["BindingRawScore"] / 16
# Round scores to 2 decimals
df_BT_scores$BindingRawScore <- round(df_BT_scores$BindingRawScore, 2)

# Remove practice trials from Multimodal Span
df_MS_filtered <- df_MS %>% filter(practice != "yes")
# Select max Multimodal Span score per participant
df_MS_scores <- df_MS_filtered %>% 
  group_by(subject_nr) %>% 
  slice_max(MultimodalSpanScore, n = 1, with_ties = FALSE)
# Convert to proportion correct
df_MS_scores <- df_MS_scores["MultimodalSpanScore"] / 11
# Round scores to 2 decimals
df_MS_scores$MultimodalSpanScore <- round(df_MS_scores$MultimodalSpanScore, 2)

# Select max Symmetry Span score per participant
df_SS_scores <- df_SS %>% 
  group_by(subject_nr) %>% 
  slice_max(score_symmetry_span, n = 1, with_ties = FALSE)
# Keep score column only
df_SS_scores <- df_SS_scores["score_symmetry_span"]
# Round scores to 2 decimals
df_SS_scores$score_symmetry_span <- round(df_SS_scores$score_symmetry_span, 2)

# Select max WM Updating score per participant
df_UT_scores <- df_UT %>% 
  group_by(subject_nr) %>% 
  slice_max(WMUExperimentalScore, n = 1, with_ties = FALSE)
# Convert to proportion correct
df_UT_scores <- df_UT_scores["WMUExperimentalScore"] / 36
# Round scores to 2 decimals
df_UT_scores$WMUExperimentalScore <- round(df_UT_scores$WMUExperimentalScore, 2)

#### FLUID INTELLIGENCE (Gf) TASKS

# Select max Number Series score per participant
df_NS_scores <- df_NS %>% 
  group_by(subject_nr) %>% 
  slice_max(ExpermentalNumberSeriesScore, n = 1, with_ties = FALSE)
# Convert to proportion correct
df_NS_scores <- df_NS_scores["ExpermentalNumberSeriesScore"] / 15
# Round scores to 2 decimals
df_NS_scores$ExpermentalNumberSeriesScore <- round(df_NS_scores$ExpermentalNumberSeriesScore, 2)

# Select max RAPM score per participant
df_RAPM_scores <- df_RAPM %>% 
  group_by(subject_nr) %>% 
  slice_max(ExperimentalRAPMScore, n = 1, with_ties = FALSE)
# Convert to proportion correct
df_RAPM_scores <- df_RAPM_scores["ExperimentalRAPMScore"] / 18
# Round scores to 2 decimals
df_RAPM_scores$ExperimentalRAPMScore <- round(df_RAPM_scores$ExperimentalRAPMScore, 2)

# Select max Letter Series score per participant
df_LS_scores <- df_LS %>% 
  group_by(subject_nr) %>% 
  slice_max(ExperimentalLetterSeriesScore, n = 1, with_ties = FALSE)
# Convert to proportion correct
df_LS_scores <- df_LS_scores["ExperimentalLetterSeriesScore"] / 15
# Round scores to 2 decimals
df_LS_scores$ExperimentalLetterSeriesScore <- round(df_LS_scores$ExperimentalLetterSeriesScore, 2)

#### COMBINE SOCIODEMOGRAPHICS + WM SCORES

# Extract unique sociodemographic variables
df_sociodemographic <- distinct(df_RS[c(1:10)])
# Merge sociodemographics with WM task scores
df_WM_Tasks <- cbind(df_sociodemographic, df_RS_scores, df_NB_scores, df_OS_scores, df_BT_scores, df_MS_scores, df_SS_scores, df_UT_scores)
# Rename WM task columns
colnames(df_WM_Tasks)[c(11:17)] <- c("Reading_Span", "NBack", "Operation_Span", "Binding_Task", "Multimodal_Span", "Symmetry_Span", "WM_Updating_Task")
# Display WM dataframe interactively
datatable(df_WM_Tasks, 
          caption = "Df with sociodemographics and WM performance data", 
          options = list(
            scrollX = TRUE, 
            pageLength = 10, 
            searching = FALSE, 
            dom = 'tp'
            ), 
          class = 'display compact stripe hover'
          )
#### COMBINE SOCIODEMOGRAPHICS + Gf SCORES

# Merge sociodemographics with Gf task scores
df_Gf_Tasks <- cbind(df_sociodemographic, df_NS_scores, df_RAPM_scores, df_LS_scores)
# Rename Gf task columns
colnames(df_Gf_Tasks)[c(11:13)] <- c("Number_Series", "RAPM", "Letter_Series")
# Display Gf dataframe interactively
datatable( df_Gf_Tasks, 
           caption = "Df with sociodemographics and Gf performance data", 
           options = list( 
             scrollX = TRUE, 
             pageLength = 10, 
             searching = FALSE, 
             dom = 'tp' 
             ), 
           class = 'display compact stripe hover' 
           )

2. Applying Inclusion and Exclusion Criteria

Participation in this study was restricted to Portuguese citizens between 18 and 35 years of age who had completed at least a secondary education degree. Participants who obtained a score of zero in more than one cognitive task were also excluded from the analyses.

In the following subsections, we apply these inclusion and exclusion criteria to define the final sample. Specifically, we remove participants who did not meet nationality, age, or education requirements, and exclude those whose performance data suggest invalid or noncompliant responses (i.e., zero scores in multiple tasks). Before applying these criteria, the complete dataset contained the following number of participants:

print(paste("number of participants:", length(df_WM_Tasks$subject_nr)))
## [1] "number of participants: 169"

2.1. Filtering by Portuguese Citizenship

# Keep only Portuguese participants (various input spellings accounted for)
df_WM_Tasks <- df_WM_Tasks[df_WM_Tasks$nacionalidade %in% c("portuguesa", "Portuguesa", "portugues", "Portugues", 
                                                            "PORTUGUESA", "PORTUGUES", "pt", "PT", "Pt", 
                                                            "portuguese", "Portuguese", "Portugal", "portugal", 
                                                            "PORTUGAL", "prtuguesa", "Portugu\\êsa"), ]

# Apply same nationality filter to Gf dataset
df_Gf_Tasks <- df_Gf_Tasks[df_Gf_Tasks$nacionalidade %in% c("portuguesa", "Portuguesa", "portugues", "Portugues", 
                                                            "PORTUGUESA", "PORTUGUES", "pt", "PT", "Pt", 
                                                            "portuguese", "Portuguese", "Portugal", "portugal", 
                                                            "PORTUGAL", "prtuguesa", "Portugu\\êsa"), ]

# Print updated sample size
print(paste("number of participants:", length(df_WM_Tasks$subject_nr)))
## [1] "number of participants: 166"

2.2. Filtering by Age (18 – 35 Years)

# Keep participants aged between 18 and 35 in WM dataset
df_WM_Tasks <- df_WM_Tasks[df_WM_Tasks$idade >= 18 & df_WM_Tasks$idade <= 35, ]

# Apply same age filter to Gf dataset
df_Gf_Tasks <- df_Gf_Tasks[df_Gf_Tasks$idade >= 18 & df_Gf_Tasks$idade <= 35, ]

# Print updated sample size
print(paste("number of participants:", length(df_WM_Tasks$subject_nr)))
## [1] "number of participants: 164"

2.3. Filtering by Educational Level

# Exclude participants without completed secondary education in WM dataset
df_WM_Tasks <- df_WM_Tasks %>%
  filter(!(freqEnsSup == "no" & HabNaoEstudante == "Nenhuma das opções anteriores se aplica"))

# Apply same education filter to Gf dataset
df_Gf_Tasks <- df_Gf_Tasks %>%
  filter(!(freqEnsSup == "no" & HabNaoEstudante == "Nenhuma das opções anteriores se aplica"))

# Print updated sample size
print(paste("number of participants:", length(df_WM_Tasks$subject_nr)))
## [1] "number of participants: 164"

2.4. Excluding Participants with Zero Scores in More Than One Task

# Keep only ID and task score columns in WM dataset
df_WM_Tasks <- df_WM_Tasks[c(1,11:17)]

# Keep only ID and task score columns in Gf dataset
df_Gf_Tasks <- df_Gf_Tasks[c(1,11:13)]

# Combine WM and Gf datasets for zero-checking
df_concat_WM_Gf <- cbind(df_WM_Tasks, df_Gf_Tasks[-c(1)])

# Identify participants with zero scores in more than one task
subjects_with_multiple_zeros <- df_concat_WM_Gf %>%
  mutate(num_zeros = rowSums(across(-subject_nr, ~ . == 0))) %>%
  filter(num_zeros > 1) %>%
  distinct(subject_nr)

# Remove those participants from WM dataset
df_WM_Tasks <- df_WM_Tasks %>%
  filter(!subject_nr %in% subjects_with_multiple_zeros$subject_nr)

# Remove those participants from Gf dataset
df_Gf_Tasks <- df_Gf_Tasks %>%
  filter(!subject_nr %in% subjects_with_multiple_zeros$subject_nr)

# Print final sample size
print(paste("number of participants:", length(df_WM_Tasks$subject_nr)))
## [1] "number of participants: 163"

3. Identification and Treatment of Outliers and Missing Values

In this section we screened data for univariate extreme values, missing data, and multivariate anomalies. First, univariate outliers (defined as scores exceeding ±3 standard deviations from the mean) were identified and replaced with missing values. Next, missing data were imputed using multiple imputation to generate plausible replacements while preserving the covariance structure of the measures. Finally, multivariate outliers were detected using Mahalanobis distance and excluded from all imputed datasets.

3.1. Detection and Removal of Univariate Outliers (Mean ±3 SD Criterion)

# Convert zeros scores to missing values
df_WM_Tasks[df_WM_Tasks == 0] <- NA

# Replace values beyond ±3 SD of the mean by NA
df_WM_Tasks <- df_WM_Tasks %>%
  mutate(across(
    -subject_nr,
    ~ ifelse(. < (mean(., na.rm = TRUE) - 3 * sd(., na.rm = TRUE)) |
               . > (mean(., na.rm = TRUE) + 3 * sd(., na.rm = TRUE)),
             NA, .)
  ))


# Apply same procedure to Gf data
df_Gf_Tasks[df_Gf_Tasks == 0] <- NA
df_Gf_Tasks <- df_Gf_Tasks %>%
  mutate(across(
    -subject_nr,
    ~ ifelse(. < (mean(., na.rm = TRUE) - 3 * sd(., na.rm = TRUE)) |
               . > (mean(., na.rm = TRUE) + 3 * sd(., na.rm = TRUE)),
             NA, .)
  ))

#Displays table containning information regrading the sociodemographic information and the socres in the WM tasks of the participants 
datatable(
  df_WM_Tasks,
  caption = "df WM Tasks",
  options = list(
    scrollX = TRUE,       
    pageLength = 10,      
    searching = FALSE,    
    dom = 'tp'            
  ),
  class = 'display compact stripe hover'
)
#Displays table containning information regrading the sociodemographic information and the socres in the WM tasks of the participants
datatable(
  df_Gf_Tasks,
  caption = "df Gf Tasks",
  options = list(
    scrollX = TRUE,       
    pageLength = 10,      
    searching = FALSE,    
    dom = 'tp'            
  ),
  class = 'display compact stripe hover'
)

3.2. Replacement of Missing Values via Multiple Imputation

###Creates the predictor matrix that will be used to impute the missing values of each WM task. 
##The predictor matrix specifies which variables are used in the linear regression of each of the imputation models. 
##A value of 1 specifies that the variable given in the column name is used in the model to impute the variable given in the row name. 
##A value of  0 specifies that this variable is not used in this model.
imp_gen_df_WM_Tasks <- mice(df_WM_Tasks, print = FALSE)
pred1 <- imp_gen_df_WM_Tasks$predictorMatrix 
pred1[,"subject_nr"] <- 0
pred1["NBack","Multimodal_Span"] <- 0

##Creates the predictor matrix that will be used to impute the missing values of each Gf task. 
##The predictor matrix specifies which variables are used in the linear regression of each of the imputation models. 
##A value of 1 specifies that the variable given in the column name is used in the model to impute the variable given in the row name. 
##A value of  0 specifies that this variable is not used in this model.
imp_gen_df_Gf_Tasks <- mice(df_Gf_Tasks, print = FALSE)
pred2 <- imp_gen_df_Gf_Tasks$predictorMatrix 
pred2[,"subject_nr"] <- 0

###This section of the script uses multiple imputation to generate 20 databases in which the missing values of the WM tasks will be replaced by plausible values.
###The missing values replaced through multiple imputation will be slightly different each time you run this function. 
imp_gen_df_WM_Tasks <- mice(data=df_WM_Tasks,
                     method = 'pmm',
                     m=20,
                     maxit=10,
                     predictorMatrix = pred1,
                     printFlag = FALSE)

###This section of the script uses multiple imputation to generate 20 databases in which the missing values of the Gf tasks will be replaced by plausible values.
###The missing values replaced through multiple imputation will be slightly different each time you run this function. 
imp_gen_df_Gf_Tasks <- mice(data=df_Gf_Tasks,
                     method = 'pmm',
                     m=20,
                     maxit=10,
                     predictorMatrix = pred2,
                     printFlag = FALSE)


###Creates the variables mice.imp_DB and fills them with the DB that were generated through multiple imputation.
mice.imp_df_WM_Tasks <- NULL
mice.imp_df_Gf_Tasks <- NULL
mice.imp_df_WM_Gf_Tasks <- NULL
ParaMaha <- NULL
for(i in 1:20) mice.imp_df_WM_Tasks[[i]] <- complete(imp_gen_df_WM_Tasks, action=i, inc=FALSE)
subject_id <- mice.imp_df_WM_Tasks[[1]][,1]
for(i in 1:20) mice.imp_df_WM_Tasks[[i]] <- mice.imp_df_WM_Tasks[[i]][,-c(1)]
for(i in 1:20) ParaMaha[[i]] <- complete(imp_gen_df_WM_Tasks, action=i, inc=FALSE)
for(i in 1:20) ParaMaha[[i]] <- ParaMaha[[i]][,-c(1)]
for(i in 1:20) mice.imp_df_Gf_Tasks[[i]] <- complete(imp_gen_df_Gf_Tasks, action=i, inc=FALSE)
for(i in 1:20) mice.imp_df_Gf_Tasks[[i]] <- mice.imp_df_Gf_Tasks[[i]][,-c(1)]
for(i in 1:20) mice.imp_df_WM_Gf_Tasks[[i]] <- add_column(mice.imp_df_WM_Tasks[[i]],mice.imp_df_Gf_Tasks[[i]],.after = 8)
for(i in 1:20) mice.imp_df_WM_Gf_Tasks[[i]] <- data.frame(mice.imp_df_WM_Gf_Tasks[[i]])

3.3. Identification and Removal of Multivariate Outliers Using Mahalanobis Distance

###Calculates the Mahalonobis Distance and the respective p-value for each participant.
for(i in 1:20) mice.imp_df_WM_Gf_Tasks[[i]]$Mahalanobis_WMC_Gf <- mahalanobis(mice.imp_df_WM_Gf_Tasks[[i]],colMeans(mice.imp_df_WM_Gf_Tasks[[i]]),cov(mice.imp_df_WM_Gf_Tasks[[i]]))
### Compute corresponding p-values from the chi-square distribution (df = 9)
for(i in 1:20) mice.imp_df_WM_Gf_Tasks[[i]]$pvalue_WMC_Gf <- pchisq(mice.imp_df_WM_Gf_Tasks[[i]]$Mahalanobis, df=9, lower.tail=FALSE)
### Reinsert the subject ID column at the beginning of each dataset
for(i in 1:20) mice.imp_df_WM_Gf_Tasks[[i]] <- add_column(mice.imp_df_WM_Gf_Tasks[[i]],subject_id,.before = 1)

#############################################################################################################
# Extract only subject_id and pvalue_WMC_Gf columns from each imputed dataset
pvals_list <- lapply(mice.imp_df_WM_Gf_Tasks, function(df) df %>% select(subject_id, pvalue_WMC_Gf))

# Merge them all into one dataframe by subject_id
pvals_merged <- reduce(pvals_list, left_join, by = "subject_id")

# Compute the average p-value across the 20 datasets for each subject
pvals_merged <- pvals_merged %>%
  mutate(mean_pvalue = rowMeans(select(., starts_with("pvalue_WMC_Gf")), na.rm = TRUE))

# Identify subjects to exclude (average p < .05)
subjects_to_exclude <- pvals_merged %>%
  filter(mean_pvalue < 0.05) %>%
  pull(subject_id)

# Exclude those subjects from all 20 imputed datasets
for (i in 1:20) {
  mice.imp_df_WM_Gf_Tasks[[i]] <- mice.imp_df_WM_Gf_Tasks[[i]] %>%
    filter(!subject_id %in% subjects_to_exclude)
}

print(paste("number of participants:", nrow(mice.imp_df_WM_Gf_Tasks[[1]])))
## [1] "number of participants: 156"
#Example of an table containing imputed values of WM and Gf scores 
datatable(
  mice.imp_df_WM_Gf_Tasks[[1]],
  caption = "Example imputed df_WM_Gf_Tasks",
  options = list(
    scrollX = TRUE,       
    pageLength = 10,      
    searching = FALSE,    
    dom = 'tp'            
  ),
  class = 'display compact stripe hover'
)
###################################################################################################################################################################
##Generates a DB that includes averaged estimates across the 20 DBs.

### Create empty lists for storing pooled WM and Gf datasets
para_merge_WM_Tasks <- NULL
para_merge_Gf_Tasks <- NULL

### Extract subject IDs from the first imputed dataset
subject_id <- mice.imp_df_WM_Gf_Tasks[[1]][, c(1)]

### Select WM task columns (2–8) from each imputed dataset
for (i in 1:20) para_merge_WM_Tasks[[i]] <- mice.imp_df_WM_Gf_Tasks[[i]][, c(2:8)]

### Select Gf task columns (9–11) from each imputed dataset
for (i in 1:20) para_merge_Gf_Tasks[[i]] <- mice.imp_df_WM_Gf_Tasks[[i]][, c(9:11)]

### Keep only valid data frames in the WM task list
para_merge_WM_Tasks <- para_merge_WM_Tasks[vapply(para_merge_WM_Tasks, is.data.frame, logical(1))]

### Identify numeric columns common to all WM datasets
num_cols <- Reduce(
  intersect,
  lapply(para_merge_WM_Tasks, function(d) names(d)[vapply(d, is.numeric, logical(1))])
)

### Convert WM datasets to numeric matrices and combine into an array
arr <- simplify2array(lapply(para_merge_WM_Tasks, function(d) as.matrix(d[num_cols])))

### Compute mean scores across the 20 imputed WM datasets
para_merge_WM_Tasks <- as.data.frame(apply(arr, c(1, 2), mean, na.rm = TRUE))

### Keep only valid data frames in the Gf task list
para_merge_Gf_Tasks <- para_merge_Gf_Tasks[vapply(para_merge_Gf_Tasks, is.data.frame, logical(1))]

### Identify numeric columns common to all Gf datasets
num_cols <- Reduce(
  intersect,
  lapply(para_merge_Gf_Tasks, function(d) names(d)[vapply(d, is.numeric, logical(1))])
)

### Convert Gf datasets to numeric matrices and combine into an array
arr <- simplify2array(lapply(para_merge_Gf_Tasks, function(d) as.matrix(d[num_cols])))

### Compute mean scores across the 20 imputed Gf datasets
para_merge_Gf_Tasks <- as.data.frame(apply(arr, c(1, 2), mean, na.rm = TRUE))

### Combine participant IDs, WM, and Gf task means into one dataframe
df_combined_WM_Gf_tasks <- cbind(subject_id, para_merge_WM_Tasks, para_merge_Gf_Tasks)

### Round all numeric variables to two decimal places
df_combined_WM_Gf_tasks <- df_combined_WM_Gf_tasks %>%
  mutate(across(where(is.numeric), ~ round(., 2)))

4. Assessment of Univariate and Multivariate Normality

To verify that the data met the statistical assumptions required for the subsequent analyses, we examined univariate and multivariate normality across all 20 imputed datasets. Univariate normality was evaluated based on the criteria of |Skewness| < 2 and |Kurtosis| < 7. Multivariate normality was assessed using Mardia’s skewness and kurtosis coefficients, with values of Kurtosis < 3 and nonsignificant skewness (p > .05) indicating acceptable fit to the multivariate normal distribution. For each imputed dataset, Mardia’s test was computed separately, and pooled estimates of skewness and kurtosis were then obtained from the combined dataset.

The results of these analyses are reported in Table 1, presented in Section 6 (Descriptive Analysis).

###############################################################################################################################################################
paraNomr <- NULL
MultNormalVar <- NULL
###Generated 20 DBs that only include the scores of the WM and the Gf tasks.
for(i in 1:20) paraNomr[[i]] <- mice.imp_df_WM_Gf_Tasks[[i]][,-c(1,12,13)]
###Computes Mardia's Skewness coefficients and assess the skewness and kurtosis values of each of the 20 BDs. 
for(i in 1:20) MultNormalVar[[i]]<- mvn(paraNomr[[i]], mvnTest = "mardia")

##Generates two columns displaying pooled estimates of skewness and kurtosis for the distribution of each WM and Gf task
Megred_MultNormVar <- df_combined_WM_Gf_tasks[,-c(1)]
Megred_MultNormVar <- mvn(df_combined_WM_Gf_tasks,mvnTest = "mardia")
Megred_MultNormVar <- data.frame(Megred_MultNormVar$Descriptives)
Megred_MultNormVar <- Megred_MultNormVar[-c(1),c(9,10)]

5. Internal Consistency Analysis (Cronbach’s Alpha and McDonalds’ Omega)

This section evaluates the internal consistency of all WM and Gf tasks using Cronbach’s α and McDonald’s ω. Trial-level data for each task, stored in the file DB_alphas_omegas_tasks.xlsx, were coded as 1 (correct) or 0 (incorrect) to compute reliability coefficients at the item level.

Following conventional psychometric standards, internal consistency was considered acceptable when both α and ω values were equal to or greater than .70. For illustration, we display the Reading Span dataset as an example of the item-level structure used for these analyses.

The resulting α and ω coefficients for all tasks are summarized in Table 1, presented in Section 6 (Descriptive Analysis).

###Reading Span
original_RS <- read_excel("DB_alphas_omegas_tasks.xlsx",sheet = "Reading_Span")
original_RS <- original_RS %>%
  filter(!subject_nr %in% subjects_to_exclude)
original_RS <- original_RS[,-c(1)]
original_RS <- original_RS[-c(82),]
#Example of an table containing imputed values of WM and Gf scores 
datatable(
  original_RS,
  caption = "Df Correct/Incorrect by Item Reading Span",
  options = list(
    scrollX = TRUE,       
    pageLength = 10,      
    searching = FALSE,    
    dom = 'tp'            
  ),
  class = 'display compact stripe hover'
)
###Cronbach's Alpha
alpha_original_RS <- alpha(original_RS)
alpha_original_RS_CIs <- alpha.ci(original_RS,n.obs = nrow(df_combined_WM_Gf_tasks),n.var = 60,digits = 2)

###McDonald's W
omega_original_RS <- omega(original_RS,nfactors = 1,fm="ml")

#####################################################################################################################################################################
###Number Series
original_NSer <- read_excel("DB_alphas_omegas_tasks.xlsx",sheet = "Number_Series")
original_NSer <- original_NSer %>%
  filter(!subject_nr %in% subjects_to_exclude)
original_NSer <- original_NSer[,-c(1)]
original_NSer <- original_NSer[-c(82),]
original_NSer_sItem2 <- original_NSer[,-c(2)]

###Cronbach's Alpha
alpha_original_NSer_sItem2 <- alpha(original_NSer_sItem2)
alpha_original_NSer_sItem2_CIs <- alpha.ci(original_NSer_sItem2,n.obs = nrow(df_combined_WM_Gf_tasks),n.var = 15,digits = 2)

###McDonald's W
omega_original_sItem2 <- omega(original_NSer_sItem2,nfactors = 1,fm="ml")

#####################################################################################################################################################################
###NBack
original_NBack <- read_excel("DB_alphas_omegas_tasks.xlsx",sheet = "N_Back")
original_NBack <- original_NBack %>%
  filter(!subject_nr %in% subjects_to_exclude)
original_NBack <- original_NBack[,-c(1)]
original_NBack <- original_NBack[-c(82),]

###Cronbach's Alpha
alpha_original_NBack <- alpha(original_NBack)
alpha_original_NBack_CIs <- alpha.ci(original_NBack,n.obs = nrow(df_combined_WM_Gf_tasks),n.var = 12,digits = 2)

###McDonald's W
omega_original_NBack <- omega(original_NBack,nfactors = 1,fm="ml")

#####################################################################################################################################################################
###RAPM
original_RAPM <- read_excel("DB_alphas_omegas_tasks.xlsx",sheet = "RAPM")
original_RAPM <- original_RAPM %>%
  filter(!subject_nr %in% subjects_to_exclude)
original_RAPM <- original_RAPM[,-c(1)]
original_RAPM <- original_RAPM[-c(82),]

###Cronbach's Alpha
alpha_original_RAPM <- suppressMessages(
  suppressWarnings(
    psych::alpha(original_RAPM, check.keys = TRUE)
  )
)
alpha_original_RAPM_CIs <- alpha.ci(original_RAPM,n.obs = nrow(df_combined_WM_Gf_tasks),n.var = 18,digits = 2)

###McDonald's W
omega_original_RAPM <- omega(original_RAPM,nfactors = 1,fm="ml")

#####################################################################################################################################################################
###Operation Span
original_OS <- read_excel("DB_alphas_omegas_tasks.xlsx",sheet = "Operation_Span")
original_OS <- original_OS %>%
  filter(!subject_nr %in% subjects_to_exclude)
original_OS <- original_OS[,-c(1)]
original_OS <- original_OS[-c(82),]

###Cronbach's Alpha
alpha_original_OS <- alpha(original_OS)
alpha_original_OS_CIs <- alpha.ci(original_OS,n.obs = nrow(df_combined_WM_Gf_tasks),n.var = 60,digits = 2)

###McDonald's W
omega_original_OS <- omega(original_OS,nfactors = 1,fm="ml")

#####################################################################################################################################################################
###Binding Task
original_BT <- read_excel("DB_alphas_omegas_tasks.xlsx",sheet = "Binding_Task")
original_BT <- original_BT %>%
  filter(!subject_nr %in% subjects_to_exclude)
original_BT <- original_BT[,-c(1)]
original_BT <- original_BT[-c(82),]

###Cronbach's Alpha
alpha_original_BT <- alpha(original_BT)
alpha_original_BT_CIs <- alpha.ci(original_BT,n.obs = nrow(df_combined_WM_Gf_tasks),n.var = 16,digits = 2)

###McDonald's W
omega_original_BT <- omega(original_BT,nfactors = 1,fm="ml")

#####################################################################################################################################################################
###Multimodal Span
original_MS <- read_excel("DB_alphas_omegas_tasks.xlsx",sheet = "Multimodal_Span")
original_MS <- original_MS %>%
  filter(!subject_nr %in% subjects_to_exclude)
original_MS <- original_MS[,-c(1)]
original_MS <- original_MS[-c(82),]
original_MS <- original_MS[,-c(5)]

###Cronbach's Alpha
alpha_original_MS <- alpha(original_MS)
alpha_original_MS_CIs <- alpha.ci(original_MS,n.obs = nrow(df_combined_WM_Gf_tasks),n.var = 5,digits = 2)

###McDonald's W
omega_original_MS <- omega(original_MS,nfactors = 1,fm="ml")

#####################################################################################################################################################################
###Letter Series
original_LSer <- read_excel("DB_alphas_omegas_tasks.xlsx",sheet = "Letter_Series")
original_LSer <- original_LSer %>%
  filter(!subject_nr %in% subjects_to_exclude)
original_LSer <- original_LSer[,-c(1)]
original_LSer <- original_LSer[-c(82),]

###Cronbach's Alpha
alpha_original_LSer <- alpha(original_LSer)
alpha_original_LSer_CIs <- alpha.ci(original_LSer,n.obs = nrow(df_combined_WM_Gf_tasks),n.var = 15,digits = 2)

###McDonald's W
omega_original_LSer <- omega(original_LSer,nfactors = 1,fm="ml")

#####################################################################################################################################################################
###Symmetry Span
original_SS <- read_excel("DB_alphas_omegas_tasks.xlsx",sheet = "Symmetry Span")
original_SS <- original_SS %>%
  filter(!subject_nr %in% subjects_to_exclude)
original_SS <- original_SS[,-c(1)]
original_SS <- original_SS[-c(82),]
#here

###Cronbach's Alpha
alpha_original_SS <- alpha(original_SS)
alpha_original_SS_CIs <- alpha.ci(original_SS,n.obs = nrow(df_combined_WM_Gf_tasks),n.var = 60,digits = 2)

###McDonald's W
omega_original_SS <- omega(original_SS,nfactors = 1,fm="ml")

#####################################################################################################################################################################
###WM Updating Task
original_WMU <- read_excel("DB_alphas_omegas_tasks.xlsx",sheet = "WMU Task")
original_WMU <- original_WMU %>%
  filter(!subject_nr %in% subjects_to_exclude)
original_WMU <- original_WMU[,-c(1)]
original_WMU <- original_WMU[-c(82),]


###Cronbach's Alpha
alpha_original_WMU <- alpha(original_WMU)
alpha_original_WMU_CIs <- alpha.ci(original_WMU,n.obs = nrow(df_combined_WM_Gf_tasks),n.var = 12,digits = 2)

###McDonald's W
omega_original_WMU <- omega(original_WMU,nfactors = 1,fm="ml")

####################################################################################################################################################################
##Generates the columns with the values of Cronbach's alpha and McDonald's Omega that will be used appended to the table with descriptive statistics.

##Column containing the Cronbach's alpha for each task.
aRS <- alpha_original_RS$total$raw_alpha
aNB <- alpha_original_NBack$total$raw_alpha
aOS <- alpha_original_OS$total$raw_alpha
aBT <- alpha_original_BT$total$raw_alpha
aMS <- alpha_original_MS$total$raw_alpha
aSS <- alpha_original_SS$total$raw_alpha
aWMU <- alpha_original_WMU$total$raw_alpha
aNSer <- alpha_original_NSer_sItem2$total$raw_alpha
aRAPM <- alpha_original_RAPM$total$raw_alpha
aLSer <- alpha_original_LSer$total$raw_alpha

col_alphas <- data.frame(c(aRS,aNB,aOS,aBT,aMS,aSS,aWMU,aNSer,aRAPM,aLSer))

##Column containing McDonalds' Omega for each task.
wRS <- omega_original_RS$omega.tot
wNB <- omega_original_NBack$omega.tot
wOS <- omega_original_OS$omega.tot
wBT <- omega_original_BT$omega.tot
wMS <- omega_original_MS$omega.tot
wSS <- omega_original_SS$omega.tot
wWMU <- omega_original_WMU$omega.tot
wNSer <- omega_original_sItem2$omega.tot
wRAPM <- omega_original_RAPM$omega.tot
wLSer <- omega_original_LSer$omega.tot

col_omegas <- data.frame(c(wRS,wNB,wOS,wBT,wMS,wSS,wWMU,wNSer,wRAPM,wLSer))

6. Descriptive Statistics

This section presents the descriptive statistics of all working memory (WM) and fluid intelligence (Gf) tasks based on pooled estimates obtained from the 20 imputed datasets. For each measure, we computed the mean, standard deviation (SD), standard error (SE), variance, range, skewness, and kurtosis. In addition, Cronbach’s α and McDonald’s ω coefficients are reported to provide a comprehensive view of the reliability of each task.

The resulting summary of all variables is displayed in Table 1.

###Calculates the mean of the WM and Gf tasks and the respective covariance matrix.
descriptive_stats <- fmi(paraNomr)

### Calculates the mean and covariance matrix across imputed datasets
descriptive_stats <- fmi(paraNomr)

### Initialize empty dataframe for diagonal extraction
aux_calc_dev <- data.frame()

### Extract diagonal values (variances) from covariance matrix
for(i in 1:10) aux_calc_dev[i,1] <- descriptive_stats$Covariances$coef[i,i]

### Compute standard deviation (SD) for each task
for(i in 1:10) aux_calc_dev[i,2] <- sqrt(aux_calc_dev[i,1])

### Compute standard error (SE) for each task (n = 156)
for(i in 1:10) aux_calc_dev[i,3] <- aux_calc_dev[i,2]/sqrt(156)

### Create descriptive statistics dataframe (means + variability indices)
descriptive_stats <- descriptive_stats$Means$coef
descriptive_stats <- data.frame(descriptive_stats)
descriptive_stats <- add_column(descriptive_stats, aux_calc_dev[,2], .after = 1)
descriptive_stats <- add_column(descriptive_stats, aux_calc_dev[,3], .after = 2)
descriptive_stats <- add_column(descriptive_stats, aux_calc_dev[,1], .after = 3)

### Add task names for labeling
TaskNames <- c("Reading Span","N-Back Task","Operation Span","Binding Task",
               "Multimodal Span","Symmetry Span","WM Updating Task",
               "Number Series","RAPM","Letter Series")
descriptive_stats <- add_column(descriptive_stats, TaskNames, .before = 1)
colnames(descriptive_stats) <- c("Task_Names","Mean","SD","SE","Var")

### Compute maximum values for each task across all imputations
MaxMatrix <- matrix(data = NA, ncol = 10, nrow = 20)
for(i in 1:20) for(j in 1:10) MaxMatrix[i,j] <- max(paraNomr[[i]][,j])
MaxColumn <- NULL
for(i in 1:10) MaxColumn[i] <- max(MaxMatrix[,i])
MaxColumn <- data.frame(MaxColumn)
for(i in 1:10) MaxColumn[i,] <- round(MaxColumn[i,], digits = 2)
for(i in 1:10) MaxColumn[i,] <- toString(MaxColumn[i,])

### Compute minimum values for each task across all imputations
MinMatrix <- matrix(data = NA, ncol = 10, nrow = 20)
for(i in 1:20) for(j in 1:10) MinMatrix[i,j] <- min(paraNomr[[i]][,j])
MinColumn <- NULL
for(i in 1:10) MinColumn[i] <- min(MinMatrix[,i])
MinColumn <- data.frame(MinColumn)
for(i in 1:10) MinColumn[i,] <- round(MinColumn[i,], digits = 2)
for(i in 1:10) MinColumn[i,] <- toString(MinColumn[i,])

### Create a combined range column for each task
RangExpTasks <- MaxColumn
for(i in 1:10) RangExpTasks[i,] <- paste(MinColumn[i,], "-", MaxColumn[i,])

### Remove redundant columns (variance, SE) before merging
descriptive_stats <- descriptive_stats[,-c(4,5)]

### Add range, skewness, and kurtosis from Mardia’s analysis
descriptive_stats <- add_column(descriptive_stats, Megred_MultNormVar, .after = 3)
descriptive_stats <- add_column(descriptive_stats, RangExpTasks, .after = 3)

### Add Cronbach’s alpha and McDonald’s omega values
descriptive_stats <- add_column(descriptive_stats, col_alphas, .after = 6)
descriptive_stats <- add_column(descriptive_stats, col_omegas, .after = 7)

### Round numeric values to two decimals
descriptive_stats <- descriptive_stats %>%
  mutate(across(where(is.numeric), ~ round(., 2)))

### Rename and reorder columns for clarity
colnames(descriptive_stats) <- c("Task Names","Mean","SD","Range",
                                 "Skewness","Kurtosis","Alpha","Omega")
descriptive_stats <- descriptive_stats[c(1,3,6,2,7,5,4,10,8,9),]
rownames(descriptive_stats) <- NULL

### Reorder task names for final presentation
descriptive_stats[1] <- c("Reading Span","Operation Span","Symmetry Span",
                          "N-back Task","Updating Task","Multimodal Span",
                          "Binding Task","Letter Series","Number Series",
                          "Raven Advanced Progressive Matrices")

### Display final descriptive statistics table
knitr::kable(descriptive_stats, caption = "<b><span style='color:black;'>Table 1. Descriptive statistics (Mean, SD, Range, Skewness, Kurtosis, Alpha, Omega) for the Working Memory and Fluid Intelligence tasks</span></b>")
Table 1. Descriptive statistics (Mean, SD, Range, Skewness, Kurtosis, Alpha, Omega) for the Working Memory and Fluid Intelligence tasks
Task Names Mean SD Range Skewness Kurtosis Alpha Omega
Reading Span 0.77 0.16 0.27 - 1 -0.99 0.56 0.92 0.92
Operation Span 0.76 0.14 0.37 - 1 -0.45 -0.37 0.90 0.90
Symmetry Span 0.65 0.17 0.23 - 1 -0.28 -0.31 0.89 0.89
N-back Task 0.62 0.24 0.08 - 1 -0.21 -0.85 0.83 0.83
Updating Task 0.56 0.27 0.06 - 1 0.04 -1.28 0.93 0.93
Multimodal Span 0.38 0.07 0.27 - 0.55 0.45 0.00 0.51 0.55
Binding Task 0.81 0.15 0.31 - 1 -0.80 -0.13 0.74 0.74
Letter Series 0.37 0.13 0.07 - 0.73 0.13 -0.13 0.69 0.74
Number Series 0.41 0.11 0.2 - 0.8 0.63 -0.05 0.64 0.71
Raven Advanced Progressive Matrices 0.48 0.16 0.11 - 0.89 0.26 -0.41 0.63 0.64

Table 1 presents α and ω values for each task. Reliability was acceptable across all WM tasks except the multimodal span (α = .51; ω = .55). The low estimates likely reflect the unique structure of this task: the multimodal span ended when participants failed to reproduce two consecutive sequences of the same length. Scores were also the lowest among all WM tasks (.27–.55), indicating participants could replicate at most four sequences. Consequently, reliability was computed from only four trials, which may explain the reduced estimates. Considering this, the results concerning this task should be interpreted with caution.

Gf task reliability estimates also fell marginally below the conventional threshold for acceptable internal consistency (α and ω ≥ .70; Adadan & Savasci, 2012; McDonald, 1999). However, both these tasks are well-established Gf measures, typically showing high reliability (Buehner et al., 2005; Wiley et al., 2011). The relatively lower estimates observed here likely reflect the time limits imposed — participants had only five or ten minutes to complete as many items as possible. Comparable α and ω values have been reported in studies using similar constraints (Wilhelm et al., 2013).

7. Percentils

To facilitate the interpretation of individual differences in WM performance, this section presents percentile-based reference values for each WM task. Percentiles were computed from the pooled dataset using quartile, tercile, and median cutoffs (5th, 25th, 33rd, 50th, 66th, 75th, and 95th percentiles). These values provide an empirical reference for categorizing participant performance levels across tasks (e.g., lower, average, and higher performers).

The resulting table displays percentile scores for each WM measure based on the final pooled estimates.

### Compute percentiles (5th, 25th, 33rd, 50th, 66th, 75th, 95th) for all variables
percentils <- NULL
percentils <- data.frame(apply(df_combined_WM_Gf_tasks[], 2, quantile, 
                               probs = c(0.05, 0.25, 0.33, 0.50, 0.66, 0.75, 0.95)))

### Extract task names and add them as a separate column
Task_Names_1 <- rownames(percentils)
percentils <- add_columns(percentils, Task_Names_1, .before = 1)

### Remove redundant column
percentils <- percentils[,-c(2)]

### Reorder columns to desired task sequence
percentils <- percentils[, c(1,2,4,7,3,8,6,5,11,9,10)]
percentils <- data.frame(percentils)

### Keep only relevant columns for WM tasks
percentils <- percentils[, c(1,5,8,6,3,10,4,7)]

### Rename columns with task names and percentile label
TaskNames <- c("Percentile", "Reading Span", "Operation Span", "Symmetry Span", 
               "N-Back Task", "WM Updating Task", "Multimodal Span", "Binding Task")
colnames(percentils) <- TaskNames
percentils <- percentils %>%
  mutate(across(where(is.numeric), ~ round(., 2)))

### Display final percentile table
knitr::kable(percentils,  caption = "<b><span style='color:black;'>Table 2. Percentiles for the WM tasks</span></b>")
Table 2. Percentiles for the WM tasks
Percentile Reading Span Operation Span Symmetry Span N-Back Task WM Updating Task Multimodal Span Binding Task
5% 0.42 0.51 0.35 0.25 0.16 0.27 0.50
25% 0.68 0.67 0.57 0.42 0.33 0.36 0.75
33% 0.75 0.70 0.58 0.50 0.39 0.36 0.75
50% 0.78 0.78 0.67 0.67 0.53 0.36 0.88
66% 0.87 0.84 0.73 0.75 0.72 0.36 0.88
75% 0.88 0.88 0.75 0.83 0.81 0.45 0.94
95% 0.97 0.97 0.92 1.00 0.97 0.55 1.00

8. Correlation Matrix

This section computes the correlation matrix among the WM tasks included in the OpenWMB battery. The resulting correlation matrix summarizes the strength and direction of associations between all WM measures, serving as the empirical foundation for the subsequent factor analyses. Correlations were rounded to two decimal places for clarity.

### Initialize empty objects for correlation header and p-values
head1 <- NULL
pval1 <- NULL

### Compute pooled correlation and covariance estimates across 20 imputations
corr_WMC <- miceadds::micombine.cor(mi.res = paraNomr, variables = c(1:10))

### Extract p-values for the correlations
pval1 <- corr_WMC[, c(9)]

### Remove unneeded columns (retain only correlations between WM tasks)
corr_WMC <- corr_WMC[, -c(4:11)]

### Reshape correlation data into a wide matrix format
corr_WMC <- spread(corr_WMC, variable1, r)

### Fill diagonal with 1.00 (perfect self-correlation)
for(i in 1:10) corr_WMC[c(i), c(i + 1)] <- 1.00

### Reorder columns for logical presentation of tasks
corr_WMC <- corr_WMC[, c(1,9,7,10,5,11,4,2)]

### Reorder rows to match task order
corr_WMC <- corr_WMC[c(8,6,9,4,10,3,1),]

### Round correlation coefficients to two decimal places
corr_WMC <- corr_WMC %>%
  mutate(across(where(is.numeric), ~ round(., 2)))

### Remove row names for display consistency
rownames(corr_WMC) <- NULL

### Assign column names for tasks
colnames(corr_WMC) <- c("Task","Reading Span","Operation Span","Symmetry Span",
                        "N-Back Task","WM Updating Task","Multimodal Span","Binding Task")

### Assign task names to the first column
corr_WMC[, c(1)] <- c("Reading Span","Operation Span","Symmetry Span",
                      "N-Back Task","WM Updating Task","Multimodal Span","Binding Task")

### Display the final correlation matrix
knitr::kable(corr_WMC, caption = "<b><span style='color:black;'>Table 3. Correlation matrix of the WM tasks</span></b>")
Table 3. Correlation matrix of the WM tasks
Task Reading Span Operation Span Symmetry Span N-Back Task WM Updating Task Multimodal Span Binding Task
Reading Span 1.00 0.54 0.50 0.38 0.48 0.24 0.28
Operation Span 0.54 1.00 0.42 0.36 0.58 0.30 0.27
Symmetry Span 0.50 0.42 1.00 0.28 0.54 0.32 0.32
N-Back Task 0.38 0.36 0.28 1.00 0.40 0.14 0.36
WM Updating Task 0.48 0.58 0.54 0.40 1.00 0.39 0.31
Multimodal Span 0.24 0.30 0.32 0.14 0.39 1.00 0.31
Binding Task 0.28 0.27 0.32 0.36 0.31 0.31 1.00

All correlations were positive (most rs > .30) and significant, except between the n-back task and the multimodal span. These magnitudes were consistent with previous findings (Lewandowsky et al., 2010; Oswald et al., 2015; Schmiedek et al., 2014). The nonsignificant correlation between these two tasks was not considered problematic, as they tap distinct aspects of working memory (updating and binding, respectively).

9. Exploratory Factor Analysis

This section investigates the latent structure underlying performance across the WM tasks using Exploratory Factor Analysis (EFA). Prior to extraction, data suitability was assessed using the Kaiser–Meyer–Olkin (KMO) measure of sampling adequacy and Bartlett’s test of sphericity. Following standard recommendations, KMO values above .50 (Kaiser & Rice, 1974) and a significant Bartlett’s test (p < .05) were considered indicative of adequate factorability.

The number of factors to retain was determined based on Kaiser’s criterion (eigenvalues > 1), the scree plot, and parallel analysis. A single-factor solution was then extracted using the Maximum Likelihood (ML) method with Promax rotation, providing estimates of each task’s loading on the latent WM factor. Model adequacy was further examined through the proportion of explained variance and internal consistency of the factor.

##Creates a df with the merged imputed values of the WM tasks
merged_imp_gen_for_efa <- df_combined_WM_Gf_tasks[2:8]

##Calculates the Kaiser-Meyer-Olkin measure (KMO). 
##A KMO above .50 suggest that the characteristics of the data are suitable for factorial analysis (Kaiser & Rice, 1974). 
##The result of the KMO (.83) suggested that our data was adequate to conduct an EFA. 
KMO(merged_imp_gen_for_efa)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = merged_imp_gen_for_efa)
## Overall MSA =  0.83
## MSA for each item = 
##     Reading_Span            NBack   Operation_Span     Binding_Task 
##             0.84             0.82             0.83             0.82 
##  Multimodal_Span    Symmetry_Span WM_Updating_Task 
##             0.83             0.84             0.82
##Calculates Bartlett's test of sphericity. 
##A significant Bartlett's test of sphericity suggests that the characteristics of the data are suitable for factorial analysis.
##The results of Bartlett’s test of sphericity (χ2(21) = 301.53, p < .001) suggested that our data was adequate to conduct this analysis.
cortest.bartlett(merged_imp_gen_for_efa)
## $chisq
## [1] 304.7483
## 
## $p.value
## [1] 3.434143e-52
## 
## $df
## [1] 21
##Calculates eigenvalues
##We applied Kaiser’s criterion (1970) to decide how many factors to retain (eigenvalues > 1).
ev <- eigen(cor(merged_imp_gen_for_efa))
round(ev$values,2)
## [1] 3.29 0.89 0.86 0.62 0.52 0.49 0.35
##Computes a scree test to evaluate how many factors should be extracted.
scree(merged_imp_gen_for_efa, pc=FALSE)

##Computes a parallel analysis to evaluate how many factors should be extracted.
fa.parallel(merged_imp_gen_for_efa, fa="pc")

## Parallel analysis suggests that the number of factors =  NA  and the number of components =  1
##The results of the scree test, the parallel analysis, and Kaiser’s criterion (1970) suggested that a single factor was enough to accommodate all WM tasks.

##Assesses the fit of the unifactorial structure of WMC and estimates the loadings of each task in this factor.
##Method of extraction = Maximum Likelihood, Rotation = Promax. These extraction methods were selected because our data was normally distributed. 
fit <- factanal(merged_imp_gen_for_efa, factors = 1, rotation="promax")
### Display factor loadings (cutoff = .40) and model summary
print(fit, digits=3, cutoff=0.4, sort=TRUE)
## 
## Call:
## factanal(x = merged_imp_gen_for_efa, factors = 1, rotation = "promax")
## 
## Uniquenesses:
##     Reading_Span            NBack   Operation_Span     Binding_Task 
##            0.530            0.723            0.480            0.794 
##  Multimodal_Span    Symmetry_Span WM_Updating_Task 
##            0.799            0.561            0.393 
## 
## Loadings:
##                  Factor1
## Reading_Span     0.685  
## NBack            0.526  
## Operation_Span   0.721  
## Symmetry_Span    0.662  
## WM_Updating_Task 0.779  
## Binding_Task     0.454  
## Multimodal_Span  0.449  
## 
##                Factor1
## SS loadings      2.720
## Proportion Var   0.389
## 
## Test of the hypothesis that 1 factor is sufficient.
## The chi square statistic is 25.8 on 14 degrees of freedom.
## The p-value is 0.0274
### Extract and visualize factor loadings diagram
loads <- fit$loadings
fa.diagram(loads)     

##The unifactorial structure accounted for 38% of the variance in the WM tasks. 
##All tasks presented acceptable factor loadings (> .40).

##Computes Cronbach's alpha for the unifactorial structure of WMC.
alpha(merged_imp_gen_for_efa)
## 
## Reliability analysis   
## Call: alpha(x = merged_imp_gen_for_efa)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
##       0.79      0.81     0.8      0.37 4.2 0.023 0.65 0.12     0.37
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.73  0.79  0.83
## Duhachek  0.74  0.79  0.83
## 
##  Reliability if an item is dropped:
##                  raw_alpha std.alpha G6(smc) average_r S/N alpha se  var.r
## Reading_Span          0.75      0.77    0.76      0.36 3.4    0.027 0.0113
## NBack                 0.77      0.79    0.78      0.39 3.8    0.024 0.0128
## Operation_Span        0.75      0.77    0.75      0.36 3.3    0.027 0.0104
## Binding_Task          0.77      0.80    0.78      0.40 3.9    0.024 0.0148
## Multimodal_Span       0.79      0.81    0.79      0.41 4.1    0.024 0.0105
## Symmetry_Span         0.75      0.77    0.76      0.36 3.4    0.027 0.0131
## WM_Updating_Task      0.74      0.76    0.74      0.34 3.1    0.030 0.0098
##                  med.r
## Reading_Span      0.32
## NBack             0.32
## Operation_Span    0.32
## Binding_Task      0.40
## Multimodal_Span   0.40
## Symmetry_Span     0.37
## WM_Updating_Task  0.32
## 
##  Item statistics 
##                    n raw.r std.r r.cor r.drop mean    sd
## Reading_Span     156  0.72  0.73  0.67   0.61 0.77 0.161
## NBack            156  0.69  0.63  0.54   0.49 0.62 0.232
## Operation_Span   156  0.73  0.74  0.69   0.63 0.76 0.144
## Binding_Task     156  0.57  0.61  0.50   0.44 0.81 0.147
## Multimodal_Span  156  0.47  0.57  0.45   0.40 0.38 0.072
## Symmetry_Span    156  0.70  0.71  0.65   0.58 0.65 0.166
## WM_Updating_Task 156  0.83  0.78  0.76   0.67 0.56 0.273

The KMO value (.83) and Bartlett’s test of sphericity (χ²(21) = 305.26, p < .001) indicated that the data were suitable for factor analysis. The scree test, parallel analysis, and Kaiser’s criterion (1970) all supported a single-factor solution (eigenvalue = 3.29, evidencing convergent validity. This factor explained 39% of the variance in WM measures, a large effect size (ƒ² = .64). All tasks showed acceptable loadings (> .40; Field, 2017), and the factor was interpreted as an indirect estimate of WMC.

10. Confirmatory Factor Analysis

This section tests the adequacy of the unifactorial model of WM identified in the exploratory factor analysis (EFA). A confirmatory factor analysis (CFA) was performed to evaluate whether a single latent WMC factor could adequately represent the shared variance across the seven WM tasks.

The model was estimated using maximum likelihood (ML) on the 20 imputed datasets. Model fit was evaluated using several standard indices: the chi-square to degrees of freedom ratio (χ²/df; < 3 = acceptable), the Comparative Fit Index (CFI; > .90 = acceptable), the Root Mean Square Error of Approximation (RMSEA; < .08 = acceptable), and the Standardized Root Mean Square Residual (SRMR; < .08 = acceptable).

The analysis also generated the implied covariance and correlation matrices, standardized residuals, and a path diagram displaying standardized factor loadings, squared multiple correlations (R²), and residual variances for each observed indicator.

#Unifactorial model of WMC
model1WMCGeneralFactor <-
  'WMC =~ NA*Reading_Span + Operation_Span + Symmetry_Span + NBack + WM_Updating_Task + Multimodal_Span + Binding_Task
  WMC ~~ 1*WMC'

#The function bellow fits the unifactorial model to the data.
fit_model1WMCGeneralFactor <- runMI(model1WMCGeneralFactor, 
                                    data=mice.imp_df_WM_Gf_Tasks,
                                    fun="cfa")

##Computes the fit indexes of the model (e. g., chi-square, CFI, RMSEA,SRMR).
aaa <- summary(fit_model1WMCGeneralFactor,fit.measures = TRUE, standardized = TRUE,rsquare = TRUE)
## lavaan.mi object based on 20 imputed data sets. 
## See class?lavaan.mi help page for available methods. 
## 
## Convergence information:
## The model converged on 20 imputed data sets 
## 
## Rubin's (1987) rules were used to pool point and SE estimates across 20 imputed data sets, and to calculate degrees of freedom for each parameter's t test and CI.
## 
## Model Test User Model:
## 
##   Test statistic                                22.938
##   Degrees of freedom                                14
##   P-value                                        0.061
## 
## Model Test Baseline Model:
## 
##   Test statistic                               284.017
##   Degrees of freedom                                21
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.966
##   Tucker-Lewis Index (TLI)                       0.949
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.064
##   Confidence interval - lower                    0.000
##   Confidence interval - upper                    0.109
##   P-value H_0: RMSEA <= 0.05                     0.281
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.048
##########################
#Computes the implied pooled covariance matrix.
model_implied_1WMCGeneralFactor_cov <- fitted(fit_model1WMCGeneralFactor, omit.imps = c("no.conv", "no.se"))
model_implied_1WMCGeneralFactor_cov <- model_implied_1WMCGeneralFactor_cov$cov
if (is.list(model_implied_1WMCGeneralFactor_cov)) {
  model_implied_1WMCGeneralFactor_cov <- model_implied_1WMCGeneralFactor_cov[[1]]
}
model_implied_1WMCGeneralFactor_cov <- round(model_implied_1WMCGeneralFactor_cov, 2)
knitr::kable(model_implied_1WMCGeneralFactor_cov, caption = "<b><span style='color:black;'>Table 4. Implied Covariance Matrix of Model 1</span></b>")
Table 4. Implied Covariance Matrix of Model 1
Reading_Span Operation_Span Symmetry_Span NBack WM_Updating_Task Multimodal_Span Binding_Task
Reading_Span 0.03 0.01 0.01 0.01 0.02 0.00 0.01
Operation_Span 0.01 0.02 0.01 0.01 0.02 0.00 0.01
Symmetry_Span 0.01 0.01 0.03 0.01 0.02 0.00 0.01
NBack 0.01 0.01 0.01 0.06 0.03 0.00 0.01
WM_Updating_Task 0.02 0.02 0.02 0.03 0.07 0.01 0.01
Multimodal_Span 0.00 0.00 0.00 0.00 0.01 0.01 0.00
Binding_Task 0.01 0.01 0.01 0.01 0.01 0.00 0.02 
#Computes the implied pooled correlation matrix.
model_implied_1WMCGeneralFactor_corr <- data.frame(model_implied_1WMCGeneralFactor_cov)
colnames(model_implied_1WMCGeneralFactor_corr) <- c("Reading_Span","Operation_Span","Symmetry_span","NBack","WM_Updating_Task","Multimodal_Span","Binding_Task")
model_implied_1WMCGeneralFactor_corr <- data.matrix(model_implied_1WMCGeneralFactor_corr)
model_implied_1WMCGeneralFactor_corr <- cov2cor(model_implied_1WMCGeneralFactor_corr)
model_implied_1WMCGeneralFactor_corr <- round(model_implied_1WMCGeneralFactor_corr, 2)
knitr::kable(model_implied_1WMCGeneralFactor_corr, caption = "<b><span style='color:black;'>Table 5. Implied Correlation Matrix of Model 1</span></b>")
Table 5. Implied Correlation Matrix of Model 1
Reading_Span Operation_Span Symmetry_span NBack WM_Updating_Task Multimodal_Span Binding_Task
Reading_Span 1.00 0.41 0.33 0.24 0.44 0.00 0.41
Operation_Span 0.41 1.00 0.41 0.29 0.53 0.00 0.50
Symmetry_Span 0.33 0.41 1.00 0.24 0.44 0.00 0.41
NBack 0.24 0.29 0.24 1.00 0.46 0.00 0.29
WM_Updating_Task 0.44 0.53 0.44 0.46 1.00 0.38 0.27
Multimodal_Span 0.00 0.00 0.00 0.00 0.38 1.00 0.00
Binding_Task 0.41 0.50 0.41 0.29 0.27 0.00 1.00 
##Computes the standardized residual matrix of covariance. 
Std_residual_matrix_1WMCGenFactor <- cfa(model1WMCGeneralFactor,data = df_combined_WM_Gf_tasks, estimator = "ML")
Std_residual_matrix_1WMCGenFactor <- resid(Std_residual_matrix_1WMCGenFactor,type="standardized")
Std_residual_matrix_1WMCGenFactor <- Std_residual_matrix_1WMCGenFactor$cov
if (is.matrix(Std_residual_matrix_1WMCGenFactor)) {
  Std_residual_matrix_1WMCGenFactor <- Std_residual_matrix_1WMCGenFactor
} else if (is.list(Std_residual_matrix_1WMCGenFactor)) {
  Std_residual_matrix_1WMCGenFactor <- Std_residual_matrix_1WMCGenFactor[[1]]
} else {
  Std_residual_matrix_1WMCGenFactor <- as.matrix(Std_residual_matrix_1WMCGenFactor)
}
Std_residual_matrix_1WMCGenFactor <- round(Std_residual_matrix_1WMCGenFactor, 2)
knitr::kable(Std_residual_matrix_1WMCGenFactor, caption = "<b><span style='color:black;'>Table 6. Standardized Residual Matrix of Covariance of Model 1</span></b>")
Table 6. Standardized Residual Matrix of Covariance of Model 1
Reading_Span Operation_Span Symmetry_Span NBack WM_Updating_Task Multimodal_Span Binding_Task
Reading_Span 0.00 1.49 1.26 0.90 -2.48 -1.45 -0.54
Operation_Span 1.49 0.00 -1.90 -0.18 1.12 -0.49 -1.34
Symmetry_Span 1.26 -1.90 0.00 -1.45 0.82 0.53 0.42
NBack 0.90 -0.18 -1.45 0.00 0.05 -1.55 2.29
WM_Updating_Task -2.48 1.12 0.82 0.05 0.00 1.14 -1.09
Multimodal_Span -1.45 -0.49 0.53 -1.55 1.14 0.00 1.77
Binding_Task -0.54 -1.34 0.42 2.29 -1.09 1.77 0.00 
##########################
#Generates a plot displays the standardized factor loadings, squared multiple correlations, and standardized error terms of model 1.
fit_model_test1 <- cfa(model1WMCGeneralFactor, 
                       data=mice.imp_df_WM_Gf_Tasks[[1]],estimator="ML")
SEMPLOT <- semPlot::semPlotModel(fit_model_test1)
aaaa <- data.frame(summary(fit_model1WMCGeneralFactor,fit.measures = TRUE, standardized = TRUE,rsquare = TRUE))
aaaa <- aaaa$std.all[-c(16:22)]
SEMPLOT@Pars$std <- aaaa
Rsq_1WMC <- summary(fit_model1WMCGeneralFactor,fit.measures = TRUE, standardized = TRUE,rsquare = TRUE)
Rsq_1WMC <- data.frame(Rsq_1WMC)
for(i in 1:length(Rsq_1WMC[,"std.all"])) Rsq_1WMC[i,"RSq"] <- as.character(round(Rsq_1WMC[i,"std.all"]^2,2))
vecNames1WMC <- c("RS","OS","SS","NB","UT","MS","BT") 
vecRsq1WMC <- Rsq_1WMC[1:7,"RSq"]
vecIndi1WMC <- paste(vecNames1WMC," (",vecRsq1WMC,")")
vecIndi1WMC[8] <- "WMC"
vecIndi1WMC[5] <- "UT ( 0.60 )"
vecIndi1WMC <- c("RS (.47)","OS (.48)","SS (.46)","NB (.23)","UT (.60)","MS (.18)","BT (.18)", "WMC")
Path1WMC <- semPaths(SEMPLOT,'std',intercepts = TRUE, residuals = TRUE, layout = "tree", edge.color = "black",rotation = 4,sizeMan = 8,sizeLat = 10,nodeLabels = vecIndi1WMC, edge.width=0.7, fade = FALSE, fixedStyle = c("black"),freeStyle = c("black"),cut=1, ,edge.label.cex = 0.75,mar=c(4,4,4,4))

Model 1 showed a good fit to the data, χ²(14) = 23.26, p = .06; χ²/df = 1.67; CFI = 0.97; RMSEA = 0.065; SRMR = 0.048. All factor loadings were significant (p < .001) and acceptable, indicating strong factor reliability. Standardized loadings, squared multiple correlations, and error terms are shown in the path diagram. Additionally, Model 1 presented good reliability (ω = .79).

11. Structural Equation Model

This section tests the relationship between the unifactorial WM model and the Fluid Intelligence (Gf) factor using a Structural Equation Model (SEM). The model specifies two latent variables—one for WMC, defined by seven working memory tasks, and one for Gf, defined by three reasoning tasks.

Model estimation was performed across the 20 imputed datasets using Maximum Likelihood (ML) estimation within a multiple imputation framework. Model adequacy was evaluated using standard fit indices (χ², CFI, RMSEA, and SRMR). We also inspected the implied covariance and correlation matrices, and standardized residuals to assess local fit and interpret the structural relationships. A final path diagram summarizes standardized loadings, R² values, and error variances for each indicator.

##Model that was used to assess the correlation between the unifactorial model of WMC and the Gf factor.
model_1GfFactor_1WMCFactor <- 'WMC =~ NA*Reading_Span + Operation_Span + Symmetry_Span + NBack + WM_Updating_Task + Multimodal_Span + Binding_Task
              Gf =~ NA*RAPM + Letter_Series + Number_Series
              Gf ~ WMC
              WMC ~~ 1*WMC
              Gf ~~ 1*Gf'

#The function bellow fits the model to the data.
fit_model_1GfFactor_1WMCFactor <- runMI(model_1GfFactor_1WMCFactor, 
                                        data=mice.imp_df_WM_Gf_Tasks,
                                        fun="sem")

#Computes the fit indexes of the model (e. g., chi-square, CFI, RMSEA,SRMR).
aaa <- summary(fit_model_1GfFactor_1WMCFactor,fit.measures = TRUE, standardized = TRUE,rsquare = TRUE,test="D2")
## lavaan.mi object based on 20 imputed data sets. 
## See class?lavaan.mi help page for available methods. 
## 
## Convergence information:
## The model converged on 20 imputed data sets 
## 
## Rubin's (1987) rules were used to pool point and SE estimates across 20 imputed data sets, and to calculate degrees of freedom for each parameter's t test and CI.
## 
## Model Test User Model:
## 
##   Test statistic                                75.412
##   Degrees of freedom                                34
##   P-value                                        0.000
## 
## Model Test Baseline Model:
## 
##   Test statistic                               508.406
##   Degrees of freedom                                45
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.911
##   Tucker-Lewis Index (TLI)                       0.882
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.088
##   Confidence interval - lower                    0.061
##   Confidence interval - upper                    0.115
##   P-value H_0: RMSEA <= 0.05                     0.012
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.058
##########################
##Computes the implied pooled covariance matrix.
model_implied_1GfFactor_1WMCFactor_cov <- fitted(fit_model_1GfFactor_1WMCFactor, omit.imps = c("no.conv", "no.se"))
if (is.list(model_implied_1GfFactor_1WMCFactor_cov)) {
  model_implied_1GfFactor_1WMCFactor_cov <- model_implied_1GfFactor_1WMCFactor_cov[[1]]
}
model_implied_1GfFactor_1WMCFactor_cov <- round(model_implied_1GfFactor_1WMCFactor_cov, 2)
knitr::kable(model_implied_1GfFactor_1WMCFactor_cov, caption = "<b><span style='color:black;'>Table 7. Implied Covariance Matrix of Model 2</span></b>")
Table 7. Implied Covariance Matrix of Model 2
Reading_Span Operation_Span Symmetry_Span NBack WM_Updating_Task Multimodal_Span Binding_Task RAPM Letter_Series Number_Series
Reading_Span 0.03 0.01 0.01 0.01 0.02 0.00 0.01 0.01 0.01 0.01
Operation_Span 0.01 0.02 0.01 0.01 0.02 0.00 0.01 0.01 0.01 0.01
Symmetry_Span 0.01 0.01 0.03 0.01 0.02 0.00 0.01 0.01 0.01 0.01
NBack 0.01 0.01 0.01 0.06 0.03 0.00 0.01 0.01 0.01 0.01
WM_Updating_Task 0.02 0.02 0.02 0.03 0.07 0.01 0.02 0.02 0.02 0.01
Multimodal_Span 0.00 0.00 0.00 0.00 0.01 0.01 0.00 0.00 0.00 0.00
Binding_Task 0.01 0.01 0.01 0.01 0.02 0.00 0.02 0.01 0.01 0.00
RAPM 0.01 0.01 0.01 0.01 0.02 0.00 0.01 0.02 0.01 0.01
Letter_Series 0.01 0.01 0.01 0.01 0.02 0.00 0.01 0.01 0.02 0.01
Number_Series 0.01 0.01 0.01 0.01 0.01 0.00 0.00 0.01 0.01 0.01 
#Computes the implied pooled correlation matrix.
model_implied_1GfFactor_1WMCFactor_corr <- data.frame(model_implied_1GfFactor_1WMCFactor_cov)
colnames(model_implied_1GfFactor_1WMCFactor_corr) <- c("Reading_Span","Operation_Span","Symmetry_span","NBack","WM_Updating_Task","Multimodal_Span","Binding_Task","RAPM","Letter_Series","Number_Series")
model_implied_1GfFactor_1WMCFactor_corr <- data.matrix(model_implied_1GfFactor_1WMCFactor_corr)
model_implied_1GfFactor_1WMCFactor_corr <- cov2cor(model_implied_1GfFactor_1WMCFactor_corr)
model_implied_1GfFactor_1WMCFactor_corr <- round(model_implied_1GfFactor_1WMCFactor_corr, 2)
knitr::kable(model_implied_1GfFactor_1WMCFactor_corr, caption = "<b><span style='color:black;'>Table 8. Implied Correlation Matrix of Model 2</span></b>")
Table 8. Implied Correlation Matrix of Model 2
Reading_Span Operation_Span Symmetry_span NBack WM_Updating_Task Multimodal_Span Binding_Task RAPM Letter_Series Number_Series
Reading_Span 1.00 0.41 0.33 0.24 0.44 0.00 0.41 0.41 0.41 0.58
Operation_Span 0.41 1.00 0.41 0.29 0.53 0.00 0.50 0.50 0.50 0.71
Symmetry_Span 0.33 0.41 1.00 0.24 0.44 0.00 0.41 0.41 0.41 0.58
NBack 0.24 0.29 0.24 1.00 0.46 0.00 0.29 0.29 0.29 0.41
WM_Updating_Task 0.44 0.53 0.44 0.46 1.00 0.38 0.53 0.53 0.53 0.38
Multimodal_Span 0.00 0.00 0.00 0.00 0.38 1.00 0.00 0.00 0.00 0.00
Binding_Task 0.41 0.50 0.41 0.29 0.53 0.00 1.00 0.50 0.50 0.00
RAPM 0.41 0.50 0.41 0.29 0.53 0.00 0.50 1.00 0.50 0.71
Letter_Series 0.41 0.50 0.41 0.29 0.53 0.00 0.50 0.50 1.00 0.71
Number_Series 0.58 0.71 0.58 0.41 0.38 0.00 0.00 0.71 0.71 1.00 
###Computes the standardized residual matrix of covariance.
Std_residual_matrix_1GfFactor_1WMCFactor <- cfa(model_1GfFactor_1WMCFactor,data = df_combined_WM_Gf_tasks, estimator = "ML")
Std_residual_matrix_1GfFactor_1WMCFactor <- resid(Std_residual_matrix_1GfFactor_1WMCFactor,type="standardized")
Std_residual_matrix_1GfFactor_1WMCFactor <- Std_residual_matrix_1GfFactor_1WMCFactor$cov
if (is.matrix(Std_residual_matrix_1GfFactor_1WMCFactor)) {
  Std_residual_matrix_1GfFactor_1WMCFactor <- Std_residual_matrix_1GfFactor_1WMCFactor
} else if (is.list(Std_residual_matrix_1GfFactor_1WMCFactor)) {
  Std_residual_matrix_1GfFactor_1WMCFactor <- Std_residual_matrix_1GfFactor_1WMCFactor[[1]]
} else {
  Std_residual_matrix_1GfFactor_1WMCFactor <- as.matrix(Std_residual_matrix_1GfFactor_1WMCFactor)
}
Std_residual_matrix_1GfFactor_1WMCFactor <- round(Std_residual_matrix_1GfFactor_1WMCFactor, 2)
knitr::kable(Std_residual_matrix_1GfFactor_1WMCFactor, caption = "<b><span style='color:black;'>Table 9. Standardized Residual Matrix of Covariance of Model 2</span></b>")
Table 9. Standardized Residual Matrix of Covariance of Model 2
Reading_Span Operation_Span Symmetry_Span NBack WM_Updating_Task Multimodal_Span Binding_Task RAPM Letter_Series Number_Series
Reading_Span 0.00 2.04 1.90 1.17 -1.84 -1.37 -0.56 -1.71 -1.52 -0.08
Operation_Span 2.04 0.00 -0.93 -0.05 0.72 -0.74 -1.62 0.40 -2.16 1.29
Symmetry_Span 1.90 -0.93 0.00 -1.09 0.86 0.42 0.27 1.02 -2.28 -1.35
NBack 1.17 -0.05 -1.09 0.00 -0.38 -1.78 2.09 1.25 -0.24 -0.55
WM_Updating_Task -1.84 0.72 0.86 -0.38 0.00 0.39 -2.21 -0.07 -0.95 3.53
Multimodal_Span -1.37 -0.74 0.42 -1.78 0.39 0.00 1.46 1.39 1.72 -0.83
Binding_Task -0.56 -1.62 0.27 2.09 -2.21 1.46 0.00 0.97 1.94 0.79
RAPM -1.71 0.40 1.02 1.25 -0.07 1.39 0.97 0.00 2.11 -3.73
Letter_Series -1.52 -2.16 -2.28 -0.24 -0.95 1.72 1.94 2.11 0.00 0.72
Number_Series -0.08 1.29 -1.35 -0.55 3.53 -0.83 0.79 -3.73 0.72 0.00 
##########################
#Generates a plot displays the standardized factor loadings, squared multiple correlations, and standardized error terms of model 2.
fit_model_test3 <- sem(model_1GfFactor_1WMCFactor, 
                       data=mice.imp_df_WM_Gf_Tasks[[1]],estimator="ML")
SEMPLOT <- semPlot::semPlotModel(fit_model_test3)
aaaa<- data.frame(summary(fit_model_1GfFactor_1WMCFactor,fit.measures = TRUE, standardized = TRUE,rsquare = TRUE))
aaaa <- aaaa$std.all[-c(24:34)]
aaaa2 <- aaaa
aaaa2[c(1:7)] <- rev(aaaa2[c(1:7)])
aaaa2[c(9:10)] <- rev(aaaa2[c(9:10)])
aaaa2[c(14:20)] <- rev(aaaa2[c(14:20)])
aaaa2[c(22:23)] <- rev(aaaa2[c(22:23)])
SEMPLOT@Pars$std <- aaaa2
Rsq_1WMC_1Gf <- summary(fit_model_1GfFactor_1WMCFactor,fit.measures = TRUE, standardized = TRUE,rsquare = TRUE)
Rsq_1WMC_1Gf <- data.frame(Rsq_1WMC_1Gf)
for(i in 1:length(Rsq_1WMC_1Gf[,"std.all"])) Rsq_1WMC_1Gf[i,"RSq"] <- as.character(round(Rsq_1WMC_1Gf[i,"std.all"]^2,2))
vecNames1WMC_1Gf <- c("RS","OS","SS","NB","UT","MS","BT","RAPM","L_Ser","N_Ser","Gf") 
vecRsq_1WMC_1Gf <- Rsq_1WMC_1Gf[1:11,"RSq"]
vecIndi1WMC_1Gf <- paste(vecNames1WMC_1Gf," (",vecRsq_1WMC_1Gf,")")
vecIndi1WMC_1Gf[12] <- "WMC"
vecIndi1WMC_1Gf[6] <- "MS ( 0.20 )"
vecIndi1WMC_1Gf <- c("BT (.21)","MS (.20)","UT (.65)","NB (.23)","SS (.43)","OS (.46)","RS (.42)","RAPM (.37)","NSer (.45)","LSer (.49)","WMC","Gf")
SEMpath <- semPaths(SEMPLOT,'std',intercepts = TRUE, residuals = TRUE, layout = "tree", edge.color = "black",rotation = 2,sizeMan = 9,sizeLat = 11, optimizeLatRes = TRUE,
         nodeLabels = vecIndi1WMC_1Gf, edge.width=0.7, fade = FALSE, fixedStyle = c("black"),freeStyle = c("black"),cut=1, edge.label.cex = 0.75, mar = c(5,5,5,5))

Model 2 showed an acceptable fit to the data, (χ²(34) = 74.94, p = .001), χ²/df = 2.20, CFI = 0.91, RMSEA = 0.088, and SRMR = 0.058. All factor loadings were significant (p < .001), indicating good factor reliability. The WMC and Gf factors were strongly correlated (r = .81, p < .001), consistent with previous findings (Hicks et al., 2016; Schmiedek et al., 2009, 2014; Wilhelm et al., 2013). Standardized loadings, squared multiple correlations, and error terms are shown in the path diagram. The Gf factor also demonstrated adequate reliability (ω = .70).

12. Discussion

The main goal of this analysis was to evaluate the psychometric properties of the OpenWMB, an open-source, automated battery of heterogeneous WM tasks. Overall, the battery and most individual WM tasks demonstrated good internal consistency (most indexes > .70).

The OpenWMB also showed adequate convergent validity, evidenced by moderate to strong positive correlations among WM tasks (most rs > .30) and by the results of the EFA and CFA.

All tasks loaded onto a single latent factor representing a common source of variance — likely individual differences in WM. This factor accounted for 38% of the variance in WM performance, corresponding to a large effect size (ƒ² = .64). Moreover, the general WM factor correlated strongly with a latent Gf factor (r = .81, p < .001), replicating the well-established association between WMC and Gf and supporting the battery’s predictive validity.

In sum, the OpenWMB demonstrated strong psychometric qualities, including internal consistency, convergent validity, and predictive validity. Thus, the battery provides a reliable and valid estimate of general WM by capturing several of its functional components (e.g., simultaneous storage and processing, updating, binding), while minimizing paradigm- and task-specific variance.