p-Wasserstein Linear Projections
Description
![[Experimental]](data:image/svg+xml;base64,PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHhtbG5zOnhsaW5rPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5L3hsaW5rIiB3aWR0aD0iMTM4IiBoZWlnaHQ9IjIwIiByb2xlPSJpbWciIGFyaWEtbGFiZWw9ImxpZmVjeWNsZTogZXhwZXJpbWVudGFsIj4KICAgIDx0aXRsZT5saWZlY3ljbGU6IGV4cGVyaW1lbnRhbDwvdGl0bGU+CiAgICA8bGluZWFyR3JhZGllbnQgaWQ9InMiIHgyPSIwIiB5Mj0iMTAwJSI+CiAgICAgICAgPHN0b3Agb2Zmc2V0PSIwIiBzdG9wLWNvbG9yPSIjYmJiIiBzdG9wLW9wYWNpdHk9Ii4xIiAvPgogICAgICAgIDxzdG9wIG9mZnNldD0iMSIgc3RvcC1vcGFjaXR5PSIuMSIgLz4KICAgIDwvbGluZWFyR3JhZGllbnQ+CiAgICA8Y2xpcFBhdGggaWQ9InIiPgogICAgICAgIDxyZWN0IHdpZHRoPSIxMzgiIGhlaWdodD0iMjAiIHJ4PSIzIiBmaWxsPSIjZmZmIiAvPgogICAgPC9jbGlwUGF0aD4KICAgIDxnIGNsaXAtcGF0aD0idXJsKCNyKSI+CiAgICAgICAgPHJlY3Qgd2lkdGg9IjU1IiBoZWlnaHQ9IjIwIiBmaWxsPSIjNTU1IiAvPgogICAgICAgIDxyZWN0IHg9IjU1IiB3aWR0aD0iODMiIGhlaWdodD0iMjAiIGZpbGw9IiNmZTdkMzciIC8+CiAgICAgICAgPHJlY3Qgd2lkdGg9IjEzOCIgaGVpZ2h0PSIyMCIgZmlsbD0idXJsKCNzKSIgLz4KICAgIDwvZz4KICAgIDxnIGZpbGw9IiNmZmYiIHRleHQtYW5jaG9yPSJtaWRkbGUiIGZvbnQtZmFtaWx5PSJWZXJkYW5hLEdlbmV2YSxEZWphVnUgU2FucyxzYW5zLXNlcmlmIiB0ZXh0LXJlbmRlcmluZz0iZ2VvbWV0cmljUHJlY2lzaW9uIiBmb250LXNpemU9IjExMCI+CiAgICAgICAgPHRleHQgYXJpYS1oaWRkZW49InRydWUiIHg9IjI4NSIgeT0iMTUwIiBmaWxsPSIjMDEwMTAxIiBmaWxsLW9wYWNpdHk9Ii4zIiB0cmFuc2Zvcm09InNjYWxlKC4xKSIgdGV4dExlbmd0aD0iNDUwIj5saWZlY3ljbGU8L3RleHQ+CiAgICAgICAgPHRleHQgeD0iMjg1IiB5PSIxNDAiIHRyYW5zZm9ybT0ic2NhbGUoLjEpIiBmaWxsPSIjZmZmIiB0ZXh0TGVuZ3RoPSI0NTAiPmxpZmVjeWNsZTwvdGV4dD4KICAgICAgICA8dGV4dCBhcmlhLWhpZGRlbj0idHJ1ZSIgeD0iOTU1IiB5PSIxNTAiIGZpbGw9IiMwMTAxMDEiIGZpbGwtb3BhY2l0eT0iLjMiIHRyYW5zZm9ybT0ic2NhbGUoLjEpIiB0ZXh0TGVuZ3RoPSI3MzAiPmV4cGVyaW1lbnRhbDwvdGV4dD4KICAgICAgICA8dGV4dCB4PSI5NTUiIHk9IjE0MCIgdHJhbnNmb3JtPSJzY2FsZSguMSkiIGZpbGw9IiNmZmYiIHRleHRMZW5ndGg9IjczMCI+ZXhwZXJpbWVudGFsPC90ZXh0PgogICAgPC9nPgo8L3N2Zz4K)
This function will calculate linear projections from a set of predictions into the space of the covariates in terms of the p-Wasserstein distance.
Usage
WpProj(
X,
eta = NULL,
theta = NULL,
power = 2,
method = c("L1", "binary program", "stepwise", "simulated annealing", "L0"),
solver = c("lasso", "ecos", "lpsolve", "mosek"),
options = NULL
)
WpProj(
X,
eta = NULL,
theta = NULL,
power = 2,
method = c("L1", "binary program", "stepwise", "simulated annealing", "L0"),
solver = c("lasso", "ecos", "lpsolve", "mosek"),
options = NULL
)
Arguments
X |
An n×p matrix of covariates
|
eta |
An n×s matrix of predictions from a model
|
theta |
An optional An p×s parameter matrix for selection methods. Only makes sense if the original model is a linear model.
|
power |
The power of the Wasserstein distance to use. Must be >= 1.0. Will default to 2.0.
|
method |
The algorithm to calculate the Wasserstein projections. One of "L1", "binary program", "IP", "stepwise","simulated annealing", or "L0". Will default to "L1" if not provided. See details for more information.
|
solver |
Which solver to use? One of "lasso", "ecos", "lpsolve", or "mosek". See details for more information
|
options |
Options passed to the particular method and desired solver. See details for more information.
|
Details
Methods
The WpProj function is a wrapper for the various Wasserstein projection methods. It is designed to be a one-stop shop for all Wasserstein projection methods. It will automatically choose the correct method and solver based on the arguments provided. It will also return a standardized output for all methods. Each method has its own set of options that can be passed to it. See the documentation for each method for more information.
For the L1 methods, see L1_method_options() for more information. For the binary program methods, see binary_program_method_options() for more information. For the stepwise methods, see stepwise_method_options() for more information. For the simulated annealing methods, see simulated_annealing_method_options() for more information.
In most cases, we recommend using the L1 methods or binary program methods. The L1 methods are the fastest and applicable to Wasserstein powers of any value greater than 1 and function as direct linear projections into the space of the covariates. The binary program methods instead preserve the coefficients of the original model if this is of interest, such as when the original model was already a linear model. The binary program will instead function as a way of turning on and off certain coefficients in a way that minimizes the Wasserstein distance between reduced and original models. Of note, we also have available an approximate binary program method using a lasso solver. This method is faster than the exact binary program method but is not guaranteed to find the optimal solution. It is recommended to use the exact binary program method if possible. See binary_program_method_options() for more information on how to set up the approximate method as some arguments for the lasso solver should be specified. For more information on how this works, please also see the referenced paper.
The stepwise, simulated annealing, and L0 methods also select covariates like the binary program methods but they can be slower. They are presented merely for comparison purposes given they were used in the original paper.
Wasserstein distances and powers
The Wasserstein distance is a measure of distance between two probability distributions. It is defined as:
Wp(μ,ν)=(infπ∈Π(μ,ν)∫Rd×Rd∥x−y∥pdπ(x,y))1/p,
where Π(μ,ν) is the set of all joint distributions with marginals μ and ν. The Wasserstein distance is a generalization of the Euclidean distance, which is the case when p=2. In our function we have argument power that corresponds to the p of the equation above. The default power is 2.0 but any value greater than or equal to 1.0 is allowed. For more information, see the references.
The particular implementation of the Wasserstein distance is as follows. If μ is the original prediction from the original model, then we seek to find a new prediction ν that minimizes the Wasserstein distance between the two: argminνWp(μ,ν).
Value
object of class WpProj, which is a list with the following slots:
-
call: The call to the function
-
theta: A list of the final parameter matrices for each returned model
-
fitted.values: A list of the fitted values for each returned model
-
power: The power of the Wasserstein distance used
-
method: The method used to calculate the Wasserstein projections
-
solver: The solver used to calculate the Wasserstein projections
-
niter: The number of iterations used to calculate the Wasserstein projections. Not all methods return a number of iterations so this may be NULL
-
nzero: The number of non zero coefficients in the final models
References
Dunipace, Eric and Lorenzo Trippa (2020) https://arxiv.org/abs/2012.09999.
Examples
if(rlang::is_installed("stats")) {
# note we don't generate believable data with real posteriors
# these examples are just to show how to use the function
n <- 32
p <- 10
s <- 21
# covariates and coefficients
x <- matrix( stats::rnorm( p * n ), nrow = n, ncol = p )
beta <- (1:10)/10
#outcome
y <- x %*% beta + stats::rnorm(n)
# fake posterior
post_beta <- matrix(beta, nrow=p, ncol=s) + stats::rnorm(p*s, 0, 0.1)
post_mu <- x %*% post_beta #posterior predictive distributions
# fit models
## L1 model
fit.p2 <- WpProj(X=x, eta=post_mu, power = 2.0,
method = "L1", #default
solver = "lasso" #default
)
## approximate binary program
fit.p2.bp <- WpProj(X=x, eta=post_mu, theta = post_beta, power = 2.0,
method = "binary program",
solver = "lasso" #default because approximate algorithm is faster
)
## compare performance by measuring distance from full model
dc <- distCompare(models = list("L1" = fit.p2, "BP" = fit.p2.bp))
if(rlang::is_installed(c("ggplot2","ggsci"))) {
plot(dc)
}
## compare performance by measuring the relative distance between a null model
## and the predictions of interest as a pseudo R^2
r2.expect <- WPR2(predictions = post_mu, projected_model = dc) # can have negative values
r2.null <- WPR2(projected_model = dc) # should be between 0 and 1
if(rlang::is_installed(c("ggplot2","ggsci"))) {
plot(r2.null)
}
## we can also examine how predictions change in the models for individual observations
if(rlang::is_installed(c("ggplot2","ggsci","ggridges"))) {
ridgePlot(fit.p2, index = 21, minCoef = 0, maxCoef = 10)
}
}
if(rlang::is_installed("stats")) {
n <- 32
p <- 10
s <- 21
x <- matrix( stats::rnorm( p * n ), nrow = n, ncol = p )
beta <- (1:10)/10
y <- x %*% beta + stats::rnorm(n)
post_beta <- matrix(beta, nrow=p, ncol=s) + stats::rnorm(p*s, 0, 0.1)
post_mu <- x %*% post_beta
fit.p2 <- WpProj(X=x, eta=post_mu, power = 2.0,
method = "L1",
solver = "lasso"
)
fit.p2.bp <- WpProj(X=x, eta=post_mu, theta = post_beta, power = 2.0,
method = "binary program",
solver = "lasso"
)
dc <- distCompare(models = list("L1" = fit.p2, "BP" = fit.p2.bp))
if(rlang::is_installed(c("ggplot2","ggsci"))) {
plot(dc)
}
r2.expect <- WPR2(predictions = post_mu, projected_model = dc)
r2.null <- WPR2(projected_model = dc)
if(rlang::is_installed(c("ggplot2","ggsci"))) {
plot(r2.null)
}
if(rlang::is_installed(c("ggplot2","ggsci","ggridges"))) {
ridgePlot(fit.p2, index = 21, minCoef = 0, maxCoef = 10)
}
}
,
Clemens Schmid [ctb] (<https://orcid.org/0000-0003-3448-5715>, ETA
progres bar is adapted from their code),
Espen Bernton [ctb] ('Hilbert Sort' adapted from their code),
Mathieu Gerber [ctb] ('Hilbert Sort' adapted from their code),
Pierre Jacob [ctb] ('Hilbert Sort' adapted from their code),
Bin Dai [ctb] (W2 projections adapted from their 'OEM' code),
Jared Huling [ctb] (<https://orcid.org/0000-0003-0670-4845>, W2
projections adapted from their 'OEM' code),
Yixuan Qiu [ctb] (W2 projections adapted from their 'OEM' code),
Dominic Schuhmacher [ctb] ('Shortsimplex 'optimal transport method
adapted from their code),
Nicolas Bonneel [ctb] ('Network Simplex' algorithm adapted from their
code)