Addressing the problem
Each subsection below runs one usual tool and the mixture on the same data and compares the results. The mixture has \(K\) components. The prediction of \(y\) at a value of \(x\) is the mean of \(y\) once \(x\) is fixed, written \(E[y \mid x]\).
## Mean of the first variable (y) when the second (x) is fixed at each value
## of xv: the component means of y, weighted by how likely each component is
## at that x.
cond_mean <- function(fit, xv) {
vapply(xv, function(xx) {
g <- gmm_conditionalise(fit, given = c(NA, xx))
sum(g@weights * vapply(g@means, function(m) m[1L], numeric(1L)))
}, numeric(1L))
}Regression: lm()
The usual tool is lm(), which fits one straight line.
The mixture route fits a mixture to the pairs \((y, x)\) and conditions on \(x\). With one component, the conditional
mean is exactly the least-squares line. With three components, the line
can bend.
set.seed(20260617)
n <- 400L
x <- runif(n, -3, 3)
y <- 0.3 * x + 1.2 * pmax(x, 0) + rnorm(n, sd = 0.4) # bends at x = 0
dat <- data.frame(y = y, x = x)
joint <- gmm_target_from_samples(cbind(y, x))
fit1 <- fit_proxymix(joint, N = 1L, regime = "moment", ridge_eps = 0)
fit3 <- fit_proxymix(joint, N = 3L, regime = "sample", max_iter = 150L)
## With one component, the slope of E[y | x] should equal the lm slope.
slope_mix <- gmm_conditionalise(fit1, given = c(NA, 1))@means[[1L]] -
gmm_conditionalise(fit1, given = c(NA, 0))@means[[1L]]
slope_lm <- unname(coef(lm(y ~ x, dat))["x"])
diff_reg <- abs(slope_mix - slope_lm)grid_reg <- data.frame(x = seq(-3, 3, length.out = 200L))
grid_reg$ols <- as.numeric(predict(lm(y ~ x, dat), newdata = grid_reg))
grid_reg$mix <- cond_mean(fit3, grid_reg$x)
ggplot2::ggplot() +
ggplot2::geom_point(data = dat, ggplot2::aes(x, y), colour = "grey60",
alpha = 0.4, size = 0.7) +
ggplot2::geom_line(data = grid_reg,
ggplot2::aes(x, ols, colour = "lm (K = 1)"),
linewidth = 0.9) +
ggplot2::geom_line(data = grid_reg,
ggplot2::aes(x, mix, colour = "mixture (K = 3)"),
linewidth = 0.9) +
ggplot2::scale_colour_manual(
name = NULL,
values = c("lm (K = 1)" = "#0072B2", "mixture (K = 3)" = "#D55E00")
) +
ggplot2::labs(
x = "x", y = "y",
title = "Regression: a straight line and a conditioned mixture"
) +
ggplot2::theme_minimal(base_size = 11) +
ggplot2::theme(legend.position = "top")
The least-squares line (one component) and the conditional mean of a three-component mixture, on data whose true relationship bends at zero. The mixture follows the bend, and the straight line does not.
Each component gives its own straight line for \(y\) given \(x\). The probability that a point at \(x\) belongs to component \(k\), written \(\pi_k(x)\), changes with \(x\). The prediction \(\sum_k \pi_k(x)\, \mu_k^{y \mid x}\) therefore switches smoothly from one line to another, which produces the bend. This method is known as Gaussian-mixture regression (de Veaux, 1989).
Kernel regression: Nadaraya-Watson
A kernel smoother predicts \(y\) at
\(x\) as a weighted average of the
observed \(y\) values, with weights
that fall off with distance from \(x\).
The Nadaraya-Watson estimator (Nadaraya, 1964; Watson, 1964) uses
normal-shaped weights, and the width of those weights, \(h\), is called the bandwidth. It is
available in R as stats::ksmooth() and in the
np package. The mixture route places one small normal
distribution on every data point, a kernel density estimate of \((y, x)\), and conditions on \(x\). The conditional mean is then exactly
the Nadaraya-Watson estimate.
h <- 0.4 # bandwidth
nw <- function(xq) {
vapply(xq, function(q) {
w <- dnorm(q, x, h) # Nadaraya-Watson weights
sum(w * y) / sum(w)
}, numeric(1L))
}
## One normal component per data point, then condition on x.
kde <- gmm(weights = rep(1 / n, n),
means = lapply(seq_len(n), function(i) c(y[i], x[i])),
covariances = rep(list(diag(c(h^2, h^2))), n))
xq <- seq(-2.5, 2.5, length.out = 21L)
diff_nw <- max(abs(nw(xq) - cond_mean(kde, xq)))grid_nw <- data.frame(x = seq(-3, 3, length.out = 200L))
grid_nw$ols <- as.numeric(predict(lm(y ~ x, dat), newdata = grid_nw))
grid_nw$nw <- nw(grid_nw$x)
ggplot2::ggplot() +
ggplot2::geom_point(data = dat, ggplot2::aes(x, y), colour = "grey60",
alpha = 0.4, size = 0.7) +
ggplot2::geom_line(data = grid_nw,
ggplot2::aes(x, ols, colour = "least squares (K = 1)"),
linewidth = 0.9) +
ggplot2::geom_line(data = grid_nw,
ggplot2::aes(x, nw, colour = "kernel (K = n)"),
linewidth = 0.9) +
ggplot2::scale_colour_manual(
name = NULL,
values = c("least squares (K = 1)" = "#0072B2",
"kernel (K = n)" = "#D55E00")
) +
ggplot2::labs(
x = "x", y = "y",
title = "From a straight line to a kernel smoother"
) +
ggplot2::theme_minimal(base_size = 11) +
ggplot2::theme(legend.position = "top")
The two ends of one scale: the least-squares line (one component) and the Nadaraya-Watson smoother (one component per data point). The same conditioning step produces both.
After conditioning, component \(i\)
has mean \(y_i\) and weight \(\pi_i(x) \propto \mathcal{N}(x \mid x_i,
h^2)\). The conditional mean \(\sum_i
\pi_i(x)\, y_i\) is the Nadaraya-Watson ratio, and the bandwidth
\(h\) is the spread of each component.
Least squares uses \(K = 1\) and the
kernel smoother uses \(K = n\).
from_kde() gives mixtures with any number of components in
between.
Clustering: kmeans()
kmeans() splits the rows into groups. The
mclust package fits a Gaussian mixture for the same purpose
(Fraley and Raftery, 2002). In the mixture route, each component is a
cluster. Each row gets a probability of belonging to each component,
called its responsibility, instead of a single label.
set.seed(20260617)
x_clust <- rbind(
mvnfast::rmvn(150L, c(-2, -1), 0.5 * diag(2)),
mvnfast::rmvn(150L, c(2, 0), matrix(c(0.6, 0.3, 0.3, 0.4), 2L)),
mvnfast::rmvn(150L, c(0, 2.5), 0.3 * diag(2))
)
colnames(x_clust) <- c("V1", "V2")
target_clust <- gmm_target_from_samples(x_clust)
fitc <- fit_proxymix(target_clust, N = 3L, regime = "sample",
max_iter = 150L)
## Responsibility of each component for each row.
responsibilities <- function(fit, xx) {
comp <- vapply(seq_len(gmm_n_components(fit)), function(k) {
fit@weights[k] *
mvnfast::dmvn(xx, mu = fit@means[[k]], sigma = fit@covariances[[k]])
}, numeric(nrow(xx)))
comp / rowSums(comp)
}
resp <- responsibilities(fitc, x_clust)
mean_confidence <- mean(apply(resp, 1L, max))| component 1 | component 2 | component 3 |
|---|---|---|
| 1 | 0 | 0 |
| 1 | 0 | 0 |
| 1 | 0 | 0 |
| 1 | 0 | 0 |
The responsibility of component \(k\) for a point \(x\) is \(\pi_k\,\mathcal{N}(x \mid \mu_k, \Sigma_k) / \sum_j \pi_j\,\mathcal{N}(x \mid \mu_j, \Sigma_j)\), where \(\pi_k\) is the weight, \(\mu_k\) the mean and \(\Sigma_k\) the covariance of component \(k\).
Principal components: prcomp()
prcomp() finds the directions in which the data vary
most. These directions are the eigenvectors of the covariance matrix
(Jolliffe, 2002). The mixture route fits one component, whose covariance
equals the sample covariance, and takes its eigenvectors. A single
normal distribution fitted this way is the probabilistic form of
principal components (Tipping and Bishop, 1999).
fit_pca <- fit_proxymix(target_clust, N = 1L, regime = "moment",
ridge_eps = 0)
ev <- eigen(fit_pca@covariances[[1L]])$vectors
pr <- prcomp(x_clust)$rotation
## Each direction may point either way, so compare absolute values.
diff_pca <- max(abs(abs(ev) - abs(unname(pr))))vals <- eigen(fit_pca@covariances[[1L]])$values
mu_pca <- unname(fit_pca@means[[1L]])
axes <- data.frame(
x = mu_pca[1L], y = mu_pca[2L],
xend = mu_pca[1L] + 2 * sqrt(vals) * ev[1L, ],
yend = mu_pca[2L] + 2 * sqrt(vals) * ev[2L, ]
)
pts <- data.frame(x_clust, cluster = factor(max.col(resp)))
ggplot2::ggplot() +
ggplot2::geom_point(data = pts,
ggplot2::aes(V1, V2, colour = cluster),
alpha = 0.6, size = 0.9) +
ggplot2::geom_segment(
data = axes, ggplot2::aes(x = x, y = y, xend = xend, yend = yend),
arrow = grid::arrow(length = grid::unit(0.2, "cm")), linewidth = 0.8
) +
ggplot2::scale_colour_viridis_d(name = "cluster", end = 0.85) +
ggplot2::coord_equal() +
ggplot2::labs(x = expression(x[1]), y = expression(x[2]),
title = "Clusters and principal directions") +
ggplot2::theme_minimal(base_size = 11)
Clusters (colour) from a three-component fit and principal directions (arrows) from a one-component fit to the same data.
With more than one component, the eigenvectors of each component’s covariance give the main directions of variation within that cluster.
Penalised regression: ridge
When predictors are strongly correlated, least-squares slopes become
unstable. Ridge regression (Hoerl and Kennard, 1970) shrinks the slopes
toward zero with a penalty of size \(\lambda\). It is available in
glmnet and MASS::lm.ridge(). In the mixture
route, ridge_eps = lambda adds \(\lambda\) to the variance of each variable
before conditioning.
lambda <- c(0, 0.5, 2, 8)
slope_pm <- vapply(lambda, function(lam) {
f <- fit_proxymix(joint, N = 1L, regime = "moment", ridge_eps = lam)
gmm_conditionalise(f, given = c(NA, 1))@means[[1L]] -
gmm_conditionalise(f, given = c(NA, 0))@means[[1L]]
}, numeric(1L))
slope_formula <- cov(x, y) / (var(x) + lambda)
diff_ridge <- max(abs(slope_pm - slope_formula))| Penalty (lambda) | proxymix slope | Ridge formula |
|---|---|---|
| 0.0 | 0.8968 | 0.8968 |
| 0.5 | 0.7680 | 0.7680 |
| 2.0 | 0.5367 | 0.5367 |
| 8.0 | 0.2434 | 0.2434 |
Conditioning gives the slope \(\beta = \Sigma_{xx}^{-1}\Sigma_{xy}\), where \(\Sigma_{xx}\) is the covariance matrix of the predictors and \(\Sigma_{xy}\) their covariance with \(y\). Adding \(\lambda I\) to \(\Sigma_{xx}\) gives \((\Sigma_{xx} + \lambda I)^{-1}\Sigma_{xy}\), which is the ridge estimate.
Treatment effects: a separate model for each arm
The treatment effect at \(x\) is the
difference in the mean outcome between treated and untreated units with
that value of \(x\). It is called the
conditional average treatment effect. A common way to estimate it is the
two-model learner, or T-learner: fit one regression to the treated
units, another to the untreated units, and subtract. Causal forests
(Wager and Athey, 2018; package grf) and double machine
learning (package DoubleML) are other options. The mixture
route fits one mixture over the outcome, the treatment and the covariate
with fit_uplift(). Every later step uses this one fit
without refitting.
The data below are simulated with a known effect, \(\tau(x) = 0.5 + x\), so the estimate can be compared with the truth. The treatment is assigned at random with probability 0.5.
set.seed(20260902)
n_up <- 600L
x_up <- rnorm(n_up)
t_up <- rbinom(n_up, 1L, 0.5)
tau_true <- function(v) 0.5 + v
y_up <- 1 + tau_true(x_up) * t_up + rnorm(n_up, sd = 0.5)
dat_up <- data.frame(y = y_up, t = t_up, x = x_up)
model <- fit_uplift(dat_up, "y", "t", "x", N = 2L, regime = "sample",
max_iter = 80L, seed = 1L)
model
#> <uplift_model>: K = 2 regimes, assume = "ignorability"
#> outcome : y (continuous)
#> treatment : t (arms 0 vs 1)
#> covariates: x
#> trained on: 600 unitsIn this printed output, and in the report below, “regimes” means the mixture’s components.
proxy_cate() returns the estimated effect with a
standard error and a 95 per cent interval for each unit supplied. The
check below compares it with the T-learner built from lm().
In the model y ~ t * x, the coefficients of t
and t:x are the differences in intercept and slope between
the two arms’ least-squares lines.
grid_up <- data.frame(x = seq(-2, 2, length.out = 41L))
cate <- proxy_cate(model, grid_up)
err_cate <- max(abs(cate$tau - tau_true(grid_up$x)))
## T-learner: one least-squares line per arm, fitted as one model.
fit_arms <- lm(y ~ t * x, data = dat_up)
b_arms <- coef(fit_arms)[c("t", "t:x")]
diff_arms <- max(abs(cate$tau - (b_arms[[1L]] + b_arms[[2L]] * grid_up$x)))
a_arms <- cbind(0, 1, 0, grid_up$x) # picks out t + x * t:x
se_arms <- sqrt(rowSums((a_arms %*% vcov(fit_arms)) * a_arms))
se_ratio <- range(cate$se / se_arms)
## Distance of the fitted intercept and slope from 0.5 and 1, in
## standard errors.
z_arms <- (b_arms - c(0.5, 1)) / sqrt(diag(vcov(fit_arms)))[c("t", "t:x")]
## Treatment value at the centre of each mixture component.
arm_of_component <- vapply(model@fit@means, function(m) {
m[model@roles$treatment]
}, numeric(1L))cate_df <- as.data.frame(cate)
cate_df$x <- grid_up$x
cate_df$truth <- tau_true(grid_up$x)
ggplot2::ggplot(cate_df, ggplot2::aes(x)) +
ggplot2::geom_ribbon(ggplot2::aes(ymin = ci_lo, ymax = ci_hi),
fill = "#56B4E9", alpha = 0.35) +
ggplot2::geom_line(ggplot2::aes(y = tau, colour = "proxymix estimate"),
linewidth = 0.9) +
ggplot2::geom_line(ggplot2::aes(y = truth, colour = "true effect"),
linewidth = 0.9, linetype = "dashed") +
ggplot2::geom_hline(yintercept = 0, colour = "grey60",
linewidth = 0.3) +
ggplot2::scale_colour_manual(
name = NULL,
values = c("proxymix estimate" = "#0072B2",
"true effect" = "#D55E00")
) +
ggplot2::labs(
x = "covariate x", y = "treatment effect on y",
title = "Treatment effect from one mixture fit"
) +
ggplot2::theme_minimal(base_size = 11) +
ggplot2::theme(legend.position = "top")
The estimated treatment effect with its 95 per cent interval, and the true effect 0.5 + x. The estimate is a straight line, the difference between the two arms’ least-squares lines.
proxy_decide() recommends treatment where the estimated
effect, multiplied by the value of one unit of outcome, exceeds the cost
of treating.
decision <- proxy_decide(model, grid_up, value = 1, cost = 0.5)
switch_x <- grid_up$x[min(which(decision$action == 1L))]| Covariate x | Estimated effect | True effect | Recommended arm | Net value |
|---|---|---|---|---|
| -2 | -1.665 | -1.5 | 0 | -2.165 |
| -1 | -0.615 | -0.5 | 0 | -1.115 |
| 0 | 0.435 | 0.5 | 0 | -0.065 |
| 1 | 1.485 | 1.5 | 1 | 0.985 |
| 2 | 2.534 | 2.5 | 1 | 2.034 |
proxy_cate() accepts only the treatment values seen in
the data, here 0 and 1, and refuses others, as shown below.
refusal <- tryCatch(
proxy_cate(model, grid_up, t1 = 100, t0 = 0),
error = function(e) e
)
msg_lines <- strsplit(conditionMessage(refusal), "\n")[[1L]]
writeLines(strwrap(msg_lines, width = 70L, exdent = 2L))
#> `t1`/`t0` do not match the treatment levels observed at fit time.
#> ℹ Observed treatment levels: 0 and 1.
#> ℹ Supplied: t1 = 100, t0 = 0.
#> ℹ Omit `t1`/`t0` to use the observed arms, or pass one of the
#> observed levels.The identification report states what the estimate is, the assumption it rests on, how well the treated and untreated units overlap in \(x\), and how far the estimate could move if a hidden group affected both treatment and outcome. It also lists what is not identified, meaning what the data cannot determine however many units are observed. An example is the counterfactual outcome of one unit: its outcome under the treatment it did not receive.
proxy_identification_report(model, grid_up)
#> Warning: The treatment is almost constant within every regime, so the within-regime
#> effect is zero and the gap equals the whole estimated effect.
#> ℹ The gap carries no information about confounding for this model.
#> == Identification report ==================================
#> Estimand : CATE / uplift: E[Y | do(T=t1), X] - E[Y | do(T=t0), X]
#> Assumption : ignorability
#> requires (Y(0), Y(1)) independent of T given X.
#> Regimes : K = 2 Outcome scale: continuous
#> Units : 41
#> Overlap : 100.0% of units adequately supported
#> Confounding gap (value at risk if a latent regime confounds):
#> mean |Delta| = 1.12, max |Delta| = 2.534
#> NOT identified : the individual counterfactual law
#> (its variance and tail probabilities).
#> ===========================================================The estimated effect is the difference of two conditional means of the same mixture, \(E[y \mid t = 1, x]\) minus \(E[y \mid t = 0, x]\). It is the same conditioning step as in the regression above. Reading this difference as a causal effect requires ignorability: given \(x\), the treatment must be unrelated to anything else that affects the outcome. Random assignment, as here, satisfies it.
The report carries a warning because, on data like these, the confounding gap equals the whole estimated effect. Its largest value, 2.534, is the estimated effect at \(x = 2\). Each component holds one arm, so a hidden group that matched the arms could account for the entire difference between them. The gap is therefore not a separate check.
Other uses
Missing values can be filled in with
gmm_conditionalise() on the observed variables, either with
the conditional mean or with a random draw from rgmm().
Density estimation is the basic use of a mixture.
from_kde() compresses a kernel density estimate into a
mixture with few components. gmm_observe() performs the
Kalman filter update for a mixture.