PCA and PLS charts: synthetic examples and held-out prediction

First inspect the geometry

The data lie on the max-normalized simplex: the surface of points with nonnegative coordinates of which the largest equals one. Distance on this surface is the length of the shortest path that stays on it. We place points on two adjacent faces. A latent coordinate t runs from -0.2 to 0.2 and is mapped to (1, 1+t) on one face and to (1-t, 1) on the other, and the remaining coordinates stay below one.

t.line <- c(seq(-.2, -.005, length.out=30), seq(.005, .2, length.out=30))
seam <- cbind(ifelse(t.line < 0, 1, 1-t.line),
              ifelse(t.line < 0, 1+t.line, 1), .2)
flat <- unfold.simplex.chart(seam, c(1, .95, .2))
stopifnot(flat$certified)
old.par <- par(mfrow=c(1,2), mar=c(4,4,2,1))
plot(seam[,1:2], asp=1, pch=16, col="#0072B2", xlim=c(.76,1.24),
     ylim=c(.76,1.24), xlab="First coordinate", ylab="Second coordinate",
     main="Ambient coordinates")
plot(flat$coordinates[,1:2], asp=1, pch=16, col="#009E73",
     xlim=c(.76,1.24), ylim=c(.76,1.24), xlab="First coordinate",
     ylab="Unfolded second coordinate", main="Unfolded coordinates")

A straight latent coordinate appears bent in ambient coordinates. Unfolding rotates the adjacent face into the query face; both axes use the same physical scale.

par(old.par)
variance.one <- function(x) {
    values <- svd(scale(x, center=TRUE, scale=FALSE), nu=0, nv=0)$d^2
    values[1]/sum(values)
}
knitr::kable(data.frame(Chart=c("Ambient", "Unfolded"),
    `Variance retained by one PC (%)`=100*c(variance.one(seam),
        variance.one(flat$coordinates)), check.names=FALSE), digits=1)

That the unfolded chart retains more variance is a statement about the geometry. It is not evidence that the outcome is predicted better; the examples below examine prediction separately. The charts with and without unfolding use the same neighbors and the same weights. Unfolding is refused near a corner, when two coordinates are tied for the maximum, or when a neighborhood extends over more than two faces. Such cases are counted and reported; they are not forced into a chart.

Define examples before comparing methods

Each observation has six coordinates. The seam coordinate t and a broad coordinate u3 vary substantially. A narrow coordinate u4 has a much smaller range; the remaining two coordinates are nuisance features. The table gives the true log odds. The true probability is plogis(eta) and the observed response is a Bernoulli draw with that probability.

The narrow coordinate varies only from 0.10 to 0.14, yet in the low-variance example its effect on the log odds is made as large as the effect of the broad coordinate in the high-variance example. This is a deliberately hard case for dimension reduction based on variance. It is not a claim that such signals are common in real data. The predictors are not standardized before the directions are computed, because standardizing would change the geometry and remove part of the contrast under study. The component scores are scaled within the local logistic fit.

scenarios <- c("Seam", "High variance", "Low variance", "Constant", "Nonlinear")
make.example <- function(seed, scenario, n=360L) {
    set.seed(seed)
    t <- runif(n, -.2, .2)
    X <- cbind(ifelse(t < 0, 1, 1-t), ifelse(t < 0, 1+t, 1),
               runif(n, .1, .5), runif(n, .10, .14),
               runif(n, .02, .06), runif(n, .02, .06))
    eta <- switch(scenario,
        "Seam"=-1+2*t/.2,
        "High variance"=-1+2*(X[,3]-.3)/.2,
        "Low variance"=-1+2*(X[,4]-.12)/.02,
        "Constant"=rep(-1,n),
        "Nonlinear"=-1+2.5*((t/.2)^2-1/3))
    probability <- plogis(eta)
    list(X=X, y=rbinom(n,1,probability), probability=probability, t=t)
}

