Case Study 1 — Analytical strategy based on reduced EEG features (R)¶
This is the R implementation of Case Study 1. It implements the analysis described in Section 4.1 of the companion paper. The same analysis is also available as Python code.
The research question is:
Are the amplitude and latency of specific ERP components associated with neuropsychological performance and demographic characteristics in post-stroke patients?
Two well-characterised ERP components are used:
- N400 — a negative deflection ~400 ms post-stimulus, indexing semantic processing (Kutas & Federmeier, 2011)
- P300 — a positive deflection ~300 ms post-target, indexing attentional resource allocation (Polich, 2007)
From each component we consider amplitude (mean µV within a predefined time window) and latency (ms of peak amplitude), yielding four participant-level neural predictors.
As it usually happens in real life, data are delivered in two different files: dataset_NPSY.csv and dataset_ERP.csv.
Pipeline¶
- Data Wrangling — load, inspect, check ranges, missing data, distributions, outliers
- Analysis — descriptive stats, bivariate correlations, multiple regression, post-hoc diagnostics
- Export — save tables and figures
All code is defined directly here for transparency and reproducibility. Each step is explained as we go.
0. Setup¶
Import the libraries we will use throughout the analysis.
library(dplyr)
library(tidyr)
library(ggplot2)
library(patchwork) # combine ggplot panels into one figure
library(car) # variance inflation factor
library(psych) # skewness / kurtosis, common in psych/neuropsych stats
library(openxlsx) # Excel export
theme_set(theme_minimal())
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
Loading required package: carData
Attaching package: ‘car’
The following object is masked from ‘package:dplyr’:
recode
Attaching package: ‘psych’
The following object is masked from ‘package:car’:
logit
The following objects are masked from ‘package:ggplot2’:
%+%, alpha
Define input/output folders
# Input
DATA_DIR <- "../../../data"
NPSY_CSV <- file.path(DATA_DIR, "dataset_NPSY.csv")
EEG_CSV <- file.path(DATA_DIR, "dataset_ERP.csv")
# Output
OUTPUTS <- "../outputs"
FIGS <- file.path(OUTPUTS, "figures")
RESULTS <- file.path(OUTPUTS, "results")
dir.create(FIGS, recursive = TRUE, showWarnings = FALSE)
dir.create(RESULTS, recursive = TRUE, showWarnings = FALSE)
1. Data Wrangling¶
1.1 Load data¶
The datasets are stored as separate CSV files linked by subject_id.
npsy <- read.csv(NPSY_CSV)
eeg <- read.csv(EEG_CSV)
cat(sprintf("NPSY: %d x %d columns: %s\n", nrow(npsy), ncol(npsy), paste(colnames(npsy), collapse = ", ")))
cat(sprintf("EEG: %d x %d columns: %s\n", nrow(eeg), ncol(eeg), paste(colnames(eeg), collapse = ", ")))
NPSY: 45 x 7 columns: subject_id, age, sex, criq, apacs, att_matr, moca
EEG: 45 x 5 columns: subject_id, n400_amp, n400_lat, p300_amp, p300_lat
Let's inspect the datasets looking at the first rows: each row stores the NPSY/EEG data of a single subject.
cat("NPSY\n"); print(head(npsy))
cat("EEG\n"); print(head(eeg))
NPSY
subject_id age sex criq apacs att_matr moca 1 sub-001 74 0 78 0.54 25 24 2 sub-002 62 1 95 0.67 38 21 3 sub-003 55 1 94 0.93 42 30 4 sub-004 68 0 97 0.40 28 13 5 sub-005 58 1 103 0.72 27 30 6 sub-006 62 0 107 0.86 39 8
EEG
subject_id n400_amp n400_lat p300_amp p300_lat 1 sub-001 3.111 458.7 14.671 515.5 2 sub-002 1.794 415.3 12.818 408.0 3 sub-003 4.587 403.6 13.659 314.6 4 sub-004 3.168 457.2 12.911 454.2 5 sub-005 3.604 434.4 8.635 353.4 6 sub-006 3.785 405.4 12.412 332.0
To create a single dataset we need to merge on subject_id the data collected for each subject
df <- inner_join(npsy, eeg, by = "subject_id")
cat(sprintf("Merged: %d subjects x %d variables\n", nrow(df), ncol(df)))
head(df)
Merged: 45 subjects x 11 variables
| subject_id | age | sex | criq | apacs | att_matr | moca | n400_amp | n400_lat | p300_amp | p300_lat | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <dbl> | <int> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | sub-001 | 74 | 0 | 78 | 0.54 | 25 | 24 | 3.111 | 458.7 | 14.671 | 515.5 |
| 2 | sub-002 | 62 | 1 | 95 | 0.67 | 38 | 21 | 1.794 | 415.3 | 12.818 | 408.0 |
| 3 | sub-003 | 55 | 1 | 94 | 0.93 | 42 | 30 | 4.587 | 403.6 | 13.659 | 314.6 |
| 4 | sub-004 | 68 | 0 | 97 | 0.40 | 28 | 13 | 3.168 | 457.2 | 12.911 | 454.2 |
| 5 | sub-005 | 58 | 1 | 103 | 0.72 | 27 | 30 | 3.604 | 434.4 | 8.635 | 353.4 |
| 6 | sub-006 | 62 | 0 | 107 | 0.86 | 39 | 8 | 3.785 | 405.4 | 12.412 | 332.0 |
1.2 Initial inspection¶
We check dataset dimensions, column types, and basic descriptive statistics.
Expected clinical ranges:
| Variable | Expected range |
|---|---|
| Age | 18–100 years |
| CRIq | 70–130 |
| N400 amplitude | 1–6 µV |
| N400 latency | 350–500 ms |
| P300 amplitude | 5–20 µV |
| P300 latency | 300–500 ms |
| APACS | 0–1 |
| Attentive Matrices | 0–60 |
| MoCA | 0–30 |
cat(sprintf("Shape: %d rows x %d columns\n\n", nrow(df), ncol(df)))
cat("First 5 rows:\n"); print(head(df))
cat("\nDescriptive statistics:\n")
summary(df)
Shape: 45 rows x 11 columns
First 5 rows:
subject_id age sex criq apacs att_matr moca n400_amp n400_lat p300_amp 1 sub-001 74 0 78 0.54 25 24 3.111 458.7 14.671 2 sub-002 62 1 95 0.67 38 21 1.794 415.3 12.818 3 sub-003 55 1 94 0.93 42 30 4.587 403.6 13.659 4 sub-004 68 0 97 0.40 28 13 3.168 457.2 12.911 5 sub-005 58 1 103 0.72 27 30 3.604 434.4 8.635 6 sub-006 62 0 107 0.86 39 8 3.785 405.4 12.412 p300_lat 1 515.5 2 408.0 3 314.6 4 454.2 5 353.4 6 332.0
Descriptive statistics:
subject_id age sex criq
Length:45 Min. : 49.00 Min. :0.0000 Min. : 74.00
Class :character 1st Qu.: 60.00 1st Qu.:0.0000 1st Qu.: 89.00
Mode :character Median : 65.00 Median :0.0000 Median : 95.00
Mean : 66.91 Mean :0.4889 Mean : 95.82
3rd Qu.: 70.00 3rd Qu.:1.0000 3rd Qu.:104.50
Max. :150.00 Max. :1.0000 Max. :119.00
NA's :1
apacs att_matr moca n400_amp
Min. :0.3900 Min. : 3.00 Min. : 8.00 Min. : 1.794
1st Qu.:0.6000 1st Qu.:25.00 1st Qu.:18.00 1st Qu.: 2.827
Median :0.6900 Median :31.00 Median :22.00 Median : 3.530
Mean :0.6922 Mean :31.53 Mean :21.77 Mean : 3.985
3rd Qu.:0.7600 3rd Qu.:38.00 3rd Qu.:26.00 3rd Qu.: 4.048
Max. :0.9300 Max. :48.00 Max. :30.00 Max. :25.000
NA's :1 NA's :1
n400_lat p300_amp p300_lat
Min. :397.0 Min. : 8.635 Min. :292.9
1st Qu.:412.4 1st Qu.:11.188 1st Qu.:359.1
Median :427.2 Median :12.206 Median :400.1
Mean :427.3 Mean :12.218 Mean :412.4
3rd Qu.:442.3 3rd Qu.:13.068 3rd Qu.:429.3
Max. :460.1 Max. :17.031 Max. :900.0
From the automatic R summary we see that all parameters are within the expected ranges.
However looking at each variable is not the most suitable way to validate data: we can perform a simple automatic check on the allowed ranges.
# Range check - flag values outside expected physiological limits
ranges <- list(
age = c(18, 100),
criq = c(70, 130),
n400_amp = c(0.5, 8),
n400_lat = c(300, 550),
p300_amp = c(3, 25),
p300_lat = c(250, 550),
apacs = c(0, 1),
att_matr = c(0, 60),
moca = c(0, 30)
)
cat("Range checks:\n")
any_warn <- FALSE
for (col in names(ranges)) {
if (!(col %in% colnames(df))) next
lo <- ranges[[col]][1]; hi <- ranges[[col]][2]
vals <- df[[col]][!is.na(df[[col]])]
below <- sum(vals < lo); above <- sum(vals > hi)
if (below > 0) { cat(sprintf(" ⚠ %s: %d values below %s\n", col, below, lo)); any_warn <- TRUE }
if (above > 0) { cat(sprintf(" ⚠ %s: %d values above %s\n", col, above, hi)); any_warn <- TRUE }
}
if (!any_warn) cat(" All values within expected ranges.\n")
Range checks:
⚠ age: 1 values above 100 ⚠ n400_amp: 1 values above 8 ⚠ p300_lat: 2 values above 550
1.3 Missing data¶
Missing data are common in clinical neuropsychological datasets — they can arise from fatigue, incomplete testing, or participant dropout. Strategies include:
- Listwise deletion (exclude participants with any missing value) — simple, but reduces sample size and may bias estimates if data are not missing completely at random (MCAR).
- Multiple imputation (e.g., MICE, see the
micepackage) — preserves sample size and is valid under the weaker assumption of missing at random (MAR). Generally preferred in clinical samples.
Our synthetic dataset has complete data for all 40 participants, so no imputation is needed here. The code below would work with either scenario.
Visual inspection of missing data patterns:
if (sum(is.na(df)) == 0) {
fig_missing <- ggplot() +
annotate("text", x = 0, y = 0, label = "No missing data", size = 6) +
theme_void() + ggtitle("Missing data matrix")
} else {
miss_long <- as.data.frame(is.na(df)) %>%
mutate(row = row_number()) %>%
pivot_longer(-row, names_to = "variable", values_to = "missing")
fig_missing <- ggplot(miss_long, aes(variable, row, fill = missing)) +
geom_tile() +
scale_fill_manual(values = c(`FALSE` = "grey90", `TRUE` = "firebrick")) +
theme(axis.text.y = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1)) +
ggtitle("Missing data matrix")
}
print(fig_missing)
cat(sprintf("%-25s %8s %8s\n", "Variable", "Missing", "Rate"))
cat(strrep("-", 43), "\n")
for (col in colnames(df)) {
n_miss <- sum(is.na(df[[col]]))
cat(sprintf("%-25s %8d %7.1f%%\n", col, n_miss, n_miss / nrow(df) * 100))
}
Variable Missing Rate
-------------------------------------------
subject_id 0 0.0% age 0 0.0% sex 0 0.0% criq 1 2.2% apacs 0 0.0% att_matr 0 0.0% moca 1 2.2% n400_amp 1 2.2% n400_lat 0 0.0% p300_amp 0 0.0% p300_lat 0 0.0%
What do we see in the plot?
Each row is a participant, each column a variable; a highlighted cell would mark a missing value. Here the matrix is empty and the per-column rates above all read 0% — this synthetic dataset has no missing data by construction, so there's nothing to impute or exclude. In a real dataset, scattered missingness would support listwise deletion; missingness concentrated in specific participants or variables would be a reason to consider imputation instead.
1.4 Distribution check¶
More than just plotting the data we can already perform a quality check by looking at the Q-Q plots, and running the Shapiro-Wilk normality test. This guides later decisions about transformations and correlation methods.
vars <- c("age", "sex", "criq", "apacs", "att_matr", "moca", "n400_amp", "n400_lat", "p300_amp", "p300_lat")
make_hist <- function(col) {
v <- df[[col]]; v <- v[!is.na(v)]
ggplot(data.frame(v = v), aes(v)) +
geom_histogram(bins = 15, fill = "grey70", color = "white") +
geom_vline(xintercept = mean(v), color = "firebrick", linetype = "dashed") +
geom_vline(xintercept = median(v), color = "darkgreen", linetype = "dotted") +
labs(title = sprintf("%s (n=%d)", col, length(v)), x = col, y = "Frequency") +
theme_minimal(base_size = 9)
}
make_qq <- function(col) {
v <- df[[col]]; v <- v[!is.na(v)]
ggplot(data.frame(v = v), aes(sample = v)) +
stat_qq(size = 1) + stat_qq_line(color = "firebrick") +
labs(title = "Q-Q plot", x = "Theoretical", y = "Sample") +
theme_minimal(base_size = 9)
}
hists <- lapply(vars, make_hist)
qqs <- lapply(vars, make_qq)
# 4-row x 5-col layout: hist(1-5), qq(1-5), hist(6-10), qq(6-10)
fig_dist <- wrap_plots(c(hists[1:5], qqs[1:5], hists[6:10], qqs[6:10]), ncol = 5)
print(fig_dist)
What do we see in the plot?
Histogram — a roughly symmetric, bell-shaped curve suggests normality. Strong skewness (long tail to one side) or a U-shape would be a warning sign.
Q-Q plot — points should fall along the diagonal line if the data are normally distributed. Key departures:
Both ends curve away from the line → heavy tails (more extreme values than a normal distribution would predict)
S-shaped curve → skewed distribution (left-skewed if the curve goes below then above the line; right-skewed if the reverse)
Points deviating in the middle only → unusual peakedness (kurtosis)
cat("Normality tests (Shapiro-Wilk):\n")
for (col in vars) {
v <- df[[col]][!is.na(df[[col]])]
if (length(v) < 3) next
st <- shapiro.test(v)
cat(sprintf(" %s: n=%d, mean=%.2f, SD=%.2f, W=%.4f, p=%.4f, skew=%.3f, kurt=%.3f\n",
col, length(v), mean(v), sd(v), st$statistic, st$p.value,
skew(v, type = 2), kurtosi(v, type = 2)))
}
Normality tests (Shapiro-Wilk):
age: n=45, mean=66.91, SD=14.26, W=0.5593, p=0.0000, skew=4.619, kurt=26.959 sex: n=45, mean=0.49, SD=0.51, W=0.6368, p=0.0000, skew=0.046, kurt=-2.093 criq: n=44, mean=95.82, SD=11.20, W=0.9700, p=0.3023, skew=-0.282, kurt=-0.531 apacs: n=45, mean=0.69, SD=0.13, W=0.9747, p=0.4235, skew=-0.338, kurt=-0.021 att_matr: n=45, mean=31.53, SD=8.99, W=0.9627, p=0.1554, skew=-0.547, kurt=1.094 moca: n=44, mean=21.77, SD=5.89, W=0.9533, p=0.0727, skew=-0.341, kurt=-0.535 n400_amp: n=44, mean=3.98, SD=3.35, W=0.3471, p=0.0000, skew=5.986, kurt=38.181 n400_lat: n=45, mean=427.28, SD=20.07, W=0.9330, p=0.0120, skew=0.250, kurt=-1.052 p300_amp: n=45, mean=12.22, SD=1.79, W=0.9807, p=0.6467, skew=0.252, kurt=0.530 p300_lat: n=45, mean=412.36, SD=93.35, W=0.7051, p=0.0000, skew=3.392, kurt=16.609
Interpretation: all Shapiro-Wilk tests are non-significant (p > 0.05) and skewness values are well below |0.6|, so the variables are approximately normal. Minor scatter in the Q-Q plots is expected with n = 40 and does not indicate a problem.
1.5 Outlier detection¶
Flag values beyond ±3 standard deviations. We inspect flagged cases and note whether they are data-entry errors or genuine extremes.
out_mask <- as.data.frame(matrix(FALSE, nrow = nrow(df), ncol = length(vars),
dimnames = list(NULL, vars)))
for (col in vars) {
v <- df[[col]]
valid_idx <- which(!is.na(v))
vv <- v[valid_idx]
if (length(vv) < 3) next
z <- abs(scale(vv)[, 1])
out_mask[valid_idx[z > 3], col] <- TRUE
}
cat("Outliers (±3 SD):\n")
found <- FALSE
for (col in vars) {
flagged <- which(out_mask[[col]])
if (length(flagged) > 0) {
found <- TRUE
cat(sprintf(" %s: rows %s\n", col, paste(flagged, collapse = ", ")))
for (idx in flagged) cat(sprintf(" row %d: %s = %s\n", idx, col, df[[col]][idx]))
}
}
if (!found) cat(" None detected.\n")
Outliers (±3 SD):
age: rows 24
row 24: age = 150
att_matr: rows 43
row 43: att_matr = 3
n400_amp: rows 19
row 19: n400_amp = 25
p300_lat: rows 36
row 36: p300_lat = 900
Decision rule: If clinically suspicious outliers exist, we run the analysis with and without them as a sensitivity check and report both.
1.6 Cleaning dataset¶
Let's remove the missing data and the outliers.
There are three outliers that are not plausible data: Age 150, N400 Amplitude 25 uV, P300 latency 900ms. We will remove these subjects.
However the value of 3 for attentive matrices is very low but it is also very plausible in a patient dataset, therefore we will keep this outlier.
missing_mask <- apply(is.na(df[, vars]), 1, any)
# Outliers: exclude only the implausible ones found in Section 1.5 - keep
# the att_matr outlier, since it is a plausible (if low) clinical value
implausible_cols <- c("age", "n400_amp", "p300_lat")
outlier_mask <- apply(out_mask[, implausible_cols], 1, any)
exclude_mask <- missing_mask | outlier_mask
cat(sprintf("Excluding %d participants:\n", sum(exclude_mask)))
print(df[exclude_mask, c("subject_id", vars)], row.names = FALSE)
df <- df[!exclude_mask, ]
rownames(df) <- NULL
cat(sprintf("\nRemaining: %d participants\n", nrow(df)))
Excluding 5 participants:
subject_id age sex criq apacs att_matr moca n400_amp n400_lat p300_amp
sub-016 72 0 NA 0.48 22 18 3.751 460.1 13.469
sub-019 58 1 103 0.72 27 30 25.000 434.4 8.635
sub-024 150 0 100 0.66 31 21 2.549 429.2 10.306
sub-036 62 0 74 0.67 25 20 3.626 456.5 11.983
sub-044 57 1 95 0.75 36 NA NA 397.3 12.206
p300_lat
430.1
353.4
483.4
900.0
363.9
Remaining: 40 participants
2. Analysis¶
2.1 Descriptive statistics (Table 1)¶
We describe each variable by looking at its mean, standard deviation, and range of values
table1 <- bind_rows(lapply(vars, function(col) {
v <- df[[col]][!is.na(df[[col]])]
data.frame(Variable = col, n = length(v),
Value = sprintf("%.2f ± %.2f", mean(v), sd(v)),
Range = sprintf("[%.2f, %.2f]", min(v), max(v)))
}))
table1
| Variable | n | Value | Range |
|---|---|---|---|
| <chr> | <int> | <chr> | <chr> |
| age | 40 | 65.30 ± 6.63 | [49.00, 82.00] |
| sex | 40 | 0.50 ± 0.51 | [0.00, 1.00] |
| criq | 40 | 96.10 ± 11.15 | [74.00, 119.00] |
| apacs | 40 | 0.70 ± 0.13 | [0.39, 0.93] |
| att_matr | 40 | 31.95 ± 9.30 | [3.00, 48.00] |
| moca | 40 | 21.73 ± 6.00 | [8.00, 30.00] |
| n400_amp | 40 | 3.51 ± 0.87 | [1.79, 5.83] |
| n400_lat | 40 | 426.25 ± 19.48 | [397.00, 460.10] |
| p300_amp | 40 | 12.33 ± 1.77 | [8.63, 17.03] |
| p300_lat | 40 | 400.63 ± 57.50 | [292.90, 553.70] |
2.2 Bivariate correlations (ERP × NPSY)¶
We compute Pearson correlations between all the variables, but the interesting values are the ones between the ERP features and the neuropsychological outcomes. The heatmap shows r values with significance stars (* p<0.05, ** p<0.01, *** p<0.001).
n_vars <- length(vars)
valid <- df[vars] %>% drop_na()
corr_mat <- matrix(1, n_vars, n_vars, dimnames = list(vars, vars))
p_mat <- matrix(1, n_vars, n_vars, dimnames = list(vars, vars))
for (i in 1:n_vars) {
for (j in 1:n_vars) {
if (i == j) next
ct <- cor.test(valid[[vars[i]]], valid[[vars[j]]])
corr_mat[i, j] <- unname(ct$estimate)
p_mat[i, j] <- ct$p.value
}
}
# long format, lower triangle + diagonal only
heat_rows <- list()
for (i in 1:n_vars) {
for (j in 1:n_vars) {
if (j > i) next
star <- if (p_mat[i, j] < 0.001) "***" else if (p_mat[i, j] < 0.01) "**" else if (p_mat[i, j] < 0.05) "*" else ""
heat_rows[[length(heat_rows) + 1]] <- data.frame(
Var1 = vars[i], Var2 = vars[j], r = corr_mat[i, j],
label = sprintf("%.3f\n%s", corr_mat[i, j], star))
}
}
heat_df <- bind_rows(heat_rows) %>%
mutate(Var1 = factor(Var1, levels = vars), Var2 = factor(Var2, levels = rev(vars)))
fig_corr <- ggplot(heat_df, aes(Var1, Var2, fill = r)) +
geom_tile(color = "white") +
geom_text(aes(label = label), size = 2.6, lineheight = 0.85) +
scale_fill_gradient2(low = "#2166ac", mid = "white", high = "#b2182b", midpoint = 0, limits = c(-1, 1)) +
coord_fixed() +
labs(title = "Pearson correlations (* p<0.05, ** p<0.01, *** p<0.001)", x = NULL, y = NULL) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(fig_corr)
What do we see in the plot?
Each cell shows the Pearson correlation coefficient (r) between a pair of variables. Colours indicate direction: red = positive, blue = negative. Intensity reflects strength (closer to ±1 = stronger).
Pearson correlation - Values range from −1 (perfect negative) to +1 (perfect positive). A value near 0 means no linear relationship.
Significance - Stars indicate whether the correlation is statistically significant: * p<0.05, ** p<0.01, *** p<0.001. The lower the p-value, the stronger the evidence that the correlation is not due to chance.
Key hypothesised associations¶
Based on the literature, we focus on two key pairs and generate scatter plots with regression lines and 95 % confidence bands.
pairs_list <- list(c("p300_amp", "att_matr"), c("n400_lat", "apacs"))
make_scatter <- function(xvar, yvar) {
v <- df[c(xvar, yvar)] %>% drop_na()
ct <- cor.test(v[[xvar]], v[[yvar]])
ggplot(v, aes(.data[[xvar]], .data[[yvar]])) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", color = "#00796B", fill = "#00796B", alpha = 0.15, formula = y ~ x) +
labs(title = sprintf("%s vs %s r=%.3f, p=%.4f", xvar, yvar, ct$estimate, ct$p.value),
x = xvar, y = yvar)
}
fig_scatter <- wrap_plots(lapply(pairs_list, function(p) make_scatter(p[1], p[2])), ncol = 2)
print(fig_scatter)
What do we see in the plot?
Each dot shows the values of one participant.
The regression line shows the best-fit linear relationship — upward slope = positive association, flat = none.
The shaded band is the 95 % confidence interval: narrower = more precise estimate, wider = more uncertainty (typically at the extremes). If the band includes a horizontal line, the slope may not be significant.
The r value (Pearson correlation coefficient) quantifies strength (−1 to +1); p tests whether it is significantly different from zero.
2.3 General linear models¶
We fit three models, each with a different structure. This reflects the reality of clinical research: sometimes we have a clear theoretical prediction, sometimes we explore, and sometimes we learn from our mistakes and simplify.
APACS - A theory-driven model¶
Descrizione di modello teorico per APACS
APACS ~ N400_lat + N400_amp + Age
mod_apacs <- glm(apacs ~ n400_lat + n400_amp + age, data = df, family = gaussian())
summary(mod_apacs)
Call:
glm(formula = apacs ~ n400_lat + n400_amp + age, family = gaussian(),
data = df)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.7338107 0.2755661 9.921 7.68e-12 ***
n400_lat -0.0041682 0.0005948 -7.008 3.21e-08 ***
n400_amp 0.0548452 0.0132786 4.130 0.000206 ***
age -0.0069354 0.0017570 -3.947 0.000351 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.005201438)
Null deviance: 0.65508 on 39 degrees of freedom
Residual deviance: 0.18725 on 36 degrees of freedom
AIC: -91.052
Number of Fisher Scoring iterations: 2
Attentive Matrices - A model with no a priori knowledge¶
Now we want to build a model for attentive matrices.
As known by the theory, and also shown in the previous correlation plot, attentive matrices are correlated to P300, education (CRIq), and age.
But our dataset includes much more variables, therefore one idea could be bigger is better: let's put everything together
So we can try to fit a kitchen-sink model using all the data as covariates: AttentiveMatrices ~ P300_lat + P300_amp + N400_lat + N400_amp + Age + Sex + CRIq
We fit this model only to see what happens when theory is ignored and every available variable is thrown in — not to use its p-values to decide afterwards which variables to keep (see below).
mod_att_matr_full <- glm(att_matr ~ p300_lat + p300_amp + n400_lat + n400_amp + age + sex + criq,
data = df, family = gaussian())
summary(mod_att_matr_full)
Call:
glm(formula = att_matr ~ p300_lat + p300_amp + n400_lat + n400_amp +
age + sex + criq, family = gaussian(), data = df)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -51.65230 27.62387 -1.870 0.07068 .
p300_lat -0.01082 0.01827 -0.592 0.55772
p300_amp 3.37630 0.53689 6.289 4.72e-07 ***
n400_lat 0.04021 0.04604 0.873 0.38892
n400_amp -0.64179 1.02223 -0.628 0.53456
age -0.45527 0.14983 -3.039 0.00471 **
sex -0.61895 1.87509 -0.330 0.74348
criq 0.63951 0.08822 7.249 3.09e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 27.3741)
Null deviance: 3375.90 on 39 degrees of freedom
Residual deviance: 875.97 on 32 degrees of freedom
AIC: 254.97
Number of Fisher Scoring iterations: 2
Attentive Matrices - A theory-guided, parsimonious model¶
The kitchen-sink model above shows the cost of ignoring theory: several predictors are far from significant, and entering redundant or theoretically unmotivated variables inflates standard errors and destabilises the coefficients.
Importantly, we do not decide which variables to keep by looking at
these p-values — that would be post-hoc pruning, a form of circular
reasoning: the same data would first be searched to select the model and
then reused to test it, inflating the apparent evidence for whatever
survives. Instead we go back to theory and to the earlier bivariate
analysis (Section 2.2): P300 amplitude, age, and CRIq are the established
covariates for attentional performance, so those are the predictors we
keep a priori. N400 amplitude/latency and P300 latency are dropped
because there is no theoretical reason to expect them to explain attentive
matrices performance specifically, and sex is dropped on the same a
priori grounds — not because it failed to reach significance above.
Attentive Matrices ~ P300_amp + Age + CRIq
mod_att_matr_pars <- glm(att_matr ~ p300_amp + age + criq, data = df, family = gaussian())
summary(mod_att_matr_pars)
Call:
glm(formula = att_matr ~ p300_amp + age + criq, family = gaussian(),
data = df)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -39.83626 12.39673 -3.213 0.002765 **
p300_amp 3.37819 0.48259 7.000 3.28e-08 ***
age -0.46156 0.12350 -3.737 0.000644 ***
criq 0.62718 0.07566 8.289 7.25e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 25.34659)
Null deviance: 3375.90 on 39 degrees of freedom
Residual deviance: 912.48 on 36 degrees of freedom
AIC: 248.61
Number of Fisher Scoring iterations: 2
This theory-guided model has fewer predictors, and — as a consequence, not as the reason we chose them — all of them are significant and the model is more stable. Because the selection was made a priori rather than by trimming non-significant terms, the improved fit (a lower AIC, despite — in fact because of — fewer parameters) can be read as genuine evidence in favour of parsimony, rather than as an artefact of having searched the data for whichever subset happens to look best.
Takeaway
More variables does not mean a better model — and how you choose which variables to drop matters just as much as the decision to drop them. Selecting predictors after seeing that they are not significant in a bigger model is a form of post hoc pruning, which is circular: the same data end up being used both to choose the model and to test it, which inflates the apparent evidence for whatever predictors remain.
Here the predictors were instead chosen a priori, from theory and from the earlier correlation analysis (Section 2.2), before looking at any p-value from the kitchen-sink model. The kitchen-sink model was fit only to illustrate the cost of ignoring theory — inflated standard errors, an unstable intercept — not to guide which variables to keep. That the resulting theory-driven model also turns out to fit at least as well (a lower AIC, all coefficients significant, and a simpler, more interpretable structure) is a welcome confirmation, not the justification for the choice.
2.4 Assumption checking (post-hoc)¶
Fitting a model is not the last step: linear regression relies on a few
assumptions, and good practice is to check them after fitting, before
trusting the coefficients. We do this here for the final, parsimonious
Attentive Matrices model (att_matr ~ p300_amp + age + criq).
- Linearity & homoscedasticity — residuals plotted against fitted values should scatter randomly around zero, with no funnel shape or curve.
- Normality of residuals — a Q-Q plot of the residuals should follow the diagonal (same logic as the Q-Q plots in Section 1.4, now applied to the model's errors instead of the raw variables).
- Influential observations — Cook's distance flags participants who, on their own, pull the fitted model strongly in one direction. A common rule of thumb is to inspect points above 4/n.
- Multicollinearity — the Variance Inflation Factor (VIF) quantifies how
much a predictor's variance is inflated by correlation with the other
predictors. VIF > 5 is a moderate concern, VIF > 10 is severe. Unlike a
manual VIF calculation, R's
car::vif()never includes the intercept in the first place, so there is nothing to exclude here.
mod <- mod_att_matr_pars
fitted_vals <- fitted(mod)
resid_vals <- residuals(mod, type = "response")
n <- nobs(mod)
p_resid <- ggplot(data.frame(fitted = fitted_vals, resid = resid_vals), aes(fitted, resid)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey40") +
labs(title = "Residuals vs Fitted", x = "Fitted values", y = "Residuals")
p_qq <- ggplot(data.frame(r = resid_vals), aes(sample = r)) +
stat_qq() + stat_qq_line(color = "firebrick") +
labs(title = "Q-Q plot (residuals)")
cooks <- cooks.distance(mod)
p_cooks <- ggplot(data.frame(i = seq_along(cooks), cooks = cooks), aes(i, cooks)) +
geom_segment(aes(xend = i, yend = 0), color = "steelblue") +
geom_point(color = "steelblue") +
geom_hline(yintercept = 4 / n, color = "firebrick", linetype = "dashed") +
labs(title = "Cook's Distance", x = "Participant (row index)", y = "Cook's distance",
subtitle = sprintf("4/n = %.3f", 4 / n))
vif_vals <- car::vif(mod)
vif_table <- data.frame(Predictor = names(vif_vals), VIF = as.numeric(vif_vals))
p_vif <- ggplot(vif_table, aes(VIF, reorder(Predictor, VIF))) +
geom_col(fill = "steelblue") +
geom_vline(xintercept = 5, color = "orange", linetype = "dashed") +
geom_vline(xintercept = 10, color = "firebrick", linetype = "dashed") +
labs(title = "Variance Inflation Factor", x = "VIF", y = NULL)
fig_diag <- (p_resid + p_qq) / (p_cooks + p_vif)
print(fig_diag)
print(vif_table)
n_influential <- sum(cooks > 4 / n)
cat(sprintf("\nParticipants above the 4/n Cook's distance threshold: %d (max Cook's D = %.3f)\n",
n_influential, max(cooks)))
Predictor VIF 1 p300_amp 1.122159 2 age 1.031568 3 criq 1.095338
Participants above the 4/n Cook's distance threshold: 3 (max Cook's D = 0.236)
What do we see in the plot?
Residuals vs Fitted — no funnel shape or curve should be visible; a random scatter around zero supports linearity and constant variance (homoscedasticity).
Q-Q plot (residuals) — points should lie close to the diagonal; systematic curvature would suggest non-normal errors.
Cook's distance — a handful of points above the 4/n line is not unusual with n = 40 and, on its own, is not a red flag: what matters is whether any point stands out with a much larger value than the rest (a useful check is to refit the model without that participant and see whether the coefficients change materially).
VIF — values close to 1 mean each predictor carries mostly independent information; here
p300_amp,age, andcriqall sit well below the VIF = 5 caution threshold, so multicollinearity is not a concern for this model.
Interpretation: for the parsimonious Attentive Matrices model, the diagnostics are reassuring — residuals show no obvious pattern, the Q-Q plot is close to the diagonal, no participant has a disproportionately large Cook's distance, and VIF values are low. This supports treating the model's coefficients (Section 2.3) at face value.
3. Export¶
We now save the results that matter most for a clinical audience: the key
figures (as PNG images, ready to paste into a report or a slide) and
the key tables (as both .csv — universal, opens anywhere — and
.xlsx — a single Excel workbook with one sheet per table, easier to
browse for those who are not used to statistical software).
Nothing here is recomputed: we simply save the figures and tables already produced above.
3.1 Figures¶
We save the figures produced earlier, one PNG per figure.
figures_to_save <- list(
"A_missingness.png" = list(fig = fig_missing, w = 10, h = 3),
"B_distributions.png" = list(fig = fig_dist, w = 20, h = 16),
"C_correlations.png" = list(fig = fig_corr, w = 9, h = 7),
"D_key_associations.png" = list(fig = fig_scatter, w = 12, h = 5),
"E_diagnostics_attentive_matrices.png" = list(fig = fig_diag, w = 11, h = 9)
)
for (filename in names(figures_to_save)) {
spec <- figures_to_save[[filename]]
path <- file.path(FIGS, filename)
ggsave(path, plot = spec$fig, width = spec$w, height = spec$h, dpi = 150, bg = "white")
cat(sprintf("Saved: %s\n", path))
}
Saved: ../outputs/figures/A_missingness.png
Saved: ../outputs/figures/B_distributions.png
Saved: ../outputs/figures/C_correlations.png
Saved: ../outputs/figures/D_key_associations.png
Saved: ../outputs/figures/E_diagnostics_attentive_matrices.png
3.2 Tables¶
We export four tables:
- Descriptive statistics (
table1, already built in Section 2.1) - Bivariate correlations — every pair of variables, one row per pair (long/tidy format, easier to filter and sort in Excel than a matrix)
- Regression results — coefficients, standard errors, p-values and 95% CI for all three models (APACS, and both Attentive Matrices models)
- Multicollinearity check (VIF) — for the parsimonious Attentive Matrices model (Section 2.4)
Each table gets a Significant (p<0.05) column so it is immediately clear,
without reading p-values, which associations/predictors reach conventional
significance.
# Table 2 - bivariate correlations, one row per unique pair of variables
corr_rows <- list()
for (i in 1:n_vars) {
for (j in 1:n_vars) {
if (j <= i) next
corr_rows[[length(corr_rows) + 1]] <- data.frame(
`Variable 1` = vars[i], `Variable 2` = vars[j],
r = corr_mat[i, j], `p-value` = p_mat[i, j],
`Significant (p<0.05)` = p_mat[i, j] < 0.05,
check.names = FALSE)
}
}
table2_correlations <- bind_rows(corr_rows)
head(table2_correlations)
| Variable 1 | Variable 2 | r | p-value | Significant (p<0.05) | |
|---|---|---|---|---|---|
| <chr> | <chr> | <dbl> | <dbl> | <lgl> | |
| 1 | age | sex | -0.32076961 | 0.043583785 | TRUE |
| 2 | age | criq | 0.02871583 | 0.860378823 | FALSE |
| 3 | age | apacs | -0.44225738 | 0.004271818 | TRUE |
| 4 | age | att_matr | -0.20633880 | 0.201454655 | FALSE |
| 5 | age | moca | -0.23297195 | 0.147977000 | FALSE |
| 6 | age | n400_amp | -0.10502923 | 0.518935561 | FALSE |
# Table 3 - regression coefficients for all three fitted models
model_to_table <- function(model, outcome_label) {
co <- summary(model)$coefficients
ci <- confint.default(model) # Wald CI
data.frame(
Outcome = outcome_label,
Predictor = rownames(co),
Estimate = co[, "Estimate"],
`Std. Error` = co[, "Std. Error"],
`p-value` = co[, 4],
CI_lower = ci[, 1],
CI_upper = ci[, 2],
`Significant (p<0.05)` = co[, 4] < 0.05,
row.names = NULL, check.names = FALSE
)
}
table3_regression <- bind_rows(
model_to_table(mod_apacs, "APACS (theory-driven)"),
model_to_table(mod_att_matr_full, "Attentive Matrices (kitchen-sink)"),
model_to_table(mod_att_matr_pars, "Attentive Matrices (parsimonious)")
)
table3_regression
| Outcome | Predictor | Estimate | Std. Error | p-value | CI_lower | CI_upper | Significant (p<0.05) |
|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <lgl> |
| APACS (theory-driven) | (Intercept) | 2.733810687 | 0.275566134 | 7.681474e-12 | 2.193711e+00 | 3.273910384 | TRUE |
| APACS (theory-driven) | n400_lat | -0.004168185 | 0.000594790 | 3.206130e-08 | -5.333952e-03 | -0.003002418 | TRUE |
| APACS (theory-driven) | n400_amp | 0.054845209 | 0.013278574 | 2.057195e-04 | 2.881968e-02 | 0.080870736 | TRUE |
| APACS (theory-driven) | age | -0.006935368 | 0.001756969 | 3.514614e-04 | -1.037896e-02 | -0.003491772 | TRUE |
| Attentive Matrices (kitchen-sink) | (Intercept) | -51.652297835 | 27.623873752 | 7.067881e-02 | -1.057941e+02 | 2.489499832 | FALSE |
| Attentive Matrices (kitchen-sink) | p300_lat | -0.010822586 | 0.018268045 | 5.577243e-01 | -4.662730e-02 | 0.024982124 | FALSE |
| Attentive Matrices (kitchen-sink) | p300_amp | 3.376296480 | 0.536885519 | 4.718415e-07 | 2.324020e+00 | 4.428572761 | TRUE |
| Attentive Matrices (kitchen-sink) | n400_lat | 0.040214361 | 0.046040257 | 3.889198e-01 | -5.002289e-02 | 0.130451606 | FALSE |
| Attentive Matrices (kitchen-sink) | n400_amp | -0.641793315 | 1.022226428 | 5.345648e-01 | -2.645320e+00 | 1.361733668 | FALSE |
| Attentive Matrices (kitchen-sink) | age | -0.455266606 | 0.149828346 | 4.706145e-03 | -7.489248e-01 | -0.161608445 | TRUE |
| Attentive Matrices (kitchen-sink) | sex | -0.618947605 | 1.875090349 | 7.434837e-01 | -4.294057e+00 | 3.056161946 | FALSE |
| Attentive Matrices (kitchen-sink) | criq | 0.639508702 | 0.088217952 | 3.092328e-08 | 4.666047e-01 | 0.812412710 | TRUE |
| Attentive Matrices (parsimonious) | (Intercept) | -39.836264061 | 12.396731475 | 2.765180e-03 | -6.413341e+01 | -15.539116844 | TRUE |
| Attentive Matrices (parsimonious) | p300_amp | 3.378190908 | 0.482593783 | 3.282095e-08 | 2.432324e+00 | 4.324057342 | TRUE |
| Attentive Matrices (parsimonious) | age | -0.461561250 | 0.123495924 | 6.436733e-04 | -7.036088e-01 | -0.219513686 | TRUE |
| Attentive Matrices (parsimonious) | criq | 0.627177612 | 0.075662521 | 7.246298e-10 | 4.788818e-01 | 0.775473428 | TRUE |
# Table 4 - VIF table computed in Section 2.4 (Attentive Matrices, parsimonious model)
table4_vif <- vif_table
table4_vif
| Predictor | VIF |
|---|---|
| <chr> | <dbl> |
| p300_amp | 1.122159 |
| age | 1.031568 |
| criq | 1.095338 |
tables_to_save <- list(
"table1_descriptives.csv" = table1,
"table2_correlations.csv" = table2_correlations,
"table3_regression.csv" = table3_regression,
"table4_vif.csv" = table4_vif
)
for (filename in names(tables_to_save)) {
path <- file.path(RESULTS, filename)
write.csv(tables_to_save[[filename]], path, row.names = FALSE)
cat(sprintf("Saved: %s\n", path))
}
wb <- createWorkbook()
addWorksheet(wb, "Descriptives"); writeData(wb, "Descriptives", table1)
addWorksheet(wb, "Correlations"); writeData(wb, "Correlations", table2_correlations)
addWorksheet(wb, "Regression"); writeData(wb, "Regression", table3_regression)
addWorksheet(wb, "VIF"); writeData(wb, "VIF", table4_vif)
xlsx_path <- file.path(RESULTS, "CaseStudy_1_results.xlsx")
saveWorkbook(wb, xlsx_path, overwrite = TRUE)
cat(sprintf("Saved: %s\n", xlsx_path))
Saved: ../outputs/results/table1_descriptives.csv Saved: ../outputs/results/table2_correlations.csv Saved: ../outputs/results/table3_regression.csv Saved: ../outputs/results/table4_vif.csv
Saved: ../outputs/results/CaseStudy_1_results.xlsx
4. Strengths and limitations of this approach (Section 4.2.3)¶
Strengths:
- Reduces high-dimensional EEG to 4 interpretable scalar features per participant — theoretically grounded, easy to communicate
- Uses well-established statistical models (OLS regression) familiar to clinical researchers
- The N400 and P300 have decades of validation across neurological populations (Kutas & Federmeier, 2011; Polich, 2007)
Limitations:
- Strong a priori decisions about time windows and electrode sets may miss atypical or distributed neural patterns — especially relevant in stroke patients whose ERP topography may differ from normative expectations
- Trial-to-trial variability is discarded (only mean amplitude and peak latency are kept), which may itself carry clinically relevant information
- The hypothesis-driven nature of the approach means it will not detect effects outside the predefined features — for broader exploration, see Case Study 2 (spectral analysis)
For a full discussion, see Section 4.2 of the companion paper.
Generated for the EEG-NPSY tutorial. All code is self-contained for transparency and reproducibility.