Using Shrinkage-Based Methods for Covariate Selection

Why Use Shrinkage Before Refitting?

A covariate-rich NLME model can contain many plausible relationships, some of which are weakly supported by the data. Estimating all candidate effects without constraint can leave noisy coefficients that are difficult to distinguish. A shrinkage-based strategy fits a model with a penalty that pulls covariate coefficients toward zero, uses those penalized estimates to form a shortlist, and then refits the shortlist without the penalty. The purpose of the first fit is to prioritize candidates, not to provide the final effect estimates.

This approach concerns fixed-effect covariate coefficients. It is distinct from using empirical Bayes estimates of individual random effects as screening tools: the penalized quantities here describe population-level covariate relationships, not individual deviations. The earlier material on empirical Bayes estimates is therefore relevant background, but the two kinds of estimates answer different questions.

Definition: A shrinkage-based covariate selection strategy fits a model whose objective includes a penalty on candidate covariate coefficients, uses the resulting shrunken coefficients to identify a reduced set of candidates, and refits that reduced model without the penalty. In a ridge-style version, the penalty is proportional to the sum of squared candidate coefficients. Ridge shrinkage generally makes coefficients smaller but does not make them exactly zero.

Model and Assumptions

Let \(y_{ij}\) be the observation for individual \(i\) at time \(t_{ij}\), and let \(x_i\) denote that individual's covariates. A general NLME model can be written as

y_{ij}=f(t_{ij},\theta_i)+\epsilon_{ij},\qquad \eta_i\sim N(0,\Omega).

Here \(f\) is the structural model, \(\theta_i\) contains individual parameters, \(\eta_i\) represents individual random effects, and \(\epsilon_{ij}\) follows the specified residual-error model. For example, a positive parameter \(\theta_{ik}\) may have a log-linear covariate model,

\log(\theta_{ik})=\beta_{k0}+\sum_{m=1}^{M}\beta_{km}z_{im}+\eta_{ik}.

In this expression, \(z_{im}\) is candidate covariate \(m\), \(\beta_{k0}\) is the baseline log parameter, and \(\beta_{km}\) is its covariate coefficient. The selection procedure assumes that the structural and residual models are adequate for the question, the likelihood can be fit reliably, and candidate covariates are coded meaningfully. It does not correct misspecification of those models, confounding, or inadequate covariate support.

For a ridge-style screen, a common penalized objective is

Q_{\lambda}(\beta,\psi)=-2\ell(\beta,\psi)+\lambda\sum_{m=1}^{M}w_m\beta_m^2,\qquad \lambda\geq 0.

Here \(\ell\) is the NLME log likelihood, \(\psi\) collects other model parameters, \(\beta_m\) are the candidate covariate coefficients being screened, and \(w_m\) are fixed penalty weights. Baseline parameters and other nuisance parameters are typically not penalized. The precise objective and optimization method should be held consistent across the screen. Because coefficient scales depend on covariate coding, continuous covariates are commonly centered and scaled before fitting; categorical covariates require a declared reference level and coding. Otherwise, the penalty may reflect measurement units rather than evidence for a relationship.

What the Penalty Does

The shrinkage effect can be seen in a local quadratic approximation to the likelihood. Consider one coefficient \(\beta\), with other parameters treated as fixed for this justification. Let \(\hat{\beta}\) be its unpenalized estimate, and let \(I\) be the curvature of \(-2\ell\) at that estimate. If the negative twice log likelihood is approximately quadratic near \(\hat{\beta}\), then

-2\ell(\beta)\approx C+I(\beta-\hat{\beta})^2,

where \(C\) does not depend on \(\beta\). Adding the ridge penalty gives \(C+I(\beta-\hat{\beta})^2+\lambda\beta^2\). Differentiating with respect to \(\beta\), setting the derivative to zero, and solving yields

2I(\beta-\hat{\beta})+2\lambda\beta=0 \quad\Longrightarrow\quad \tilde{\beta}_{\lambda}=\frac{I}{I+\lambda}\hat{\beta}.

