3 Usage
The argument loss selects the estimator;
penalty is one of "grSCAD",
"grLasso" or "grMCP"; lambda
controls sparsity and is best chosen by cross-validation (see
below).
3.1 Regression
set.seed(1)
n <- 1000; p <- 10
B <- matrix(0, p, 2); B[1, 1] <- B[2, 2] <- 1
x <- matrix(rnorm(n * p), n, p)
y <- (x %*% B[, 1]) / (0.5 + (x %*% B[, 2] + 1)^2) + 0.2 * rnorm(n)
## penalized principal least-squares SVM (P^2LSM)
fit <- ppm(x, y, H = 10, C = 1, loss = "lssvm", penalty = "grSCAD", lambda = 0.01)
round(fit$evectors[, 1:2], 3)
#> [,1] [,2]
#> [1,] 1 0
#> [2,] 0 1
#> [3,] 0 0
#> [4,] 0 0
#> [5,] 0 0
#> [6,] 0 0
#> [7,] 0 0
#> [8,] 0 0
#> [9,] 0 0
#> [10,] 0 0
summary(fit, d = 2)
#> Penalized Principal Machine (P2M) summary
#> Loss: lssvm Penalty: grSCAD lambda = 0.01
#> Working dimension d = 2
#>
#> Estimated basis of the central subspace:
#> Dir1 Dir2
#> x1 1e+00 -1e-04
#> x2 1e-04 1e+00
#> x3 0e+00 0e+00
#> x4 0e+00 0e+00
#> x5 0e+00 0e+00
#> x6 0e+00 0e+00
#> x7 0e+00 0e+00
#> x8 0e+00 0e+00
#> x9 0e+00 0e+00
#> x10 0e+00 0e+00
#>
#> Selected variables (2 of 10): x1, x2The same interface fits any other estimator. For example, the penalized principal asymmetric least squares regression (\(\mathrm{P}^2\mathrm{AR}\)) examines the response at several asymmetry levels \(\tau_k\):
## penalized principal asymmetric least squares regression (P^2AR)
fit_ar <- ppm(x, y, H = 10, C = 1, loss = "asls", penalty = "grSCAD", lambda = 0.05)
round(fit_ar$evectors[, 1:2], 3)
#> [,1] [,2]
#> [1,] 1.00 -0.02
#> [2,] 0.02 1.00
#> [3,] 0.00 0.00
#> [4,] 0.00 0.00
#> [5,] 0.00 0.00
#> [6,] 0.00 0.00
#> [7,] 0.00 0.00
#> [8,] 0.00 0.00
#> [9,] 0.00 0.00
#> [10,] 0.00 0.003.2 Binary classification
For a binary response the weighted losses ("wlssvm",
"wlogit", "wl2svm", "wsvm") code
the classes internally as \(\{-1,
+1\}\), and \(\mathrm{P}^2\mathrm{AR}\) applies to both
response types. The example below generates a binary response whose
central subspace is again spanned by the first two predictors.
yb <- sign(x[, 1] + x[, 2]^3 / 3 + 0.2 * rnorm(n))
## penalized principal asymmetric least squares (P^2AR)
fit_arb <- ppm(x, yb, loss = "asls", penalty = "grSCAD", lambda = 0.2)
round(fit_arb$evectors[, 1:2], 3)
#> [,1] [,2]
#> [1,] 0.927 -0.376
#> [2,] 0.376 0.927
#> [3,] 0.000 0.000
#> [4,] 0.000 0.000
#> [5,] 0.000 0.000
#> [6,] 0.000 0.000
#> [7,] 0.000 0.000
#> [8,] 0.000 0.000
#> [9,] 0.000 0.000
#> [10,] 0.000 0.000
## penalized principal weighted logistic regression (P^2WLR)
fit_wlr <- ppm(x, yb, loss = "wlogit", penalty = "grSCAD", lambda = 0.005)
summary(fit_wlr, d = 2)
#> Penalized Principal Machine (P2M) summary
#> Loss: wlogit Penalty: grSCAD lambda = 0.005
#> Algorithm: iterative GCD (line search on), 7 iterations, converged
#> Working dimension d = 2
#>
#> Estimated basis of the central subspace:
#> Dir1 Dir2
#> x1 0.9450 -0.3271
#> x2 0.3271 0.9450
#> x3 0.0000 0.0000
#> x4 0.0000 0.0000
#> x5 0.0000 0.0000
#> x6 0.0000 0.0000
#> x7 0.0000 0.0000
#> x8 0.0000 0.0000
#> x9 0.0000 0.0000
#> x10 0.0000 0.0000
#>
#> Selected variables (2 of 10): x1, x23.3 Monitoring convergence with the line search
For the iterative GCD losses the fitted object records the number of
outer iterations (iter), the termination status
(status), and, when the line search is on, the objective
path (objective) and the number of halvings \(m_t\) at each iteration
(n.halving). For \(\mathrm{P}^2\mathrm{LR}\) the safeguard is
on by default:
## penalized principal logistic regression (P^2LR)
fit_lr <- ppm(x, y, loss = "logit", penalty = "grSCAD", lambda = 0.01)
print(fit_lr)
#> Penalized Principal Machine (P2M) for Sufficient Dimension Reduction
#>
#> Call:
#> ppm(x = x, y = y, loss = "logit", lambda = 0.01, penalty = "grSCAD")
#>
#> Loss: logit Penalty: grSCAD lambda = 0.01 gamma = 3.7 C = 1 H = 10
#> Data: n = 1000, p = 10 (continuous response)
#> Algorithm: iterative GCD, step-halving line search on (max.halving = 10; 0 of 7 updates damped)
#> Iterations: 7 (converged)
#> Leading eigenvalues of the working matrix M:
#> 0.00040482, 5.05e-06, 0, 0, 0
fit_lr$n.halving # m_t at each iteration (0 = plain update)
#> [1] 0 0 0 0 0 0 0
all(diff(fit_lr$objective) <= 0) # Q never increases
#> [1] TRUE
## the plain updates give the same fit whenever no step is damped
fit_lr0 <- ppm(x, y, loss = "logit", penalty = "grSCAD", lambda = 0.01,
line.search = FALSE)
all.equal(fit_lr$M, fit_lr0$M)
#> [1] TRUEFor \(\mathrm{P}^2\mathrm{AR}\) the plain updates are the default, and the safeguard can be requested explicitly. Here it damps an update, terminates earlier, and returns the same basis as the plain updates.
fit_ar_ls <- ppm(x, y, loss = "asls", penalty = "grSCAD", lambda = 0.05,
line.search = TRUE)
c(plain = fit_ar$iter, safeguarded = fit_ar_ls$iter)
#> plain safeguarded
#> 88 56
table(fit_ar_ls$n.halving) # how often each m_t was used
#>
#> 0 5
#> 55 1
round(fit_ar_ls$evectors[, 1:2], 3)
#> [,1] [,2]
#> [1,] 1.00 -0.02
#> [2,] 0.02 1.00
#> [3,] 0.00 0.00
#> [4,] 0.00 0.00
#> [5,] 0.00 0.00
#> [6,] 0.00 0.00
#> [7,] 0.00 0.00
#> [8,] 0.00 0.00
#> [9,] 0.00 0.00
#> [10,] 0.00 0.00
plot(fit_ar_ls$objective, type = "l",
xlab = "iteration t", ylab = expression(Q(theta^(t))))
3.4 Choosing
lambda by cross-validation
ppm_tune() performs \(K\)-fold cross-validation, selecting the
lambda that maximizes the held-out distance correlation
between the response and the estimated sufficient predictors. Call
set.seed() beforehand for reproducible folds. The
line-search settings are passed to every cross-validation fit.
set.seed(1)
cv <- ppm_tune(x, y, loss = "lssvm", d = 2, n.fold = 5,
nlambda = 10, lambda.max = 0.02)
cv$opt.lambda
#> [1] 0.009283178
summary(cv$fit, d = 2)
#> Penalized Principal Machine (P2M) summary
#> Loss: lssvm Penalty: grSCAD lambda = 0.00928318
#> Working dimension d = 2
#>
#> Estimated basis of the central subspace:
#> Dir1 Dir2
#> x1 1e+00 -3e-04
#> x2 3e-04 1e+00
#> x3 0e+00 0e+00
#> x4 0e+00 0e+00
#> x5 0e+00 0e+00
#> x6 0e+00 0e+00
#> x7 0e+00 0e+00
#> x8 0e+00 0e+00
#> x9 0e+00 0e+00
#> x10 0e+00 0e+00
#>
#> Selected variables (2 of 10): x1, x2