Principal Component Analysis (PCA) transforms a set of correlated signals (e.g. EEG channels) into the same number of uncorrelated signals, the principal components (PCs), ordered by the variance they explain: the first PC is the combination of the channels with the largest variance, the second the combination with the largest remaining variance that is uncorrelated with the first, and so on.
EEG channels are strongly correlated (volume conduction, common reference), so a few PCs usually explain most of the variance. PCA is used to:
reduce the dimensionality of the data, e.g. before ICA or machine learning,
remove components dominated by noise or artifacts (with caution, see below),
detect rank deficiency (channels that are linear combinations of other channels).
How PCA Works
The channels are centered (their means are subtracted).
The covariance matrix of the channels is computed.
Its eigenvectors are the PC weights (the direction of each PC in the channel space) and its eigenvalues the variances of the PCs.
The PCs are the projections of the data on the eigenvectors: \(\mathrm{PC} = W^T (x - \bar{x})\).
The data are not scaled: channels with larger amplitudes contribute more to the first PCs. The data can be reconstructed from the PCs: \(x = W \, \mathrm{PC} + \bar{x}\), exactly with all PCs and approximately with the first \(k\).
Decomposition
pca_decompose() fits one PCA model to all epochs of the selected channels and returns the first n PCs:
pc, pcv, pcm, pc_model =pca_decompose(eeg; ch ="eeg", n =19)
pc: PCs × samples × epochs,
pcv: variance of each PC in each epoch, in % of the total variance of the epoch,
pcm: channel means,
pc_model: the fitted model (MultivariateStats.PCA), used to reconstruct the signal.
cv =cumsum(pcv[:, 1])for k in [1, 3, 5, 10]println("first ", rpad(k, 2), " PCs: ", round(cv[k]; digits =1), "%")end
first 1 PCs: 65.1%
first 3 PCs: 87.2%
first 5 PCs: 92.4%
first 10 PCs: 97.8%
The first PC explains about two thirds of the variance, the first 5 PCs more than 90%: the channels are strongly correlated.
PC weights
The weights of a PC (a column of the projection matrix) show how much each channel contributes to it; as topographic maps:
W = NeuroAnalyzer.MultivariateStats.projection(pc_model)plot_topo(eeg; data = W[:, 1], ch ="eeg", title ="PC 1", ps =:m, cb_title ="Weight", nmethod =:none)
plot_topo(eeg; data = W[:, 2], ch ="eeg", title ="PC 2", ps =:m, cb_title ="Weight", nmethod =:none)
PC 1 has weights of the same sign in almost all channels: it is activity common to the channels (e.g. activity at the reference electrode). Fp2, F4 and Fz, which correlate poorly with the other channels in this recording (see Bad channels and epochs), hardly contribute to it. The sign of a PC is arbitrary (the PC and its weights can both be multiplied by −1).
Reconstruction
pca_reconstruct() reconstructs the signal from the PCs and the model. Keeping only the first \(k\) PCs (setting the others to zero) gives the best approximation of the data with \(k\) dimensions:
x = eeg.data[get_channel(eeg; ch ="eeg"), :, :]for k in [1, 3, 5, 10, 19] pc_k =copy(pc) pc_k[(k +1):end, :, :] .=0 e_k =pca_reconstruct(eeg; ch ="eeg", pc = pc_k, pc_model = pc_model) x_k = e_k.data[get_channel(eeg; ch ="eeg"), :, :]println(rpad(k, 3), "PCs: ",round(100* (1-var(x .- x_k) /var(x)); digits =1),"% of the variance reconstructed")end
1 PCs: 65.1% of the variance reconstructed
3 PCs: 87.2% of the variance reconstructed
5 PCs: 92.4% of the variance reconstructed
10 PCs: 97.8% of the variance reconstructed
19 PCs: 100.0% of the variance reconstructed
Removing a component
Setting a PC to zero and reconstructing removes its activity from all channels. PCs are defined by variance only, so a PC usually contains both artifacts and brain activity; removing it removes both. Artifacts are better separated by ICA, which uses statistical independence instead of variance.
pc_r =copy(pc)# remove PC 1pc_r[1, :, :] .=0eeg_r =pca_reconstruct(eeg; ch ="eeg", pc = pc_r, pc_model = pc_model)
Rank Deficiency
If a channel is a linear combination of other channels, the data have fewer independent dimensions than channels (the data are rank deficient), and the last PCs have (almost) zero variance. This happens after re-referencing to the average reference (the channels sum to zero), after interpolating channels and after removing ICA components:
eeg_avg =reference_avg(eeg)_, pcv_avg, _, _ =pca_decompose(eeg_avg; ch ="eeg", n =19)println("variance of the last PC, recorded reference: ",round(pcv[end, 1]; sigdigits =2), "%")println("variance of the last PC, average reference: ",round(pcv_avg[end, 1]; sigdigits =2), "%")
variance of the last PC, recorded reference: 0.097%
variance of the last PC, average reference: 3.5e-29%
After average referencing, one PC has no variance: the 19 channels contain only 18 independent signals. Methods that need full-rank data, such as ICA, should be given (at most) as many components as the rank of the data, e.g. after reducing the data with PCA.
