library("gamlss2")
## longley data example with
## strongly correlated predictors
data("longley")
## fit all predictors linearly
m1 <- gamlss2(
Employed ~ lin(~ GNP.deflator + GNP + Unemployed +
Armed.Forces + Population + Year),
data = longley,
family = NO
)
## fit the same model using a ridge-penalized block,
## ridge() standardizes the model-matrix columns before applying the
## common ridge penalty.
m2 <- gamlss2(
Employed ~ ridge(~ GNP.deflator + GNP + Unemployed +
Armed.Forces + Population + Year),
data = longley,
family = NO
)
## Compare model summaries.
summary(m1)
summary(m2)
## extract coefficients
cm <- cbind(
"linear" = specials(m1, model = "mu", elements = "coefficients"),
"ridge" = specials(m2, model = "mu", elements = "coefficients")
)
print(cm)
## ridge terms combined with ordinary linear effects,
## Year is left unpenalized, while the remaining correlated
## economic variables are shrunk jointly
m3 <- gamlss2(
Employed ~ Year +
ridge(~ GNP.deflator + GNP + Unemployed +
Armed.Forces + Population),
data = longley,
family = NO
)
summary(m3)
## interactions
m4 <- gamlss2(
Employed ~ ridge(~ GNP + Unemployed + GNP:Unemployed),
data = longley,
family = NO
)
specials(m4, model = "mu", elements = "coefficients")Bundled Linear and Ridge-Penalized Model Terms
Description
Constructs a model-matrix block for use as a single model term in a gamlss2 formula. Function lin() fits the block with only a negligible numerical penalty, whereas ridge() standardizes its columns and estimates a ridge penalty for the complete block.
Usage
lin(x, ..., ridge = FALSE, scale = FALSE)
ridge(...)
Arguments
x
|
A one-sided formula or the first unquoted predictor expression. A formula is recommended when transformations, factors, or interactions are used, for example |
…
|
In |
ridge
|
Logical. Should the identity-penalty smoothing parameter be estimated? This is primarily a low-level argument; |
scale
|
Logical. Should the columns of the constructed design matrix be centered and divided by their sample standard deviations? Columns with zero or non-finite standard deviation are left unchanged. Function |
Details
Model-matrix construction. The supplied predictors are converted to a design matrix using model.matrix. Any column named “(Intercept)” is removed because the distributional predictor has its own intercept. Apart from that removal, standard R formula rules and contrasts are used. Thus, ~ x1 * x2 expands to two main effects and their interaction, while ~ x1:x2 produces a pure interaction without adding the marginal terms. Pure interactions are therefore supported. For a factor f, ~ f uses the current contrast coding (normally treatment coding), whereas ~ 0 + f requests one indicator column per level.
Supplying several predictors without a formula is a convenience shortcut: lin(x1, x2) is equivalent to lin(~ x1 + x2), and similarly for ridge(). A formula should be used whenever operators such as :, *, -, or 0 + need to be interpreted explicitly.
Unpenalized linear blocks. Function lin() fixes the identity-penalty multiplier at a negligible value. It therefore behaves like weighted linear regression for the current working response, while entering all columns as one model term rather than as separate ordinary linear terms. This can be useful for keeping a block of columns together and for using the same model-matrix and prediction machinery as penalized terms. Ordinary linear effects can still be specified directly in the surrounding gamlss2 formula.
Ridge-penalized blocks. For ridge(), coefficients are estimated from the current weighted working response using an identity penalty. Ignoring the small numerical stabilization used by the fitting algorithm, the update minimizes
\((z - X\beta)^\mathsf{T} W (z - X\beta) + \lambda\,\beta^\mathsf{T}\beta,\)
where \(z\) and \(W\) are the current working response and weights. A large \(\lambda\) shrinks the whole coefficient block toward zero.
Before fitting, every nonconstant column of \(X\) is standardized using its training-data mean and standard deviation. The stored transformation is reused for prediction. Consequently, the identity penalty is applied on a comparable scale across numeric predictors, indicator columns, and interaction columns.
The smoothing parameter \(\lambda\) is selected during backfitting using the smooth-term criterion supplied to gamlss2(). The default is corrected AIC (criterion = “aicc”); supported alternatives include “gcv”, “aic”, “gaic”, and “bic”. For example, use gamlss2(…, criterion = “bic”) to select the ridge penalty by BIC. The previous estimate is used to warm-start subsequent backfitting updates.
The effective degrees of freedom are computed from the penalized weighted hat matrix.
Value
lin() and ridge() return a smooth specification object of class “lin.smooth.spec”, which is completed to an object of class “lin.effect” during model setup.
See Also
gamlss2, specials, model.matrix, smooth.construct