Fit on training observations; choose dimension on tuning observations

Each of five independent repetitions per example has 180 training, 60 tuning and 120 test observations. The observations are independent in this simulation. For repeated specimens or matched sets, splitting must respect those groups instead.

We compare three charts: PCA in the original coordinates without unfolding (ambient PCA), PCA after unfolding, and PLS after unfolding. All three use 60 training neighbors, Gaussian weights, a ridge penalty of 1 on the slopes, and pseudo-counts of half a success and half a failure for the intercept. The neighborhood size and the ridge penalty are held fixed so that differences can be attributed to the charts; no claim is made that these values are optimal.

For each chart we report two results: one with a chart of dimension one, and one with the dimension chosen among 1, 2 and 3 as the one with the lowest Brier loss on the tuning observations (the smaller dimension wins a tie). The predictions for the tuning observations are made from the training observations only. The predictions for the test observations come from the same fits; the tuning observations are not added to the training set. Neither the test outcomes nor the true probabilities play any part in choosing a model.

methods <- c("Ambient PCA", "Unfolded PCA", "Unfolded PLS")
fit.chart <- function(X, y, query, method, dimension, degree=1L) {
    fit.local.likelihood(X, y, X.eval=query, likelihood.family="bernoulli",
        chart.geometry="linf.simplex", unfold.faces=method != "Ambient PCA",
        coordinate.method=if (method == "Unfolded PLS") "local.pls" else "local.pca",
        support.size=60L, chart.dim=dimension, degree=degree,
        kernel="gaussian", lambda.ridge=1, intercept.mass=.5,
        max.iter=70L, tol=1e-7, fallback="degree0")
}
train <- 1:180
tune <- 181:240
test <- 241:360

Here is a complete fit for the low-variance example. Replacing "Unfolded PLS" with "Unfolded PCA" changes the directions while retaining the same geometry and likelihood settings.

dat <- make.example(271301, "Low variance")
fit <- fit.chart(dat$X[train,], dat$y[train], dat$X[test,],
                 "Unfolded PLS", dimension=1L)
head(fit$fitted.values)
fit$diagnostics$status.counts

The loop below carries out the comparison. It includes two baselines. The local constant uses the same neighbors and weights as the charts but no directions, and the global constant is the smoothed proportion of successes in the whole training set. Both are useful for reference, especially in the example whose true probability is constant.

The results shown are from all 25 repetitions (five examples, five repetitions each). To keep the vignette quick to build, they are read from a file stored with the package instead of being recomputed. The last section explains how to rerun the loop with the same seeds, sample sizes and settings. The geometry demonstration and the probability slices are computed each time the vignette is built.