For \(I>0\) and \(\lambda\geq0\), the factor \(I/(I+\lambda)\) lies in \((0,1]\). For nonzero \(\hat{\beta}\), the approximate penalized estimate has the same sign and no greater magnitude; its magnitude is strictly smaller when \(\lambda>0\). At \(\lambda=0\), it is unchanged. If \(\hat{\beta}=0\), the penalized estimate is also zero, and describing the estimates as having the same sign is not meaningful. In a full nonlinear mixed-effects fit, the approximation is only local: parameters interact, and numerical estimates need not obey this one-coefficient formula exactly.

Worked Example: Calculating Approximate Shrinkage

Suppose a hypothetical unpenalized fit estimates a covariate coefficient as \(\hat{\beta}=0.30\) per 10-unit increase in a covariate, and the local curvature is \(I=25\). Use a ridge penalty with \(\lambda=25\). The local approximation gives

\tilde{\beta}_{\lambda} =\frac{I}{I+\lambda}\hat{\beta} =\frac{25}{25+25}(0.30) =0.50(0.30)=0.15.

The estimate is reduced from 0.30 to 0.15 log-parameter units per 10-unit increase. Under the log-linear parameterization, the corresponding approximate multiplicative change is \(\exp(0.15)=1.162\), or about a 16.2% increase in the parameter per 10-unit increase. This is a screening estimate under the penalty, not the final unpenalized estimate. The result also shows why coefficients with less information may be shrunk more for the same \(\lambda\): the information-to-penalty ratio is smaller.

A Practical Screening and Refit Workflow

The screen is most useful when it is treated as a disciplined way to reduce a candidate set, not as an automatic test of whether an effect exists. The full covariate model should be specified from scientifically plausible candidates, with coding and scaling decided before interpreting coefficient magnitudes. Then fit the penalized model, inspect the shrunken effects, and define a shortlist rule in advance or at least document it. The rule may use standardized coefficient magnitude together with direction, uncertainty, and scientific plausibility; a numerical cutoff is a screening convention, not a significance threshold.

Specify and prepare candidates.
Use a defensible structural and residual model, encode covariates consistently, and scale continuous candidates so their penalization is comparable.

Fit the penalized model.
Apply the penalty to the candidate coefficients, not indiscriminately to baseline or nuisance parameters. Check convergence and whether estimates are sensitive to the chosen penalty strength.

Form a shortlist.
Use a stated screening rule to identify candidates worth retaining. Do not treat a nonzero ridge estimate as proof of an effect: ridge generally leaves every coefficient nonzero.

Refit without the penalty.
Fit the reduced NLME model using the original likelihood, then assess estimates, uncertainty, model adequacy, and whether the retained relationships make sense in context.

Worked Example: Screening Three Standardized Covariates

Suppose a hypothetical model has three continuous candidate covariates, each standardized to mean 0 and standard deviation 1. The following are approximate unpenalized coefficients and local curvatures from the full model. Use \(\lambda=20\) and a screening rule that retains a candidate when the absolute approximate penalized coefficient is at least 0.08 log-parameter units per standard deviation.

Candidate	\(\hat{\beta}\)	\(I\)	Approximate factor \(I/(I+\lambda)\)	\(\tilde{\beta}_{\lambda}\)
A	0.24	40	40/60 = 0.667	0.160
B	0.18	10	10/30 = 0.333	0.060
C	-0.12	30	30/50 = 0.600	-0.072

For A, the calculation is \((40/60)(0.24)=0.160\), which exceeds 0.08 in absolute value. For B it is \((10/30)(0.18)=0.060\), below the cutoff. For C it is \((30/50)(-0.12)=-0.072\), also below the cutoff. Under this explicitly chosen rule, only A is shortlisted. This result is not a hypothesis test: the cutoff has no automatic false-positive rate, and changing \(\lambda\), scaling, or the rule can change the shortlist.

Suppose the reduced unpenalized refit estimates A as 0.21 with standard error 0.09. The coefficient is per one standard deviation increase in A. Its approximate 95% Wald interval on the log scale is \(0.21\pm1.96(0.09)=[0.0336,0.3864]\). Exponentiating gives \([1.034,1.472]\), so the fitted parameter is estimated to be about 3.4% to 47.2% higher per standard deviation increase in A, under this model. Because A was selected using the same data, this ordinary interval does not account for the selection step and should not be presented as fully selection-adjusted uncertainty.

