Sparse PCA for mixed data

Introduction

The function sparsePCAmix in the sparsePCA package implements sparse principal component analysis (PCA) for mixed data combining numerical and categorical variables. This implementation, as well as that of the function sparsePCA, which performs sparse and group-sparse PCA for numerical data, relies on either a block optimization algorithm or a deflation-based algorithm. The variance of the possibly non-orthogonal principal components is computed using the optimal projected variance.

library(sparsePCA)

Arguments

Several arguments may be passed to sparsePCAmix, but the main ones are:

X is a data frame containing numerical and/or categorical columns.
m is the number of principal components to compute.
lambda is the vector of sparsity parameters for the \(m\) components. In practice, the same sparsity parameter in \([0, 1]\) is chosen for all dimensions (here, the sparsity parameters are normalized to promote more balanced sparsification across components).
block is a logical value. If block = TRUE (the default), a block optimization algorithm is used to compute the m components simultaneously. Otherwise, a deflation-based algorithm is used to compute the components recursively.
mu is a vector of \(m\) strictly positive weights (required only for the block algorithm). Strictly decreasing weights recover the standard PCA solution when all sparsity parameters are set to zero. With constant weights, the standard PCA solution is not necessarily recovered; instead, an orthogonal rotation of it may be obtained. By default, \(\mu_k = 1/k, \, k=1,\dots,m\).

Numerical comparisons of the three approaches (deflation, block with equal weights, and block with different weights) showed better and more robust results for the block algorithm with different weights. This is why these are the default settings of the function arguments.

Output

The sparsePCAmix function returns an object of class sparsePCA. This object mainly contains the matrix of sparse loadings and the matrix of principal component scores. Two methods are associated with this class: the plot method, which displays the observations according to their principal component scores, together with the proportion of variance explained by each component, and the predict method, which predicts the principal component scores for new observations.

Heart disease data

The HDdata dataset is a simplified version of the heart disease database available from the UCI Machine Learning Repository¹. It contains data on 270 patients described by a mixture of numerical and categorical variables. The six numerical variables are:

age
resting blood pressure in mm Hg
serum cholesterol in mg/dl
maximum heart rate achieved
oldpeak = ST depression induced by exercise relative to rest
number of major vessels (0-3) colored by fluoroscopy

The seven categorical variables are:

sex (0 = female, 1 = male)
chest pain type (1 = typical angina, 2 = atypical angina, 3 = non-anginal pain, 4 = asymptomatic)
fasting blood sugar > 120 mg/dl (1 = True, 0 = False)
resting electrocardiographic results (0 = normal, 1 = ST-T wave abnormality, 2 = elevation or depression of > 0.05 mV)
exercise-induced angina (1 = Yes, 0 = No)
slope of the peak exercise ST segment (1 = upsloping, 2 = flat, 3 = downsloping)
thal: 3 = normal; 6 = fixed defect; 7 = reversible defect
heart disease: 1 = absent, 2 = present

The last binary variable, heart disease, is the variable of interest and is therefore not used to compute the principal components.

data("HDdata")

PCAmix and the choice of the number of components

The method PCAmix, implemented in the R package PCAmixdata, is first applied to the heart disease dataset. This method is designed for PCA of mixed data (Chavent et al. 2022) and has the interesting property of reducing to standard PCA when all variables are numerical and to MCA (Multiple Correspondence Analysis) when all variables are categorical.

library(PCAmixdata)

# The binary variable `heart_disease` is removed

X <- HDdata[, -14]

# The dataset is split into numerical and categorical datasets

X.quanti <- splitmix(X)$X.quanti
X.quali <- splitmix(X)$X.quali

# PCAmix is applied to these two datasets

out <- PCAmix(X.quanti, X.quali,
              rename.level = TRUE, graph = FALSE)

barplot(out$eig[, 1], names = 1:18,
        xlab = "dim", ylab = "eigenvalue")

The eigenvalues are the variances of the principal components and the barplot of the eigenvalues suggests to choose \(m=3\) components.

# Proportion of explained variance.
m <- 3
out$eig[1:m, "Proportion"]

##     dim 1     dim 2     dim 3 
## 17.868093  9.281815  8.259613

# Loadings 
out$V[, 1:m]