rows <- list()
illustrations <- list()
append.result <- function(scenario, repetition, method, mode, dimension,
                          prediction, truth, response, fallbacks, tuning.fallbacks=NA_integer_) {
    stopifnot(all(is.finite(prediction)), all(prediction > 0 & prediction < 1))
    data.frame(scenario=scenario, repetition=repetition, method=method,
        mode=mode, dimension=dimension,
        probability.mse=mean((prediction-truth)^2),
        brier=mean((prediction-response)^2),
        log.loss=-mean(response*log(prediction)+(1-response)*log1p(-prediction)),
        fallbacks=fallbacks, tuning.fallbacks=tuning.fallbacks)
}
for (s in seq_along(scenarios)) for (repetition in 1:5) {
    scenario <- scenarios[s]
    dat <- make.example(271000+100*s+repetition, scenario)
    example.predictions <- list()
    for (method in methods) {
        fits <- lapply(1:3, function(d) fit.chart(dat$X[train,], dat$y[train],
            dat$X[c(tune,test),], method, d))
        tuning.loss <- vapply(fits, function(f)
            mean((f$fitted.values[1:60]-dat$y[tune])^2), numeric(1))
        best <- which.min(tuning.loss)
        for (mode in c("One dimension", "Tuned dimension")) {
            d <- if (mode == "One dimension") 1L else best
            chosen <- fits[[d]]
            prediction <- chosen$fitted.values[61:180]
            rows[[length(rows)+1L]] <- append.result(scenario, repetition, method,
                mode, d, prediction, dat$probability[test], dat$y[test],
                sum(chosen$diagnostics$per.eval$fallback[61:180]),
                sum(chosen$diagnostics$per.eval$fallback[1:60]))
            if (repetition == 1L) example.predictions[[paste(method,mode)]] <- prediction
        }
    }
    local <- fit.chart(dat$X[train,], dat$y[train], dat$X[test,],
                       "Ambient PCA", dimension=1L, degree=0L)
    constants <- list("Local constant"=local$fitted.values,
        "Global constant"=rep((sum(dat$y[train])+.5)/(length(train)+1), length(test)))
    for (method in names(constants)) {
        rows[[length(rows)+1L]] <- append.result(scenario, repetition, method,
            "Constant", 0L, constants[[method]], dat$probability[test], dat$y[test], 0L)
        if (repetition == 1L) example.predictions[[method]] <- constants[[method]]
    }
    if (repetition == 1L) illustrations[[scenario]] <- list(data=dat,
        predictions=example.predictions)
}
results <- do.call(rbind, rows)
stopifnot(nrow(results) == 200L)

Measure estimation and prediction separately

Two measures are reported, and lower is better for both. The probability mean squared error (MSE) compares the predicted probability with the true probability, which is known here, and so measures the accuracy of estimation directly. The Brier loss compares the predicted probability with the observed binary response, and therefore also contains the variation of the response itself, which no method can remove. The log loss on the test observations is stored in results as well. A low MSE, or an optimizer that converged, does not show that the predicted probabilities are well calibrated.

The table gives means over the five repetitions. The figures show each repetition as well as the mean, so that a favorable average cannot conceal repetitions in which the order of the methods is reversed. Five repetitions are enough for an illustration. They do not make a simulation study with adequate power, and they do not rank the methods in general.

summary.table <- aggregate(cbind(probability.mse,brier) ~ scenario+method+mode,
                           results, mean)
summary.table <- summary.table[order(match(summary.table$scenario, scenarios),
    match(summary.table$mode,c("One dimension","Tuned dimension","Constant")),
    match(summary.table$method,c(methods,"Local constant","Global constant"))),]
knitr::kable(summary.table, digits=4, row.names=FALSE,
    col.names=c("Example","Method","Chart","Probability MSE","Brier loss"))

colors <- c("Ambient PCA"="#0072B2", "Unfolded PCA"="#009E73",
    "Unfolded PLS"="#D55E00", "Local constant"="#666666", "Global constant"="#AA4499")
labels <- c("PCA", "Unfolded\nPCA", "Unfolded\nPLS", "Local\nconstant", "Global\nconstant")
old.par <- par(mfrow=c(5,2), mar=c(3.2,4,2,1), oma=c(0,0,1,0), cex=.8)
for (scenario in scenarios) for (mode in c("One dimension","Tuned dimension")) {
    tab <- results[results$scenario==scenario & results$mode %in% c(mode,"Constant"),]
    order.methods <- names(colors)
    plot(NA, xlim=c(.7,5.3), ylim=c(0, max(tab$probability.mse)*1.08+1e-6),
         xaxt="n", xlab="", ylab="Probability MSE", main=paste(scenario,mode,sep=" — "))
    axis(1, at=1:5, labels=labels, cex.axis=.8)
    for (j in 1:5) {
        values <- tab$probability.mse[tab$method==order.methods[j]]
        points(j+seq(-.12,.12,length.out=5), values, pch=16, cex=.7,
               col=adjustcolor(colors[j],alpha.f=.6))
        points(j, mean(values), pch=18, cex=1.3, col=colors[j])
    }
}

