Skip to content

PLS-DA

Partial Least Squares Discriminant Analysis with double cross-validation.

Functions

motco.stats.pls.plsda_doubleCV(X, y, cv1_splits=7, cv2_splits=8, n_repeats=30, max_components=50, random_state=1203, n_jobs=1, progress=True)

Run canonical double-nested cross-validation for PLS-DA.

For each outer fold, the inner CV averages AUROC across its V folds and selects the n_LV with the highest mean. The outer fold's selected n_LV is then evaluated on its held-out test set. Per repeat, the K outer-fold test AUROCs are aggregated by mean (and sample std), and the per-fold n_LV choices are aggregated by mode (parsimony tie-break). One final model is refit per repeat on the full input using the modal n_LV.

Parameters:

Name Type Description Default
X DataFrame

The predictor variables.

 required
y Union[DataFrame, Series]

The outcome variable.

 required
cv1_splits int

Number of folds in the CV1 (inner) loop. Default: 7.

 7
cv2_splits int

Number of folds in the CV2 (outer) loop. Default: 8.

 8
n_repeats int

Number of repeats of the K-fold outer CV. Default: 30.

 30
max_components int

Maximum number of LV to test (candidates are 1..max_components-1). Default: 50.

 50
random_state int

For reproducibility. Default: 1203.

 1203
n_jobs int

Number of parallel workers for the inner CV loop. Use -1 for all available CPUs. Default: 1 (serial).

 1
progress bool

Whether to display a tqdm progress bar over outer folds. Default: True.

 True

Returns:

Type Description
dict with keys:
  • "table": pd.DataFrame with one row per repeat and four columns "rep" (int), "LV" (int, modal across K outer folds with parsimony tie-break), "AUROC" (float, mean across K outer folds), "AUROC_std" (float, sample std across K outer folds; NaN if K < 2).
  • "models": list of length n_repeats. Each entry is a PLSRegression refit on the full input (X, one-hot(y)) with the corresponding row's modal n_LV.
