reml_grad_from_dV: Evaluate the REML gradient with respect to a single...

In lgspline: Lagrangian Multiplier Smoothing Splines for Smooth Function Estimation

Evaluate the REML gradient with respect to a single correlation parameter

Description

Computes the derivative of the negative REML objective with respect to a scalar correlation parameter, given the matrix derivative \partial \mathbf{V}/\partial \rho.

Usage

reml_grad_from_dV(dV, model_fit, glm_weight_function, ...)

Arguments

dV	N \times N numeric matrix giving \partial\mathbf{V}/\partial\rho for one scalar parameter \rho, with the chain rule through any reparameterization already applied.
model_fit	A fitted lgspline object; see lgspline.
glm_weight_function	The GLM weight function used during model fitting; see the glm_weight_function argument of lgspline.
...	Additional arguments forwarded to glm_weight_function.

Details

Notation

\mathbf{D} = \mathrm{diag}(d_1,\ldots,d_N) is the diagonal matrix of observation weights (observation_weights), and \mathbf{W} = \mathrm{diag}(w_1,\ldots,w_N) is the diagonal matrix of GLM working weights evaluated at the current fitted values via glm_weight_function, with the observation-weight contribution carried separately by \mathbf{D}. For canonical GLM families, w_i is the usual IRWLS/Fisher-scoring weight on the mean–variance scale; for example, logistic regression gives w_i = \mu_i(1-\mu_i). The combined weighting entering the information matrix is \mathbf{W}\mathbf{D}. In Gaussian identity models both reduce to scalar multiples of the identity.

\mathbf{V} is the N \times N correlation matrix implied by VhalfInv, with \mathbf{V}^{-1} = (\mathbf{V}^{-1/2})^\top \mathbf{V}^{-1/2}.

The penalized observed information at the current iterate is

\mathbf{M} = (\mathbf{X}^*)^\top \mathbf{V}^{-1}\mathbf{W}\mathbf{D}\,\mathbf{X}^* + \mathbf{U}^\top\boldsymbol{\Lambda}\mathbf{U},

where \mathbf{X}^* = \mathbf{X}\mathbf{U} is the constrained design (N \times P) and the first term is the quadratic approximation to the penalized log-likelihood Hessian at the current \boldsymbol{\mu}. For non-Gaussian families this is a local approximation (the IRWLS/Fisher scoring Hessian), not an exact GLS information matrix.

The constraint projection is \mathbf{U} = \mathbf{I} - \mathbf{G}\mathbf{A}(\mathbf{A}^\top\mathbf{G}\mathbf{A})^{-1} \mathbf{A}^\top with \mathbf{G} = \mathbf{M}^{-1}. \mathbf{U} is idempotent (\mathbf{U}^2 = \mathbf{U}) but not symmetric, so \mathbf{U}^\top\boldsymbol{\Lambda}\mathbf{U} \neq \mathbf{U}\boldsymbol{\Lambda}\mathbf{U}^\top.

The U stored in model_fit$U is on the expansion- and response-standardised scale. The rescaled version used here is

\tilde{\mathbf{U}} = \mathbf{U} \cdot \mathrm{diag}(1, s_1, \ldots, s_{p-1}, 1, s_1, \ldots) / \hat{\sigma}_y,

where s_j are the expansion scales and \hat{\sigma}_y standardises the response. All quantities below use \tilde{\mathbf{U}} in place of \mathbf{U}.

REML objective and its gradient

The REML objective is constructed by integrating out the fixed effects from the penalized log-likelihood, using a Laplace approximation to the marginal likelihood for non-Gaussian families. This approximation is exact for Gaussian identity models and is the standard extension used in restricted maximum likelihood estimation for GLMMs.

The REML correction term is -\tfrac{1}{2}\log|\mathbf{M}|, where \mathbf{M} is the penalized observed information defined above. Differentiating with respect to \rho (noting only \mathbf{V} depends on \rho) gives the REML correction gradient

-\frac{1}{2}\mathrm{tr}\!\Bigl( \mathbf{M}^{-1} (\mathbf{X}^*)^\top\mathbf{V}^{-1} \frac{\partial\mathbf{V}}{\partial\rho} \mathbf{V}^{-1}\mathbf{W}\mathbf{D}\,\mathbf{X}^* \Bigr).

Full gradient

\frac{\partial(-\ell_R)}{\partial\rho} = \frac{1}{N}\Biggl[ \underbrace{ \frac{1}{2}\mathrm{tr}\!\Bigl( \mathbf{V}^{-1}\frac{\partial\mathbf{V}}{\partial\rho} \Bigr) }_{\text{(i) log-det of }\mathbf{V}} \underbrace{ -\frac{1}{2\tilde{\sigma}^2} \mathbf{r}^\top\frac{\partial\mathbf{V}}{\partial\rho}\mathbf{r} }_{\text{(ii) residual quadratic form}} \underbrace{ -\frac{1}{2}\mathrm{tr}\!\Bigl( \mathbf{M}^{-1}(\mathbf{X}^*)^\top\mathbf{V}^{-1} \frac{\partial\mathbf{V}}{\partial\rho} \mathbf{V}^{-1}\mathbf{W}\mathbf{D}\,\mathbf{X}^* \Bigr) }_{\text{(iii) REML correction}} \Biggr],

where the whitened residual is

\mathbf{r} = \mathrm{diag}\!\left(\sqrt{d_i}/\sqrt{w_i}\right) \mathbf{V}^{-1/2}(\mathbf{y} - \boldsymbol{\mu}),

and r_i = [\mathbf{V}^{-1/2}(\mathbf{y}-\boldsymbol{\mu})]_i \sqrt{d_i}/\sqrt{w_i}.

Term (i) does not involve \mathbf{D} or \mathbf{W}; the log-determinant of \mathbf{V} depends only on the correlation structure. For non-Gaussian families terms (ii) and (iii) are evaluated at the current IRWLS iterate and constitute a local approximation.

Value

A scalar: the gradient of the negative REML objective with respect to \rho, divided by N.

See Also

lgspline, lgspline.fit

lgspline documentation built on Aug. 5, 2026, 1:10 a.m.