Held-out probability error for all five repetitions. Small points are independent repetitions; diamonds are their means. Each panel uses its own labelled vertical scale. Lower is better. The constant baselines are identical in the two columns.

par(old.par)

Inspect the first repetition, without choosing favorable seeds

The next plot shows predicted against true probabilities for the test observations of the first repetition of each example. The first repetition was designated in advance and was not chosen for its results. Each point is one test observation. On the diagonal, the prediction equals the true probability; the diagonal is not a fitted line. In the constant example all observations have the same true probability, so their vertical spread shows how much each method follows the noise.

old.par <- par(mfrow=c(5,2), mar=c(3,3.3,2,1), cex=.8)
for (scenario in scenarios) for (method in c("Unfolded PCA","Unfolded PLS")) {
    item <- illustrations[[scenario]]
    prediction <- item$predictions[[paste(method,"One dimension")]]
    plot(item$data$probability[test], prediction, xlim=c(0,1), ylim=c(0,1),
         xlab="True probability", ylab="Predicted probability", pch=16, cex=.65,
         col=adjustcolor(colors[method],alpha.f=.55), main=paste(scenario,method,sep=" — "))
    abline(0,1,lty=2,col="grey35")
}

Predicted versus known true probabilities for first-repetition test observations. Columns compare unfolded PCA and PLS, each with one dimension. Points near the diagonal recover the true probabilities; all panels share the same axes.

par(old.par)

View the estimated probability functions

The curves below are fitted to the training observations of the first repetition. In each panel the coordinate that carries the signal varies over its range while the other coordinates are fixed at the centers of their ranges. The curves are the fitted probabilities at those points; they are not averages of the points in the scatterplots above. The black curve is the true probability. All methods use charts of dimension one, and nothing was tuned again for these plots.

slices <- list()
old.par <- par(mfrow=c(3,2), mar=c(4,4,2,1), cex=.85)
for (scenario in scenarios) {
    grid <- seq(-1,1,length.out=60)
    t <- if (scenario %in% c("Seam","Nonlinear","Constant")) .2*grid else rep(.08,60)
    X.slice <- cbind(ifelse(t < 0,1,1-t), ifelse(t < 0,1+t,1), .3,.12,.04,.04)
    if (scenario=="High variance") X.slice[,3] <- .3+.2*grid
    if (scenario=="Low variance") X.slice[,4] <- .12+.02*grid
    eta <- if (scenario=="Constant") rep(-1,60) else if (scenario=="Nonlinear")
        -1+2.5*(grid^2-1/3) else -1+2*grid
    truth <- plogis(eta)
    dat <- illustrations[[scenario]]$data
    prediction <- vapply(methods, function(method)
        fit.chart(dat$X[train,],dat$y[train],X.slice,method,1L)$fitted.values,
        numeric(60))
    axis.value <- if (scenario=="High variance") X.slice[,3] else if (scenario=="Low variance")
        X.slice[,4] else t
    axis.label <- if (scenario=="High variance") "Broad coordinate u3" else if (scenario=="Low variance")
        "Narrow coordinate u4" else "Seam coordinate t"
    plot(axis.value,truth,type="l",lwd=2,col="black",ylim=c(0,1),
         xlab=axis.label,ylab="Probability",main=scenario)
    for (j in seq_along(methods)) lines(axis.value,prediction[,j],col=colors[j],lwd=1.5,lty=j)
    slices[[scenario]] <- list(X=X.slice,truth=truth,prediction=prediction)
}
plot.new()
legend("center",c("Truth",methods),col=c("black",colors[methods]),
       lty=c(1,1,2,3),lwd=c(2,1.5,1.5,1.5),bty="n",cex=.9)

Conditional probability slices through the same first-repetition training fits. Nuisance coordinates are fixed at their generating centers. Black is the known truth; colored curves use one-dimensional charts. All probability axes run from zero to one.

par(old.par)

Check selected dimensions and fallbacks