Worked Example: Why Correlated Candidates Can Be Hard to Separate

Consider a hypothetical limiting case in which two standardized candidate covariates are identical for every individual, \(z_{i1}=z_{i2}\). Suppose the likelihood in the relevant direction depends only on their sum, \(\beta_1+\beta_2=0.30\). With equal ridge penalties, the penalty among coefficient pairs with that sum is proportional to \(\beta_1^2+\beta_2^2\). Substitute \(\beta_2=0.30-\beta_1\):

\beta_1^2+(0.30-\beta_1)^2 =2(\beta_1-0.15)^2+0.045.

This expression is minimized at \(\beta_1=0.15\), so \(\beta_2=0.15\) as well. The penalized fit divides the combined signal equally, but the data cannot distinguish the two covariate relationships because the covariates are identical. A shortlist that keeps both coefficients based on their sizes would not establish two separate effects. With highly correlated rather than identical covariates, the allocation may still be unstable. The screen must therefore be interpreted alongside the candidate support and collinearity checks, not as a way to resolve information the study does not contain.

Diagnostics, Alternatives, and Limitations

Before using the shortlist, check whether the penalized optimization converged, whether the unpenalized and penalized fits use the same data and likelihood specification, and whether the fitted model gives plausible predictions and residual behavior. Large changes under modest changes in penalty strength, implausible signs, or coefficients that are extremely sensitive to covariate coding are warnings that the ranking is not robust. Sparse support, strong correlation, confounding, and model misspecification can all undermine the interpretation even when the optimizer reports convergence.

The main alternatives answer related but different questions. Likelihood-ratio or Wald-based covariate procedures assess specified effects in unpenalized models, while a shrinkage screen prioritizes candidates by estimates deliberately pulled toward zero. Direct unpenalized forward or backward searches make inclusion decisions through a sequence of model comparisons; a ridge screen instead estimates many candidates together under a penalty and then requires a separate rule to reduce the model. These methods do not eliminate the need to assess the final model, and their results need not agree.

A ridge penalty is useful when the aim is to damp noisy coefficients and create a manageable shortlist, but it has important limits. It does not itself provide a calibrated inclusion probability, a conventional hypothesis test, or a confidence interval that accounts for selection. It may retain weak effects as small nonzero estimates, and correlated candidates can share or redistribute apparent signal. The final unpenalized refit restores coefficients to the original model scale, but it does not undo uncertainty introduced by having selected candidates using the same data. Report the penalty, scaling, shortlist rule, and refit clearly, and regard the refit as conditional on the selection procedure.

Common pitfall: A learner may conclude that a covariate is important because its ridge estimate is nonzero, or may interpret the shrunken coefficient as the final effect estimate. That reasoning is incorrect: ridge generally leaves coefficients nonzero, and the penalty intentionally changes their magnitudes. Use the penalized fit for screening under a stated rule, then estimate and interpret retained effects from the unpenalized reduced model while acknowledging selection-related uncertainty.

Key takeaway: Shrinkage-based covariate selection is a two-stage strategy: penalize candidate coefficients to prioritize them, then refit a reduced NLME model without the penalty. The screen is a ranking aid, not proof of a relationship; the retained effects require contextual assessment and careful uncertainty interpretation.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

In a ridge-style NLME screen, which quantities are penalized, and why are baseline parameters commonly left unpenalized?
If \(I=12\), \(\lambda=8\), and \(\hat{\beta}=-0.25\), calculate the approximate penalized estimate.
Why should continuous covariates be scaled before applying a common ridge penalty?
What does a nonzero ridge coefficient establish about the covariate effect?
Why does refitting without the penalty not remove all uncertainty caused by screening?

Solutions. (1) Candidate covariate coefficients are penalized; leaving baseline parameters unpenalized avoids shrinking the reference level or overall scale toward zero. (2) The factor is \(12/(12+8)=0.600\), so the estimate is \(0.600(-0.25)=-0.150\). (3) Scaling makes coefficient magnitudes more comparable, so the penalty is less dependent on arbitrary measurement units. (4) Nothing definitive about existence: ridge generally produces nonzero estimates, and a screening rule is needed. (5) The same data helped determine which covariates entered the refit, so ordinary refit uncertainty does not account for that selection process.