##                              dim 1        dim 2       dim 3
## age                    0.271129454  0.387807928 -0.11624299
## blood_pressure         0.169782824  0.397578481  0.25529942
## serum_cholestoral      0.092568756  0.374719232 -0.27641166
## max_heart_rate        -0.381266791  0.032870017  0.20574799
## oldpeak                0.381499745  0.032754807  0.21668747
## number_vessels         0.282099424  0.073924728 -0.13878254
## sex=0                 -0.061117826  0.217212679 -0.12943637
## sex=1                  0.061117826 -0.217212679  0.12943637
## chest_pain=1           0.003678278  0.029554954  0.11641504
## chest_pain=2          -0.085709257 -0.008834936 -0.02068031
## chest_pain=3          -0.076300669  0.076972971  0.03483134
## chest_pain=4           0.158331648 -0.097692989 -0.13056606
## fasting_blood_sugar=0 -0.012744094 -0.078519107 -0.11968921
## fasting_blood_sugar=1  0.012744094  0.078519107  0.11968921
## resting_results=0     -0.078118847 -0.100240552 -0.01962458
## resting_results=1      0.004554606  0.020241835 -0.01189911
## resting_results=2      0.073564240  0.079998716  0.03152369
## induced_angina=0      -0.155537784  0.092568507  0.05440910
## induced_angina=1       0.155537784 -0.092568507 -0.05440910
## slope=1               -0.182767198  0.018092463 -0.01504295
## slope=2                0.149997574 -0.043115563 -0.10912250
## slope=3                0.032769624  0.025023099  0.12416545
## thal=3                -0.164223908  0.142491801 -0.06352670
## thal=6                 0.020544628 -0.026394079  0.02980747
## thal=7                 0.143679280 -0.116097722  0.03371923

The objective is now to use sparsePCAmix to construct principal components with sparse loadings while preserving as much of the explained variance as possible. Here, the first three principal components explain 17.8%, 9.2%, and 8.2% of the total variance, respectively.

Choice of the common sparsity parameters and regularization paths

Once the number of components, \(m\), has been selected, a common sparsity parameter \(\lambda\) must be chosen. For that purpose, we jointly examine the regularization path of the proportion of explained variance and the regularization path of the number of selected variables.

A one-dimensional grid of common values for the sparsity parameter is first built: \[0 \leq \lambda = \lambda_1 = \lambda_2 = \lambda_3 \leq 1\] by letting \(\lambda\) vary from 0 to 1 in steps of 0.01.

grid <- seq(0, 1, by = 0.01) # grid of common lambda values

The loading matrices are computed over this one-dimensional grid, together with the corresponding proportions of variance explained (computed using the optimal projected variance) and numbers of selected variables for each component.

m <- 3
pev <- matrix(NA, length(grid), m) # proportion of explained variance
nsel <- matrix(NA, length(grid), m) # number of selected variables


for (k in 1:length(grid))
{
  res <- sparsePCAmix(X = X,
                      lambda = rep(grid[k], m),
                      m = m, 
                      rename.level = TRUE)
  pev[k,] <-res$pev
  nsel[k,] <- res$nsel
}
colnames(pev) <- colnames(nsel) <- paste("dim", 1:m, sep = "")

The figure below plots the proportion of variance explained by the principal components (pev) as a function of the common sparsity parameter \(\lambda\).

k <- 33
matplot(grid, pev,
        xlab = "lambda",
        ylab = "pev",
        type = "l", 
        col = 1:m, 
        lty = 1) 

legend("topright",
       legend = colnames(pev),
       lty = 1,
       col = 1:m,
       bty = "n")

abline(v = grid[k], lty = 2)

The next figure plots the number of selected variables (nsel) as a function of \(\lambda\).

matplot(grid, nsel,
        ylab = "nsel",
        xlab = "lambda",
        type="l", 
        col=1:m,
        lty=1) 

legend("topright", 
       legend = colnames(pev),
       lty  =1,
       col = 1:m,
       bty = "n")

abline(v = grid[k], lty = 2)

These two regularization paths suggest to choose \(\lambda=0.32\) (vertical dotted line).

Some results of sparsePCAmix

The function sparsePCAmix is now applied using a common sparsity parameter, \(\lambda = 0.32\), across the three dimensions.

# sparsePCAmix with lambda = c(0.32, 0.32, 0.32)
m <- 3
out <- sparsePCAmix(X = X,
                    lambda = rep(0.32, m),
                    m = m, 
                    rename.level = TRUE)

The levels of a categorical variable form a group, and their loadings are set to zero simultaneously. For instance, the categorical variable fasting_blood_sugar is not selected in the first dimension because the loadings of all its levels are equal to zero. In contrast, the loadings of the levels of the variable chest_pain are non-zero, and this categorical variable is therefore selected in the first dimension.

# The selected numerical or categorical variables
out$sel

## $dim1
## [1] "max_heart_rate" "oldpeak"        "chest_pain"     "induced_angina"
## [5] "slope"          "thal"          
## 
## $dim2
## [1] "age"               "blood_pressure"    "serum_cholestoral" "sex"              
## 
## $dim3
## [1] "fasting_blood_sugar" "slope"

The four principal components built from these sparse loadings explain 29.2% of the variance, compared with 35.4% for PCAmix.

# Proportion of explained variance
out$pev

##       dim1       dim2       dim3 
## 0.14927345 0.08064836 0.06292405

The first principal component, built from two numerical variables (max_heart_rate and oldpeak) and four categorical variables (chest_pain, induced_angina, slope, and thal), still discriminates between the presence and the absence of heart disease in the 270 patients.

# Plot of the patients 
plot(out,
     labels = TRUE,
     color_by = HDdata$heart_disease,
     col = c("steelblue", "tomato"))

A method predict uses the coefficients of the principal components of sparsePCAmix, to predict the scores of new observations.