tuned <- results[results$mode=="Tuned dimension",]
selection <- as.data.frame(xtabs(~scenario+method+dimension, tuned))
selection <- selection[selection$Freq>0,]
knitr::kable(selection, row.names=FALSE,
    col.names=c("Example","Method","Selected dimension","Repetitions"))

knitr::kable(aggregate(cbind(fallbacks,tuning.fallbacks) ~ scenario+method+mode,
                      results, sum), row.names=FALSE,
    col.names=c("Example","Method","Chart","Test fallbacks","Tuning fallbacks"))

A fallback occurs when a chart cannot be certified or the optimizer fails; the prediction is then the smoothed local proportion. These predictions are kept in the scores above and are not removed to make a method look better. The table counts fallbacks only among the fits that are reported, namely those with dimension one and those with the selected dimension. Some fits are counted under both headings, so the two counts should not be added. results also contains the selected dimension and the test log loss of each fit.

Reading the comparison

The low-variance example shows that using the response to choose directions is not sufficient in itself. With one dimension, unfolded PCA and PLS have mean probability MSEs of 0.0472 and 0.0473, respectively. After selecting dimension on the tuning observations, these become 0.0442 and 0.0162. Weighted PLS is based on covariances, so it remains sensitive to the scale of the predictors, to the noise in a binary response, and to the number of directions kept. In this example its advantage for the narrow signal appears only after the dimension has been selected; its first direction alone does not provide it.

In the constant example the global constant has the lowest probability MSE. In the high-variance example PCA with a tuned dimension has a lower mean probability MSE than PLS with a tuned dimension. Tuning the dimension sometimes makes the test results worse than the fit of dimension one, because 60 tuning observations give a noisy criterion for the choice. This is why the vignette shows the baselines, both the fixed and the selected dimension, and every repetition.

Each example was designed to isolate one issue. The seam example separates the correction of a fold in the geometry from the learning of a direction related to the response. The high-variance example suits the criterion of PCA. The low-variance example puts the information about the response in a direction that PCA is likely to discard; a chart of higher dimension lets PCA recover it, at the price of more variable estimates. PLS can help at low dimension, provided its directions are computed from the training data alone.

The constant example is a necessary control, since it contains no direction to be found. The nonlinear example tests whether a nonlinear function is close enough to linear within a neighborhood. Where the signal is symmetric about a point, the local covariance on which PLS relies can be weak even though the response depends strongly on position. None of the methods is expected to be best throughout, and the variation among repetitions matters as much as the averages.

Before relying on such a comparison in a scientific analysis, extend the simulation to realistic sample sizes, class balance, dependence among observations, missing data and geometry. Tune the neighborhood size and the penalties within the training data, assess the calibration of the probabilities on a separate hold-out sample, and examine the results near corners of the simplex and where the data are sparse. Under case-control sampling, the probabilities estimated in the selected sample do not give the prevalence in the population unless the sampling fractions are known.

Reproduce and adapt

The stored results were produced by the code shown above. They contain all 25 seeds, all 200 rows of results, and the predictions of the first repetition of each example, and nothing else: no real data and no reference predictions from other software. They also record the package versions and the SHA-256 checksum of this vignette. When the vignette is built, it checks that the stored results were produced from the present version of its source.

To recompute every comparison while rendering, use:

options(geosmooth.recompute.vignettes=TRUE)
rmarkdown::render("vignettes/pca-pls-chart-comparison.Rmd")

After changing this vignette, regenerate the stored results by running Rscript vignettes/precompute-pca-pls.R from the package source directory. The computation runs serially and needs only single-threaded BLAS and OpenMP. The package versions used for the stored results are listed below. ?fit.local.likelihood and ?unfold.simplex.chart describe the fitting options and the diagnostics.

#> geosmooth   dgraphs      Rcpp    digest 
#>   "0.2.0"   "0.4.1"   "1.1.2"  "0.6.39"