Source code in src/motco/stats/pls.py 
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
def plsda_doubleCV(
    X: pd.DataFrame,
    y: Union[pd.DataFrame, pd.Series],
    cv1_splits: int = 7,
    cv2_splits: int = 8,
    n_repeats: int = 30,
    max_components: int = 50,
    random_state: int = 1203,
    n_jobs: int = 1,
    progress: bool = True,
) -> PLSDAResult:
    """
    Run canonical double-nested cross-validation for PLS-DA.

    For each outer fold, the inner CV averages AUROC across its V folds and
    selects the n_LV with the highest mean. The outer fold's selected n_LV is
    then evaluated on its held-out test set. Per repeat, the K outer-fold test
    AUROCs are aggregated by mean (and sample std), and the per-fold n_LV
    choices are aggregated by mode (parsimony tie-break). One final model is
    refit per repeat on the full input using the modal n_LV.

    Parameters
    ----------
    X: pd.DataFrame
        The predictor variables.
    y: Union[pd.DataFrame, pd.Series]
        The outcome variable.
    cv1_splits: int
        Number of folds in the CV1 (inner) loop. Default: 7.
    cv2_splits: int
        Number of folds in the CV2 (outer) loop. Default: 8.
    n_repeats: int
        Number of repeats of the K-fold outer CV. Default: 30.
    max_components: int
        Maximum number of LV to test (candidates are 1..max_components-1). Default: 50.
    random_state: int
        For reproducibility. Default: 1203.
    n_jobs: int
        Number of parallel workers for the inner CV loop. Use -1 for all
        available CPUs. Default: 1 (serial).
    progress: bool
        Whether to display a tqdm progress bar over outer folds. Default: True.

    Returns
    -------
    dict with keys:
      - "table": pd.DataFrame with one row per repeat and four columns
            "rep" (int), "LV" (int, modal across K outer folds with parsimony
            tie-break), "AUROC" (float, mean across K outer folds),
            "AUROC_std" (float, sample std across K outer folds; NaN if K < 2).
      - "models": list of length n_repeats. Each entry is a PLSRegression
            refit on the full input (X, one-hot(y)) with the corresponding
            row's modal n_LV.
    """
    _X_arr = np.asarray(X, dtype=float)
    _y_arr = np.asarray(y)
    _n_x = _X_arr.shape[0]
    _n_y = _y_arr.shape[0]
    if _n_x != _n_y:
        raise ValueError(
            f"X has {_n_x} rows but y has {_n_y} rows — number of rows must match."
        )
    _n_classes = len(np.unique(_y_arr))
    if _n_classes < 2:
        raise ValueError(
            f"y has {_n_classes} unique class(es); at least 2 are required."
        )
    if max_components > _X_arr.shape[1]:
        raise ValueError(
            f"max_components={max_components} exceeds the number of features in X "
            f"({_X_arr.shape[1]}); reduce max_components."
        )
    if not np.all(np.isfinite(_X_arr)):
        raise ValueError("X contains NaN or Inf values.")
    encoder = OneHotEncoder(sparse_output=False)
    yd_arr = encoder.fit_transform(np.array(y).reshape(-1, 1))
    yd = pd.DataFrame(yd_arr)
    cv2 = RepeatedStratifiedKFold(
        n_splits=cv2_splits, n_repeats=n_repeats, random_state=random_state
    )
    cv1 = StratifiedKFold(n_splits=cv1_splits)

    n_lv_candidates = list(range(1, max_components))
    n_candidates = len(n_lv_candidates)

    rep_means: list[float] = []
    rep_stds: list[float] = []
    rep_lvs: list[int] = []

    outer_aurocs: list[float] = []
    outer_n_lv: list[int] = []
    fold_in_repeat = 0

    cv2_iter = tqdm(
        cv2.split(X, y),
        total=cv2_splits * n_repeats,
        desc="PLS-DA CV2",
        unit="fold",
        disable=not progress,
    )
    for rest, test in cv2_iter:
        X_rest = X.iloc[rest, :]
        X_test = X.iloc[test, :]
        y_rest = y.iloc[rest]
        yd_rest = yd.iloc[rest, :]
        yd_test = yd.iloc[test, :]

        # Inner CV: average AUROC across V folds for each candidate n_LV.
        inner_aurocs = np.zeros((cv1_splits, n_candidates))
        for v_idx, (train, validation) in enumerate(cv1.split(X_rest, y_rest)):
            X_train = X_rest.iloc[train, :]
            yd_train = yd_rest.iloc[train, :]
            X_val = X_rest.iloc[validation, :]
            yd_val = yd_rest.iloc[validation, :]
            args = zip(
                n_lv_candidates,
                repeat(X_train),
                repeat(yd_train),
                repeat(X_val),
                repeat(yd_val),
            )
            fold_aurocs: list[float]
            if n_jobs == 1:
                fold_aurocs = [_plsda_auroc(*a) for a in args]  # type: ignore[misc]
            else:
                n_workers = (
                    (multiprocessing.cpu_count() or 1) if n_jobs == -1 else max(1, n_jobs)
                )
                with multiprocessing.Pool(processes=n_workers) as pool:
                    fold_aurocs = pool.starmap(_plsda_auroc, args)  # type: ignore[arg-type]
            inner_aurocs[v_idx, :] = fold_aurocs

        mean_per_n_lv = inner_aurocs.mean(axis=0)
        # np.argmax returns the first occurrence on ties → parsimony.
        n_lv_star = int(np.argmax(mean_per_n_lv)) + 1

        outer_score = float(
            _plsda_auroc(n_lv_star, X_rest, yd_rest, X_test, yd_test)  # type: ignore[arg-type]
        )
        outer_aurocs.append(outer_score)
        outer_n_lv.append(n_lv_star)

        fold_in_repeat += 1
        if fold_in_repeat == cv2_splits:
            mean_auroc = float(np.mean(outer_aurocs))
            std_auroc = (
                float(np.std(outer_aurocs, ddof=1))
                if len(outer_aurocs) > 1
                else float("nan")
            )
            mode_lv = _modal_int_with_parsimony(outer_n_lv)
            rep_means.append(mean_auroc)
            rep_stds.append(std_auroc)
            rep_lvs.append(mode_lv)
            outer_aurocs = []
            outer_n_lv = []
            fold_in_repeat = 0

    model_table = pd.DataFrame(
        {
            "rep": list(range(1, n_repeats + 1)),
            "LV": rep_lvs,
            "AUROC": rep_means,
            "AUROC_std": rep_stds,
        }
    )
    model_table["rep"] = model_table["rep"].astype(int)
    model_table["LV"] = model_table["LV"].astype(int)
    model_table["AUROC"] = model_table["AUROC"].astype(float)
    model_table["AUROC_std"] = model_table["AUROC_std"].astype(float)

    best_models: list[PLSRegression] = []
    for lv in rep_lvs:
        model = PLSRegression(
            n_components=lv, scale=True, max_iter=1000
        ).fit(_X_arr, yd_arr)
        best_models.append(model)

    return {"models": best_models, "table": model_table}

motco.stats.pls.calculate_vips(model, components=None)

Estimates Variable Importance in Projection (VIP) in Partial Least Squares (PLS)

Parameters:

Name Type Description Default
model

model generated from the PLSRegression function

 required
components Union[None, list[int]]

if not None, a list of integers indicating the components to compute the VIPs from. If None, all components are taken into account. Default None.

 None

Returns:

Name Type Description
vips array

variable importance in projection for each variable
Source code in src/motco/stats/pls.py 
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
def calculate_vips(
    model,
    components: Union[None, list[int]] = None,
) -> np.ndarray:
    """
    Estimates Variable Importance in Projection (VIP)
    in Partial Least Squares (PLS)

    Parameters
    ----------
    model: PLSRegression
        model generated from the PLSRegression function
    components: Union[None, list[int]]
        if not None, a list of integers indicating the components to compute
        the VIPs from. If None, all components are taken into account.
        Default None.

    Returns
    -------
    vips: np.array
        variable importance in projection for each variable
    """
    if components is not None:
        t = model.x_scores_[:, components]
        w = model.x_weights_[:, components]
        q = model.y_loadings_[:, components]
        p, h = w.shape
    else:
        w = model.x_weights_
        t = model.x_scores_
        q = model.y_loadings_
        p, h = w.shape

    vips = np.zeros((p,))
    s = np.diag(np.matmul(np.matmul(np.matmul(t.T, t), q.T), q)).reshape(h, -1)
    total_s = np.sum(s)
    for i in range(p):
        weight = np.array([(w[i, j] / np.linalg.norm(w[:, j])) ** 2 for j in range(h)])
        vips[i] = np.sqrt(p * np.matmul(s.T, weight)[0] / total_s)
    return vips