---
title: "Worked example: msPCA on mtcars"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Worked example: msPCA on mtcars}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

## Overview

This vignette shows the basic workflow of `msPCA` on the built-in `mtcars`
dataset. We compute sparse principal components, inspect the solution with the
`print()` and `summary()` S3 methods, and compare the sparse result with dense
PCA.

## Install and load

Install the package directly from CRAN.

```{r install, eval=FALSE}
install.packages("msPCA")
```

```{r load}
library(msPCA)
```

## Fit two sparse PCs

We work with the correlation matrix of `mtcars` and ask for two 4-sparse
principal components under the default orthogonality constraint.

```{r fit}
Sigma <- cor(mtcars)

set.seed(42)
res <- mspca(Sigma, r = 2, ks = c(4, 4), verbose = FALSE)
```

`print()` shows the sparse loading matrix restricted to the union of all
active variables, together with the percentage of variance explained and
the number of non-zero loadings per component.

```{r print}
print(res)
```

`summary()` gives a fuller breakdown: a per-PC table of variance explained,
sparsity, and each component's largest violation against any other
component, followed by the full pairwise violation matrix, which shows how
well the constraint is satisfied for each pair of components. The header
states which constraint definition is reported: the one used to fit.

```{r summary}
summary(res)
```

## Working from the raw data matrix

By default (`type = "Sigma"`) the first argument is a covariance/correlation
matrix. Set `type = "X"` to pass the raw data matrix instead (rows are
observations, columns are variables). With `type = "X"`, msPCA applies the
algorithm to the data directly: each matrix--vector product
$\boldsymbol{\Sigma}\boldsymbol{\beta} = \boldsymbol{X}^\top(\boldsymbol{X}\boldsymbol{\beta})/(n - 1)$
is computed without ever forming the $p \times p$ matrix. This is
mathematically equivalent but more scalable when $p \gg n$.

The preprocessing arguments control which matrix is implicitly used:

- `center` (default `TRUE`) subtracts column means.
- `scale` (default `TRUE`) divides by column
  standard deviations; with `scale = TRUE` the algorithm operates on the
  correlation matrix, with `scale = FALSE` on the covariance matrix.
- `divisor` selects the normalization: `"n-1"` (default, matching `cov`/`cor`)
  or `"n"`.

With `scale = TRUE` and `divisor = "n-1"`, the raw-data call targets exactly
the same problem as the correlation-matrix call above, and returns the same
solution here.

```{r fit_X}
X <- as.matrix(mtcars)

set.seed(42)
res_X <- mspca(X, r = 2, ks = c(4, 4), type = "X", scale = TRUE, verbose = FALSE)
print(res_X)
```

The same dual interface is available for the single-component `tpm()`.

## Orthogonality versus zero correlation

Sparse loading vectors are not automatically non-redundant, so `mspca()`
imposes an explicit constraint between components. The
default (`feasibilityConstraintType = 0`) enforces orthogonality of the loading
vectors. Setting `feasibilityConstraintType = 1` instead enforces zero pairwise
correlation between the resulting scores. The choice can lead to different
solutions when the variables are strongly correlated.

```{r fit_corr}
set.seed(42)
res_corr <- mspca(Sigma, r = 2, ks = c(4, 4),
                  feasibilityConstraintType = 1, verbose = FALSE)
print(res_corr)
summary(res_corr)
```

Both sets of violations are computed at fit time and stored, so the solution
can be inspected under the other definition without a refit:

```{r nonredundancy}
res_corr$nonredundancy$uncorrelatedness   # the constraint that was enforced
res_corr$nonredundancy$orthogonality      # the one that was not
```

## Diagnostics

The utility functions `feasibility_violation_off()` and
`fraction_variance_explained()` can be called directly for custom reporting or
for comparing solutions across methods — in particular for scoring loadings
that did not come from `mspca()`, such as those of a competing package, which
carry no stored diagnostics.

```{r diagnostics}
# Orthogonality and zero-correlation violations for the default solution
feasibility_violation_off(Sigma, res$x_best, feasibilityConstraintType = 0)
feasibility_violation_off(Sigma, res$x_best, feasibilityConstraintType = 1)

# The same two numbers, already stored on the fitted object
sum(res$nonredundancy$orthogonality,    na.rm = TRUE)
sum(res$nonredundancy$uncorrelatedness, na.rm = TRUE)

# Total and per-PC fraction of variance explained
fraction_variance_explained(Sigma, res$x_best)
fraction_variance_explained_perPC(Sigma, res$x_best)
```

The zero-correlation violation is normalized by the total variance
$\mathrm{tr}(\Sigma)$, so that each pairwise term is a fraction of the total
variance and the measure is comparable between a covariance matrix and the
corresponding correlation matrix; the orthogonality violation needs no such
normalization. Note that `res$feasibility_violation`, the quantity the
solver compares against `feasibilityTolerance`, also includes the diagonal norm
terms $\bigl|\,\|\boldsymbol{u}_t\|_2^2 - 1\bigr|$, and so is at least as
large as the off-diagonal diagnostic above.

## Comparison with dense PCA

The first two dense principal components explain more variance, but all
variables receive non-zero loadings.

```{r dense_pca}
pca_res <- prcomp(mtcars, scale. = TRUE)
fraction_variance_explained(Sigma, pca_res$rotation[, 1:2])
```

Sparse PCA trades explained variance for a more interpretable loading
pattern.

## Where to go next

- `vignette("case-study-snp500", package = "msPCA")` applies the same workflow
  to a 423-stock correlation matrix shipped with the package, and shows how the
  two constraint types can yield noticeably different factor structures.
- `vignette("algorithm-and-implementation", package = "msPCA")` documents the
  optimization problem, both algorithms, the implementation, and how to set
  `ks`, `feasibilityConstraintType` and the iteration budgets.
- A comparison against eight competing functions, drawn from seven packages, on
  four datasets is available on the package website, under
  [Benchmarking against other sparse PCA packages](https://jeanpauphilet.github.io/msPCA/articles/benchmarking.html).